Journal of Contaminant Hydrology 107 (2009) 162–174

Contents lists available at ScienceDirect

Journal of Contaminant Hydrology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j c o n h yd

Analytical solution of two-dimensional solute transport in an aquifer–aquitard system Hongbin Zhan a,⁎, Zhang Wen b,c, Guanhua Huang b,c, Dongmin Sun d a b c d

Department of Geology and Geophysics, Texas A&M University, College Station, TX 77843-3115, United States Department of Irrigation and Drainage, College of Water Conservancy and Civil Engineering, China Agricultural University, Beijing, 100083, PR China Chinese–Israeli International Center for Research in Agriculture, China Agricultural University, Beijing, 100083, PR China Department of Environmental Sciences, University of Houston-Clear Lake, Houston, TX 77058-1098, United States

a r t i c l e

i n f o

Article history: Received 6 March 2008 Received in revised form 21 April 2009 Accepted 28 April 2009 Available online 7 May 2009 Keywords: Aquitard diffusion Aquitard advection Mass conservation Solute transport Laplace transform

a b s t r a c t This study deals with two-dimensional solute transport in an aquifer–aquitard system by maintaining rigorous mass conservation at the aquifer–aquitard interface. Advection, longitudinal dispersion, and transverse vertical dispersion are considered in the aquifer. Vertical advection and diffusion are considered in the aquitards. The first-type and the thirdtype boundary conditions are considered in the aquifer. This study differs from the commonly used averaged approximation (AA) method that treats the mass flux between the aquifer and aquitard as an averaged volumetric source/sink term in the governing equation of transport in the aquifer. Analytical solutions of concentrations in the aquitards and aquifer and mass transported between the aquifer and upper or lower aquitard are obtained in the Laplace domain, and are subsequently inverted numerically to yield results in the real time domain (the Zhan method). The breakthrough curves (BTCs) and distribution profiles in the aquifer obtained in this study are drastically different from those obtained using the AA method. Comparison of the numerical simulation using the model MT3DMS and the Zhan method indicates that the numerical result differs from that of the Zhan method for an asymmetric case when aquitard advections are at the same direction. The AA method overestimates the mass transported into the upper aquitard when an upward advection exists in the upper aquitard. The mass transported between the aquifer and the aquitard is sensitive to the aquitard Peclet number, but less sensitive to the aquitard diffusion coefficient. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Interaction between aquifers and aquitards is an important process affecting flow and transport in subsurface flow systems. Most aquitards consist of silt and clay and are well capable of storing water and solute, due to their large values of porosity. Thus when a solute in an aquifer contacts a previously solutefree aquitard, a concentration gradient exists across the aquiferaquitard interface; and molecular diffusion will drive the solute into the aquitard. Furthermore, leakage often exists across the aquitard, thus advection in the aquitard will be another important mechanism for solute transport there. ⁎ Corresponding author. Tel.: +1 979 862 7961; fax: +1 979 845 6162. E-mail address: [email protected] (H. Zhan). 0169-7722/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2009.04.010

Diffusion at the aquifer–aquitard interface is somewhat similar to diffusion at the matrix-fracture boundary (Tang et al., 1981; Sudicky and Frind, 1982, 1984; Fujikawa and Fukui, 1990; Liu et al., 2004). The difference is that the aperture of a fracture is much smaller than the aquifer thickness and the flow velocity in the fracture is often much greater than that in the aquifer under the same hydraulic gradient. Many studies on fractured media have shown that matrix diffusion is the primary factor for retarding contaminants in the fractures (e.g. Neretnieks, 1980; Rasmuson and Neretnieks, 1981; Neretnieks et al., 1982; Moreno et al., 1985; Liu et al., 2004). Aquitard diffusion was shown to be a controlling factor affecting solute transport in laboratory experiments by Sudicky et al. (1985), Starr et al. (1985), Young and Ball (1998), and in numerous field aquifer studies such as Gillham et al. (1984), Johnson et al. (1989), Ball

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

et al. (1997a,b), Liu and Ball (1999), Hendry et al. (2003), Hunkeler et al. (2004), Parker et al. (2004) and others. After passing of the solute front in the aquifer, back diffusion from the aquitard to the aquifer is the primary cause of the tailing effect observed in the aquifer, which has caused great difficulty for contaminant remediation (Liu and Ball, 2002). Another important reason causing the tailing effect in real aquifer–aquitard systems is diffusion into and out of low-permeability lenses within the aquifer materials (Gillham et al., 1984). Sudicky et al. (1985) have investigated the aquitard diffusion effect in an artificial sandy aquifer whose thickness is about 0.02 to 0.03 m in the laboratory. Because the aquifer thickness is so thin in respect to the horizontal transport distance, they could approximate the diffusive mass flux at the aquifer–aquitard interface as a volumetric source/sink term in the governing equation of solute transport in the aquifer. Chen (1985) and Tang and Aral (1992a,b) have also adopted this approach to study dispersion–diffusion in an aquifer–aquitard system with radial and uniform flows, respectively. The implication of this approach is that the transverse mixing across the aquifer thickness is so rapid that a thickness-averaged concentration can be used. The same approximation has been broadly adopted in studying matrix diffusion in fractured media (Neretnieks, 1980; Rasmuson and Neretnieks, 1981; Tang et al., 1981; Neretnieks et al., 1982; Sudicky and Frind, 1982, 1984; Fujikawa and Fukui, 1990; Johns and Roberts, 1991; Liu et al., 2004). Sudicky et al. (1985) have realized that in real aquifers, the transverse mixing is probably not always rapid enough to warrant the usage of a thickness-averaged approach. Therefore, they proposed an alternative method of treating the diffusive flux at the aquifer–aquitard interface as a boundary condition rather than a source/sink term in the governing equation of transport in the aquifer. But to make the work amenable to analytical solution, Sudicky et al. (1985) have neglected the longitudinal dispersion and only considered the transverse vertical dispersion in the aquifer. Such a treatment might provide a satisfactory interpretation of conservative solute transport in a very thin artificial sandy aquifer used in the experiment of Sudicky et al. (1985). However, for more realistic thicker aquifers, Starr et al. (1985) showed that longitudinal dispersion is important and should be taken into account. In a different paper, Johns and Roberts (1991) have proposed a model for investigating solute transport in largeaperture fractures by considering lateral dispersion to the small aperture regions and diffusion to the rock matrix. Longitudinal dispersion is not considered in that study. Advection in the aquitard is rarely considered in previous analytical solutions although it can be dealt with in numerical simulations (see for example Bester et al., 2005). Most confined aquifer systems are recharged by leakage through the overlying aquitard, which means, in case of transport, that there will be advection across the interface. Advective mass flux in the aquitard may be small, but may be at a similar or greater scale when compared to diffusive flux across the interface. Therefore, to make the analytical solution useful, both advective and diffusive fluxes in the aquitard are considered in this study. It is also important to compare the relative magnitudes of advective versus diffusive fluxes across the interface under realistic conditions. In this study, we will solve for two-dimensional transport in the aquifer and one-dimensional advection–diffusion in

163

the aquitard simultaneously for a fully penetrating, horizontally infinite source without using the averaged approximation employed by Sudicky et al. (1985), Chen (1985), Tang and Aral (1992a,b), and others. Mass balance requirement is rigorously maintained. An example of such a two-dimensional transport scenario is the leaking of toxic materials that are buried in a long trench or a large landfill site. There is no doubt that three-dimensional numerical simulations of flow and transport in complex aquifer–aquitard systems can be carried out with the present-day's computational power. For instance, Martin and Frind (1998) have carried out a threedimensional numerical simulation of groundwater flow and capture zone description for a multiple layer aquifer–aquitard system in the Waterloo Moraine. Bester et al. (2005) have carried out a numerical simulation of road salt impact on an urban wellfield located in an aquifer–aquitard system. The purpose of this paper is to illustrate the essence of the aquitard effect on solute transport from an analytical perspective. Such analytical solutions may serve the purpose of validating the numerical simulations which may suffer from various types of numerical errors. For instance, the numerical errors in numerical models tend to be the largest at the aquifer–aquitard interfaces (Martin and Frind, 1998; Bester et al., 2005). For the sake of simplicity of presentation, we only discuss a conservative solute. 2. Conceptual and mathematical models 2.1. Conceptual model The system investigated is an aquifer bounded at the top and bottom by two aquitards, or an aquifer bounded at the top by an aquitard and at the bottom by impermeable bedrock. The aquifer–aquitard and the bedrock–aquifer boundaries are assumed to be horizontal. The aquifer is homogeneous and horizontally isotropic with constant longitudinal and transverse vertical dispersivities. The aquitards are also homogeneous and sufficiently thick so that solute diffusion is not affected by their thicknesses. This assumption may be reasonable because the penetration depths of solute into the aquitards via diffusion are often limited; meaning that only those regions of the aquitards that are very close to the aquifer will be affected by solute transport. However, if one is dealing with a very thin aquitard, it is possible that solute can penetrate through the entire thickness of the aquitard into the adjacent aquifer. If this is the case, one must consider the affected adjacent aquifer as well. Such a circumstance is of interest but will not be discussed in this article. The aquifer, the aquitards, and the bedrock extend horizontally to infinity. Fig. 1A–B show the schematic diagrams of two cases that will be investigated. They are an aquifer bounded by upper and lower aquitards, and an aquifer bounded by an upper aquitard and lower impermeable bedrock, respectively. The big arrows there show the groundwater flow directions. The hydraulic conductivity of the aquitard is assumed to be a few orders of magnitude smaller than that of the aquifer, thus advective flow in the aquitard is nearly perpendicular to the aquifer–aquitard interface. Flow in the aquifer is nearly horizontal. We set up the coordinate system as follows. The x- and z-axes are along the horizontal and vertical directions, respectively with the origin at the left boundary. We choose the

164

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

tration can be easily calculated via the resident concentration using the method of van Genuchten and Parker (1984). The solution is based on the assumption of a constant horizontal velocity in the aquifer and constant vertical velocities in the aquitards. This will lead to constant dispersion coefficients in the aquifer and the aquitards. In reality, velocities in the aquifer and the aquitards could be variable and so are the dispersion coefficients. A numerical simulation is probably more useful when variable velocities are encountered. It is assumed that the aquifer and the aquitard are free from solutes at the start of transport. Based on the conceptual model, the following mathematical model is established. 2.2. Solute transport in an aquifer bounded by upper and lower aquitards The governing equations and the initial and boundary conditions for the aquifer and aquitards of this case are as follows. For the aquifer,

Fig. 1. Schematic diagram of: A) an aquifer bounded at top and bottom by aquitards; B) an aquifer bounded at top by an aquitard and at the bottom by impermeable bedrock.

origin of the coordinate system at the center of the left boundary of the aquifer in Fig. 1A, whereas it is at the bottom of the left boundary of the aquifer in Fig. 1B. We will deal with two kinds of boundary conditions, namely the first-type and the third-type conditions for the source. At first, a constant concentration (the first-type) is maintained at the left boundary for the entire thickness of the aquifer. The length of the source in the direction transverse to the x–z plane is assumed to be sufficiently large that dispersion in that direction is not considered. The first-type boundary condition has been commonly used in many contaminant transport problems described by Huyakorn et al. (1987), Leij et al. (1991), Batu (1996), and others. However, some investigators such as Kreft and Zuber (1978, 1979) and van Genuchten and Parker (1984) have pointed out that the first-type boundary condition might not correctly predict the resident concentration, although it can correctly predict the flux concentration. They have shown that the third-type or the flux-type boundary condition is probably more suitable for predicting the resident concentration in a column test. The difference between the first-type and the third-type boundary conditions is expected to be small when the concentration gradient near the boundary is small. In real aquifer contamination problems, both boundaries seem possible, depending on the site-specific conditions. In this article, we will derive solutions for both types of boundary conditions. Unless stated otherwise, the concentration is the resident concentration. If needed, the flux concen-

AC A2 C A2 C AC = Dx 2 + Dz 2 − v ; At Ax Ax Az

ð1Þ

C ðx = 0; z; t Þ = C0 ;

ð2Þ

C ðx = + ∞; z; t Þ = 0;

ð3Þ

C ðx; z; t = 0Þ = 0; for x N 0;

ð4Þ

and for the upper aquitard, AC1 A2 C1 AC = D01 − v1z 1 ; At Az Az2

ð5Þ

C1 ðx; z = B; t Þ = C ðx; z = B; t Þ;

ð6Þ

θ1 v1z C1 ðx; z = B; t Þ − θ1 D01 = − θDz

AC1 ðx; z = B; t Þ Az

ð7Þ

AC ðx; z = B; t Þ ; Az

C1 ðx; z = + ∞; t Þ = 0;

ð8Þ

C1 ðx; z; t = 0Þ = 0;

ð9Þ

and for the lower aquitard, AC2 A2 C2 AC = D02 − v2z 2 ; At Az Az2

ð10Þ

C2 ðx; z = − B; t Þ = C ðx; z = − B; t Þ;

ð11Þ

θ2 v2z C2 ðx; z = − B; t Þ − θ2 D02 = − θDz

AC2 ðx; z = − B; t Þ Az

ð12Þ

AC ðx; z = − B; t Þ ; Az

C2 ðx; z = − ∞; t Þ = 0;

ð13Þ

C2 ðx; z; t = 0Þ = 0;

ð14Þ

where C, C1, C2 represent the resident concentrations in the aquifer, the upper aquitard, and the lower aquitard respectively;

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

C0 is the constant concentration at x = 0; v is the average pore velocity of groundwater flow in the aquifer, v1z and v2z are the vertical pore velocities of groundwater flow in the upper and lower aquitards, respectively and are positive for upward flow; Dx and Dz are the longitudinal and vertical hydrodynamic dispersion coefficients respectively for the aquifer, and Dx = αxv +D0⁎, Dz =αzv + D0⁎, where αx and αz are the longitudinal and transverse vertical dispersivities, respectively, D0⁎is the effective molecular diffusion coefficient of the solute in the aquifer (Bear, 1972); D01 and D02 are the effective molecular diffusion coefficients in the upper and lower aquitards, respectively; B is the half thickness of the aquifer; θ, θ1, and θ2 are the porosities of the aquifer, the upper aquitard, and the lower aquitard, respectively; and t is time. The effective molecular diffusion coefficients in the upper and lower aquitards are D01 =τ1D0, and D02 =τ2D0, respectively, where τ1 and τ2 are tortuosities of the upper and lower aquitards, respectively, and D0 is the free molecular diffusion coefficient of the solute. Eqs. (6) and (7) are the continuities of concentration and mass flux at the upper aquifer–aquitard interface at z = B, respectively, whereas Eqs. (11) and (12) are the continuities of concentration and mass flux at the lower aquifer–aquitard interface at z = − B, respectively. Both aquitards are thought to be sufficiently thick for the transport problem discussed here. If using the third-type boundary condition, Eq. (2) will be replaced by   AC + = vC : j vC − Dx 0 Ax x = 0

ð15Þ

A point to note is that diffusion along the x-axis is not included in the governing equations of transport in the aquitards (see Eqs. (5) and (10)). This is based on two considerations. First, one purpose of this study is to compare the result with those obtained with the averaged approximation which only consider diffusion along the z-axis in the aquitard as well (Sudicky et al., 1985; Chen, 1985; Tang and Aral, 1992a,b). Second, adding the diffusion term along the x-axis in Eqs. (5) and (10) will make the analytical modeling challenging and complex. From a rigorous physical point of view, neglecting diffusion along the x-axis in the aquitard will certainly introduce some errors to the modeling results. Now the question is: how large are such errors? At the moving front of a contaminant plume where the concentration gradient along the x-axis could be very large, the diffusion term along the x-axis could be comparable with other terms in the equations, thus one is expected to see some non-negligible errors. On the other hand, at regions where the concentration gradients along the x-axis are nearly zeroes, such as at regions far behind the moving front of a contaminant front, the errors of neglecting the diffusion term along the x-axis should be minimized and negligible. Fortunately, many problems associated with aquitard contamination often involve time scales of many decades and contaminant moving fronts that have traveled long distances from persistent sources of contaminant, thus the result developed in this study will be applicable for regions that are close to those sources (or far behind the moving fronts). Nevertheless, one still hopes that an analytical or semi-

165

analytical solution can be developed to consider diffusions in both the x- and z-axes in the future. Defining the dimensionless terms in Table 1 where the subscript “D” represents the dimensionless term, and Pe, Pe1, Pe2 are the Peclet numbers of the aquifer, the upper aquitard, and the lower aquitard, respectively (Bear, 1972). Pe1 or Pe2 will be positive for upward advection and negative for downward advection. After transforming Eqs. (1)–(9) into dimensionless form, applying the Laplace transform to the governing equations and boundary conditions will result in the following equation groups. pC D =

A2 C D A2 C D AC + − Pe D ; 2 AxD AxD Az2D

ð16Þ

C D ðxD = 0; zD ; pÞ = 1 = p;

ð17Þ

C D ðxD = + ∞; zD ; pÞ = 0;

ð18Þ

and for the upper aquitard, pC 1D = β1

A2 C 1D AC − Pe1 1D ; AzD Az2D

ð19Þ

C 1D ðxD ; zD = 1; pÞ = C D ðxD ; zD = 1; pÞ;

ð20Þ

½

AC 1D ðxD ; zD = 1; pÞ − Pe1 C 1D ðxD ; zD = 1; pÞ AzD

=

AC D ðxD ; zD = 1; pÞ ; AzD

σ 1 β1



ð21Þ

C 1D ðxD ; zD = + ∞; pÞ = 0;

ð22Þ

and for the lower aquitard, pC 2D = β2

A2 C 2D AC − Pe2 2D ; AzD Az2D

ð23Þ

C 2D ðxD ; zD = − 1; pÞ = C D ðxD ; zD = − 1; pÞ; "

AC 2D ðxD ; zD = − 1; pÞ − Pe2 C 2D ðxD ; zD = − 1; pÞ AzD

σ 2 β2 =

ð24Þ # ð25Þ

AC D ðxD ; zD = − 1; pÞ ; AzD

C 2D ðxD ; zD = − ∞; pÞ = 0;

ð26Þ

̅ , C 2D ̅ are the Laplace transforms of CD, C1D, C2D, where C D̅ , C 1D respectively, and p is the Laplace transform parameter in respective to the dimensionless time. All the associated terms are explained in Table 1.

Table 1 Dimensionless variables used in this study. qffiffiffiffi β1 = xD = Bx DDxz ; zD = Bz ; tD = DB2z t CD =

C C0

M1D =

; C1D = C0 B2

C1 C0

; C2D =

Mffiffiffiffiffiffiffiffiffiffiffi 1 p ; Dx = Dz

C2 C0

M2D =

;

σ1 = C0 B2

D01 Dz θ1 θ

; β2 =

; σ2 =

D02 Dz θ2 θ

;

Mffiffiffiffiffiffiffiffiffiffiffi 2 p ffiffiffiffiffiffiffiffi ; Pe1 = Pe = pvB Dx = Dz

Dx Dz

;

v1z B Dz ; Pe2

=

v2z B Dz

166

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

The procedure of solving the equation group of Eqs. (16)– ̅ , (22) is given in Appendix A. The derived solutions for C D̅ , C 1D ̅ are: and C 2D 2

∞ X

CD =

An exp4−

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   Pe2 + 4 p + ω2n − Pe 2

n=0

3 xD 5 cosðωn zD + μ n Þ;

ð27Þ

− 1 V zD V 1;

C 1D =

½

∞ X

An exp − 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðPe1 =β 1 Þ2 + 4p = β 1 − Pe1 = β 1

n=0



qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   Pe2 + 4 p + ω2n − Pe

2



ðzD − 1Þ cosðωn + μ n Þ;

zD z 1:

C 2D =

½

∞ X

+

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   Pe2 + 4 p + ω2n − Pe

An exp − 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðPe2 =β 2 Þ2 + 4p = β 2 + Pe2 = β2

n=0

2

ð28Þ

xD

ð29Þ

xD



ðzD + 1Þ cosðωn − μ n Þ;

zD V − 1:

2.4. Numerical consideration

where ωn and μn are the frequency and the phase terms that are determined via Eqs. (A8) and (A9), respectively in Appendix A, and An =

4F1 ðωn ; μ n Þ ; pF2 ðωn ; μ n Þ

ð30Þ

where F1(ωn, μn) and F2(ωn, μn) are given in Eq. (A11). If the third-type boundary condition is used at xD = 0, Eq. (17) will be replaced by 1 AC D CD − Pe AxD

!

j

= 1 = p:

ð31Þ

xD = 0 +

The solutions for the third-type boundary condition are identical to Eqs. (27)–(29) except that An is given by An =

4F1 ðωn ; μ n Þ ; pF2 ðωn ; μ n ÞGðωn Þ

ð32Þ

where G(ωn) is: Gðωn Þ =

with negligible diffusion into the bedrock (Fig. 1B), the concentrations in the aquifer and aquitard can be easily obtained by modifying the solutions derived in Section 2.2. This can be done as follows. One can use the x axis or (z = 0 line) as a symmetric line to make an image aquifer and an image lower aquitard from z = 0 to −∞. The image aquifer and aquitard have the same hydrological parameters as their counterparts, respectively. Furthermore, the advective velocities in the upper aquitard and the lower image aquitard have the same magnitude but opposite direction. The image aquifer and the aquifer are now combined into an equivalent hypothetical aquifer with thickness 2B. The system generated in such a way is a special case of Section 2.2 as mentioned after Eq. (33). Therefore, the problem shown in Fig. 1B can be regarded as the upper half-plane of a hypothetical aquifer of thickness 2B with identical upper and lower aquitard parameters and the same magnitude but opposite aquitard advective velocities. With such a modification, the solutions derived in Section 2.2 can be directly applied for the case of Fig. 1B.

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    1 + 4 p + ω2n = Pe2 + 1 = 2:

ð33Þ

We call the results of Eqs. (27)–(29) the Zhan method hereinafter. For the special case in which the upper and lower aquitards have identical hydrologic parameters, μ n = 0. Furthermore, if the advective velocities in the upper and lower aquitards have the same magnitude but opposite direction, then the system investigated becomes symmetric about the line of z = 0 which can be treated as a no-flux boundary. This can also be seen easily from Eqs. (28) and (29). 2.3. Solute transport in an aquifer bounded by an upper aquitard and lower impermeable bedrock If the aquifer is bounded by an upper aquitard and lower impermeable bedrock, which is treated as a no-flux boundary

Now the solutions are derived in the Laplace domain. Analytical inverse Laplace transform is often difficult, if not impossible for the problems discussed here. Previous studies by Tang et al. (1981), Sudicky et al. (1985), and many others have shown that even under simplified conditions, analytical inverse Laplace transform of the transport problem might end in multiple integrations that can only be calculated numerically. Therefore, we will use numerical inverse Laplace transform to find the concentrations in the real time domain. The numerical inverse Laplace transform can be carried out using the methods developed either by Stehfest (1970), or Talbot (1979), or de Hoog et al. (1982), among others. The Laplace transform Galerkin technique (LGT) proposed by Sudicky (1989) is another highly efficient and accurate numerical inverse Laplace transform method for solving flow and transport. Algorithms by Talbot (1979), de Hoog et al. (1982), and Sudicky (1989) involve complex variables, which are difficult to use in our solutions involving a frequency term ωn and a phase term μn. The Stehfest (1970) algorithm, on the other hand, only requires real values of p. It has been used successfully by several hydrologists in similar problems (Moench and Ogata, 1981; Moench, 1991), but it sometimes suffers from oscillation and convergence problems. Fortunately, this algorithm is found to work well for the problems investigated here when comparing the numerical solutions with analytical solutions under simplified situations. The number of terms used in the Stehfest method is N = 12, which seems to yield the best result for the problems investigated here. Under the special condition that transport in the aquitard is negligible, closed-form analytical solution exists. The results obtained from the numerical inverse Laplace transform method is then tested against the closed-form analytical solution and agreement has been reached. Nevertheless, such an agreement may not be enough to prove that the numerical inverse Laplace transform method is always accurate for the general case that the aquitard transport is non-negligible; it gives us some confidence of using the numerical inverse Laplace transform.

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

167

3. Comparison with previous studies 3.1. Comparison with analytical solutions It can be shown that for no flow and no mass flux across the aquifer–aquitard interfaces, and constant source concentration, the solution simplifies to the well-known one-dimensional transport solution by Ogata and Banks (1961). Without the interface flux, β1 = β2 = 0 and ωn = nπ, where n = 0, 1, 2, ....from Eq. (A6). The only non-zero An term from Eq. (A7) is A0 = 1/p. After this, it is easy to verify that the inverse Laplace transform of (A4) will result in the following solution in real time domain: CD =

   

1 x − Pet x + Pet erfc D pffiffiffiffiffi D + expðλxD Þerfc D pffiffiffiffiffi D ; 2 tD 2 tD 2

ð34Þ

where erfc ( ) is the complementary error function. Eq. (34) is the well-known one-dimensional solution derived by Ogata and Banks (1961). For the third-type boundary condition at xD = 0, A0 = 1/[pG(0)], and An = 0 for n ≠ 0, where G(0) is given by Eq. (33) by setting ωn = 0. It is not difficult to verify that the inverse Laplace transform of the general solution (A4) under this circumstance will become:

CD =

" # !1     1 x − Pet 1 Pe2 tD 2 ðx −PetD Þ2 1 x + Pet 2 − exp − D erfc D pffiffiffiffiffi D + 1 + PexD + Pe tD expðPexD Þerfc D pffiffiffiffiffi D 2 tD π 4tD 2 tD 2 2 4

ð35Þ

Eq. (35) is the well known one-dimensional solution for the third-type boundary derived by van Genuchten (1981) and several others. Detailed derivation of the inverse Laplace transform can be found from van Genuchten (1981) or from the author upon request. 3.2. Averaged approximation of aquitard advection and diffusion In many previous studies, aquitard diffusion is approximated as a volumetric source/sink term in the governing equation of transport in the aquifer, and aquitard advection is not taken into account (Neretnieks, 1980; Tang et al., 1981; Chen, 1985; Sudicky et al., 1985; Fujikawa and Fukui, 1990; Tang and Aral, 1992a,b). In the following, we are going to include advection as part of the mass transport mechanism. The governing Eq. (1) under this approximation is modified as: AC A2 C AC C C = Dx 2 − v − 1 + 2 ; 2θB 2θB At Ax Ax

ð36Þ

where Γ1 and Γ2 are the mass fluxes across the upper and lower interfaces of the aquifer and they include both advective flux and diffusive flux: C1 =

    AC AC θ1 v1z C1 − θ1 D01 1 ; C2 = θ2 v2z C2 −θ2 D02 2 : Az z = B Az z = − B

j

j

ð37Þ

The governing equations of the aquitards are still the same as Eqs. (5) and (10). C and C1 are continuous at the upper aquifer– aquitard interface, whereas C and C2 are continuous at the lower aquifer–aquitard interface. Changing above equations into their dimensionless forms and applying Laplace transforms, one can obtain the following solutions in the Laplace domain for the firsttype boundary condition as: 2 6 6 1 C D = exp6 6− p 4

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Pe21 + 4pβ1 + Pe1 + σ 2 Pe22 + 4pβ2 − Pe2 − Pe Pe2 + 4p + σ 1 2 2

C 1D = C D ðxD ; zD = 1; pÞ × exp4−

C 2D

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðPe1 =β1 Þ2 + 4p = β1 − Pe1 = β1 2

3 7 7 xD 7 7; 5

ð38Þ

3 ðzD − 1Þ5;

2qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðPe2 =β2 Þ2 + 4p = β2 + Pe2 = β2 = C D ðxD ; zD = − 1; pÞ × exp4 ðzD + 1Þ5: 2

ð39Þ

ð40Þ

If ignoring the advective fluxes in the aquitards and assuming the upper and lower aquitards have identical parameters, then the solutions come back to those derived by Tang et al. (1981), Fujikawa and Fukui (1990), and others. The same procedures can be

168

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

applied for the third-type boundary condition and will not be repeated here. Inverse Laplace transform of Eqs. (38)–(40) will yield the concentrations in the real time domain. We call the results of Eqs. (38)–(40) the averaged approximation (AA) method hereinafter. Now we compare the AA method with the Zhan method derived in Section 2.2. 4. Results analysis The results of this study are formulated in dimensionless forms which can be applied to any suitable real parameters. For the sake of helping readers to understand the range of those dimensionless variables, a table of realistically possible values is provided in Table 2 and is briefly illustrated. The effective molecular diffusion coefficients chosen for the aquitard (D01 and D02) and the aquifer (D0⁎) are equivalent to a dilute NaCl solution with D01 = D02 = D0⁎ = 1.16 × 10− 9 m2/s. The average pore velocity v is chosen to be 0.10 m/day or 1.16×10− 6 m/s. The same values of effective molecular diffusion coefficient and average pore velocity were also used in Sudicky et al. (1985). For a field-scale dispersion problem, we choose the longitudinal dispersivity αx =2 m. The transverse vertical dispersivity is expected to be 0.0095 of the longitudinal one as αz =0.019 m. Such a choice is consistent with field-scale dispersion (Domenico and Schwartz, 1998; Fetter, 1999) but is greater than the local dispersivity used in the laboratory experiment of Sudicky et al. (1985) and the field experiment of Sudicky (1986), which is 0.001 m. Therefore, the corresponding dispersion coefficients are Dx =αxv +D0⁎ = 2.32 × 10− 6 m2/s, and Dz =αzv +D0⁎ = 2.32× 10− 8 m2/s. The aquifer thickness is 2B = 4 m. The corresponding Peclet number Pe is 10. Two different leakage velocities in the aquitard are v1z =1.16×10− 10 m/s, and 1.16×10− 9 m/s, respectively, corresponding to Pe1 of 0.1, and 1. The positive sign of v1z (or v2z) indicates that water is moving upward. The porosity difference among the aquitards and aquifer is expected to be a secondary effect, thus the identical value of porosity 0.36 is assigned for all units. The same porosity is used for the aquitard in Sudicky et al. (1985). The first-type boundary condition is used in the following discussion.

4.1. Transport in the aquifer–aquitard system 4.1.1. Concentration profiles in the aquifer–aquitard system Fig. 2A–C show the vertical concentration profiles in the aquifer–aquitard system under conditions of σ1 = σ2 = 1, Pe = 10, xD = 10, and tD = 2. Fig. 2A is for a case with different hydrologic parameters of the upper and lower aquitards with β1 = 0.05, β2 = 0.1, Pe1 = 0.5, and Pe2 = 1. Fig. 2B is for a case Table 2 The default values used in this study. Parameter name

Symbol

Default value

Half aquifer thickness in Fig. 1A or full aquifer thickness in Fig. 1B Effective molecular diffusion coefficients of the aquitard Longitudinal dispersion coefficient of the aquifer Transverse dispersion coefficient of the aquifer Average pore velocity in the aquifer Average pore velocity in the aquitard

B

2m

D01, D02

1.16 × 10− 9 m2/s

Dx

2.32 × 10− 6 m2/s

Dz

2.32 × 10− 8 m2/s

v v1z, v2z

1.16 × 10− 6 m/s 1.16 × 10− 10 m/s 1.16 × 10− 9 m/s 0.36

Porosity

θ = θ1 = θ2

with identical hydrologic parameters of the upper and lower aquitards with β1 = β2 = 0.1, and Pe1 = Pe2 = 1. A point to note is that advective velocities in the upper and lower aquitards are both upward in Fig. 2A and B. Fig. 2C is for a case that is identical to that of Fig. 2B except that the advective velocities in the upper and lower aquitards have opposite direction, i.e., Pe1 = 1 but Pe2 = −1. Fig. 2C is the special case discussed before in Section 2.2 after Eq. (33) with a symmetric system about z = 0. The case shown in Fig. 2C for z N 0 would also be the solution applying for the conceptual model of Fig. 1B. Several interesting points can be made from these three figures. Firstly, as can be seen from the figures, Fig. 2A and B have asymmetric shapes of concentration distributions along the z-axis since the advective velocities in the upper and lower aquitards are both along the upward direction. However, Fig. 2C has a symmetric shape of concentration distribution along the z-axis as expected. For the asymmetric cases of Fig. 2A and B, greater discrepancy exists between the result of the AA method and the Zhan method in the upper aquitard than that in the lower aquitard. Secondly, the AA method only yields a z-independent concentration within the aquifer and a declined concentration with z in the aquitard. Notice that zD = 1 is the interface between the aquifer and the upper aquitard. The Zhan method, however, provides a concentration profile that continuously varies with z. Because of the difference between the dispersion coefficient of the aquifer and the diffusion coefficient of the aquitard, concentration gradients are not continuous at the aquifer–aquitard interface. Fig. 2A–B show that the AA method considerably overestimates the concentration in the upper aquitard up to a fairly large distances from the aquifer–aquitard interface. However, for the symmetric case of Fig. 2C, the degree of such overestimation is much less. Fig. 2C shows that for the symmetric case and for the case with the impermeable bedrock, the difference between the AA solution and our solution is actually quite small. Thirdly, the AA method and the Zhan method seem to converge to negligible concentration at approximately the same distances from the aquifer–aquitard interface. For instance, at zD ≈ 3 in Fig. 2A–C, the concentrations calculated from the Zhan method and the AA method will both become nearly negligible.

4.1.2. Breakthrough curves (BTCs) in the aquifer Fig. 3 shows the comparison of BTCs in the aquifer calculated by the Zhan method at zD = 0 and zD = 0.5 under conditions of β1 = β2 = 0.1, σ1 = σ2 = 1, Pe = 10, Pe1 = Pe2 = 1, and xD = 10. The result from the AA method, which is independent of the z-axis, is also shown in this Figure. Several observations can be made. First, the concentration at zD = 0 is always greater than that at zD = 0.5. This simply reflects the fact that a concentration gradient exists to drive mass upward into the upper aquitard. Second, the AA method overestimates the aquifer concentration at early times but underestimates the aquifer concentration at late times. Third, even under relatively long times, aquifer concentrations cannot reach the source concentration C0 (or CD

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

169

than the source concentration. This phenomenon is somewhat similar to that observed for a radioactively decaying solute. 4.1.3. Comparison of the semi-analytical solutions and the numerical simulations To test the semi-analytical solutions obtained here, we have conducted numerical simulations using the model MT3DMS developed by Zheng and Bennett (2002). The numerical simulations are straightforward, thus are only briefly outlined here. A two-dimensional aquifer–aquitard system is generated and the parameters used are the same as those outlined in Section 4. The idea is to reproduce the situation of β1 = β2 = 0.1, σ1 = σ2 = 1, Pe = 10, and Pe1 = Pe2 = 1. The simulated domain includes an aquifer that is 4 m thick and 1000 m long, and an upper and a lower aquitard. The thickness of each aquitard is 50 m, which is sufficiently large to satisfy the boundary conditions of Eqs. (8) and (13). The hydraulic conductivity of the aquitards is approximately three orders of magnitude lower than that in the aquifer, so that flow is essentially vertical in the aquitards. The vertical flow velocity in the aquitards is 1/100 of the horizontal flow velocity in the aquifer. The flow velocities in the aquifer and aquitards are obtained through Modflow after appropriately specifying the boundary condition of flow and these velocities are used by MT3DMS to generate the concentration field (Zheng and Bennett, 2002). The grid spacing was made sufficiently fine to reduce the numerical errors. After trying different grid spacings, we decided to use a spacing of 0.2 m for the aquifer and the aquitards. Further refining the grid space to 0.1 m did not produce noticeable changes to the output. Different solver options of MT3DMS had been tried except the explicit finite difference method, and similar results were obtained. The hybrid Method of Characteristic (HMOM) was eventually chosen because of its relative short computational time. A constant–concentration boundary condition was used on the left-hand side of the aquifer in the numerical simulation. However, one must be careful with the nature of comparison of the analytical and numerical modeling results.

Fig. 2. The vertical concentration profiles in the aquifer–aquitard system under conditions of σ1 = σ2 = 1, Pe = 10, xD = 10, and tD = 2. (A) β1 = 0.05, β2 = 0.1, Pe1 = 0.5, and Pe2 = 1; (B) β1 =β2 = 0.1, and Pe1 =Pe2 = 1; (C) β1 = β2 = 0.1, Pe1 = 1, and Pe2 =− 1.

of 1). That is because in addition to the horizontal mass flux in the aquifer, the vertical mass flux across the aquifer–aquitard boundary serves as a sink to keep aquifer concentration lower

Fig. 3. The breakthrough curves (BTC) at zD = 0, and zD = 0.5 under conditions of β1 = β2 = 0.1, σ1 = σ2 = 1, Pe = 10, Pe1 = Pe2 = 1, and xD = 10 from the Zhan method. The result of the AA method, which is independent of the z axis, is also shown in this figure.

170

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

Such a comparison is by no means to be exact in a quantitative sense because of a few reasons. First, diffusion along the xaxis is not considered in the governing equations of transport in the aquitard (see Eqs. (5) and (10)), but it is included in MT3DMS. In this sense, the numerical modeling seems to have an advantage over the analytical modeling of this study. Unfortunately, our numerical exercises indicate that MT3DMS performs poorly near the aquifer–aquitard interfaces. Such poor performance is caused by the abrupt change of aquifer parameters and concentration gradients near the aquifer– aquitard interfaces. This issue cannot be simply resolved by refining the grid blocks near those interfaces. New numerical modeling strategies should be used to accurately model the transport behavior near the aquifer–aquitard interfaces. The finding of our numerical exercises seems to agree with previous numerical simulation which indicated that the numerical errors tended to be the largest at the aquifer– aquitard interfaces (Martin and Frind, 1998; Bester et al., 2005). In this regard, the analytical modeling seems to have some advantages over the numerical modeling. The idea of comparing the analytical and numerical modeling results is not to show which one is better; rather, it is to show how different are those results in a qualitative sense. The default parameters used in the comparison come from Table 2. As an example, Fig. 4 shows the concentration distribution in the system obtained from the numerical simulation against that obtained from the Zhan method at xD = 2.5 (or x = 50 m) at tD = 2 or (t = 3991 day). The advections in both aquitards are upward thus Fig. 4 is for an asymmetric case. The horizontal flow velocity in the aquifer is 0.1 m/day. This implies that the averaged travel distance of the solute is x = 3991 × 0.1 m = 399.1 m, which is considerably larger than the coordinate of the observation point in Fig. 4 (which is x = 50 m). Thus the result in Fig. 4 can be regarded as the one far behind the front of the solute. As can be seen from this figure, the result of the numerical simulation is close to that of the Zhan method

in the lower aquitard but it differs considerably from that of the Zhan method in the upper aquitard and in the upper portion of the aquifer. Large differences between the analytical and numerical modeling results exist at both the upper and lower aquifer–aquitard interfaces. The asymmetric shapes of the concentration distribution in Fig. 4 are apparently due to the asymmetric aquitard advections which are upward in both aquitards here. 4.1.4. Mass transported into the aquitard It is interesting to know how much mass was stored in the aquitard due to advective and diffusive transport across the aquifer–aquitard boundary. The total masses per unit width along the direction perpendicular to the xz plane in the upper and lower aquitards are denoted as M1 and M2, respectively. Their dimensionless forms, denoted as M1D and M2D are defined in Table 1. The Laplace transforms of M1D and M2D, denoted as M ̅1D and M ̅2D respectively, can be obtained through the following integration: Z M1D =



0

Z dxD



1

Z C 1D dzD ;

M2D =

∞ 0

Z dxD

−1 −∞

C 2D dzD ð41Þ

This can be done straightforwardly. Using Eqs. (28) and (29), one has M 1D =

∞ X n=0

M2D =

∞ X n=0

4An cosðωn + μ n Þ ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   2 2 Pe + 4 p + ωn − Pe ðPe1 =β 1 Þ2 + 4p = β 1 − Pe1 = β 1 4An cosðωn + μ n Þ ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   Pe2 + 4 p + ω2n − Pe ðPe2 = β 2 Þ2 + 4p = β 2 + Pe2 = β2

ð42Þ where An is given by Eq. (30) for the first-type boundary condition at x = 0 or Eq. (32) for the third-type boundary ̅ and condition at x = 0. Inverse Laplace transforms of M 1D ̅ in Eq. (42) will result in the solution in the real time M 2D domain for the Zhan method. If using Eqs. (39) and (40) for the AA method, one can ̅ as: ̅ and M 2D obtain M 1D 4= p M 1D = sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Pe2 + 4p + σ 1

Pe21 + 4pβ1 + Pe1

+ σ2

Pe22 + 4pβ2 − Pe2

− Pe

1 × qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðPe1 = β1 Þ2 + 4p = β1 − Pe1 = β1 4= p M 2D = sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 2 Pe1 + 4pβ1 + Pe1 + σ 2 Pe22 + 4pβ2 − Pe2 − Pe Pe + 4p + σ 1 1 × qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ðPe2 = β2 Þ2 + 4p = β 2 + Pe2 = β2

Fig. 4. Comparison of the vertical concentration profiles in the aquifer– aquitard system obtained from the numerical simulation using MT3DMS and the Zhan method. β1 = β2 = 0.1, σ1 = σ2 = 1, Pe = 10, Pe1 = Pe2 = 1, tD = 2. The profiles at xD = 2.5 is shown here. zD = 1 and − 1 are the upper and lower aquifer–aquitard interfaces. The first-type solute transport boundary condition is used on the left hand side of the aquifer.

ð43Þ

̅ in Eq. (43) Inverse Laplace transforms of M ̅1D and M 2D will result in the solution in the real time domain for the AA method. As an example, a comparison of the M1D values calculated from Eq. (42) (the solid line) and Eq. (43) (the dashed line) is shown in Fig. 5 for β1 = β2 = 0.05, σ1 = σ2 = 1, Pe = 10, Pe1 = Pe2 = 0.1 and 1. Several observations can be made from Fig. 5.

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

171

Fig. 5. The dimensionless mass transported to the aquitard per unit width (M1D) calculated from the Zhan method and that from the AA method versus the dimensionless time (tD) under the conditions of β1 = β2 = 0.05, σ1 = σ2 = 1, Pe = 10, Pe1 = Pe2 = 0.1 and 1.

Fig. 6. The dimensionless mass transported to the aquitard per unit width (M1D) calculated from the Zhan method and that from the AA method versus the dimensionless time (tD) under the conditions of β1 = β2 = 0.1, σ1 = σ2 = 1, Pe = 10, Pe1 = Pe2 = 0.1 and 1.

First, the AA method always overestimates the mass transported to the aquitard, and the error increases with the Peclet numbers. For example, at tD = 2 and Pe1 = Pe2 = 0.1, the error is about 18%, while at tD = 2 and Pe1 = Pe2 = 1, the error is about 144%. This is understandable. In the AA method, the aquitard diffusion is regarded as a volumetric source/sink term in the governing equation of transport in the aquifer, thus every point from zD = 0 to zD = 1 within the aquifer is capable of losing mass, resulting in a greater loss of mass to the aquitard. This is contradictory to the realistic mass transfer that only occurs at the aquifer–aquitard interface. The Zhan method treats the mass flux to the aquitard as an interface term rather than a volumetric term, thus it results in less total mass transported. Second, the value of the aquitard Peclet number greatly influences the M1D values and a greater Peclet number results in a greater M1D value when time is sufficiently large (tD N 1). The reason is quite straightforward since a greater Peclet number indicates greater upward advective transport from the aquifer to the upper aquitard, thus more mass ends up in the upper aquitard. Third, the functionality of M1D versus tD is nonlinear. In fact, the slope of M1D versus tD increases with tD when tD is relatively small and eventually approach a constant when tD is close to 2, indicating an advective transport mechanism at the late time. A tD of 2 is equivalent to a real time of t = 4000 days or close to 11 years (see Tables 1 and 2). Therefore, both diffusion and advection will affect the earlytime mass transport, but only advection will affect the late time transport between the aquifer and the aquitard. Fig. 6 is similar to Fig. 5 except that both the β1 and β2 values are increased to 0.1. It is interesting to see that when the β1 and β1 values are doubled, indicating that the aquitard diffusion becomes stronger, the mass transported to the upper aquitard increases as well, as expected. However, the percentage of increase is greater for a case with a smaller aquitard Peclet number. For instance, at tD of 2 and Pe1 = Pe2 = 0.1, the M1D value increases 20% from Figs. 5 and 6, while at tD of 2 and Pe1 = Pe2 = 1, the M1D value increases

only 0.4% from Figs. 5 and 6. This finding implies that the mass transported between the aquifer and the aquitard is relatively insensitive to the aquitard diffusion when the aquitard Peclet number is greater than unity. Fig. 7 is similar to Figs. 5 and 6 except that the upper and lower aquitards have different diffusion coefficients. In this case, the upper aquitard diffusion coefficient is one half of the lower one. One can tell that by reducing the upper aquitard diffusion, the M1D value is slightly reduced. But in general, the M1D value is less sensitive to the difference of the aquitard diffusion coefficients. A great discrepancy has been observed between the Zhan method and the AA method when Pe1 = Pe2 = 1. Fig. 8 shows the results for different aquitard Peclet numbers. In this example, the Peclet number of the upper aquitard is 1/10 of that of the lower aquitard. As can be seen obviously, the M1D value is sensitive to the difference of the

Fig. 7. The dimensionless mass transported to the aquitard per unit width (M1D) calculated from the Zhan method and that from the AA method versus the dimensionless time (tD) under the conditions of β1 = 0.05, β2 = 0.1, σ1 = σ2 = 1, Pe = 10, Pe1 = Pe2 = 0.1 and 1.

172

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

Peclet numbers between the upper and lower aquitards, particularly at late times.

generally overestimates the concentration in the upper aquitard for this asymmetric case. When flow velocities in the upper and lower aquitards have the same magnitude but opposite directions, a symmetric concentration distribution along the vertical direction will be obtained. For this case, the concentrations in the aquifer obtained by the two methods differ considerably. However, much less discrepancy has been observed in the aquitards. Comparison of the numerical simulation using MT3DMS and the Zhan method indicates that the result of the numerical simulation could differ considerably from that of the Zhan method for an asymmetric case where aquitard advection is in the same direction. More specifically, in the aquitard where the advection is towards the aquifer–aquitard interface, the numerical result agrees with the result of the Zhan method reasonably well; in the aquitard where the advection is away from the aquifer–aquitard interface, the numerical result differs considerably from the result of the Zhan method. Compared to the Zhan method, the AA method overestimates mass transported into the upper aquitard when an upward advection exists. The mass transported between the aquifer and the aquitard is sensitive to the aquitard Peclet number, but less sensitive to the aquitard diffusion coefficient, particularly at late times.

5. Summary and conclusions

Acknowledgements

This study deals with two-dimensional solute transport in an aquifer–aquitard system by maintaining rigorous mass conservation at the aquifer–aquitard interface. Advection, longitudinal dispersion, and transverse vertical dispersion are considered in the aquifer. Vertical advection and diffusion are considered in the aquitards. The first-type and the third-type boundary conditions are considered in the aquifer. This study differs from previous AA method that treats the mass flux between the aquifer and aquitard as an averaged volumetric source/sink term in the governing equation of transport in the aquifer. Instead, it solves the transport governing equations in the aquifer and the aquitards explicitly by maintaining continuity of concentration and vertical mass flux at the aquifer–aquitard interfaces. Analytical solutions of concentrations in the aquitards and aquifers are obtained in the Laplace domain and are inverted numerically to yield concentrations in the real time domain. The following conclusions can be drawn from this study. The BTCs and concentration profiles obtained in the Zhan method are drastically different from those obtained using the AA method. The AA method shows a sharp change of concentration profile near the aquifer–aquitard interface whereas the Zhan method shows a continuously decreasing concentration near that interface. Concentration distributions in the aquitards are closely related to advection in the aquitards. When flow velocities in the upper and lower aquitards have the same magnitude and both are along the upward direction, an asymmetric concentration distribution along the vertical direction will result. The AA method overestimates the aquifer concentration at early times but underestimates the aquifer concentration at late times. However, a greater discrepancy exists between the result of the AA method and the Zhan method in the upper aquitard than that in the lower aquitard. The AA method

This study is partially supported by National Science Foundation of China (#50428907), the CONACYT-TAMU Collaborative Research Program, and the Advanced Research Program (ARP) from the Texas Higher Education Coordination Board, USA. We thank two anonymous reviewers for their constructive comments. In particular, we thank the Editor-inChief E. O. Frind for his critical and detailed suggestions of this manuscript. Dr. Frind's comments have helped us improve this manuscript substantially.

Fig. 8. The dimensionless mass transported to the aquitard per unit width (M1D) calculated from the Zhan method and that from the AA method versus the dimensionless time (tD) under the conditions of β1 = β2 = 0.05, and β1 = β2 = 0.1, σ1 = σ2 = 1, Pe = 10, Pe1 = 0.1 and Pe2 = 1.

Appendix A. Solute transport in an aquifer bounded by two aquitards The following solution is proposed for C D̅ : CD =

∞ X

An C xD ðxD ; p; nÞ cosðωn zD + μ n Þ;

ðA1Þ

n=0

̅ where An is the coefficient that needs to be determined, C xD (xD, p, n) is the part related to the x coordinate, and ωn is the frequency of the Fourier transform along the z-axis. Substituting Eq. (A1) into Eq. (16) results in the following equa̅ (xD, p, n): tion for C xD d2 C xD PedC xD 2 − − p + ωn C xD = 0: 2 dxD dxD

ðA2Þ

Considering the boundary condition of Eq. (18) will yield the general solution of Eq. (A2): 2 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 3   Pe2 + 4 p + ω2n − Pe AxD 5: ðA3Þ C xD ðxD ; p; nÞ = exp4−@ 2

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

Therefore, the proposed solution for the aquifer becomes: 2 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 3   ∞ Pe2 + 4 p + ω2n − Pe X 4 @ AxD 5 cosðωn zD + μ Þ; CD = An exp − n 2 n=0 − V zD V 1

ðA4Þ

The solution to Eq. (19) with the boundary conditions of Eqs. (20) and (22) is: 2 C 1D = C D ðxD ; zD = 1; pÞexp4−

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðPe1 =β 1 Þ2 + 4p = β1 − Pe1 = β1 2

zD z 1:

3 ðzD − 1Þ5;

ðA5Þ

singularities of ωn =π/4+nπ/2 separate the entire domain into many sections with an equal interval of π/2. It is worthwhile to pffiffiffi point out that near the singular point of ωn = b, there are alwaysptwo roots of ωn within that π/2 section including ffiffiffi ωn = b. Only one root of ωn exists within each of the rest π/2 sections. After determining ωn and μn, one can determine An from the boundary condition Eq. (17): ∞ X

An cosðωn zD + μ n Þ = 1 = p:

C 2D

zD V − 1:

ðA6Þ

Substituting Eq. (A4) into Eqs. (21) and (25), and considering Eqs. (A5) and (A6) will lead to: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ωn tanðωn + μ n Þ = σ 1 Pe21 = 4 + pβ1 + σ 1 Pe1 = 2; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ωn tanðωn − μ n Þ = σ 2 Pe22 = 4 + pβ2 − σ 2 Pe2 = 2:

ðA7Þ

Conducting inverse Fourier transform of Eq. (A10) results in:

×

2 n

− σ 1σ 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Pe21 = 4 + pβ1 + Pe1 = 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Pe22 = 4 + pβ2 − Pe2 = 2

= ωn

½

ðA8Þ

Þ tanð2ω Þ n

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  σ 1 Pe21 = 4 + pβ1 + σ 1 Pe1 = 2



 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  + σ 2 Pe22 = 4 + pβ2 − σ 2 Pe2 = 2 ;

ð

2

ωn + σ 1 σ 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Pe21 = 4 + pβ1 + Pe1 = 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  × Pe22 = 4 + pβ2 − Pe2 = 2 = ωn

½

An =

4 F ðω ; μ Þ × 1 n n ; p F2 ðωn ; μ n Þ

ðA11Þ

where F1 ðωn ; μ n Þ = sinðωn Þ cosðμ n Þ;

ðA12Þ

F2 ðωn ; μ n Þ = 2ωn + sinð2ωn Þ cosð2μ n Þ: ̅ , and C 2D ̅ are derived after substituting Therefore, C D̅ , C 1D Eq. (A11) into Eqs. (A4), (A5), and (A6), respectively. For the special case in which the upper and lower aquitards have identical parameters, μn = 0. References

From Eq. (A7), one has:

ðω

ðA10Þ

n=0

Similarly, 2qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðPe2 = β2 Þ2 + 4p = β2 + Pe2 = β 2 = C D ðxD ; zD = − 1; pÞexp4 ðzD + 1Þ5; 2

173

ðA9Þ

Þ tanð2μ Þ n

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  σ 1 Pe21 = 4 + pβ1 + σ 1 Pe1 = 2



 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  − σ 2 Pe22 = 4 + pβ2 − σ 2 Pe2 = 2 ; ωn is determined from Eq. (A8). After that, μn is calculated from Eq. (A9). For the numerical calculation, accurate determination of ωn and μn are important since these parameters have been used in solutions (see Eqs. (27)–(33)). As can be seen from Eq. (A8), ωn is the root of a function tan(2ωn)=αωn/(ωn2 −b), where a and b are two dummy variables associated with all the other 2 parameters of Eq. (A8). pffiffiffi Notice that function αωn/(ωn −b)has a singularity at ωn = b, and the function of tan(2ωn) has infinite numbers of singularities at ωn =π/4+nπ/2, where n is an integer varying from negative infinity to positive infinity. The

Ball, W.P., Liu, C., Xia, G., Young, D.F., 1997a. A diffusion-based interpretation of tetrachloroethene and trichloroethene concentration profiles in a groundwater aquitard. Water Resour. Res. 33 (12), 2741–2757. Ball, W.P., Xia, G., Durfee, D.P., Wilson, R.D., Brown, M.J., Mackay, D.M., 1997b. Hot methanol extraction for the analysis of volatile organic chemicals in subsurface core samples from Dover Air Force Base, Delaware. Ground Water Monit. Remediat. 17 (1), 104–121. Batu, V., 1996. A generalized three-dimensional analytical solute transport model for multiple rectangular first-type sources. J. Hydrol. 174, 57–82. Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York, USA. Bester, M., Frind, E.O., Molson, J.W., Rudolph, D.L., 2005. Numerical investigation of road salt impact on an urban well field. Ground Water 44 (2), 165–175. Chen, C.S., 1985. Analytical and approximate solutions to radial dispersion from an injection well to a geological unit with simultaneous diffusion into adjacent strata. Water Resour. Res. 21 (8), 1069–1076. de Hoog, F.R., Knight, J.H., Stokes, A.N.,1982. An improved method for numerical inversion of Laplace transforms. SIAM J. Sci. Statist. Comput. 3 (3), 357–366. Domenico, P.A., Schwartz, F.W., 1998. Physical and Chemical Hydrogeology, 2nd ed. Wiley, New York, USA. Fetter, C.W., 1999. Contaminant Hydrogeology, 2nd edition. Prentice-Hall, Upper Saddle River, NJ, USA. Fujikawa, Y., Fukui, M., 1990. Adsorptive solute transport in fractured rock: analytical solutions for delta-type source conditions. J. Contam. Hydrol. 6, 85–102. Gillham, R.W., Sudicky, E.A., Cherry, J.A., Frind, E.O., 1984. An advective diffusive concept for solute transport in heterogeneous unconsolidated geological deposits. Water Resour. Res. 20 (3), 369–378. Hendry, M.J., Ranville, J.R., Boldt-Leppin, B.E.J., Wassenaar, L.I., 2003. Geochemical and transport properties of dissolved organic carbon in a clay-rich aquitard. Water Resour. Res. 39 (7), 1194. doi:10.1029/2002WR001943. Hunkeler, D., Chollet, N., Pittet, X., Aravena, R., Cherry, J.A., Parker, B.L., 2004. Effect of source variability and transport processes on carbon isotope ratios of TCE and PCE in two sandy aquifers. J. Contam. Hydrol. 74, 265–282. Huyakorn, P.S., Ungs, M.J., Mulkey, L.A., Sudicky, E.A., 1987. A threedimensional analytical method for predicting leachate migration. Ground Water 25 (5), 588–598. Johns, R.A., Roberts, P.V., 1991. A solute transport model for channelized flow in a fracture. Water Resour. Res. 27 (8), 1797–1808.

174

H. Zhan et al. / Journal of Contaminant Hydrology 107 (2009) 162–174

Johnson, R.L., Cherry, J.A., Pankow, J.F., 1989. Diffusive contaminant transport in natural clay: a field example and implications for clay-lined waste disposal sites. Environ. Sci. Technol. 23, 340–349. Kreft, A., Zuber, A., 1978. On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions. Chem. Eng. Sci. 33, 1471–1480. Kreft, A., Zuber, A., 1979. On the use of the dispersion model of fluid flow. Int. J. Appl. Radiat. Isot. 30, 705–708. Leij, F.J., Skaggs, T.H., van Genuchten, M.Th., 1991. Analytical solutions for solute transport in three-dimensional semi-infinite porous media. Water Resour. Res. 27 (10), 2719–2733. Liu, C.X., Ball, W.P., 1999. Application of inverse methods to contaminant source identification from aquitard diffusion profiles at Dover AFB, Delaware. Water Resour. Res. 35 (7), 1975–1985. Liu, C.X., Ball, W.P., 2002. Back diffusion of chlorinated solvent contaminants from a natural aquitard to a remediated aquifer under well-controlled field conditions: predictions and measurements. Ground Water 40 (2), 175–184. Liu, H.H., Bodvarsson, G.S., Zhang, G., 2004. The scale-dependency of the effective matrix diffusion coefficient. Vadose Zone J. 3, 312–315. Martin, P.J., Frind, E.O., 1998. Modeling a complex multi-aquifer system: the Waterloo moraine. Ground Water 36 (4), 679–690. Moench, A.F., 1991. Convergent radial dispersion: a note on evaluation of the Laplace transform solution. Water Resour. Res. 27 (12), 3261–3264. Moench, A.F., Ogata, A., 1981. A numerical inversion of the Laplace transform solution to radial dispersion in a porous medium. Water Resour. Res. 17 (1), 250–252. Moreno, L., Neretnieks, I., Eriksen, T., 1985. Analysis of some laboratory tracer runs in natural fissures. Water Resour. Res. 21 (7), 951–958. Neretnieks, I., 1980. Diffusion in the rock matrix: an important factor in radionuclide retardation? J. Geophys. Res. 85 (B8), 4379–4397. Neretnieks, I., Eriksen, T., Tähtinen, P., 1982. Tracer movement in a single fissure in granitic rock: some experimental results and their interpretation. Water Resour. Res. 18 (4), 849–858. Ogata, A., Banks, R.B., 1961. A solution of differential equation of longitudinal dispersion in porous media. U.S. Geol. Surv. Prof. Pap. 411, A1–A7. Parker, B.L., Cherry, J.A., Chapman, S.W., 2004. Field study of TCE diffusion profiles below DNAPL to assess aquitard integrity. J. Contam. Hydrol. 74, 197–230. Rasmuson, A., Neretnieks, I., 1981. Migration of radionuclides in fissured rock: the influence of micropore diffusion and longitudinal dispersion. J. Geophys. Res. 86 (B5), 3749–3758.

Stehfest, H., 1970. Numerical inversion of Laplace transforms. Commun. ACM 13 (1), 47–49. Sudicky, E.A., 1986. A natural gradient experiment on solute transport in a sand aquifer: spatial variability of hydraulic conductivity and its role in the dispersion process. Water Resour. Res. 22 (13), 2069–2082. Sudicky, E.A., 1989. The Laplace transform Galerkin technique: a timecontinuous finite element theory and application to mass transport in groundwater. Water Resour. Res. 25 (8), 1833–1846. Sudicky, E.A., Frind, E.O., 1982. Contaminant transport in fractured porous media: analytical solutions for a system of parallel fractures. Water Resour. Res. 18 (6), 1634–1642. Sudicky, E.A., Frind, E.O., 1984. Contaminant transport in fractured porous media: analytical solutions for a two-member decay chain in a single fracture. Water Resour. Res. 20 (7), 1021–1029. Sudicky, E.A., Gillham, R.W., Frind, E.O., 1985. Experimental investigation of solute transport in stratified porous media, 1. The nonreactive case. Water Resour. Res. 21 (7), 1035–1041. Starr, R.C., Gillham, R.W., Sudicky, E.A., 1985. Experimental investigation of solute transport in stratified porous media, 2. The reactive case. Water Resour. Res. 21 (7), 1043–1050. Talbot, A., 1979. The accurate numerical inversion of the Laplace transforms. J. Inst. Math. Appl. 23, 97–120. Tang, Y., Aral, M.M., 1992a. Contaminant transport in layered porous media, 1. General solution. Water Resour. Res. 28 (5), 1389–1397. Tang, Y., Aral, M.M., 1992b. Contaminant transport in layered porous media, 2. Applications. Water Resour. Res. 28 (5), 1399–1406. Tang, D.H., Frind, E.O., Sudicky, E.A., 1981. Contaminant transport in fractured porous media: analytical solution for a single fracture. Water Resour. Res. 17 (3), 555–564. van Genuchten, M., 1981. Analytical solutions for chemical transport with simultaneous adsorption, zero-order production and first-order decay. J. Hydrol. 49, 213–233. van Genuchten, M.Th., Parker, J.C., 1984. Boundary conditions for displacement experiments through short laboratory soil. Soil Sci. Soc. Am. J. 48, 703–708. Young, D.F., Ball, W.P., 1998. Estimating diffusion coefficients in lowpermeability porous media using a macropore column. Environ. Sci. Technol. 32 (17), 2578–2584. Zheng, C., Bennett, G.D., 2002. Applied Contaminant Transport Modeling, 2nd ed. John Wiley & Sons, New York, USA.