Chapter 23 Modeling Contaminant Transport 1 Chapter 23 Modeling Contaminant Transport Chapter Contents 23.1 Derivation of the Mass Transport Equations 23.2 Advection-Dispersion-Reaction Equation 23.3 Boundary and Initial conditions 23.4 Analytical Solutions of Mass Transport Equations in One Dimension 23.5 Analytical Solutions of Mass Transport Equations in Three Dimensions 23.6 Other Useful Analytical Solutions 23.7 Particle Tracking Methods 23.8 Numerical Solutions of Mass Transport Equations 23.9 Introduction to CONTz.exe Mass transport models play a vitally important role in the analysis of problems of contamination. The main application is in forward modeling - predicting how contaminant plumes spread away from a source as a function of time. This approach is used in practice * to explain how a plume of contamination evolved and how it will behave in the future, * to estimate mass transport parameters (for example, dispersivity values or Kd values) through a trial-and-error process that involves fitting the model to some observed concentration distribution, and * to design and to evaluate the likely behavior of systems for contaminant remediation. The modeling process is similar to that developed for ground-water flow and well hydraulics. We start with a partial differential equation that represents the mass transport and mass transfer processes, such as advection, dispersion, and kinetic decay, and the dimensionality of the system. Accordingly, in the first section of this chapter, we develop various mass transport equations and solve them subject to boundary and initial conditions in order to provide concentration distributions as a function of space and time. Chapter 23 Modeling Contaminant Transport 2 As was the case with well hydraulics, two main approaches may be used to solve mass transport equations. The analytical approach is based on classical mathematical procedures for solving differential equations, which in many cases do not require a computer. They are similar in this respect to the Theis solution. Not unexpectedly, the analytical models are limited to simple mass transport problems with constant parameter values. Similarly, they are usually restricted to modeling the transport from only one to perhaps several dissolved constituents. This chapter presents several analytic solutions that provide basic tools for looking at contamination problems. Sometimes the objective of a study is to represent a complex hydrogeological setting, coupled reactions among constituents, and heterogeneities in system parameters rigorously. In this case, numerical models provide the appropriately powerful modeling approach. Industry standard codes like MT3D are part of the MODFLOW family and find wide use in the evaluation of practical hydrogeologic problems. In this chapter, we touch on numerical models with an application that involves the simulation of a large plume on Cape Cod, Massachusetts. 23.1 Derivation of the Mass Transport Equations This section develops a general differential equation, describing the transport of a dissolved constituent, subject to key physical and chemical transport processes. As was the case with flow, conservation principles, here solute mass conservation, are applied in relation to a control volume (Figure 23.1). In words, solute is conserved in the system such that Solute inflow rate ­ Solute outflow rate = Change of solute storage with time (23.1) Chapter 23 Modeling Contaminant Transport 3 Figure 23.1 Mass flux balance in a control volume: Net Flux In ­ Net Flux Out = Storage Change. The derivation of the differential equation proceeds by determining the mass inflow rate across the face ABCD (x = x - x/2) and the mass outflow rate across the face EFGH (x = x + x/2). Taylor's expansion series of Jx at the outflow face EFGH (x= x + x/2) is given by Thomas (1972) as L+) 2 x ( x J +J=J xx x xxXxxxX ==+= ||| 2/ (23.2) Similarly, the Taylor's expansion series at the inflow face ABCD (x = x - x/2) is given by L+) 2 x ( x J J=J xx x xxXxxxX - ==-= ||| 2/ (23.3) Chapter 23 Modeling Contaminant Transport 4 The solute mass flux out of the control volume in the x-direction is expressed as ( ) zyx x J zyJJ= xx x xxxxxxxxx -=- = +=-= 2/2/ |||rateinflownetSolute (23.4) Similarly, the net inflow rates for solutes in y- and z-directions in the control volume are zyx y J = yy y y - = |rateinflownetSolute (23.5) zyx z J = zz z z - = |rateinflownetSolute (23.6) The solute mass change rate in the control volume is zyx t nC = )( ratechangestorageSolute (23.7) Inserting Eqs. (23.4) through (23.7) into Eq. (23.1), we get t nC z J y J x J zz z yy y xx x = - - - === )( (23.8) Inserting Eqs. (19.20) through (19.22) into Eq. (23.8), we obtain + + z C nD xy C nD xx C nD x xzxyx 22 + + + z C nD zz C nD yy C nD y zzyzy 2 t nC z Cnv y Cnv x Cnv izyx = - - - )()()()( (23.9) Eq. (23.9) can be written in an equivalent but more compact form as t nC nvCn =- )( )()D( (23.10) Chapter 23 Modeling Contaminant Transport 5 Eq. (23.10) states that the contaminant net inflow rate in a control volume is equal to the contaminant mass storage change rate. 23.2 Advection-Dispersion-Reaction Equation The chemical reactions of a dissolved species (often contaminants) can be incorporated into Eq. (23.10) by adding a reaction term on the left side of the equation t nC rnvCCn =+- )( )()D( (23.11) where r is the reaction rate defined as mass produced from chemical reactions in a unit volume per unit time. For a sorbing species, another term, accounting for the change in mass storage, needs to be added t s t nC rCnCn b + =+- )()( )v()D( (23.12) Assuming that sorption is modeled as a linear equilibrium, reversible reaction, Eq. (23.12) is rewritten as t nC RrCnCn f =+- )( )v()D( (23.13) where n K R db f += 1 (23.14) where Rf is the retardation factor, b is the bulk density of the soil or rock, and Kd is the distribution coefficient of a species in the soil or rock. Additional chemical reactions can be incorporated via the reaction term, r. For a decaying species, described by a first-order kinetic reaction, Eq. (23.13) is written as Chapter 23 Modeling Contaminant Transport 6 t nC RCnRCnCn ff =-- )( )v()D( (23.15) where is the decay rate constant. Recall that the decay rate is related to the half-life of the reaction by )2(ln 2/1 =t (23.16) With additional modifications, we can derive other useful forms of the transport equations. For example, with velocity |v| = 0, Eq. (23.15) becomes a diffusion equation. t nC RCnRCn ff =- )( )D( (23.17) The solution of Eq. (23.17) is similar to that of Eq. (5.14) or (5.22), which is called a diffusion equation. It can be applied to a variety of practical problems. For example, the transport of a volatile organic contaminant source in unsaturated zones, vapor concentration caused by the source may be treated as a diffusion problem. For a unidirectional flow (vx = vL, vy = 0, and vz = 0) in a homogeneous medium, the dispersion tensor D can be simplified, as shown by Eqs. (19.23) through (19.26). mxLLx DvDD +== (23.18) mxTHTHy DvDD +== (23.19) mxTVTVz DvDD +== (23.20) 0====== zyyzzxxzyxxy DDDDDD (23.21) Thus, Eq. (23.17), in turn, reduces to t nC RCR x C v z C D y C D x C D ffxzyx =- - + + )( 2 2 2 2 2 2 (23.22) Chapter 23 Modeling Contaminant Transport 7 23.3 Boundary and Initial conditions Three types of boundary conditions specify the mass exchange between a simulation domain and its surrounding regions: First-type (or Dirichlet-type) boundary conditions prescribe the concentration along a boundary segment (Figure 23.2). ),,,(),,,( 0 tzyxctzyxc = for (x, y, z)1 (23.23) This boundary condition is analogous to the constant-head boundary conditions that were developed for problems of flow. Figure 23.2. Possible boundary conditions in mass transport problems. Third-type (or Cauchy-type) boundary conditions prescribe the solute flux along a boundary segment (Figure 23.2). 0cnqcnqn x c nD iiiii j ij =+ - for (x, y, z)3 (23.24) Second-type (or Neumann-type) boundary conditions describe a situation in which there is no dispersive or diffusive flux across the boundary (Figure 23.2). Chapter 23 Modeling Contaminant Transport 8 0= - i j ij n x c nD for (x, y, z)2 (23.25) The initial condition for a mass transport problem specifies the initial concentration of dissolved mass within the simulation domain and is expressed as ),,()0,,,( zyxczyxc init= 23.26) where Cinit is the concentration distribution in the simulation domain at the beginning of the simulation (Figure 23.2). 23.4 Analytical Solutions of Mass Transport Equations in One Dimension This section introduces a few basic analytical solutions for one-dimensional mass transport problems. The mass transport equation in one-dimensional space is written as CR x C v x C D t C R fxxf - - = 2 2 (23.27) where vx is the linear ground-water velocity in the x direction, Rf is the retardation factor, and Dx is the dispersion coefficient. Instantaneous Point Source Figure 23.3. A simple one-dimensional problem of mass transport from a point source. Chapter 23 Modeling Contaminant Transport 9 The first case is one-dimensional transport due to an instantaneous point source (Figure 23.3). The boundary conditions for mass transport in one dimension with an instantaneous point source (Figure 23.3) are expressed as 0= x C , at x = 0 and (23.28) with the initial condition (at t = 0): 0MM = , at x = 0 (23.29) 0=C , at 0 < x < (23.30) The mass M0 ( = C0V0) represents the product of the source concentration C0 and the source volume V0. The solution for mass transport in one dimension with an instantaneous point source is given by ),(),( 1 0 txX AnR M txC f = (23.31) t RtD Rtvx fx fx fx e RtD txX - - - = /4 )/( 1 2 /4 1 ),( (23.32) where M0 is the contaminant mass, and A is the cross-sectional area of the one-dimensional model. Example 23.1: The contaminant mass at the source is 10 g at x = 0 cm in a column with a crosssectional area of 1 cm2 . The longitudinal dispersivity is 1 cm. A water flux of 0.35 cm3 /min is maintained continuously at x = 0. The porosity of the medium is 0.35, and the molecular diffusion coefficient of the contaminant is 1.0 x 10-5 cm2 /min. Assuming an instantaneous source, calculate the contaminant concentration profiles along the column at 100 minutes and the breakthrough curve at x = 25 cm for decay constants = 0 and 0.01 min-1 . Solution: The flow velocity is Chapter 23 Modeling Contaminant Transport 10 min/0.1 )1(35.0 min/35.0 2 3 cm cm cm nA Q vx === The tortuosity of the porous medium is calculated according to Eq. (19.7) 7.035.0 3/1 == The dispersion coefficient is calculated by Eq. (19.19) min/0.1min)/101(7.0)1min)(/0.1( 225 cmcmcmcmDx =×+= The concentration is calculated using the Microsoft Excel file "analytic1.xls" located in directory Chap23 on the accompanying CD. There are two calculation sheets in the file "analyitic1.xls". The calculation sheet "1d-point-source" applies to this example. Concentration profiles along the column are shown in Figure 23.4 for the different decay constants. Notice how decay reduces the total amount of mass in the system. The breakthrough curves for these two cases are shown in Figure 23.5 at 25 cm. Concentration along the flow direction 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0 50 100 150 200 x (cm) Concentration Figure 23.4. Concentration profiles in the direction of flow at 100 minutes from a point source. Chapter 23 Modeling Contaminant Transport 11 Breakthrough curve 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0 10 20 30 40 50 60 Time (min) Concentration Figure 23.5. Breakthrough curve at x = 25 cm for a point source. Semi-infinite Systems with Continuous First-Type Boundary Conditions at the Source Figure 23.6. A problem of one-dimensional mass transport with a continuous source at constant concentration. Here is an analytic solution for a source with a constant concentration. The boundary conditions for a semi-infinite system with continuous first-type source boundary conditions (Figure 23.6) are expressed as 0CC = , for x = 0 (23.33) 0= x C , for x = (23.34) Chapter 23 Modeling Contaminant Transport 12 with the initial condition (at t = 0): 0=C , for 0 < x < (23.35) It is assumed that flow is in the x-direction at a constant velocity vx and a longitudinal dispersion coefficient (Dx), defined as Lvx + Dm ( is the tortuosity of a medium and Dm is the molecular diffusion coefficient). The solution of Eq. (23.27) with boundary conditions (23.33) through (23.35) was given by van Genuchten and Alves (1982) + + - = +- fx f Uv D x fx f Uv D x RtD RUtx erfce RtD RUtx erfce C txC x x x X /2 / /2 / 2 ),( )( 2 )( 20 (23.36) where U = xx Dv 42 + , and erfc() is the complementary error function. The complementary error function erfc() is related to the error function erf() by )(1)( erferfc +=- (23.37) )(1)( erferfc -= (23.38) )()( erferf -=- (23.39) Values of erf() and erfc() are presented in Table 23.1. For = 0 (no decay), Eq. (23.36) reduces to the well-known Ogata and Banks (1961) solution + + - = fx fxD xv fx fx RtD Rtvx erfce RtD Rtvx erfc C txC x x /2 / /2 / 2 ),( 0 (23.40) Chapter 23 Modeling Contaminant Transport 13 Table 23.1. Values of erf() and erfc() x erf(x) erfc(x) -3.0 -1.000000 2.000000 -2.8 -0.999925 1.999925 -2.6 -0.999764 1.999764 -2.4 -0.999311 1.999311 -2.2 -0.998137 1.998137 -2.0 -0.995323 1.995323 -1.8 -0.989091 1.989091 -1.6 -0.976348 1.976348 -1.4 -0.952285 1.952285 -1.2 -0.910314 1.910314 -1.0 -0.842701 1.842701 -0.8 -0.742101 1.742101 -0.6 -0.603856 1.603856 -0.4 -0.428392 1.428392 -0.2 -0.222703 1.222703 0.0 0.000000 1.000000 0.2 0.222703 0.777297 0.4 0.428392 0.571608 0.6 0.603856 0.396144 0.8 0.742101 0.257899 1.0 0.842701 0.157299 1.2 0.910314 0.089686 1.4 0.952285 0.047715 1.6 0.976348 0.023652 1.8 0.989091 0.010909 2.0 0.995323 0.004677 2.2 0.998137 0.001863 2.4 0.999311 0.000689 2.6 0.999764 0.000236 2.8 0.999925 0.000075 3.0 1.000000 0.000000 The second term in Eq. (23.40) is often negligible because erfc() is near zero when xvx/Dx is large giving - = - fx f Uv D x RtD RUtx erfce C txC x X /2 / 2 ),( )( 20 (23.41) The following example illustrates the use of analytical soltution for the calculation of a breakthrough curve and a concentration profile along a column. Chapter 23 Modeling Contaminant Transport 14 Example 23.2: The concentration of a contaminant is maintained at 10 g/cm3 with time at x = 0 cm in a column with a cross-sectional area of 1 cm2 . The longitudinal dispersivity is 1 cm. A water flux of 0.35 cm3 /min is supplied continuously at x = 0. The porosity of the medium is 0.35, and the molecular diffusion coefficient of the contaminant is 1.0 x 10-5 cm2 /min. Calculate the contaminant concentration in the direction of flow at 100 min and the breakthrough curve at x = 25 cm for decay constants = 0.0 and 0.01 min-1 . Solution: The calculation of the concentration is similar to that for an instantaneous point source. The concentration is calculated using the Microsoft Excel file "analytic1.xls" located in directory chap23 in the accompanying CD. There are two calculation sheets in the file "analytic1.xls". The calculation sheet "1d-continuous-source" is designed for this example. Concentration along the ground-water flow direction is shown in Figure 23.7. As indicated in the figure, contaminant decay reduces the total amount of mass in the system. The breakthrough curves are shown in Figure 23.8. Concentration along the flow direction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 50 100 150 200 250 x (cm) Concentration Figure 23.7. Concentration profiles in the direction of flow for a continuous source. Chapter 23 Modeling Contaminant Transport 15 Breakthrough Curve 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 5 10 15 20 25 30 Time (min) Concentration Figure 23.8. Breakthrough curves at x = 25 cm for a continuous source. Semi-infinite System with a Third-type Source Boundary Condition Figure 23.9. A problem of one-dimensional mass transport with a continuous source that provides a constant mass flux. This problem involves mass transport in a one-dimensional semi-infinite system with a third-type source boundary condition (Figure 23.9). Mathematically, these boundary conditions are written as x C DCvCv xxx -=0 , for x = 0 (23.42) Chapter 23 Modeling Contaminant Transport 16 0= x C , for x = (23.43) The initial condition (at t = 0) is: 0=C , for 0 < x < (23.44) The solution of Eq. (23.27) for conditions (23.42) through (23.44) was given by van Genuchten and Alves (1982) as - + + + = -- fx f Uv D x xfx fxR t v D x x x RtD RUtx erfce Uv U RtD Rtvx erfce D v CtxC x xf x x /2 / /2 / 2 ),( )( 2 )(2 0 + - + + fx f Uv D x x RtD RUtx erfce Uv U x x /2 /)( 2 (23.45) For = 0, the solution (23.41) reduces to + - = - - fx fx t/RD t/Rvx fx x fx fx e RD tv t/RD t/Rvx erfcCtxC 4 )( 2 0 2 22 1 ),( + ++- fx fxD xv fx x x x t/RD t/Rvx erfce RD tv D xv x x 2 1 2 1 2 (23.46) as given by Lindstrom et al.(1967). The calculated concentrations by Eqs. (23.40) and (23.46) will be different at early times. 23.5 Analytical Solutions of Mass Transport Equations in Three Dimensions For contaminant transport in three-dimensional space with a constant linear flow velocity in the x-direction, the mass transport equation is CR x C v z C D y C D x C D t C R fxzyxf - - + + = 2 2 2 2 2 2 (23.47) Chapter 23 Modeling Contaminant Transport 17 The sections that follow present various solutions to this equation. Instantaneous Point Source Figure 23.10. Geometry for a point source in a three-dimensional domain. The next solution applies to mass transport in three dimensions with an instantaneous point source. The boundary conditions, as depicted in Figure 23.10 are 0MM = , At x = 0, y = ys and z = zs (23.48) 0,0 = = x C C , for x = (23.49) with the initial condition: 0=C , for 0 < x < and t = 0 (23.50) The solution is 111 0 ),,,( ZYX nR M tzyxC f = (23.51) with fy s RtD yy fy e RtD tyY /4 )( 1 2 /4 1 ),( - - = (23.52) Chapter 23 Modeling Contaminant Transport 18 += + - - - fz s fz s RtD zz RtD zz fz ee RtD tzZ /4 )( /4 )( 1 22 /4 1 ),( (23.53) where M0 is the initial contaminant mass at x = 0, y = ys and z = zs, x is the distance from the source to a calculation point down gradient, y is the lateral distance from the source to the calculation point, z is the distance from the water table to the calculation point, and X1(x, t) is given by Eq. (23.32). For a case with complete vertical mixing of dissolved mass in an aquifer or an average concentration in the vertical direction, Eq. (23.51) simplifies to this two-dimensional form: 11 0 ),,,( YX bnR M tzyxC f = (23.54) where b is the aquifer thickness. For a point source, the maximum concentration at the center of the plume where x = vt, y = ys, and z = zs is expressed as t zyxff ss e DDDRtnR M tzzyyvtxC ==== 2/12/3 0 )()/(4 ),,,( (23.55) As time increases, the maximum concentration is gradually reduced due to dispersion (Figure 23.11). The spreading due to dispersion is 3 )2( 3 )2( , 3 )2( 2/12/12/1 tD and tDtD z z y y x x === (23.56) where is the standard deviation so that x, y, and z represent three spreading lengths within which about 99.7% of the mass is contained. Chapter 23 Modeling Contaminant Transport 19 Figure 23.11. Plan view of a plume developed from a point source at several different times (from Domenico and Schwartz, 1998. Physical and Chemical Hydrogeology). Copyright 1998 by John Wiley & Sons, Inc. Reprinted by permission of John Wiley & Sons, Inc. Example 23.3: The contaminant source mass is 10 kg at x = 0, y = 0, and z = 0. The values of longitudinal, transverse, and vertical dispersivities are 1, 0.1, and 0.01 m, respectively. The linear ground water velocity is 0.35 m/day. The porosity is 0.35, and the molecular diffusion coefficient is 1.0 x 10-6 m2 /day. Calculate the contaminant concentration along and transverse to the direction of flow at 100 days, assuming the aquifer to be thick (that is, complete 3-D spreading). If the thickness is assumed to be 5 or 50 m, calculate the contaminant concentration in the aquifer. Solution: The concentration is calculated using the Microsoft Excel file "analytic2.xls" located in directory Chap23 on the accompanying CD. There are two calculation sheets in the file "analytic2.xls". The Excel sheet "3d point source" is based on the formulation (23.51) while the Excel sheet "2d point source" is based on the Eq. (23.54). The calculated concentrations in the direction of flow are shown in Figure 23.12. Concentration profiles transverse to the direction of flow are shown in Figure 23.13. Apparently, the average concentration in an aquifer of finite Chapter 23 Modeling Contaminant Transport 20 thickness may be larger or smaller than that in a thick aquifer depending on the thickness of the aquifer. Vertical dispersion is an important factor for diluting the concentration. 0.0E+00 2.0E-03 4.0E-03 6.0E-03 8.0E-03 1.0E-02 1.2E-02 1.4E-02 1.6E-02 0 50 100 150 200 x (m) Concentration(kg/m3) b=5 m b =50 m 3d Figure 23.12. Concentration profiles in the direction of flow for two- and three-dimensional point sources. 0.0E+00 2.0E-03 4.0E-03 6.0E-03 8.0E-03 1.0E-02 1.2E-02 1.4E-02 1.6E-02 -60 -40 -20 0 20 40 60 Y (m) Concentration(kg/m3) b=5 m b =50 km 3d Figure 23.13. Concentration profiles transverse to the direction of flow for two- and threedimensional point sources. Chapter 23 Modeling Contaminant Transport 21 Domenico-Robbins Solution Figure 23.14. Geometry for a patch source at the water table with unidirectional flow (modified (from Domenico and Robins, 1985). Domenico (1987) derived an approximate analytical solution for mass transport in a constant velocity field (x-direction) for a continuous source of regular geometry (Figure 23.14). The contaminant source starts at the water table and extends downward. It has a width Y and a height Z with a concentration of C0. The solution is 222),,,( ZYXtzyxC = (23.57) with ( ) ( ) +- = +- fx xxfx vD D xv RtD vDRtvx erfce C txX xx x x /2 /41/ 2 ),( 2/12/411 20 2 2/12 (23.58) - - + = 2/12/12 )/(2 2/ )/(2 2/ 2 1 ),( xyxy vxD Yy erf vxD Yy erftyY (23.59) Chapter 23 Modeling Contaminant Transport 22 - - + = 2/12/12 )/(2)/(22 1 ),( xzxz vxD Zz erf vxD Zz erftzZ (23.60) where vx is the linear ground-water velocity. For = 0, Eq. (23.57) as shown in Domenico and Robbins (1985) reduces to 223),,,( ZYXtzyxC = (23.61) with - = fx fx RtD Rtvx erfc C txX /2 / 2 ),( 0 3 (23.62) For two-dimensional problems (source height being same as the aquifer thickness), the term Z2(z,t) is set to 1 in Eqs. (23.57) and (23.61). Similarly, for a one-dimensional problem, the terms Y2(y,t) and Z2(z,t) are set to 1 in Eqs. (23.57) and (23.61). An example calculation will be given for the Domenico-Robbins solution in section 23.9 using the code CONTz.exe that is introduced in that section. AT123D solution (Yeh, 1981) The solution to Eq. 23.22 can be derived by the superposition of the Green's functions. The Green function of a mass transport problem is the solution of mass transport for an instantaneous point source, with unit source mass flux under the same boundary and initial conditions. The solution to the mass transport problem is = 1 0 ),,,;,,,(),,,( t o R ddddtzyxMGtzyxC (23.63) where M is the strength of the source, G is the Green function of the problem, t is the time, t1 is the source duration, R0 is the source region, and x, y, and z are the coordinates of the system. Chapter 23 Modeling Contaminant Transport 23 Figure 23.15. Source geometry and aquifer specifications used in the code AT123D. Yeh (1981) used the Green function scheme to derive a three-dimensional solution for sources with a regular geometry (Figure 23.15). The contaminant source is assumed to be a rectangular block with a beginning coordinate of (L1, B1, H1) and an ending coordinate of (L2, B2, H2). The thickness and width of the aquifer are B and H, respectively. The geometric integration was calculated analytically in the solution. Thus, numerical integration is only required over the time coordinate for a non-instantaneous source. The solution is > < = sourceousinstantaenan),0,,()0,,()0,,( Tt,),,(),,(),,( Tt,),,(),,(),,( t)z,y,C(x, 0 1 1 tzFtyFtxF nR M dtzFtyFtxF nR M dtzFtyFtxF nR M zyx f T o zyx f t o zyx f (23.64) Chapter 23 Modeling Contaminant Transport 24 where M1 is the mass release rate (MT-1 ), M0 is the contaminant mass, n is the porosity, Rf is the retardation factor, T is the source duration, and Fx, Fy, and Fz are the integration of the Green functions within the source in the x, y, and z directions, respectively. The solution is programmed in the code AT123D (Yeh, 1981). This code can be applied to a variety of mass transport problems involving a homogeneous medium, various source geometries, and unidirectional flow. Comparatively, the AT123D solution is much more flexible than other solutions that we have introduced. AT123D has been incorporated into the computer program (CONTz.exe) and will be introduced in section 23.9. 23.6 Other Useful Analytical Solutions Contaminant Transport with Retardation Most of the solutions thus far have included the effects of sorption. What we have not shown, however, is how easily any analytical solution to an advection-dispersion equation can be modified to include sorption. With sorption described in terms of a distribution coefficient, Kd, parameter values are adjusted using the retardation factor, Rf. Because Rf is a constant (Rf = 1 + Kd b/n), Eq. (23.13) can be rewritten as t C nR r C R v C R D fff =+- )()( or t C rCvCD =+- ')'()'( (23.65) where fff nR r r R v R === ', v ', D 'D (23.66) Chapter 23 Modeling Contaminant Transport 25 The resulting Eq. (23.65) is a simplified form of the transport equation where the effects of sorption are incorporated neatly in the parameter values. Figure 23.16 shows the transport of a retarded species as compared with an unretarded species. Figure 23.16.The transport of a retarded species as compared to an unretarded species (from Domenico and Schwartz, 1998. Physical and Chemical Hydrogeology). Copyright 1998 by John Wiley & Sons, Inc. Reprinted by permission of John Wiley & Sons, Inc. Contaminant Transport under Unsaturated Flow Conditions Under unsaturated conditions, the porosity (n) in Eq. 23.13 is replaced by soil moisture ( = nSw). In addition, the ground-water flow velocity and retardation factor need to be calculated differently. The equation for the linear ground-water velocity vector V is rewritten as ( ) - - -== z h nS K y h nS K x h nS K vvvv w zz w yy w xx zyx ,,,, (23.67) The equation for retardation is rewritten as w db f nS K R += 1 (23.68) where b is the bulk density of the soil or rock, Kd is the distribution coefficient of a chemical species in the soil or rock, n is the porosity, and Sw is the water saturation. Example 23.3: The total organic carbon content in the soil is 0.01% and the partition coefficient of the contaminant between organic carbon and water (Koc) is 200 cm3 /g. The bulk density of the Chapter 23 Modeling Contaminant Transport 26 aquifer material is 1.65 g/cm3 ; the average soil moisture content is 0.8; the soil porosity is 0.25; the specific discharge is 0.01 m/day; and the longitudinal dispersivity is 1 m. Calculate the ground-water velocity, the advective velocity of the retarded contaminant, and the retarded dispersion coefficient. Solution: The distribution coefficient is gcmgcmfKK ococd /2%)01.0)(/200( 33 === The retardation factor is 5.17 )8.0)(25.0( ))/2)(/65.1( 1 33 =+= gcmcmg Rf The linear ground-water velocity is daym daym nS q v w /05.0 )8.0)(25.0( /01.0 === The advective velocity of the retarded contaminant is daym daym R v v f /0029.0 5.17 /05.0 ' === The longitudinal dispersion coefficient is daymmdaymvD Lx /05.0)1)(/05.0( 2 === The retarded longitudinal dispersion coefficient (D) is daym daym R D D x x /0029.0 5.17 /05.0 2 2 ' === Chapter 23 Modeling Contaminant Transport 27 Contaminant Transport with Multiple Pulse Source Conditions Figure 23.17. Representation of multiple pulse sources in time. If the source concentration varies with multiple pulses (Figure 23.17), the analytical solution is ),,,()(),,,( 1 1 0,10, - = - --= i i j jj TtzyxGCCtzyxC (23.69) where Ci,0 is the source concentration for the time interval i (i=1,2,...,I), G is the analytical solution of the problem with a unit source concentration, Ti is the duration of the time interval i, and t is time. It is assumed that T0 = 0 and C00 = 0. For a single-pulse source, the solution is written as [ ] >-- < = 110 10 ,),,,(),,,( ),,,,( ),( TtTtzyxGtzyxGC TttzyxGC txC (23.70) Chapter 23 Modeling Contaminant Transport 28 Eq. (23.70) indicates that the analytical solution for a single-pulse source at times after the source becomes inactive is the superposition of two continuous sources with opposite strengths. A simple procedure is to calculate the concentration for t > T1 using a spreadsheet program: (1) enter time in the first column; (2) calculate the concentration values using C0G(x,y,z,t) in the second column; (3) copy the column into the third column and shift the starting time to T1 to get C0G(x, y, z, t-T1); and (4) subtract the concentration in the third column from the concentration values in the second column to get C0[G(x,y,z,t)-G(x,y,z,t-T1)]. Example 23.4: Calculate the breakthrough curve at x = 25 cm for example 23.2 if the source is active for only 24 min. Solution: This problem gives I = 2, T1 = 25 min, C10 = 1 g/cm3 , T2 = 54 min, and C20 = 1 g/cm3 . The analytical solution is written as > -- - - < - = - - 1 1 )( 2 1 )( 20 , /2 /)( /2 / , /2 / 2 ),( Tt RtD RTtUx erfc RtD RUtx erfce Tt RtD RUtx erfce C txC fx f fx f Uv D x fx f Uv D x x X x X (23.71) A MS Excel file "analytic3.xls", located in the directory Chap23 on the accompanying CD, is used to calculate the breakthrough curve. The results of this calculation are shown in Table 23.2 and displayed in Figure 23.18. The curve labeled c(g/cm3 ) is the breakthrough curve for the pulse source, whereas the curve labeled c2(g/cm3 ) is the breakthrough curve for a continuous source. Chapter 23 Modeling Contaminant Transport 29 Table 23.2. Breakthrough curves for a pulse source as calculated by Eq. (23.71) t(min) c2(g/cm3 ) c1(g/cm3 ) c(g/cm3 ) 8 0.00E+00 0.00E+00 0.00E+00 10 3.98E-04 0.00E+00 3.98E-04 12 3.98E-03 0.00E+00 3.98E-03 14 1.88E-02 0.00E+00 1.88E-02 16 5.58E-02 0.00E+00 5.58E-02 18 1.22E-01 0.00E+00 1.22E-01 20 2.15E-01 0.00E+00 2.15E-01 22 3.26E-01 0.00E+00 3.26E-01 24 4.43E-01 0.00E+00 4.43E-01 26 5.55E-01 0.00E+00 5.55E-01 28 6.56E-01 0.00E+00 6.56E-01 30 7.41E-01 0.00E+00 7.41E-01 32 8.09E-01 3.98E-04 8.09E-01 34 8.62E-01 3.98E-03 8.58E-01 36 9.03E-01 1.88E-02 8.84E-01 38 9.32E-01 5.58E-02 8.76E-01 40 9.53E-01 1.22E-01 8.32E-01 42 9.68E-01 2.15E-01 7.54E-01 44 9.79E-01 3.26E-01 6.53E-01 46 9.86E-01 4.43E-01 5.43E-01 48 9.91E-01 5.55E-01 4.35E-01 50 9.94E-01 6.56E-01 3.38E-01 52 9.96E-01 7.41E-01 2.55E-01 54 9.97E-01 8.09E-01 1.88E-01 56 9.98E-01 8.62E-01 1.36E-01 58 0.998908 0.902574 0.096334 60 0.999301 0.932045 0.067256 60 0.999301 0.932045 0.067256 Chapter 23 Modeling Contaminant Transport 30 Breakthrough curve at x = 25 cm 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 10 20 30 40 50 60 t (min) Concentration c2(g/cm3) c(g/cm3) Figure 23.18. Examples of breakthrough curves for a pulse source. Contaminant Transport in a Fractured Medium In fractured media, several processes influence concentration distributions: sorption on to the fracture surface, diffusion into the matrix, and advective and dispersive transport along the fracture. Most theories for contaminant transport in fractured media assume that there is no advective mass transport in the matrix. Chapter 23 Modeling Contaminant Transport 31 Figure 23.19. Mass transport in a set of parallel fractures with diffusion into the matrix. A simple parallel fracture model by Sudicky and Frind (1981) illustrates the basic mass transport mechanisms in fractured media (Figure 23.19). The differential equation for solute transport along a fracture is written as b q CR x C v x C D t C R fxxf -- - = 2 2 (23.72) where C is the contaminant concentration, Rf is the retardation factor representing sorption onto the fracture wall, b is half of the fracture aperture, x is the distance along the fracture in the direction of flow, and q is the diffusive flux in or out of the matrix . The retardation in the fracture is expressed as b K R f f += 1 (23.73) where Kf is the distribution coefficient for the fracture surface, defined as the ratio of the solute mass adsorbed per unit area of fracture-wall surface to the solute concentration in the solution. The differential equation for mass transport in the matrix is Chapter 23 Modeling Contaminant Transport 32 '' '' ' 2 2 * CR y C D t C R - = (23.74) where R is the retardation factor in the matrix, C is the concentration in the matrix, D* is the diffusion coefficient in the matrix, and y is the direction perpendicular to the flow in the fracture. The retardation factor and the diffusion coefficient in the matrix are ' 1' n K R bm += (23.75) mDD =* (23.76) where Km is the distribution coefficient for sorption in the matrix, b is the bulk density of the matrix material, n is the matrix porosity, is the tortuosity of the matrix, and Dm is the molecular diffusion coefficient. The matrix diffusion term q can be expressed in terms of Fick's first law as y C Dnq -= ' ' * (23.77) The effect of matrix diffusion is similar to retardation where the transport velocity of the contaminant along the fractures is less than the linear velocity of the water. Simulation results (CRWMS M&O, 1998) concluded that (1) as the diffusion coefficient and matrix porosity of the matrix blocks increase, contaminant spreading in the fracture is reduced; (2) as the fracture spacing is increased under a same specific discharge in the fractured medium, matrix diffusion is less effective in attenuating the spread of the contaminant. Here, we present an analytical solution by Tang et al. (1981) for mass transport in a fracture and a matrix with infinite matrix spacing, neglecting the effect of dispersion in the fracture. It turns out that this solution is sufficiently straightforward for a spreadsheet calculation. The equation is written as Chapter 23 Modeling Contaminant Transport 33 > ++ - < = - - 0', '2'22 0',0 ),,( 2/12/10 2/12/1 TW T W erfceW T W erfcee C T tyxC WWv Rz (23.78) where C(x,y,t) is the concentration in the matrix. T and W are expressed as 2/1 ' -= v zR tT f (23.79) and )( ')'(' * 2/1* bx D R b DRn v z W -+= (23.80) Concentration in fractures may be determined from Eq. (23.78) by setting x = b. Example 23.4: It is known that C0 =1, Rf = 1, v = 600 m/yr, 2b = 0.001 m, n'=0.25, R' = 1, D* = 0.00315576 m2 /yr, = 0 yr-1 , x = 0.00005 m. Calculate: (1) the breakthrough curve for z = 100 m, t = 7.5 yr, nt = 30; (2) the concentration along a fracture for t = 1000 yr, z = 100 m, and nz = 30. Solution: A Microsoft Excel file "analytic4.xls", located in the directory Chap23 on the accompanying CD, is used to calculate the breakthrough curve and the concentration along the fracture. The symbols in the worksheet are the same as those in Eqs. (23.78) through (23.80). The variables in yellow, blue, and green background can be changed for new calculations. When users change the values of these variables, click the calculate button to update the values in result cells and charts. The results of the calculation are shown in Figures 23.20a and 23.20b, respectively. Without matrix diffusion, there is only advective transport in the fracture, so that the value of concentration should be 1 on the two curves. Chapter 23 Modeling Contaminant Transport 34 Breakthrough curve 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 Time (yr) Concentration (a) Concentration alonga fracture 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 500 1000 1500 2000 2500 3000 3500 Distance along a fracture (m) Concentration (b) Figure 23.20. Results of a calculation with the simplified solution of Tang et al. (1981). Panel (a) shows a breakthrough curve at 100 m; Panel (b) shows the concentration profile along the fractures. 23.7 Particle Tracking Methods Numerical approaches have the potential to address problems of mass transport in complex natural systems. Unfortunately, there are significant disadvantages in using these codes. One problem is the steep learning curve in actually learning to set up and run cases, and to provide accurate results. Another is the significant burden of field data that are required to Chapter 23 Modeling Contaminant Transport 35 describe a problem fully. Particle-tracking approaches, which have been used to advantage in many practical problems, provide a useful compromise between power and complexity. For many problems one really doesn't need a full-blown mass transport simulation. The question may be simply to identify the position of a contaminant front and how it moves in space and time, to determine why production wells are being contaminated, or to trace water originating from different sources. Particle tracking describes the path that hypothetical particles follow through time when added at a source. This application assumes that advection is the main mass transport process that influences the migration of the particles. The track of the particle with time defines a pathline through the system. In other words, the particle travels in the direction of ground-water flow with the linear ground-water velocity. A particle-tracking code is normally run together with a ground-water flow code (for example, MODFLOW). The solution to the flow code defines a hydraulic-head distribution that with a simple application of Darcy's equation can be transformed into a two-dimensional or three-dimensional velocity field. This velocity field provides information to the particles as to where and how fast they should move. For a typical problem, one would start by adding reference particles at some source. Given the details of the velocity field and some small elapsed time, the particles are stepped down the pathlines. Generally, this simple scheme neglects key chemical processes. However, making simple adjustments to the velocity values can include the effects of retardation. Recall that the velocity of the contaminant equals the velocity of the water divided by the retardation factor. Here is the algorithm that lets one define the evolution of a contaminant plume in time and space. New coordinates of particles can be calculated from old coordinates. ),,(1 nnnx f nn zyxv R t xx +=+ (23.81) Chapter 23 Modeling Contaminant Transport 36 ),,(1 nnny f nn zyxv R t yy +=+ (23.82) ),,(1 nnnz f nn zyxv R t zz +=+ (23.83) where (xn, yn, zn) is the old coordinates of a particle, t is the time step, (vx, vy, vz) is the velocity along the flow path, and (xn+1, yn+1, zn+1) is the new coordinates of the particle. Improvement was made to calculate the velocity more accurately (Pollock, 1994; Zheng, 1990). The best known of the particle-tracking codes is MODPATH (Pollock, 1994). This code functions as a post-processor to MODFLOW and with minimal additional data provides a useful mass transport capability. A particle-tracking code like MODPATH could be used to determine what the source of water is for a well emplaced next to a stream (Figure 23.21) (Franke et al., 1998). Before the well was installed, all of the ground water discharged into the stream. This ground water originated from recharge from the land surface and lateral inflow from the valley walls. Once the well was installed, these sources were portioned between the well and the stream according to the pumping rate. Simply adding the pathlines to the cross section (Figure 23.21) makes it possible to differentiate between flow to the well and flow to the stream. Chapter 23 Modeling Contaminant Transport 37 Figure 23.21. Hypothetical model-calculated ground-water flow paths at equilibrium from the water table to a nearby stream and to a discharging well screened at the bottom of the aquifer (from Franke et al, 1998a). An example from Masterson et al. (1997) illustrates how MODFLOW could be used to predict possible patterns of contaminant spreading. Particles were added at a source that coincided with an actual landfill and tracked through two slightly different versions of the same hydrogeologic system (Figure 23.22a,b). In one case, the horizontal hydraulic conductivity of the Chapter 23 Modeling Contaminant Transport 38 moraine sediments was assigned a hydraulic conductivity of 50 ft per day. In the second model, they were assigned a hydraulic conductivity of 150 ft per day. These changes to the hydrogeologic setting had relatively minor effects on the simulation of the water-table configuration (compare Figures 23.22a and b). However, the simulated patterns of contaminant Figure 23.22. Simulated effects of increasing horizontal hydraulic conductivity of moraine sediments on the water-table configuration and ground-water flowpaths near the Landfill-1 contaminant plume, western Cape Cod (from Masterson et al, 1997). Chapter 23 Modeling Contaminant Transport 39 migration differed appreciably (compare Figures 23.22a and b). The second of the simulation trials (Figure 23.22b) best matched the known configuration of the landfill plume. 23.8 Numerical Solutions of Mass-Transport Equations The numerical approaches are a family of computer-based techniques for solving contaminant transport equations. They approximate forms of the advection-dispersion equation as a system of algebraic equations, or alternatively simulate transport through the spread of a large number of moving reference particles. Whatever the procedure used, it invariably has to be coded for solution on a high-speed computer. Numerical approaches easily accommodate variability in flow and transport parameters (e.g., hydraulic conductivity, porosity, dispersivity, cation exchange capacity, and etc.). This flexibility in representing parameters facilitates modeling of layering or other more complex patterns of variability in two and three dimensions. One can simulate the complex plumes that occur with real problems. Thus, numerical approaches are readily adapted to site-specific problems and thus are particularly useful in practice. Other important features of the numerical approaches are the flexibility in implementing complex boundary conditions and in accounting for a variety of important mass transport processes. Codes are now available to model (1) radionuclide chains in which radioisotopes and their daughters decay, each with a variable half-life; (2) equilibrium reactions such as precipitation, dissolution, and complexation; and (3) competitive sorption or ion exchange. Developing the increased sophistication to deal with reactions is necessary to model contaminant groups such as trace metals or radionuclides accurately. Chapter 23 Modeling Contaminant Transport 40 The modeling process can be examined independently of a particular solution technique. The starting place in nearly all cases is with one of the various forms of the advection-dispersion equation. Because of the power of the numerical approaches, the transport equation is commonly formulated in two or three dimensions. Although our emphasis here is on mass transport involving a single contaminant, modeling approaches have evolved to the point where a coupled set of differential equations can be formulated for problems involving several reacting species. One of the first steps in developing a computer model is to subdivide the region in terms of cells or elements. This process makes it possible to account for varying parameter values within the domain. Furthermore, subdividing the region also serves to define the nodes at which concentrations are calculated. The way in which nodes are defined and how the domain is subdivided (e.g., squares, rectangles, triangles, stream tube segments, and etc.) depends on the specific numerical technique. The Otis Air Base case study presented later shows how one goes about subdividing a region and assigning boundary conditions. There is really no difference in the boundary conditions and transport parameters required in the numerical approaches as compared to the analytical approaches. Complicating the situation with numerical models is the need to provide boundary conditions at a large number of node points or nodal blocks (areas around the nodes). For example, values of concentration or loading rates defining various boundary conditions now need to be specified for all nodes located along the boundary of the domain. Initial conditions and transport parameters are specified for all nodes except in some cases for nodes with constant concentrations. It is no real problem to assign transport parameters to nodal areas, other than perhaps not knowing what the values are. Velocity values are an important exception because they are more difficult to specify. In most real systems with pumping and injection wells, and a variable Chapter 23 Modeling Contaminant Transport 41 hydraulic conductivity field, velocities can be extremely variable. Both the magnitude and the direction of flow can change in space and time. In terms of modeling, measured or guessed velocity values are not adequate. Continuity considerations in all numerical solutions require a smooth and accurate representation of the velocity field. This kind of field can only be obtained by simulation with a flow model. Velocity values are calculated by applying the Darcy equation with calculated values of hydraulic and known parameters, or directly calculating velocity values at the nodes. Figure 23.23 illustrates how flow and mass transport models are used together to predict contaminant distributions. Input necessary for the flow model could include the hydraulic conductivity distribution, the water-table configuration, and other boundary conditions. The transport model is coupled to the flow model by the velocity terms as shown in Figure 23.23. Figure 23.23. Illustration of numerical solution of flow and mass transport equations. It is beyond the scope of our simplified treatment here to discuss the procedures for solving mass transport equations numerically. With experience, one can learn to use various public domain numerical mass transport codes. The most commonly used codes include Chapter 23 Modeling Contaminant Transport 42 MT3DMS (Zheng and Wang, 1998), SUTRA (Voss, 1984), HST3D (Kipp, 1997), and SWMS_3D (Šimůnek et al., 1995). MT3DMS is the modular three-dimensional transport model (MT3D) with significantly expanded capabilities (http://hydro.geo.ua.edu/mt3d). SUTRA is a two-dimensional solute or energy transport finite-element code for saturated or unsaturated flow (http://h2o.usgs.gov/software). The HST3D is a heat and transport three-dimensional code with integrated finite difference scheme (http://h2o.usgs.gov/software). The SWMS_3D is a threedimensional water unsaturated flow and solute transport code developed by U. S. Salinity Laboratory (http://www.ussl.ars.usda.gov/MODELS/MODELS.HTM). Case Study in the Application of a Numerical Model A model study of contamination at Otis Air Base at Cape Cod, Massachusetts (Figure 23.24) (LeBlanc, 1984a), illustrates some of the steps in using a mass transport model. At Otis Air Base, approximately 1740 m3 /day of treated wastewater has been disposed of in the subsurface since 1936. The disposal unit is an unconfined sand and gravel aquifer approximately 35 m thick. A zone of contamination has developed that is approximately 915 m long, 23 m thick, and 3.35 km long. This site is being studied by the U.S. Geological Survey as part of the Toxic-Waste Ground-Water Contamination Program (Franks, 1987; LeBlanc, 1984b; Ragone, ed., 1988). The treated sewage effluent contains above-background concentrations of Na+ , Cl-, ammonium, nitrate, phosphate, detergents, and several different volatile organic compounds. Boron, one of the trace metals in the effluent, was selected by LeBlanc (1984a) for a detailed model study because (1) its concentration has remained relatively constant in the effluent at 500 mg/L, (2) it has a relatively low background level in the native ground water (30 mg/L), and (3) it tends not to react chemically during transport. The objective in modeling was to guide the collection of data from the site and to test hypotheses concerning the character of contaminant migration. Chapter 23 Modeling Contaminant Transport 43 Figure 23.24. Sewage effluent plume at Otis Air Base on Cape Cod, Massachusetts (from LeBlanc, 1984a). The first step in developing the flow and boron transport models was to define the region of interest and to establish boundary conditions for flow. LeBlanc (1984a) used the following guidelines, which should apply generally to most studies. The domain had to be large enough to include all of the existing plume (Figure 23.24) as well as providing room for future spreading. The side and bottom boundaries downgradient of the disposal area are ponds, rivers, and saltwater bays. These features are a good choice for boundaries because in the absence of detailed hydraulic head data, they provide the best places to estimate boundary conditions. In effect by assuming ground-water discharges at these locations, boundaries for ground-water flow can be estimated in terms of constant head nodes or leakage fluxes. North of Coonamessett and Johns Ponds, there are no natural hydrogeologic boundaries. In this area, flow lines formed the side boundaries (Figure 23.24). The flow lines are imaginary boundaries, located so as not to intersect the contaminant plume. The 60-ft equipotential line arbitrarily defined the top or northern boundary. This boundary could have been placed anywhere north of the site with the proviso again that contaminants from the site could not intersect this boundary. The Konikow and Bredehoeft (1978) model requires that the site be subdivided into a regular finite difference grid. In this example, the grid consisted of 40 rows and 36 columns. The nodal blocks are square, with dimensions chosen so that there is not an unmanageably large Chapter 23 Modeling Contaminant Transport 44 number of cells (an upper limit might be 50 or so cells in the row/column directions). However, the blocks are sufficiently small to ensure that the plume is not localized in just a few cells. The plume in this case has a width of about seven cells. Figure 23.25 shows how the model replicates the region shape and how the boundary conditions for flow are included. Figure 23.25. Simulation grid with hydraulic properties and boundary conditions (from LeBlanc, 1984a). LeBlanc modeled the pattern of ground-water flow first. Because the density and viscosity of the contaminated ground water are nearly the same as the uncontaminated water, the distribution of hydraulic head and hence the velocity field are unaffected by the migration of the plume. The flow equation therefore can be solved first independently of the mass transport equation. Observations at the site indicate further that the gradients in hydraulic head do not change significantly with time. Thus ground-water flow is steady, which means that the groundwater velocities need to be calculated only once. Figure 23.26 compares the observed configuration of the water table with the "best fit" from a series of simulations. Such trials are designed to calibrate the model--a process of selecting the set of parameters that produces the best simulation of a known history. The Chapter 23 Modeling Contaminant Transport 45 calibration here took initial estimates of model parameters (for example, hydraulic conductivity, recharge rates) and adjusted them by trial and error until the model successfully reproduced the observed configuration of the water table. Because of the lack of data, no attempt was made to reproduce spatial variability in hydraulic conductivity or recharge rates. Nevertheless, the simulated and observed results compare well (Figure 23.26). The model as represented by the parameters and boundary conditions is not unique. Another set of different parameters and boundary conditions could give the same or a better fit. By keeping the hydraulic conductivity close to the measured values, the simulated hydraulic heads and resulting velocity field should be a reasonable representation of the flow system. Figure 23.26. Comparison of the observed (dashed lines) and simulated (solid lines) water-table elevations for November 1979 (from LeBlanc, 1984a). The only information required for the transport simulation with a conservative species like boron is the longitudinal and transverse dispersivities, plus the loading function at the source. In the absence of actual dispersivity data, values were selected from the literature for similar types of geologic materials. Here, this uncertainty in choosing a value is not a problem Chapter 23 Modeling Contaminant Transport 46 because analysis shows that with active recharge, simulated concentrations are not sensitive to the particular choice of dispersivity value. A lack of data may not be so easily dealt with in every model study. Contaminant loading was approximated by fixing the boron concentration at 500 mg/L at the four cells representing the infiltration beds. Tests with the model attempted to reproduce the historical spread of contaminants over a 40-year period between 1940 and 1978-79. Figures 23.27a and b compares the simulated and observed distributions of boron after 40 years. As LeBlanc (1984b) points out, the simulated path of the plume agrees reasonably well with that observed in 1978-79, particularly at concentrations above 50 mg/L. However, at lower concentrations, the center line of the simulated plume was located east of the observed center line. The simulated plume is also somewhat longer and wider than the observed plume. LeBlanc (1984b) explained the difference between observed and simulated plumes as follows. The observed plume could in fact be larger than is represented in Figure 23.27b because of uncertainty in establishing the edges of the plume at concentrations below 100 mg/L. Inaccuracies in the simulated flow field could have existed, producing somewhat more divergent flow than actually occurs. This problem could be related to the complex interaction between ground waters and surface waters. Finally, the model could contain undetected numerical dispersion. This study went on to evaluate the suitability of management options for the contamination. With reasonable confidence in the ability of this model to predict boron distributions, LeBlanc (1984b) examined the consequences of operating the site in the future as it is now or eliminating it. Chapter 23 Modeling Contaminant Transport 47 Figure 23.27. A comparison between (a) the simulated boron concentration 40 years after disposal began and (b) Observed boron concentration 1978 - 79 (from LeBlanc, 1984a). 23.9 Introduction to CONTz.exe A Window 32-bit based simulation, CONTz.exe, was implemented for modeling mass transport. The code is located in Chap23/program on the accompanying CD-ROM. When the program is installed, a program icon will be added to the Start/Programs Menu. CONTz.exe is Chapter 23 Modeling Contaminant Transport 48 activated by double-clicking the icon (Figure 23.28). The code includes six menus: File, Edit, Model, View, Options, and Help. The Edit menu is not implemented in the current version of the code. The File menu contains the commands Open, Save, Save As, Import Backdrop Map (Dib or BMP file), Print, Print Preview, Print Setup, and Exit. The imported backdrop map provides a base map for the output. The size of the backdrop map is specified using the left bottom and the top right coordinates of the rectangle containing the map. Figure 23.28. The CONTz menu showing the options for model calculations. The Model menu contains three options: Domenico-Robbins, AT123D, and Measured concentration. Selecting the Domenico-Robbins option prompts a dialog box, as shown in Figure 23.29. The dialog box is self-explanatory with reference to Eqs. (23.53) through (23.59). The units of length and time should be consistent in the input. For example, if the unit of length is m and the unit of time is day, the unit of ground-water velocity should be m/day. The unit of concentration can be any unit of mass per volume (ppm, ppb, and kg/m3 ). The code will calculate concentration values in an x-y plane for a depth of z from the water table. When the parameters Chapter 23 Modeling Contaminant Transport 49 in the dialog box are assigned, the simulation results are displayed by clicking the "OK" button. The simulation results are displayed as contour lines by default. Figure 23.29. Dialog box providing input parameters for the Domenico-Robins solution. The dialog box for the input parameters to AT123D is activated by clicking AT123D (Figure 23.30). Both continuous and instantaneous sources can be simulated with AT123D. In the current code, variable source strength is not implemented, although the original AT123D has that option. Specify the mass release rate by clicking the button "Mass release rate" (Figure 23.31). If the mass transport type is chemical, the output concentration will be mass/time*1000. For example, the output concentration will be g/m3 if the mass release rate is in kg/day or the total mass released is in kg and the time is in day. The soil and waste properties are activated by clicking the button "Properties" (Figure 23.32). The aquifer and source geometry are activated by clicking the button "Geometry" (Figure 23.33). Chapter 23 Modeling Contaminant Transport 50 Figure 23.30. Dialog box providing various control parameters for the AT123D solution. Figure 23.31. Dialog box for specifying the mass release rate for the AT123D solution. Chapter 23 Modeling Contaminant Transport 51 Figure 23.32. Dialog box for specifying soil and water properties in the AT123D solution. Figure 23.33. Dialog box for defining the aquifer and source sizes in the AT123D solution. In the View menu, the viewable area of the simulation domain can be controlled by clicking the command Zoom in, Zoom out, Zoom, Window, and Grid. Zoom in will shrink the viewable area by half in each direction. Zoom out will increase the viewable area by two times in each direction. Users can increase or decrease the viewable area by specifying a multiple factor of the dimension in the x- and y-direction. The Window command allows users to specify a viewable area by drawing a window on the screen. The Grid command returns you to the Chapter 23 Modeling Contaminant Transport 52 original size of the simulation domain. The Report command in the View menu displays the simulation report in Notepad or Wordpad. Three types of maps can be used to display the simulation results. The Contour command displays the simulation results with contour lines; the Colormap command adds a color to each calculation point; and the ColorContour command fills space between the contour lines with colors. In the Options menu, the Font command selects the Font for any text type in the display window, although the command will not change the size of the text. The Contour-lines command prompts a dialog box for specifying the properties of contour lines (Figure 23.34). The height of the contour label and the distance between labels are specified using the simulation scale. A maximum of 20 contour values can be entered in the table. The Save report command will save the simulation results in one file and concentration on the grid in SURFER format. Figure 23.34. Contour options for CONTz. Chapter 23 Modeling Contaminant Transport 53 Example 23.5: The simulation domain for a contamination problem has the dimensions 0 < x < 750 m and -50 m < y < 50m. The source is located within x = 0, -12.5m < y < 12.5m, and 0