## Abstract

Protecting groundwater resources plays an important role in watershed management. For this purpose, studies on groundwater flow dynamics incorporating surface water–groundwater interactions have been conducted including analytical, numerical, and experimental models. In this research, a stream–aquifer system was considered to understand the physical behavior of surface water–groundwater interactions. Interactions in a stream–aquifer system were incorporated into the mathematical modeling by defining the stream head as a boundary condition for the groundwater flow equation. This boundary was chosen as a sloping stream boundary, which is an approach in representing the natural conditions of the stream and may be used to define continuous interactions between stream and aquifer. A semi-analytical solution for transient 2D groundwater flow was developed for the considered problem. Isotropic, homogeneous, and finite aquifer assumptions were made in order to define the aquifer characteristics. Then, a series of laboratory experiments was conducted to simulate this stream–aquifer system. Finally, a numerical model was developed by using Visual MODFLOW to verify analytical and experimental results. Numerical results matched with both analytical solutions and the experimental observations.

## INTRODUCTION

Protecting and managing groundwater resources play an important role in watershed management. Therefore, the groundwater flow dynamics needs to be investigated and understood in detail. A great deal of analytical, numerical, and experimental studies are conducted on solutions of groundwater flow equation. The analytical solutions of this equation involve some simplifications including mostly linearization of the partial differential equation. The hydrologic features of the chosen system and the boundary conditions are other factors which are effective on the solution methods (Spanoudaki *et al.* 2010). One of the major assumptions is made on the dimension of the flow. Solutions of a 1D equation are commonly developed among the models which exist for both confined and unconfined aquifers. For a confined and homogeneous aquifer, solutions of a 1D groundwater flow equation have been developed by Hantush (1961), Ferris (1963), Cooper & Rarabough (1963), Hall & Moench (1972), and Oakes & Wilkinson (1972). The most comprehensive 1D solution was done by Singh (2003, 2004) for a confined aquifer. The solution was done for semi-pervious stream banks with partial penetration in an infinite, homogeneous, and isotropic aquifer. In addition, four different types of boundary conditions were applied to the equation: sinusoidal, linear, unit step change, and symmetrical flood waves. Also, Townley (1995) developed a 1D solution of groundwater flow equation by defining a time dependency for the source/sink term.

For an unconfined aquifer, the 1D flow equation is known as the Boussinesq equation. The linearization of this equation depends on the assumptions of the system. Several solutions have been developed under different assumptions, such as fully or partially penetrating stream with or without pervious banks; homogeneous, isotropic, and finite or infinite aquifer (Workman *et al.* 1997; Serrano & Workman 1998; Ostfeld *et al.* 1999). The most comprehensive solution for an unconfined aquifer was developed by Serrano & Workman (1998). 2D problems have also been considered with similar specifications as in 1D problems by Higgins (1980), Neuman (1981), Van de Giesen *et al.* (1994), Moench & Barlow (2000), and Barlow *et al.* (2000). Boyraz & Kazezyılmaz-Alhan (2014) developed a solution under a sloping stream boundary in a 2D horizontal domain. Serrano (2003) solved the 3D flow equation by decomposition method in a finite aquifer which is bounded by two different fully penetrating streams. The water tables of streams have been defined by analytical functions. Elastic storage was neglected and the nonlinear kinematic free surface boundary was used. As well as these approaches, which determine the stream–aquifer interactions by assuming the stream level as a boundary condition, other approaches are still valid that integrate the surface water flow equations with the groundwater flow equation. The integrated hydrodynamic solutions of the surface water and the groundwater systems have been accomplished by Hunt (1990), Hantush *et al.* (2002), Hantush (2005), and Lal (2001). The most comprehensive solution was developed by Lal (2001), who presented a solution of the coupled linearized 2D Boussinesq equation with the 1D diffusion equation.

The laboratory studies on stream–groundwater interactions have mostly been conducted for the case where the interactions occur in the x-z plane. Thus, the system involves a stream along x-direction lying on a vertically oriented aquifer. Open channel flume and tanks are used to be able to model these interactions (Rowe 1960; Elliott & Brooks 1997; Packman *et al.* 2004; Salehin *et al.* 2004; Schillig *et al.* 2014). Different stream conditions and aquifer properties have been considered and different approaches of solution methods employed. Rivière *et al.* (2014) determined the stream infiltration rate and capacity in a homogeneous sandy aquifer by using a fiberglass tank. By considering a disconnected case of stream–aquifer interactions, it is concluded that the infiltration rate of the interactions is not affected by the hydraulic head variance between the river and groundwater. Marion *et al.* (2008) used an open channel flume, which was filled with gravel and sand mix, to understand the streambed interface dynamics for different sediment beds, multi-layer aquifer, and different sloped streams. They found that the interaction has an important role on accumulation of dissolved solid material. Consequently, many studies have investigated the behavior of interactions in a vertical direction and its effects on contaminant transport (Saha 2009; Jazaei *et al.* 2014). Moreover, different shapes of stream bed have been addressed as the drivers of interactions and dye injected models applied to determine the interaction flow direction (Tonina & Buffington 2009). Interactions have also been investigated in field studies with different approaches (Loheide & Gorelick 2006; Novakowski *et al.* 2006; Gleeson *et al.* 2014; Batlle-Aguilar *et al.* 2015; Cai *et al.* 2016; Becher Quinodoz *et al.* 2017; Nakamura *et al.* 2017). These studies emphasized the crucial role of the interactions in a watershed and the impact of the interactions on water resources. Despite all of these studies, the interaction flow effects on groundwater flow and head distribution in a horizontal plane of a stream–aquifer system need further investigations as studies on this kind of a system are rare in the literature.

In order to contribute to the aforementioned studies, the objectives of this study were: (i) to determine the groundwater head distribution analytically in a stream–aquifer region where a sloped water table is defined as stream boundary; and (ii) to show the applicability of the analytical solution by using a laboratory-scale model. A 2D semi-analytical solution was derived under the assumption of isotropic, homogeneous, and finite aquifer. Then, an experimental setup was developed to simulate the groundwater flow in a stream–aquifer region with sloping stream boundary. Sand with different particle size distributions were used in the experiments and the observations were done for both steady and transient flow cases. Finally, the semi-analytical solution and experimental data were verified with a finite difference numerical method via employing the Visual MODFLOW 2010.

## METHODOLOGY

*h*is the groundwater head (L),

*K*is the hydraulic conductivity of the aquifer (L/T),

*S*is the specific storage (1/L), and

_{s}*W*is the source/sink term (1/T).

### Analytical solution

*h*is small with respect to magnitude of

*h*, Equation (1) reduces to the following: The problem was considered in a finite domain. The stream was located along the y-axis of the aquifer. The streambed was assumed as fully saturated (Figure 1). The following boundary conditions were applied to the problem: where

*J*is the streambed slope (

*L/L*),

*H*

_{0}is the initial head of the aquifer (

*L*), and

*L*and

_{x}*L*are the lengths of the aquifer along

_{y}*x*and

*y*directions (

*L*), respectively. The boundaries and the stream–aquifer system are shown in Figure 1.

*p*is the Laplace variable. Second, Fourier cosine transform was applied to Equation (8). All boundary conditions were also transformed into the Fourier cosine domain. For a finite domain, the Fourier cosine transform is given as follows: After applying Equation (9) to Equation (8), the following is obtained:

*et al.*1982) is applied to Equation (12) in order to take the inverse Laplace by a numerical inversion method, as shown below:

*h(x,y,t)*is the groundwater head solution in time domain which is obtained by taking the inverse Laplace integrals (Equation (13)).

*h*(

_{1}*x,y,t*),

*h*(

_{2}*x,y,t*), and

*h*(

_{3}*x,y,t*) are the numerical forms of inverse Laplace algorithm (Equation (14)) which are obtained by applying the trapezoidal rule to inverse Laplace integral where the real and imaginary parts of Laplace parameter

*p*are separated as .

### Experimental setup

An experimental setup was developed to observe the groundwater flow in the stream–aquifer system with sloping stream boundary. The same initial and boundary conditions of the analytical problem were created within the experimental setup (Figure 2(a)–2(c)). A plexiglass box with dimensions of 50 cm × 50 cm × 10 cm was used for the aquifer region. The box was placed on a movable table which could be adjusted to different slope values. Inflow and outflow pools were placed at the inlet and outlet points of the stream to adjust the stream head. On one side of the aquifer, a sloping stream boundary was employed by placing a stream with a length of 50 cm in the y-direction. The width of the stream channel was 1.5 cm. To adjust the streambed slope, a wooden stick (50 cm × 1.5 cm × 0.4 cm) was used as shown in Figure 2(b). The rest of the plexiglass box was filled with sand. A special filter was used to avoid sedimentation. Several wells were located in the aquifer to measure the groundwater head. Each well was of 1 cm diameter. The groundwater level was measured by using thin wooden sticks manually (Figure 2(c)).

### Numerical model

In order to verify the analytical solution and experimental results, a numerical model was developed using Visual MODFLOW. Visual MODFLOW, which is widely used to simulate the groundwater flow and transport, is used for the numerical solutions of groundwater problems. The groundwater flow equation is solved by MODFLOW (1988, 2000, 2005) package (Hill *et al.* 2000) of Visual MODFLOW using the finite difference method. Moreover, several packages are integrated into the Visual MODFLOW to model different systems which interact with the aquifer, such as wells, streams, and wetlands. In addition, Visual MODFLOW also has the capability of simulating contaminant transport processes in the aquifer.

In this study, the constant head boundary (CHB) package of Visual MODFLOW was employed to simulate the sloping stream head boundary. This boundary type fixes the hydraulic head value in the selected region during the selected analysis period (stress periods). It acts as an infinite source of water to represent the head. Linear interpolation of the specific heads in time between the beginning and end of each stress period was carried out by this package. No-flux boundaries were automatically defined to the border cells if no other boundary condition was introduced. Aquifer type and properties should be defined to complete the model. There are several numerical engines incorporated into the Visual MODFLOW and each of them needs an initial value to start for marching in time and obtaining the numerical solutions. Preconditioned conjugate gradient (PCG) solver was used for the solution of the stream–aquifer system discussed in this study.

## RESULTS

Three synthetic examples were solved using the analytical solution presented in the previous section. In the first example, the analytical solution results were compared with the Visual MODFLOW numerical model results. In the second and third examples, the experimental results were compared with the Visual MODFLOW numerical model results for steady state and transient cases, respectively. In addition, the analytical solutions for the case of same geometry and properties of the experimental system were also obtained and included in the comparison graphs.

### Example 1

In this example, an aquifer domain of 2,000 m × 2,000 m and a stream with a streambed slope of 0.0015 m/m at the boundary were selected. The initial depth of groundwater was chosen as *H*_{0} = 6 m. Figure 3(a)–3(c) show the transient and Figure 3(d) the steady state groundwater head distributions obtained by both analytical and numerical methods. As seen from this figure, the analytical solution results were very close to the Visual MODFLOW results for both transient and steady state cases. Very small differences between the two solution methods were observed due to the neglected source/sink term *W* in the groundwater flow equation when developing the analytical solution, whereas MODFLOW takes this term into account in the numerical solution. A symmetrical groundwater head distribution with respect to the x-axis passing through y = 1,000 m was observed in the case of steady state solution (Figure 3(d)). Consequently, the comparison of the models showed that the analytical solution gives accurate results.

### Example 2

In this example, an aquifer domain of 50 cm × 50 cm and a stream with a streambed slope of 0.028 m/m at the boundary were selected. Coarse sand was used as the soil type of the aquifer, where the grain size changed between 0.5 mm and 1 mm. This type of sand allows for completion of the experiment in a short period of time and thus is suitable for steady state observations. The constant head permeability test was conducted to determine the hydraulic conductivity of the sand and was found as 5 × 10^{−6} m/s.

In the experimental setup, the constant hydraulic head changed from 8.3 cm at upstream to 6.9 cm at downstream. After 2 minutes, the system reached the steady state conditions. The measurements were done after 10 minutes of the start of the experiment to make sure that the steady state condition was achieved in the stream–aquifer system. The measured groundwater level data and the corresponding analytical solution are shown in Figure 4. The groundwater level was measured at 18 points. The comparison shows that the analytical solution results are mostly compatible with the measurements. In addition, the numerical solution of the example is also presented in the same figure. As can be seen, both analytical solution and experimental results matched well with the numerical solution.

### Example 3

In this example, the transient groundwater head distribution was examined with the experimental setup. The hydraulic head slope was 0.025 and the stream water depth decreased from 8.8 cm at upstream to 7.55 cm at downstream. Fine sand was used as the soil type of the aquifer. The hydraulic conductivity of the fine sand was measured as 5 × 10^{−7} m/s with the constant hydraulic head test. The number of the observation wells was reduced to four to be able to measure the groundwater level in short time intervals. Figure 5 shows the groundwater head distribution after the first minute the experiment started. Experimental observations matched relatively well with the analytical and numerical solutions. In fact, many trials have been done to catch better agreements between observed and calculated groundwater head values. During the repetitive experiments, the difficulty of observation for the transient case was noticed. This situation was attributed to the scale of the experimental setup and isotropic and homogeneous aquifer assumptions.

## DISCUSSION

The problem considered in this study showed the stream–aquifer relation in a system where the groundwater was fed only by the stream. In addition, this system allowed for observation of the effects of the sloped stream head on groundwater head distribution. The considered problem in this study had two important assumptions on boundaries. The first assumption was the definition of a sloping stream boundary which is a good approximation of a natural stream–aquifer system. This type of boundary can also be used in interaction modeling, e.g., flooding. Incorporating the sloping stream as a boundary condition helps in simulating the real conditions of the stream–aquifer system better than a non-sloped stream boundary. Considering variable hydraulic head in stream–aquifer systems led us to improve the existing solutions of groundwater problems. Especially, the effects of surface water–groundwater interaction on groundwater head distribution differ with this boundary condition and it is worth investigating this issue further in future studies. The second assumption was the no-flow boundary condition which helps in modeling the stream–aquifer interaction analytically in groundwater flow equation. No-flow boundary may also be a good representation for aquifers lying in a large domain. The analytical solution in this study has been developed under the assumption of a neglected source/sink term. The effect of the interaction flow rate as a source/sink term on groundwater head distribution is also worth investigating in groundwater flow studies.

In order to observe the hydraulic behavior of the stream–aquifer system considered in this study, a series of experiments were conducted. It is a reduced scale physical model of the problem domain. For the steady state condition, the results of the analytical solution and the experimental measurements matched quite well. However, the experimental results for the transient condition were not as good as the ones for the steady state case. This situation may be attributed to the complexity of preparing a homogeneous and isotropic aquifer. Unless the aquifer is prepared exactly the same as the one defined in the analytical problem, the observed groundwater head values at the time steps presenting in the unsteady case do not match well with the calculations. In addition, employing a small scale physical model increases this said effect on the observations.

The Visual MODFLOW numerical solution performs well for a wide and long aquifer. However, as the dimensions of the model decrease, the accuracy of the numerical model also decreases. Time steps are also important for conditionally stable numerical methods. In Visual MODFLOW, several numerical methods are available which can be chosen according to the physical problem with various scales. For a small scale model such as the one used in the experimental part of this study, PCG solver was determined as the best numerical method after testing different numerical methods. This solver solves the groundwater flow equation by using residual and head change criteria in inner and outer iterations. The numerical results matched both with measurements and analytical solution results within an acceptable range.

## CONCLUSIONS

In this study, a semi-analytical solution was developed for a 2D groundwater flow equation in an isotropic and homogeneous aquifer. The selected boundaries allowed for simulation of the stream–aquifer interactions in the system. Then, an experimental setup was developed to have a better understanding of the hydraulic behavior of the stream–aquifer system. Finally, a numerical model was developed using Visual MODFLOW to evaluate and verify both analytical and experimental results. The analytical results matched very well with the numerical results for both steady and unsteady state cases. The experimental observations performed also quite well for the steady state case, whereas there is room for improvement in unsteady state experimental results. All comparisons showed that the analytical solution worked well and provides a physical insight into the groundwater problems which involve stream–aquifer interactions. The presented solution not only contributes to the analytical and numerical studies of groundwater but is also useful in field studies and forecasting groundwater hydraulics.

## ACKNOWLEDGEMENTS

The first author would like to thank the Scientific and Technological Research Council of Turkey (TUBITAK) for the financial support of his PhD studies. Also, the authors would like to express their gratitude to the anonymous reviewers and the Editor for their excellent suggestions, which strengthened the paper.

## REFERENCES

*MS Thesis*