# A Fourier-Laplace transform finite element method (FLTFEM) for the analysis of contaminant transport in porous media

код для вставкиСкачатьINTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) A FOURIER-LAPLACE TRANSFORM FINITE ELEMENT METHOD (FLTFEM) FOR THE ANALYSIS OF CONTAMINANT TRANSPORT IN POROUS MEDIA J. C. WANG* AND J. R. BOOKER Department of Civil Engineering, University of Sydney, NSW 2006, Australia SUMMARY In many practical situations, a simpli"ed two-dimensional model is often su$ciently accurate for practical purposes, however on some occasions it is still necessary to consider the full three-dimensional situation. Conventional numerical techniques such as the "nite di!erence and "nite element methods are able to deal with contaminant transport under 3D conditions, however the computational e!ort involved in using them is often very large. Under some circumstances it is possible to perform 3D analyses more e$ciently. A particular case is when the geometry is 2D but the problem is actually 3D. In such cases the analysis can be simpli"ed by application of a Fourier transform. In this paper, the dimensions of the contaminant transport problem will be reduced by one and the transient mass transport equation will be transformed into a time-independent mass transport equation through the use of a combination of integral transforms (Laplace and Fourier transforms) and a "nite element technique. The mathematical model developed here is called Fourier}Laplace transform "nite element method (FLTFEM). By using the FLTFEM a three-dimensional transient problem can be analysed as a two-dimensional time-independent problem and the computational cost will thus greatly be reduced. Copyright ( 1999 John Wiley & Sons, Ltd. KEY WORDS: contaminant transport; Fourier transform; Laplace transform INTRODUCTION In many important practical situations, a simpli"ed two-dimensional model is often su$ciently accurate for practical purposes; however on some occasions it is still necessary to consider the full three-dimensional situation. Conventional numerical techniques such as the "nite di!erence and "nite element methods are able to deal with contaminant transport under 3D conditions; however the computational e!ort involved in using them is often very large. Under some circumstances it is possible to perform 3D analyses more e$ciently. A particular case is when the geometry is 2D but the problem is actually 3D. In such cases the analysis can be simpli"ed by application of a Fourier transforms. Although the combined use of integral transforms and approximate numerical methods to solve linear di!erential equations can be seen in the literature,1}7 it appears most of them are only based on the application of a Laplace transform. Booker and Leo8 were probably the "rst to use * Correspondence to: J. C. Wang, Department of Civil Engineering, University of Sydney, Sydney, NSW 2006, Australia. E-mail: jc.wang@civil.usyd.edu.au CCC 0363}9061/99/141763}34$17.50 Copyright ( 1999 John Wiley & Sons, Ltd. Received October 1998 Revised 26 October 1998 1764 J. C. WANG AND J. R. BOOKER a combination of both Laplace transform and Fourier transform and boundary element methods to study the e!ects of contaminant migration from waste repositories. In this paper, the dimensions of the contaminant transport problem will be reduced by one and the transient mass transport equation will be transformed into a time-independent mass transport equation through the use of a combination of integral transforms (Laplace and Fourier transforms) and a "nite element technique. The mathematical model developed here is called Fourier}Laplace transform "nite element method (FLTFEM). By using the FLTFEM, a threedimensional transient problem can be analysed as a two-dimensional time-independent problem and the computational cost will thus greatly be reduced. The numerical method (FLTFEM) developed here may be applied to assess the impact of contamination transport (see Figure 1) due to a surface waste disposal repository or a prismatic waste repository which the geometry being considered is the same along one spatial dimension. One restriction of FLTFEM is that the material properties must not vary in one spatial dimension, although they may vary in the remaining two spatial dimensions. For example, in the case of contamination due to a crack on the bottom of a long channel (Figure 1(b)), the geometry and the material properties are the same along the y direction, however there may be layered soils with di!erent properties in the z direction (Figure 1(b)) or the x direction. Figure 1. Contamination problems due to a prismatic waste repository and a surface land"ll. (a) A surface land"ll; (b) a case of a prismatic waste repository Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1765 ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA In the FLTFEM, a Gauss}Lengendre quadrature formula is used to invert the Fourier transform while in the case of the Laplace transform, inversion is performed using a convenient technique developed by Stehfest.9, 10 2. DEVELOPMENT OF FLTFEM FOR CONTAMINANT TRANSPORT PROBLEMS 2.1. Governing equations The theory of a regular two- or three-dimensional fracture network within fractured media developed by Rowe and Booker11, 12 o!ers a theoretical base to apply "nite element techniques to analyse the contaminant migration in fractured media. Hence, the general form of the governing mathematical equation of contaminant transport in nonfractured media and in fractured media based on the theory of Rowe and Booker11,12 can be written as $ ' (D ' $c)!V ' $c"n ! ! Lc #g#q Lt (1) where c is the concentration of the contaminant in the pore water of non-fractured media or the concentration in the fractures of fractured media, D is the apparent hydrodynamic dispersion ! tensor, V is the vector of the components of the apparent groundwater velocity (Darcy velocity), ! $ is the gradient operator, n is the porosity of the soil in the non-fractured media or the fracture porsity in the fractured media, g is the sum of the rates at which contaminant is lost from the groundwater due to sorption onto the solid matrix and the e!ects of radioactive decay and biodegradation and q is the rate per unit volume at which contaminant is being transported into the matrix material by matrix di!usion in the fractured media or is omitted in the non-fractured media. The source term q in the Laplace domain has been established by Rowe and Booker11}13 and is given in the appendix. If the principle axes of hydrodynamic dispersion are chosen to be parallel to the co-ordinates axes, the contaminant transport equation for three-dimensional homogeneous non-fractured or fractured media under steady groundwater #ow condition is given by the following equation: A B L2c L2c L2c Lc Lc Lc Lc D #D #D !< !< !< "nR #j*c #q !xx Lx2 !yy Ly2 !zz Lz2 !x Lx !y Ly !z Lz Lz (2) where c is the concentration of the contaminant in the pore water of non-fractured media or the concentration in the fractures of fractured media, n is the porosity of the soil in the non-fractured media or the fracture porosity in the fractured media, R is the retardation factor (R"R "1#oK /n) in the non-fractured media or R"R "1#b (K /n ), in the fractured $ $ & & & & media in which b , K , and n are de"ned by Rowe and Booker),11, 12 j*, j*"j/R , in which j is & & & $ the decay constant, in the non-fractured media or j*"j/R in the fractured media, < , < , < & !x !y !z are the components of the apparent groundwater velocity (Darcy velocity) in the x, y, and z directions and D , D , D are the apparent hydrodynamic dispersion coe$cients. !xx !yy !zz 2.2. Transformations of the governing equations Both Fourier and Laplace transformations are utilized in the Fourier}Laplace transform "nite element method (FLTFEM). The Fourier integral transform is "rst applied to the governing Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1766 J. C. WANG AND J. R. BOOKER equations in order to make a reduction of the dimensionality of the contaminant transport problem by one. In applying such a transform, the assumption is made that the material properties of the soil and geometry do not vary in the y direction although they may change in the xz plane and then a Fourier transform = P~= c (x, y, z, t) e~*uy dy 1 = c" C (x, u, z, t)`*yu du 2n P ~= C" (3) (4) of the governing equation (2) becomes A B L2C L2C LC LC LC D !(u2 D #iu< ) C#D !< !< "nR #j* C #Q !xx Lx2 !yy !y !zz Lz2 !x Lx !z Lz Lt (5) where C is the Fourier transform of concentration (c), Q is the Fourier transform of q and u is the Fourier transform parameter. Following the application of a Fourier integral transformation, a Laplace transform is introduced: = P 0 C(x, u, z, t) e~st dt CM " (6) Then the contaminant transport equation in Fourier}Laplace transform space can be written as D L2CM L2 CM LCM L CM #D !< #< "(u2 D #iu< #h) CM !nRC !xx Lx2 !zz Lz2 !x Lx !z Lz !yy !y 0 (7) where s is the parameter of Laplace transform, CM is the concentration in Fourier}Laplace transform space, C is the initial concentration in Fourier transforms space, 0 = P~= c0 e~*uy dy C " 0 G h" n R (s#j*) in a non-fractured material $ n R (s#j*)#g6 (s#j* ) in a fractured material . f & in which g6 and j are given in the appendix. . It is obvious that the transformed problem in the Fourier}Laplace transform space becomes 2D problem provided that the boundary conditions can be expressed in terms of C. This will always be the case if the boundary conditions are of the same type along any line x"constant and z"constant. For the special case of 3D contaminant migration from a surface land"ll through a horizontally layered soil media in which the soil properties are the same at any horizontal location within the layer, if a second Fourier transform is taken along the x-axis, then the method reduces to the "nite layered method.14 Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1767 2.3. Boundary and initial conditions In order to completely de"ne the problem, it is necessary to specify the boundary and initial conditions. Although these conditions depend on the particular problems being treated, the ones which occur often will be described here. More complicated situations such as the speci"c #ux boundary condition treated in the section on leakage will be discussed in detail when the related problem is introduced. 2.3.1. Specixed concentration at boundary. As previously discussed the FLTFEM can be applied when the boundary conditions are of the same type and the material properties do not vary in one spatial dimension. A particular case of a boundary condition which is often considered in this paper is introduced here. Considering a surface land"ll in which the leachate in the land"ll has a constant concentration (c ) and zero concentration is assumed elsewhere. The schematic B diagram of this situation is shown in Figure 2. The boundary condition of this surface land"ll may be written as G c , c" B 0, D y D(H/2, D x D(¸/2, z"0 D y D'H/2, D x D'¸/2, z"0 Introducing a Fourier transform, the transformed surface boundary condition can be expressed as = P~= C " B H@2 P~H@2 c e~*uy dy" A B 2c uH c e~*uy dy" B sin B u 2 (8) Suppose that the concentration c is constant with time. It follows that taking a Laplace B transform of the boundary conditions in Fourier}Laplace transform space leads to G 2c B sin su CM " B 0, A B uH , 2 D x D(¸/2, z"0 D x D'¸/2, z"0 (9) where c is the constant concentration on the boundary in real space and CM is the transformed B B concentration in Fourier}Laplace transform space. Figure 2. Boundary condition of a surface land"ll Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1768 J. C. WANG AND J. R. BOOKER Figure 3. Schematic of an initially contaminated block If the material properties of the soil do not vary in the y direction, then the above transformed boundary conditions can be applied in conjunction with the FLTFEM to determine the contaminant migration from a 3D land"ll by solving a 2D time-independent transport (7). 2.3.2. Initial concentration in the contaminated zone. Considering the case of a block of contaminated soil (Figure 3) in which the initial concentration (c ) is assumed to be spatially constant 0 within the contaminated zone and zero initial concentration elsewhere; then the initial condition may be written as G c , D y D(H/2, D x D(¸ /2, D z D(¸ /2 x z c" 0 i 0 D y D'H/2, D x D'¸ /2, D z D'¸ /2 x z Taking the Fourier transform of this relation it is found that G 2c :H@2 c e~*uy dy" 0 sin (uH/2) when D x D(¸ /2; D z D(¸ /2 x z u ~H@2 0 CM " 0 0 when D x D'¸ /2; D z D'¸ /2 x z (10) where c is the initial concentration in real space and C is the initial concentration in Fourier 0 0 transform space. If the material properties of the soil do not vary in the y direction and the boundary conditions are of the same type in the y direction, then the above transformed initial conditions can be incorporated with equation (7) in the FLTFEM to examine the extent of contamination due to this initially contaminated soil 2.4. The xnite element scheme in transformed space The "nite element approximation of the transformed contaminant transport equation (7) together with its associated boundary and initial conditions of particular problems can be obtained by using the Galerkin "nite element technique. By using the method of weighted residuals in the "nite element scheme, the equation for the contaminant transport problem can be Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1769 ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA expressed as PP) C = (x, z) D L2 CI L2 CI L CI L CI #D !< !< !xx Lx2 !zz Lz2 !x Lx !z Lz D !(u2D #iu< #h) CI #n R C d)"0 !yy !y 0 (11) where = (x, z) are the weighting functions (i.e. interpolation functions) and ) represents the domain of a typical section. If CI denotes the approximate solution in Fourier}Laplace transform space, the approximate solution can be written as m n (12) CI (x, z, u, s)" + + Ne (x, z) CM e (u, s) i i e/1 i/1 where Ne (x, z) are interpolation functions, CM e (u, s) are the unknown nodal values of the "eld i i variable in Fourier}Laplace transform space, n is the number of nodes within the element, m is the number of elements, and the superscript e indicates a particular element. Using integration by parts and substituting equation (12) into equation (11), the intergral in equation (11) can be evaluated and leads to a system of linear equations which may be cast in matrix form as [M] MC1 N"MF N#MF N (13) * & where [M] is the global matrix including advection}dispersion, dispersion, sorption, decay matrices which may be obtained by combining the element matrices and is given by equation (14), [MC1 N is a vector containing the unknown nodal concentrations in Fourier}Laplace transform space, MF N is a vector that depends on the known initial conditions which is given by equation * (15), and MF N is a vector due to the prescribed #ux ( f ) which is the speci"ed outward normal #ux & n along the boundaries is given by equation (16). If no speci"ed #ux is encountered, the #ux is zero and this term will vanish: m [M]" + e/1 P m !+ e/1 m #+ e/1 C DC DC D C DC DC D L Ne LNe 1 1 Lx Lz F F L Ne LNe n n Lz Lx A% P L Ne LNe 1 1 Lx Lz F F L Ne LNe n n Lx Lz A% P Ae CD Copyright ( 1999 John Wiley & Sons, Ltd. De !xx 0 ve !x 0 De !zz 0 L Ne LNe 1 2 n Lx Lz dA L Ne LNe 1 2 n Lz Lz Ne 2 Ne n 1 dA 0 ve !z Ne 1 2 Ne n Ne 1 F [he#De u2#iu<e ] [Ne 2 Ne] dA !yy !y 1 n Ne n (14) Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1770 J. C. WANG AND J. R. BOOKER m MF N" + * e/1 P Ae m MF N" + & e/1 CD PC D Ne 1 F [ne Re Ce ] dA 0 Ne n (15) Ne 1 F [FM ] d! n Ne n (16) !e The above equations describe the "nite element approach for approximating the equation governing contaminant migration by using the FLTFEM. The initial condition for a contaminant transport problem is incorporated into the transport equation (7), when a Laplace transform is introduced and is then incorporated into the analysis through equations (13) and (15). The concentration is also often speci"ed on a portion of the boundary, this is easily incorporated into a "nite element analysis by specifying the transformed concentration at those boundary nodes. Alternatively, the normal #ux entering the boundary may be speci"ed. It is then necessary to transform the speci"ed normal #ux, this will be described in detail later. The condition of continuity of #ux and concentration at boundaries separating di!erent materials is implicit in the "nite element procedure. 2.5. Numerical inversion methods The solutions obtained directly from the "nite element method are in the Fourier}Laplace transform space. Hence, numerical inversion techniques are necessary to invert those solutions in transformed space into the solutions in real time and real space. The inversion of the Laplace transform can be performed using a convenient algorithm developed by Stehfest9,10 while in the case of Fourier transform inversion, a conventional numerical integration technique, Gauss}Legendre quadrature is used. The inversion of Laplace transform has been given in detail in References 9, 10 and 15. Hence, only the inversion of Fourier transform will be discussed here. The concentration of contaminant (c) is related to the transformed concentration (C) by equation (4) which is restated here for each of reference: 1 c" 2n = P~= C (x, u, z, t) e`*yu du (17) The integral in equation (17) may be evaluated numerically by the use of Gauss quadrature as follows. First the in"nite integral is truncated and then broken into a number of subintervals so that 1 c" 2n = P~= C (x, u, z, t) e`*yu du : 1 k/`p~1 + C (x, u, z, t) e`*yu du 2n K/~p (18) where ¸ is the length of the subinterval. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1771 Each of the integration on these subintervals is then approximated using Gaussian quadrature so that m C (x, u, z, t) e`*yu du : + = C (x, u , z, t) e`*yuK (19) K K k¸ K/1 where = are the Gauss weights, u are the Gauss co-ordinates and m is the total number of K K Gauss points. Details of Gauss}Lengendre quadrature may be found in References 16}18. Generally, using more Gauss points will produce more accurate results. In order to assure the accuracy of numerical integration technique the numerical solutions obtained in this paper was corroborated by using large integration range and more Gauss point. Experience in dealing with the contaminant migration problems considered here shows that both the range of integration and the number of Gauss points required to achieve a speci"ed accuracy of solution depend on the problem being considered and the time at which the solution is being evaluated. Generally, a larger range of integration and number of Gauss points will be required when the solution is evaluated at larger times. P (k#1)¸ 3. VERIFICATION The technique (FLTFEM) developed in this paper was veri"ed by evaluating four test examples including two 2D and two 3D cases of contaminant migration in the non-fractured media. In all test examples the soil is assumed to have the same geometry and material properties in the y direction and thus the Fourier transforms can be taken along the y-axis to reduce the dimensionality of problems being treated. As mentioned earlier, in the FLTFEM, a Gauss}Lengendre quadrature is used to invert the Fourier transform and the Stehfest inversion method is used to invert the Laplace transform. The numerical solutions obtained by FLTFEM were compared with the analytical solutions derived by Cleary and Ungs,19 Leij and Dane,20 and Booker.21 3.1. Two-dimensional contaminant transport problems The FLTFEM is "rst applied to contaminant migration under two-dimensional conditions. As was mentioned previously, through the application of FLTFEM, these 2D contaminant transport problems can be solved by using a 1D time-independent transformed mass transport equation in Fourier}Laplace space. 3.1.1. Example for constant concentration surface source. Suppose that a surface land"ll (Figure 4(a)) of a "nite width of ¸ is the y direction but with a relatively larger length in the x direction overlays a homogeneous, isotropic porous medium. Under these conditions, there will be no variation of concentration in the x direction. A undirectional steady #ow is also assumed to have a vertical seepage velocity < in the z direction and the values of the dispersion coe$cients in that ; direction and orthogonal to it are respectively denoted as D and D . In this example, the z y contaminant is assumed to be non-reactive and non-decaying (R "1 and j"0). The concentra$ tion of the contaminant is assumed to be spatially and temporally constant in the land"ll and initially zero elsewhere. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1772 J. C. WANG AND J. R. BOOKER Figure 4. Schematic diagram for 2D example of constant concentration surface source. (a) 2D problem; (b) Transformed 1D problem In order to be speci"c, the width of the land"ll source (¸) is assumed to be 100 m. The transport parameters D "1 m2/day, D "0.1 m2/day, and < "0.1 m/day are used. Since z y z the geometry and material properties of soil are assumed to be the same in the y direction, the Fourier transform can be taken along the y-axis in the numerical approach and this twodimensional contaminant transport problem is then reduced to a problem involving only one spatial dimension (see Figure 4(b)). A one-dimensional "nite element mesh consisted of 111 linear bar elements and each element of length of 2 m was used to analyse this two-dimensional problem. The transformed concentration, which was given by equation (9) in Section 2.3, is speci"ed at the boundary node of the ground surface because a constant concentration (c ) is assumed at the B surface boundary. A no-#ux condition was assumed at the lower boundary in the FEM. However, the remote boundary (222 m from the ground surface) was found to be deep enough to avoid any noticeable boundary e!ects. The numerical results of FLTFEM were compared with the analytical solution of Cleary and Ungs.19 By using FLTFEM, the concentration of a contaminant at particular positions and at speci"ed times can be quickly evaluated. In this case it was found that 60 G points, one subinterval and ¸"4 (the range of integration of each subinterval used in Gauss}Lengendre quadrature) was adequate to obtain an accurate solution because the solutions examined were found to be almost the same when even more Gauss points, a larger value of ¸ or more subintervals were used. Comparisons of FLTFEM numerical results and the analytical solutions are shown in Figure 5. It can be seen that numerical results of concentration distributions are nearly identical to those obtained analytically. 3.1.2. Example for an initially contaminated rectangular zone. The purpose of this example is to test the ability of the FLTFEM to model contaminant transport from an initially contaminated region. For this reason, groundwater contamination due to waste liquid seeping from the bottom of a waste repository into a homogeneous aquifer is considered. It is assumed that this gives rise to an initially contaminated area of a spatially constant concentration in the aquifer and that there is no migration above and below the aquifer (see Figure 6). An uniform steady groundwater #ow with seepage velocity < "1 m/y is also assumed in the x direction of the x-axis and the transport parameters are taken to be D "10 m2/y and x Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1773 Figure 5. Two-dimensional examples of constant concentration surface source Figure 6. A schematic showing a vertical section of the aquifer and the surface impoundment D "1 m2/y. Both the length and the width of the initially contaminated area are assumed to be y 10 m. A plan of the details are shown schematically in Figure 7(a). Since the material properties of the aquifer are assumed to be constant in the y direction the FLTFEM converts the 2D governing equation of contaminant transport into a 1D transport equation. Figure 7(b) shows the one-dimensional "nite element mesh of 240 linear bar elements, which have an equal length of 1 m, used to analyse this problem. A mesh of "nite extent has to be used since it is not possible to model an in"nite boundary, no #ux boundaries are assumed at the boundaries of "nite element domain. However the 1D "nite element mesh (total length of 240 m) used to analyse this problem was found to be su$ciently large to avoid any noticeable boundary e!ects. The treatment of the initial condition is incorporated into the analysis through equations (7), (10), (13) and (15). In this case, 60 Gauss points and two subintervals with ¸"4 (the range of integration of each subinterval) were found to be adequate in the Gauss}Lengendre quadrature since no signi"cant Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1774 J. C. WANG AND J. R. BOOKER Figure 7. Schematic diagram for 2D example of initial concentration rectangle source. (a) 2D problem; (b) Transformed 1D problem Figure 8. Two-dimensional example of initial concentration rectangle source di!erence was found in the solutions even when more Gauss points, a larger value of ¸ or more subintervals were used. Figure 8 shows the comparison of the numerical solutions obtained by the FLTFEM and the analytical solutions from Booker.21 It can be seen that the numerical results obtained using the FLTFEM provide a very good "t to the analytical solutions as well. 3.2. Three-dimensional contaminant transport problems The accuracy of the FLTFEM has been demonstrated by the accuracy of the numerical solutions it provided to the 2D contaminant transport test examples. On some occasions, Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1775 however, it is necessary to consider contaminant migration under three-dimensional condition. Therefore, it is desirable to extend the veri"cation of the FLTFEM to some 3D contaminant transport problems. As previously discussed, by applying the Fourier transform together with the Laplace transform, a transformed 2D time-independent mass transport equation can be used to solve some particular three-dimensional transient contaminant migration problems. This leads to signi"cant savings in computational cost. 3.2.1 3D case of a landxll with a leachate which has a constant concentration. In order to further verify the FLTFEM under 3D conditions, a waste disposal impoundment (land"ll) constructed overlaying a homogeneous soil, shown in Figure 9, is considered. It is assumed that the waste is continuously dumped into the land"ll and that the concentration within the land"ll remains constant over the period of interest. It is also assumed that the groundwater #ow is in the vertical direction and has a constant seepage velocity < . This 3D contaminant transport problem has z been solved analytically by Leij and Dane.20 The contaminant is also assumed to be non-reactive and non-decaying (R "1 and j"0). The $ Fourier transform is taken along the y-axis as it was in the previous 2D problems. Thus this three-dimensional consideration of contaminant transport from a constant concentration land"ll source can then be reduced to a problem involving only two spatial dimensions (shown in Figure 10). Figure 9. A surface land"ll Figure 10. Schematic diagrams of 3D case of a land"ll with constant concentration condition and the transformed 2D problem. (a) 3D real problem; (b) Transformed 2D problem Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1776 J. C. WANG AND J. R. BOOKER The transformed 2D problems is symmetric about the z-axis. Therefore, a simplifed half plane (xz) domain (Figure 10(b)) was used to analyse this problem and a no-#ux boundary was assumed as the z-axis (x"0 m) because the problem is symmetric about it. It was assumed that there was a constant concentration (c ) in the land"ll because the land"ll contained a very large amount of B contaminant. Zero concentration was also assumed on the ground surface since it was assumed that periodic rainfall #ushed any surface contamination away. Since it is not possible to model an in"nite boundary a mesh of "nite extent was used. A zero #ux was assumed at the other boundaries of the "nite element domain. In order to be speci"c, the transport parameters D "1 z m2/day, D "D "0.1 m2/day, and < "0.1 m/day were used in this 3D test example and both x y z ¸ and ¸ were assumed to be 10 m. x y In the "nite element approach, a two-dimensional rectangular "nite element mesh with total number of 300 linear quadrilateral elements (Figure 11) was used. The mesh was re"ned at the edge of the land"ll in order to model the high concentration gradients which occur there. The transformed concentration (CM ) which was discussed in Section 2.3 was speci"ed at the boundary B nodes of land"ll (z"0) and is found to be G A B 2c u¸ B sin y su 2 CM " B 0 when 0(x(5 m elsewhere No perceptible boundary e!ects were found during the period of interest indicating that the positioning of the boundary was satisfactory. In this case, the solutions obtained by using 60 G points, two subintervals and ¸"4 (the range of integration of each subinterval) in the Gauss}Lengendre quadrature were found to be satisfactory because the solutions obtained were virtually the same when more Gauss points, a larger value of ¸ or more subintervals were used. Comparisons of the numerical result obtained by using the FLTFEM and the analytical solution are shown in Figures 12(a)}14(a). Excellent agreement is observed between the numerical and Figure 11. The "nite element mesh used in FLTFEM for 3D constant concentration surface source Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1777 Figure 12. Three-dimensional example of constant concentration surface source at y"5 m plane. (a) Concentration distribution curves at di!erent values of x; (b) concentration contours at y"5 m plane (time"100 days) analytic solutions at di!erent values of y. Figures 12(b)}14(b) also show the concentration contours of half of the xz plane but at di!erent values of y. Inspecting these "gures, it can be seen that the concentration pro"les at diferent values of x are obvious di!erent. This is because of land"ll was considered to be of "nite length in the x direction. A signi"cant di!erence is also observed at the same value of x (x"4 or 5 m) but di!erent values of y (i.e. y"5 m in Figure 12, y"0 m in Figure 13 and y"6 m in Figure 14) since the land"ll considered also has a "nite length in the y direction. 3.2.2. Case of a block of initially contaminated soil. The next example considers the extent of contamination due to a three-dimensional block of contaminated soil which was assumed to be at Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1778 J. C. WANG AND J. R. BOOKER Figure 13. Three-dimensional example of constant concentration surface source at y"0 m plane. (a) Concentration distribution curves at di!erent values of x; (b) concentration contours at y"0 m plane (time"300 days) a certain depth (h) beneath the ground surface in a deep homogeneous soil (see Figure 15). A spatially uniform initial concentration was assumed within this contaminated block and the initial concentration is assumed to be zero elsewhere. It is also assumed that there is zero concentration on the surface of ground because periodic rainfall #ushes any surface contamination away. The dimension of the contaminated soil and the speci"c values of parameters used are shown in Figure 15. A uniform steady groundwater #ow of seepage velocity < "10 m/y was assumed in x the x direction and the values of the dispersion coe$cients adopted were D "10 m2/y and xx Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1779 Figure 14. Three-dimensional example of constant concentration surface source at y"6 m plane. (a) Concentration distribution curves at di!erent values of x; (b) concentration contours at y"6 m plane (time"300 days) D "D "1 m2/y, respectively. The length (¸ ) and width (¸ ) of this initially contaminated soil yy zz x y are 10 m and the height (¸ ) of it is assumed to be 5 m. z By introducing a Laplace transform and a Fourier transform which was taken along the y direction, this 3D problem was then transformed into a simpler 2D time-independent problem which is shown schematically in Figure 16. A two-dimensional "nite element mesh shown in Figure 17 with a total number of 528 linear quadrilateral elements was used. The nodal concentrations are zero at the surface and a no-#ux condition was used at the other boundaries of "nite element domain. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1780 J. C. WANG AND J. R. BOOKER Figure 15. Schematic diagram of 3D case of a block of contaminated soil with initial concentration condition Figure 16. 2D domain solved in FLTFEM for a 3D problem of a block of initially contaminated soil The treatment of initial condition has been discussed previously and the transformed initial concentration (C ) in Fourier transform space, is given below: 0 2c u¸ 0 sin y when !¸ /2(x(¸ /2,!¸ /2(z(¸ /2 x x z z 2 C " u 0 0 elsewhere G A B This condition is incorporated into the transport equation (7) when a Laplace transform is introduced and is then incorporated into the analysis through equations (13) and (15). Numerical results shows that using 60 G points, three subintervals and ¸"4 (the range of integration of each subinterval) are su$cient to obtain an accurate solution in this case considering that nearly identical solutions were obtained when more Gauss points, a larger value of ¸ or Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1781 Figure 17. The "nite element mesh used in FLTFEM for a 3D problem of a block of initially contaminated soil Figure 18. Three-dimensional problem of a block of initially contaminated soil more subintervals were used. It is also found that there is no noticeable e!ect due to the zero-#ux boundaries used. Comparisons of numerical results and analytical solutions from Booker21 are given in Figures 18}20. It can be seen that at di!erent times and di!erent values of y the simulation results obtained by the FLTFEM are in very close agreement with those obtained analytically. 4. LEAKAGE FROM A WASTE REPOSITORY Often, leachate migration from sanitary land"lls is an important source of groundwater pollution. As well, other types of land disposal units including domestic septic tanks, waste disposal through wells and deeply buried were repository, etc., are potential sources of contamination. Thus, it is Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1782 J. C. WANG AND J. R. BOOKER Figure 19. Three-dimensional problem of a block of initially contaminated soil Figure 20. Three-dimensional problem of a block of initially contaminated soil Figure 21. Schematic diagram of leakage problem in a waste repository necessary to be able to predict the likely impact of the contamination due to such buried water repository and to assess if remediation action is required. Suppose leakage is observed from a long deeply buried cylindrical pipeline (shown schematically in Figure 21) which was constructed for conveying chemical liquid or waste substances. The Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1783 long cylindrical pipe is assumed to be buried deep in the non-fractured porous media. It is assumed that the pipe develops a break of width (H) at one of the pipe joints and that a constant rate of contaminant mass #ux ( f ) is released from the circumference of the pipeline into the 0 surrounding soil. Since the pipe is assumed to be very long and that the material properties of soil do not vary in one spatial dimension (the y direction in this case) thus the FLTFEM can be applied to assess the impact of contamination due to this leakage. Two cases of leakage from this cylindrical pipeline are investigated in this section. The "rst one is the case of leaking of the contaminant from the pipe occurring over a long period (¹ PR). $ The second is the case of leakage which occurs over a period of "nite duration (¹ ). The governing $ equations and boundary condition for the problem are described below. The governing equations: The governing equation is equation (2). As in the cases considered previously, the assumptions that the same geometry and material properties along the y-axis are made and then, by introducing the Fourier and Laplace transforms, the transport equation in transformed space is given by equation (7). Boundary condition at the section of leakage. The mass #ux of contaminant #owing out from the cylindrical pipeline is assumed to be temporally and spatially constant during the period (¹ ) of $ leaking. In other words, there is a constant rate of #ux ( f ), for the duration ¹ , over the section 0 $ (H) where the leakage occurs (see Figure 21). This boundary condition can be expressed as: G f" f 0 0 D y D(H/2, r"R and if 0(t(¹ over the section of the leakage $ Dy D'H/2, r"R where r is the radial distance from the centre of the pipe (r"Jx2#z2), and R is the radius of the pipe. The contaminant #ux on the circumference of the pipe, in Fourier}Laplace transform space, then becomes A B 2f uH FM " 0 sin (1!e~sT$) su 2 (20) where FM is the radial contaminant mass #ux in Fourier}Laplace transform space, f is the radial 0 contaminant mass #ux in real space, s is the Laplace transform parameter, u is the Fourier transform parameter, and ¹ is the duration of leakage. $ If the leaking of the contaminant from the pipe occurs over a long period (¹ PR), the $ contaminant #ux emanating from the circumference of the pipe in Fourier}Laplace transform space is given by A B 2f uH FM " 0 sin su 2 (21) 4.1. Case 1: leakage occurred over an inxnite duration In this case, it is assumed that a long cylindrical pipeline having a diameter of 2 m is constructed in an idealized isotropic clayey soil and that is develops a break of 0)2 m wide. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1784 J. C. WANG AND J. R. BOOKER A constant rate of contaminant mass #ux ( f "1000 g/m2/a) is released into the surrounding soil 0 over the section where the leakage occurs. The lacking of contaminant is assumed to occur over an in"nite duration (¹ "R). $ Also it is assumed that there is a uniform advection velocity (< "!0)02 m/a) in the x x direction and the hydrodynamic dispersion coe$cient is 0)05 m2/a in all directions because the groundwater velocity in clayey soil is relatively small. The porosity of soil is assumed to be 0)4 and the contaminant is conservative (K "0) and does not decay (j"0). $ The material properties of soil and the geometry considered do not vary along the y direction and thus a Fourier transform can be taken along the y-axis. This enables the three-dimensional problem of leakage from a long cylindrical pipe to be transformed into a 2D time-independent problem which can be evaluated using a 2D "nite element mesh. Figure 22 shows the actual 3D situation and the 2D cross-section of the domain. A half-plane 2D "nite element domain was used in the FLTFEM. The 2D problem is symmetric about the x-axis and thus a non-#ux condition was assumed at the boundary z"0 m. The 2D "nite element mesh of 440 linear quadrilateral elements shown in Figure 23 was used to analyse this problem. Figure 22. Schematic diagram of 3D leakage problem and its 2D domain used in the FLTFEM Figure 23. The 2D "nite element mesh used in FLTFEM for 3D leakage problem from a cylindrical waste repository Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1785 A constant #ux ( f ) was released into soil from the circumference of the cylindrical pipe. Thus, it is 0 necessary to transform the constant #ux into the nodal #uxes through equations (16) and (21) these nodal #uxes are applied at the boundary nodes of pipe. A no #ux boundary condition was also adopted at the other boundaries in this case since it is not possible to model an in"nite domain. The simulation results were examined to ensure that there was no noticeable boundary e!ect due to the "nite size of the mesh adopted. In this case a satisfactory solution was obtained by using 100 G points, 1 subinterval and ¸"100 (the range of integration of each subinterval) in the FLTFEM since no signi"cant di!erence in the solution was found by using a larger value of ¸, more Gauss points or more subintervals. The concentration distributions in the radial direction at di!erent angles for di!erent values of y are shown in Figures 24 and 25. The semi-analytical solutions for this problem developed by Booker and Leo8 are also presented in the "gures. It can be seen the numerical solutions are almost identical to those of semi-analytical solutions. Inspection of the simulation results (Figures 24 and 25) shows that the highest peak concentration is found in the direction of h"03 (the upstream side of pipe) and followed by the h"903 (the top side of pipe) at the two cross-sections (y"0 and 0.1 m). The smallest peak concentration is Figure 24. Solutions of leakage problem in a waste repository for in"nite leaking duration (¹ "R) at y"0 m plane $ Figure 25. Solutions of leakage problem in a waste repository for "nite leaking duration (¹ "R) at y"0.1 m plane $ Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1786 J. C. WANG AND J. R. BOOKER found along the direction of h"1803 (the downstream side of pipe) which has the same direction as that (negative x) of groundwater #ow. This is because the migration of the contaminant in the opposite direction of the #ow (along the direction of h"03) was retarded by the groundwater #ow and the contaminant is taken away along the direction of h"1803 from the pipe. Thus, the contaminant is transported further away from the waste repository in the direction of h"1803 (the downstream side of pipe) than in the direction of h"03 (the upsteam side of pipe). No matter what direction (h) is considered, the peak concentration is always found at the boundary of pipe. 4.2. Case 2: leakage occurred over a period of xnite duration Now, exactly same leakage problem as discussed above (see also Figure 22) was again considered, except that the leaking of contaminant is now assumed to occur at a constant rate of #ux ( f "1000 g/m2/a) but over a "nite period of time (¹ "10 y). The transport parameters 0 $ adopted in this case are identical to those used in the previous case and the "nite element mesh shown in Figure 23 was also used. As would be expected, the simulation results at 10 y time are found to be exactly same as the results obtained in the previous case because during this period (10 y) the two situations are identical. The solution was evaluated at both 12 and 20 y which are some time after the leaking ceased and both numerical and semi-analytic results are given in Figures 26 and 27. Similar behaviour to that found in the previous case is also observed in di!erent radial directions. The highest peak concentration is found in the direction of h"03 (the upstream side of pipe) and followed by the h"903 (the top side of pipe). The smallest peak concentration is found along the direction of h"1803 (the downstream side of pipe) which has the same direction as that of groundwater #ow. It is also observed that the contaminant is transported further away in the direction of h"1803 (the downstream side of pipe) and the migration of the contaminant in the opposite direction of the #ow is retarded by the groundwater #ow. After an additional 10 y (t"20 y), Figures 26 and 27 show that the concentration of the contaminant becomes smaller as time increases and the peak concentration in the direction of h"1803 (the downstream side of pipe) is found not to be at the boundary of pipe. This is because Figure 26. Solutions of leakage problem in a waste repository for "nite leaking duration (¹ "10 yr) at y"0 m plane $ Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1787 Figure 27. Solutions of leakage problem in a waste repository for "nite leaking duration (¹ "10 yr) at y"0.1 m plane $ leaking of the contaminant ceased 10 y previously and after that time the concentration at the boundary of the pipe starts to decrease since contaminant is transported away from the pipe and is not replenished by the leak. Based on the results presented, it is evident that good agreement is obtained between FLTFEM solutions and the semi-analytical solutions of Booker and Leo8. It has been demonstrated by this case that the FLTFEM can also be used to analyse a contaminant transport problem which involves a change of the environmental conditions by using the multi-stage approach.22 5. 3D WASTE REPOSITORY IN FRACTURED MEDIA The migration of contaminant in non-fractured media under two- or three-dimensional situation has been discussed in previous sections. Sometimes the migration of a contaminant occurs in a "ssured porous media. For example, if a waste disposal repository was constructed overlaying a "ssured media, it may cause contamination of the "ssured material beneath it sometime after its construction because of aging, design shortcomings or acts of nature. Therefore, a 3D waste disposal repository placed over a fractured material, shown schematically in Figure 28, is examined in this section by using FLTFEM. Although, in many important practical situations such as land"ll, a simpli"ed 2D or 1D model is often su$ciently accurate for practical purposes because the longitudinal dimension of the contaminant source under consideration is often relatively large compared with the other dimensions or the depth of penetration of contaminant from a land"ll into the soil beneath it over the period of interest. It is also interesting to know what the di!erence will be between a 3D analysis and a simpli"ed 2D or 1D analysis of a surface waste repository. It is assumed that waste is dumped into a waste repository, shown in Figure 28(a), continuously and that this leads to a spatially uniform and temporally constant concentration of contaminant within the repository. A constant concentration (c ) was thus assumed in this waste repository. B Initially, it is assumed that there is no contamination of the "ssured material. Also zero concentation was assumed on the ground surface because it was assumed that periodic rainfall #ushed any surface contamination away. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1788 J. C. WANG AND J. R. BOOKER Figure 28. Schematic diagrams of a waste repository under 3D and 2D considerations. (a) 3D surface waste disposal repository; (b) 3D real problem; (c) 2D problem Figure 29. 2D fracture system The fracture system begin considered is shown in Figure 29 and the material properties of fractured media are listed below. two sets of fractures: set 1 and 3 fracture spacings: H "H "0)05 m 1 3 h "h "0)001 m 1 3 < "n < "0 m/a !x pxx x < "n < "0 m/a !y pyy y < "n < "0)01 m/a !z pzz z fracture openings: Darcy velocity: Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA Dispersivity: a "1 m, a "0)1 m L T fracture di!usion coe$cient: D "0)01 m2/a o matrix di!usion coe$cient: D "0)0036 m2/a . matrix porosity: n "0)1 . 1789 The 3D model of this waste repository is shown in Figure 28(b). A simpli"ed 2D model in the xz plane shown in Figure 28(c) and an idealized 1D model were also used to analyse this problem. The FLTFEM was used to implement the 3D analysis while the numerical model (LTFEM) developed by Booker and Wang23 for the fractured porous media was used to analyse the pollutant migration under 2D and 1D conditions. The problem domain (see Figure 28(c), and 2D "nite element mesh (see Figure 30) which has 1800 linear quadrilateral elements were used in LTFEM for the analysis of contaminant migration under 2D condition and a 1D mesh consisted of 100 linear bar elements each of length of 1 m was used for the analysis under 1D condition. The material properties do not vary in the y direction, so by applying Fourier transform along the y-axis, the 3D contamination problem can be treated as a 2D problem. Thus, the exactly same problem domain (see Figure 28(c)) and the 2D "nite element mesh (see Figure 30) can be used in FLTFEM for the analysis of contaminant migration under 3D conditions. As mentioned above, a contant concentration was assumed in the waste repository and zero concentration was assumed on the ground surface. A no-#ux condition was assumed at the other boundaries of the "nite element domain. However, when the numerical results were examined to perceptible boundary e!ects were found. The simulation results show that 60 G points, "ve subintervals and ¸"4 (the range of integration of each subinterval) are satisfactory to obtain an accurate solution because no signi"cant di!erence was found in the solution by using more Gauss points, a larger value of ¸, or more subintervals. Figure 30. Finite element mesh Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1790 J. C. WANG AND J. R. BOOKER Figure 31. Concentration contours at time"50 yr of 1D analysis Figure 32. Concentration contours at time"50 yr of 2D analysis The concentation contours at time"50 years for the analysis using both 2D and 1D models are presented in Figures 31 and 32. Also Figures 33 and 34 show the concentration contours obtained for the analysis using the 3D model at di!erent locations of y"0 m, (the centre of the waste repository) and y"5 m (the edge of the waste repository). Figures 35 and 36 show the concentration contours at time"100 y obtained by using both 2D and 1D models. Also the concentration contours obtained by the use of the 3D model at di!erent locations of y"0 m (the centre of the waste repository) and y"5 m (the edge of the waste repository) are presented in Figures 37 and 38. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1791 Figure 33. Concentration contours at y"0 m (¹"50 yr) of 3D analysis Figure 34. Concentration contours at y"5 m (¹"50 yr) of 3D analysis Examining Figures 31}33 shows that the contaminant plume at the plane of y"0 m obtained by using the 3D analysis is very close to the results obtained by using the 2D analysis and the simulation results of the contaminant migration along the depth by the use of 1D model is also close to those of the analyses of using 2D and 3D models at the position near the centre of the waste repository. This indicates that the extent of the contaminant beneath the centre line (y"0 m) of the waste repository may be analysed by using a simpli"ed 2D model or even using a simple 1D model without sacri"cing the accuracy in this case at 50 yr. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1792 J. C. WANG AND J. R. BOOKER Figure 35. Concentration contours at time"100 yr of 1D analysis Figure 36. Concentration contours at time"100 yr of 2D analysis However, examination of the contaminant plume on the plane of y"5 m (the edge of the waste repository) shows that there are signi"cant three-dimensional e!ects. The edge e!ects can be seen when comparing results of the 1D and 2D analyses (see Figures 31 and 32). As would be expected, investigation of the simulation results (shown in Figures 35}38) at a larger time (100 y) of these three models shows that the di!erence between a 3D model and a simpli"ed 2D and 1D model becomes greater at locations further from the centre of the waste repository and the edge e!ects become more signi"cant as time increases (see Figures 32 and 36 or Figures 33 and 37). Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA 1793 Figure 37. Concentration contours at y"0 m (time"100 yr) of 3D analysis Figure 38. Concentration contours at y"5 m (time"100 yr) of 3D analysis 6. CONCLUSIONS The Fourier}Laplace "nite element method (FLTFEM) developed in this paper has been demonstrated to be an accurate method for analysing some particular mulit-dimensional contaminant migration problems such as those involving a surface waste disposal repository and a long prismatic waste repository in which the geometry and material properties do not vary along one spatial dimension. The analysis uses a Fourier transform to reduce the dimensions of the multi-demensional contaminant transport problem by one, while using Laplace transform to Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1794 J. C. WANG AND J. R. BOOKER simplify a transient contaminant transport equation to a time-independent di!erential equation involved only spatial derivatives. The "nite element scheme is used to approximate the transformed governing equations. Once the solutions in Fourier}Laplace transform space are obtained, the solutions in real space and time can be attained through the numerical inversion techniques. The inversion of Laplace transform can be performed using a convenient algorithm developed by Stehfest 9,10 while in the case of the inversion of Fourier transform, a conventional numerical integation technique, Grass}Lengendre is used. With the combination use of Fourier and Laplace transforms, the Fourier}Laplace transform "nite element method (FLTFEM) greatly simpli"es a 3D contaminant transport problem and this certainly results in a computational bene"t. APPENDIX: EXPRESSION OF THE LAPLACE TRANSFORM OF q AND g FUNCTION The Laplace transform of the function g is used in the theory of Rowe and Booker11,12 to simulate the di!usive transport of contaminant from the fractures into the solid matrix in the fractured media. The expression for g6 and qN , which is the Laplace transform of the rate (q) at which the contaminant migrates from the fractures into the solid matrix per unit volume of fracture-matrix system, are brie#y introduced below. Considering that the molecular di!usion is the dominating process in the solid matrix and the advection is neglected, it is found that the concentration (c ) in the matrix satis"ed the equation . Lc n D +2c "(n #o K ) .#n j c . .& . . . . Lt m . . (22) where c is concentration in the solid matrix, D is the di!usion coe$cient of the solid matrix, n . .& . is the porosity of the solid matrix, o is the dry density of the solid matrix, K is the distribution . . coe$cient of the solid matrix, and j is the sum of the radioactive decay and biological . degradation constants in the solid matrix. The concentration of the surface of the matrix block will be equal to that in the "ssure system (c ), so & G c on the surface of the solid matrix c " & . 0 initially By taking Laplace tansform, equation (22) becomes n D +2 cN "n R (s#j* ) cN . . . .& . . . (23) where R "1#o K /n is the retardation coe$cient of the solid matrix, j* "j /R and . . . . . . . cN "cN on the surface of the matrix. . & Now, the Laplace transform of q may be de"ned as : FM dS . qN "! ./ : d< . (24) where !F denotes the normal component of #ux entering the matrix, : FM dS is the integral ./ ./ . along the surface S of the fracture, and : d< is an integral over the volume < of the solid . . . Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1795 ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA matrix. Using the divergence theorem, it is found : n D +2 cN d< n R (s#j* ) : cN d< . ." . . . . . qN " . .& (25) : d< : d< . . It was shown by Rowe and Booker11,12 that the Laplace transform of q can always be expressed in the form (26) qN "(s#j* ) g6 cN . & Referring to the equations (A.4) and (A.5), it is found that the Laplace transform of g is given by : cN d< . . g6"n R (27) . . : cN d< & . Case 1: A single set of xssures: Assuming there is only one set of "ssures (set 1) and considering di!usive transport into the matrix in the x direction, the concentration within the matrix satisi"es the equation L2cN ."n R (s#j* ) cN n D . .& Lx2 . . . . (28) where cN "cN . & x"$H 1 It it thus found that cosh kx cN "cN . & cosh kH 1 (29) tah kH 1 kH 1 (30) and g6"n R . . where k2"R (s#j* )/D . . .& . Also g6 may be expressed as a series expansion, C D 1 = (s#j*) m g6"n R 1!2 + (31) . . (s#j* #a2 (D /R ) (a H )2 i .& . i 1 . i/1 Case 2. Two sets of xssures: Assuming there are two sets of "ssures (sets 1 and 3) and considering di!usive transport matrix in the x and z directions, Laplace transform of g is given by Rowe and Booker.11, 12 C D 1 1 = = (s#j* ) . g6"n R 1!4 + + . . (s#j* #(a2#c2) (D /R ) (a H )2 (c H )2 k .& . i . i 1 k 3 i/1 k/1 where a "(i!1/2)n/H and c "(k!1/2) n/H . i 1 k 3 Copyright ( 1999 John Wiley & Sons, Ltd. (32) Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999) 1796 J. C. WANG AND J. R. BOOKER Case 3. Three sets of xssures: Finally, if there are three sets of "ssures (sets 1, 2 and 3) and it is observed that the transport distances of interest along the fracture system are large compared to the dimensions of the block, then Laplace transform of g is found by Rowe and Booker:11,12 C D 1 1 1 = = = (s#j* ) . g6"n R 1!8 + + + . . [s#j* #(a2#b2#c2) (D /R )] (a H )2 (b H )2 (c H )2 . i j k .& . i 1 j 2 k 3 i/1 k/1 j/1 (33) where a "(i!1/2)n/H , c "(k!1/2) n/H and b "( j!1/2)n/H . i 1 k 3 j 2 REFERENCES 1. J. A. Liggett and L. F. Liu, ¹he Boundary Integral Equation Method for Porous Media Flow, George Allen & Unwin., U.K., 1983. 2. H. T. Chen, T. M. Chen and C. K. Chen, &Hybrid Laplace transform/"nite element method for one-dimensional transient heat conduction problems' Comput. Meth. Appl. Mech. Engng., 63(1), 83}95 (1987). 3. Y. S. Yu, M. Heidari and W. T. Wang, &A semi-analytical method for solving unsteady #ow equations', Paper presented at Proc. ASCE National Conf. on Hydraulic Eng., Am. Soc. Civ. Eng., Williamsburg, 1987. 4. H. T. Chen and C. K. Chen &Application of hybrid Laplace transform/"nite di!erence method to transient heat conduction problems', Numer. Heat ¹ransfer, 14, 343}356 (1988). 5. E. A. Sudicky, &The Laplace transform Galerkin technique: a time-continuous "nite element theory and application to mass transport in groundwater', =ater Resour. Res., 25(8), 1833}1846 (1989). 6. C. K. Chen, T. M. Chen and J. W. Cleaver, &New hybrid laplace transform/"nite element method for two-dimensional transient heat conduction problem', Numer. Heat ¹ransfer, 20, 191}205 (1991). 7. G. J. Moridis and D. J. Reddel, &The Laplace transform "nite di!erence method for simulation of #ow through porous media', =ater Resour. Res., 27(8), 1873}1884, 1991. 8. J. R. Booker and C. J. Leo, &A FLT boundary integral equation method for analysis of contaminant leakage', in Valliappan, Pulmano, Tin-Loi (eds), Computational Mechanics- from Concepts to Computations, Proc. 2nd AsianPaci,c Conf. on Computational Mechanics, Sydney, Australia, 1993. 9. H. Stehfest, &Algorithm 368, Numerical inversion of Laplace transforms', Commun. Assoc. Comput. Mach. 13(1), 47}49 (1970). 10. H. Stehfest, &Remark on Algorithm 368, Numerical inversion of Laplace transforms', Commun. Assoc. Comput. Mach. 13(10), 624}625 (1970). 11. R. K. Rowe and J. R. Booker &A semi-analytic model for contaminant migration in a regular two- or threedimensional fractured network: conservative contaminants', Int. J. Numer. Anal. Meth. Geomech., 13(5), 531}550 (1989). 12. R. K. Rowe and J. R. Booker &Contaminant migration in a two- or three-dimensional fractured network: reactive contaminants', Int. J. Numer. Anal. Meth. Geomech., 14(6), 401}425 (1990). 13. R. K. Rowe and J. R. Booker, &Modelling of contaminant movement through fractured or jointed media with parallel fractures', Proc. of 6th int. Conf. on Numerical Methods in Geomechnics, 1988, pp. 855}862. 14. R. K. Rowe and J. R. Booker, &A "nite layer technique for calculating three-dimensional pollutant migration in soil', Geotechnique, 36(2), 205}214 (1986). 15. J. C. Wang and J. R. Booker, &A Laplace transform "nite element method (LTFEM) for the analysis of contaminant transport in porous media', Research Report (R756), Department of Civil Engineering, University of Sydney, 1997. 16. M. Abramowitz and I. A. Stegun, &Handbook of Mathematical Function, Dover Publications, New York, 1965. 17. B. Carnahan, H. A. Luther and J. O. Wilkes, 0Applied Numerical Methods', Wiley, New York, 1969. 18. A. Ralston and P. Rabinowitz, A First Cource in Numerical Analysis, 2nd edr, McGraw-Hill, New York, 1978. 19. R. W. Cleary and M. J. Ungs &Groundwater pollution and hydrology, mathematical modles and computer programs', Rep. 79-=R-15, Water Resour. Program, Princeton Univ., Princeton, NJ, 1978. 20. F. J. Leij and J. H. Dane, &Analytical solutions of the one-dimensional advection equation and two- or threedimensional dispersion equation', =ater Resour. Res., 26(7), 1475}1482 (1990). 21. J. R. Booker, &Analytic and semi analytic methods for the analysis of contaminant migration in soils', Proc. 8th Int. Conf. on Computer Methods and Adcances in Geomechanics, Morgantown, West Virginia, USA, 1994, pp. 3}19. 22. J. C. Wang, &A study of contaminant transport in porous media', Ph.D. ¹hesis, Department of Civil Engineering, University of Sydney, 1997. 23. J. R. Booker and J. C. Wang &A Laplace transform "nite element method (LTFEM) for the analysis of contaminant transport in fractured porous media', Int. J. Numer. Anal. Meth. Geomech., (1999) submitted. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)

1/--страниц