Efficient conceptual model for simulating the effect of aquifer heterogeneity on natural groundwater discharge to rivers

Advances in Water Resources 34 (2011) 1377–1389

Contents lists available at SciVerse ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Efficient conceptual model for simulating the effect of aquifer heterogeneity on natural groundwater discharge to rivers David Pulido-Velazquez a,⇑, Carlos Llopis-Albert b,1, Salvador Peña-Haro c, Manuel Pulido-Velazquez b,2 a

Instituto Geológico y Minero de España (IGME), Urb. Alcázar del Genil, 4-Edif. Zulema, Bajo. 18006 - Granada, Spain Research Institute of Water and Environmental Engineering, Universitat Politècnica de València, Spain c Institute of Environmental Engineering, ETH Zürich, Zürich, Switzerland b

a r t i c l e

i n f o

Article history: Received 11 January 2011 Received in revised form 20 July 2011 Accepted 20 July 2011 Available online 31 July 2011 Keywords: Natural discharge Heterogeneity Eigenvalue techniques Conceptual model Linear reservoir

a b s t r a c t When linearity can be assumed (linear response of heads to stresses), stream–aquifer flow exchange can be simulated as the drainage of a number of independent linear reservoirs. This conceptual model, which can be mathematically deduced in a univocal way from an eigenvalue solution of the linear groundwater flow problem, facilitates the understanding of the physical phenomenon and the analysis of influencing factors. The number of reservoirs required to simulate stream depletion in some ideal homogeneous cases of stream–aquifer connection was analyzed in detail in a previous investigation using analytical eigenvalue solutions [16]. However, most aquifers are heterogeneous in nature and numerical solutions must be employed to analyze whether they could also be simulated using few reservoirs. This paper presents a stochastic analysis of the influence of heterogeneity on the simulation of natural groundwater discharges in aquifers connected to rivers, as a series of linear reservoirs. A Monte-Carlo approach was employed to perform this study. The results show that, on a monthly time scale, many cases (even heterogeneous aquifers) can be simulated using just a few reservoirs with sufficient accuracy and at minimum computational cost. Therefore, this modeling technique can be useful to efficiently simulate the integrated management of complex water resources systems at the basin scale (with many aquifers, reservoirs, demands, etc.) that need to simultaneously consider surface and groundwater flow and stream–aquifer interaction. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction The interest in modeling processes and the dynamics of groundwater (GW)–surface water (SW) interactions using different approaches has grown steadily over the last two decades. Many GW–SW interaction models have been developed for use in analysing water quantity management problems at different spatial– temporal scales. In these cases, different techniques have been increasingly used to explore hypotheses and to develop new conceptual models [5]. This is the case of eigenvalue techniques. Eigenvalue techniques have been used to derive analytical (e.g. [10,18,16,3]) and numerical solutions of groundwater flow and stream–aquifer interaction problems (e.g. [7,21,1,22,13,14]) that have helped to develop conceptual models. When linearity can be accepted (linear response of heads to stresses), the solution of ⇑ Corresponding author. Tel.: +34 963 943 474; fax: +34 963 944 436. E-mail addresses: [email protected] (D. Pulido-Velazquez), [email protected] (C. Llopis-Albert), [email protected] (S. Peña-Haro), [email protected] (M. Pulido-Velazquez). 1 Tel.: +34 963877058; fax: +34 963877618. 2 Tel.: +34 963879616; fax: +34 963877618. 0309-1708/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2011.07.010

the stream–aquifer flow exchange over a stress period can be formulated using an explicit and continuous in time (time is not discretized) state equation whose structure shows that stream– aquifer interaction problem can be simulated as the drainage of a finite number of linear reservoirs (LR). This conceptual approach, called the embedded multireservoir model (EM model) [17,15], presents computational advantages (whose importance increases as the number of reservoirs required for the simulation reduces) and facilitates the understanding of the physical phenomenon and its influencing factors. Analytical EM solutions are obtained without either discretizing the spatial domain nor the stress periods, and they solve the partial differential equation of the GW–SW problem explicitly (neither time nor space are discretized), which is only possible for ideal cases (homogeneous aquifers with simple geometry and boundary conditions) [16]. Meanwhile, numerical EM solutions can be obtained assuming a spatial discretization in order to approach the spatial derivatives of the governing equation, and solving explicitly with respect to the time variable but not for the spatial domain. Spatial discretization allows approximation of the solution of the groundwater flow equation for more complex aquifers (heterogeneous aquifers with irregular domains and more complex

1378

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

boundary conditions). The most common way of discretizing the spatial domain is using a grid, as in classic finite difference (FD) models. It has been demonstrated that the EM model allows simulation of the GW–SW interactions (two-way interactions) obtained from any FD linear flow model, providing practically the same results, since both solve the same spatial approximation to the flow problem (same geometry, boundary conditions, hydrodynamic parameters and stresses defined over the grid) but without discretizing the stress periods [15]. A disadvantage of EM models with respect to FD models is that they cannot solve non-linear flow problems (problems governed by a non-linear partial equation and/or with non-linear boundary conditions). The eigenvalue technique employed to define the EM models allows the flow problem to be solved in a different way to the FD models, and it yields a solution whose structure (a state equation) provides certain computational advantages [15]. These computational advantages have been stated and evaluated in various previous papers [1,14]. They are mainly due to: – FD models require division of the stress periods (i.e. the time interval over which all the external stresses remain constant) into smaller time steps and obtain a separate solution for each, while the EM model calculates the solution for a particular stress period directly. – The FD model requires the hydraulic head in each cell to be obtained for each time step, while the EM model allows GW–SW flow variables to be calculated directly, without obtaining the hydraulic head in each cell. The computational advantages increase as the number of terms (number of LRs) required to obtain the EM solution decreases. The possibility of simulating stream–aquifer interactions using a reduced number of LRs could be useful to simulate or optimize the management of complex water resources systems at basin scale (with several aquifers, reservoirs, demands, etc.) that simultaneously consider surface and groundwater flow and the interaction between them for long time horizons and monthly stress periods [2,17]. The number of LRs required by the EM models for some ideal homogeneous cases (stream depletion produced by pumping or recharge in a rectangular homogeneous aquifer perfectly connected with a river along one of its sides) has been analyzed using an analytical solution [16]. Analytical solutions are often useful for a preliminary assessment and, in this case, they demonstrate that stream–aquifer interaction can be simulated with a reduced number of LRs in some homogeneous cases (the exact solution is closed-form and we do not need to use the infinite terms of the exact solution). But numerical solutions are required to analyze more complex problems, as in the case of heterogeneous aquifers. In a previous investigation [15], numerical EM models were employed to simulate stream–aquifer depletion produced by constant pumping in certain ideal cases. However, most aquifers are heterogeneous and, from a practical point of view, it would be interesting to study whether GW–SW problems can also be approached using a small number of LRs, regardless of how complex the heterogeneity. A detailed analysis of the effect of the heterogeneity is needed to extract conclusions in this respect. This paper focuses on an analysis of the impact of heterogeneity on an EM simulation of a particularly interesting GW–SW problem relating to the natural aquifer discharges to a river. It is a significant problem, as highlighted by the great number of papers published about this issue [4,12,23]. In this paper, in order to assess whether this heterogeneous aquifer problem can be approached using a reduced number of LRs, we develop a stochastic study using a Monte-Carlo approach, in which we analyse the effect of recharge events that are evenly (homogeneously) spread over the

aquifer (a hypothesis made to represent cases of natural aquifer recharge) for multiple synthetic case studies generated using different degrees of heterogeneity. Although the study assumes uniform recharge (in order to solve problems of discharge under a natural regime), the linearity of the problem and the solution proposed (EM model) can be applied to modeling pumped abstractions that are evenly spread over the aquifer (which are equivalent but with the opposite sign). Thus, the conclusions are applicable to cases with two-way interactions between groundwater and stream. In order to generate the heterogeneous case studies we focused on the hydraulic conductivity (K), since this is usually the parameter with the most significant spatial variation. In fact, it can vary spatially by several orders of magnitude. For instance, the aquifer at the Columbus Air Force Base in Mississippi, commonly known as the Macrodispersion Experiment (MADE) site, is a strongly heterogeneous system with a variance of the natural logarithm of K of nearly 4.5 [19,20,9]. In addition, fractured rock aquifers are also characterized by high variance in conductivity [8]. In this paper, we wanted to study a wide range of heterogeneities and analyse the sensitivity of the solution to the hypothesis adopted in the definition of our case study. For this reason, we consider that the use of synthetic case studies provides more flexibility in achieving the study’s objective (testing the efficiency of the solution over a wide range of heterogeneity). Nevertheless, we have attempted to use the analysis performed with the selected synthetic cases to provide a preliminary idea of the necessities for any real problem. The paper has been organized into the following sections: Methods (with description of the Monte-Carlo methodology and the cases studies adopted to allow the work to be reproducible); Theory (where we include an extension of the background, with a summary of the EM formulations and we comment on the foundation for future research); Results and discussion (including the results obtained to test the accuracy of the EM models by comparing EM solutions with those obtained using MODFLOW [11] and we explore their significance and comment on their limitations and possible future research); and, lastly, the Conclusions reached. 2. Methodology: the Monte-Carlo approach A stochastic Monte-Carlo approach is used to analyze the influence of heterogeneity in the definition of an efficient conceptual– numerical EM model to simulate natural discharges of aquifers connected to a river. The approach is built systematically using the following steps: first, multiples case studies are defined in order to cover a wide range of degrees of heterogeneity, using a constant recharge scenario to approximate their behaviour under a natural regime (2.1); then, MODFLOW is used to simulate groundwater flow and stream–aquifer interaction for these case studies (2.2); different EM approaches are applied to simulate each case (2.3); finally, the natural discharges obtained with MODFLOW and with the EM models are compared. 2.1. Definition of the heterogeneous case studies In order to analysis the effect of heterogeneity, multiple case studies are defined with different levels of heterogeneity (null [homogeneous cases], intermediate and extreme) for the geometry, discretization and boundary conditions (BC) represented in Fig. 1. A no-flow boundary condition (Neumann BC) is selected to define the aquifer’s limits, while a head-dependent flow boundary condition (a Cauchy BC) is used to simulate the aquifer discharge. This BC is formulated as the product of the stream–aquifer hydraulic connection (defined by a conductance parameter [L2/T]) multiplied by the hydraulic head gradient (difference between the stream head and the aquifer head in the places where the BC is defined).

1379

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389 Table 1 Heterogeneous case studies considered. Perfect stream–aquifer connection (20,000 m2/day) Heterogeneous cases High heterogeneity level (var[LnT] = 4.5) Mean heterogeneity level (var[LnT] = 0.5) Partial stream–aquifer connections (200 m2/day) High heterogeneity level (var[LnT] = 4.5)

Fig. 1. Discretization and boundary conditions.

An homogeneous storage coefficient equal to 0.0001 is considered in all the cases. The transmissivity, usually the parameter with most spatial variation in an aquifer [19,20,9], was defined for different levels of heterogeneity. Each case study with a specific heterogeneity level requires the generation of multiple transmissivity (T) fields to carry out a stochastic analysis of the results obtained by using a Monte-Carlo approach. One hundred unconditional T fields were generated for each case study by means of sequential Gaussian simulation using the computer code GCOSIM3D [6]. It is assumed that this set is large enough to provide a significant representation of the variability of this parameter. The same stochastic structure was considered for all simulated T fields; in this way, we simplify our analysis and avoid uncertainty about the stochastic structure. Consequently, all T fields are equally likely realizations and therefore, plausible representations of reality, since they display the same degree of spatial variability. The stochastic

Aquifer diffusivity (T/S m2/day) 10^7

10^5

10^3

10

Case I

Case II

Case III

Case IV

Case I’

Case II’

CASE Ip

structure was defined using a spherical variogram with a range approximately equal to 1/5 of the aquifer size. The different degrees of heterogeneity were generated by defining two different variances of the natural logarithm of T (equal to 0.5 for the intermediate degree of heterogeneity and 4.5 for the extreme one). Previously developed research using analytical solutions showed that, in simple homogeneous aquifers under homogeneous distributed stress, the number of LRs required to simulate stream– aquifer interaction depends on the stream–aquifer hydraulic connection and the aquifer diffusivity [L2/T], which is defined as the transmissivity divided by the storage coefficient. For this reason, in order to analyse the influence of diffusivity on the results, we considered different values of mean transmissivity (1000 m2/day; 10 m2/day; 0.1 m2/day and 0.001 m2/day; obtained for each case study as the mean transmissivity in each cell), which (with constant homogeneous storage equal to 0.0001) define cases with different diffusivity (107 m2/day; 105 m2/day; 103 m2/day; and 10 m2/ day). The homogeneous case studied also showed that the cases requiring a higher number of LRs are related to perfect stream– aquifer connection [16]. It seems logical that this also happens for heterogeneous cases (cases with higher stream–aquifer connection will be related to a need for a higher number of LRs). Therefore, we assumed perfect GW–SW connection (stream–aquifer conductance equal to 20,000 m2/day) to define the synthetic case studies. In order to check the sensitivity of the results to the GW–SW hydraulic connection we also solved a case for partial stream–aquifer connection (conductance equal to 200 m2/day). The synthetic case studies finally considered are summarized in Table 1. Cases I, II, III, and IV (see Table 1) were defined with high heterogeneity (variance of the natural logarithm of T equal to 4.5) and various mean diffusivity values in order to assess the influence of the mean diffusivity. To analyze the sensitivity of the results to the degree of heterogeneity, two further cases were defined with an intermediate degree of heterogeneity (variance of the natural logarithm of T equal to 0.5): case I’ and II’ (see Table 1). A case study with partial connection was also analyzed to analyse the sensitivity to stream–aquifer connection (case Ip that differs from case I in the stream–aquifer connection, which is equal to 200 m2/ day). In addition, we also solved a number of homogenous cases to try to identify a relationship between the number of terms required in the heterogeneous cases and in the homogeneous ones. In order to compare results with the heterogeneous cases we considered the same constant storage equal to 0.0001 and homogeneous transmissivity values equal to the mean heterogeneous transmissivities (1000 m2/day; 10 m2/day; 0.1 m2/day; 0.001 m2/ day and 0.00001 m2/day; obtained for each case study as the mean transmissivity value in each cell), which (using the constant homogeneous storage coefficient equal to 0.0001) define cases with the following diffusivity values: 107 m2/day; 105 m2/day; 103 m2/day; 10 m2/day and 0.01 m2/day.

1380

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

2.2. Simulations: variables, models and scenarios

Table 2 Hydrodynamic parameters for the homogeneous case studies.

In our study, we focussed on the variable instantaneous stream–aquifer flow exchange rate at the end of each stress period. The analysis was done on a monthly scale, which is the most common temporal scale employed in river basin management models). Since the objective is to analyze the accuracy of the EM simulation, we compared our EM results with results from MODFLOW, since the accuracy of this finite difference method has been widely demonstrated. For this reason, we defined a variable to study the relative differences (E(%)SOLUCION(j)) between the results obtained from MODFLOW [11] (CVMODFLOW), and those obtained with EM models (CVSOLUCION(j)) over a stress period t, which is evaluated using the expression

E%ðtÞSOLUCIONðjÞ ¼

CVðtÞMODFLOW  CVðtÞSOLUCIONðjÞ CVðtÞMODFLOW

 100

ð1Þ

For each field of these case studies we simulated a scenario that represents the influence of a change in natural recharge. Since the problem is linear, and we intended to analyze relative errors in the simulation of changes in natural recharge for different EM approaches, a simplified scenario with constant uniform recharge, zero initial conditions and zero head along the river was simulated. The simulation was conducted over a period long enough to reach a steady state. This scenario was simulated with MODFLOW for the multiple fields of heterogeneous cases defined (see Section 2.1). The stream–aquifer relationship was modelled using the General-Head Boundary (GHB) package. MODFLOW is the most commonly used code in groundwater flow, and its accuracy to simulate discrete groundwater flow problems defined over finite difference grids has been widely demonstrated. Therefore, results provided by MODFLOW were considered appropriate to estimate the relative errors that appear when other approaches are used. We also simulated the same scenario (for each field of the case studies) using different EM numerical approaches, which are defined as the summation of different numbers of LR (from just one up to the number of active cells in the spatial discretization of the aquifer geometry). The parameters and initial conditions of each LR are mathematically defined in a univocal way from the original finite difference model, preserving its geometry, boundary conditions, hydrodynamic parameters (heterogeneity) and spatial distribution of stresses (see summary of the formulation in Section 3). The stream–aquifer evolution obtained with MODFLOW and the EM approaches was then employed to estimate the relative differences between the two solutions, as explained above (Eq. (1)). 2.3. Analysis of results For each field of the case study we obtained the number of LR required to maintain the relative error (E(%)SOLUCION(j), see Eq. (1)) below a certain threshold (1% and 5%) during the simulation. Statistical analysis of the results (number of LRs required) for all the fields in each case study identifies the number of LRs required to simulate heterogeneous cases, whilst maintaining the error below a certain threshold (which is the main objective of this research, to check if the heterogeneous cases could be also simulated with a reduced number of LRs). A sensitivity analysis (to stream–aquifer connection, mean diffusivity parameters and heterogeneity) was done using the results from the different case studies. We also solved a number of homogenous cases in order to compare results with the heterogeneous cases, in an attempt to identify an experimental relationship between the number of LR required in homogeneous vs. heterogeneous cases. The homogeneous cases studied (defined by modifying the diffusivity and the stream–aquifer connection) are included in Table 2.

Stream–aquifer connection (m2/day) Homogeneous cases 20,000 (‘‘perfect connection’’) 200

Aquifer diffusivity (T/S m2/day) 10^7

10^5

10^3

10

0.1

Hom 1

Hom 2

Hom 3

Hom 4

Hom 5

Hom 6

Hom 7

Hom 8

Hom 9

Hom 10

In a previous article using analytical solutions [16], relative error was defined and studied as a percentage of the intensity of the applied stress: CVðtÞMODFLOW CVðtÞSOLUCIONðjÞ ðE%ðtÞSOLUCIONðjÞ ¼  100, where R is the reR charge rate). In this study, we are interested in limiting the error (defined in accordance with Eq. (1)) in the simulation of natural aquifer discharges with respect to an ‘‘exact solution’’ (obtained by simulation with MODFLOW) in each of the monthly time steps (denoted by j in Eq. (1)). For this reason, the number of terms required in this analysis will be higher than in that previous research [16]. This stochastic analysis has been described in detail in Section 5. 3. Theory: EM numerical–conceptual models The conceptual–numerical EM model employed to study the influence of heterogeneity on the number of LRs required to simulate with sufficient accuracy, was deduced in a previous investigation [15] and its formulation is summarized in this section. It was deduced by comparing the conceptual–analytical eigenvalue solution (the EM model, see [16]) of stream–aquifer interaction and the numerical eigenvalue method solution for distributed-parameter groundwater-flow problems [1]. It allows solution of the GW–SW interaction (two-way interactions) obtained using a calibrated finite difference (FD) model through the drainage of some independent LRs (see Fig. 2) in the same way that an analytical multireservoir model does [16]. The parameters of the multi-reservoir model (discharge coefficients and allocation factors) can be estimated in a univocal way from the original calibrated FD model. The model preserves the aquifer characteristics considered in the distributed model used to define it (geometry, boundary conditions, hydrodynamic parameters, and spatial distribution of stresses). Its formulation is summarized as follows. The linear 2D groundwater flow equation, defined with timeinvariant transmissivity, can be expressed as

LðhÞ þ w ¼ S 

@h ; @t

ð2Þ

where h = h(x, y, t) is the hydraulic head, [L]; S(x, y) is the dimensionless storage coefficient; w(x, y, t) is the external stress function  @  @  (pumping or recharge rate per unit area), [L/T]; LðÞ ¼ @x T x @x þ    @ @ T is a linear operator in which T (x, y) and T (x, y) are the y x y @y @y main components of the transmissivity tensor [L2/T]. These components can be considered invariant over time for aquifers with linear behaviour (confined aquifers or unconfined aquifers with small hydraulic head changes compared to their saturated thickness). When the spatial domain is discretized using a finite difference grid, the linear groundwater flow problem defined from Eq. (2) and linear boundary conditions can be formulated with the following matrix expression [21]

T  fHg þ fWg ¼ SF 

dfHg ; dt

ð3Þ

1381

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

Fig. 2. Conceptualization of the EM model for a single stress (modified from Pulido-Velazquez et al. [15]).

where {H} is the hydraulic head vector, with (N  1) dimensions, N is the number of active cells in the aquifer discretization [L]; T is a symmetrical matrix (N  N) whose elements depend on the transmissivity values, space discretization and boundary conditions [L2/ T]; SF is a diagonal matrix (N  N) whose elements depend on storage coefficients and space discretization [L2]; {W} is an invariant vector (N  1) in each stress period that depends on the external stresses and non-zero Dirichlet or Cauchy boundary conditions [L3/T]. The solution over a stress period s(t  1, t) of duration Dt during which the external stress w(x, y) remains invariant can be obtained by solving the eigenproblem [21]:

T  A ¼ a  SF  A;

ð4Þ

where a and A are the diagonal eigenvalue matrix (N  N), [T1] and the corresponding eigenvector matrix (N  N) of the eigenproblem [L1]. The stresses can be expressed as a linear combination of some unit stresses (basic stresses) with time-invariant spatial distribution:

fWðsÞg ¼ D  fRðsÞg;

ð5Þ

where D is a matrix whose columns are defined with the unitary percentage of the stresses applied among the discretization cells (N  Na), and Na is the number of basic stresses defined; {R(s)} is the intensity at which the basic vectors are applied (Na  1) during the stress period s that ends at time t [L3/T]. The surface water body-aquifer flow exchange rate can be approximated with the drainage of p + 1 independent LRs formulated as follows, with truncated matrices of p + 1 components (considering only the terms corresponding to the p + 1 smallest eigenvalues, where p + 1 < N, and N is the number of active cells in the spatial discretization), which has been denoted with the superscript r.

QðtÞ ¼ far g  fVðtÞg; fVðtÞg ¼ fEðDtÞgr  fVðt  1Þgr þ



ðIr  EðDtÞr Þ



ar

 fB  fRðsÞgg ¼ fT1g  fVðt  1Þgr þ fT2g  fBr  fRðsÞgg; r

r

ð6Þ

where {V(t)}r is a vector with p + 1 components representing the storage in each LR at time t. {T1}r = {E(Dt)}r is a vector of (p + 1) components that represents the decline of the initial conditions, whose ith component is defined by T1(i) = ei = exp(a(i)  Dt). ar, Ir, E(Dt)r are the reduced diagonal matrices defined by the first p + 1 components of the diagonal matrices a (see Eq. (6)) that represent the discharge coefficients of the LR, I (identity matrix) and E(Dt) (whose no nulls components are eii = exp (a(i)  Dt)). {T2}r is a vector of (p + 1) components employed to obtain the effect produced by the applied stresses. Its ith component is defined by T2(i) = (1  exp( a(i)  Dt))/(ai). {V(t  1)}r is a vector with p + 1 components representing the storage in each LR at time t  1. The initial storage condition should be exactly reproduced by the model and so it must be estimated using the following expression

V ri ð0Þ ¼ fi  li ð0Þ if i < p;

V rp ð0Þ ¼ Vð0Þ 

p1 X

fi  li ð0Þ;

ð7Þ

i¼0

where

fir

¼

" N X

# SF j ðaji Þ

for i 6 p

ð8Þ

j¼1

and aji represents the jth component of the eigenvector {ai}, and the fir terms represent the storage beneath the surface defined by each eigenvector {ai} that corresponds to the p + 1 smallest eigenvalues. r

li ð0Þ ¼

" N X

#

ðaij Þ  SF jj  hj

for i 6 p:

ð9Þ

j¼1 m

Br is the matrix of allocation factors bi that represent the percentage of stresses applied in each LR. The components in each column (there are as many columns as stresses) represent the allocation factors that correspond to each stress, m, and can be estimated in a univocal way from a calibrated FD model as explained below. The fluid mass balance must be preserved as in the truncated conceptual analytical solution. As a consequence, the stresses must

1382

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

be totally distributed among the LRs considered in the conceptual Pp m model (p + 1). This is mathematically formulated as i¼0 bi ¼ 1. Therefore, the last row in the matrix of allocation factors has to be obtained by imposing this requirement. The other components, m bi , will be defined by the following equation:

Br ¼ Fr  AT;r  D if

m

i < p and bp ¼ 1 

p1 X

m

bi ;

ð10Þ

i¼0

where Fr is a diagonal matrix whose ith diagonal component is obtained applying Eq. (8), AT,r ((p + 1)  N) is the transpose matrix of the truncated eigenvalue matrix, which is defined with the eigenvectors related to the p + 1 smallest eigenvalues. This formulation can be applied to simulate two-way interactions of GW–SW for any kind of stress applied (point pumped abstraction or recharge and distributed pumping or recharge) whenever the flow problem is linear (the governing equation and the BC are linear). Nevertheless, in order to limit this research to a manageable problem, as explained in the introduction, we focused on cases with homogenous distributed stresses usually applied to solve natural recharge in an aquifer. Although, in this paper, we consider solely discharges under a natural regime – assumed to occur by evenly-distributed recharge, the linearity of the problem and its proposed solution (using EM models), means that the same conclusions are applicable to modeling abstractions that are evenly-distributed over the aquifer (equivalent, though with the opposite sign). Therefore, the conclusions are applicable

to problems with two-way interactions between groundwater and stream (The river could recharge the aquifer in cases with distributed pumping.) Previous research has shown that, for homogeneous cases [16] where the problem is linear, point action (pumping from or recharge into wells) can also be modelled using EM models, but the effect of heterogeneity on the number of terms is unknown [15]. A further IMP investigation to study the impact of heterogeneity in these cases and we consider that this is of sufficient importance to warrant further investigation.

4. Results and discussion 4.1. Heterogeneous case studies For each heterogeneous case study, the designed Monte-Carlo approach (see Section 2 for details about the case studies, the simulations performed [variables, models employed and scenarios simulated], and results selected for this analysis) was applied to obtain the cumulative density function (CDF) of the variable ‘‘number of LRs’’ required to simulate natural aquifer discharges, whilst maintaining the maximum relative error (Eq. (1)) beneath a certain threshold (5% and 1%). The CDF obtained for the heterogeneous fields defined under different hypotheses (cases summarized in Table 1) are represented in Fig. 3a (In Fig. 3b the x-axis scale has been changed to show the CDF of cases I, I’, Ip, II and II’ in greater detail).

Linear reservoirs required to simulate aquifer discharge rate at monthly scale with maximum relative error of 5%

Cumulative density function

1

CASE I CASE I' CASE Ip CASE II CASE II' CASE III CASE IV

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000

Number of linear reservoirs Linear reservoirs required to simulate aquifer discharge rate at monthly scale with maximum relative error of 1%

Cumulative density function

1

CASE I CASE I' CASE Ip CASE II CASE II' CASE III CASE IV

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000

Number of linear reservoirs Fig. 3a. CDF of the number of LRs required to simulate aquifer discharge on a monthly scale using EM models, for maximum relative errors of 1% & 5% (cases I, I’, Ip, II, II’, III, IV).

1383

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389 Linear reservoirs required to simulate aquifer discharge rate at monthly scale with maximum relative error of 5%

Cumulative density function

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 CASE I

CASE I'

CASE II

CASE II'

CASE Ip

0.1 0 0

10

20

30

40

50

60

70

80

90

100

110

Number of linear reservoirs Linear reservoirs required to simulate aquifer discharge rate at monthly scale with maximum relative error of 1%

Cumulative density function

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

CASE I

CASE I'

CASE Ip

CASE II

CASE II'

0 0

10

20

30

40

50

60

70

80

90

100

110

Number of linear reservoirs Fig. 3b. CDF of the number of LRs required to simulate aquifer discharge on a monthly scale using EM models, for maximum relative errors of 1% & 5% (cases I, I’, Ip, II, II’).

The analysis of selected results allows us to draw the following conclusions about the sensitivity of results to the predefined admissible error, mean aquifer diffusivity, heterogeneity, and stream–aquifer connection: – If we reduce the admissible error, the number of LRs required is higher and its dispersion (the degree of scatter of data considered as the magnitude of the range over which the variable takes values) decreases (see Fig. 3 for all case studies with maximum relative errors of 1% and 5%). This was also shown for other thresholds in case studies II and II’ (see Fig. 4). This conclusion, confirmed with the Monte Carlo approach, could be expected a priori: the more accurate the solution required, the greater the number of LRs used in the EM approach. As explained in the theory section, EM solutions are obtained as a summation of LRs whose parameters (discharge coefficients and allocation factors) can be obtained in a univocal way from a FD model. In this summation, the LRs related to the smallest discharge coefficients have greatest influence on the final solution. The terms T1 and T2 that depend on the temporal scale and the discharge coefficients [T1[a(i), Dt] = exp( a(i)  Dt) and T2[a(i), Dt] = (1  exp (a(i)  Dt))/(ai)] decrease as the discharge coefficients increase. In order to show this aspect for one field of case I and case IV, Fig. 5 shows the value adopted by functions T1 and T2 for the first ten LRs at different temporal scales (monthly and daily). Fig. 5 shows how the smallest discharge coefficients are related to the highest values of the functions T1 and T2, and, therefore, the LRs related to the smallest discharge coefficients will have the greatest influence

on the solution. Although the relative importance of the LRs related to the smallest eigenvalues is small, all of them add some information (the exact solution of the problem is obtained using the summation of all the LRs, as many as the number of active cells). Therefore, if we demand a more accurate solution, it should be related with an increment in the number of LRs required. – The cases defined with mean diffusivity values equal to or higher than 105 m2/day significantly reduce the number of LRs required to simulate with sufficient accuracy (see Fig. 3). We can analyze this conclusion by considering the formulation of the EM solutions. As described in the theory section, the LRs related to the smallest discharge coefficients have the greatest influence on the final solution (the T1 and T2 terms decrease as the discharge coefficients decreases, see Section 3). In homogeneous cases [16], it has been demonstrated that high diffusivity cases correspond to high values of the smaller discharge coefficient (a, which is proportional to the diffusivity), which is related to accurate solutions defined using a reduced number of LR. This must be due to the relative importance of the LR, which depends on the terms T1 and T2, producing a faster relative reduction in cases with higher discharge coefficients (see Fig. 5). As can be observed in Fig. 5, in cases with high diffusivity, using few LRs clearly has a greater influence on the final solution, while in cases with low diffusivity, T1 and T2 adopt almost the same values for all the represented LRs (all of them will be important in the final solution). This helps to understand the tremendous dispersion of the number of LRs required, depending on diffusivity.

1384

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

Linear reservoirs required to simulate natural aquifer discharges at monthly scale with various maximum error thresholds (case II)

1

Cumulative density function

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

Error <1%

Error <5%

Error <10%

Error <15%

Error <20%

0 0

10

20

30

40

50

60

70

80

90

100

110

Number of linear reservoirs Linear reservoirs required to simulate natural aquifer discharges at monthly scale with various maximum error thresholds (case II')

1

Cumulative density function

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2

Error <1%

Error <5%

Error <10%

Error <15%

0.1 Error <20%

0 0

10

20

30

40

50

60

70

80

90

100

110

Number of linear reservoirs Fig. 4. CDF of the number of LR required to simulate stream-aquifer discharge on a monthly scale with EM models for a maximum relative error of 1%, 5%, 10%, 15% and 20%. Case II and II’.

A more detailed analysis of cases I and II (cases with diffusivity higher than 105 m2/day) will be performed (see Figs. 3b, 4, and Table 3), since the interest in this methodology will grow if the number of mathematical operations required to simulate stream–aquifer flow can be significantly reduced. These high values of diffusivity are more common in confined aquifers, due to their smaller storage capacity. On the other hand, the assumption of linearity required by this methodology (linear response of head to stresses) is only strictly fulfilled in confined aquifers, although it could be applied also to unconfined aquifers when hydraulic head changes are negligible in comparison with the saturated thickness. – The most heterogeneous cases (greater K variance) show higher dispersion and higher mean and maximum values of the variable ‘‘number of LRs required’’ to simulate aquifer discharge on a monthly scale with EM models for maximum relative errors of 1% and 5%, as seen in Fig. 3b (comparing case I and I’, and case II and II’). Intuitively, we would expect that greater complexity or wider spatial variability in flow associated with highly heterogeneous conditions would oblige the use of a more complete EM model (with a higher number of LRs) to solve the problem, but this Monte-Carlo approach has allowed this idea to be demonstrated and quantified. – When the stream is partially connected, the dispersion of the distribution of the number of LRs required to simulate hetero-

geneous fields increases (as can be deduced comparing results of case I and Ip in Fig. 3b). The results of the Monte-Carlo approach seem to indicate that the effects of heterogeneity are amplified where there is partial connectivity. It seems that partial connection adds complexity to the original heterogeneous flow solutions and would force the use of more complete EM solutions (solutions with a higher number of components [LRs]) to obtain an accurate approximation. For each heterogeneous case defined with diffusivity equal to or greater than 105 m2/day, we also estimated (Table 3) the mean value of the maximum relative error obtained for each of the 100 generated fields when simulating the discharge with different EM solutions (derived with different numbers of LR: 1, 3, 7, 13, 20 and 30). The main conclusions about the sensitivity of this solution (mean maximum error of EM solutions defined with a fixed small number of LRs) to the parameters employed to define the case studies are: – Higher heterogeneity levels correspond to higher mean maximum errors. This would accord with the intuitive expectation described above, since more LRs should be required due to the increased complexity of flow modelled. – Higher mean maximum errors are associated with lower mean diffusivity. This would be due to the LRs related to smaller eigenvalues having a lower relative weight in the final solution than in cases including greater diffusivity.

1385

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389 30.00

T2( . t=month); dimensionless

T1( . t=month); dimensionless

1.20

1.00 0.80

Case I (T/S =10^7 m2/day) 0.60

Case IV (T/S =10 m2/day)

0.40

T1( , t=month) 0.20

25.00

Case I (T/S =10^7 m2/day) 20.00

Case IV (T/S =10 m2/day)

15.00 10.00

0.00

T2( , t=month) 5.00

0.00

1

2

3

4

5

6

7

8

9

10

1

2

Eigenvalue rank (from smaller to higher)

3

4

5

6

7

8

9

10

Eigenvalue rank (from smaller to higher)

T2( . t=month); dimensionless

1.00

Case I (T/S =10^7 m2/day) 0.80

0.60

0.40

T2( , t=month)

0.20

0.00 1

2

3

4

5

6

7

8

9

10

1.1

1.1

1.0

1.0

T2( . t=day); dimensionless

T1( . t=day); dimensionless

Eigenvalue rank (from smaller to higher)

0.9

T1( , t=day)

0.8

0.7 0.6 0.5 0.4

Case I (T/S =10^7 m2/day)

0.3

Case IV (T/S =10 m2/day)

0.2 0.1

0.9

T2( , t=day)

0.8

0.7

Case I (T/S =10^7 m2/day)

0.6

Case IV (T/S =10 m2/day)

0.5 0.4 0.3 0.2 0.1

0.0

0.0

1

2

3

4

5

6

7

8

9

10

1

2

Eigenvalue rank (from smaller to higher)

3

4

5

6

7

8

9

10

Eigenvalue rank (from smaller to higher)

T1( . t=day); dimensionless

0.40

Case I (T/S =10^7 m2/day)

0.35 0.30

T1( , t=day)

0.25 0.20 0.15

0.10 0.05 0.00 1

2

3

4

5

6

7

8

9

10

Eigenvalue rank (from smaller to higher) Fig. 5. T1[a(i), Dt] and T2[a(i), Dt] obtained for the 10 first linear reservoirs of one field of case I and case IV.

– The heterogeneous cases with a partial stream–aquifer connection are related to higher mean maximum errors than the same case with perfect stream–aquifer connection, which can be explained based on the hypothesis that the effects of heterogeneity are amplified when there is partial connection.

Fig. 6 depicts the monthly evolution of relative errors obtained when simulating the described scenarios with various EM models for one transmissivity field of case study II and II’. The figure shows that the maximum relative error appears at the beginning of the simulation. Although only the results for one field of cases II and

1386

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

II’ are included, the same pattern was observed in other simulated cases, and this behaviour is observed since the fastest changes in the solution appear at the beginning of the stress period (because we have assumed null initial condition [level 0 over the entire aquifer] and the most rapid modification of levels and flows will be produced in the first time step during which recharge is made). The figure also shows that higher relative errors are associated with greater heterogeneity (field of case II) with the same number of LRs. The heterogeneous fields also need more time to reach relative errors near to zero.

Table 3 Mean value of the maximum relative errors that would appear when simulating natural aquifer discharges on a monthly scale with EM solutions for various values of LR.

Number of linear reservoir = 1 Number of linear reservoir = 3 Number of linear reservoir = 7 Number of linear reservoir = 13 Number of linear reservoir = 20 Number of linear reservoir = 30

Case I

Case I’

Case Ip

Case II

Case II’

3.2 0.5 0.1 0.0

0.3 0.1 0.0 0.0

19.3 7.8 3.3 1.5

48.7 34.2 22.1 14.4

14.8 7.6 3.0 1.8

0.0

0.0

0.5

10.0

1.7

0.0

0.0

0.2

7.1

1.7

A separate question was to study sensitivity to the temporal scale, and to this end, cases II’ and II were reanalyzed on a daily timescale. Fig. 7 shows the CDF of the variable ‘‘number of LRs required’’ to simulate natural discharges on a daily timescale, maintaining a maximum relative error (Eq. (1)) below certain thresholds (1%, 5%, 10%, 15% and 20%). The number of LRs required is significantly higher when the temporal scale decreases from monthly to daily periods. This is because the relative contribution to the solution of LRs related to smaller eigenvalues will decrease as the temporal scale decreases. It is due to the exponential terms T1[a(i), Dt] = exp(a(i)  Dt) and T2[a(i), Dt] = (1  exp (a(i)  Dt))/(ai) which decrease faster when we increase the temporal scale Dt (see Fig. 5). Therefore, a higher number of LRs will be required for the final solution. The dispersion of the variable ‘‘number of LRs’’ required also shows a significant increment when the temporal scale decreases, which also seems logical when we consider the faster reduction in the terms T1 and T2. 4.2. Homogeneous case studies Lastly, we analyzed the results obtained for the homogenous cases defined in the methodology in order to identify any relationship between the number of LRs required in the heterogeneous cases and in the homogeneous ones. For each of these homogenous cases (Table 2), we obtained the number of terms required to

CASE II- Field 1: T mean = 10 m2/día; var [Ln(T)]=4.5 Variable simulated: instantaneous aquifer discharge flow rate at monthly scale

Relative error [E(%)]

55

NLR =1

NLR = 3

50

NLR = 7

NLR = 13

45

NLR = 20

40 35 30 25 20 15 10 5 0 0

2

4

6

8

10

Time (months) CASE II'- Field 1: T mean = 10 m2/día; var [Ln(T)]=0.5 Variable simulated: instantaneous aquifer discharge flow rate at monthly scale

Relative error [E(%)]

55

NLR =1

NLR = 3

50

NLR = 7

NLR = 13

45

NLR = 20

40 35 30 25 20 15 10 5 0 0

2

4

6

8

10

Time (months) Fig. 6. Time evolution of the relative errors obtained when simulating the described scenario with EM models defined with different number of linear reservoirs (NLR). Case II and II’.

1387

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389 Linear reservoirs required to simulate aquifer discharge rate at daily scale with various maximum error thresholds (CASE II)

Cumulative density function

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

Error <20%

Error <5%

Error <15%

Error <10%

Error <1%

0 0

50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800

Number of linear reservoirs Linear reservoirs required to simulate aquifer discharge rate at daily scale with various maximum error thresholds (CASE II')

Cumulative density function

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

Error <20%

Error <5%

Error <15%

Error <10%

Error <1%

0 0

50

100

150

200

250

300

350

400

450

500

Number of linear reservoirs Fig. 7. CDF of the number of LR required to simulate natural discharges on a daily timescale with EM models for different maximum relative error thresholds. Solutions for case II and II’.

Table 4 LR required to simulate aquifer discharge on a monthly scale in homogeneous cases (with 1000 active cells) with maximum relative error thresholds of 1% and 5% using EM models.

The analysis of these results allows us to deduce the following conclusions:

Aquifer diffusivity (T/S m2/day) 10^7 10^5 10^3 10 0.1 Homogeneous cases. Error <1% on a monthly scale 20,000 (‘‘perfect connection’’) 1 81 352 978 981 200 1 1 335 977 976

– The most demanding homogeneous cases (need the most number of LRs) are those defined with a perfect stream–aquifer connection, in accordance with the conclusions obtained from the analytical solution [16]. – A qualitative relationship between the diffusivity value in homogeneous cases and the mean diffusivity in heterogeneous cases could be deduced from these results. The higher diffusivity in homogeneous cases and the higher mean diffusivity in heterogeneous cases correspond to the less demanding cases (fewer LRs required, see Table 4 and Fig. 3).

simulate the described problem of natural aquifer discharges on the monthly scale with maximum errors of 5% and 1% with regards to the exact solution (see Table 4). Note that the results obtained for the number of terms required (Table 4) are higher than those obtained in a previous study [16] using a analytical solution. This is due to a different definition of the relative error. In the previous study [16], relative error was defined as a percentage of the intensity of the applied stress CVðtÞMODFLOW CVðtÞSOLUCIONðjÞ (E%ðtÞSOLUCIONðjÞ ¼  100), where R is the R recharge rate, whereas in this paper (Eq. (1)), we consider relative error with respect to an ‘‘exact solution’’.

Table 5 shows the number of terms required by homogeneous cases compared to heterogeneous problems (homogeneous cases with diffusivity equal to the mean diffusivity of the heterogeneous cases), and the probability of needing a smaller number of terms to simulate fields with the same mean diffusivity and different heterogeneity level. The first column of this table (Table 5) shows that the probability that fewer LRs are required to simulate a heterogeneous field (for the cases case I, case I’, case II, case II’, case III and case IV) than to simulate a homogeneous one (cases Hom 1, Hom2, Hom3, Hom4 with constant transmissivity equal to the mean transmissivity of the heterogeneous cases) decreases as heterogeneity increases (see Table 5). This tallies with an intuitive expectation that, as heterogeneity increases, so does the spatial complexity of the flow

Stream–aquifer connection (m2/day)

Aquifer diffusivity (T/S m2/day) 10^7

10^5

Homogeneous cases. Error <5% on a monthly scale 20,000 (‘‘perfect connection’’) 1 70 200 1 1

10^3

10

0.1

99 75

973 972

974 973

Stream–aquifer connection (m2/day)

1388

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

Table 5 Number of linear reservoirs (NLR) required by the selected homogeneous cases with maximum relative error of 1% & 5%. Probability (P) of needing fewer LR for heterogeneous cases. Homogeneous

Quantile defined by the ‘‘number of LR’’ (NLR) required in the homogeneous case (in the corresponding heterogeneous CDF)

Cases

Quantile defined by adding 15 to the ‘‘number of LR’’ (NLR) required in the homogeneous case (in the corresponding heterogeneous CDF)

Error <5%

Error <1%

Error <5%

Error <1%

Hom 1

NLR required by Hom1 (NLR_Hom1) P(NLR_Hom1 6 NLR_CASE I) P(NLR_Hom1 6 NLR_CASE I’)

1 98% 98%

1 25% 98%

15 + NLR required by Hom1 (15 + NLR_Hom1) P(15 + NLR_Hom1 6 NLR_CASE I) P(15 + NLR_Hom1 6 NLR_CASE I’)

16 100% 100%

16 100% 100%

Hom 2

NLR required by Hom2 (NLR_Hom2) P(NLR_Hom2 6 NLR_CASE II) P(NLR_Hom2 6 NLR_CASE II’)

70 97% 100%

81 90% 100%

15 + NLR required by Hom2 (15 + NLR_Hom2) P(15 + NLR_Hom2 6 NLR_CASE II) P(15 + NLR_Hom2 6 NLR_CASE II’)

85 98% 100%

96 97% 100%

Hom 3

NLR required by Hom3 (NLR_Hom3) P(NLR_Hom3 6 NLR_CASE III)

99 3%

352 3%

15 + NLR required by Hom3 (15 + NLR_Hom3) P(15 + NLR_Hom3 6 NLR_CASE III)

114 2%

367 4%

Hom 4

NLR required by Hom4 (NLR_Hom4) P(NLR_Hom4 6 NLR_CASE IV)

973 34%

978 30%

15 + NLR required by Hom4 (15 + NLR_Hom4) P(15 + NLR_Hom4 6 NLR_CASE IV)

988 78%

993 67%

and that the probability of needing a greater number of LRs increases to give an accurate simulation. Though it was not a key objective of the paper, we also tried to identify an experimental rule to use the homogenous solution with perfect GW–SW connection (the most demanding homogeneous cases) to identifying and upper bound of LRs required to simulate some heterogeneous cases. The rule is deduced experimentally (we have not found a physical explanation for it) and is based on the solutions obtained for the particular cases studies selected in this investigation and, therefore, is suggested with some reservations (differences in this rule could be observed in cases with different geometries and BC,. . .). The third column of Table 5 shows that, if the diffusivity is equal to or higher than 105 m2/day, practically all the fields generated (more than 97%) for a heterogeneous cases with perfect stream–aquifer connection, could be simulated if we use a number of LRs equal to 15 plus the number required by the homogeneous cases (see second column of Table 5: ‘‘quantile defined adding 15 to the number of LR required in the homogeneous case’’). Therefore, in these cases the homogenous solution could be employed to define an upper bound of the LR required to simulate heterogeneous cases with perfect stream–aquifer connection by adding 15 LRs to the number of terms required by the homogeneous field with the same mean diffusivity. For diffusivity values equal to or smaller than 103 m2/day a simple rule cannot be identified to use the homogenous solutions for an upper bound of the LR required for the heterogeneous cases, due to the higher diffusivity and maximum and mean errors observed in these heterogeneous cases. 4.3. Limitations and future research We conclude this section (analysis of results) by describing the limitations and possible future research. – Although the EM can be applied to simulate two-way interactions of GW–SW whenever the flow problem is linear (linear governing equation and BC) and for any kind of stress configuration (either point pumping-recharge or distributed pumpingrecharge), our analysis focuses on cases with homogenous distributed stresses. In any event, although the present study considers solely discharges under a natural regime (which frequently approximates to a uniform distribution), the linearity of the problem and proposed solution (using EM models), means that the same conclusions could be applied to modelling pumped abstractions made evenly across the aquifer (the same but with the opposite sign). Thus the conclusions are applicable to problems with two-way interactions between groundwater

and stream (The river could recharge the aquifer in cases with distributed pumping.) – The effect of heterogeneity is unknown for practical cases involving point recharge or point abstractions. Previous research has shown that, for homogeneous cases [16] where the problem is linear, these point events are modelable using EM models. We consider that this concept is sufficiently important to form the object of future investigation. – The present investigation has demonstrated that, in many heterogeneous cases, a reduced number of LRs can be used to simulate discharges from an aquifer to a river, and it has highlighted the factors that influence this. In addition, it has given a qualitative idea of the way these influences operate, based on the quantitative results obtained from the case studies designed. Nevertheless, the relationship between the number of terms required for homogeneous and heterogeneous cases has been defined in an experimental way and based only on the results of the present study. No physical explanation has so far been found to account for this phenomenon. 5. Conclusions This research analyses the relationship between the error and the number of linear reservoirs (LR) for natural aquifer discharges to a river in heterogeneous aquifers. The study identifies an upper bound on the number of LRs required to simulate heterogeneous cases, while maintaining the error under a certain threshold. A stochastic analysis of the influence of heterogeneity was performed using a Monte-Carlo approach. On a monthly scale, the case studies show that only a few reservoirs are required to accurately simulate the numerically-computed aquifer discharge in cases with high diffusivity (equal to or higher than 105 m2/day). Therefore, for these cases it is of interest to use this technique to simulate conjunctive-use alternatives for complex water resources systems at river basin scale (where there are many aquifers, reservoirs, demands, etc., which should simultaneously consider surface and groundwater flow components as well as the hydraulic interaction between them). Other conclusions derived from the sensitivity analysis performed are: – The most demanding homogeneous cases (those which most LRs) are those with perfect stream connection and low diffusivity values. This is due to a smaller relative increment of the discharge coefficients, and the influence of these on the final solution, defined by negative exponential functions, does not decrease so quickly.

D. Pulido-Velazquez et al. / Advances in Water Resources 34 (2011) 1377–1389

– The most heterogeneous cases (greatest K variance) show higher dispersion and higher mean and maximum number of LRs required, as expected, due to the higher complexity or spatial variability that forces the use of a higher number of components (LR) to give a sufficiently accurate solution. – The analysis for heterogeneous cases with low stream–aquifer connection shows higher dispersion and mean and maximum errors (the partial connection seems to amplify the effects of heterogeneity). – As can be expected, if we reduce the admissible error, the number of LRs required is greater and its dispersion decreases. – The cases defined with high mean diffusivity values (equal to or higher than 105 m2/day) allow a significant reduction in the number of LRs required to simulate with sufficient accuracy, due to a faster relative reduction of the solution’s components associated with the smaller eigenvalues (discharge coefficients). – The maximum relative errors appear at the beginning of the stress periods, as expected, since this is the period during which changes in groundwater discharges are relatively greater. – When the temporal scale is reduced to a daily scale, significantly more linear reservoirs are required, since the relative contribution of the LR related to the smaller eigenvalues on the solution decreases as the temporal scale decreases.

Acknowledgements The authors thank the reviewers for their comments, which have helped to improve significantly the clarity of the original manuscript. This research has been supported by the SAWARES project (GESHYDRO CGL2009-13238-C02-01 & HYDROECOCLIMATE CGL2009-13238-C02-02) of the Spanish Ministry of Science and Innovation (Plan Nacional I+D+I 2008-2011) and the European Community 7th Framework Project GENESIS (226536) on groundwater systems. References [1] Andreu J, Sahuquillo A. Efficient aquifer simulation in complex systems. J Water Resour Plan Manage 1987;113(1):110–29. [2] Andreu J, Capilla J, Sanchís E. AQUATOOL, a generalized decision support system for water-resources planning and management. J Hydrol 1996;177:269–91. [3] Capilla JE, Pulido-Velazquez D, Sahuquillo A, Andreu J. Solving the steady state groundwater flow equation for finite linear aquifer using a generalized Fourier series approach in two dimensional domains. Int J Numer Methods Eng 2009;79:179–204.

1389

[4] Dewandela B, Lachassagneb P, Bakalowiczc M, Weng Ph, Al-Malkid A. Evaluation of aquifer thickness by analysing recession hydrographs. Application to the Oman ophiolite hard-rock aquifer. J Hydrol 2003;274:248–69. [5] Fleckenstein JH, Krause S, Hannah DM, Boano F. Groundwater–surface water interactions: new methods and models to improve understanding of processes and dynamics. Adv Water Resour 2010. doi:10.1016/j.advwatres.2010.09.01. [6] Gómez-Hernández JJ, Journel AG. Joint simulation of multi Gaussian random variables. In: Soares A. editor. Geostatistics Tróia, vol. 92(1), Dordrecht: Kluwer; 1993. pp. 85–94. [7] Kuiper LK. Analytic solution of spatially discretized groundwater flow equations. Water Resour Res 1973;9(4):1094–7. [8] Llopis-Albert C, Capilla JE. Stochastic inverse modelling of hydraulic conductivity fields taking into account independent stochastic structures: a 3D case study. J Hydrol 2010;391:277–88. doi:10.1016/j.jhydrol.2010.07.028. [9] Llopis-Albert C, Capilla JE. Gradual conditioning of non-gaussian transmissivity fields to flow and mass transport data. 3 Application to the macrodispersion experiment (MADE-2) site, on Columbus air force base in mississippi (USA). J Hydrol 2009;371. doi:10.1016/j.jhydrol.2009.03.016. [10] Manglik A, Rai SN. Modelling of water table fluctuations in response to time varying recharge and withdrawal. Water Resour Manag 2000;14(5):339–47. [11] McDonald MG, Harbaugh AW. A modular three dimensional finite difference ground water flow model. US Geol. Surv., Washington, DC. Open – file report 83-875; 1988. [12] Moore RD. Storage-outflow modeling of stream flow recessions, with application to shallow-soil forested catchment. J Hydrol 1997;198:260–70. [13] Pulido-Velazquez D, Sahuquillo A, Andreu J. A two-step explicit solution of the Boussinesq equation for efficient simulation of unconfined aquifers in conjunctive-use models. Water Resour Res 2006;42:W05423. doi:10.1029/ 2005WR004473. [14] Pulido-Velazquez D, Sahuquillo A, Andreu J, Pulido-Velazquez M. A general methodology to simulate groundwater flow of unconfined aquifers with a reduced computational cost. J Hydrol 2007;338:42–56. doi:10.1016/ j.jhydrol.2007.02.009. [15] Pulido-Velazquez D, Sahuquillo A, Andreu J, Pulido-Velazquez M. An efficient conceptual model to simulate surface water body–aquifer interaction in conjunctive use management models. Water Resour Res 2007;43:W07407. doi:10.1029/2006WR005064. [16] Pulido-Velazquez MA, Sahuquillo-Herraiz A, Ochoa-Rivera JC, PulidoVelazquez D. Modeling of stream–aquifer interaction: the embedded multireservoir model. J Hydrol 2005;313(3–4):166–81. [17] Pulido-Velázquez M, Sahuquillo A, Andreu J. Economic optimization of conjunctive use of surface and groundwater at the basin scale. J Water Resour Plan Manage 2006;132(6):454–67. [18] Rai SN, Ramana DV, Singh RN. On the prediction of ground-water mound formation in response to transient recharge from a circular basin. Water Resour Manage 1998;12:271–84. [19] Rehfeldt KR, Boggs JM, Gelhar LW. Field study of dispersion in a heterogeneous aquifer 3. Geostatistical analysis of hydraulic conductivity. Water Resour Res 1992;28(12):3309–24. [20] Salamon P, Fernández-Garcia D, Gómez-Hernández JJ. Modeling tracer transport at the MADE site: the importance of heterogeneity. Water Resour Res 2007;43:W08404. doi:10.1029/2006WR005522. [21] Sahuquillo A. An eigenvalue numerical technique for solving unsteady groundwater continuously in time. Water Resour Res 1983;19(1):87–93. [22] Sloan WT. A physics-based function for modeling transient groundwater discharge at the watershed scale. Water Resour Res 2000;36(1):225–41. [23] Van de Griend AA, de Vries JJ, Seyhan E. Groundwater discharge form areas with a variable specific drainage resistence. J Hydrol 2002;259:203–20.