Development of 3 D Finite Element Method for Non-Aqueous Phase Liquid Transport in Groundwater as Well as Verification

Groundwater contamination previously occurred at a broad range of locations in present-day China. There are thousands of kinds of contaminants which can be divided into soluble and insoluble categories in groundwater. In recent years, the non-aqueous phase liquid (NAPL) pollution that belongs to the multi-phase seepage flow phenomenon has become an increasingly prominent topic due to the challenge brought by groundwater purification and its treatment. Migrating with seepage flow and moving into the potable water sources, these contaminants directly endanger people’s health. Therefore, it is necessary to research how these contaminants not only migrate, but also are then accordingly remedied. First, as an analysis means, an effective numerical method is necessary to be built. A three-dimensional finite element method program for analyzing two-phase flow in porous media, which can be applied to the immiscible contaminant transport problem in subsurface flow has been developed in this paper. The fundamental theory and numerical discretization formulations are elaborated. The numerical difficulty brought about by the distinct non-linearity of the temporal evolution of saturation-dependent variables is overcome by the mixed-form formulation. The effectiveness of simultaneous solution (SS) method and its improvement in efficiency are explained. Finally, two computational examples are given for verifying the correctness and demonstrating the preliminary applicability. In addition, the function of two-phase immiscible flow, especially in Fast Lagrangian Analysis of Continua (FLAC) is used to simulate the same examples and the results are compared to further verify the correctness of the numerical development.


Introduction
Different from the traditional of subsurface contaminant problems dominated by diffusion dissolved in water and convection miscible in fluid transport, in recent years, pollution of non-aqueous phase liquids (NAPLs), which belongs to multi-phase seepage flow phenomenon, is becoming increasingly prominent and complicated in the groundwater purification and treatment.Such liquid pollutants are immiscible in water.Most of them are organic pollutants such as hydrocarbons, petroleum products, chlorinated organic solvents, etc.According to the density of NAPLs compared with water, they are divided into light non-aqueous liquids (LNAPLs, lighter than water) and heavy non-aqueous liquids (DNAPLs, heavier than water).During actual industry production, they are commonly seen and exposed in many places [1][2][3]: Strong organic chlorine solvents, such as trichloroethylene and tetrachloroethylene, are widely applicable to dry cleaning of metal, integrated circuits, and electronic components to remove paints, oil, and fat stains [4].Additionally, a large number of gas stations and oil storage facilities have been newly built in many places in the world since the development of petroleum energy and chemicals play an important role in human industrial production activities.Due to improper disposal or poor management, or accidental breakage of storage containers or pipelines, they leak out and seep into the ground [5,6].As shown in Figure 1, these leaking petroleum products will pass through surface soil to aquifers to the aquitard, and they will sneak into the groundwater and expand laterally on the aquitard top boundary.Additionally, the volatile component will quickly spread far away, and even spillover out of the surface and pollute the air.NAPLs can be chronically bounded in the pores of media, and then accumulated there in large quantities.They are extremely difficult to degrade, finally becoming a persistent source of contamination.It retains in the environment, able to biologically accumulate in ecosystems, has substantial negative impacts on human health and the environment, and even directly affects the price of land.
Processes 2018, 6, x FOR PEER REVIEW 2 of 20 Additionally, a large number of gas stations and oil storage facilities have been newly built in many places in the world since the development of petroleum energy and chemicals play an important role in human industrial production activities.Due to improper disposal or poor management, or accidental breakage of storage containers or pipelines, they leak out and seep into the ground [5,6].As shown in Figure 1, these leaking petroleum products will pass through surface soil to aquifers to the aquitard, and they will sneak into the groundwater and expand laterally on the aquitard top boundary.Additionally, the volatile component will quickly spread far away, and even spillover out of the surface and pollute the air.NAPLs can be chronically bounded in the pores of media, and then accumulated there in large quantities.They are extremely difficult to degrade, finally becoming a persistent source of contamination.It retains in the environment, able to biologically accumulate in ecosystems, has substantial negative impacts on human health and the environment, and even directly affects the price of land.In these cases, it is necessary to analyze the mutual dependence and simultaneous flowing of the non-miscible phase liquid and the ground water in the infiltration field.Compared with the miscible fluid, the difference of NAPL in subsurface water is that the different liquids are separated by the interfaces existing in the pores of the medium.The liquid with higher wettability tends to adhere to pore walls.Due to the interface between the liquids is bent by interfacial tension, it is recessed to the liquid with the higher wettability.To balance this effect, the less humidified liquid has a higher pressure in the pores than the more humidified liquid.This usually not negligible pressure difference is called capillary pressure.When the pore size is kept constant, capillary pressure depends on the percentage of the pore space occupied by different liquid phase contents, that is, saturation.The effective relative permeability of each different liquid also depends on the degree of the phase corresponding saturation.These saturation-dependent parameters exhibit nonlinear characteristics in time evolution, and sometimes even relatively strong nonlinearities, posing a challenge to the rigor of numerical solutions.

Numerical Developing Status on NAPL Contaminant Analysis
The migration of incompatible pollutants in groundwater was only noticed internationally in the late 1980s [7][8][9].In the past 40 years, one scholar after another conducted comparatively systematic research.For most cases of the soil contaminated with NAPLs, generally speaking, the problem of two-phase flow is involved in the water-saturated domain of interest, while that of threephase flow considering condensation and gasification phenomena should be taken into account in In these cases, it is necessary to analyze the mutual dependence and simultaneous flowing of the non-miscible phase liquid and the ground water in the infiltration field.Compared with the miscible fluid, the difference of NAPL in subsurface water is that the different liquids are separated by the interfaces existing in the pores of the medium.The liquid with higher wettability tends to adhere to pore walls.Due to the interface between the liquids is bent by interfacial tension, it is recessed to the liquid with the higher wettability.To balance this effect, the less humidified liquid has a higher pressure in the pores than the more humidified liquid.This usually not negligible pressure difference is called capillary pressure.When the pore size is kept constant, capillary pressure depends on the percentage of the pore space occupied by different liquid phase contents, that is, saturation.The effective relative permeability of each different liquid also depends on the degree of the phase corresponding saturation.These saturation-dependent parameters exhibit nonlinear characteristics in time evolution, and sometimes even relatively strong nonlinearities, posing a challenge to the rigor of numerical solutions.

Numerical Developing Status on NAPL Contaminant Analysis
The migration of incompatible pollutants in groundwater was only noticed internationally in the late 1980s [7][8][9].In the past 40 years, one scholar after another conducted comparatively systematic research.For most cases of the soil contaminated with NAPLs, generally speaking, the problem of two-phase flow is involved in the water-saturated domain of interest, while that of three-phase flow considering condensation and gasification phenomena should be taken into account in the vadose zone.Conceptually, NAPLs are presented in small pores, pore wedges, bypassed pores as films or lenses on water or solid surfaces.They might not be drained from pores in vadose zones after long drainage periods in a strongly water-wet porous medium [10].White, Oostrom, and Lenhard implemented a numerical model on the flow of a nonvolatile nonaqueous phase liquid (NAPL) and aqueous phases that accounts for mobile, entrapped, and residual NAPL in variably-saturated water-wet porous media.The results were also compared against those from detailed laboratory experiments [11].Comparisons between the numerical simulations and experiments demonstrated the necessity to include the residual NAPL formation process in multiphase flow simulators.Xue et al. established a coupling model for analyzing the transport of organic contaminants in soil and water environments [12].Wu and Wang researched the relationship between the riser oil holding rate and the volume fraction with regards to oil measured by the coaxial conductivity sensor, and designed an oil content measuring system based on coaxial conductivity sensor [13].Chen, Yang, and Tian conducted a study on microscopic mechanism of non-aqueous phase liquids (NAPLs) migration in porous media [14].Kleinknecht and Braun undertook not only a large column experiment, but also its numerical simulation to study the density-driven migration of DNAPL gas (carbon disulfide (CS 2 ) vapor) in the vadose zone [15].Kikumoto and Nakamura developed a comprehensive numerical method for simulating transport of non-aqueous phase liquids (NAPLs) in unsaturated subsurface domains [16].Javanbakht et al. conducted a study where X-ray microtomography experiments were performed to investigate the impact of surfactants and microemulsions on the mobilization and resolvability of NAPL in heterogeneous rocks [17].Xie et al. investigated the dynamic behavior, such as flow rate and multi-scale time irreversibility, of different flow patterns based on the measurement signals obtained from oil-gas-water three-phase and oil-water two-phase flow experiments [18,19].Tan et al. studied the flow patterns of horizontal oil-water two-phase pipe flow with water holdup fluctuations provided by a set of conductivity and capacitance sensors [20].Picchi and Battiato, tackled the problem of the limitations of Darcy's law in properly modeling the flow at the continuum scale by proposing a set of upscaled equations based on pore-scale flow regimes, that is, the topology of flowing phases [21].Balasuriya et al. performed a study on detailed field characterization of elemental mercury DNAPL distribution with depth, together with two-phase flow modelling by using STOMP [22].Li et al. developed an energy demodulation algorithm for flow velocity measurement of oil-gas-water three-phase flow [23].Kacem and Benadda built a model to simulate the multiphase extraction (MPE) applied to soil polluted by toluene.The transport and transfer between three phases by using the capillarity equations were simulated [24].These studies involve various aspects of multiphase flows in subsurface water purification.However, in most practical problems, because of the heterogeneity of the considered domain, the irregular shapes of its boundaries, the existence of flow turbulence and phase interfacial interactions, and the liquid mixture in porous medium often exhibit complex behaviors.It is not possible to solve these mathematical models analytically.Instead, the mathematical model is transformed into a numerical one that can be solved by means of computer programs [25].Although there is a large amount of literature on finite difference or finite volume methods for multi-phase flow, the flexible geometry and intrinsic boundary adaptation ability were usually not considered sufficiently.The literature on finite element methods is particularly needed when conducting such issues.However, among these studies and studies, those that discuss the finite element method in this field were rarely seen, especially in the three-dimensional FEM literature, which introduce numerical programs for their counterparts.Descriptions in such detail with operable procedures is hardly found.
This paper introduces the program development of the three-dimensional finite element method for non-miscible and incompressible two-phase flow and preliminary verification with two examples, which belong to LNAPL (oil) and DNAPL (trichloroethylene), respectively.

Fundamental Theory of Two-Phase Flow in Porous Media
In order to use mathematical means to depict the motion disciplines of fluids in porous media, the continuum medium method is introduced [26].It is needed to define the parameters of fluids and porous medium at arbitrary point of domain in this method, such as the velocity of flow, porosity, driving pressure, etc.Therefore, it is hypothesized that a medium is continuously filled over the entire field of research and the actual multiphase seepage porous flow microscopic structure is replaced by average macroscopic meaning, that is, the continuum medium mechanics which at least contains two phases-fluid phase and porous media phase.Both medium and fluid are continuously filled with the entire domain of interest.The parameters of the porous medium, fluid and motion such as porosity, permeability, density, flow rate, concentration, etc., can be defined at any point.This method to describe the migration of fluids and contaminants in porous media is called the continuum medium method, and it has avoided the difficulty of conducting the law of fluid mass point motion in a single pore, studying its parameters using their macroscopic average values instead.These obtained parameters such as flow rate, pressure and concentration are the best approximation of the actual flow, which meets the actual demands.

Governing Equations
In a porous medium, an infinitesimal elemental volume centered at a point of media is taken, which is a volume that can represent the average physical properties near the point.It is called a representative elemental volume (REV, sometimes called the control volume).Similar to the partial differential equation for establishing single-phase seepage flow of porous media, as shown in Figure 2, a cuboid REV is considered.
Processes 2018, 6, x FOR PEER REVIEW 4 of 20 porous medium at arbitrary point of domain in this method, such as the velocity of flow, porosity, driving pressure, etc.Therefore, it is hypothesized that a medium is continuously filled over the entire field of research and the actual multiphase seepage porous flow microscopic structure is replaced by average macroscopic meaning, that is, the continuum medium mechanics which at least contains two phases-fluid phase and porous media phase.Both medium and fluid are continuously filled with the entire domain of interest.The parameters of the porous medium, fluid and motion such as porosity, permeability, density, flow rate, concentration, etc., can be defined at any point.This method to describe the migration of fluids and contaminants in porous media is called the continuum medium method, and it has avoided the difficulty of conducting the law of fluid mass point motion in a single pore, studying its parameters using their macroscopic average values instead.These obtained parameters such as flow rate, pressure and concentration are the best approximation of the actual flow, which meets the actual demands.

Governing Equations
In a porous medium, an infinitesimal elemental volume centered at a point of media is taken, which is a volume that can represent the average physical properties near the point.It is called a representative elemental volume (REV, sometimes called the control volume).Similar to the partial differential equation for establishing single-phase seepage flow of porous media, as shown in Figure 2, a cuboid REV is considered.
In the concept of continuum mechanics, it must be small enough to approximate macroscopic continuity and large enough to be equivalent to microscopic statistical averaging.The flow rate (Q ) flowing into the elemental volume can be expressed as: where ρ is the mass density of the liquid; u, v, and w are the average flow rates through the left, front, and bottom of the element (positive and negative, respectively, represent in and out).
In direction x, if the mass flux flowing into the element is ρu, after the same time interval ∆x, the mass flux flow out of the element can be expressed as Taylor series.
If high-order terms of a small variable are omitted, the corresponding flow (Q ) flowing out of the cuboid per unit time can be expressed as: On the other hand, the cumulative storage amount (Q ) in the cuboid after per unit time can be expressed by the following formula: Mass conservation of flows passing through a small rectangular parallelepiped.
In the concept of continuum mechanics, it must be small enough to approximate macroscopic continuity and large enough to be equivalent to microscopic statistical averaging.
The flow rate (Q i ) flowing into the elemental volume can be expressed as: where ρ is the mass density of the liquid; u, v, and w are the average flow rates through the left, front, and bottom of the element (positive and negative, respectively, represent in and out).
In direction x, if the mass flux flowing into the element is ρu, after the same time interval ∆x, the mass flux flow out of the element can be expressed as Taylor series. (ρu If high-order terms of a small variable are omitted, the corresponding flow (Q t ) flowing out of the cuboid per unit time can be expressed as: On the other hand, the cumulative storage amount (Q a ) in the cuboid after per unit time can be expressed by the following formula: where φ is the porosity of the media, φ = V v /V, V v and V are the volume of the pores and medium, respectively; S is the saturation of the flow phase, S = V f /V v , V f is the flow phase volume.For a flow phase, according to the mass conservation principle 1), (3), and (4) yields: Introducing Einstein's summation convention, Equation (5a) can be written as follow: As the simplified form of conservation of momentum in the flow phase, Darcy's law applies and can be also expressed as: where, k is relative permeability co-efficient of flow phase, µ is dynamic viscosity of flow phase, K ij is absolutely permeability tensor, p is the pressure of flow phase, g is gravitational acceleration, H i is the coordinate in the direction of each coordinate axis, γ is relative density of fluid, and δ ji is the sign of the Kronecker delta.The item ρgH i = 0 only when H i along the gravity direction (assumed direction 3).Thus, Equation (6a) can be also written as: As already mentioned in the above, the flow phases discussed in this paper are all approximated as incompressible, and the pores of the porous media are also approximated as rigid.
Substituting Equation (6b) into Equation (5b), then phase control equations yield as follows: where ∇• is the divergence operator of a vector field; ∇ is the gradient operator of a scalar field; p α , k α , ρ α , S α , and µ α are pressure, relative permeability coefficient, mass density, saturation, and dynamic viscosity of each phase; κ is the intrinsic or absolute permeability.z is the vertical coordinate of the porous medium.In this paper, the intrinsic or absolute permeability and porosity are only functions of position."w" and "o" are representing water and non-aqueous liquids.
As for two-phase flow that contains water and non-aqueous only, there are: and: Noticing Equation (8a), Equation ( 7) can be rewritten as the following equation group: where S is saturation of water, S = S w ; µ r is relative viscosity of water; K is hydraulic conductivity, K = κρ w g/µ w ; ψ α (α = w, o) is pressure expressed in the water head of phase α, ψ α = p α /(ρ w g); ρ r is the relative density, ρ r = ρ o /ρ w .

Finite Differential Discretization in Time
The finite difference discretization is used for time evolution of Equation ( 9) and the following equations are obtained: where the superscript indicates "n" is the known value of the variable at the end of the previous time step, ψ n = ψ(t), "n + 1" indicates the value to be evaluated at the end of the current time step, ψ n+1 = ψ(t + ∆t); and the superscript indicates "v" is the value obtained by the current time step in the previous iteration, "v + 1" indicates the value obtained by the current time step variable at the completion of this iteration; θ is a parameter that indicates the type of the difference method; when θ = 1, the backward difference format is adopted, and when θ = 1/2, the central difference format is adopted.In order to solve the equations, the pressure head of each phase ψ w = ψ w (x, y, z, t) and ψ o = ψ o (x, y, z, t). are selected to be the basic variables, and the following equations are used: Based on these, the saturation of water can be expressed as the function of ψ w or ψ o .In two-phase flow, the void space is completely filled by the two fluids.One of the fluids (the wetting fluid) wets the porous medium more than the other (the non-wetting fluid).As a result, the pressure in the non-wetting fluid will be higher than the pressure in the wetting fluid.The pressure difference is the capillary pressure ψ c , which is a function of saturation.As shown in Figure 3a, the saturation of water phase can be expressed as a function of the capillary pressure head, S = S(ψ c ), where capillary pressure is defined as ψ c = ψ o − ψ w .Water phase saturation can be obtained by the following finite-difference equation: As shown in Figure 3b, the relative permeability coefficient of each phase is a function of water phase saturation, and it can also be further expressed as a function of ψ w and ψ o .For example, permeability coefficient of water phase can be expressed as: Processes The Taylor series of water phase saturation expands to: where: This expansion corresponds to the backward difference format, which has a first-order precision after omitting the high-order terms.If the Taylor series of saturation expands as follows: where: adding the first and the second formula of Equation (14b), the follow equation can be obtained: Equation ( 15) corresponds to central difference format, which has second-order precision after omitting higher-order terms.
If the specific water content function is defined as C(ψ) = ϕ ∂S/ ∂ψ, which represents the change in water content caused by per unit change in capillary pressure head.Then, the evolution of The Taylor series of water phase saturation expands to: where: This expansion corresponds to the backward difference format, which has a first-order precision after omitting the high-order terms.If the Taylor series of saturation expands as follows: where: adding the first and the second formula of Equation (14b), the follow equation can be obtained: Equation ( 15) corresponds to central difference format, which has second-order precision after omitting higher-order terms.
If the specific water content function is defined as C(ψ) = φ∂S/∂ψ, which represents the change in water content caused by per unit change in capillary pressure head.Then, the evolution of saturation can be expressed as: where: It should be pointed out that ∆ t S = S n+1,v+1 − S n+1,v should be able to converge to 0 in the gradual iteration.In theory, Equation ( 16) can be rewritten into: 2.3.Mixed-Form Formulation Adopted in FEM Program Development However, Equation ( 17) does not consistently meet the requirements for the conservation of mass.As shown in Figure 4, because of the nonlinearity of the water content characteristic curve, there exists such a case: where the actual saturation change is large corresponding to a small capillary pressure head increment.Additionally, if Equation ( 17) is substituted into Equation (10) to eliminate the saturation term, and equations with the pressure head as the unknown quantity will be directly formed.These equations whose linear values correctly approximate the saturation require unusually high convergence accuracy, so the efficiency is extremely low or not feasible at all.If this is not recognized and general convergence accuracy is set, the calculation process may be unstable or the results will be totally unacceptable.It should be pointed out that Δ S = S , − S , should be able to converge to 0 in the gradual iteration.In theory, Equation ( 16) can be rewritten into:

Mixed-Form Formulation Adopted in FEM Program Development
However, Equation ( 17) does not consistently meet the requirements for the conservation of mass.As shown in Figure 4, because of the nonlinearity of the water content characteristic curve, there exists such a case: where the actual saturation change is large corresponding to a small capillary pressure head increment.Additionally, if Equation ( 17) is substituted into Equation (10) to eliminate the saturation term, and equations with the pressure head as the unknown quantity will be directly formed.These equations whose linear values correctly approximate the saturation require unusually high convergence accuracy, so the efficiency is extremely low or not feasible at all.If this is not recognized and general convergence accuracy is set, the calculation process may be unstable or the results will be totally unacceptable.The main performance is that the conservation of mass is not guaranteed, and the position of the depth of the infiltration is estimated incorrectly.Therefore, it is necessary to adopt the mixed-form formula advocated by Celia et al. [27].Equation ( 10) is expressed as follows: Incorrect calculated saturation increment t S Correct saturation increment  S  t ψc t ψ t The main performance is that the conservation of mass is not guaranteed, and the position of the depth of the infiltration is estimated incorrectly.Therefore, it is necessary to adopt the mixed-form formula advocated by Celia et al. [27].Equation ( 10) is expressed as follows: Processes 2019, 7, 116 9 of 20 Equation (10) and Equation ( 18) are in hybrid form because some items are based on the pressure head, and some items are based on the saturation.In Equation ( 18), the right side of the equations are known values of the variable at the last time step or the end of the previous iteration.The left items of the equations are expressed as the pressure head that will increase in the current iteration, and the quantity is the unknown value to be sought.

Numerical Discretization Formulations Using the Finite Element Method
In space, the shape of the variable to be sought is discretized as follows: We adopt an eight-node hexahedron iso-parametric element, where: is an interpolation function.P 0 (x ξ , y ξ , z ξ represents the local coordinates of the corner points in the eight-node hexahedral element; ξ is the local grid point number for the corresponding element, ξ = 1, 2, . . ., 8.
According to the standard Galerkin finite element method format, the discrete residual function should be equal to 0 in the global domain, so there are: where e. is element sequence number, e = 1, 2, . . ., N; m is sequence number for the grid point in the whole domain, m = 1, 2, . . ., NP.
Applying partial integral and Gaussian divergence theorem to Equation ( 20), after settlement it yields: where: ∂x j dV e (21d) when m = (ξ), δ m(ξ) = 1, when m = (ξ), δ m(ξ) = 0; (ξ) represents the global node number corresponding to the local node number of the node in the element.n i is normal direction vector The number of equations or variables to be solved in Equation (21a) is double to the total number of nodes because the two-phase flow itself has two degrees of freedom.Noticing that in this form the coefficient matrix will remain symmetrical, an unsteady iterative solution can be used to solve the linear equations formed during each iteration of each time step, such as the conjugate gradient method, which is very efficient, stable, fast, and accurate.As for the iterative steps for solving a system of nonlinear equations, as expressed by Equation ( 18) and ( 21), if an incremental form of the unknown pressure head is used, the general iterative method is equivalent to the Newton-Raphson iterative solution when the incremental form is not used.

Water Flooding Oil in the Pillar Rock Sample (LNAPL)
The first example is to simulate a rock sample water displacing oil experiment to test the numerical method presented in this paper.As shown in Figure 4, the static rock sample is filled with oil, such as the initial saturation condition of the water is zero.The initial value of the capillary pressure of the sample for oil and water is measured in advance, for example, the result is water head pressure 25 cm.The side is sealed and the top boundary is kept as a certain fixed pressure.At the beginning of the experiment, water is injected from the bottom boundary with a certain, but slightly higher, pressure to gradually displace most of the oil from the top.The specific pressure values at the top and bottom boundaries are shown in Figure 5.
The capillary pressure head and the relative permeability coefficient depend on water saturation using the formulas of Van Genucheten [28] and Parker et al. [29][30][31]: where, S we is the effective water saturation, S we = (S w − S wr )/(1 − S w − S wr ), S wr and S or are the residual saturations of water and non-aqueous liquid, respectively; k w and k o are the relative permeability coefficients of water and non-aqueous liquid; α and n are both van Genuchten constants determined by experimental data; and r = 1 − 1/n.In this example, λ = 0.431 cm −1 , n = 3.7.The residual saturation of oil, and water are both 0. The required samples and other physical parameters of the liquids are shown in Table 1.q = q = 0 Boundary condition: Boundary condition: Ψ = 4.5cm, q = 0 Mesh generation into 1 106 The residual saturation of oil, and water are both 0. The required samples and other physical parameters of the liquids are shown in Table 1.Different from the numerical method above, FLAC2D, version 7.0, is calculating with finite-difference method.Noticed that the sample is a quadrangular, and FLAC considers only a plane model (the third direction is processed by unit length).In order to testify the influence of the sample thickness in this example, a unit thickness was set to the quadrangle and it was found that this is a complete one-dimensional problem.Therefore, the difference in thickness can be ignored.In FLAC, the increments of the nodal pore pressure and saturation are expressed as the following equations [28,32]: where, by definition: and ψ w , ψ o , S w , S o are nodal pressure, and saturation of water and non-wetting fluid.K w , K o are fluid bulk moduli, V is the nodal volume, Q is the nodal flow rate, β is the undrained coefficient, which is a constant and equal to one for mechanical coupling and zero for a stand-alone fluid flow calculation.ψ c is the derivative of the capillary curve.Differentiating f equations with respect to S w gives the following Equation (30): Visualizing the calculated date results from self-development by visualization software, and then generating the saturation distribution contours, the water saturation distribution in the rock sample after 5000 s is shown in Figure 6.The more detailed water distribution information is shown in Figure 7.The water content distribution at that time is shown as Figure 8.The specific water capacity refers to a change in the volumetric water content corresponding to the unit change of pressure head.The distribution of the relative permeability coefficients of water and oil after 5000 s is shown in Figure 9.
As shown in Figures 6-9, the results which were calculated by self-developed program are essentially consistent with those correspondingly by FLAC.The discrepancy yet existing between the outcomes from source code and FLAC is due to the difference in numerical formulations between the implicit FEM and explicit finite difference method.Concretely speaking, the latter is an overlaid quadrilateral (four triangular elements for each) 'element' as its feature.Additionally, the striking divergence between the theory described before and mentioned in this chapter will also contribute to this subtle difference.It lies in whether the compressibility of fluids is considered.
then generating the saturation distribution contours, the water saturation distribution in the rock sample after 5000 s is shown in Figure 6.The more detailed water distribution information is shown in Figure 7.The water content distribution at that time is shown as Figure 8.The specific water capacity refers to a change in the volumetric water content corresponding to the unit change of pressure head.The distribution of the relative permeability coefficients of water and oil after 5000 s is shown in Figure 9.As shown in Figures 6-9, the results which were calculated by self-developed program are essentially consistent with those correspondingly by FLAC.The discrepancy yet existing between the outcomes from source code and FLAC is due to the difference in numerical formulations between the implicit FEM and explicit finite difference method.Concretely speaking, the latter is an overlaid quadrilateral (four triangular elements for each) 'element' as its feature.Additionally, the striking divergence between the theory described before and mentioned in this chapter will also contribute to this subtle difference.It lies in whether the compressibility of fluids is considered.As shown in Figures 6-9, the results which were calculated by self-developed program are essentially consistent with those correspondingly by FLAC.The discrepancy yet existing between the outcomes from source code and FLAC is due to the difference in numerical formulations between the implicit FEM and explicit finite difference method.Concretely speaking, the latter is an overlaid quadrilateral (four triangular elements for each) 'element' as its feature.Additionally, the striking divergence between the theory described before and mentioned in this chapter will also contribute to this subtle difference.It lies in whether the compressibility of fluids is considered.Therefore, before explaining the specific reason of this discrepancy, it is necessary to clarify the impact of the difference based on the theory first.To simplify the comparison process, a twodimensional Darcy's flow was elaborated.In a saturated fluid steady horizontal flow in the ground, if considering the compressibility of the fluid, the governing equations are the follows: while ignoring compressibility of fluid: Comparing Equation (31) with Equation ( 32), it can be noticed that, in fact, FLAC's built-in algorithm has considered redundant items ∂ρ/(∂x) and ∂ρ/ ∂y (fluid density ρ can be regarded as the function of moduli K α , while a constant in self-development code), which contributes to the main difference.So as to cut down the deviations to the lowest, the most direct and simple solution is to Therefore, before explaining the specific reason of this discrepancy, it is necessary to clarify the impact of the difference based on the theory first.To simplify the comparison process, a two-dimensional Darcy's flow was elaborated.In a saturated fluid steady horizontal flow in the ground, if considering the compressibility of the fluid, the governing equations are the follows: while ignoring compressibility of fluid: Comparing Equation (31) with Equation ( 32), it can be noticed that, in fact, FLAC's built-in algorithm has considered redundant items ∂ρ/(∂x) and ∂ρ/∂y (fluid density ρ can be regarded as the function of moduli K α , while a constant in self-development code), which contributes to the main difference.So as to cut down the deviations to the lowest, the most direct and simple solution is to increase the value K α .However, in realistic calculations, oversized K α (almost the actual value) will result in a significant increase in the calculation time, as shown in Equation (33), and even terminating the process in order to ensure the convergence for transit seepage flow: where L z is the smallest zone size in the simulation, and k o is the saturated mobility coefficient for the non-wetting fluid (k o = k w µ w /µ o ).Although it is difficult to precisely position the interface between water and oil, it is believed that the required accuracy can still be attained by finer gridding and longer computational time.

Simulation of Trichloroethylene Transport in the Aquifer (DNAPL)
The second example is described below.It is used to verify the correctness of the developed program by simulating a typical contamination transport problem.As shown in Figure 10, there is a non-aqueous liquid (trichloroethylene) pressure head with a 0.1 m source at the midpoint of the aquifer model and a low permeability rock formation in the center of the aqueous model (0.667 m × 0.100 m).The top and bottom of the model are impervious borders, while the left and right sides are fixed hydrostatic pressure boundaries based on initial saturation conditions.In the horizontal direction, the water is driven by a 2% pressure gradient and slowly seeps from right to left.
Although it is difficult to precisely position the interface between water and oil, it is believed that the required accuracy can still be attained by finer gridding and longer computational time.

Simulation of Trichloroethylene Transport in the Aquifer (DNAPL)
The second example is described below.It is used to verify the correctness of the developed program by simulating a typical contamination transport problem.As shown in Figure 10, there is a non-aqueous liquid (trichloroethylene) pressure head with a 0.1 m source at the midpoint of the aquifer model and a low permeability rock formation in the center of the aqueous model (0.667 m × 0.100 m).The top and bottom of the model are impervious borders, while the left and right sides are fixed hydrostatic pressure boundaries based on initial saturation conditions.In the horizontal direction, the water is driven by a 2% pressure gradient and slowly seeps from right to left.
The physical properties of the medium are shown in Table 2.For model aquifers and lowpermeability aquifers, the residual saturation of water is 5%, while the residual saturation of nonaqueous liquids is 0.5.The physical properties of the medium are shown in Table 2.For model aquifers and low-permeability aquifers, the residual saturation of water is 5%, while the residual saturation of non-aqueous liquids is 0.5.The calculation results can be obtained by discretizing the fields and performing numerical simulations using a computer program written in accordance with the finite element formula described above.Meanwhile, the same parameters according to Table 2 in FLAC were set and then this simulation ran.To further verify the correctness of the source program, the numerical simulation result contour figures of the two algorithms after 360 min, 720 min, 1080 min, 1440 min, 1800 min, 2160 min, 2520 min, and 2880 min are compared.As an example, the distribution and variation of trichloroethylene saturation over time in the model domain is shown in Figure 10.

Media
During the infiltration process, under the action of gravity, DNAPL (trichloroethylene) moves in three directions in this model simultaneously and generates an approximately circular pollution area, the closer the place is away from the injection source, the higher the concentration of contaminant.Meanwhile, horizontal flow driving force makes the circular area gradually transform into an elliptical shape, as shown in Figure 11a.
As the infiltration process continues, DNAPL has reached the aquitard and accumulated at the top surface of there, shown in Figure 11b.With the infiltration time increased, DNAPL migration states are shown in the next series of contour figures (Figure 11c-h).
According to these comparison results of the above saturation contours, it can be concluded that either calculated by the self-developed source code or FLAC, the range scale of transport of trichloroethylene in the aquifer after the same time interval is basically regarded as consistent.
These two typical examples have been well simulated using two difference methods: the self-developed finite element method and the advanced finite differential method in FLAC as the existing commercial simulator.The calculation results had both been obtained ideally.They suggest that the comparisons were essentially in consistency, which verifies the correctness of the developed source program, and also show its preliminary practical applicability.The calculation results can be obtained by discretizing the fields and performing numerical simulations using a computer program written in accordance with the finite element formula described above.Meanwhile, the same parameters according to Table 2 in FLAC were set and then this simulation ran.To further verify the correctness of the source program, the numerical simulation result contour figures of the two algorithms after 360 min, 720 min, 1080 min, 1440 min, 1800 min, 2160 min, 2520 min, and 2880 min are compared.As an example, the distribution and variation of trichloroethylene saturation over time in the model domain is shown in Figure 10.
During the infiltration process, under the action of gravity, DNAPL (trichloroethylene) moves in three directions in this model simultaneously and generates an approximately circular pollution area, the closer the place is away from the injection source, the higher the concentration of contaminant.Meanwhile, horizontal flow driving force makes the circular area gradually transform into an elliptical shape, as shown in Figure 11a.

Figure 1 .
Figure 1.Schematic illustration of the subsurface transport mechanism of NAPL contaminants.

Figure 1 .
Figure 1.Schematic illustration of the subsurface transport mechanism of NAPL contaminants.

Figure 2 .
Figure 2. Mass conservation of flows passing through a small rectangular parallelepiped.

Figure 3 .
Figure 3. Nonlinear character of water and NAPL two-phase flow: (a) Two-phase flow characteristic curve of the water-non-aqueous liquid; and (b) the relative permeability coefficient of each phase depends on the change of the water phase saturation.

Figure 3 .
Figure 3. Nonlinear character of water and NAPL two-phase flow: (a) Two-phase flow characteristic curve of the water-non-aqueous liquid; and (b) the relative permeability coefficient of each phase depends on the change of the water phase saturation.
curve of water-NAPL two phase Imbibition Drainage

Figure 4 .
Figure 4. Illustration of inability of using single form formulation to model immiscible flow.

Figure 4 .
Figure 4. Illustration of inability of using single form formulation to model immiscible flow.

Figure 5 .
Figure 5. Displacing test thorough a rock sample: (a) Model size and boundary conditions; and (b) model meshing.

Figure 6 .
Figure 6.Contrast of simulated water saturation distribution after 5000s: (a) Instantaneous water saturation distribution calculated by self-developed source codes; and (b) instantaneous water saturation distribution calculated by FLAC.

Figure 6 .Figure 6 .
Figure 6.Contrast of simulated water saturation distribution after 5000s: (a) Instantaneous water saturation distribution calculated by self-developed source codes; and (b) instantaneous water saturation distribution calculated by FLAC.

Figure 8 .
Figure 8. Simulated specific volumetric water content distribution after 5000 s.

Figure 10 .
Figure 10.Aqueous model test: (a) Model size and boundary conditions; and (b) model meshing.

Figure 10 .
Figure 10.Aqueous model test: (a) Model size and boundary conditions; and (b) model meshing.

Figure 11 .
Figure 11.Distribution of saturation with regards to trichloroethylene and its change with time simulated using self-developed source code (left) and FLAC (right).

Table 1 .
Physical properties of rock sample and fluids (24 °C).

Table 1 .
Physical properties of rock sample and fluids (24 • C).

µ
Dynamic viscosity of flow phase µ w Dynamic viscosity of water µ o Dynamic viscosity of NAPL µ r Relative viscosity of NAPL to water (µ r = µ o /µ w ) K Hydraulic conductivity ψ w Pressure in water head of water (ψ w = p w /(ρ w g)) ψ o Pressure in water head of NAPL ( ψ o = p o /(ρ w g)) ψ c Capillary pressure in water head (ψ c = ψ o − ψ w )