Extension of the Henry semi-analytical solution for saltwater intrusion in stratified domains

The Henry semi-analytical solution is developed for stratified aquifers with exponential and Gaussian permeability–depth relationships and small dispersion. The semi-analytical solution is sought by expanding the stream function and the concentration in an infinite Fourier series truncated at given orders. Due to the heterogeneity, the semi-analytical solution contains additional terms and quickly becomes impractical because of the high truncation orders required to yield stable results. An appropriate evaluation of the most expensive term, involving fairly complex summations, is proposed to render the computation of the semi-analytical solution affordable. Three configurations with the same transmissivity as that for the homogeneous Henry problem are investigated. The first configuration considers an exponential decay of the permeability (k) with depth. The second and third configurations consider an exponential increase and a Gaussian high-low-high distribution of k with depth, respectively. The three configurations are also investigated numerically using an accurate numerical code based on the method of lines and advanced spatial discretization schemes. An excellent agreement is obtained between the semi-analytical and the numerical solutions for all test cases. The results show that the amount of saltwater intrusion is strongly dependent on the heterogeneity. Further, the heterogeneous Henry problem is found to be more suitable for benchmarking density-driven numerical codes because the numerical solution is more sensitive to the grid size than for the homogeneous problem.     

双层饱和软土地基一维大应变固结研究

研究了双层软土地基一维非线性大应变固结问题,同时考虑了土体受沉积作用影响的自重分布以及渗透性在固结过程中的变化等因素。对此复杂问题,利用半解析法进行了求解。通过与已有的大应变固结解析解对比,检验了半解析法在大应变固结分析中的有效性。最后较为详细地分析了双层软土地基一维大应变固结性状,并与相应条件下的小应变固结作了比较。     

Application of empirical and semi-analytical algorithms to MERIS data for estimating chlorophyll a in Case 2 waters of Lake Erie

Lake Erie is biologically the most active lake among the Great Lakes of North America, experiencing seasonal harmful algal blooms (HABs). The early detection of HABs in the Western Basin of Lake Erie (WBLE) requires a more efficient and accurate monitoring tool. Remote sensing is an efficient tool with high spatial and temporal coverage that can allow accurate and timely detection of the HABs. The WBLE is heavily influenced by the surrounding terrestrial ecosystem via rivers such as the Sandusky River and the Maumee River. As a result, the optical properties of the WBLE are influenced by multiple color producing agents (CPAs) such as phytoplankton, colored dissolved organic matter (CDOM), organic detritus, and terrigenous inorganic particles. The diversity of the CPAs and their non-linear interactions makes these waters optically complex, and the task of optical remote sensing for retrieving estimates of CPAs more challenging. Chlorophyll a, which is the primary light harvesting pigment in all phytoplankton, is used as a proxy for algal biomass. In this study, several published remote sensing algorithms and band ratio models were applied to the reflectance data from the full resolution MERIS sensor to remotely estimate chlorophyll a concentrations in the WBLE. Efficiency of the sensor and the algorithms performance were tested through a least squares regression and residual analysis. The results indicate that, among the suite of existing bio-optical models, the Simis semi-analytical algorithm provided the best model results for measures of algal biomass in the optically complex WBLE with R ² of 0.65, RMSE 0.85 μg/l, (n = 71, P < 0.05). The superior results of this model in detecting chlorophyll a are attributed to several factors including optimizing spectral regions that are less sensitive to CDOM and the incorporation of correction factors such as absorption effects due to pure water (a _w), backscatter (b _b) from suspended matter and interference due to phycocyanin (δ), a major accessory pigment in the WBLE.     

A general semi-analytical solution for consolidation around an expanded cylindrical and spherical cavity in modified Cam Clay

This paper presents a general semi-analytical solution for undrained cylindrical and spherical cavity expansion in Modified Cam Clay (MCC) and subsequent consolidation. The undrained cylindrical and spherical cavity expansion response in MCC model is obtained through the similarity solution technique. Then, the subsequent consolidation process around the cavity is governed by the classical partial differential equation for consolidation. Finite Difference Method (FDM) is selected for solving the consolidation equation numerically. The proposed semi-analytical solution is validated by comparing the prediction of the dissipations of the pore pressure with Randolph’s closed-form solution for elastic-perfectly plastic soil. Parametric study shows that G₀/p₀′, R and M have significant influence on the cavity wall excess pore pressure dissipation curve, while it is not sensitive to the value of ν′. It is also found that the negative pore pressure generates around the expanded cylindrical and spherical cavity wall during the consolidation process when R > 5 for typical Boston blue clay. The developed solution has potential applications in geotechnical problems, such as the pile foundation, in-situ test, tunnel construction, compaction grouting, and so forth.     

Nonlinear consolidation in randomly heterogeneous highly compressible aquitards

Severe land subsidence due to groundwater extraction may occur in multiaquifer systems where highly compressible aquitards are present. The highly compressible nature of the aquitards leads to nonlinear consolidation where the groundwater flow parameters are stress-dependent. The case is further complicated by the heterogeneity of the hydrogeologic and geotechnical properties of the aquitards. The effect of realistic vertical heterogeneity of hydrogeologic and geotechnical parameters on the consolidation of highly compressible aquitards is investigated by means of one-dimensional Monte Carlo numerical simulations where the lower boundary represents the effect of an instant drop in hydraulic head due to groundwater pumping. Two thousand realizations are generated for each of the following parameters: hydraulic conductivity (K), compression index (C _c), void ratio (e) and m (an empirical parameter relating hydraulic conductivity and void ratio). The correlation structure, the mean and the variance for each parameter were obtained from a literature review about field studies in the lacustrine sediments of Mexico City. The results indicate that among the parameters considered, random K has the largest effect on the ensemble average behavior of the system when compared to a nonlinear consolidation model with deterministic initial parameters. The deterministic solution underestimates the ensemble average of total settlement when initial K is random. In addition, random K leads to the largest variance (and therefore largest uncertainty) of total settlement, groundwater flux and time to reach steady-state conditions.     

Spatial analysis of the frequency–magnitude distribution of aftershock activity of 2003 Bam earthquake: southeast Iran

The Bam earthquake (2003 December 26, M _W = 6.6) was one of the largest earthquakes that occurred in southeast of Iran during last century. It took place along an N–S trending right-lateral strike-slip fault, almost near the southern end of Nyband–Gowk fault. In this study, we mapped the frequency–magnitude distribution of aftershock events spatially across the Bam aftershock zone. The b-value varies between 0.6 and 1.1 across the Bam rupture zone. The overall depth distribution of b-value in Bam aftershock zone reveals two distinct increases in b-value: (1) at depths of 8–10 km and (2) shallower than 4 km beneath the Bam city. There is no correlation between high b- value anomalies found in this study and the region of largest slip, whereas the spatial correlation between high b-value anomalies and the zone of low V _s and high σ (in earlier tomography study) is obvious. This correlation reveals that material properties and increasing heterogeneity are more important in controlling b-value distribution in Bam earthquake rupture zone. The high b-value anomaly near the surface of northern part of rupture zone may be related to unconsolidated and water-rich quaternary alluvial sediments and probable low-strength rocks beneath them. The high b-value anomaly at depth range 8–10 km can be correlated with fractured and fluid-filled mass, which may result from the movement of magma during Eocene volcanism in the Bam area. In this study, the induced changes in pore fluid pressure due to main shock are suggested as a mechanism for aftershock generation.     

The influence of porosity gradients on mixing coefficients in sediments

The surface layer of aquatic sediments is a zone characterized by both porosity gradients and intensive mixing. In the standard approach, porosity gradients are ignored when estimating mixing intensity. Here, model formulations with both constant and varying porosity are contrasted to estimate mixing coefficients D_b from tracer depth profiles. Complementing the well-known exponential solution of the constant-porosity model, we present a general solution to the variable-porosity model in terms of hypergeometric functions. When using these models in a forward way, the tracer activities predicted by the variable-porosity model are higher than those generated by the constant-porosity model. Similarly, when inverse modelling, D_b values estimated by the variable-porosity model are systematically higher than those derived from the constant-porosity model. Still, differences in D_b values remain relatively small. When applying both mixing models to excess ²¹⁰Pb data profiles from slope sediments, a maximal difference of 30% is obtained between D_b values, the average deviation being 16%. A systematic exploration of parameter space predicts a maximal underestimation of 60% when deriving D_b values from the constant-porosity mixing model. Given the uncertainty imposed by other model assumptions underlying the diffusive mixing model, the influence of porosity gradients on D_b values must be classified as rather modest. Hence, the current mixing coefficient database is not biased by the constant porosity approximation.     

An application of Green’s function technique for computing well inflow without radial flow assumption

Well modeling plays an important role in numerical reservoir simulation. The main difficulty in well modeling is the difference in scale between the wellbore radius and well gridblock dimension used in the simulation. The Peaceman equation is widely used in reservoir simulation to match gridblock pressure to the local solution of the diffusivity equation describing the flow near the well. However, this approach was developed under the assumption of radial flow. At the same time, the well inflow equation can be solved within the Green’s function (GF) formalism which allows the solution to be obtained without the assumption of radial flow. The GF solution can be presented as a series over the eigenvalues of the Laplace differential operator. However, this series converges conditionally and its direct summation is time-consuming. In Posvyanskii et al. (2008), a method for fast summation of such a series was proposed and successfully applied for analyzing the pressure build up curves. In this paper, we adopt the same technique for calculating the well indices for horizontal, slanted and partially penetrated wells. Additionally, the role of different boundary conditions is considered. The semi-analytical expressions for well indices are obtained and compared to the solution of the Peaceman equation. It is shown that in some cases, the difference between these solutions can be significant. The use of the obtained expression in numerical flow simulation allows well inflow to be modeled with high accuracy even on a coarse grid.     

Design and performance of a novel multi-function debris flow mitigation system in Wenjia Gully,Sichuan

The post-earthquake debris flows in the Wenjia Gully led to the exposure of the shortcomings in the design of the original conventional debris flow mitigation system. A predicament for the Wenjia mitigation system is a large amount of loose material (est. 50 × 10⁶ m³) that has been deposited in the gully by the co-seismic landslide, providing abundant source material for debris flows under saturation. A novel design solution for the replacement mitigation system was proposed and constructed, and has exhibited excellent performance and resilience in subsequent debris flows. The design was governed by the three-phase philosophy of controlling water, sediment, and erosion. An Early Warning System (EWS) for debris flow that uses real-time field data was developed; it issues alerts based on the probabilistic and empirical correlations between rainfall and debris flows. This two-fold solution reduces energy of the debris flow by combining different mitigation measures while minimizing the impact through event forecasting and rapid public information sharing. Declines in the number and size of debris flows in the gully, with increased corresponding rainfall thresholds and mean rainfall intensity-duration (I-D) thresholds, indicate the high efficacy of the new mitigation system and a lowered debris flow susceptibility. This paper reports the design of the mitigation system and analyzes the characteristics of rainfall and debris flow events that occurred before and after implementation of the system; it evaluates the effectiveness of one of the most advanced debris flow mitigation systems in China.     

A novel nonlinear solution for the polygon scaled boundary finite element method and its application to geotechnical structures

The polygon scaled boundary finite element method is semi-analytical and known for its high precision. However, the material nonlinearity cannot be maintained because this method uses an analytical solution in the radial direction. In this paper, a novel nonlinear algorithm is developed by introducing internal Gaussian points over a subdomain. The response of nonlinearity for a concrete-faced rockfill dam is modeled. The results correspond well with the results from finite element modelling, which demonstrates the method can be used to describe the nonlinear characteristics of geomaterials. Furthermore, this method offers promising flexibility for analyzing complex geometries without decreasing the precision.     

Modeling of shear failure in fractured reservoirs with a porous matrix

A finite volume-based numerical modeling framework using a hierarchical fracture representation (HFR) has been developed to compute flow-induced shear failure. To accurately capture the mechanics near fracture manifolds, discontinuous basis functions are employed which ensure continuity of the displacement gradient across fractures. With these special basis functions, traction and compressive forces on the fracture segment can be calculated without any additional constraints, which is extremely useful for estimating the irreversible displacement along the fracture (slip) based on a constitutive friction law. The method is further extended to include slip-dependent hydraulic aperture change and grid convergent results are obtained. Further, the change in hydraulic aperture is modeled using an asymptotic representation which respects the experimentally observed behavior of pore volume dilation due to shear slip. The model allows the initial rapid increase in hydraulic aperture due to shear slip and asymptotically approaches a finite value after repeated shearing of a fracture segment. This aperture increase is the only feedback for mechanics into the fluid flow for a linear elastic mechanics problem. The same model is also extended to include poroelastic relations between flow and mechanics solver. The grid convergence result in the case of poroelastic flow-mechanics coupling for flow-induced shear failure is also obtained. This proves the robustness of the numerical and analytical modeling of fracture and friction in the extended finite volume method (XFVM) set-up. Finally, a grid convergent result for seismic moment magnitude for single fracture and fracture network with random initial hydraulic and friction properties is also obtained. The b-value, which represents the slope of seismic moment occurrence frequency decay vs seismic moment magnitude, which is approximately constant in a semi-logarithmic plot, is estimated. The numerical method leads to converged b-values for both single fracture and fracture network simulations, as grid and time resolutions are increased. For the resulting linear system, a sequential approach is used, that is, first, the flow and then the mechanics problems are solved. The new modeling framework is very useful to predict seismicity, permeability, and flow evolution in geological reservoirs. This is demonstrated with numerical simulations of enhancing a geothermal system.     

指数形式渗流定律下软土非线性大变形固结分析

考虑土中指数形式渗流定律以及土体的非线性固结特性,以超静孔隙水压力为变量在拉格朗日坐标系内建立了软土一维大变形固结问题的控制方程及其求解条件,并运用有限差分法获取其数值解答。在指数形式渗流定律退化为达西定律下,通过将差分解与已有的半解析解进行对比,验证了数值计算的可靠性。最后对指数形式渗流定律下软土一维非线性大变形固结性状进行计算分析,结果表明： 1时,软土的非线性大变形固结速率会随外载增大而减慢; 1时,软土的非线性大变形固结速率会随着外荷载的增加而加快;软土非线性大变形固结速率要比非线性小变形固结速率快,且差别会随荷载增大而加剧;此外,大变形固结理论的最终沉降值要小于小变形固结理论,且差别会随着荷载的增大而加剧。     

Isotopic effects in fracture-dominated reactive fluid-rock systems

A mathematical model is presented that describes the effects of pore fluid aqueous diffusion and reaction rate on the isotopic exchange between fluids and rocks in reactive geo-hydrological systems where flow is primarily through fractures. The model describes a simple system with parallel equidistant fractures, and chemical transport in the matrix slabs between fractures by aqueous diffusion through a stagnant pore fluid. The solid matrix exchanges isotopes with pore fluid by solution-precipitation at a rate characterized by a time constant, R (yr⁻¹), which is an adjustable parameter. The effects of reaction on the isotopes of a particular element in the fracture fluid are shown to depend on the ratio of the diffusive reaction length for that element (L) to the fracture spacing (b). The reaction length depends on the solid-fluid exchange rate within the matrix, the partitioning of the element between the matrix pore fluid and the matrix solid phase, the porosity and density of the matrix, and the aqueous diffusivity. For L/b < 0.3, fluid-rock isotopic exchange is effectively reduced by a factor of 2L/b relative to a standard porous flow (single porosity) model. For L/b > 1, the parallel fracture model is no different from a porous flow model. If isotopic data are available for two or more elements with different L values, it may be possible to use the model with appropriate isotopic measurements to estimate the spacing of the primary fluid-carrying fractures in natural fluid-rock systems. Examples are given using Sr and O isotopic data from mid-ocean ridge (MOR) hydrothermal vent fluids and Sr isotopes in groundwater aquifers hosted by fractured basalt. The available data for MOR systems are consistent with average fracture spacing of 1-4 m. The groundwater data suggest larger effective fracture spacing, in the range 50-500 m. In general, for fractured rock systems, the effects of fracture-matrix diffusive exchange must be considered when comparing isotopic exchange effects for different elements, as well as for estimating water age using radioactive and cosmogenic isotopes.     

Thermal expansion,Debye temperature and Grüneisen parameter of synthetic (Fe,Mg)SiO3 orthopyroxenes

Thermal expansion properties of synthetic orthopyroxenes (Fe_0.20Mg_0.80)SiO₃, (Fe_0.40Mg_0.60)SiO₃, (Fe_0.50Mg_0.50)SiO₃, (Fe_0.75Mg_0.25)SiO₃ and (Fe_0.83Mg_0.17)SiO₃ were systematically studied by means of single-crystal x-ray diffraction in the temperature range from 296 to 1300 K. The measurements of unit cell dimensions as a function of temperature reveal that the a and c dimensions and the unit cell volume V increase nonlinearly with a positive curvature with rising temperature, whereas the b dimension behaves differently, depending on the total Fe content. For Mg-rich orthopyroxenes (Fe/(Fe+Mg)<30%), the b dimension expands similarly as the a and c dimensions, but it exhibits a nonlinear increase with a negative curvature for orthopyroxenes with Fe/(Fe+Mg)>30%. Together with the high temperature neutron diffraction data on enstatite (MgSiO₃) (McMullan, Haga and Ghose, unpublished) and x-ray diffraction data on ferrosilite (FeSiO₃) (Sueno et al. 1976), the measured unit cell dimensions were analyzed in terms of the Grüneisen theory of thermal expansion. The linear thermal expansion coefficients α _a and α _c both increase as temperature is elevated, with α _c increasing faster, while α _b changes gradually from increasing for Mg-rich orthopyroxenes to decreasing for Fe-rich orthopyroxenes. The relative magnitudes of linear thermal expansion coefficients are always in the order α _b >α _c >α _a between 300 and 500 K, but at higher temperatures, the order changes to α _c >α _b >α _a for Mg-rich orthopyroxenes and α _c >α _a >α _b for Fe-rich ones. The linear thermal expansion behavior is interpreted on the basis of the structural mechanical model of Weidner and Vaughan (1982). The anomalous behavior of α _b is mainly attributed to the changes in the Fe₂₊ population at the M2 site and the relative stiffness of the M2(Fe₂₊)-O bonds compared to the M2(Mg₂₊)-O bonds. The volume thermal expansion coefficients are nonlinear functions of temperature and lie between 23 and 49×10^?6/K. The previously reported results of mean volume thermal expansion coefficients appear to represent the α _V values characteristic of higher temperatures compared to our results. The thermal Debye temperatures are composition-dependent, decreasing linearly from 812 (MgSiO₃) to 561 K (FeSiO₃), and are systematically higher than the corresponding acoustic Debye temperatures. The Grüneisen parameters range from 0.85 to 0.89 and do not seem to vary with composition. The linear compressibilities derived from thermal expansion and elastic moduli data agree very well. The pressure derivatives of the isothermal bulk modulus (dK₀/dP) are also composition-dependent and decrease from 11.2 (MgSiO₃) to 8.77 (FeSiO₃). Such large values indicate possible anomalous elastic behavior of orthopyroxenes at high pressures in the Earth's upper mantle.     

Impulse response of viscous remanent magnetization: laboratory measurements by a pulse induction system

Transient electromagnetic responses measured in the field or in the laboratory may bear effects of viscous remanent magnetization (VRM) associated with magnetic relaxation of ultrafine grains of ferrimagnetic minerals or superparamagnetism. The behavior of VRM can be studied in time or frequency domain, TDEM measurements being advantageous because they are done in the absence of primary field and owing to broad time range providing high accuracy of VRM parameters. Another advantage is that the rate of viscous decay measured as voltage decay does not need to be corrected for stable and/or slowly decaying viscous component of total remanence. Time-dependent transient responses of viscous decay follow the power law a ⋅ t^–b, where a is the initial emf signal (varying in a broad range) and b is the exponent approaching 1. Laboratory tests with a pulse induction coil system reveal a strong linear correlation of the parameter a with frequency-dependent magnetic susceptibility Δκ used commonly for constraining the relative abundances of superparamagnetic particles. Such systems are thus suitable for quick measurements of the large number of samples for detection of superparamagnetic (SP) particles and quantifying their contribution. The difference of b from 1, though being minor, exceeds markedly its error in estimates from measured data. Simulated TDEM responses of a superparamagnetic ground show both parameters (a and b) to depend on particle volume distribution, which is prerequisite for inversion of time-domain transients to magnetic properties of rocks and soils.     

A model for silicate melt viscosity in the system CaMgSi2O6-CaAl2Si2O8-NaAlSi3O8

Five hundred eighty-five viscosity measurements on 40 melt compositions from the ternary system CaMgSi₂O₆ (Di)-CaAl₂Si₂O₈ (An)-NaAlSi₃O₈ (Ab) have been compiled to create an experimental database spanning a wide range of temperatures (660-2175°C). The melts within this ternary system show near-Arrhenian to strongly non-Arrhenian properties, and in this regard are comparable to natural melts. The database is used to produce a chemical model for the compositional and temperature dependence of melt viscosity in the Di-An-Ab system. We use the Vogel-Fulcher-Tammann equation (VFT: log η = A + B/(T − C)) to account for the temperature dependence of melt viscosity. We also assume that all silicate melts converge to a common viscosity at high temperature. Thus, A is independent of composition, and all compositional dependence resides in the parameters B and C. The best estimate for A is −5.06, which implies a high-temperature limit to viscosity of 10^-5.06 Pa s. The compositional dependence of B and C is expressed by 12 coefficients (b_i=1,2.6, c_j=1,2..6) representing linear (e.g., b_i=1:3) and higher order, nonlinear (e.g., b_i=4:6) contributions. Our results suggest a near-linear compositional dependence for B (<10% nonlinear) and C (<7% nonlinear). We use the model to predict model VFT functions and to demonstrate the systematic variations in viscosity due to changes in melt composition. Despite the near linear compositional dependence of B and C, the model reproduces the pronounced nonlinearities shown by the original data, including the crossing of VFT functions for different melt compositions. We also calculate values of T_g for melts across the Di-An-Ab ternary system and show that intermediate melt compositions have T_g values that are depressed by up to 100°C relative to the end-members Di-An-Ab. Our non-Arrhenian viscosity model accurately reproduces the original database, allows for continuous variations in rheological properties, and has a demonstrated capacity for extrapolation beyond the original data.     

X-ray and electron-optical investigation of synthetic high-temperature plagioclases

High-temperature plagioclases NaAlSi₃O₈(Ab)-CaAl₂Si₂O₈(An) have been prepared by dry devitrification of glasses. X-ray powder photographs were taken to follow the time-temperature dependence of indicators of structural state and to determine the lattice parameters of the plagioclases as well as their K-exchanged equivalents. Transmission electron microscopy (TEM) has been used to study the anorthite-type antiphase domain textures. The results are as follows:

1. TEM and X-ray powder data suggest to subdivide the high-temperature plagioclase series into three structural regions: (a) An0 to An30–40: monalbite/analbite to high albite-type structures, (b) An60–70 to An100: anorthite-type structures, (c) An30–40 to An60–70: transitional structures assumed to represent the early stages in the development of the structural complexities found in chemically intermediate low temperature plagioclases.
2. b-Antiphase domains typical for the anorthite-type structure could be imaged in samples ranging between An100 and An67.5, whereas c-antiphase domains have only been found in An90 and An100. The b-domains developed during isothermal crystallization, rather than cooling. The transition from the anorthite-type structure (S.G.I \(\bar 1\) ) to an intermediate C \(\bar 1\) -type structure at An60–70 appears to be a gradual one.
3. The relative contribution of various structural and chemical influences on the cell edges and angles is discussed in detail. It is shown that the nonlinear variation of γ with An content is due to incomplete Al,Si disorder in the high series.

     

分布热源作用下裂隙岩体渗流-传热的拉氏变换-格林函数半解析计算方法

以裂隙岩体高放射性核废物地下处置库性能评估为目标,提出了分布热源作用下单裂隙岩体渗流-传热的简化概念模型、控制微分方程和拉氏变换-格林函数半解析法,为进一步采用半解析法计算分布热源作用下多裂隙岩体的渗流-传热问题奠定了基础。针对单裂隙岩体的渗流-传热问题,建立考虑岩石内热源和二维热传导的控制微分方程,利用拉氏变换域微分方程的基本解建立格林函数积分方程,采用解析法处理其中的奇点,通过数值积分和拉氏数值逆变换求解,计算任意时刻裂隙水和岩石的温度分布。通过算例,与基于岩石一维热传导假定的解析解进行了对比,并计算分析了分布热源作用下单裂隙岩体的渗流-传热特征及其对裂隙开度、岩石热传导系数和热流集度的敏感度。算例表明,（1）就裂隙水温度而言,由于考虑了岩石的二维热传导,拉氏变换-格林函数半解析解小于基于岩石一维热传导假定的解析解;（2）裂隙水温度和岩石温度对裂隙开度和热流集度的敏感度较大,对岩石热传导系数的敏感度较小。     

Hybrid discretization of multi-phase flow in porous media in the presence of viscous,gravitational, and capillary forces

Multi-phase flow in porous media in the presence of viscous, gravitational, and capillary forces is described by advection diffusion equations with nonlinear parameters of relative permeability and capillary pressures. The conventional numerical method employs a fully implicit finite volume formulation. The phase-potential-based upwind direction is commonly used in computing the transport terms between two adjacent cells. The numerical method, however, often experiences non-convergence in a nonlinear iterative solution due to the discontinuity of transmissibilities, especially in transition between co-current and counter-current flows. Recently, Lee et al. (Adv. Wat. Res. 82, 27–38, 2015) proposed a hybrid upwinding method for the two-phase transport equation that comprises viscous and gravitational fluxes. The viscous part is a co-current flow with a one-point upwinding based on the total velocity and the buoyancy part is modeled by a counter-current flow with zero total velocity. The hybrid scheme yields C¹-continuous discretization for the transport equation and improves numerical convergence in the Newton nonlinear solver. Lee and Efendiev (Adv. Wat. Res. 96, 209–224, 2016) extended the hybrid upwind method to three-phase flow in the presence of gravity. In this paper, we present the hybrid-upwind formula in a generalized form that describes two- and three-phase flows with viscous, gravity, and capillary forces. In the derivation of the hybrid scheme for capillarity, we note that there is a strong similarity in mathematical formulation between gravity and capillarity. We thus greatly utilize the previous derivation of the hybrid upwind scheme for gravitational force in deriving that for capillary force. Furthermore, we also discuss some mathematical issues related to heterogeneous capillary domains and propose a simple discretization model by adapting multi-valued capillary pressures at the end points of capillary pressure curves. We demonstrate this new model always admits a consistent solution that is within the discretization error. This new generalized hybrid scheme yields a discretization method that improves numerical stability in reservoir simulation.     

Application of the Fourier method to the numerical solution of moving boundary problem in heat conduction

An algorithm for the solution of a nonlinear problem of phase boundary movement and evolution of temperature distribution due to the perturbation in the basal heat flux has been discussed. The reduction of the problem to a system of nonlinear ordinary differential equations with the help of a Fourier series method leads to a stiff system. This stiffness is taken care of by the use of a modified Euler’s method. Various cases of basal heat flow variation have been considered to show the performance and stability of the technique for such a nonlinear system. The first case of step-wise function is taken to analyse the performance of the technique, and the study has been extended to other general cases of linear increase, periodic variation, and box and triangular function type variations in the heat flux. In the step-wise case the phase boundary attains a constant position rapidly if the supplied heat flux is sufficiently large. The effect of periodicity in the heat flow is clearly depicted in the phase boundary movement, where the phase boundary oscillates about the mean position at large times. The absence of any constant level in the case of linear increase in heat flux is due to a very large value of heat flux. In the cases of box car and triangular heat flux the boundary starts moving downward after the cessation of excess heat flux but does not immediately return to its original preperturbation state, instead approaches it at large times. This technique may be applied to more general cases of heat flow variation.