Analytical and Numerical Groundwater Flow Solutions for the FEMME-Modeling Environment

Abstract

Simple analytical and numerical solutions for confined and unconfined groundwater-surface water interaction in one and two dimensions were developed in the STRIVE package (stream river ecosystem) as part of FEMME (flexible environment for mathematically modelling the environment). Analytical and numerical solutions for interaction between one-dimensional confined and unconfined aquifers and rivers were used to study the effects of a 0.5 m sudden rise in the river water level for 24 h. Furthermore, a two-dimensional groundwater model for an unconfined aquifer was developed and coupled with a one-dimensional hydrodynamic model. This model was applied on a 1 km long reach of the Aa River, Belgium. Two different types of river water level conditions were tested. A MODFLOW model was set up for these different types of water level condition in order to compare the results with the models implemented in STRIVE. The results of the analytical solutions for confined and unconfined aquifers were in good agreement with the numerical results. The results of the two-dimensional groundwater model developed in STRIVE also showed that there is a good agreement with the MODFLOW solutions. It is concluded that the facilities of STRIVE can be used to improve the understanding of groundwater-surface water interaction and to couple the groundwater module with other modules developed for STRIVE. With these new models STRIVE proves to be a powerful example as a development and testing environment for integrated water modeling.

1. Introduction

There is a need to evaluate groundwater-surface water (GW-SW) interaction for water and ecosystem management. This is essential because linkages and feedback between groundwater and surface water systems affect both the quantity and quality of available water required by humans and ecosystems [1,2,3,4]. Therefore, the research topic of GW-SW interaction has gained importance in the last two decades, because of its role in conjunctive use, riparian zone management, and ecohydrology. In addition, understanding the interaction between groundwater and surface water can be important in the determination of migration pathways for contaminants [5]. The hydraulics of groundwater interaction with adjoining streams, canals, and drains is an important aspect of many hydrogeologic systems. Examples are the support of groundwater discharge to stream flow, bank storage attenuation of flood waves, and how groundwater discharge to streams lowers water tables maintains favorable root-zone salinity levels and prevents water logging of soil [6].

A variety of investigation methods have been used to study the hydraulic interaction of stream-aquifer systems including analytical, numerical, chemical, and field methods [7,8,9,10,11]. Recent examples of improved capabilities of MODFLOW [12,13,14] for stream-aquifer interaction are GSFLOW [15], HYDRUS [16,17,18] and unsaturated-zone flow (UZF1) packages [19], MIN3P [20], PARFLOW [21], the integrated water flow model (IWFM) [22], and SWAT-MODFLOW [23]. Numerical, chemical, and field methods have been widely used in different regions [24,25,26,27,28,29,30]. To study the interaction of groundwater and surface water flow in a river and adjacent areas, analytical solutions are often advantageous because of their simplicity. They are more general than site-specific field experiments but yet easier to implement for a particular site than numerical models [6]. In fact, several analytical solutions have been published for evaluation of the interaction of groundwater systems and hydraulically connected surface water bodies such as streams, lakes, reservoirs, drains, and canals. These solutions can be useful for understanding base flow processes, determining aquifer hydraulic properties, and predicting the response of aquifers to changing stream stage. Most of the solutions have been developed for confined and unconfined aquifers, such as [6,31,32,33,34,35,36,37].

The goal of this paper is to develop and test modules for groundwater-surface water interaction as part of the STRIVE (stream river ecosystem) package within the flexible environment for mathematically modelling the environment (FEMME) software [38]. Both numerical and analytical solutions have been developed to evaluate hydraulic interaction of river-aquifer systems. The analytical solutions from [31,32] for an unconfined aquifer and from [33] for a confined aquifer to calculate groundwater heads and discharges of the aquifer are described. The numerical solutions are based on the explicit finite difference approximation of the transient flow equation in a saturated, homogeneous, and isotropic aquifer [39]. A two-dimensional groundwater module for an unconfined aquifer is coupled to a one-dimensional hydrodynamic model [40]. This model was tested for a part of the Aa River, Belgium. Inter-model comparison is performed with MODFLOW. We present the modeling methodologies (FEMME, STRIVE package, hydrodynamic model, analytical and numerical solutions) as well as applications and comparison between analytical and numerical solutions, coupling with a hydrodynamic model, and comparison between STRIVE and MODFLOW.

2. Materials and Methods

2.1. FEMME Modeling Environment and STRIVE Package

FEMME was developed by the Netherlands Institute of Ecology (NIOO) [38]. FEMME is a modeling environment for the development and application of ecological time-dependent processes by using numerical integration of time-dependent differential equations. FEMME is constructed for ecosystem modeling and enables the simulation of different physical, biogeochemical, and transport processes of river ecosystems, like retention or exchange of matter [38,41]. The program is written in FORTRAN; it is designed to implement 0- to multi-dimensional, time-dependent models. FEMME is open source and facilitates the use of pre-defined integration tools in a modular FORTRAN environment [38]. FEMME consists of a wide range of numerical functions as integration, forcing, and calibration functions, as well as data manipulation functions. These technical possibilities allow the user to focus on the scientific part of model development rather than having to deal with complex programming issues. Hence, the environment allows rapid and efficient code development for environmental applications.

The STRIVE-package is a set of modules incorporated in the FEMME environment. It consists of different modules for macrophyte growth, water quality, and hyporheic processes, which can be coupled to form numerical models for specific research questions regarding river ecosystems. Hence, the module for hydrodynamic flow can, e.g., interact with a macrophyte growth routine, which influences the Manning coefficient and therefore the flow simulation. A heat transport module was implemented [42,43] in the STRIVE package for studying the vertical interaction in the hyporheic zone between rivers and aquifers. In this article, the STRIVE package is extended with groundwater modules, which are necessary to understand the interaction between a river and its riparian margin.

2.2. Hydrodynamic Module in STRIVE Package

A hydrodynamic surface water flow module [40] was developed for the STRIVE package and applied for a part of the Aa River, Belgium. The module is based on the Saint-Venant equations for one dimensional unsteady open channel flow [44]. Based on the capabilities of STRIVE, we hypothesize that the implementation of simple analytical and numerical groundwater flow solutions coupled with the hydrodynamic surface water module will allow the investigation of groundwater-river exchange processes in more detail.

2.3. Analytical Solutions for Groundwater-Surface Water Interaction

In this section, analytical solutions [31,32,33] for the interaction between surface water and unconfined and confined aquifers are presented.

2.3.1. Edelman Analytical Solution

We describe here the interaction between a river and a one-dimensional homogeneous semi-infinite and unconfined aquifer, which is bounded at one side by a fully-penetrating stream and below by an impermeable stratum (Figure 1). Under steady-state conditions, the water table in the aquifer and the water level in the river coincide, and there is no flow in or out of the aquifer. A sudden rise in the water level of the river induces a flow from the river towards the aquifer. As a result, the water table in the aquifer starts rising until it reaches the same level as that in the river. The flow from the river to the aquifer is unsteady and one-dimensional.

If the Dupuit approximation, i.e., (a) the equipotential lines are vertical and its consequence; (b) the slope of the groundwater table is equivalent to the hydraulic gradient and is invariant with depth, holds [45], then the flow problem can be described by the partial differential equation:

$$\frac{\partial h}{\partial t}=\frac{kb}{S_y}\frac{{\partial }^{2}h}{\partial x^2}$$

(1)

where kb = T is the transmissivity of the homogeneous aquifer [L²T⁻¹], h is the hydraulic head in the aquifer [L], t is the time [T], x is the distance from the river bank [L], and S_y is the specific yield (−). A general solution to Equation (1) does not exist and integration is possible only for specific boundary conditions [46].

Edelman [31] presented a solution for the case of a sudden change of the water level in the river and a constant water level thereafter, h(x, t) = head.

$$h(x,t)=a(1-\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^u{e}^{-{\mu }^2}}d\mu )$$

(2)

where h(x, t) is the head in the aquifer [L] at the horizontal coordinate x [L] and time t [T], a is the sudden change in the water level of the river [L], and u is a dimensionless auxiliary variable.

$$u=\sqrt{x^2{S}_y/4Kbt}$$

(3)

$$erf(u)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^u{e}^{-{\mu }^2}}d\mu$$

(4)

where erf(u) is the error function, and 1-erf(u) = erfc(u) is the complementary error function. Hence, the head in the aquifer can be formulated:

$$h(x,t)=a\cdot erfc(u)$$

(5)

The flux in the aquifer per unit length of river at distance x is:

$$q(x,t)=\frac{a}{\sqrt{\pi }}\frac{\sqrt{Kb{S}_y}}{\sqrt{t}}{e}^{-u^2}$$

(6)

Equation (6) gives the discharge from one side of the river. This equation can also be used if the water level in the river suddenly drops, inducing a flow from the aquifer to the river, resulting in a fall of the water table in the aquifer.

2.3.2. Lockington Analytical Solution

Lockington [32] derived simple analytical solutions for the one-dimensional Boussinesq equation using a weighted residual method. His approach can be applied to both a recharging and dewatering semi-infinite unconfined, homogeneous aquifer with a fully penetrating trench. Only the analytical solution for a recharging aquifer is discussed here (Figure 1) [33].

$$h=h_0+(h_1-h_0){\left(1-\frac{x}{\lambda }\sqrt{\frac{S_y}{Kt}}\right)}^{\frac{1}{\mu }}$$

(7)

where h is piezometric head [L], h₁ is the water level in the river [L], h₀ is the initial water level in the river and aquifer [L], K is the hydraulic conductivity of the aquifer [LT⁻¹], Sy is the specific yield of the aquifer [–], and λ and µ are parameters, which are defined as:

$$\mu =-\frac{3}{4}(1+N)+\frac{N}{(2-A)}+\frac{{\left[{(2-A)}^2(1+2N)+N^2{(2+A)}^2\right]}^{\frac{1}{2}}}{4(2-A)}$$

(8)

$${\lambda }^2=\frac{(1+\mu )(1+2\mu )}{2{\mu }^2}(h_0+h_1)$$

(9)

where A is a constant defined as:

$$A=\frac{4\left[h_0+(1+N)h_1\right]}{(1+N)(2+N)(h_1+h_0)}$$

(10)

In which the exponent N is estimated as in Parlange et al. [47]:

$$N=2.27932-\frac{3h_0}{(h_1+2h_0)}$$

(11)

The flux from one side of the river is given by:

$$q=\frac{C_r(h_1-h_0)\sqrt{K{S}_y}}{2\sqrt{t}}$$

(12)

where C_r is a recharge coefficient, given by:

$$C_r^2=\frac{(1+2\mu )(h_1+h_0)}{2(1+\mu )}$$

(13)

2.3.3. Bruggeman Analytical Solution

Bruggeman [33] derived the following general analytical solution:

$$h(x,t)=a\cdot 2^n\Gamma (1+\frac{n}{2})t^{\frac{n}{2}}erfc(u)$$

(14)

where:

$$u=\sqrt{x^{2}S/4Kbt}$$

(15)

S is the storage coefficient of the aquifer (−), and n = 0, 1, 2, 3… depends on the type of the water level change in the river. For n = 0, the change in the water level is assumed to be a sudden change. Similarly, an n value of 1 and 2 indicate a linear and parabolic water level change, respectively. Using the following relationship:

$$i^{n}erfc(0)=\frac{1}{2^n\Gamma (1+\frac{n}{2})},n=-1,0,1,2,\dots$$

(16)

Equation (14) can be simplified:

$$h(x,t)=a\cdot t^{\frac{n}{2}}\frac{i^{n}erfc(u)}{i^{n}erfc(0)}$$

(17)

The flux is calculated based on Darcy’s law:

$$q(x,t)=\frac{a}{2}\cdot t^{\frac{n-1}{2}}\sqrt{KbS}\frac{i^{n-1}erfc(u)}{i^{n}erfc(0)}$$

(18)

where:

$$i^{n}erfc(u)=-\frac{u}{n}{i}^{n-1}erfc(u)+\frac{1}{2n}{i}^{n-2}erfc(u),n=1,2,3,\dots$$

(19)

with:

$$i^{0}erfc(u)=erfc(u),$$

(20)

$$i^{-1}erfc(u)=\frac{2}{\sqrt{\pi }}{e}^{-u^2}$$

(21)

For n = 0, a sudden change in the surface water level, Equation (17) simplifies to:

$$h(x,t)=a\cdot erfc(u)$$

(22)

and Equation (18) for the flux from one side of the aquifer to the river reduces to:

$$q(x,t)=\frac{a}{\sqrt{\pi }}\frac{\sqrt{KbS}}{\sqrt{t}}{e}^{-u^2}$$

(23)

2.4. Numerical Solutions for Groundwater-Surface Water Interaction

Numerical solutions for transient groundwater flow (Equation (24)) are implemented in STRIVE. The solutions are based on the explicit finite difference approximation of transient flow in a saturated, homogeneous, incompressible, and isotropic aquifer [39].

$$\frac{{\partial }^{2}h}{\partial x^2}+\frac{{\partial }^{2}h}{\partial y^2}=\frac{S}{T}\frac{\partial h}{\partial t}-\frac{R(x,y,t)}{T}$$

(24)

where x and y are the spatial coordinates [L], and R is recharge [L].

In order to solve this diffusion equation (Equation (24)), it is necessary to prescribe boundary and initial conditions [39]. Finite difference approximations for different unsteady-state groundwater flow problems are developed in the following sections.

2.4.1. One-Dimensional Flow in a Confined Aquifer

The explicit finite difference approximation for the head for one-dimensional flow in a confined aquifer can be calculated from the heads at time moment n.

$$h_i^{n+1}=h_i^n+F(h_{i+1}^n-2h_i^n+h_{i-1}^n)$$

(25)

with:

$$F=T\Delta t/S{(\Delta x)}^2$$

(26)

In order to implement Equation (25) in STRIVE, the equation should be written as the rate of change in head:

$$\frac{d{h}_i}{dt}=G(h_{i+1}^{}-2h_i^{}+h_{i-1}^{})$$

(27)

where G is defined as:

$$G=T/S{(\Delta x)}^2$$

(28)

The criterion for stability of the numerical solution of Equation (25) is F < 0.25. ∆t is the time step [T], ∆x is the size of the grid cell [L], $h_i^n$ and $h_i^{n+1}$ are the heads in the center of grid cell i at respectively time step n and n + 1, $h_{i+1}^n$, $h_{i-1}^n$ are the heads at times step n in the centers of respectively grid cell i+1 and i−1. The discharge from the river to the aquifer or vice versa is calculated by Darcy’s law:

$$Q=-KA\frac{dh}{dl}$$

(29)

where Q is the discharge from or into the aquifer [L³T⁻¹], K is the hydraulic conductivity [LT⁻¹], and A is the cross-sectional area normal to the flow direction [L²].

2.4.2. One and Two-Dimensional Flow in an Unconfined Aquifer

The explicit finite difference approximation for the head at time step n+1 in terms of υ for one-dimensional flow in unconfined aquifers is [39]

$${\upsilon }_i^{n+1}={\upsilon }_i^n+F\sqrt{{\upsilon }_i^n}({\upsilon }_{i+1}^n-2{\upsilon }_i^n+{\upsilon }_{i-1}^n)$$

(30)

with:

$$\upsilon =h^2$$

(31)

where F is defined as in Equation (26) but with Sy instead of S.

In order to implement Equation (30) in STRIVE, the equation is written as a time derivative:

$$\frac{d{\upsilon }_i}{dt}=G\sqrt{{\upsilon }_i^{}}({\upsilon }_{i+1}^{}-2{\upsilon }_i^{}+{\upsilon }_{i-1}^{})$$

(32)

where G is defined as in Equation (28) but with Sy instead of S. Considering a finite set of points on a regularly spaced grid (Figure 2), the explicit finite difference approximation for unconfined two-dimensional flow is:

$${\upsilon }_{ij}^{n+1}={\upsilon }_{ij}^n+F\sqrt{{\upsilon }_{ij}^n}({\upsilon }_{i+1,j}^n+{\upsilon }_{i,j+1}^n-4{\upsilon }_{ij}^n+{\upsilon }_{i-1,j}^n+{\upsilon }_{i,j-1}^n)$$

(33)

Here, F is equal to:

$$F=K\Delta t/S_y{a}^2$$

(34)

where a = Δx = Δy.

The formulation in STRIVE is:

$$\frac{d{\upsilon }_{ij}}{dt}=\frac{K\sqrt{{\upsilon }_{ij}^{}}}{S_y{a}^2}({\upsilon }_{i+1,j}^{}+{\upsilon }_{i,j+1}^{}-4{\upsilon }_{ij}^{}+{\upsilon }_{i-1,j}^{}+{\upsilon }_{i,j-1}^{})$$

(35)

3. Application and Discussion

3.1. One-Dimensional Analytical and Numerical Solutions for Confined and Unconfined Aquifers

The implemented one-dimensional analytical and numerical solutions for groundwater heads and discharges in unconfined or confined aquifers are compared as a function of time and distance. The solutions were set up in the STRIVE package for a domain perpendicular to the river to calculate the rise in the head and flow in the unconfined or confined aquifer at a distance x from the river after a sudden water level rise of 0.5 m.

Table 1. Parameters and dimensions of the unconfined and confined aquifer.

Parameter	Value	Dimension
Hydraulic conductivity (K)	10	m d−1
Storage coefficient (S)	0.2	(−)
Thickness of the aquifer (b)	10	m
Specific yield (Sy )	0.2	(−)
Initial head everywhere in the aquifer	10.4	m
Head in the river	10.9	m

specifies the details for the unconfined and confined aquifer.

The heads are calculated at the center of cells, and the discharges are calculated using Darcy’s equation at the interface of the cells based on the head calculations for a one-day simulation with a time step of 0.001 day. The results of this application are presented below as a comparison between the analytical and numerical solutions.

Figure 3 shows the numerical solution for one-dimensional transient flow through an unconfined aquifer and the analytical solutions of [31,32] for the heads at different lateral distances from the river at 1.5, 12, and 24 h The differences between the solutions of [31,32] and the numerical one are negligible; the maximum difference is 0.021 m between the Edelman [31] and the numerical solutions, and the root mean squared error (RMSE) is 0.0075 m after 1.5 h. The maximum difference between the Lockington [32] and the numerical solutions was 0.015 m, and the RMSE was 0.0048 m.

Figure 4 shows the match between the three solutions for the discharges at the river-aquifer boundary for a simulation period of 1 day. The maximum difference between the analytical (Edelman [31] and Lockington [32]) and the numerical solutions for the groundwater flux per meter of river length is 0.17 m² d⁻¹, while the RMSE is 0.1 m² d⁻¹ at 1.5 h. The difference decreases with time.

The groundwater heads and fluxes simulated with the numerical solution for one-dimensional transient flow in a confined aquifer and the analytical solution of Bruggeman [35] are presented in Figure 5 and Figure 6. Figure 5 shows the relationship between groundwater heads as a function of distance from the river at different times (1.5, 12, 24 h). It is observed that the effect of a 0.5 m sudden rise in the river stage is negligible beyond a distance of 80 m from the river-aquifer boundary. The distance increases with time; at 1.5 h the effect disappears after 50 m and at 24 h it is negligible after 80 m from the river-aquifer boundary. The maximum change in head in the aquifer is attained after one day for locations further than 10 m away from the river. The maximum difference between the numerical and Bruggeman [33] solutions was 0.01 m at 1.5 h.

Figure 6 shows the groundwater flux as a function of distance from the river at different times (1.5, 12, 24 h). It is noticed that the maximum discharges are attained within the first 1.5 h and range from 4.13 to 1.25 m² d⁻¹ per meter of river length for a location 5 m away from the river-aquifer boundary. The maximum difference between the numerical and Bruggeman solutions was 0.3 m² d⁻¹ and the RMSE was 0.0145 m at 1.5 h. The accuracy of the numerical solution is mainly determined by the space and time discretization. It is concluded that the results of the analytical solutions of Edelman [31], Lockington [32], and Bruggeman [33] are in a good agreement with the numerical solution.

3.2. Two-Dimensional Numerical Solution in an Unconfined Aquifer

The two-dimensional groundwater module for an unconfined aquifer is applied to the Aa River aquifer, Belgium (based on data from [42,43]) (Figure 7).

The simulated area is simplified to be as Figure 8; it has a length of 1000 m (between weirs 3 and 4 in Figure 7), a width of 200 m, a thickness of 20 m, a hydraulic conductivity K_x and K_y of 10 m d⁻¹, a specific yield S_y of 0.2, and a grid cell size of 10 m. To simplify implementation of the problem in STRIVE, we simulated the reach of 1000 m of the Aa River as a straight line, as shown in Figure 8. The groundwater module developed in STRIVE was applied for two cases of the river boundary conditions. In case 1, the river boundary is estimated by interpolation based on the two head values assigned to the upstream and downstream points of the river. The model was run in steady state by using boundary conditions based on interpolation of heads assigned at the corner cells (a = 10.5 m, b = 10.2 m, c = 10.5 m, d = 10.8 m) (Figure 8). The boundary condition heads are interpolated linearly based on the heads assigned at the corner cells. The initial head everywhere in the aquifer is specified as 10.4 m. A time step of 10 seconds was used for a duration of 1000 days.

The simulated groundwater head and the lateral flows are presented in Figure 9 and Figure 10 respectively. From Figure 9, it is observed that the hydraulic gradient is directed towards the river boundary, due to the selected boundary conditions.

Figure 10 shows in two dimensions the lateral groundwater fluxes obtained from the STRIVE model. The fluxes are directed towards the river (values are negative); the maximum and minimum discharges from the aquifer to the river along the river are respectively 1.49 and 0.74 m³ d⁻¹ at 900 m and 200 m viewed from the upstream boundary respectively.

In the second case, the river boundary is based on water levels from a hydrodynamic surface water model, which is also integrated in STRIVE [40]. The water level in the hydrodynamic surface water model was calculated based on an upstream discharge condition, formulated as a half sine wave (Figure 11). Other boundary conditions and initial conditions are the same as those applied in case 1. The model was run in transient for a period of 30 days with a time step of 1 min and output interval of 1 h. The river boundary values in this case were obtained from the output of the hydrodynamic model.

Figure 12 presents the two-dimensional groundwater heads as simulated with STRIVE for this second case. The groundwater heads are presented at two times, in order to study the effect of the river stages on the groundwater heads in the aquifer and to test the response of the interaction between the river and the aquifer. Figure 12a shows the groundwater heads in the aquifer at 24 h (the beginning of the pulse of the upstream discharge in the river). It is observed that the hydraulic gradient is directed towards the river boundary. Figure 12b presents the groundwater heads in the aquifer at 32 h (the maximum effect of the pulse of the upstream discharge in the river). It is observed that at the river-aquifer boundary, the hydraulic gradient is directed from the river boundary to the aquifer.

In Figure 13, the effect of the river pulse on groundwater heads at different times (24, 32, and 480 h) and at different distances along the river at 100 m (near to the upstream model boundary) are shown. At 24 and 480 h (before the beginning and after the ending of the pulse), there is no effect of the river pulse on the groundwater heads, and the flow is directed towards the river boundary. However, at 32 h, the effect of the river pulse on the groundwater heads in the model appears clearly at the edge of the river-aquifer boundary (30 m) and disappears after a distance of 30 m from the river-aquifer boundary.

The river-aquifer fluxes along the river at 24, 32, and 480 h are shown in Figure 14. At the beginning of the pulse, the Aa River gains water from the aquifer, while at 32 h, the aquifer gains water from the Aa River. The Aa River comes back again to gain water from the aquifer after 480 h. The direction and the rate of flow depend on the difference between the river stage and the groundwater heads.

3.3. Comparison between STRIVE and MODFLOW Results

The groundwater flow model MODFLOW-2000 [12] was used for model comparison of the STRIVE two-dimensional groundwater flow implementation. A simple MODFLOW model was set up using the same data as was used for the second case example in STRIVE. The PCG2 solver was used and both steady-state and transient state simulations were performed in order to compare these results with STRIVE results.

In order to see the agreement between the STRIVE and MODFLOW results, the groundwater heads in lateral direction from the river-aquifer boundary at different distances along the river at 5 m (upstream river), 505 m (the middle of the river), and 995 m (downstream river) are shown in Figure 15. Figure 16 shows the groundwater heads in lateral direction from the river-aquifer boundary at distance of 505 m along the river at different simulation times (1, 6, 12, 24 h). The maximum head difference between STRIVE and MODFLOW was 0.004 m, and the RMSE is 7.1 × 10⁻⁶ m and 1.8 × 10⁻⁵ m for steady-state and transient state simulations, respectively.

4. Conclusions

The facilities of STRIVE were tested for the interaction between groundwater flow in confined and unconfined aquifers with streams in one and two dimensions. The analytical solutions implemented in STRIVE include Edelman [31] and Lockington [32] for unconfined aquifers and Bruggeman [33] for confined aquifers. Groundwater heads and discharges in the aquifers were calculated based on a sudden change in the river stage. The results of the analytical solutions were compared with one-dimensional numerical solutions which were implemented in STRIVE for confined and unconfined aquifers. The results of the analytical solutions for confined and unconfined aquifers were in good agreement with the numerical results.

In addition, a two-dimensional groundwater model for an unconfined aquifer was developed in STRIVE and coupled with a one-dimensional hydrodynamic model in STRIVE. This model was applied on a 1 km long reach of the Aa River. The model was tested for two cases with different boundary conditions. In the first case, the river boundary was interpolated, and the model was simulated in steady-state. In the second case, the river boundary was linked with the water levels from a hydrodynamic surface water model. MODFLOW models were set up for these cases as well, in order to check the implementation in STRIVE. The results of the two-dimensional groundwater model developed in STRIVE showed that there is a very good agreement with MODFLOW.

It is concluded that analytical and numerical solutions for groundwater-surface water interaction for unconfined and confined aquifers have been successfully implemented in STRIVE. Hence, STRIVE is extended in terms of modeling facilities for groundwater-surface water interaction, but also due to these implemented sub-models, new integration possibilities with existing modules such as the hydrodynamic, hyporheic zone and sediment appear. The flexibility of STRIVE has proven to be a major advantage in developing these new sub-modules. With these new models, STRIVE increases its capabilities without becoming a dedicated type of model with a graphical user interface. Rather, it is a powerful example of a development and testing environment for integrated water modeling.