An efficient model for solving density driven groundwater flow problems based on the network simulation method

Journal of Hydrology (2007) 339, 39– 53

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/jhydrol

An efficient model for solving density driven groundwater flow problems based on the network simulation method ´lez Ferna ´ndez A. Soto Meca *, F. Alhama, C.F. Gonza Applied Physics Department, Polytechnic University of Cartagena, ETSII Campus Muralla, 30203 Cartagena, Spain Received 17 July 2006; received in revised form 27 February 2007; accepted 7 March 2007

KEYWORDS Saltwater intrusion; Density-dependent flow; Groundwatermodelling; Network simulation method

Summary The Henry and Elder problems are once more numerically studied using an efficient model based on the Network Simulation Method, which takes advantage of the powerful algorithms implemented in modern circuit simulation software. The network model of the volume element, which is directly deduced from the finite-difference differential equations of the spatially discretized governing equations, under the streamfunction formulation, is electrically connected to adjacent networks to conform the whole model of the medium to which the boundary conditions are added using adequate electrical devices. Coupling between equations is directly implemented in the model. Very few, simple rules are needed to design the model, which is run in a circuit simulation code to obtain the results with no added mathematical manipulations. Different versions of the Henry problem, as well as the Elder problem, are simulated and the solutions are successfully compared with the analytical and numerical solutions of other authors or codes. A grid convergence study for the Henry problem was also carried out to determine the grid size with negligible numerical dispersion, while a similar study was carried out with the Elder problem in order to compare the patterns of the solution with those of other authors. Computing times are relatively small for this kind of problem. ª 2007 Elsevier B.V. All rights reserved.

Introduction The problem of saline intrusion in aquifers, of great interest nowadays, has been studied by a large numbers of research-

* Corresponding author. Tel.: +34 968325512; fax: +34 968325337. E-mail address: [email protected] (A. Soto Meca).

ers. Besides the books of Segol (1994), Holzbecher (1998) and Bear et al. (1999), many fundamental papers treat the problem, providing approximate analytical or numerical solutions. In the field of density driven groundwater flow, the Henry and the Elder problems (e.g., Henry, 1964; Elder, 1967b; Voss and Souza, 1987) have become two of the standard problems used as a benchmark for variable-density flow and transport tests.

0022-1694/$ - see front matter ª 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2007.03.003

40

A. Soto Meca et al.

Nomenclature a b c C d D Df g G j l M, N n Pe p q Q R Ra t v V x, y

discharge parameter (dimensionless) inverse seepage Peclet number (dimensionless) salt concentration (kg m3) capacitor (F) height of the aquifer (m) salt dispersion coefficient (m2 s1) molecular diffusivity (m2 s1) gravity acceleration controlled current source generalized flux (salt density flux or water density flux) also, electric current density (A/m2) length of the aquifer (m) number of the volume elements along the y and x coordinate, respectively porosity (dimensionless) Peclet number, Pe = Q/nD (dimensionless) pressure (Pa) specific discharge (m s1) recharge per unit with (m2 s1) resistor (X) Rayleigh number, Ra = kDqgd/Dfl time (s) velocity (m s1) electric potential (V) spatial coordinates (m)

With regards to the Henry problem, it was Henry who first developed a semianalytical solution for the steadystate salt distribution in a cross-section of a confined coastal aquifer, preserving the miscible displacement. The first authors to numerically reproduce the semianalytical solution proposed by Henry were Pinder and Cooper (1970), who obtained a transient solution by the method of characteristics. However, unlike Henry, they used a different form of the solute mass balance equation since they divided the specific discharge by the porosity. Four years later, Lee and Cheng (1974) solved the steady Henry problem by the finite element method, using a streamfunction formulation. Lee and Cheng compared their isochlor (c = 0.5) steadystate solution with those reported by Henry and by Pinder and Cooper. The three solutions are in relatively good agreement along the aquifer, but present certain differences at the toe of saltwater circulation. Another contribution to the unsteady solution of the Henry problem was provided by Segol et al. (1975) who, in addition to assuming Pinder’s mass balance equation, introduced a mixed boundary condition on the top portion of the seaward side to overcome the numerical difficulty caused by the physically unrealistic Henry boundary condition. They solved the problem by a Galerkin-finite element method, in which the pressure and velocities are obtained simultaneously in order to guarantee continuity in the velocity field. A numerical model for a low cost simulation was introduced by Frind (1982) for long-term transient response using the equivalent freshwater head formulation and the Galerkin-finite element method. Frind uses Pinder’s mass balance

j l q e n w Dx Dy

permeability of the porous medium dynamic viscosity (kg m1 s1) density (kg m3) density-dependent coefficient (m3/kg) aspect ratio (dimensionless) streamfunction (m2 s1) half length of the volume element (m) half thickness of the volume element (m)

Subscripts 1, 2 used to identify the controlled current source, G c related to salt concentration i, j centre of the dimensionless volume element i ± Dx right and left ends of the volume element j ± Dy upper and bottom ends of the volume element in input of the volume element 0 related to freshwater out output of the volume element s related to saltwater x, y related to spatial coordinates w related to streamfunction Superscripts 0 related to dimensionless quantity

equation and Segol’s boundary conditions, and studies the apparent differences obtained when using a constant or a velocity-dependent dispersion coefficient. A three-dimensional finite element model based on this formulation (with Pinder’s mass balance equation and Segol’s boundary conditions) was presented by Huyakorn et al. (1987), who studied several types of aquifer, including Henry’s. Huyarkon compared his isochlor (c 0 = 0.5) solution with those of Henry, Lee et al. and Frind, finding the solutions to be very close. Voss and Souza (1987) claimed that the discrepancies in the published papers were due to the use of two different dispersion coefficients. These discrepancies arose from Pinder and Cooper’s study (1970), since they did not use the effective recharge, as did Henry (1964), in the mass balance equation but this quantity was divided by the porosity. Senger and Fogg (1990) illustrated the use of streamfunction for modelling the Henry problem, while Evans and Raffensperger (1992) demonstrated that these functions must be defined in terms of mass flux (rather than volume flux) in order to diminish the significant errors that appear with their use, especially when flow parallels large density gradients. Galeati et al. (1992) investigated the extent to which the dependence of fluid density on salt concentration affects the numerical solution by solving the flow and advection– dispersion equations in coupled, partially coupled and completely decoupled modes. Croucher and O’Sullivan (1995) presented a grid converge study to reduce the numerical dispersion and used the streamfunction formulation, while Kolditz et al. (1998) commented that the Boussinesq approximation requires

An efficient model for solving density driven groundwater flow problems based on the network simulation method perpendicular Darcy velocity vectors and density gradients. Both authors also present a comparison of their solutions (isochlor, c 0 = 0.5) with previous results. More recently, Simpson and Clement (2003), who compare their results with those of Frind, provided a coupled versus uncoupled analysis, by means of a Galerkin-finite element model, to show that the saline water distribution in the Henry problem is mainly dictated by boundary-forcing and not by density-dependent effects. Gotovac et al. (2003) used the streamfunction formulation and a numerical collocation method (successively dividing the domain by half in each direction toward the discharge zone for the critical parameter b = 0.035), to study the different versions of the Henry problem. Their inland distance solutions, 1.37 and 1.16, for b = 0.1 and 0.035, respectively, coincides with the estimations of Croucher and O’Sullivan. In 1994, Segol provided a more precise semianalytical solution of the original Henry problem, enlarging the truncation to a scheme with 183-terms. Segol compared this solution with a numerical solution of the modified Henry problem (with Segol boundary condition or mixed condition). As regards this comparison, Simpson and Clement (2004) emphasized that, while this solution is the first to show agreement with other numerical solutions, the comparison is not appropriate. The authors re-evaluated both the semianalytical solution (assuming a scheme with 203terms) and the numerical solution for the standard Henry problem, and found an excellent degree of agreement between them. In addition, as an alternative option for increasing the influence of the density effects and to improve the worthiness of the Henry problem, Simpson and Clement obtained a new semianalytical solution for half an effective freshwater recharge. The Elder problem (Elder, 1967b), which treats heat-driven convection in a small domain, was physically solved by Elder using a Hele-Shaw cell to verify his numerical analysis based on the streamfunction formulation for a two-dimensional cross-section. This serves as an example of free convection phenomena, whereby the fluid flow within the domain is initiated by vertical temperature gradients which, in turn, induce a complex pattern of fingering to mix

through the domain. One justification for the use of Elder’s problem as a benchmark is the existence of laboratory results for comparison. Indeed, the Elder experiment was selected as a test case for model validation in the international Hydrocoin project (HYDROCOIN, 1990). Diersch (1981), and later Voss and Souza (1987), were the first to use Elder’s results in the context of solute transport in a porous medium. Since then, the solute analogy of the Elder problem has been widely adapted for use by the hydrology community. The Elder problem is representative of a number of hydrogeological systems where unstable variable-density flow behaviour is likely to be important, e.g., salt lake systems, saline disposal basins, dense contaminant and leachate plumes, geothermal reservoirs, etc. Simpson and Clement (2003) stated that the results of the Elder problem can only be matched in a qualitative sense because the problem is highly sensitive to discretization. They used a coupled versus uncoupled strategy to show that the Elder problem was a better test of variable- density flow phenomena than the Henry problem. In this paper, we present a powerful numerical model, based on the network simulation method (Gonza ´lez-Ferna ´ndez and Alhama, 2002), which has been demonstrated as accurate and efficient in other engineering fields (Zueco et al., 2004; Horno et al., 1995), to solve the 2-D salt intrusion unsteady problems of Henry, (in its different versions) and Elder, using the streamfunction formulation. The model is directly derived from the finite-difference differential equations based on spatial discretization of the governing equations. Once the model is designed, the user does not need to realize mathematical manipulations since the software used to simulate the model does this. Furthermore, as the conservation laws are inherent in the software (Kirchhoff’s current law), the numerical treatment related to these laws (both for salt and water density fluxes) does not need to be taken into account. This is an important advantage of the model, which also needs relatively small computing times. Solutions provided by the proposed method are compared successfully with those of other authors.

1.0

0.8

y'

0.6

0.4

0.2

0.0 0.0

0.5

1.0

x'

Figure 1

41

Geometry of the Henry problem.

1.5

2.0

42

A. Soto Meca et al.

The governing equations of the Henry and Elder problems Henry’s problem concerns a uniform, isotropic, rectangular aquifer, with initial uniformly zero salt concentration and impermeable top and bottom domain boundaries. The aquifer receives a constant freshwater recharge from the left side, while the right side is exposed to a stationary body of seawater. As time progresses, saltwater enters the aquifer at the base of the right boundary (Fig. 1) and mixes with freshwater (dispersion process), inducing a circulation from the floor of the sea into the dispersion zone from where it flows back into the sea. After a lengthy time, a final dynamic equilibrium state is reached. Formulation of the Henry problem (Croucher and O’Sullivan, 1995) involves Darcy’s law, the mass conservation (or continuity) equation and the (salt) transport equation. These are, respectively: j ~ gÞ ð1Þ q ¼  ðrp  q~ l rðq~ qÞ ¼ 0 ð2Þ oc þ r  ðc~ v Þ  r  ðDrcÞ ¼ 0 ð3Þ ot In addition, the concentration dependence of the saltwater density is assumed as q ¼ q0 ð1 þ ecÞ

ð4Þ

Since the dimensionless parameter e = (q  q0)/q is of the order of magnitude 0.025, the Boussisneq approximation hypothesis (q ffi q0) is assumed, and Eq. (2) reduces to r~ q ¼ 0. By introducing the new dimensionless variables 0

x ¼ x=l;

0

y ¼ y=d;

0

0

c ¼ c=cs ; 0

q ¼ ðq  q0 Þ=ðqs  q0 Þ ¼ c ;

0

q ¼ qðd=Q Þ;

t ¼ ðD=d 2 Þt 0

ow0 oy 0

q0y ¼

c0 ð0; y 0 Þ ¼ 0

ow0 ox 0

oc0 ðx 0 ; 0Þ oc0 ðx 0 ; 1Þ ¼ ¼0 oy 0 oy 0 w0 ðx0 ; 0Þ ¼ 0 w0 ðx 0 ; 1Þ ¼ 1 ow0 ð0; y 0 Þ ow0 ð2; y 0 Þ ¼ ¼0 ox0 ox 0

Eqs. (2) and (3) may be written in the dimensionless form (Lee and Cheng, 1974; Croucher and O’Sullivan, 1995; Holzbecher, 1998 and Gotovac et al., 2003), as ð5Þ ð6Þ

Figure 2

ð7aÞ ð7bÞ ð7cÞ ð7dÞ

while the initial condition is given by c0 ðx 0 ; y 0 ; t0 ¼ 0Þ ¼ 0

ð7eÞ

Elder’s problem concerns a uniform, isotropic, rectangular aquifer. No fluid flows in or out of any of the boundaries, except in the two upper corners of the domain held at zero pressure head. The pressure is initially hydrostatic and no solute is present. The solute source located along the middle half of the top boundary is specified as a constant unit value, while the bottom boundary is maintained at a constant zero value. For the streamfunction, a Dirichlet condition is imposed on all the edges. As regards solute transport, the vertical edges obey the Neumann condition. The normalized concentration has a Dirichlet condition c 0 = 1 for the solute source and c 0 = 0 at the bottom boundary. At the upper boundary with no solute source a Neumann condition prevails. The geometry of the Elder problem is given in Fig. 2. The governing equations, in dimensionless form, using the streamfunction formulation are o2 w0 o2 w0 oc0 þ ¼ Ra ox 02 oy 02 ox 0  0 0  2 0 2 0 oc oc ow oc ow0 oc0 oc0 ¼ 0 þ   02 02 0 0 0 0 ox oy oy ox ox oy ot

x 0 ¼ x=d;

o2 w0 o2 w0 1 oc0 þ ¼ ox 02 oy 02 a ox 0   o2 c0 o2 c0 1 ow0 oc0 ow0 oc0 oc0 ¼ 0 þ   02 02 0 0 0 0 ox oy ox oy ot b oy ox

c0 ð2; y 0 Þ ¼ 1

ð8Þ ð9Þ

with Ra the Rayleigh number (Holzbecher, 1998). The dimensionless variables are:

and the streamfunction (that comes from r~ q ¼ 0) q0x ¼

where a ¼ jgdðqlQS q0 Þ and b ¼ nD are the discharge parameter Q and the inverse of the seepage Peclet number (b = 1/Pe), respectively (Lee and Cheng, 1974). The boundary conditions expressed as a function of the dimensionless variables, for an aspect ratio n = l/d = 2, are

y 0 ¼ y=d;

q0y ¼ qy ðd=Df Þ; w0 ¼ w=Df ;

q0x ¼ qx ðd=Df Þ;

t0 ¼ tðDf =d 2 Þ;

c0 ¼ ðc  co Þ=ðcs  co Þ;

The boundary and initial conditions are given by the equations: c0 ðx 0 ; 0Þ ¼ 0 c0 ð1 6 x 0 6 3; 1Þ ¼ 1 oc0 ð0; y 0 Þ oc0 ð4; y 0 Þ ¼ ¼0 oy 0 oy 0

Geometry of the Elder problem.

ð10Þ ð11Þ

An efficient model for solving density driven groundwater flow problems based on the network simulation method w0 ðx 0 ; 0Þ ¼ 0 w0 ðx 0 ; 1Þ ¼ 0 w0 ð0; y 0 Þ ¼ 0 w0 ð4; y 0 Þ ¼ 0 ð12Þ oc0 ðx 0 ; 1Þ ¼ 0 0 6 x 0 < 1 and 3 < x 0 6 1 ð13Þ oy 0 c0 ðx 0 ; y 0 ; t ¼ 0Þ ¼ 0 ð14Þ

The network model The network simulation method (Gonza ´lez-Ferna ´ndez and Alhama, 2002) is based on the electric analogy of the transport processes, in which the following equivalence between variables is generally defined: j (electric current, W/m2)  flow variable of the problem V (electric potential, V)  potential variable of the problem Classical thermoelectric analogy has nothing to do with the network simulation method because it always refers to linear problems and uncoupled equations. If the mathematical model is defined by coupled equations, such as occurs in this work, separate networks are implemented in the model, one for each pair of conjugate variables. The coupling of the equations is easily implemented in the network model by using a special type of device that is integrated in the libraries of the simulation codes, called a dependent current source. This source is able to provide a current that is a function of both the flow and/or potential variables at any node of the medium. The linear terms of the equation are also generally implemented by simple capacitors (storage terms) and resistors (dissipative terms). Since only three elements make up the model (resistors, capacitors and sources), very few programming rules are required for designing the model. Each term of the equations is, in this way, implemented by a sole electric element, which is interconnected with others according to the algebraic sign in the equation. Finally, each network of the volume element is directly connected to adjoining networks to make up the whole

Figure 3

43

medium and, finally, boundary conditions are also implemented by suitable elements. Once the complete network model has been designed, its simulation is carried out using a standard circuit analysis code such as Pspice (PSPICE 6.0, 1994), without any other mathematical requirements. This provides the nearly exact solution of the model. Pspice was initially developed by the Integrated Circuit Group of the Electronic Research Laboratory at the University of California, Berkeley, and the first version was finished in 1992. Since that time, thousands of copies of the new versions of the program have been given to universities and companies in the electronic industry. The widespread use of Pspice and its new versions attests to the applicability of the program to a large variety of circuit simulation problems and provides a valuable base of experience that demonstrates the advantages and disadvantages of the powerful, efficient and reliable numerical algorithms that are employed in the program. The system of equations that describes the complete circuit is determined by the model equations for each element and topological constraint, which are determined by interconnecting the elements, reflecting Kirchhoff’s Current and Voltage Laws (Kisrchhoff’s Current law is a conservation principle for the flow and solute, while Kirchhoff’s Voltage law ensures the uniqueness of the concentration and streamfunction variables). For transient analysis, Pspice maintains an internal timestep, which is continuously adjusted to maintain accuracy, while not performing unnecessary steps. The time-step is reduced by the code during the simulation so that the integrated charges and currents are sufficiently accurate. During periods of inactivity, the internal time-step is increased, while during active periods, it is decreased. The maximum internal step size can be controlled by using software to specify it. The minimum time-step is the overall run time divided by 1E12, while the initial time-step is fixed as a function of both the total time required for the simulation and the value of the relative accuracy, a data (parameter) also specified by the user. If convergence is not

Nomenclature of the volume element.

44

A. Soto Meca et al. Table 1 problem

Properties and parameters used for the Henry

Symbol

6.6 · 10 (Henry, 1964) 2.31· 106 (Pinder and Cooper, 1970) 9.81 1.020408 · 109 6.6 · 105 0.35 1.0 · 103 1000 1025 2 1

g k q n l q0 qs l d

reached, the initial time-step is successively reduced in the simulation until convergence disappears. The circuit equations are a system of algebraic-differential equations of the form F(x, x 0 , t) = 0, where x is the vector of the unknown circuit variables, x 0 is the time derivative of x and F is, in general, a nonlinear operator. The equation formulation algorithm based on a combination of Cutset and Loop Analysis is a modified version of the clas-

Table 2

m2 s1

m s2 m2 m s1 – kg m1 s1 kg m3 kg m3 m m

sical Nodal Analysis. This method provides the same generality as other formulation methods, but produces a nearsymmetric system of equations that are solved with an amount of computational effort that is comparable to that needed for the simplest Nodal Methods. The Markowitz algorithm is used for reordering the system of equations. Once the equations F(x, x 0 , t) = 0 are reordered, direct elimination methods, such as LU factorization and sparse-matrix, are used to obtain the solution. If the circuit contains elements that are modelled by nonlinear equations, then the solution is obtained by an iterative sequence of linearized solutions. The modified Newton–Raphson algorithm, which approximates each nonlinearity by a Taylor series that is truncated after the first order term, is the most common method of linearization. The modified Newton–Raphson method contains the reliable and efficient simple-limiting of Colon (algorithm) that requires the lowest number of iterations for convergence. For the implicit transient analysis method of numerical integration, stiffly-stable algorithms are used, such as Trapezoidal integration. The local truncation error of the integration is proportional to the time-step, which is successively reduced to make the error negligible. Design of the network model starts from the finite-difference differential equation derived from the spatial discretization of the governing Eqs. (5) and (6), retaining time as a continuous variable. According to the nomenclature of the square volume element of Fig. 3, the resulting discretized equations are

Minimum and maximum time-step for different grid sizes b = 0.1

20 · 10 40 · 20 80 · 40

Unit 6

D

Figure 4 Network model of the volume element: (a) streamfunction variable and (b) concentration.

Value

b = 0.035

b = 0.2

Time-step minimum (s)

Time-step maximum (s)

Time-step minimum (s)

Time-step maximum (s)

Time-step minimum (s)

Time-step maximum (s)

8.35E6 6.55E6 6.55E6

0.02 0.02 0.02

8.35E6 3.80E6 3.80E6

0.02 0.02 0.02

5.18E6 3.80E6 3.80E6

0.02 0.02 0.02

An efficient model for solving density driven groundwater flow problems based on the network simulation method  0   0  D Dw Dw  Dx 0 out Dx 0 in 2Dx 0  0   0  D Dw Dw 1 Dc0 ¼   2Dy 0 a 2Dx 0 Dy 0 out Dy 0 in dc0 1 Dw0 Dc0 1 Dw0 Dc0 D þ   dt0 b 2Dy 0 2Dx 0 b 2Dx 0 2Dy 0 2Dx 0  0   0  D Dc Dc   ¼0 Dy 0 out Dy 0 in 2Dy 0



Dc0 Dx 0

ð15Þ 

 

out

Dc0 Dx 0

  in

ð16Þ

Each term of these equations may be considered as a transport variable – water flow for Eq. (15) and salt for Eq. (16) – and the equations themselves can be considered as the balance of these variables within the volume element. Defining

Figure 5

Figure 6

w0iþDx0 ;j  w0i;j iþDx0 ;j Dx0 w0i;j  w0iDx0 ;j jw0 0 ¼ iDx ;j Dx0 w0i;jþDy0  w0i;j jw0 0 ¼ i;jþDy Dy 0 0 wi;j  w0i;jDy 0 jw0 0 ¼ i;jDy Dy 0 0 ciþDx0 ;j  c0iDx0 ;j jw0i;j ¼ a c0iþDx0 ;j  c0i;j jc0 0 ¼ iþDx ;j Dx 0 c0i;j  c0iDx0 ;j jc0 0 ¼ iDx ;j Dx 0 jw0

¼

Simulation results. Salt concentration isolines for the original Henry problem, b = 0.1.

Simulation results. Salt concentration isolines for the modified Henry problem (Pinder version), b = 0.035.

45 ð17aÞ ð17bÞ ð17cÞ ð17dÞ ð17eÞ ð18aÞ ð18bÞ

46

A. Soto Meca et al. c0i;jþDy 0  c0i;j Dy 0

jc0

¼

jc0

c0i;j  c0i;jDy 0 ¼ Dy 0

i;jþDy 0

i;jDy 0

j1;c0 ¼ i;j

0 1 ðwi;jþDy0  wi;jDy 0 Þðc0iþDx0 ;j  c0iDx0 ;j Þ 2Dx 0 b

0 1 ðwiþDx0 ;j  wiDx0 ;j Þðc0i;jþDy 0  c0i;jDy 0 Þ i;j 2Dx 0 b oc0 jC0i;j ¼ 2Dy 0 ot

j2;c0 ¼

Eqs. (15) and (16) can be written as

ð18cÞ

jw0

 jw0

þ jw0

jc0

 jc0

þ jc0

iþDx 0 ;j

iþDx 0 ;j

ð18dÞ ð18eÞ ð18fÞ ð18gÞ

iDx 0 ;j

iDx 0 ;j

i;jþDy 0

i;jþDy 0

 jw0

i;jDy 0

 jc0

i;jDy 0

 jw0i;j ¼ 0

ð19Þ

þ j1;c0  j2;c0 þ jc0 ¼ 0 i;j

i;j

i;j

ð20Þ

From the point of view of the network model, the last two equations can be considered as Kirchhoff’s current law at the node (i, j) for the variables water density flux and salt density flux, respectively. On the one hand, the terms of equations (17a,b)–(18a,b) and (17c,d)–(18c,d) are easily implemented in the model by means of simple resistors of value R = 2Dx 0 and R = 2Dy 0 , respectively (Fig. 4). These resistors are of the same value since Dx 0 = Dy 0 . On the other hand, the term (18g) is associated with a condenser of value 2Dy 0 (F), while the nonlinear terms (17e), (18e) and (18f) are

Figure 7 Comparison of the solutions provided by different authors. Isoline c 0 = 0.5. Dashed line: Henry; dash-dot line: Lee; dotted line: Gotovac; triangle up (m): semianalytical, Clement; solid line with triangle down (.): Voss and Souza (SUTRA); solid line: present analysis.

Figure 8 Concentration c 0 = 0.5 (b = 0.035). Dotted line: Pinder and Cooper; dash-dot line: Segol et al.; dashed line: Voss and Souza (SUTRA); triangle up (m): Gotovac et al; solid line: present analysis.

An efficient model for solving density driven groundwater flow problems based on the network simulation method implemented by three controlled current sources, Gw0i;j , G1;c0i;j and G2;c0i;j . In this kind of source, the current may be defined by programming as a function of the potential variables (c 0 and/or w 0 ) of the other nodes of the model. The initial condition, Eq. (7e), is implemented by fixing the initial voltage of the condensers at the value zero. A total number of N · M volume elements are connected in series to represent the whole aquifer. Finally, the boundary conditions are implemented by means of simple devices. For example, Eqs. (7a) and (7c), which are related to the concentration and streamfunction at the boundaries (left and right and bottom and top, respectively), are implemented by means of constant voltage sources, while for Eqs. (7b) and (7d), simple resistors of a very high (infinite) value connected at the boundaries are needed to assume these adiabatic conditions.

47

The network model for the Elder problem is the same as that for the Henry problem, except for the boundary conditions, which are applied according to Eqs. (10)–(13).

Simulation and results Henry problem The general properties of the Henry problem are listed in Table 1. For the original Henry (1964) n = 2, a = 0.263 and b = 0.1; for the Pinder version (Pinder and Cooper, 1970) n = 2, a = 0.263 and b = 0.035, and for the modified Henry problem introduced by Simpson and Clement (2004) n = 2, a = 0.1315 and b = 0.2. Boundary conditions for the three cases are shown in Fig. 1. The large numerical dispersion that

Figure 9 Comparison of isochlors. c 0 = 0.25, 0.5 and 0.75. Solid line: present study, triangle up (m): semianalytical Simpson and Clement solution.

Table 3

x-Toe position of the isochlor c 0 = 0.5 for different authors and codes

b = 0.1

b = 0.035

b = 0.2

1.089 Henry (1964) 1.242 Lee and Cheng (1974) 1.409 Voss and Souza (1987) 1.400 Sutra (Segol, 1994) 1.375 Croucher and O’Sullivan (1995) 1.371 Gotovac et al. (2003) 1.393 Simpson and Clement (2004) Semianalytical 1.394 Seawat (Langevin and Guo, 1999) 1.373 Present study

1.220 Pinder and Cooper (1970) 1.245 Segol et al. (1975) 1.220 Voss and Souza (1987) 1.155 Sutra (Segol, 1994) 1.154 Gotovac et al. (2003) 1.158 Seawat (Langevin and Guo, 1999) 1.158 Present study –

1.074 Seawat (Langevin and Guo, 1999) 1.078 Simpson and Clement (2004) Semianalytical 1.074 Simpson and Clement (2004) Numerical 1.059 Present study –

–

–

– – –

48

A. Soto Meca et al.

Figure 10

Figure 11

Comparison of steady-state streamlines, b = 0.1. Solid line: present study, triangle up (m): Gotovac et al.

Comparison of steady-state streamlines, b = 0.035. Solid line: present study, triangle up (m): Gotovac et al.

Distance intruded (m)

1.4

0.25

1.2 1

0.5

0.8

0.75 0.6 0.4 0.2 0 0

100

200

300

400

Time (minutes)

Figure 12 Transient position of the intersection of the isochlors c 0 = 0.25, 0.5 and 0.75 with the base of the aquifer. Triangle up (dashed line): present study, continuous line: SEAWAT, Langevin and Guo (2006).

appears in the outflow region for the Pinder version makes this case very interesting to study. All the simulations were carried out for a 120 · 60 (equal volume elements) grid. Table 2 shows the minimum and maximum time-step used in the simulation for different grid sizes. This parameter always starts from an initial value (200E6 s for all cases in the table), which is defined as a small fraction of the total simulation time required. Minimum time-step occurs at the beginning of the simulation (periods of activity), where large changes in the variables take place, while the timestep increases for periods of inactivity. Fig. 5 shows the isolines of the salt concentration in the aquifer in the steady state for the original Henry problem, while Fig. 6 shows these curves for Pinder’s version of Henry’s problem. Isolines of both figures were depicted by interpolating the numerical values provided by simulation of the model, using the software MatLab 6.5.1. Curves move

An efficient model for solving density driven groundwater flow problems based on the network simulation method

49

Table 4 Computing times, number of nodes and toe for the isochlor c 0 = 0.5 as a function of the grid size (number of volume elements) Grid size

Computing time (s) b = 0.1

Computing time (s) b = 0.035

Computing time (s) b = 0.2

Number of nodes

Toe shoe (m) b = 0.1

Toe shoe (m) b = 0.035

Toe shoe (m) b = 0.2

20 · 10 40 · 20 80 · 40 120 · 60

2.160 55.59 1789 8621

3.36 84.42 2050 1028E+01

2.14 78.86 1829 8646

630 2460 9720 21,780

1.386 1.375 1.373 1.373

1.186 1.160 1.151 1.151

1.072 1.061 1.059 1.059

landward in the transition zone for both figures but slight differences between them can be appreciated due to the dispersion value (or, what is the same, due to the recharge rate), which is larger (smaller recharge rate) for the original Henry problem. As a consequence, the toe of the isochlors of the lower dispersion factor (Fig. 5) intersects the base of the aquifer at locations farther away from the seawater side, while the upper part of the curves contract in the seepage (discharge) area. For the original Henry problem, the steady condition is reached around a dimensionless time of 0.12, a value that comes from the solution of the transient concentration. The value used by Frind was 0.103 while that used by Pinder and Cooper, Segol et al. and Voss and Souza was 0.0396. We ensured steady results by extending all the simulations to a dimensionless time of 0.2. For the isochlor c 0 = 0.5, Fig. 7 compares the simulation results of the original Henry problem, which is quite similar to Croucher and O’Sullivan (1995) solution, with those of other authors (Henry, 1964; Lee and Cheng, 1974; Gotovac et al., 2003) and the semianalytical solution of Simpson and Clement (2004). Voss and Souza’s (1987) results, which used the Segol boundary condition and did not assume the Boussinesq approximation, are included in the figure. Our solutions are quite close to the solutions of Gotovac et al. and Simpson and Clement, authors who studied the Henry original problem in detail. Gotovac et al. solved the problem numerically from the streamfunction formulation, using a domain reduction technique that reduced the numerical dispersion, while Simpson and Clement provided a semianalytical solution. As is known, the Henry solution moved away from the above solutions (except for the discharge zone) because the Fourier series were truncated, eliminating high order terms. The same occurs for the numerical solution of Lee et al., for which, according to Frind (1982), some kind of numerical dispersion arises. Finally, the solution of Voss and Souza is clearly affected by the mixed boundary condition used and cannot be directly compared with the other solutions. For the modified Henry problem (Pinder version), Fig. 8 compares the isochlor, c 0 = 0.5, with the solutions of different authors. Again, the present solution, which is quite similar to the solution of Croucher and O’Sullivan (1995), is very close to the recent solution of Gotovac et al. (2003). However, since the solutions of Segol et al. (1975) and Voss and Souza (1987) used a mixed boundary condition, comparison with these authors are not suitable at the seepage. In the rest of the dispersion zone, the isochlor differs slightly from the present solution and that of Gotovac et al.

In Simpson and Clement (2004) the worthiness of the Henry problem as a benchmark problem is discussed by means of coupled versus uncoupled analysis. The objective of their work was to increase the importance of the densitydependent effects compared with the boundary forces. They proved that by decreasing the freshwater recharge by half, the worthiness of the Henry problem increases considerably and showed that the relative importance of density-dependent effects are spatially distributed within the plume. They presented new semianalytical and numerical solutions for this modified problem (which must not be confused with the one that uses the mixed boundary condition on the seaward boundary). Fig. 9 compares the above semianalytical results with the numerical solution of the present analysis which, in turn, is quite close to the numerical solution of Simpson and Clement (2004). The x-toe position given by different authors and codes for the three cases, b = 0.1, 0.035 and 0.2, is shown in Table 3. As seen, the agreement between the results of the present study (120 · 60 grid size) and those of the more recent authors is reasonable. Figs. 10 and 11 show the steady-state streamfunction for b = 0.1 (original Henry problem) and 0.035 (Pinder et al. version), respectively, and their comparison with the results of Gotovac et al. (2003), which are very close. The w 0 = 0 streamline separates the region of saltwater recirculation (w 0 < 0), which contains streamlines that begin and end at the sea boundary, from the upper region of non-recirculation (w 0 > 0), which contains fresh water through-flow, and meets the aquifer base at the stagnation point x 0 = 1.059, 0.925 and 0.410 (b = 0.1, b = 0.035 and b = 0.2, respectively). Figs. 10 and 11 include a comparison with the results of Gotovac et al., which are nearly the same. Finally, Fig. 12

Table 5 Symbol k n l g q0 Dq Df l d Ra

Properties and parameters of the Elder problem Value

Unit 13

4.845 · 10 0.1 1.0 · 103 9.807 1000 200 3.565 · 107 600 150 400

m2 – Pa s m s2 kg m3 kg m3 m2 s1 m m –

50

A. Soto Meca et al.

Figure 13 Left column: computed salinities, 0.1–0.9 isolines. Right column: streamline patterns for six dimensionless times, t 0 = 0.005, 0.01, 0.02, 0.05, 0.075 and 0.1.

An efficient model for solving density driven groundwater flow problems based on the network simulation method Table 6

51

Flow direction in the central section with respect to the mesh discretization

Discretization number of unknowns

Very coarse < 800

Coarse 1000–2000

Fine 3500–5000

Very fine 6000–10000

Olteam and Bue ´s (2001) Oldenburg and Pruess (1995) Kolditz et al. (1998) Diersch and Kolditz (2002) Frolkovic and De Schepper (2001) Present study

#" – – – – #

#" # # # # #

" " " " " "

" " " – – "

Figure 14 Coarse grid of 40 · 20, t 0 = 0.02. (a): computed salinity, 0.1–0.9 isolines; (b): streamlines pattern.

shows an additional test of the proposed model. This figure refers to the modified Henry problem of Simpson and Clement (2004) and successfully compares the solutions of the

present method (for a 40 · 20 grid size) and SEAWAT (Langevin and Guo, 1999) for the transient movement of the intersection of isochlors c 0 = 0.25, 0.5 and 0.75 with the base of the aquifer. A study of the grid size provided the results given in Table 4. Numerical dispersion nearly disappears for a grid size of 80(horizontal) · 40(vertical), except for the case b = 0.035 (Fig. 6) where a negligible dispersion still exists at the outflow boundary. This is also true for a 120 · 60 grid size (the dispersion disappears, for example, if the ten columns near the seaward boundary are changed to twenty). The network simulation method is able to deal with any kind of non-uniform grid as required by the problem. According to the requirements of the volume element techniques, the interface between the regions of large and small volume elements is generally formed by volume elements of suitable geometry – triangular, trapezoidal and other configurations can easily be implemented in the network method. Table 4 also shows the relatively low computation time for small 20 · 10 and 40 · 20 grids, and the longer time needed for large grids. The case b = 0.035, where numerical difficulties caused by the dominance of the convective terms are critical, needed the longest computation time. As regards the x-toe position for different grids (also shown in Table 4), it can be seen that all the grids provide almost identified results (errors less than 1% for grids of 40 · 20 and above).

Figure 15 Position of the isochlors c 0 = 0.2 and 0.6. Coarse grid of 40 · 20, t 0 = 0.05. Dashed line: Elder (1967b), dash-dot line: Voss and Souza (1987), continuous line: present study.

52

A. Soto Meca et al.

Elder problem

Conclusions

The values of the parameters and physical characteristics of the Elder problem are given in Table 5. A maximum of 1800 (60 horizontal · 30 vertical) rectangular volume elements were used for a half domain model, taking advantage of the symmetry of the problem. The movement of the solute in the closed domain is initiated exclusively by a diffusion process. The increase in water density due to changes in solute concentration promotes a circulation process, which disperses the salt throughout the aquifer. Fig. 13 shows the isochlors and streamfunction patterns for typical dimensionless times. Small eddies of reverse circulation emerge at the two ends of the source at t 0 = 0.005 and t 0 = 0.01. In the region of positive streamfunction values, flows circulate counterclockwise around the highest value of w, while in the region of negative streamfunction values, flows circulate clockwise around the lowest value of w. Complicated streamline patterns with two large and two small eddies can be observed on both sides of the domain at t 0 = 0.02. These eddies begin to interact and to fuse, and at the time t 0 = 0.05 the two small eddies disappear. Finally, for times longer than 0.005, the patterns basically remain unchanged. Oldenburg and Pruess (1995), who take the Elder experiment to verify the Tough2 code (Pruess, 1991), refer to two possible solutions: downwelling and upwelling. Kolditz et al. (1998) showed that for fine meshes (a spatial discretization consisting of at least 4400 bilinear finite elements that correspond to 4359 nodes) there is a central upwelling flow, while for coarse meshes there is a central downwelling flow. Their results, obtained with Feflow and Rockflow simulators (Diersch, 1994; Kolditz et al., 2002), agreed with those already provided by Oldenburg and Pruess. Using a mixed hybrid finite element method, Ackerer et al. (1999) always obtained the upwelling solution. Frolkovic and De Schepper (2001) presented results for the Elder problem based on a strict mesh convergence study. These results were independently confirmed by Diersch and Kolditz (2002), who also summarize the findings of recent research, indicating the flow direction in the central section with respect to the mesh discretization. For the fine mesh of 60 · 30 volume elements (which corresponds to 5490 nodes since each volume element has five nodes) the present method shows a clear central upwelling flow, high concentrations increasing in the upper region of the central part of the domain (Fig. 13). A mesh convergence study, where the mesh size is consecutively refined, 20 · 10, 40 · 20, 60 · 30 and 80 · 40, was carried out. An important outcome regarding the mesh level is the connection between the upwelling or downwelling flow and the size of the mesh (Diersch and Kolditz, 2002). Table 6 summarizes the findings of this study, which agree with those of other authors. As a part of this study, Fig. 14a and b, which correspond to a coarse mesh level of 40 · 20, show the computed salinity and streamline patterns corresponding to a downwelling central flow. Finally, Fig. 15 shows the results of Elder (1967b), Voss and Souza (1987) and the present study for a coarse grid of 40 · 20 (2500 nodes) at time t 0 = 0.05 (10 years).

An efficient 2-D numerical model based on the network simulation method is proposed for the solution of the nonlinear conjugate unsteady problems of Henry and Elder. The basic model for a (finite) control volume is an electric network containing two independent subcircuits, one for the streamfunction variable and the other for the salt concentration variable. The coupling between variables is realised by means of simple controlled (electrical) devices defined by programming. The complete model is obtained by electrically connecting the basic models and adding the boundary conditions by simple devices. Finally, the model can then be directly simulated by suitable circuit simulation software such as Pspice, which simultaneously provides the transport and potential variables at each point for the whole transient domain. Current mathematical manipulations inherent to this kind of numerical problems, such as numerical algorithms, output data organization, convergence criteria, etc., as well as additional requirements, such as those referring to mass conservation equations, are not necessary in this method since the software chosen does this work. Simulation of the Henry and Elder problems provided the unsteady and steady solutions for isochlors and streamline patterns. Other versions of Henry were also studied (such as those of Pinder and Cooper and the modified problem of Simpson and Clement). Both for the Henry (different versions) and Elder problems studied, the results were successfully compared with those of other authors (using a 120 · 60 volume element discretization for the Henry problem and 40 · 20 for the Elder problem), demonstrating the effectiveness and robustness of the proposed method. Finally, a study of the effect of grid size for all the cases studied in this paper enabled us to compare the results with those of other authors obtained for the Henry and Elder problems.

References Ackerer, P., Younes, A., Mose ´, R., 1999. Modeling variable density flow and solute transport in porous medium: 1. Numerical model and verification. Transp. Porous Media 35 (3), 345–373. Bear, J., Cheng, A.H.-D., Sorek, S., Ouazar, D., Herrera, I., 1999. Seawater Intrusion in Coastal Aquifers – Concepts, Methods and Practices. Kluwer Academic Publishers, London. Croucher, A.E., O’Sullivan, M.J., 1995. The Henry problem for saltwater intrusion. Water Resour. Res. 31 (7), 1809–1814. Diersch, H.J., 1981. Primitive variables finite element solutions of free convection flows in porous media. Zschr. Angew. Math. Mech. (ZAMM) 61, 325–337. Diersch, H.-J., 1994. Interactive, graphics-based finite-element simulation system – FEFLOW – for modelling groundwater flow and contaminant transport processes. Berlin. Diersch, H.J., Kolditz, O., 2002. Variable-density flow and transport in porous media: approaches and challenges. Adv. Water Resour. 25, 899–7944. Elder, J.W., 1967b. Transient convection in a porous medium. J. Fluid Mech. 27, 609–623. Evans, D.G., Raffensperger, J.P., 1992. On the streamfunction for variable-density groundwater flow. Water Resour. Res. 28 (8), 2141–2145. Frind, E.O., 1982. Simulation of long-term transient density-dependent transport in ground-water. Adv. Water Resour. 5, 73–88.

An efficient model for solving density driven groundwater flow problems based on the network simulation method Frolkovic, P., De Schepper, H., 2001. Numerical modelling of convection dominated transport with density driven flow in porous media. Adv. Water Resour. 24 (1), 63–72. Galeati, G., Gambolati, G., Neuman, S.P., 1992. Coupled and partially coupled Eulerian–Lagrangian model of freshwater– seawater mixing. Water Resour. Res. 28 (1), 149–165. Gonza ´lez-Ferna ´ndez, C.F., Alhama, F., 2002. In: Horno, J. (Ed.), Heat Transfer and the Network Simulation Method. Research Signpost, Trivamdrum. Gotovac, H., Andricevic, R., Gotovac, B., Kozulic, V., Vranjes, M., 2003. An improving collocation method for solving the Henry problem. J. Contam. Hydrol. 64, 129–149. Henry, H.R., 1964. Effects of dispersion on salt encroachment in coastal aquifers. In: Sea Water in Coastal Aquifers, US Geol. Surv. Supply Pap. 1613-C, pp. 70–84. Holzbecher, E., 1998. Modelling Density-Driven Flow in Porous Media. Springer, Berlı´n. Horno, J., Gonza ´lez-Ferna ´ndez, C.F., Hayas, A.J., 1995. The network method for solutions of oscillating reaction–diffusion systems. Comput. Phys. 118, 310–319. Huyakorn, P.S., Andersen, P.F., Mercer, J.W., White Jr., H.O., 1987. Saltwater intrusion in aquifers: development and testing of a three dimensional finite element model. Water Resour. Res. 23, 293–312. HYDROCOIN, 1990. The International HYDROCOIN Project, Level 2: Model validation, OECD, Paris, 194p. Kolditz, O., Ratke, R., Diersch, H.-J.G., Zielke, W., 1998. Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models. Adv. Water Resour. 21 (1), 27–46. Kolditz, O., Habbar, A., Kaiser, R., Rother, T., Thorenz, C., 2002. ROCKFLOW—Theory and User’s Manual. Release 3.6, Institute of Fluid Mechanics. University of Hannover, Germany. Langevin, C.D., Guo, W., 1999. Improvements to SEAWAT, a variable-density modeling code [abs.] in Eos, Transactions, 80 no. 46, P. F-373. Langevin, C.D., Guo, W., 2006. MODFLOW/MTÆDMS-Based simulation of variable-density ground water flow and transport. Ground Water 44 (3), 339–351.

53

Lee, C.-H., Cheng, T.-S., 1974. On seawater encroachment in coastal aquifers. Water Resour. Res. 10 (5), 1039–1043. Oldenburg, C.M., Pruess, K., 1995. Dispersive transport dynamics in a strongly coupled groundwater-brine flow system. Water Resour. Res. 31, 289–302. Olteam, C., Bue ´s, M.A., 2001. Coupled groundwater flow and transport in porous media. A conservative or non-conservative form? Transp. Porous Media 44 (2), 219–246. Pinder, G.F., Cooper Jr., H.H., 1970. A numerical technique for calculating the transient position of the saltwater front. Water Resour. Res. 6 (3), 875–882. Pruess, K., 1991. TOUGH2 – a general-purpose numerical simulator for multiphase fluid and heat flow. Lawrence Berkeley Laboratory Report, LBL – 29400. PSPICE 6.0, 1994. Microsim Corporation Fairbanks, Irvine, CA. Segol, G., 1994. Classic Groundwater Simulations Proving and Improving Numerical Models. Prentice-Hall, Old Tappan, NJ. Segol, G., Pinder, G.F., Gray, W.G., 1975. A Galerkin-finite element technique for calculating the transient position of the saltwater front. Water Resour. Res. 11 (2), 343–347. Senger, R.K., Fogg, G.E., 1990. Streamfunctions and equivalent freshwater heads for modelling regional flow of variable-density groundwater, I, Review of theory and verification. Water Resour. Res. 26 (9), 2089–2096. Simpson, M.J., Clement, T.P., 2003. Theoretical analysis of the worthiness of the Henry and Elder problems as benchmarks of density-dependent groundwater flow models. Adv. Water Resour. 26, 17–31. Simpson, M.J., Clement, T.P., 2004. Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models. Water Resour. Res. 40, W01504. Voss, C.I., Souza, W.R., 1987. Variably density flow and solute transport simulation of regional aquifers containing a narrow fresh-water saltwater transition zone. Water Resour. Res. 23 (10), 1851–1866. Zueco, J., Alhama, F., Gonza ´lez-Ferna ´ndez, C.F., 2004. Analysis of laminar forced convection with network simulation in thermal entrance region of ducts. Int. J. Therm. Sci. 43, 443–451.