QUT Digital Repository: http://eprints.qut.edu.au/ This is the author version published as: This is the accepted version of this article. To be published as : This is the author’s version published as: Liu, Fawang and Anh, Vo and Turner, Ian (2003) A twoDimensional finite volume method for variable density flow and solute transport through saturatedunsaturated media. In: The proceeding of the International Symposium on Nonlinear Science and Application, Shanghai, 2003. Catalogue from Homo Faber 2007 Copyright 2003 [please consult the authors] A Two-Dimensional Finite Volume Method for Variable Density Flow and Solute Transport Through Saturated-Unsaturated Media F．Liu a ,b , V. a Anh b and I. Turner b Department of Mathematics, Xiamen University, Xiamen 361006, China email@example.com b School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, Qld.4001, Australia firstname.lastname@example.org, email@example.com, firstname.lastname@example.org ABSTRACT Extensive groundwater withdrawal has resulted in a severe seawater intrusion problem in the Gooburrum aquifers at Bundaberg, Queensland, Australia. Better management strategies can be implemented by understanding the seawater intrusion processes in those aquifers. To study the seawater intrusion process in the region, a two-dimensional density-dependent, saturated and unsaturated flow and transport computational model is used. The model consists of a coupled system of two non-linear partial differential equations. The first equation describes the flow of a variable-density fluid, and the second equation describes the transport of dissolved salt. A two-dimensional control volume finite element model is developed for simulating the seawater intrusion into the heterogeneous aquifer system at Gooburrum. The simulation results provide a realistic mechanism by which to study the convoluted transport phenomena evolving in this complex heterogeneous coastal aquifer. 1. Introduction Seawater intrusion is an important environmental problem for coastal aquifers of Queensland, Australia, because they supply up to 370000 megalitres of groundwater per annum for both irrigation and water supply for town and industrial purposes. Excessive demand for groundwater in coastal areas has stressed the aquifers beyond their long term yield, resulting in seawater intrusion with a substantial loss of agricultural land. The general occurrence and distribution of groundwater in the Bundaberg area is rather well understood. However, the dynamics of the circulation of saltwater is not well understood. Simulation of variable density flow and transport through saturated-unsaturated media is a difficult problem. The modelling is further compounded by the complex regional aquifer geometry in the Bundaberg area. This area is unconfined and the water surface fluctuates with time. Continued development of groundwater resource for irrigation has raised concerns of how much additional extraction can be allowed from various borehole locations without adversely affecting the water quality as a result of saltwater invasion. In this paper, a two-dimensional control volume finite element model is developed for simulating the seawater intrusion into the heterogeneous aquifer system at Gooburrum. It is demonstrated that the transport simulation can provide an excellent tool for a hydraulic investigation even when complex transition zones are involved. 2. Statement of the problem Groundwater level contours of the Bundaberg area indicate that the major flow direction of water in the upper aquifer is from the southwest to the northeast with the discharge mainly to the ocean. Under such circumstances, preliminary assessments made from two-dimensional simulations can be reasonably informative. Here, a two-dimensional density-dependent, saturated and unsaturated flow and transport computational model was used to analyse the Bundaberg groundwater system. For this study, one vertical section is modelled (refer to section A-A in Figure 1). N Pacific Ocean A Moore Park Avondale Burnett Heads describes the transport of dissolved salt. The flow equation may be written as [see Voss, 1984] Escarpment A Burnett Ri ve r Sw P C Sw Sop P t Sw C t Bundaberg Figure 1. Plan view of Gooburrum area. （1） Kk r P g QP where P ( x, z , t ) is the fluid pressure; C ( x, z , t ) is the solute concentration; QP ( x , z , t ) is the fluid mass source; is the aquifer volumetric porosity ; is the fluid dynamic viscosity; g is the gravity acceleration vector; is S op is the specific pressure storativity; ( x, z , t ) the density of fluid; Sw is the water saturation; kr is the relative permeability to fluid flow; K is the solid matrix permeability. Zone No Figure 2. Cross section along AA shown in Figure 1 with different zone permeabilities. Figure 2 shows the geological units in the vertical section interpreted using the available data from the observation bores in the area. This is a heterogeneous and anisotropic aquifer system and section A-A comprises two aquifers that are separated by a leaky layer. The zones in the diagram numbered from 1 to 8 indicate the regions of different permeabilities. The upper aquifer consists primarily of sand and gravel with minor clay, while the top surface is the rainfall recharge. Rainfall recharge infiltrates through the unsaturated portion of the upper aquifer and then reaches the water table. Thus, the upper aquifer is treated as unconfined. The lower aquifer (zones 7 and 8) consists of sand, gravel and clays. The leaky aquitard (zore 6) consists mainly of silty or sandy clays and clays. The permeability of the aquitard is estimated to be approximately 5 orders of magnitude less than that of zone 8. The inland boundary and the bottom boundary of the domain are formed by relatively impermeable sediment. The problem of seawater intrusion into coastal aquifers can be formulated in terms of two partial differential equations. The first equation describes the flow of a variable-density fluid, and the second equation Permeability m 1 1.20 1011 2 4.00 1011 3 1.00 1010 4 4.00 1010 5 5.00 1010 6 2.00 1015 7 3.50 1011 8 3.50 1010 2 Table 1. Permeability values used for the simulation Table 1 shows the permeability values corresponding to the zones numbered from 1 to 8 depicted in Figure 2. The permeability values have been assumed such that the flow in the horizontal direction is three orders of magnitude more conducting than that in the vertical direction. To describe solute transport, the following form of the Fokker-Planck equation is used [see Voss, 1984]: S wC S w D C S wCv QP C * t (2) where C * is the solute concentration of fluid sources; v ( x, z, t ) is the fluid velocity; D is the dispersion tensor. A two-dimensional xz coordinate system for the tensor is given as [Bear, 1979]: v2 v2 v2 v2 Dxx L x T z Dd Txx , Dzz T x L z Dd Tzz , v v v v vv Dxz Dzx L T x z v (3) where Dd is the molecular diffusivity coefficient; Txx the principal components of the and Tzz are tortuosity tensor; L and T are the longitudinal and transverse dispersivities respectively. The Darcy velocity vector may be expressed as Kk r v P g . Sw (4) To obtain a unique solution to (1) and (2), appropriate initial and boundary conditions must be specified. The boundary conditions used for simulating the flow in the aquifer are also shown in Figure 2. Note that the three layer system has been treated as an integral part of the dynamic aquifer system and is modelled as a single system in a cross-sectional simulation. A no-flow condition is imposed along the inland and basal boundaries of the cross-section. The upper boundary is a recharge boundary and the recharge enters variably along the top of the aquifer as exhibited in Figure 2. 3. control volume finite element method In this paper, a control volume finite element method (CVFEM) is developed for simulating seawater intrusion into the Gooburrum aquifer systems. A two-dimensional finite volume finite element method based on a triangular background interpolation mesh using linear shape functions has been proposed previously for analysing the evolution of the seawater intrusion into single and multiple coastal aquifer systems [Liu et. al, 2002]. The model formulation consisted of a ground-water flow equation (the reference hydraulic head as a variable) and a salt transport equation (the solute concentration as a variable). A two-dimensional finite volume unstructured mesh method based on a quadrilateral background interpolation mesh using bilinear shape functions has also been proposed for analysing the model system [Liu et. al, 2003]. The model formulation consists of a flow equation (the reference hydraulic head as a variable) and a modified Fokker-Planck equation (the solute concentration as a variable) [see Frind, 1982; Huyakorn et. al, 1987]. In this paper, a two-dimensional control volume finite element method based on a quadrilateral background interpolation mesh is developed for simulating seawater intrusion into aquifers with both saturated and unsaturated regions. This model formulation consists of a flow equation where the fluid pressure is used as a variable and a Fokker-Planck equation for the solute concentration. The flow and solute transport across each control surface must be determined by an integral. The CVFEM discretisation process is initiated by utilising the integrated form of Eqs. (1) and (2). Integrating the flow Eq. (1) and the transport Eq. (2) over an arbitrary control volume yields: S w Sop v S w P C dv S w dv P t C t v (5) Kk r P g QP dv v v S wC dv t (6) S w D C dv S wCv dv QP C dv * v v v Applying the Gauss divergence theorem and the lumped mass approach to Eqs. (5) and (6) gives S w dPp dC p S w Sop P dt v p S w C dt v p p p (7) Kk r P g dn QP p v p S Sw p dC p dt vp (8) Sw D C dn S wCv dn QPC v p p * S S where dn represents the components of the outward surface vector normal to the control surface S, and an anticlockwise traversal of the finite volume integration is assumed, i.e., dn can be approximated in the discrete S w dPp dC p S w S op P dt v p S w C dt v p p p sense by dn zk xi ; x and z represent the x and z components of the sub-control volumes (SCV) face; v p is the area of the control volume, and is i 1 where N PSCV is the total number of the sub-control volumes that make up the control volume associated with the node p . The integrals in Eqs. (7) and (8) are line integrals, which will be represented by the midpoint approximation for each control surface in order to attain second order accuracy. To approximate this midpoint, the argument of the integrals is required at the midpoint of the control surface and it is for these surfaces that the outward normal vector will be required. The integral in Eq. (7) can be rewritten as Kk r P g dn S j N PSCV Kk r P j j g dn j j 1 S j 1 S j 2 j N PSCV 2 K k P l , j x r zlj x j 1 l 1 l ,ups j K k P l , j z r xlj l ,ups z l 1 4 j K x kr N il , j zl j l ,ups x i 1 (12) Sw p dC p dt vp S N PSCV j 1 2 l 1 4 i 1 w l, j j N il , j j j N i D D zl , , xx l xz l p x z (13) N N il , j j j Dzxj ,l Dzzj ,l xl Ci x z l, j i N PSCV 2 4 K k j N l , j i Cl j,ups x r zl j x j 1 l 1 i 1 , l ups j K k N il , j z r xlj Pi j l ,ups z j K k j Cl .S g xlj QP C * v p z r 0 p C l ,ups The concentration (10) Cl j,S in the gravitational term within each sub-control volume is represented by the use of shape functions at the midpoint ( IPl ( xIP , zIP ), l 1,2 ) of each l l SCV, i.e., the Clj,S will be evaluated at the SCV integration point for each face: K k z r l j g xlj l ,ups 4 Cl j, S Ci j Nil , j xIPl , zIPl . j i 1 where the upstream weighting technique will be used in this work: u l j, u p s u j 1 j K k j Cl .S g xl j QP p v p z r 0 C l ,ups (9) SCVi 2 j N PSCV v [ K k N il , j z r xlj ]Pi j l ,ups z evaluated for the vertex case as vp N PSCV (14) The components of the fluid velocity will be approximated by v l vx2,l vz2,l , j p , if u l j, u p s u l j , i f v v n n C S C S 0, l 0 (11) l where CSl , (l 1, 2) is a control surface. To evaluate the terms in Eqs. (10) and (8), a local-global co-ordinate transformation is introduced. It is convenient to work in a local co-ordinate system so that each element may be treated identically, irrespective of how distorted an element may be in terms of the global co-ordinates. Using CVFEM, the final form of the discretised equations become (15) j Kk Nil , j j v x ,l x r Pi , x i 1 S w l ,ups 4 j j 4 Kk Nil , j j K z kr j v z ,l z r Pi Cl .S g . 0 x C i 1 S w l ,ups Sw l ,ups Equations (1) and (2) are transformed into a system of ordinary differential equations (ODE) of the form (7) and (8) for each node using CVFEM. In this work, we used DASSL (the differential algebraic system package) as the ODE solver [Pelzold, 1982]. This technique has been used previously with good success to solve adsorption problems involving steep gradients in a bidisperse particle [Liu and Bhatia, 1999 and 2001b], hyperbolic models of transport in bidisperse solids (Liu, Bhatia and Abarzhi, 2000), transport problems involving steep concentration gradients [Liu and Bhatia, 2001a] and modelling saltwater intrusion into coastal aquifers [Liu, Turner and Anh, 2002; Liu, Turner , Anh and Su, 2003]. 4. Numerical results and discussion The numerical techniques described in this paper have been implemented in Fortran-77. The CVFEM is used to simulate seawater intrusion into the aquifer system at Gooburrum. Figure 3 shows the non-uniform finite element mesh (724 rectangular elements and 817 nodes) employed for the chosen cross-sectional region (see Figure 1). In this application, a transient situation where pumping from both the upper and lower aquifers was activated. located at a distance of 11.85 km from the inland end of the domain at depths of 48.5, 51.5, 54.5 and 57.5 m below AHD (see Figure 2). In order to simulate the seawater intrusion into the aquifer system when the pumping from both upper and lower aquifers is activated, the steady state concentration and pressure distribution were used as the initial values. Figures 4 and 5 show the salt concentration distribution after 5 years and 10 years of pumping for a total extraction of 347 ML per annum, respectively. Here, 126 ML of water per annum is pumped from the bottom aquifer and 221 ML of water per annum is pumped from the upper aquifer. The simulation results are again similar to those observed by Bajracharya et al  using SUTRA. 60 50 40 ELEVATION, m 30 20 10 0 -10 -20 -30 -40 -50 -60 -70 0 5000 10000 15000 D IS TAN C E , m Figure 3. Nouniform finite element mesh for cross-sectional model of Gooburrum, Bundaberg, Queensland, and Australia. Figure 4. Concentration distribution after 5 years of continuous pumping at a constant rate of 347 ML per annum. Concentration units are in mass fraction. The initial pressure and concentration were also set to zero. Spatially varying recharge was used along the top boundary. The ground surface is assigned a time-independent spatially variable rainfall recharge. The recharge intensities are 60 mm/year for about 4.5 km and 90 mm/year for about 9.6 km. Note that the seawater concentrations were assigned at the specified pressure boundary condition at the right end of the problem domain. The seawater level in the offshore aquifer was fixed by specifying a constant pressure, i.e., P s g depth . In order to simulate the pumping from both the upper and lower aquifers, two wells were located in the top aquifer at a distance of 12.75 km from the inland end of the problem domain at depths of 6.5 and 9.5 m below AHD. Similarly, four wells in the bottom aquifer were Figure 5. Concentration distribution after 10 years of continuous pumping at a constant rate of 347 ML per annum. Concentration units are in mass fraction. 5. Conclusion Analysis of a freshwater-saltwater aquifer system that contains a relatively narrow transition zone is a very difficult problem. In this paper, a two-dimensional control volume finite element method for the density-dependent, saturated and unsaturated flow and transport model has been described. This model has been used to simulate seawater intrusion into the aquifer system at Gooburrum. It is demonstrated that the transport simulation can provide an excellent tool for hydraulic investigation even when complex transition zones are involved. Acknowledgements This research has been supported by the ARC SPIRT grant C10024101, the QUT small grant scheme 190752-0003/08 and the National Natural Science Foundation of China 10271098. References  K. Bajracharya, J. Arunakumaren and W. J. Huxley, “Numerical modelling of seawater intrusion in Gooburrum, Bundaberg, queensland”, Paper presented at International Groundwater Conference, International Association of Hydrogeologists, University of Melbourne, Melbourne, Australia, 1998.  J. Bear, “Hydraulics of Groundwater,” McGraw-Hill, New Your, 1979.  E. O. Frind, “Simulation of long-term trasiet density-dependent transport in groundwater,” Adv. Water Resour, 5, pp73-88, 1982.  P. S. Huyakorn, P. F. Andersen, J. W. Mercer and H. O. White Jr., “Saltwater intrusion in aquifers: Development and testing of a three-dimensional finite element model,” Water Resour. Res., 23, pp293-312, 1987.  F. Liu and S. K. Bhatia, Computationally efficient solution techniques for adsorption problems involving steep gradients in bidisperse particles,” Comp. Chem. Eng., 23, pp933-943, 1999.  F. Liu, S. K. Bhatia and I. I. Abarzhi, “Numerical solution of hyperbolic models of transport in bidisperse solids,” Comp. Chem. Eng., 24, pp1981-1995, 2000.  F. Liu and S. K. Bhatia, ”Solution techniques for transport problems involving steep concentration gradients: application to noncatalytic fluid solid Reactions,” Comp. Chem. Eng., 25, pp1159-1168, 2001a.  F. Liu and S. K. Bhatia, “Application of Petrov-Galerkin methods to transient boundary value problems in chemical engineering: adsorption with steep gradients in bidisperse solids,” Chem. Eng. Sci., 56, pp3727-3735, 2001b.  F. Liu, I. Turner and V. Anh, “An unstructured mesh finite volume method for modelling saltwater intrusion into coastal aquifers,” Korean J. Comput. Appl. Math., pp391-407, 2002.  F. Liu, I. Turner, V. Anh and Su, “A two-dimensional finite volume method for transient simulation of time- and scale-dependent transport in heterogeneous aquifer systems,” J.. Appl. Math. Comput., 11, pp215-241, 2003.  G. F. Pelzold, “A description of DASSL: a differential/algebraic equation system solver, Technical Report No. SAND 82-8637, Sandia National Laboratories, Livermore,” 1982.  C.I. Voss, “USGS SUTRA Code – History, Practical Use and Applications in Hawaii, Chapter 9 in: J. Bear et. al. (eds.) Seawater Intrusion in Coastal Aquifers – Concepts, Methods and Practices, Kluwer, Dordrech, pp249-313, 2000.
Link or Click Back
Here will be a configuration form