Modeling of Density-Dependent Flow and Transport in Anisotropic Aquifers: An Example Site of Landfill

In advection-dispersion transport processes, the density contrast between the contaminant fluid and the groundwater in the surrounding flow domain can affect the groundwater flow. Under such circumstances, the migration patterns and concentration distributions, of contaminant plumes can be changed. The primary purpose of this research is to simulate the problems of density-dependent groundwater flow and transport in homogeneous, anisotropic and inhomogeneous (or heterogeneous), anisotropic aquifers by means of a Galerkin finite element method. Secondly, this research is to show the numerical model of Galeati et al. (1992) can be applied to a landfill site. The leachate from a landfill site into the groundwater can change the contaminant concentration. Thus, there exists a density contrast between the contaminant fluid and the groundwater in the surrounding flow domain. In simulation, an example site of landfill is demonstrated. The results for contaminant sources from a landfill conclude as follows. If there does not exist a significant density contrast, the contaminant plume will spread in a shallow zone close to the water table; on the other hand, if there is such a contrast, the plume will sink downward into the groundwater flow system, and its migration pattern will only be changed after a longer transport time of six years. In a shorter time (one to three years), however, the migration pattern will not be significantly changed.

means of a Galerkin finite element method . Secondly, this research is to show the numerical model of Galeati et aL (1992) can be applied to a landfill site.
The leachate from a landfill site into the groundwater can change the conta1ninant concentration. Thus, there exists a density contrast between the contaminant fluid and the groundwater in the surrounding flow domain.
In simulation, an example site of landfill is demonstrated. The results for conta1ninant sources from a landfill conclude as follows. If there does not exist a significant density contrast, the contaminant plume will spread in a shallow zone close to the water table; on the other hand, if there is such a cont1·ast, the plume will sink downward into the groundwater flow system, . and its migration pattern will only be changed after a longer. transport time of six years. In a shorter time (one to three years), however, the migration pattern will not �e significantly changed.

. INTRODUCTION
210 TAO, Vol.6, No.2, June 1995 (4) the dispersivity, and (5) the aquifer properties (e.g., porosity, storage, hydraulic conduc tivity contrasts of aquifer materials, etc.). This research provides a more accurate prediction of the migration pattern of a contaminant plume based on the causative factors mentioned above. In detem1ining the depth and location of monitoring wells for groundwater quality, an accurate prediction of migration patterns in the vertical and horizontal section coordinates ( x, z ) and ( x, y ) in the modeling solute transport is required.
In the published literature, most density-dependent solute transport models have consid ered seawater intrusion into an aquifer (e.g., Henry, 1964;Segol et al., 1975;Kuiper, 1983;Voss and Souza, 1987;Huyakom et al. 1987;Dorgarten and Ts ang, 1991;and Galeati et al., 1992). Frind ( 1982) applied the seawater intrusion concept to simulate the solute transport of leachate resulting from a landfi ll site. Freeze and Cherry (1979) indicated that if the contaminant solution entering the flow regime has the same density as the groundwater, the contaminant plume will spread in a shallow zone close to the water table. If the contaminant solution is considerably denser than the groundwater, the plume will sink downward into the groundwater system. . Recently, Galeati et al. ( 1992) provided a more advanced continuity equation based • on different assumptions for the fluid from the one by Frind (1982). However, Galeati et al. (1992) applied their model of density-dependent transport to seawater intrusion into an aquifer, not to a landfill site. This research simulates the problems of density-dependent groundwater flow and transport in homogeneous, anisotropic and inhomogeneous, anisotropic aquifers for a landfi ll site by means of a Galerkin finite element method. An example site of landfill is demonstrated for simulation. The purpose of this research is to show the numerical model of Galeati et al. ( 1992) is applicable to a landfill site.

1 Fluid Continuity Equation
The continuity equation for the unsteady-state fl uid used in this research is the same as that used for the problem of seawater intrusion by Galeati et al. (1992). It is written as: Galeati et al. (1992) is more advanced than the one by Frind (1982). The second term on the right side of (1) which is due to the time derivative of the liquid density does not exist in the equation by Frind ( 1982).
With the density contrast (or density difference ratio) defined as: the constitutive relation is written as: where K i j is the hydraulic conductivity tensor, x j (j=l,2,3) are Cartesian coordinates (with I and 2 standing for the horizontal and 3 for the vertical directions), h is the hydraulic head, T/j is a vector such that 7J;=O along the horizontal direction and 7J ; =l along the vertical, / is the density contrast (or density difference ratio), S s is the specific storage, t is the time, </> is the porosity, Q is a fluid source tertn, po freshwater density, P max is the density corresponding to the maximum concentration, and C is a dimensionless concentration (actual divided by maximum). In (1), it is assumed that the conditions are isothe11nal, the fluid i.s incompressible and the salt does not undergo chemical reactions. It is also assumed that its concentration is not high enough to affect viscosity or to invalidate Darcy's and Fick's laws (Hassanizadeh and Leijnse, 1988). The flow equation is taken to satisfy initial and boundary conditions. For the initial condition: where h o is the initial equivalent head. For the condition on the Dirichlet boundaries r 1 : where h is the prescribed head. For the condition on the Neumann boundaries f 2: where q0 is the prescribed fluid flux normal to the boundaries r 2 (positive inward), and U; is the Darcy flux given by: For the condition on the boundaries fa (combination of the Dirichlet and Neumann boundaries):

. 2·Solute Continuity Equation
The solute continuity equation is expressed as: at (Bear, 1979), and C* Ck for injection; C* Cs for abstraction, ( 1 0) where Ck and Cs are the relative concentrations of the solute injected and pumped with the fluid sources, respectively, and D i ; is the dispersion tensor which is defined as in Galeati et al. (1992).
212 TAO, Vol.6, No.2, June 1995 The solute transport equation is taken to satisfy the initial condition: where Co is the initial concentration. For the condition on the Dirichlet boundaries r 1 : where C is the prescribed concentration. For the condition on the Neumann boundaries r 2: where q d is the prescribed dispersive flux nortnal to the boundaries r 2 (positive inward). For the condition on the Cauchy or Robin boundaries r 3: where q c is the prescribed flux of concentration across the boundaries.

. GALERKIN FINITE ELEMENT SOLUTIONS
The numerical formulations used below to solve density-dependent flow and transport in a landfi ll site follow the Galerkin finite element method.

. 1 Fluid Flow Equation
Recalling (1) for the continuity equation of the fluid, a two-dimensional system of coor dinates ( x, z) which does not necessarily coincide with the principal axes of the conductivity tensor can be written as: where i, j, m, n are nodal points for each rectangular element, and N is the basis or "' weighting function. The residuals R ( h) are then defi ned in te1·1ns of trial heads as: . R ( h )N Ldxdz 0 L = 1, 2 , · · · , n n total number of nodes.

h
( 1 8 ) where NL is the basis function at L=l,2, ... , n nodal points, n x and n z are x and z components of unit vector outward and nor1nal to the boundaries, and c m is the mean . . concentration over an element. For the sake of simplicity, C is taken to represent the spatial average within each element (Galeati, et al., 1992 (2 3) The above equation can be used to numerically solve the head value for each nodal point in inhomogeneous, anisotropic media.

Point Velocity
To avoid discontinuity across the boundary between two elements, the flow velocity for each nodal point in inhomogeneous, anisotropic media is solved from Yeh (1981) by: .
where R x , R z are residuals at x, z components, and § x , li z are trial solutions of Darcy velocities at x, z directions. With the resulting residuals set orthogonal to all the basis functions, NL over the flow domain A yields: A q z can be achieved in a similar fashion. Substituting the trial solution h in (16) into (25) and subdividing the whole domain into N rectangular elements yields: Equation (26) can be written in matrix form as: where [ N ] is the basis function matrix, ( K ] is the hydraulic conductivity matrix, and [ K /] is the hydraulic conductivity (at x, z coordinates) matrix with the density contrast. Therefore, the flow velocity for each nodal point in inhomogeneous, anisotropic media can be found.

Solute Transport Equation
It should be recalled that the solute continuity equation is expressed as: a 2 c a 2 c a 2 c a 2 c <f> Dxx 8 x2 + </>Dx z 8 x 8 z + </>D z x 8 z 8 x + </>D zz 8 z2 ac ac U x 8 x + U z 8 z - where C* =C k for injection and C* =Cs for abstraction. C k and Cs are as defined before. It is to be noted that Dxx , Dx z , D z x and D zz are allowed to vary in space. U x and U z , the Darcy fluxes, are assumed to be constant in space since the density differences as used in this research were small, and 8� and 8tfzz =O. If the density differences were high, 8Jfx" and 8Jlx1 can not be ignored, and hence, U x and U z are not constant in space. The procedure for the Galerkin finite element scheme in (28) can be achieved in a similar fashion. The trial A solution C is first defined to approximate ·the solution C such that: A where C is the concentration. The residu . als R( C) are defined in terms of concentrations as: With the approximate solution (29) substituted into (30), and the resulting residuals set orthogonal to all basis functions NL over the flow domain A:· A A R( C)N Ldxdy 0 L = 1, 2 , · · · , n n total number of nodes.
ac aNL zz az az )dx d z The Galerkin finite element scheme for solute transport in inhomogeneous, anisotropic media is similar to that of the fluid flow equation and can be written in matrix form as:

RESULTS
Since the analytical equation for (1) by Galeati et al. (1992) has not yet been refined in this research isochlors published by Frind (1982) for leachate transport are compared with those of Galeati et al. (1992), recalling that the governing equation of Frind (1982) does not include the second term as shown in equation (1).
. Additionally, the mass balance error was used to evaluate the accuracy of the numerical solution. The mass balance error can be expressed by: [Net material flux across the whole boundary] + [net rate of mass production owing to well pumping or injection] -[rate of mass storage] = 0 ( 3 4 ) Equations (9) or (28) can be written as: where Q C* is the net rate of mass production owing to well pumping or injection, and C is the trial solution of the concentration. For the Cauchy boundaries: ....
For the Neumann boundaries: For the Dirichlet boundaries: (35) can be written as: where E is the rate of mass loss, nb is the number of nodes on the whole boundary, n q is the number of pumping and/or injection wells, and nd is the total number of nodes in the whole solution region. Thus, the mass balance error (E) used as an indicator of the accuracy of the numerical solution of the transport . equation is expressed as: where pis the time level, and k is the cumulative time level. The mass balance errors shown in the simulation results are computed from equation ( 41 ).  In the simulation of this case, 275 rectangular elements and 216 nodes were used; �t was allowed to increase gradually by a . factor of 1.2 from 0.000 1 to 100 days. Figure 1 shows that the modeling results of isochlors at the transport times of 1.5, 3, and 6 years computed by Frind (1982) (in dashed lines) and this research (in solid lines) are in fair agreement. The concentration distributions . at 1.5 years of Figure l(a) with a 3.4% mass balance error (E) agree with those by Frind (1982). However, the position of 0.1· isochlor at 3 years by Frind ' (1982) has migrated to the aquifer bottom, but not in Figure l(b) from this research with a 4.6% mass balance error. This implies that the solute transport at 3 years according to Frind (1982) is faster than that in this research. The concentration distributions by Frind (1982) at a transport time of 6 years, however, agree with those of Figure 1 ( c) in this research with a 6.2% mass balance error. This indicates that the solute transport after a long time in this study and by Frind (1982) are in good agreement. The governing equation for fluid flow by Galeati et al. ( 1992) as in equation ( 1) is mostly updated and more advanced than that used by Frind (1982). Therefore, equation (1) is used to model solute transport from a landfill in this research.

SIMULATION RESULTS
The validation of the model in this research on homogeneous, anisotropic and inho mogeneous (heterogeneous), anisotropic aquifers allows for its application to other cases of interest. The following scenarios are for modeling leachate solute transport in wedge-shaped aquifers in landfill disposal sites.    nodes were used; �t was allowed to increase gradually by a factor of 1.25 from 0.000 1 to 100 days. Figure 3 shows concentrations at the transport times of 1, 3 and 6 years without density-dependent ftow (/=0). The 0.13 concentration lines in (a,b,c) of Figure 3 have not yet migrated to the aquifer bottom. Now it should be assumed that the density contrast (/) between the contaminant and the groundwater is 0.004. Figure 4 shows the concentrations at transport times of 1, 3, and 6 years with an approximate mass balance error of 2.0%. Figures 4{a· ,b) and 3(a,b) show no changes in the concentration distribution at the transport times of 1 and 3 years. However, the concentration of 0.13; at the transport time of 6 year, as shown in Figure 4( c ), has reached the bottom of the wedge-shaped aquifer, but not in Figure 3(c). When the d_ ensity contrast increases to 0.007 (/=0.007), the concentration of 0.16 at the transport time of 6 years has nearly ·reached the bottom as shown in Figure  5(c). From the comparisons of the concentration distributions· shown above, it becomes evident that the concentration distributions after the short-time transport of one year show no significant difference. However, when the transport time increases to six years, the concentration distribution and migration pattern of the contaminant plume obviously change. It is especially true that when the contaminant solution is denser than the groundwater, the plume sinks downward into the groundwater flow system and disperses away from the contaminant source.

Case2: Inhomogeneous, Anisotropic Unconfined Aquifer
. This scenario concerns an inhomogeneous, anisotropic wedge-shaped unconfined aquifer with an aquitard lens. The fi nite element discretizations are shown in Figure 6. The boundary conditions and system parameters, except for the aquitard lens, are the same as those in Case 1, and the system parameters for the lens adopted for the numerical simulation are: Kxx=0.02 ml day, Kzz=0.02 m/day, The modeling results of solute transport at 1, 3, and 6 years are shown in Figure 7, and they show that the contaminant can migrate through the aquitard lens, but its concentration contours are distorted, especially at 3 and 6 years (Figures 7b, 7c). This indicates the transport migration pattern can also be influenced by the hydraulic conductivity contrast between different geologic materials. The hydraulic conductivity of the aquifer and aquitard materials are hundred to one, i.e., two orders of magnitude different. With differences of more than two orders of magnitude in hydraulic conductivity, the transport migration pattern can be affected.

Example Site of Landfill
The numerical model was used for application in a practical way to a certain landfi ll in Taiwan. Figure 8 shows the site along with the location of its boreholes. The B-B' geological section was made as shown in Figure 9, and the hydrogeologic conditions and boundaries for the numerical simulation are known. The whole profile of the site from bottom to top is grey silty clay, grey silty sand or sandy silt (aquifer 2), grey silty clay occasionally with fine sand inclusion, and brown silty fine sand (aquifer 1). The hydrogeologic condition of this site is a homogeneous and anisotropic unconfined aquifer (i.e., aquifer I), which is underlain by an homogeneous, isotropic aquitard and a homogeneous, anisotropic confined aquifer (aquifer ..c z. x.
. . The boundary conditions and fi nite element discretizations consisting of 89 rectangular ele ments and 26 triangular elements and 125 nodal points are shown· in Figure 10.
In simulation, the time interval D.t was allowed to increas·e gradually by a factor of 1.25 from 0.001 to 100 days. The equipotential lines and flow lines for a steady-state flow in the confined aquifer are shown in Figure 11 (a). The ground ·water flow direction is from left to · right, and thus, the leachate contaminant plume is moved to the right. The concentration distributions after a 1000-day (2.7 years) and a 2000-day (5 . 5 years) transport with the density contrast of 1 =0.0071 are shown in Figures 1 l(b) and (c), respectively. In . . comparison, the plume after a 2000-day transport as in Figure ll(c) spreads more in the horizontal direction than after a 1000-day transport as in Figure 11 (b ). The plume after a 2000 -day transport has percolated and reached the confined aquifer (i.e., aquifer 2), but not i1,1 the case of the 1000-day transport. This demonstrates that the contaminant . plume sinks downward into the groundwater flow system in aquifer 2 after a longer transport time of 2000 days if the density contrast between the contaminant fluid and the groundwater is more signifi cant (i.e., /=0.0071). This means that the migration pattern of the plume is changed due to the density effect.

DISCUSSION
In the model of density-dependent flow and transport, the change in contaminant con centration may induce the velocity change of fluid flow. If the contaminant solution is considerably denser than the groundwater, the vertical fluid flow velocity increases.. Be fore the solute transport reaches a steady-state, the contaminant concentration changes with time, causing the fluid flow velocity to change with time. Therefore, the solute transport should be considered in te1111 s of an unsteadystate fl ow fi eld. It is to be noted the fluid fl ow velocity (Ui) in the advection ter·m of equations (9) and (28) is not a constant, i.e., 8 (UiC)/ 8 xi Ui( 8 C/ 8 xi) + C( 8 Ui/ 8 xi)· However, in ge�eral, Ui is dealt with as a constant (i.e., 8 ( ui C) I ox i ui ac Ix i) when the density differences are small, and aui I ox i can be ignored. This indicates that the change in the contaminant solution does not induce a distinct difference in aui I ax i. With the consideration of c aui I ox i neglected, the computation time can be reduced without any signifi cant loss of accuracy. However, If the density differences are high (I >0.1 ), ui is not constant in space and, hence, aui I ax i cannot be ignored. Otherwise, the simulated results would be subject to significant errors.
Besides the effect of density contrast, the modeling results for Cases 1 and 2 in this research indicate that the transport migration pattern can also be influenced by the hydraulic conductivity (K) contrast between the aquifer and aquitard materials, especially for differ ences of more than two orders of magnitude in hydraulic conductivity. The porosities ( </>) are also different for the aquifer and the aquitard lens. Therefore, the values of the fluid fl�w velocity (Ui) in the advection ter1ns of the solute transport equation (i.e., equation (9)) are affected. Accordingly, the concentration distribution and migration pattern is changed. Contaminants travel fast in the aquifer and become slower through the aquitard lens, thereby causing the concentration contours to be distorted, especially at longer transport times of 3 to 6 years.

800
Modeling leachate transport ·fo r the example site: (a) steady-state flow . The modeling results for the example site of landfill in this research imply that the transport migration pattern can be affected by the longitudinal dispersivity (aL) and transverse dispersivity (ar). The values of aL and ar in the aquifers I and 2 are greater than those in the aquitard layer . . Contaminants can spread to occupy a portion of the fl ow system many times larger than would be the case in the presence of advection alone. If the transverse dispersivity is very large, contaminants transported along relatively horizontal fl ow paths can . migrate deep into the flow system. Freeze and Cherry (1979) indicated that if dispersivity values are orders of magnitude larger than the values obtained from laboratory experiments, the dispersion exerts a strong influence on the contaminant transport.
From long-terrn water-quality monitoring, in a real case of a landfill site, McFarlane et al. ( 1983) found that the chemical reaction heat from solid waste decomposition was associated with the leachates and was carried into aquifers. The aquifers then became non isother1nal · due to the reaction heat, with a maximum temperature difference of 5.5°C which . can affect the viscosity (5.5°C can reduce leachate viscosity by 10°). However, this research assumes that the contaminant tran sport is isothermal with no chemical heat f or1nation during the solid waste decomposition, and therefore, the velocity change in the fluid fl.ow is induced only by the density contrast between the contaminant solution and the groundwater. This velocity change can affect the concentration distribution and migration pattern of the plume in the contarninant transport process.

. CONCLUSIONS
Because the leachate from a landfill site into the groundwater can change the con taminant concentration, there exists a density contrast between the contaminant fl uid and the groundwater in the su�ounding flow domain. This density contrast can affect the groundwater flow and, consequently, the migration patterns and concentration distribution of contaminant plumes. The density gradients can infl uence . the fl.ow pattern and hence the tran sport of solutes infiltration from a landfill.
The simulation results in this research are for the problems of density-dependent ground water flow and transport in homogeneous, anisotropic and inhomogeneous, anisotropic aquifers of a landfill site by means of a Galerkin finite element method. It is shown that if there does not exist a significant density contrast, the contaminant plume spreads in a shallow zone close to the water table; if there is such a contrast, the plume will sinks downward into the groundwater flow system, and its migration pattern is significantly changed after a longer transport time of six years. However, the migration pattern is not significantly changed in a shorter time period of one to three years.
The modeling results by Galeati et al. (1992) are in good agreement with those by Frind (1982), as shown in Figure 1. The numerical model of Galeati et al. (1992) was used to simulate the leachate transport problems including the example site. All showed <5% mass balance errors· and it can be concluded that the modeling results can also be applicable to a landfill site.