Comput Geosci (2011) 15:167–183 DOI 10.1007/s10596-010-9206-2
ORIGINAL PAPER
Computational simulation for the morphological evolution of nonaqueous phase liquid dissolution fronts in two-dimensional fluid-saturated porous media Chongbin Zhao · B. E. Hobbs · K. Regenauer-Lieb · A. Ord
Received: 18 October 2009 / Accepted: 9 August 2010 / Published online: 24 August 2010 © Springer Science+Business Media B.V. 2010
Abstract This paper deals with the computational aspects of nonaqueous phase liquid (NAPL) dissolution front instability in two-dimensional fluid-saturated porous media of finite domains. After the governing equations of an NAPL dissolution system are briefly described, a combination of the finite element and finite difference methods is proposed to solve these equations. In the proposed numerical procedure, the finite difference method is used to discretize time, while the finite element method is used to discretize space. Two benchmark problems, for which either analytical results or previous solutions are available, are used to verify the proposed numerical procedure. The related simulation results from these two benchmark problems have demonstrated that the proposed numerical procedure is useful and applicable for simulating the morphological evolution of NAPL dissolution fronts in two-dimensional fluid-saturated porous media of finite domains. As an application, the proposed numerical procedure has been used to simulate morphological evolution processes for three kinds of NAPL dissolu-
C. Zhao (B) Computational Geosciences Research Centre, Central South University, Changsha 410083, China e-mail:
[email protected] B. E. Hobbs · K. Regenauer-Lieb · A. Ord School of Earth and Environment, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia K. Regenauer-Lieb CSIRO Division of Earth Science and Resource Engineering, P. O. Box 1130, Bentley, WA 6102, Australia
tion fronts in supercritical NAPL dissolution systems. It has been recognized that: (1) if the Zhao number of an NAPL dissolution system is in the lower range of the supercritical Zhao numbers, the fundamental mode is predominant; (2) if the Zhao number is in the middle range of the supercritical Zhao numbers, the (normal) fingering mode is the predominant pattern of the NAPL dissolution front; and (3) if the Zhao number is in the higher range of the supercritical Zhao numbers, the fractal mode is predominant for the NAPL dissolution front. Keywords Nonaqueous phase liquid · Residual saturation · Computational simulation · Fingering mode · Fractal mode · Porous media
1 Introduction Effective and efficient removal of nonaqueous phase liquids (NAPLs) from contaminated sites has attracted considerable attention during the past two decades [7, 8, 10–14, 18, 23, 24, 27]. Due to their low solubility and relatively high toxicity, NAPLs consisting of immiscible fluids and organic species cannot be adequately removed from the subsurface by means of traditional pump-and-treat methods [18, 24]. To develop effective and efficient remediation techniques for NAPL-contaminated sites, a considerable amount of laboratory experiments have been conducted to understand the potential physical and chemical processes controlling the movement of the NAPL in fluidsaturated porous media. Major achievements from these laboratory experiments are: (1) the quantitative determination of mass transfer rates between an NAPL
168
and an aqueous phase liquid and (2) the recognition of fingering phenomena of NAPL dissolution fronts at the laboratory (i.e., centimeter) scale. With these achievements, both mathematical and computational models [9, 14] have been developed to simulate the morphological evolution of NAPL dissolution fronts in fluid-saturated porous media. However, due to huge requirement for computer efforts, the previous computational simulation of NAPL dissolution front evolution problems was limited to the laboratory scale [11, 14]. For the numerical simulation of reactive mass transport problems where the dissolution front instability may be unimportant, a significant amount of research has been conducted to simulate large length-scale geological and geochemical systems [1, 16, 19, 21, 22, 25, 26, 28–34]. However, when the dissolution front instability is considered, the numerical simulation is mainly focused on the generic and small length-scale problems. In this aspect, much work has been conducted on mineral dissolution front instability problems in fluidsaturated porous media [3–6, 17, 20, 35, 36, 38, 40]. Although there are some similarities between mineral dissolution front instability problems and NAPL dissolution front instability ones, the following considerable differences exist between them. First, in the case of mineral dissolution induced instability the dissolution process in a fluid-saturated porous medium is dominated by the chemical reaction process (i.e., the diagenetic process or water–rock chemical reaction process), while in the case of NAPL dissolution induced instability it is dominated by the physical elution process [14]. Second, the mineral dissolution ratio, which is defined as the ratio of the equilibrium concentration of the dissolvable mineral to the molar density of the solid mineral, is in a range of 10−3 to 10−9 for typical rocks [3], while the NAPL dissolution ratio, which is similarly defined as the ratio of the equilibrium concentration of an NAPL in the aqueous phase fluid to the mass density of the NAPL itself, is of the order of 10−3 [9]. This means that the time scale of a mineral dissolution system is much larger than that of an NAPL dissolution system. For this reason, it is possible to conduct laboratory experiments to observe NAPL dissolution front fingering phenomena [10], but it is very difficult, if not impossible, to observe the mineral dissolution front instability phenomenon through laboratory experiments. Third, the pore-fluid velocity in a mineral dissolution system is often a few centimeters per year; while the aqueous phase fluid velocity in an NAPL dissolution system may reach a few meters per day. This indicates that solute dispersion (i.e., mechanical dispersion) may be
Comput Geosci (2011) 15:167–183
neglected for a mineral dissolution system, but it must be considered in an NAPL dissolution system. Due to these differences, the computational methods developed for simulating mineral dissolution front instability problems cannot be directly used to simulate NAPL dissolution front instability problems without necessary modifications. From previous laboratory experiments and numerical simulations, only fractal fingering modes have been observed for NAPL dissolution fronts in fluid-saturated porous media. However, from the theoretical point of view, other kinds of modes, such as the fundamental mode and the (normal) fingering mode, should exist in some supercritical NAPL dissolution systems. To test this hypothesis, a numerical procedure that is a combination of the finite difference and finite element methods is proposed in this paper. After the proposed numerical procedure is verified by two benchmark problems, for which either analytical results or previous solutions are available, it is used to investigate the likelihood of occurrence for three kinds of NAPL dissolution fronts in supercritical NAPL dissolution systems. The related simulation results have demonstrated that depending upon the relative position where the Zhao number of an NAPL dissolution system is located in the range of the supercritical Zhao numbers, three kinds of modes, namely the fundamental mode, the (normal) fingering mode, and the fractal mode, can take place for NAPL dissolution fronts in two-dimensional fluidsaturated porous media. If the Zhao number of an NAPL dissolution system is in the lower range of the supercritical Zhao numbers, the fundamental mode is predominant, while if the Zhao number is in the middle range of the supercritical Zhao numbers, the (normal) fingering mode is the predominant pattern of the NAPL dissolution front. However, if the Zhao number is in the higher range of the supercritical Zhao numbers, the fractal mode is predominant for the NAPL dissolution front. Compared with a previous publication [39], the major contribution of this paper can be described as follows. In the previous paper, a detailed theoretical analysis was carried out to investigate NAPL dissolution front instability in two-dimensional fluid-saturated porous media of infinite domains. In particular, the Zhao number was proposed to represent the fundamental characteristics of an NAPL dissolution system, so that a stability criterion was established to assess the instability likelihood of an NAPL dissolution front. Although the fractal fingering mode was considered to demonstrate the usefulness and applicability of the established stability criterion, the consideration of the detailed characteristics of different mode patterns was
Comput Geosci (2011) 15:167–183
169
lacking in the previous paper [39]. The major contribution of the current paper is to apply the Zhao number and related stability criterion to better characterize NAPL dissolution fingering modes in supercritical Zhao number systems. As a result, three kinds of modes, namely the fundamental mode, the (normal) fingering mode, and the fractal mode, have been found in NAPL dissolution systems of supercritical Zhao numbers.
2 Governing equations of the NAPL dissolution problems in two-dimensional fluid-saturated porous media Based on the previous studies [9], the governing equations of an NAPL dissolution problem in a twodimensional fluid-saturated porous medium can be expressed as follows: φ φ
∂ Sn K Ceq − C =− ∂t ρn
∂ Sn φ = −∇ • ∂t
k (Sn ) K ∇ pa − Ceq − C μa ρa
Dh =
va =
(2)
(3)
where φ is the porosity; Sn is the NAPL saturation (i.e., the fraction of the void space occupied by the NAPL); ρ a and ρ n are the aqueous phase and nonaqueous phase densities; K is the mass transfer rate coefficient to express the exchange of the NAPL species from the nonaqueous phase to the aqueous phase; C is the solute concentration of the NAPL species in the bulk aqueous phase; Ceq is the equilibrium concentration of the NAPL species in the aqueous phase; Dh is the general dispersion tensor; k(Sn ) is the saturation-dependent permeability of the porous medium to aqueous phase flow; μa is the dynamic viscosity of the aqueous phase; and pa is the aqueous phase pressure. Note that Sn , C, and pa are three independent variables in these three equations so that the NAPL dissolution problem is well posed mathematically. Although most NAPLs found in natural porous media that arise from contaminant spills are mixtures, only a singlephase NAPL is assumed in this investigation.
τ Dm + α L va 0 0 τ Dm + αT va
(5)
where kf is the intrinsic permeability of the porous medium after the NAPL is completely dissolved; Sai is the irreducible saturation of the aqueous phase; τ is the tortuosity of the porous medium; Dm is the molecular diffusivity of the NAPL species in the aqueous phase; α T and α L are the transversal and longitudinal dispersivities of the NAPL species in the aqueous phase; va is the amplitude of the average linear velocity vector (va ) of the aqueous phase. According to Darcy’s law, the average linear velocity vector (va ) can be expressed as follows [2]:
(1)
∂ (1 − Sn ) C = ∇ • φ (1 − Sn ) Dh • ∇C ∂t k (Sn ) C∇ pa +∇ • μa + K Ceq − C
The following expressions can be used for the saturation-dependent permeability and the dispersion tensor, respectively: 1 − Sn − Sai 3 k (Sn ) = kf (4) 1 − Sai
−k (Sn ) ∇ pa φ (1 − Sn ) μa
(6)
where k(Sn ) is the saturation-dependent permeability of the porous medium to aqueous phase flow; μa is the dynamic viscosity of the aqueous phase; and pa is the aqueous phase pressure. The previous experimental results [9] indicate that the mass transfer rate coefficient (K) can be expressed in the following form: K = β0 Sβn1
(7)
where β 0 is a function of the porous medium, the NAPL, and the velocity, viscosity, and density of the aqueous phase fluid; β 1 is a constant. If the computational domain of an NAPL dissolution system is finite, then the boundary conditions of the problem can be expressed as follows: C = 0, C = Ceq ,
Sn = 0,
∂ pa = paxf ∂x
Sn = Sn0 ,
pa = pa0
(x = 0) (x = Lx )
(8) (9)
where Sn0 is the initial saturation of the NAPL; paxf is the pressure gradient of the aqueous phase on the upstream boundary; pa0 is the pressure of the aqueous phase on the downstream boundary; and Lx is the length of the problem domain in the x direction. Since paxf drives the aqueous phase fluid flow continuously along the positive x direction, it has a negative algebraic value (i.e., paxf < 0) for the problem under consideration.
170
Comput Geosci (2011) 15:167–183
Except for the upstream boundary, the initial conditions of the problem for the rest of the computational domain are as follows: C = Ceq , Sn = Sn0 0 < x ≤ Lx , 0 ≤ y ≤ L y (10) where L y is the length of the problem domain in the y direction. Since the use of dimensionless governing equations has some advantages in dealing with problems of multiscales and multi-processes [40], it is useful to transform the above-mentioned governing equations of the problem into the corresponding dimensionless form. For this purpose, the following dimensionless quantities have been defined: x=
x , L∗
C=
C , Ceq
pa =
pa , p∗a
αL , L∗
αT =
αT , L∗
αL =
where φ t∗ = , β0 Dh =
y=
y , L∗
Lx =
Lx L∗
t=
t ε t∗
ε=
L∗ = τ Dm t∗ ,
Dh , τ Dm
ψ (Sn ) =
Ly =
Ly L∗
(11)
(12)
Ceq << 1 ρn
(13)
φτ Dm p∗a = ψ (Snf )
ψ (Sn ) , ψ (Snf )
ψ (Sn ) =
(14) k (Sn ) μa (15)
where ψ(Snf ) is the value of ψ(Sn ) at Sn = Snf ; and Snf = 0 is the saturation of the NAPL after it is completely dissolved in the NAPL dissolution system. Inserting Eqs. 11, 12, 13, 14, and 15 into Eqs. 1, 2, and 3 yields the following dimensionless governing equations for the NAPL dissolution system: ∂ Sn + Sβn1 1 − C = 0 ∂t
(16)
∂ ε (1 − Sn ) C ∂t − ∇ • (1 − Sn ) Dh • ∇C + Cψ (Sn ) ∇ pa +
ε
∂ Sn =0 ∂t
∂ Sn ε β1 + ∇ • ψ (Sn ) ∇ pa + Sn 1 − C = 0 ρa ∂t
(17)
(18)
where ⎡
⎤ ψ (Sn ) ∂ pa 0 ⎢ 1 − αL 1 − S ∂ x ⎥ ⎥ n Dh = ⎢ ⎣ ψ (Sn ) ∂ pa ⎦ 0 1 − αT 1 − Sn ∂ x
(19)
Note that off-diagonal terms in the dispersion tensor are assumed to be zero. This is true when the longitudinal dispersivity is equal to the transversal dispersivity in the fluid-saturated porous medium. It is a good approximation when the perturbed aqueous phase (secondary) flow in the vertical direction perpendicular to the injected aqueous phase (primary) flow direction is much slower than the injected aqueous phase (primary) flow in the horizontal direction within a rectangular domain [6, 39]. Nevertheless, for NAPL dissolution problems involving irregular domains, the injected aqueous phase (primary) flow may have the same order of magnitude in both the horizontal and vertical directions, so that off-diagonal terms in the dispersion tensor need to be considered. The boundary conditions in this case can be expressed in the dimensionless form as follows: C = 0,
Sn = 0,
C = 1,
Sn = Sn0 ,
∂ pa = paxf ∂x pa = pa0
(x = 0)
x = Lx
(20) (21)
where p axf is the dimensionless pressure gradient of the aqueous phase on the upstream boundary; pa0 is the dimensionless pressure of the aqueous phase on the downstream boundary. Similarly, the initial conditions of the problem can be expressed in the following dimensionless form: 0 < x ≤ Lx , 0 ≤ y ≤ L y (22) C = 1, Sn = Sn0
3 Computational simulation of NAPL dissolution front instability problems in two-dimensional fluid-saturated porous media The main purpose of this section is to propose a numerical procedure for simulating how an NAPL dissolution front evolves from a planar shape to a complicated morphology. To verify the accuracy of the computational simulation, a benchmark problem is constructed from the theoretical analysis derived under subcritical Zhao numbers in the previous study [39]. As a result, the numerical solution obtained from the benchmark problem can be compared with the corresponding analytical solution. To further verify whether or not the proposed numerical procedure can be used to simulate the
Comput Geosci (2011) 15:167–183
171
morphological evolution of NAPL dissolution fronts in fluid-saturated porous media of supercritical Zhao numbers, a computational model that was presented by Imhoff et al. [11] on the basis of their previous laboratory experiments [10] is used in this investigation. Since computational simulation results for the complicated morphological evolution patterns of NAPL dissolution fronts compared well with those observed from the corresponding laboratory experiments [10, 14], this computational model can be used as a typical representative of the related laboratory experiments. On the other hand, because both the detailed input parameters and output results of this computational model are openly available [11], it can be served as a useful benchmark problem for verifying numerical procedures that are used to simulate the morphological evolution of NAPL dissolution fronts in fluid-saturated porous media of supercritical Zhao numbers. For this reason, the proposed numerical procedure will be verified by comparing the current results with those previously obtained from the same laboratory experiment-based computational model [11].
The dimensionless governing equations (i.e., Eqs. 16, 17, and 18) of an NAPL dissolution problem in a two-dimensional fluid-saturated porous medium can be solved using a combination of the finite element and finite difference methods. Toward this end, the finite element method is used to discretize the geometrical shape of the problem domain, while the finite difference method is used to discretize the dimensionless time. Since the dimensionless governing equations of an NAPL dissolution system are highly nonlinear, the segregated algorithm, in which Eqs. 16, 17, and 18 are solved separately in a sequential manner, is used to derive the formulation of the proposed numerical procedure. For a given dimensionless time-step, t + t, the saturation of the NAPL can be denoted by (Sn )t + t = (Sn )t + ( Sn )t + t , where (Sn )t is the saturation of the NAPL at the previous time-step and ( Sn )t + t is the saturation increment of the NAPL at the current timestep. Since Eq. 16 is a nonlinear equation, it can be linearized using the following Taylor expansion: β β1 = (Sn )t + ( Sn )t + t 1 (Sn )t + t β
β −1
( Sn )t + t
where Ct+ t is the dimensionless concentration of the NAPL at the current time-step; t is the dimensionless time increment at the current time-step. In the finite difference sense, the following relationships exist:
∂ ε (1 − Sn ) C = ε ∂t
1 − (Sn )t+ t Ct+ t
t
= − εCt+ t
(Sn )t+ t
t
Ct+ t + ε 1 − (Sn )t+ t
t ∂ Sn
(Sn )t+ t β1 = = − 1 − Ct+ t (Sn )t+ t ∂t
t
3.1 Formulation of the proposed numerical procedure for simulating the evolution of NAPL dissolution fronts
= (Sn )t 1 + β1 (Sn )t 1
Using the backward difference scheme, Eq. 16 can be written as follows. 1 1−β1 + β1 1 − Ct+ t ( Sn )t+ t (Sn )t
t = − (Sn )t 1 − Ct+ t (24)
(23)
(25)
(26)
∇ • (1 − Sn ) Dh • ∇C = ∇•
1 − (Sn )t+ t
Dh
t+ t
• ∇Ct+ t
(27)
∇ • Cψ (Sn ) ∇ pa
= C∇ • ψ (Sn ) ∇ pa + ∇ pa • ψ (Sn ) ∇C
= Ct+ t ∇ • ψ (Sn )t+ t ∇ pa t+ t
+ ∇ pa t+ t • ψ (Sn )t+ t ∇Ct+ t
(28)
Considering Eqs. 25 and 26 simultaneously yields the following equation: ε
∂S (S ) ∂ n n t+ t = 1 − εCt+ t (1 − Sn ) C + ∂t ∂t
t + ε 1 − (Sn )t+ t ×
Ct+ t
t
(29)
172
Comput Geosci (2011) 15:167–183
Since ε is much smaller than one, Eq. 29 can be approximately expressed as follows: ε
∂S ∂ n (1 − Sn ) C + ∂t ∂t Ct+ t
(Sn )t+ t + ε 1 − (Sn )t+ t
t
t β1 = − 1 − Ct+ t (Sn )t+ t =
Ct+ t + ε 1 − (Sn )t+ t
t
t+ t
− ∇ pa t+ t • ψ (Sn )t+ t ∇Ct+ t = (30)
Substituting Eqs. 27, 28, and 30 into Eq. 17 yields the following finite difference equation: ε β1 Ct+ t 1 − (Sn )t+ t + (Sn )t+ t
t
− ∇ • 1 − (Sn )t+ t Dh • ∇Ct+ t t+ t
− ∇ pa t+ t • ψ (Sn )t+ t ∇Ct+ t
ε β1 1 − (Sn )t+ t Ct + (Sn )t+ t
t
(31)
Similarly, Eq. 18 can be rewritten in the following discretized form: ∇ • ψ (Sn ) ∇ pa
(32)
Substituting Eq. 32 into Eq. 31 yields the following equation: ε β1 1 − (Sn )t+ t + (Sn )t+ t Ct+ t
t
− ∇ • 1 − (Sn )t+ t Dh • ∇Ct+ t t+ t
− ∇ pa t+ t • ψ (Sn )t+ t ∇Ct+ t − εCt+ t =
1−
1 ρa
β1 1 − Ct+ t (Sn )t+ t
ε β1 1 − (Sn )t+ t Ct + (Sn )t+ t
t
(34)
In summary, the dimensionless governing equations (i.e., Eqs. 16, 17, and 18) of an NAPL dissolution system under the condition of ε << 1 can be expressed in the finite difference form as follows: 1 1−β1 + β1 1 − Ct+ t ( Sn )t+ t (Sn )t
t = − (Sn )t 1 − Ct+ t
(35)
t+ t
− ∇ pa t+ t • ψ (Sn )t+ t ∇Ct+ t =
ε β1 1 − (Sn )t+ t Ct + (Sn )t+ t
t
(36)
∇ • ψ (Sn ) ∇ pa = ∇ • ψ (Sn )t+ t ∇ pa t+ t
= ∇ • ψ (Sn )t+ t ∇ pa t+ t 1 β1 =ε 1− 1 − Ct+ t (Sn )t+ t ρa
ε β1 1 − (Sn )t+ t Ct + (Sn )t+ t
t
ε β1 Ct+ t 1 − (Sn )t+ t + (Sn )t+ t
t
• ∇Ct+ t − ∇ • 1 − (Sn )t+ t Dh
− Ct+ t ∇ • ψ (Sn )t+ t ∇ pa t+ t =
Since ε << 1, Eq. 33 can be approximately expressed in the following form: ε β1 1 − (Sn )t+ t + (Sn )t+ t Ct+ t
t
• ∇Ct+ t − ∇ • 1 − (Sn )t+ t Dh
(33)
1 β1 = ε 1− 1 − Ct+ t (37) (Sn )t+ t ρa Note that Eqs. 16, 17, and 18 are valid for any values of ε, so that if ε is close to unity, they should be used to obtain numerical solutions instead of Eqs. 35, 36, and 37. However, the direct use of Eqs. 16, 17, and 18 may result in more computer efforts, compared with the use of Eqs. 35, 36, and 37. Since most NAPL dissolution problems may fall in the range of ε << 1 and analytical solutions are only available in this range, Eqs. 35, 36, and 37 are used to obtain numerical solutions in this investigation. Using both the proposed segregated scheme and the finite element method [39, 41], Eqs. 35, 36, and 37 are solved separately and sequentially for the saturation of the NAPL, the dimensionless concentration of the NAPL, and the dimensionless pressure of the aqueous
Comput Geosci (2011) 15:167–183
173
phase fluid at the current time-step. To consider the coupling effect between these three equations, an iteration scheme is also used in the process of finding the numerical solution. The following convergence criterion needs to be satisfied so that a convergent solution can be obtained. ⎛ N Sn 2 k k−1 E = max ⎝ − (Sn )i,t+ t , (Sn )i,t+ t i=1
∂C ∂y
Sn = 0 C = C (t)
Inflow
∂pa = p'axf ∂x
y
∂pa =0 ∂y
Sn = Sn0 C =1
Ly
x 0
NC 2 k k−1 Ci,t+ t − Ci,t+ t ,
= 0,
∂C ∂y
= 0,
∂pa =0 ∂y
Lx
i=1
N pa k p
a i,t+ t
⎞
k−1 2 ⎟ − pa i,t+ t ⎠ < E
Fig. 1 Geometry and boundary conditions for the first benchmark problem
(38)
i=1
where E and E are the maximum error at the kth iteration step and the allowable error limit; NSn , NC , and N pa are the total numbers of the degree-of-freedom for the saturation of the NAPL, the dimensionless concentration of the NAPL in the aqueous phase fluid, and the dimensionless pressure of the aqueous phase fluid, respectively; k is the index number at the current iteration step and k−1 is the index number at the previ k k k ous iteration step; (Sn )i,τ + τ , Ci,τ + τ , and pa i,τ + τ are the saturation of the NAPL, the dimensionless concentration of the NAPL, and the dimensionless pressure of the aqueous phase fluid for node i at both the curk−1 rent time-step and the current iteration step; (Sn )i,τ + τ , k−1 k−1 Ci,τ + τ , and pa i,τ + τ are the saturation of the NAPL, the dimensionless concentration of the NAPL, and the dimensionless pressure of the aqueous phase fluid for node i at the current time-step but at the previous iteration step. 3.2 Verification of the proposed numerical procedure for simulating the evolution of NAPL dissolution fronts Due to inevitable round-off errors in computation and discretization errors in temporal and spatial variables, it is necessary to verify the proposed numerical procedure so that meaningful numerical results can be obtained from a discretized computational model. For this reason, two benchmark problems, for which either the analytical results or previous solutions are available, are considered in this subsection. Figure 1 shows the geometry and boundary conditions of the first benchmark problem, in which the dimensionless-pressure gradient (i.e. paxf = −0.02) is
applied to the left boundary. This means that there is a horizontal throughflow from the left to the right of the computational model. The dimensionless height and width of the computational model are 200 and 400, respectively. Except for the left boundary, the initial residual saturation of the NAPL is 0.2, while the initial dimensionless concentration of the NAPL in the aqueous phase fluid is unity within the computational domain. The final residual saturation of a value of zero is applied to the left boundary as a boundary condition of the computational domain. Both the top and the bottom boundaries are assumed to be impermeable for the NAPL and aqueous phase fluid. The following parameters are used in the corresponding computations: the initial saturation (Sn0 ) of the NAPL is 0.2; the irreducible saturation (Sai ) of the aqueous phase fluid is 0.15; the dimensionless longitudinal dispersivity (α L ) is 0.2; the dimensionless transverse dispersivity (α T ) is 0.02; the ratio (ε) of the equilibrium concentration of the NAPL species in the aqueous phase fluid to the density of the NAPL itself is 0.001; the density ratio ρ a of the aqueous phase fluid to the NAPL is 1.0/1.46; the value of β 1 is 0.87; and the Zhao number (i.e., − paxf = 0.02) applied to the left boundary of the NAPL dissolution system is 0.02. Since the critical Zhao number of the system is approximately equal to 0.052, the NAPL dissolution system considered here is in a subcritical state. In this case, the corresponding dimensionless propagation speed of the planar NAPL dissolution front is 0.1. To appropriately simulate the propagation of the NAPL dissolution front, the whole computational domain is simulated by 19,701 four-node square elements of 20,000 nodal points in total. The initial residual saturation field of the NAPL is randomly perturbed by a small amount of 1% of the originally
174
Comput Geosci (2011) 15:167–183
input saturation of the NAPL (i.e., Sn0 = 0.2) before running the computational model. This means that the resulting initial residual saturation is of a random distribution, which has a mean value of the homogeneous residual saturation (i.e., Sn0 = 0.2) and a variation of 0.002 (i.e., 1% of Sn0 = 0.2) in the whole computational domain. Since the computational domain of the benchmark problem is of finite size, a time-dependent dimensionless-concentration boundary condition (i.e., C t = exp paxf v front t ) needs to be applied to the left boundary so that the numerical solutions can be compared with the analytical solutions. From the previous study [39] the analytical solutions for this benchmark problem can be expressed as follows: C x, t = 1, pa x, t = − pax0 Lx − x + 100,
(39)
C x, t = exp − paxf x − v front t , pa x, t = paxf x−v front t − pax0 Lx −v front t +100, Sn (x, τ ) = 0
x ≤ v front t
(40)
Note that the dimensionless pressure pa0 of the aqueous phase fluid on the downstream boundary is assumed to be 100 in the benchmark problem. Figure 2 shows the comparison of numerical solutions with analytical ones for the residual saturation of
0.2
Saturation
0.2 (t = 1000)
(t = 2000) 0.0
0.0 0
0
400
x
x
400
(A) Residual Saturation Dimensionless concentration
Fig. 2 Comparison of the numerical solutions with analytical results at different time instants
Saturation
x > v front t
1
(t = 1000) 0 0
x
400
Dimensionless concentration
Sn (x, τ ) = Sn0
the NAPL and the dimensionless concentration of the NAPL in the aqueous phase fluid at two different time instants. In this figure, thick lines show the numerical results, while thin lines show the corresponding analytical solutions. It can be observed that the overall numerical solutions agree well with the analytical solutions, except for a smooth effect on the numerically simulated propagation front of the NAPL saturation. This fact demonstrates that the proposed numerical procedure is capable of simulating the planar dissolution-front propagation within the fluid-saturated porous medium of a finite domain under subcritical Zhao number conditions. Since a finite value for the ratio (ε) of the equilibrium concentration of the NAPL species in the aqueous phase fluid to the density of the NAPL itself is used in the computational model, there is a discrepancy between the computational model used to produce the numerical results and the theoretical model used to derive the analytical solutions. Because of this discrepancy, both the sharpness and the propagation speed of the NAPL dissolution front obtained from the computational model differ slightly from those obtained from the theoretical analysis. To reduce this discrepancy, the value for the ratio (ε) of the equilibrium concentration of the NAPL species in the aqueous phase fluid to the density of the NAPL itself should be used as small as possible in the computational model. As mentioned previously, the second benchmark problem under consideration is a computational model (as shown in Fig. 3) that is based on the laboratory experiments conducted by Imhoff et al. [10]. Miller et al. [14] have successfully used this computational model to reproduce the same fingering phenomena as
1
(t = 2000) 0
0
(B) Dimensionless Concentration
x
400
Comput Geosci (2011) 15:167–183
175
∂pa = p'axf ∂x Inflow
Using the above parameters, the Reynolds number of the NAPL dissolution system is as follows:
Sn = 0, C = 0, 0 ∂C = 0 ∂y ∂pa =0 ∂y
y x
Sn = 0.14 ± 0.03, C = 1
∂C = 0 ∂y ∂pa =0 ∂y
Re = Lx
=
Ly
Fig. 3 Geometry and boundary conditions for the second benchmark problem
those observed in the corresponding laboratory experiments. The two primary purposes of considering this benchmark problem are: (1) to verify the proposed numerical procedure when it is used to simulate the morphological evolution of NAPL dissolution fronts in supercritical systems and (2) to demonstrate how to use the Zhao number and the critical Zhao number to assess the instability likelihood of an NAPL dissolution system. In order to establish this computational model, the data associated with Experiment 13 in the experimental work of Imhoff et al. [10] was used so that this computational model has a solid foundation of experiment results. The NAPL considered in the experiment is trichloroethylene (i.e., TCE). For the purpose of reproducing this experiment numerically, Imhoff et al. [11] used the following parameters in their numerical experiments: the aqueous solubility (Ceq ) of the NAPL is 1.27 kg/m3 ; the irreducible saturation (Sai ) of the aqueous phase (i.e., water) is 0.104; the molecular diffusion coefficient (Dm ) of the NAPL is 8.4 × 10−10 m2 /s; the porosity (φ) of the porous medium is 0.355; the permeability (kf ) of the porous medium in the isotropic case is 1.0 × 10−12 m2 ; the initial residual saturation (Sn0 ) of the NAPL is 0.14; the dynamic viscosity (μa ) of the aqueous phase (i.e., water) is 1.0 × 10−3 (N s/m2 ); the density (ρ a ) of the aqueous phase (i.e., water) is 1,000 kg/m3 ; the density (ρ n ) of the NAPL is 1,460 kg/m3 ; the ratio (ε) of the aqueous solubility (Ceq ) to the density (ρ n ) of the NAPL is 0.87 × 10−3 ; the Darcy velocity (Vaxf ) of the injected aqueous phase fluid (i.e., water) is 1.0 / 86,400 m/s = 1.16 × 10−5 m/s; the tortuosity (τ ) of the porous medium is 0.66; the longitudinal and transverse dispersivities (α L and α T ) of the porous medium are 7.2 × 10−4 m and 1.44 × 10−4 m, respectively; the geometrical mean particle diameter (d50 ) is 3.6 × 10−4 m; the length of the model in the injected flow direction is 0.30 m; and the length in the direction perpendicular to the injected flow is 1.0 m.
d50 ρa Vaxf φμa (1− Sn0 ) 3.6 × 10−4 × 1000 × 1.16 × 10−5 0.355 × 10−3 × (1 − 0.14)
= 1.368 × 10−2
(41)
When the correlation of Imhoff et al. [8] is used, both β 1 and β 0 involved in the mass transfer rate coefficient (K) can be evaluated using the following equations: β1 = 0.87 β0 = =
(42)
186φ 0.87 Re0.71 Dm d250 186 × 0.3550.87 × 0.013680.71 × 8.4 × 10−10 2 3.6 × 10−4
= 2.325 × 10−2 (1/s)
(43)
Thus the intrinsic characteristic time, length, and aqueous pressure of the NAPL dissolution system can be determined as follows: φ 0.355 = = 15.27 s β0 2.325 × 10−2
L∗ = τ Dm t∗
t∗ =
=
0.66 × 8.4 × 10−10 × 15.27
= 9.20 × 10−5 m p∗a =
(44)
(45)
φτ Dm φτ Dm μa = ψ (Snf ) kf
0.355 × 0.66 × 8.4 × 10−10 × 10−3 10−12 = 0.1968 N/m2 =
(46)
As stated in a previous publication [39], both the Zhao number and its critical value (i.e., the critical Zhao number) of an NAPL dissolution system are needed to assess the dissolution front instability of the NAPL dissolution system. If the Zhao number of an NAPL dissolution system is greater than the corresponding critical Zhao number of the NAPL dissolution system, then the NAPL dissolution front is unstable in the NAPL dissolution system. From the previous studies
176
Comput Geosci (2011) 15:167–183
[39], the Zhao number of the NAPL dissolution system under consideration is: " ∗ p L 1 V axf Zh = − paxf = − axf∗ = √ pa β φτ Dm 0 1.16 × 10−5 = √ 0.355 × 0.66 × 8.4 × 10−10 × 2.352 × 10−2 = 5.392
(47)
To establish a computational model for the NAPL dissolution system under consideration, other dimensionless quantities can be evaluated as follows: Lx =
Lx 0.3 = = 3,260.9, L∗ 9.2 × 10−5
αT =
(48)
αL 7.2 × 10−4 = = 7.826, L∗ 9.2 × 10−5 αT 1.44 × 10 = = 1.5652 ∗ L 9.2 × 10−5
(49)
The critical Zhao number (i.e., Zhcritical ) of the NAPL dissolution system can be determined from the following characteristic equation [39]: ⎡ ⎛ 1 ⎣ (1 + β) Zhcritical − (1 + β) − ⎝ 1 + α L Zhcritical 1 + α L Zhcritical ⎞ ⎤ 1 + α T Zhcritical ⎠ −α L m (1 − β)⎦ σ1 1 + α L Zhcritical "
+ ⎣⎝m
− α L m2 (1 + α T Zhcritical )⎠
σ1 =
β=
186φ 0.87 Dm β0 = d250
Zhcritical 1+α L Zhcritical
2
β0 = ⎣
(1 − β) ⎦=0 1 + α L Zhcritical
+
1+α L Zhcritical
+ Zhcritical 1+α L Zhcritical
≈
0.87
186φ Dm d250
(55)
d50 ρa Zh μa (1 − Sn0 )
$
"
τ Dm φ
d50 ρa Zhcritical μa (1 − Sn0 )
3,301φ 0.798 τ 0.55 D2.1 m d250
%0.71 β00.355 (56)
"
τ Dm φ
ρa Zhcritical μa (1 − Sn0 )
1 %0.71 ⎤ 0.645
⎦
1.1 (57)
Finally, L∗ is expressed as a function of β 0 . " ∗
L =
φτ Dm β0
(58)
Consideration of Eqs. 57 and 58 yields the following equation:
2
ψ (Sn0 ) k (Sn0 ) = ψ (Snf ) k (Snf )
τ Dm β0 φ
Note that in the critical state, Zh = Zhcritical , so that β 0 can be expressed as follows:
⎤
4m2 1+α T Zhcritical
$
⎡
(50) where #
"
This leads to the following equation:
1 + α T Zhcritical Zhcritical + 1 + α L Zhcritical 1 + α L Zhcritical ⎞
(54)
Using this equation, the Reynolds number of the NAPL dissolution system can be expressed as d50 ρa Zh Re = μa (1 − Sn0 )
−4
⎡⎛ "
where Vaxf is the Darcy velocity of the aqueous phase fluid after the NAPL is completely dissolved in the NAPL dissolution system. Equation 53 can be straightforwardly rewritten in the following form:
Vaxf = Zh φτ Dm β0
Ly 1.0 Ly = ∗ = = 10,869.6 L 9.2 × 10−5 αL =
Since m, α L , and α L are dependent on L∗ , it is necessary to establish a relationship between L∗ and Zhcritical below. Based on the dimensionless quantities expressed in Eqs. 14 and 15, the dimensionless Zhao number is defined as follows: " ∗ p L V 1 axf (53) Zh = − paxf = − axf∗ = √ pa β φτ Dm 0
μa (1 − Sn0 )
d250 φ 0.202 τ 0.45 3,301
1.11
(51)
Zhcritical =
(52)
For an NAPL dissolution system of a finite domain, the critical geometrical characteristic length (L y ) of the
ρa Dm (L∗ )1.818
(59)
Comput Geosci (2011) 15:167–183
177
system should be equal to half the longest wavelength (λ) of the perturbation that can grow in the system, so that the following condition exists: λ = Ly 2
(60)
Since m = 2π = Lπy , where m is the wavenumber of λ the perturbation, the dimensionless wavenumber of the perturbation can be expressed as follows: m = mL∗ =
π ∗ π L = ∗ Ly Ly
(61) &
Consideration of Eqs. 59 and 61 simultaneously yields the following equation:
Zhcritical
μa (1 − Sn0 ) π 1.818 = 1.818 ρa Dm mL y
d250 φ 0.202 τ 0.45 3,301
0
Zh =
d50 Vaxf 57.45φ 0.899 τ 0.775 D1.55 m
μa (1 − Sn0 ) ρa
1 0.55 ' 1.55
(63)
1.11 (62)
Equation 62 indicates that for an NAPL dissolution system of a finite domain, the critical Zhao number of the system is inversely proportional to the critical geometrical characteristic length (L y ) of the system. Thus, with an increase in the critical geometrical char-
Fig. 4 Comparison of current results with the previous ones due to different longitudinal and transversal dispersivities
acteristic length (L y ) of the system, there is a decrease in the critical Zhao number of the system. This implies that for the perturbation of a given wavelength, the NAPL dissolution system becomes more unstable with the increase of the critical geometrical characteristic length (L y ) of the system. Substituting Eq. 57 into Eq. 53 yields the Zhao number of an NAPL dissolution system when the correlation of Imhoff et al. [8] can be used for describing the mass transfer rate coefficient of the NAPL dissolution system.
Again, Eq. 63 indicates that the Zhao number of an NAPL dissolution system is directly proportional to the injected Darcy velocity of the aqueous phase fluid in the NAPL dissolution system. As a result, with an increase in the injected Darcy velocity of the aqueous phase fluid in an NAPL dissolution system, there is an increase in the Zhao number of the NAPL dissolution
y
x
(Previous result, α L = 7.2 × 10−4 m,αT = 1.44 × 10−4 m )
(Previous result, α L = 7.2 × 10−4 m,αT = 7.2 × 10−4 m )
(Current result, α L = 7.2 × 10−4 m,αT = 1.44 × 10−4 m )
(Current result, α L = 7.2 × 10−4 m,αT = 7.2 × 10−4 m )
(A) Vaxf
(Current result, α L = 7.2 × 10−4 m,αT = 1.44 × 10−4 m )
(B) Vaxf
= 1.16 ×10−5 m / s
(Current result, α L = 7.2 × 10−4 m,αT = 7.2 × 10−4 m ) = 2.52 × 10−5 m / s
178
system. This means that the NAPL dissolution system becomes more unstable with an increase in the injected Darcy velocity of the aqueous phase fluid within the NAPL dissolution system. On the other hand, for a given Zhao number, the injected Darcy velocity of the aqueous phase fluid within the NAPL dissolution system can be straightforwardly determined from either Eq. 54 or Eq. 63. Thus, the Zhao number of an NAPL dissolution system, rather than the injected Darcy velocity of the aqueous phase fluid within the NAPL dissolution system, is used to investigate three kinds of potential fingering modes in Section 4 later. If Eq. 57 is substituted into Eq. 58, then L∗ is expressed as a function of Zhcritical so that the critical Zhao number of the NAPL dissolution system can be determined by solving the corresponding characteristic equation (i.e., Eq. 50) of the system. As a result, the solution for the critical Zhao number of the NAPL dissolution system is equal to 0.011 approximately. Since the Zhao number of the NAPL dissolution system is about 5.4 (see Eq. 47), which is much greater than the critical Zhao number, the system under consideration is in an unstable state. Note that although both the NAPL morphological profile evolution and the NAPL morphological interface roughness evolution can be used to compare the current simulation results with the previous ones [11], the former is adopted in this investigation because it is the NAPL morphological profile evolution that can be used to show the three different kinds of fingering modes, as demonstrated later in Section 4. However, for the simulation of NAPL dissolution front instability in heterogeneous porous media, the latter may be the best way to represent the simulation results, since the finger locations are not necessarily in the same locations for the random heterogeneous porous medium. Figure 4 shows the comparison of the current simulation results with the previous solutions that can be viewed as a representative of the related laboratory experimental observations [11]. Note that in this figure, the previous solutions are presented in a black-white form, so that the color code is only applicable to the current simulation results. Two cases are considered for comparison. In the first case, the longitudinal and transverse dispersivities are 7.2 × 10−4 m and 1.44 × 10−4 m, while in the second case, both of them are equal to 7.2 × 10−4 m. To examine the possible effect of the Darcy velocity on the fingering patterns of the NAPL dissolution front, two different Darcy velocities, namely 1.16 × 10−5 m/s (equivalent to 1 m/day) and 2.52 × 10−5 m/s (equivalent to 2.18 m/day), are used to produce the related simulation results. Since the cur-
Comput Geosci (2011) 15:167–183
rent simulation results compare well with the previous solutions, it demonstrates that the proposed numerical procedure is suitable and useful for simulating the morphological evolution of NAPL dissolution fronts in supercritical systems.
4 Three different kinds of modes associated with morphological evolution of NAPL-dissolution fronts in supercritical systems After the proposed numerical procedure is verified through two benchmark problems, it is used to investigate three potential modes of NAPL dissolution fronts in the supercritical systems of finite domains. Compared with previous publications, this is the major contribution of this paper. The geometrical and boundary conditions of the computational model are the same as those used in the first benchmark problem (see Fig. 1). To examine three different modes of the NAPL dissolution fronts in supercritical systems, three cases with two different dimensionless wavenumbers, namely 0.005 and 0.001, are considered in the corresponding computations. In the first two cases, the dimensionless wavenumber is 0.005, while in the third case, the dimensionless wavenumber is 0.001. Common parameters used in three cases are as follows: the initial saturation (Sn0 ) of the NAPL is 0.2; the irreducible saturation (Sai ) of the aqueous phase fluid is 0.15; the dimensionless longitudinal and transverse dispersivities (α L and α T ) are 0.2 and 0.02, respectively; the ratio (ε) of the equilibrium concentration of the NAPL species in the aqueous phase fluid to the den sity of the NAPL itself is 0.001; the density ratio ρ a of the aqueous phase fluid to the NAPL is 1.0/1.46; the value of β 1 is 0.87. To simulate the propagation of NAPL dissolution fronts appropriately, the whole computational domain is simulated by 99,301 four-node square elements of 100,000 nodal points in total. The initial residual saturation field of the NAPL is randomly perturbed by a small amount of 1% of the originally input saturation of the NAPL (i.e., Sn0 = 0.2) before running the computational model. This means that the resulting initial residual saturation is of a random distribution, which has a mean value of the homogeneous residual saturation (i.e., Sn0 = 0.2) and a variation of 0.002 (i.e., 1% of Sn0 = 0.2) in the whole computational domain. Due to the consideration of three different ratios of the Zhao number to the critical Zhao number, three different kinds of modes, namely the fundamental mode, the fingering mode, and the fractal mode, have been obtained from the corresponding computational models.
Comput Geosci (2011) 15:167–183
179
4.1 The fundamental mode The fundamental mode is obtained in the first case. The characteristic of this case is that the Zhao number of the NAPL dissolution system is in the lower range of the supercritical Zhao numbers so that the Zhao number and its critical value have the same order of magnitude. The Zhao number used in this case is 0.06, while the critical Zhao number is 0.017. This means that the ratio of the Zhao number to its critical value is equal to 3.53. The dimensionless lengths of the computational domain are 628 and 1,570 in the x and y directions, respectively. In this case, the dimensionless time step is 11.0 in the computational simulation. Figure 5 shows the evolution process of the fundamental mode associated with the NAPL dissolution front in a lower supercritical Zhao number system. This mode is characterized by half the longest wavelength
of the perturbation that the computational model can accommodate in the y direction. Note that the fundamental mode of the NAPL dissolution front is similar to that of temperature in a convective instability system within the crust of the Earth [15, 37]. Since the Zhao number of the NAPL dissolution system is in the lower range of supercritical Zhao numbers, the amplitude of the fundamental mode grows very slowly so that it takes a considerable time for the mode to be viewed with naked eyes. For example, this mode can be clearly identified after the dimensionless time is equal to 3,850. With an increase in the dimensionless time, the amplitude of the fundamental mode grows accordingly. 4.2 The fingering mode The (normal) fingering mode is obtained in the second case. The characteristic of this case is that the Zhao
Fig. 5 The evolution processes of the fundamental mode (the first kind of mode) for an NAPL dissolution front in the supercritical system
(t
= 1100 )
(t
= 2750 )
(t
= 3850 )
(t
= 4400 )
(t
= 4950 )
(t
= 5500 )
(t
= 2200 )
(t
= 3300 )
180
Comput Geosci (2011) 15:167–183
number of the NAPL dissolution system is in the middle range of the supercritical Zhao numbers so that the Zhao number is one order of magnitude higher than its critical value. The Zhao number used in this case is 0.5, while the critical Zhao number remains 0.017 because the same dimensionless wavenumber is used for both the first and the second cases. In this situation, the ratio of the Zhao number to its critical value is equal to 29.4. The dimensionless lengths of the computational domain are the same as those used in the first case. In this situation, the dimensionless time step is 1.3 in the computational simulation. Figure 6 shows the evolution process of the fingering mode for the NAPL dissolution front in a moderate supercritical Zhao number system. This kind of mode is characterized by the pattern of a few fingers in the
NAPL dissolution system. Since the Zhao number of the NAPL dissolution system is in the middle range of supercritical Zhao numbers, the length of the fingers (i.e., the amplitude of the fingering mode) grows relatively faster, compared with that of the fundamental mode. However, the width of the finger still grows very slowly. This kind of mode can be clearly observed after the dimensionless time is equal to 260 in Fig. 6. Similarly, with an increase in the dimensionless time, the length of the finger grows gradually. 4.3 The fractal mode The fractal mode is obtained in the third case. The characteristic of this case is that the Zhao number of the NAPL dissolution system is in the higher range of
Fig. 6 The evolution processes of the fingering mode (the second kind of mode) for an NAPL dissolution front in the supercritical system
t = 130 )
(t
= 195 )
(t
= 260 )
(t
= 335 )
(t
= 390 )
(t
= 455 )
(t
= 520 )
(t
= 585 )
(
Comput Geosci (2011) 15:167–183
181
the supercritical Zhao numbers of the system so that the ratio of the Zhao number to its critical value is at least one order of magnitude higher than that associated with the second case (i.e., the normal fingering mode case). The Zhao number used in this case is 1.0, while the critical Zhao number is 3.34 × 10−3 . This implies that the ratio of the Zhao number to its critical value is equal to 283.3. The dimensionless lengths of the computational domain are 3,142 and 4,713 in the x and
y directions, respectively. In this case, the dimensionless time step is 1.9 in the computational simulation. Figure 7 shows the evolution process of the fractal mode associated with the NAPL dissolution front in a higher supercritical Zhao number system. This kind of mode has been commonly observed in the previous laboratory experiments and numerical simulations [10, 14]. It is characterized by the fractal pattern consisting of several irregular fingers. As the Zhao number of
Fig. 7 The evolution processes of the fractal mode (the third kind of mode) for an NAPL dissolution front in the supercritical system
(t
= 95 )
(t
= 285 )
(t
= 475 )
(t
= 570 )
(t
= 760 )
(t
= 950 )
(t
= 190 )
(t
= 380 )
182
the NAPL dissolution system is in the higher range of supercritical Zhao numbers, the amplitude of the fractal mode grows very fast, compared with that of the normal finger mode that is comprised of a few smooth fingers. At the early stage of the simulation, many small fingers of shorter wavelengths are generated in the computational model (see t = 95). As time goes on, these small fingers grow in length and width, sometimes merging to form longer and wider fingers of irregular fractal shapes. As shown in Fig. 7, this kind of mode can be observed after the dimensionless time is equal to 190 in the computational simulation.
5 Conclusions To investigate the potential evolution patterns of NAPL dissolution fronts in two-dimensional fluidsaturated porous media of finite domains, a numerical procedure consisting of a combination of the finite element and finite difference methods is proposed as a computational simulation tool. In the proposed numerical procedure, the finite difference method is used to discretize time, while the finite element method is used to discretize space. Two benchmark problems, for which either analytical results or previous solutions are available, are used to verify the proposed numerical procedure. The related simulation results from these two benchmark problems have demonstrated that the proposed numerical procedure is useful and applicable for simulating the morphological evolution of NAPL dissolution fronts in two-dimensional fluid-saturated porous media of finite domains. The proposed numerical procedure is used to investigate the potential morphological evolution patterns of NAPL dissolution fronts in supercritical NAPL dissolution systems. It has been recognized that depending upon the relative position where the Zhao number of an NAPL dissolution system is located in the range of the supercritical Zhao numbers, three kinds of modes, namely the fundamental mode, the (normal) fingering mode, and the fractal mode, can take place for NAPL dissolution fronts in two-dimensional fluidsaturated porous media. The related simulation results have demonstrated that (1) if the Zhao number of an NAPL dissolution system is of the same order of magnitude as the critical Zhao number of the system, then the fundamental mode is predominant; (2) if the Zhao number is one order of magnitude higher than the critical Zhao number, then the (normal) fingering mode is the predominant pattern of the NAPL dissolution front; and (3) if the Zhao number is two orders of mag-
Comput Geosci (2011) 15:167–183
nitude higher than the critical Zhao number, then the fractal mode is predominant for the NAPL dissolution front. Acknowledgements This work is financially supported by the Natural Science Foundation of China (Grant no: 10872219). The authors express their thanks to the anonymous referees for their valuable comments, which led to a significant improvement over an early version of the paper.
References 1. Alt-Epping, P., Smith, L.: Computing geochemical mass transfer and water/rock ratios in submarine hydrothermal systems: implications for estimating the vigour of convection. Geofluids 1, 163–181 (2001) 2. Bear, J.: Dynamics of Fluids in Porous Media. Elsevier, Amsterdam (2001) 3. Chadam, J., Hoff, D., Merino, E., Ortoleva, P., Sen, A.: Reactive infiltration instabilities. IMA J. Appl. Math. 36, 207–221 (1986) 4. Chadam, J., Ortoleva, P., Sen, A.: A weekly nonlinear stability analysis of the reactive infiltration interface. IMA J. Appl. Math. 48, 1362–1378 (1988) 5. Chen, J.S., Liu, C.W.: Numerical simulation of the evolution of aquifer porosity and species concentrations during reactive transport. Comput. Geosci. 28, 485–499 (2002) 6. Chen, J.S., Liu, C.W., Lai, G.X., Ni, C.F.: Effects of mechanical dispersion on the morphological evolution of a chemical dissolution front in a fluid-saturated porous medium. J. Hydrol. 373, 96–102 (2009) 7. Geller, J.T., Hunt, J.R.: Mass transfer from nonaqueous phase organic liquids in water-saturated porous media. Water Resour. Res. 29, 833–845 (1993) 8. Imhoff, P.T., Jaffe, P.R., Pinder, G.F.: An experimental study of complete dissolution of a nonaqueous phase liquid in saturated porous media. Water Resour. Res. 30, 307–320 (1994) 9. Imhoff, P.T., Miller, C.T.: Dissolution fingering during the solubilization of nonaqueous phase liquids in saturated porous media: 1. Model predictions. Water Resour. Res. 32, 1919–1928 (1996) 10. Imhoff, P.T., Thyrum, G.P., Miller, C.T.: Dissolution fingering during the solubilization of nonaqueous phase liquids in saturated porous media: 2. Experimental observations. Water Resour. Res. 32, 1929–1942 (1996) 11. Imhoff, P.T., Farthing, M.W., Gleyzer, S.N., Miller, C.T.: Evolving interface between clean and nonaqueous phase liquid (NAPL)-contaminated regions in two-dimensional porous media. Water Resour. Res. 38, 1093–1106 (2002) 12. Imhoff, P.T., Farthing, M.W., Miller, C.T.: Modeling NAPL dissolution fingering with upscaled mass transfer rate coefficients. Adv. Water Resour. 26, 1097–1111 (2003) 13. Miller, C.T., Poirier-McNeil, M.M., Mayer, A.S.: Dissolution of trapped nonaqueous phase liquids: mass transfer characteristics. Water Resour. Res. 26, 2783–2796 (1990) 14. Miller, C.T., Gleyzer, S.N., Imhoff, P.T.: Numerical modeling of NAPL dissolution fingering in porous media. In: Selim, H.M., Ma, L. (eds.) Physical Nonequilibrium in Soils: Modeling and Application, pp. 389–415. Ann Arbor Press, Michigan (1998) 15. Nield, D.A., Bejan, A.: Convection in Porous Media. Springer, New York (1992)
Comput Geosci (2011) 15:167–183 16. Ormond, A., Ortoleva, P.: Numerical modeling of reactioninduced cavities in a porous rock. J. Geophys. Res. 105, 16737–16747 (2000) 17. Ortoleva, P., Chadam, J., Merino, E., Sen, A.: Geochemical self-organization II: the reactive-infiltration instability. Am. J. Sci. 287, 1008–1040 (1987) 18. Powers, S.E., Abriola, L.M., Weber, W.J. Jr.: An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems: transient mass transfer rates. Water Resour. Res. 30, 321–332 (1994) 19. Raffensperger, J.P., Garven, G.: The formation of unconformity-type uranium ore deposits: coupled hydrochemical modeling. Am. J. Sci. 295, 639–696 (1995) 20. Renard, F., Gratier, J.P., Ortoleva, P., Brosse, E., Bazin, B.: Self-organization during reactive fluid flow in a porous medium. Geophys. Res. Lett. 25, 385–388 (1998) 21. Schafer, D., Schafer, W., Kinzelbach, W.: Simulation of reactive processes related to biodegradation in aquifers: 1. Structure of the three-dimensional reactive transport model. J. Contam. Hydrol. 31, 167–186 (1998) 22. Schafer, D., Schafer, W., Kinzelbach, W.: Simulation of reactive processes related to biodegradation in aquifers: 2. Model application to a column study on organic carbon degradation. J. Contam. Hydrol. 31, 187–209 (1998) 23. Seyedabbasi, M.A., Farthing, M.W., Imhoff, P.T., Miller, C.T.: The influence of wettability on NAPL dissolution fingering. Adv. Water Resour. 31, 1687–1696 (2008) 24. Soerens, T.S., Sabatini, D.A., Harwell, J.H.: Effects of flow bypassing and nonuniform NAPL distribution on the mass transfer characteristics of NAPL dissolution. Water Resour. Res. 34, 1657–1673 (1998) 25. Steefel, C.I., Lasaga, A.C.: Evolution of dissolution patterns: permeability change due to coupled flow and reaction. In: Melchior, D.C., Basset, R.L. (eds.) Chemical Modeling in Aqueous Systems II. Am. Chem. Soc. Symp. Ser., vol. 416, pp. 213–225 (1990) 26. Steefel, C.I., Lasaga, A.C.: A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. Am. J. Sci. 294, 529–592 (1994) 27. Willson, C.S., Hall, J.L., Miller, C.T., Imhoff, P.T.: Factors affecting bank formation during surfactant-enhanced mobilization of residual NAPL. Environ. Sci. Technol. 33, 2440– 2446 (1999) 28. Yeh, G.T., Tripathi, V.S.: A model for simulating transport of reactive multispecies components: model development and demonstration. Water Resour. Res. 27, 3075–3094 (1991)
183 29. Zhao, C., Xu, T.P., Valliappan, S.: Numerical modeling of mass transport problems in porous media: a review. Comput. Struct. 53, 849–860 (1994) 30. Zhao, C., Hobbs, B.E., Mühlhaus, H.B.: Finite element modelling of temperature gradient driven rock alteration and mineralization in porous rock masses. Comput. Methods Appl. Mech. Eng. 165, 175–186 (1998) 31. Zhao, C., Hobbs, B.E., Mühlhaus, H.B., Ord, A.: Finite element modelling of rock alteration and metamorphic process in hydrothermal systems. Commun. Numer. Methods Eng. 17, 833–843 (2001) 32. Zhao, C., Hobbs, B.E., Mühlhaus, H.B., Ord, A., Lin, G.: Finite element modeling of three-dimensional steady-state convection and lead/zinc mineralization in fluid-saturated rocks. J. Comput. Methods Sci. Eng. 3, 73–89 (2003) 33. Zhao, C., Hobbs, B.E., Ord, A., Peng, S., Mühlhaus, H.B., Liu, L.: Numerical modeling of chemical effects of magma solidification problems in porous rocks. Int. J. Numer. Methods Eng. 64, 709–728 (2005) 34. Zhao, C., Hobbs, B.E., Hornby, P., Ord, A., Peng, S.: Numerical modelling of fluids mixing, heat transfer and non-equilibrium redox chemical reactions in fluid-saturated porous rocks. Int. J. Numer. Methods Eng. 66, 1061–1078 (2006) 35. Zhao, C., Hobbs, B.E., Ord, A., Hornby, P., Peng, S.: Effect of reactive surface areas associated with different particle shapes on chemical-dissolution front instability in fluidsaturated porous rocks. Transp. Porous Media 73, 75-94 (2008) 36. Zhao, C., Hobbs, B.E., Hornby, P., Ord, A., Peng, S., Liu, L.: Theoretical and numerical analyses of chemical-dissolution front instability in fluid-saturated porous rocks. Int. J. Numer. Anal. Methods Geomech. 32, 1107–1130 (2008) 37. Zhao, C., Hobbs, B.E., Ord, A.: Convective and Advective Heat Transfer in Geological Systems. Springer, Berlin (2008) 38. Zhao, C., Hobbs, B.E., Ord, A.: Fundamentals of Computational Geoscience: Numerical Methods and Algorithms. Springer, Berlin (2009) 39. Zhao, C., Hobbs, B.E., Ord, A.: Theoretical analyses of nonaqueous-phase-liquid dissolution induced instability in two-dimensional fluid-saturated porous media. Int. J. Numer. Anal. Methods Geomech. (2009). doi:10.1002/nag.880 40. Zhao, C., Hobbs, B.E., Ord, A., Peng, S.: Effects of mineral dissolution ratios on chemical-dissolution front instability in fluid-saturated porous media. Transp. Porous Media 82, 317– 335 (2010) 41. Zienkiewicz, O.C.: The Finite Element Method. McGrawHill, London (1977)
