<figure><img src="image_1.png" /></figure>

# Steady-state semi-analytical solutions for assessing the two-dimensional groundwater hydraulic head depletion induced by river dam removal

Benoît Dewandel<sup>a,b,*</sup>, Sandra Lanini<sup>a,b</sup>, Nicolas Frissant<sup>a,b</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Patrick Lachassagne, Associate Editor

**Keywords:**

Dam removal

Impact on groundwater level

Analytical solutions

Fundamental solutions method

Groundwater hydraulics

## ABSTRACT

Dams are pervasive features of the world's river systems. The reservoir associated with the dam modifies the distribution of the hydraulic head in the aquifer and the natural groundwater flow. However, such modifications must be defined to forecast their environmental, economic and/or social impacts. Based on the method of fundamental solutions (MFS), steady-state semi-analytical solutions are proposed for evaluating the long-term spatial distribution of the rise or decline of the hydraulic head in an aquifer caused by a dam reservoir on a watercourse or its removal. MFS was chosen because of its ability to account for different boundary conditions along the banks of the watercourse. The dam reservoir can have a triangular or rectilinear geometry. The aquifer is assumed unconfined, homogeneous, inclined and finite. The solutions take into account both the reservoir-aquifer interaction and those between the stream and the aquifer up- and downstream of the reservoir. The reservoir-stream system fully penetrates the aquifer, and is separated from the aquifer by semi-pervious banks (Robin condition). Analysis of the semi-analytical solution included sensitivity tests, and assessment of the influence of several parameters: reservoir leakage parameter, width of the reservoir and width of the aquifer. The results show that the reservoir and river leakage parameters up- and downstream of the reservoir, as well as the geometry of the reservoir dam, are the main determining parameters of the extent of the induced groundwater mound. They show that low aquifer-river exchanges in the river upstream and downstream of the reservoir increases the spatial extension of the mound. The degree of error introduced by the method of linearization used to solve the governing partial differential equation is also discussed. A solution has been also developed with a pumping well, particularly to assess what the depletion of the water level at the well location might be if the dam is removed. The proposed solutions were applied at two field sites, giving satisfactory results.

These semi-analytical solutions will be useful applications for assessing the long-term spatial impact on the aquifer of the emplacement, removal or leveling of a reservoir dam.

## 1. Introduction

Dams are pervasive features of the world's river systems. The most recent edition of the World Register of Dams of ICOLD (2023) contains information on 51,325 large dams (higher than 15 m). Worldwide references on smaller dams are scarce, U.S. (USACE, 2008) referenced 2.5 million small dams (less than 1.8 m high) and 80,000 medium dams (over 1.8 m high), in Europe (without Russia), the European Rivers Network estimates the number of smaller dams (less than 15 m height) to hundreds of thousands. Worldwide, Lehner et al. (2011) estimate that smaller dams (<15 m) may be more than 16 million.

Dams provide flood control, hydroelectric power, reliable reservoirs

for surface water supply, recreational activities, etc., but nowadays there is an increasing interest for dam removal. Such interest is motivated by adverse ecological and social impacts, safety conditions associated with aging dams to prevent dam failure, and appreciation for societal values linked to healthy rivers and fisheries (e.g. World Commission on Dams, 2000; Pejchar and Warner, 2001; Bednarek, 2001; Johnson and Graber, 2002; Collins et al., 2007; Graf, 2005). Thus, removal of dams is often considered as a mean of restoring natural streamflow and sediment transport, resulting in improved riparian corridors, fishery habitats, sports fishing, recreational rafting, and reappearance of a more pastoral and riparian vegetation setting. However, dam removal may have short-term negative ecological impacts due to the increase of sediment load, which may cause suffocation and abrasion to various biota and habitats,

<footnote><sup>a</sup> BRGM, Univ. of Montpellier, Montpellier, France</footnote>

<footnote><sup>b</sup> G-eau, UMR 183, INRAE, CIRAD, IRD, AgroParisTech, Institut Agro Montpellier, BRGM, Montpellier, France</footnote>

<footnote>* Corresponding author at: BRGM, Univ. of Montpellier, Montpellier, France. E-mail address: b.dewandel@brgm.fr (B. Dewandel).</footnote>
List of symbols

a water depth in the upstream part of the reservoir at $x = -R$, in m

b water depth downstream of the reservoir at the dam location (at $x = R$), in m

$b_{dw}'$ thickness of the semi-pervious layer between the stream and the aquifer downstream of the reservoir, in m

$b_R'$ thickness of the semi-pervious layer between the reservoir and the aquifer, in m

$b_{up}'$ thickness of the semi-pervious layer between the stream and the aquifer upstream of the reservoir, in m

$c\_av$ water depth in the stream, in m

$D$ saturated thickness of the aquifer, in m

$D_0$ saturated thickness of the aquifer without watercourse (stream or reservoir), in m

$d'$ distance of image wells to the stream and reservoir banks, in m

$e^x$ exponential function

$h$ hydraulic head, in m

$h_0$ the hydraulic head in the watercourse, in m

$i$ hydraulic gradient, $i = \tan\theta$, dimensionless

$K$ hydraulic conductivity of the aquifer, in m/s

$K_0(u)$ modified Bessel function of the second kind of the zero order

$K_1(u)$ modified Bessel function of the second kind of the first order

$k_{dw}'$ hydraulic conductivity of the semi-pervious layer between the stream and the aquifer dowstream of the reservoir, in m/s

$k_R'$ hydraulic conductivity of the semi-pervious layer between the reservoir and the aquifer, in m/s

$k_{up}'$ hydraulic conductivity of the semi-pervious layer between

$L$ the stream and the aquifer upstream of the reservoir, in m/s width of the aquifer, in m

$L_R$ length of the reservoir along the y-axis, in m

$n$ number of image wells

$P$ distance, $P = 20 \cdot R$, in m

$Q_p$ pumping (or injecting) rate, in m³/s

$r$ radius, $r = \sqrt{x^2 + y^2}$, in m

$R$ half-length of the reservoir along the x-axis, in m

$R_{j,1}$, $R_{j,2}$ distances to the image wells, in m

$R_r$ reservoir leakance parameter (retardation coefficient), $R_r = \frac{K}{k_R}b_{R'}$, in m

$R_{r\_up}$ stream leakance parameter (retardation coefficient) upstream of the reservoir, $R_{r\_up} = \frac{K}{k_{r\_up}}b_{r\_up}'$, in m

$R_{r\_dw}$ reservoir-stream leakance parameter (retardation coefficient) downstream of the reservoir, $R_{r\_dw} = \frac{K}{k_{r\_dw}}b_{r\_dw}'$, in m

$s_0$ water height in the watercourse (stream or reservoir), in m

$s$ variation of the aquifer saturated thickness, $s = D-D_0$, in m

$x, y$ coordinates of a Cartesian system, in m

$x_j, y_j$ coordinates of image wells, in m

$x_w, y_w$ coordinates of the well pumping (or injecting) well, in m

$\alpha$ parameter, $\alpha = i/(2D_0)$

$\beta$ degree of saturation, dimensionless

$\gamma_j$ strength coefficient of the jth image well, m³/s

$\tilde{\gamma}_j$ strength coefficient of the jth image well, m

$\varphi$ angle between the normal to the limit watercourse-aquifer and the x-axis, in radian

$\theta$ dip angle of the impermeable basement, in radian

$\sigma$ distance between image wells, in m

$\vartheta$ angle of the diagonal of the reservoir makes with the x-axis ($\tan\theta = L_R/2R$), in radian

$\Delta\epsilon$ incremental spacing, in m

the release of contaminated sediments or the potential generation of downstream flooding (Bednarek, 2001; Roberts et al., 2006).

Although predicting impacts of the implementation or removal of dams has historically attracted the interest of researchers for several decades in various fields: i) flow regime of rivers, ii) ecosystems, iii) sediment transport, iv) geomorphic processes, v) nutrient dynamics, vi) riparian vegetation, vii) economic and viii) social impacts (e.g. Thomas, 1976; Baxter, 1977; Williams and Wolman, 1984; Shuman, 1995; Collier and Webb, 1996; Beyer, 2005; Hart and Poff, 2002; Poff and Hart, 2002; Pizzuto, 2002; Stanley and Doyle, 2002; Shafroth et al., 2002; Whitelaw and MacMullan, 2002; Johnson and Graber, 2002; Wyrick et al., 2009; Capart, 2013), little attention has been paid to their impact on groundwater (Farinacci, 2009; Berthelote, 2013; Servière, 2021; Li et al., 2023).

The reservoir associated with the dam usually creates a groundwater mounding (Fig. 1), which can extent several kilometres up- and downstream of the reservoir (Berthelote, 2013). Consequently, the removal of the dam can deplete groundwater level, reduce the productivity of tube wells or dry them up, or dry up wetlands (USDA, 2010; Berthelote, 2013; Learn, 2011 [popular newspaper]). The influence of dam reservoirs on groundwater and on surface – groundwater interactions were evaluated using field observations and numerical modelling for quantifying fluxes responses, the transport of contaminants or the effects of dam on salt-water intrusion (Girard et al., 2003; Constantz and Essaid, 2007; Ashraf et al., 2007; Berthelote, 2013; Chang et al., 2019; Aberg et al., 2022; Fang et al., 2022; Li et al., 2023). However, such studies require a large number of field observations, which are not available in most places. Therefore, in these cases, a simpler approach is needed to model changes in groundwater levels due to dam removal, which can be based on

<figure><img src="image_3.png" /><figcaption>Fig. 1. Conceptual model of the theoretical examples presented for reservoir-induced groundwater mound created by a dam along a river.</figcaption></figure>

analytical solutions.

The impact of fluctuations in the water level of a river, a lake or a sea on groundwater level is a well-known theoretical problem. Its mathematical formulation leads to a set of differential equations for which analytical solutions can be found in various textbooks and other works
(Carslaw and Jaeger, 1959; Polubarinova-Kochina, 1962, 1977, Pinder et al., 1964; de Marsily, 1986; Barlow and Moench, 1998; Bruggeman, 1999; Hayek, 2019; Xin et al., 2020; Nikoletos and Katsifarakis, 2024, etc.). Although solutions are given for a variety of hydrogeological situations: infinite to channelized aquifers, stream of infinite length, isolated rectangular ponds, partially clogged streambed, etc., no analytical solution has been proposed to evaluate the distribution of the hydraulic head induced by a reservoir created by a dam, or its changes because of its removal. To bear a maximum of resemblance to field systems, such a solution must take into account both the reservoir-aquifer interaction and those between the stream and the aquifer up- and downstream of the reservoir, and that groundwater-surface water exchanges are not necessarily perfect and the same in the reservoir and in the stream. In addition, for a closer link with the real conditions, the solution must also consider that the height of the water level in the reservoir can vary spatially. This is the subject of this study.

This study seeks steady-state semi-analytical solutions based on the method of fundamental solutions (MFS), for computing hydraulic-head disturbance created by a watercourse including a triangular reservoir of finite-size. MFS is used in the field of physical and engineering problems to overcome the difficulty of finding an analytical solution in complex systems (Chaiyo et al. 2011, Young et al., 2006). It is a meshless method based on the fundamental solutions of governing equations (Kupradze and Aleksidze 1964, Golberg, 1995), which gives it true advantages over conventional models (no model to build, no mesh generation, no discretization of physical boundaries), and allows, among other things, to consider various boundary conditions along the same discontinuity. Its basic concept is to decompose the solutions of a partial differential equation by superposition of the fundamental solution with proper strengths, depending upon the location of boundary conditions and their types. It assumes that a linear combination of the strengthened appropriate fundamental solution approximates the unknown solution sought. The accuracy of the method depends only on the distances between pairs of points where the fundamental solution is applied. In groundwater hydraulics, but not only, fundamental solutions are typically line-sink or point-source solutions. In the field of groundwater hydraulics, MFS has been rarely used (Wang and Zheng 2015). Kuo (1990), Yeung and Chakrabarty (1993) and Kuo et al. (1994) used this technique to compute drawdown caused by wells pumping inside irregularly shaped aquifers, and Huang and Yeh (2015) for estimating flux filtration on a meandering stream caused by pumping wells. Recently, this technique has been used to model the hydraulic-head disturbance created by underground dams (Dewandel et al., 2023).

In the present study, MFS is used to solve the problem of surface- groundwater flow interaction, and is particularly interesting to take into account variable groundwater-river exchanges along the reservoir and stream banks up- and downstream the reservoir. The aim of this work is to develop steady-state semi-analytical solutions for assessing the two- dimensional groundwater hydraulic head induced by a stream- reservoir system fully penetrating a semi-infinite sloping aquifer (Fig. 1). The solutions takes into account that water infiltrates from the reservoir to the aquifer and that a certain amount of this infiltrated water can return to the river. The solution also considers that the surface- groundwater exchanges, i.e. between the reservoir and the aquifer, and between the stream and the aquifer, are not necessarily perfect, and that the water height in the reservoir varies spatially. The solution is used to assess the influence of reservoir-aquifer exchanges (reservoir leakage parameter), the width of the reservoir (i.e. of the dam), and the width of the aquifer. Theoretical formulation is also given in case of a pumping/injecting well near the reservoir, to examine the behaviour of composite devices (Fig. 1). Analytical results are applied to two field examples, one from the literature (Berthelote, 2013) and the other from the Massillargues-Attuech site (France, Frissant and Ladouche, 2022). A numerical analysis of the method used, including sensitivity tests and spatial variation of strength coefficients, is shown. The degree of error introduced by the method of linearization used to solve the governing

partial differential equation (first method of linearization or h-linearization) is also discussed.

The proposed semi-analytical solutions can be implemented as operational tools for engineers assessing the impact on groundwater levels due to the implementation or removal of dams on rivers. Although these solutions cannot represent reservoir-aquifer and stream-aquifer interactions to the same degree of detail as a numerical model, they can be a useful tool for examining the influence of various factors and obtaining estimates commensurate with the level of currently available data.

## 2. Mathematical statements

### 2.1. Hydraulic head solution in the aquifer induced by a reservoir on the watercourse

This study considers an unconfined, sloping, homogeneous aquifer with a saturated thickness without watercourse $D_0$ (in m; Fig. 2). The aquifer dip and the watercourse are aligned with the x-axis direction and $\theta$ is the dip angle of the impermeable basement (in radian). The natural flow direction, and thus the hydraulic head, is parallel to the aquifer bottom, the natural (no reservoir, no stream) hydraulic gradient, $i$, being equal to tan$\theta$. The aquifer is characterized by a constant and uniform hydraulic conductivity ($K$ in m/s) and is limited in space along the y-axis by a no-flow boundary at a distance $L$ (in m) from a Cartesian coordinate system $(x, y)$, resulting in a semi-infinite aquifer. Therefore, $L$ characterizes the width of the aquifer. The reservoir created by the dam has a triangular shape and fully penetrates the aquifer. It is located at the centre of the coordinate system (Fig. 2a), and is defined by its half-length along the x-axis ($R$ in m), and its length along the y-axis ($L_R$ in m). $\theta$ is the angle that the diagonal of the reservoir makes with the x-axis (tan $\theta = L_R/2R$). The water depth in the reservoir is variable, $b$ downstream of the reservoir at the dam location (at $x = R$) and $a$ upstream (at $x = -R$), between the two the water depth changes linearly (Fig. 2b). $a$ and $b$ can take any value. A semi-pervious layer of insignificant storage capacity lies between the reservoir and the aquifer, it is defined by its hydraulic conductivity $k'_R$ (m/s) and its thickness $b'_R$ (m). The river up- and downstream of the reservoir is aligned with the x-axis and is of infinite extent, it is assumed to have a constant water height above $D_0$, $c_{av}$ (m). The stream-aquifer exchanges are also controlled by semi-pervious layers with hydraulic conductivity $k'_{up}$ (m/s) upstream and $k'_{dw}$ (m/s) downstream, and thicknesses $b'_{up}$ and $b'_{dw}$ (m) upstream and downstream respectively.

This scheme therefore considers that there are reservoir-aquifer interactions, and stream-aquifer interactions, both up- and downstream of the reservoir. As in this case the watercourse constitutes a constant hydraulic head or an imposed flux condition, the system is symmetrical according to the x-axis.

The governing steady-state groundwater flow equation corresponding to this situation is defined by the non-linear Boussinesq equation:

$$ \frac{\partial}{\partial x} \left( K D \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial y} \left( K D \frac{\partial h}{\partial y} \right) = 0 \qquad (1) $$

With $D$ the saturated thickness (in m) and $h(x,y) = D(x,y) - ix$ the hydraulic head (in m). The boundary conditions are: $D(\pm\infty, 0) = D_0$ thus $h(\pm\infty, 0) = \pm\infty$, and along the no-flow boundary, $\frac{\partial h}{\partial y}|_{y=L} = 0$. $D_0$ is the saturated thickness without reservoir and stream, and $i = \tan\theta$.

Along the reservoir and the river up- and downstream of the reservoir, the surface- groundwater interactions assume the Robin condition (or Fourier condition), therefore along the banks (see Supplemental materials):

$$ \left( \sin\varphi \frac{\partial h}{\partial x} - \cos\varphi \frac{\partial h}{\partial y} \right) K = \frac{k'_k}{b'_k} (h_0 - h) \qquad (2) $$

where $\varphi$ is the angle between the normal to the bank and the y-axis.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Sketch of the conceptual model of reservoir-stream system with a pumping well. a) Plan view with the displays of image wells; the equations correspond to the boundary condition along the reservoir and stream banks. σ: distance between image wells. b) Section view of the water height in watercourse along the x-axis. c) Section view of the water level in the aquifer parallel to the x-axis and near the pumping well (0 < y < L). d) and e) sections view parallel to the y-axis, near the stream downstream of the reservoir (d) and at the dam location (e).</figcaption></figure>

Therefore, $\varphi = \vartheta$ along the oblique bank of the reservoir ($-R \le x < R$), $\varphi = 3\pi/2$ on the downstream part of the reservoir (along $L_R$, $x = R$), and $\varphi = 0$ along the river up- and downstream of the reservoir. $k'_R$ and $b'_R$ are the hydraulic conductivities and the thicknesses of semi-pervious layers at the contact between the watercourse and the aquifer as mentioned above. $h_0$ is the hydraulic head in the watercourse, which depends on its water height ($s_0$):

$$h_0 = D_0 - ix_0 + s_0 \quad (3)$$

$x_0$ is the distance along the x-axis, which depends on the location so that $h(x,y)$ and $h_0(x_0,y)$ in Eq. (2) are along the same normal to the bank. In the stream upstream ($x < -R$) and downstream ($x > R$) of the reservoir,

$x_0 = x$, in the reservoir, along the oblique bank ($-R \le x < R$, $\varphi = \vartheta$), $x_0 = x + b'_R \cdot \sin\vartheta$, and at the dam location ($x = R$, $\varphi = 3\pi/2$), $x_0 = x - b'_R$. Concerning the water height ($s_0$), in the reservoir ($-R \le x \le R$), and assuming that the water height changes linearly along the x-axis, between the value $a$ (upstream, $x = -R$) and the value $b$ (downstream, $x = R$), $s_0 = \frac{1}{2}[x'_R(b-a)+a+b]$. In the stream upstream ($x < -R$) and downstream ($x > R$) of the reservoir, $s_0 = c \cdot av$.

Notice that the Robin condition admits two extreme cases, the Dirichlet condition (constant head, thus a perfect exchange with the aquifer), when $k'_R \to \infty$, therefore $h = h_0$ in Eq. (2), and the no-flow condition (no exchange), when $k'_R = 0$, therefore $(\cos\varphi \frac{\partial h}{\partial y} - \sin\varphi \frac{\partial h}{\partial x}) K = 0$ in
Eq. (2).

Assuming that the aquifer slope is small ($i = \tan\theta \approx \theta$; $\theta \le 0.02$), which implies that groundwater flow is parallel to the substratum, a linearized form of Eq. (1) can be found by defining $s(x, y) = D(x, y) - D_0$ and assuming that $(D - D_0)/D_0$ is small (i.e. the first method of linearization – or h-linearization – is used to solve Eq. (1); e.g. Zlotnik et al. 2017, Dewandel et al., 2023). Therefore, Eq. (1) can be rewritten as:

$$KD_0 \left( \frac{\partial^2 s}{\partial x^2} + \frac{\partial^2 s}{\partial y^2} \right) - K \cdot i \frac{\partial s}{\partial x} = 0 \quad (4)$$

$$\text{with } s(\pm\infty, 0) = 0 \text{ and } \frac{\partial s}{\partial y}|_{y=L} = 0.$$

The uncertainty introduced by this linearization method is discussed later.

Boundary conditions along the banks (Eq. (2)) becomes (see Supplemental materials):

$$K \left( \sin\varphi \frac{\partial s}{\partial x} - \cos\varphi \frac{\partial s}{\partial y} - i \cdot \sin\varphi \right) = - \frac{k'_k}{b'_k} (s - s_0 + \sin\varphi \cdot i \cdot b'_k) \quad (5)$$

where $s_0$ is the water height in the watercourse (reservoir and stream), and has been defined previously.

The flux conditions through the banks of the reservoir and the stream can be represented by applying the method of fundamental solutions. This assumes that the solution of the equation describing the hydraulic head in the sloping aquifer with a reservoir and a stream can be obtained from a linear combination of strengthened terms produced by applica-

$$H_j(x, y, x_j, y_j) = e^{\alpha(x-x_j)} [K_0(R_{j,1}\alpha) + K_0(R_{j,2}\alpha)] \quad (8)$$

the fundamental solution applied at the $j$th image wells, $x_j$ and $y_j$ the coordinates of the $j$th image wells (up to $n$), and $\gamma_j$ are strengths applied to this solution that must be evaluated to represent the reservoir-stream boundary (Fig. 2a). $R_{j,1}$ and $R_{j,2}$ are distances to the image wells

$$(R_{j,1} = \sqrt{(x-x_j)^2 + (y-y_j)^2} \quad \text{and} \quad R_{j,2} = \sqrt{(x-x_j)^2 + (2L - (y-y_j))^2}).$$

Therefore, the hydraulic head in the sloping aquifer with a reservoir-stream system is given by:

$$h(x, y) = \frac{1}{2\pi D_0 K} \sum_{j=1}^{n} \gamma_j \cdot H_j(x, y, x_j, y_j) + D_0 - i \cdot x \quad (9)$$

If image wells are regularly located within the watercourse at a distance $d'$ of its banks, the distance between them is: $\sigma = (2P + \frac{2R-d}{\cos\theta} + \tan\theta \cdot (2R - d') - d')/n$, with $P$ a distance which must be large enough not to influence the resolution of the system, i.e. greater than the impact of reservoir on the aquifer. In practice, we assume that $P = 10 \cdot (2 \cdot R)$. Note that when $\gamma_j = 0$, Eq. (9) corresponds to the hydraulic head without the reservoir-stream system.

The solution given by Eqs. (7) and (8) has to respect the boundary condition given by Eq. (5). It implies:

$$\frac{1}{2\pi D_0 K} \sum_{j=1}^{n} \gamma_j \left( \sin\varphi \frac{\partial H_j}{\partial x} - \cos\varphi \frac{\partial H_j}{\partial y} + \frac{1}{R_{r-k}} H_j \right) =$$

$$\frac{s_0}{R_{r-k}} + i \cdot \sin\varphi \left( 1 - \frac{b'_k}{R_{r-k}} \right) \quad (10)$$

F

p

tion of the appropriate fundamental solution at a series of image wells. Along the banks of the reservoir and the stream, the superposition of the perturbations induced by each image well has to respect the boundary condition given by Equation (5) (Fig. 2a).

The fundamental solution of the linearized equation (Eq. (4)) describing the steady-state hydraulic head in an inclined aquifer, bounded by a no-flow boundary along the y-axis ($\frac{\partial h_j}{\partial y}|_{y=L} = 0$) with a pumping well, can be found in Hantush (1964a,b) where an image well has been added to represent the no-flow boundary at a distance $L$ (use of the principle of superposition, see Supplemental materials). With the used coordinate system (Fig. 2a), the form of the fundamental solution is given by the saturated thickness $D$:

$$s(x, y) = D(x, y) - D_0 = \frac{Q}{2\pi D_0 K} e^{\alpha x} [K_0(r\alpha) + K_0(r_n\alpha)] = \frac{Q}{2\pi D_0 K} H(x, y) \quad (6)$$

where $Q$ is the pumping rate ($Q<0$; in m³/s), $K_0(u)$ is the modified Bessel function of the second kind of the zero order, and $e$ is the exponential function, $\alpha = i/(2D_0)$, $r = \sqrt{x^2 + y^2}$ and $r_n = \sqrt{x^2 + (2L - y)^2}$. Note that the same solution can be retrieved in Polubarinova-Kochina (1977; p.464–465) and Carlslaw and Jaeger (1959; section 10.7).

According to the method of fundamental solutions, the solution of the whole system (Eq. (4)) and boundary condition given by Eq. (5) can be expressed as:

$$s(x, y) = D(x, y) - D_0 = \sum_{j=1}^{n} \frac{\gamma_j}{2\pi D_0 K} H_j(x, y, x_j, y_j) \quad (7)$$

With where $K_1(u)$ is the modified Bessel function of the second kind of the first

where $R_{r-k}$ are reservoir-stream leakage parameters or retardation coefficients (in m), $R_r = \frac{K}{K_r} b_{R'}$ for the reservoir, and $R_{rup} = \frac{K}{K_{up}} b_{up}'$ and $R_{rdw} = \frac{K}{K_{dw}} b_{dw}'$ for the stream up- and downstream of the reservoir respectively. Note that $R_{r-k} = 0$ assumes that reservoir-stream/aquifer exchanges are perfect, and that $R_{r-k} = \infty$ assumes that there is no reservoir-stream/aquifer exchanges. Eq. (10) is a linear function of $\gamma_j$.

Along the stream, upstream ($x < -R$) and downstream ($x > R$) the reservoir, $\varphi = 0$, $y = 0$ and $s_0 = c \cdot \alpha v$. Along the oblique side of the reservoir: $\varphi = \theta$, $-R \le x < R$, $0 < y < L_R$ and $s_0 = \frac{1}{2} [x(b - \alpha) + a + b]$, and along the downstream part of the reservoir: $\varphi = 3\pi/2$, $x = R$, $0 < y \le L_R$, and $s_0 = b$. Expanded forms of Eq (10) are given in Appendix A. Eq. (10) can be decomposed into two terms: ‘F’ on the left side, and ‘p’ on the right side.

$\frac{\partial H_j}{\partial x}$ and $\frac{\partial H_j}{\partial y}$ are the derivatives of Eq. (8) according to $x$ and $y$.

According to $x$:

$$\frac{\partial H_j}{\partial x} = \alpha e^{\alpha(x-x_j)} \left\{ K_0(R_{j,1}\alpha) + K_0(R_{j,2}\alpha) - (x-x_j) \left[ \frac{K_1(R_{j,1}\alpha)}{R_{j,1}} + \frac{K_1(R_{j,2}\alpha)}{R_{j,2}} \right] \right\} \quad (11)$$

And according to $y$:

$$\frac{\partial H_j}{\partial y} = \alpha e^{\alpha(x-x_j)} \left\{ \frac{2L - (y - y_j)}{R_{j,2}} K_1(R_{j,2}\alpha) - \frac{y - y_j}{R_{j,1}} K_1(R_{j,1}\alpha) \right\} \quad (12)$$
order.

The evaluation of strength coefficients $\gamma_j$ is achieved by transposing the flow condition along the banks of the reservoir and the stream. Strength coefficients $\gamma_j$ were evaluated on the watercourse banks at the closest distance ($d'$) between the image wells and the banks of the watercourse (Fig. 2a). As the flow condition (Eq. (10)) is computed in $n$ locations, it results in a series of $n$ linear equations that have to be solved simultaneously, and which can be represented in matrix form:

$$ \frac{1}{2\pi D_0 K} \begin{bmatrix} F_{11} & F_{21} & F_{31} \cdots F_{n1} \\ F_{12} & F_{22} & F_{32} \cdots F_{n2} \\ F_{13} & F_{23} & F_{33} \cdots F_{n3} \\ \cdots & \cdots & \cdots \\ F_{1j} & F_{2j} & F_{3j} \cdots F_{nj} \\ \cdots & \cdots & \cdots \\ F_{1n} & F_{2n} & F_{3n} \cdots F_{nn} \end{bmatrix} = \begin{bmatrix} \gamma_1 \\ \gamma_2 \\ \gamma_3 \\ \cdots \\ \gamma_j \\ \cdots \\ \gamma_n \end{bmatrix} \quad \begin{bmatrix} p_1 \\ p_2 \\ p_3 \\ \cdots \\ p_j \\ \cdots \\ p_n \end{bmatrix} \quad (13) $$

where 'F' and 'p' are left and right terms of Eq. (10). Subscripts indicate the location where the calculation is performed and the number of the image well. $K$ is hydraulic conductivity (m/s).

Eq. (13) is of the form $\frac{1}{2\pi D_0 K} [\mathbf{A}] [\tilde{\mathbf{\gamma}}_j] = [\mathbf{B}]$, where $[\mathbf{A}]$ represents the equations matrix, $[\tilde{\mathbf{\gamma}}_j]$ is the unknown variable vector and $[\mathbf{B}]$ the known vector. With $[\tilde{\mathbf{\gamma}}_j] = \frac{|\tilde{\mathbf{\gamma}}_j|}{2\pi D_0 K}$ the strength coefficients were evaluated numerically using matrix calculation ($[\tilde{\mathbf{\gamma}}_j] = [\mathbf{B}] \cdot [\mathbf{A}]^{-1}$). Therefore, the hydraulic head is given by:

$$ h(x, y) = \sum_{i=1}^{n} \tilde{\gamma}_j e^{a(x-x_i)} [K_0(R_{j,1}\alpha) + K_0(R_{j,2}\alpha)] + D_0 - i \cdot x \quad (14) $$

Eq. (14) shows that the steady-state hydraulic head induced by the reservoir-stream system does not depend on the hydraulic conductivity of the aquifer but on leakance parameters between the reservoir and the aquifer ($R_r$), and those between the stream and the aquifer up- and downstream of the reservoir ($R_r$ up and $R_r$ dw).

## 2.2. Hydraulic-head solution in the aquifer induced by a reservoir on the watercourse and a pumping well

If a pumping well fully penetrating the aquifer is located near the reservoir created by the dam (Fig. 2), the same method can be used for evaluating the steady-state hydraulic head. Due to the linearity of the solution (Eqs. (4) and (7)), the principle of superposition is used for adding such source. The result is that the term corresponding to the pumping well has to be added to the second part of Eq. (7) and thus to the known vector $[\mathbf{B}]$. Therefore, $[\mathbf{B}]$ and $[\tilde{\mathbf{\gamma}}_j]$, and thus hydraulic head, now also depend on the aquifer hydraulic conductivity. The new system of linear equations obtained is also solved numerically.

For a well pumping (or injecting) at a rate $Q_p$ (in m³/s) located at $x_w$, $y_w$, the term to be added is identical to Eq. (6), except that $x$ is replaced by $x-x_w$, $y$ by $y-y_w$ and $L$ by $L-y_w$. If the well is pumping, $Q_p$ is <0; if it is injecting, $Q_p$ is >0.

Therefore, the semi-analytical steady-state hydraulic-head solution combining the reservoir-stream system with a pumping (or injecting) well is given by:

$$ h(x, y) = \sum_{j=1}^{n} \tilde{\gamma}_j e^{a(x-x_j)} [K_0(R_{j,1}\alpha) + K_0(R_{j,2}\alpha)] + D_0 - i \cdot x + \frac{Q_p}{2\pi D_0 K} e^{\alpha x_p} [K_0(r_{p1}\alpha) + K_0(r_{p2}\alpha)] \quad (15) $$

$$ \text{with } x_p = x - x_w; \quad r_{p1} = \sqrt{(x - x_w)^2 + (y - y_w)^2}, \quad r_{p2} = \sqrt{(x - x_w)^2 + (2L - (y - y_w))^2}. $$

The boundary condition (Eq. (10)) must be rewritten accordingly, and is of the form:

$$ \frac{1}{2\pi D_0 K} \sum_{j=1}^{n} \tilde{\gamma}_j \left( \sin\phi \frac{\partial H_j}{\partial x} - \cos\phi \frac{\partial H_j}{\partial y} + \frac{1}{R_{r-k}} H_j \right) = \frac{s_0 - s_p}{R_{r-k}} + i \cdot \sin\phi \left( 1 - \frac{1}{R_{r-k}} \right) - \sin\phi \frac{\partial s_p}{\partial x} + \cos\phi \frac{\partial s_p}{\partial y} \quad (16) $$

where $s_p$ is the term corresponding to the pumping well (last the term of Eq. (15), $R_{r-k}$ the leakance parameters and $s_0$ the water height in the watercourse (reservoir and stream). The pumping (or injecting) well, can be located everywhere in the aquifer. Notice that if others pumping wells need to be considered, the same superposition principle can be used. Therefore, additional terms corresponding to each well, with their own coordinates and pumping rates, must be added to Eqs. (15) and (16).

# 3. Numerical analysis of the proposed solutions

## 3.1. Accuracy of the method of fundamental solutions

The accuracy of the method used (MFS) depends on the distances between image wells ($\sigma$) and on the distance of these wells to the banks of the reservoir and the stream ($d'$, Fig. 2a). Examples here are given for the case of a reservoir and its stream (Eq. (14)) and for the solution including a pumping well (Eq. (15)).

Fig. 3 gives examples of flow computations according to Eq. (2), $F_1 = K (\sin\phi \frac{\partial h}{\partial x} - \cos\phi \frac{\partial h}{\partial y}) D \cdot \Delta\epsilon$ (flow from the aquifer, in m³/s) and $F_2 = \frac{k_p}{b_p} (h_0 - h) D \cdot \Delta\epsilon$ (flow from the watercourse, in m³/s), with $D$ the saturated thickness in the aquifer near banks (m), $\Delta\epsilon$ an incremental spacing (in m). $\Delta\epsilon = \Delta x$ along the stream up- and downstream the reservoir, $\Delta\epsilon = \Delta x/\cos\theta$ along the oblique side of the reservoir, and $\Delta\epsilon = \Delta y$ along the downstream part of the reservoir. The reservoir-stream system is identical for each case but $d'/\sigma$ ratios differ, with $50 \le n \le 3,000$ ($R = 250$ m, $L_R = 250$ m, $L = 1,250$ m, $i = 0.006$, $K = 10^{-3}$ m/s, $a = c \cdot \alpha v = 1$ m, $b = 4$ m, $R_r = 200$ m and $R_{r \text{ up}} = R_{r \text{ dw}} = 100$ m, $d' = 8$ m, without pumping well Eq. (14). Fig. 4 shows the difference between $F_1$ and $F_2$. Results are satisfactory when conditions of Eq. (2) are satisfied, therefore when $F_1$ is almost identical to $F_2$ (Fig. 3) or when their difference becomes negligible (Fig. 4). Here, for $d'/\sigma$ ratios lower than 0.9, results are clearly not satisfactory and become acceptable for higher ratios. So the question is: what ratio should you use to get satisfactory results?

Fig. 5 shows a series of sensitivity tests of the average flow balance normalized to the mean ($\sum \frac{|F_1 - F_2|}{|F_1 + F_2|}/2$), for different reservoir geometries ($250 \le R \le 750$, $100 \le L_R \le 250$, $0.3 \le a \le 1$, $1.5 \le b \le 13$), different leakance parameters ($20 \le R_r \le 250$, $10 \le R_{r \text{ up}} \le 200$ and $10 \le R_{r \text{ dw}} \le 2,000$), different $d'$ ($4 \le d' \le 15$), different pumping locations and pumping rates (at 15 and 22 m from the watercourse banks, $Q = 60$ and $150$ m³/h) and different number of image wells ($10 \le n \le 3500$). Tests were also done with different natural hydraulic gradients $i$, aquifer thicknesses $D_0$ and aquifer widths $L$ ($0.001 \le i \le 0.01$, $2.5 \le D_0 \le 20$ m, $400 \le L \le 1,500$ m). The figure shows that the accuracy of the proposed semi-analytical solution, i.e. a nil balance between $F_1$ and $F_2$, is clearly controlled by the number of image wells used, which is itself controlled by the geometry of the reservoir-stream system. Obviously, increasing $n$ increases the result accuracy, but for reasonable computation times, from a few tens of seconds to a few minutes, the number of image wells can be optimized according to $d'$. Results show that for a ratio $d'/\sigma$ ranging between 0.7 and 1.0, the semi-analytical solution gives satisfactory results. Usually, values of $d'$ between 5 and 10 m give satisfactory results. However, if pumping wells are located very close to the banks (e. g. a few meters), $d'$ must be decreased (e.g. 2–3 m) to allow a correct description of the flows, and therefore $n$ increased to maintain the $d'/\sigma$ ratio ranges.
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Sensitivity tests performed with Eq. (14) for a reservoir-stream system alone. a) Flow computed along the banks of the reservoir and the stream along the x-axis (red dotted line). b) Flow computed along the banks of the reservoir and the stream parallel to the y-axis (red dotted line, $x = R$; $0 \le y \le L_R$). $F_1$: flow from the aquifer and $F_2$: flow from the watercourse. $>0$: flows from the watercourse towards the aquifer, and $<0$: flows from the aquifer towards the watercourse. D: saturated thickness in the aquifer near banks (m), $\Delta\epsilon$: incremental spacing along the stream. Up- and downstream the reservoir $\Delta\epsilon = \Delta x = 40$ m, along the oblique side of the reservoir $\Delta\epsilon = \Delta x/\cos\theta = 44.7$ m, and $\Delta\epsilon = \Delta y = 6.3$ m along the downstream part of the reservoir. $R = 250$ m, $L_R = 250$ m, $L = 1$ 250 m, $i = 0.006$, $K = 10^{-3}$ m/s, $a = c \cdot av = 1$ m, $b = 4$ m, $R_r = 200$ m, $R_r \text{ up} = R_r \text{ dw} = 100$ m and $d' = 8$ m with $50 \le n \le 3000$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 3.2. Spatial variation of strength coefficients

Fig. 6 shows the spatial variation of strength values (Eqs. (14) and (15) with $\tilde{y}_j = y_j/(2\pi D_0 K)$, in m). In the example given, strengths are positive or negative with increasing absolute values near singularities: corners of the reservoir, at the contact between the reservoir and the stream, and near the pumping well. The strength values are positive along the oblique side of the reservoir and downstream, but negative near the pumping well, near the corner of the reservoir ($x = R, y = L_R$), and upstream of the reservoir. The sign of strengths ($>0$ or $<0$) can also vary depending on the contrasts in hydraulic conductivity between the aquifer and the banks, the location of the pumping well and its pumping rate. Generally, they decrease for banks with low permeability. Additionally, strengths values decrease with increasing number of image wells to counterbalance the number of image wells, and their absolute values tend to be larger where the hydraulic conductivity contrast is high (low permeability of the banks); see Supplementary materials.

Therefore, the spatial variation of strength values is non-linear and depends on the geometry of the system, the hydraulic conductivity contrasts between the aquifer and the banks of the watercourse and, the location of the pumping well and its rate of pumping. However, their values cannot be directly used to assess the actual flows across the watercourse banks. It is the convolution at a given location of the fundamental solution with all the strengths, which makes it possible to evaluate the hydraulic head or the flow across the banks (Eqs. (14) and (15), and Eqs. (10) and (16)).

## 4. Dependence on parameters

As the developed semi-analytical solution allows testing a variety of hydrogeological situations, the solution is used to evaluate the dependence of certain parameters ($R_r, L_R, L$) on groundwater levels due to the removal of the dam (or its emplacement). The hydraulic head are presented in term of variations ($\Delta h$), which are deduced from the
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Sensitivity tests performed with Eq. (14) for a reservoir-stream system alone, difference between the flow from the aquifer ($F_1$) and the flow from the watercourse ($F_2$). a) Along the banks of the reservoir and the stream along the x-axis (red dotted line). b) Along the banks of the reservoir and the stream parallel to the y-axis (red dotted line, $x = R$; $0 \le y \le L$). $R = 250$ m, $L_R = 250$ m, $L = 1$ 250 m, $i = 0.006$, $K = 10^{-3}$ m/s, $a = c\_av = 1$ m, $b = 4$ m, $R_r = 200$ m and $R_{r\ up} = R_{r\ dw} = 100$ m and $d' = 8$ m with $50 \le n \le 3,000$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

comparison between Eq. (14) and the situation without a reservoir ($\Delta h = h_{stream+reservoir} - h_{stream}$). The situation without a reservoir considers a stream on the entire watercourse, i.e. the water height is identical on the entire watercourse and that the leakance parameter at the location of the reservoir is equal to that of the stream upstream of the reservoir ($a = b = c\_av$ and $R_r = R_{r\ up}$ in Eq. (14)). Also, flow in and out of the aquifer across the banks (deduced from Eq. (10)) are presented in term of variations. They are also deduced from the comparison to the situation without reservoir ($\Delta Flow = Flow_{stream+reservoir} - Flow_{stream}$). $\Delta Flow$ are shown for each segment of the watercourse (reservoir and, up- and downstream of the reservoir), and were computed with an arbitrary value of hydraulic conductivity ($K$) of $10^{-3}$ m/s.

In all the cases presented, the length of the reservoir is the same ($2R = 500$ m) as well as the water heights in the reservoir and in the stream ($a = c\_av = 1$ m and $b = 4$ m), the aquifer thickness ($D_0 = 10$ m) and its gradient ($i = 0.006$). For simplicity, the stream leakance factors up- and downstream of the reservoir are set identical ($R_{r\ up} = R_{r\ dw} = 100$ m), which means that assuming a unit value for the thickness of semi-pervious layers, the streambank is 100 times less permeable than the

aquifer. The dependence of these other parameters will be discussed later.

The first two paragraphs examine the influence of the reservoir leakance factor and that of the dam width assuming an infinite aquifer width ($L = \infty$) to avoid superimposed effects.

### 4.1. Dependence on reservoir leakance factor

Fig. 7 shows the dependence on the reservoir leakance factor ($0 \le R_r \le \infty$ m) in the case of an infinite wide aquifer ($L = \infty$). The length of the reservoir along the y-axis is the same for each case ($L_R = 250$ m).

The figure shows that as long as the leakance factor is sufficiently low ($R_r < 10,000$) to allow flow towards the aquifer, the groundwater level increases in the aquifer up to several kilometres up- and downstream of the dam (Fig. 7a), but also laterally along the y-axis, perpendicularly to the stream (Fig. 7b). This increase is greater downstream for two reasons. The first is due to the greater height of water in the reservoir near the dam, and the second is due to the inclination of the aquifer, which creates greater variations downstream than upstream. This is due to the
<table><thead><tr><th>Symbols</th><th>i (-)</th><th>D0 (m)</th><th>R (m)</th><th>L<sub>R</sub> (m)</th><th>L (m)</th><th>b (m)</th><th>a (m)</th><th>c av (m)</th><th>d' (m)</th><th>Rr (m)</th><th>Rr_up (m)</th><th>Rr dw (m)</th><th>K (m/s)</th><th>xw (m)</th><th>yw (m)</th><th>Q (m3/h)</th></tr></thead><tbody><tr><td>●</td><td>0.006</td><td>2.5</td><td>250</td><td>250</td><td>875</td><td>4</td><td>1</td><td>1</td><td>8</td><td>200.0</td><td>100.0</td><td>100.0</td><td>1.0E-03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>●</td><td>0.004</td><td>10</td><td>500</td><td>250</td><td>875</td><td>5</td><td>1</td><td>1</td><td>15</td><td>200.0</td><td>100.0</td><td>100.0</td><td>1.0E-03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>●</td><td>0.01</td><td>5</td><td>250</td><td>100</td><td>1000</td><td>6</td><td>1</td><td>1</td><td>4</td><td>200.0</td><td>100.0</td><td>1000.0</td><td>1.0E-03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>●</td><td>0.01</td><td>15</td><td>400</td><td>100</td><td>400</td><td>13</td><td>1</td><td>1</td><td>10</td><td>100.0</td><td>50.0</td><td>500.0</td><td>1.0E-03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>☒</td><td>0.005</td><td>20</td><td>400</td><td>200</td><td>1500</td><td>13</td><td>1</td><td>1</td><td>7</td><td>20.0</td><td>10.0</td><td>10.0</td><td>1.0E-03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>●</td><td>0.005</td><td>20</td><td>750</td><td>200</td><td>1500</td><td>8</td><td>0.5</td><td>0.5</td><td>8</td><td>250.0</td><td>100.0</td><td>500.0</td><td>1.0E-03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>●</td><td>0.001</td><td>12</td><td>600</td><td>100</td><td>400</td><td>1.5</td><td>0.3</td><td>0.3</td><td>10</td><td>10.0</td><td>200.0</td><td>2000.0</td><td>1.0E-03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>☒</td><td>0.006</td><td>10</td><td>250</td><td>250</td><td>1250</td><td>4</td><td>1</td><td>1</td><td>5</td><td>200.0</td><td>100.0</td><td>100.0</td><td>1.0E-03</td><td>0.0</td><td>150.0</td><td>150.0</td></tr><tr><td>●</td><td>0.006</td><td>10</td><td>250</td><td>250</td><td>1250</td><td>4</td><td>1</td><td>1</td><td>5</td><td>200.0</td><td>100.0</td><td>100.0</td><td>1.0E-03</td><td>-1000</td><td>15</td><td>150.0</td></tr><tr><td>☒</td><td>0.006</td><td>10</td><td>250</td><td>250</td><td>1250</td><td>4</td><td>1</td><td>1</td><td>5</td><td>200.0</td><td>100.0</td><td>100.0</td><td>1.0E-03</td><td>1000</td><td>15</td><td>150.0</td></tr><tr><td>☐</td><td>0.006</td><td>10</td><td>250</td><td>250</td><td>1250</td><td>4</td><td>1</td><td>1</td><td>5</td><td>200.0</td><td>100.0</td><td>100.0</td><td>3.0E-04</td><td>0</td><td>150</td><td>60.0</td></tr><tr><td>☐</td><td>0.006</td><td>10</td><td>250</td><td>250</td><td>1250</td><td>4</td><td>1</td><td>1</td><td>5</td><td>200.0</td><td>100.0</td><td>1000.0</td><td>1.0E-03</td><td>272.4</td><td>125</td><td>150.0</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Sensitivity tests performed with Eqs. (14) and (15), $d'/\sigma$ vs. average flow balance normalized to the mean ($\sum |\frac{F_1-F_2}{(F_1+F_2)/2}|$), for different reservoir geometries, different leakage parameters and aquifer parameters. $250 \le R \le 750$, $100 \le L_R \le 250$, $0.3 \le a \le 1$, $1.5 \le b \le 13$, $20 \le R_r \le 250$, $10 \le R_r \le 200$ and $10 \le R_r dw \le 2,000$, $10^{-3} \le i \le 10^{-2}$, $2.5 \le D_0 \le 20$, $400 \le L \le 1,500$ m, $4 \le d' \le 15$, $Q = 60$ and $150$ m$^3$/h and $10 \le n \le 3,500$. See inserted table for pumping well locations ($x_w$, $y_w$) and prescribed K values.</figcaption></figure>

fundamental solution used, in particular because of its exponential term ($e^{\alpha x}$ giving higher values for $x > 0$).

For low $R_r$ values (<10,000), the spatial extent of the increase in hydraulic head is much higher than the reservoir area, which creates flow from the aquifer to the stream, up- and downstream of the reservoir (Fig. 7c; ΔFlow < 0). Downstream outflows are more important for the reasons mentioned above. However, when the hydraulic conductivity of the reservoir bank becomes too low and impervious $R_r \ge 10,000$, the water levels in the aquifer upstream of the reservoir continue to increase compared to the situation with a stream alone (without reservoir) as the reservoir banks, particularly its oblique side, act as barrier to ground-water flow. Downstream of the reservoir water levels decrease and become lower than that occupied by the watercourse alone, inducing a reversal of flow from stream to aquifer (Fig. 7c; ΔFlow > 0). This situation describes what can happen when the reservoir is silted with low permeability sediments.

## 4.2. Dependence on dam width

Fig. 8 presents examples of hydraulic head variation for different

reservoir geometries (ratio $L_R/R$), from $L_R/R = 2$ (i.e. $L_R = 2R$) to the case where the reservoir is rectilinear and aligned with the stream ($L_R/R = 0$). In these examples, the width of the alluvial valley is also assumed to be infinite ($L = \infty$) and the reservoir leakage parameter is constant ($R_r = 200$ m). As expected (Fig. 8a and b), the more the width $L_R$ increases, the greater the impact in terms of hydraulic head both up- and downstream of the reservoir (parallel to the x-axis) as well as on its sides (parallel to the y-axis). This is explained by the increase of outflows from the reservoir towards the aquifer (Fig. 8c). It results that the hydraulic head increases over increasingly greater distances as this ratio increases, and, consequently, that flows from the aquifer towards the stream increase, in particular downstream (Fig. 8c; ΔFlow < 0).

Note that for ratios $L_R/R \le 0.1$, variations in hydraulic head are almost identical to the case $L_R/R = 0$, which shows for practical applications that an elongated reservoir can be modelled by a rectilinear model without introducing significant errors.

## 4.3. Dependence on aquifer width

Fig. 9 presents a series of calculations of hydraulic head variations
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Plan view of strength coefficients $\tilde{y}_j$ evaluated numerically using matrix calculation (Eqs. (14) and (15)), with $\tilde{y}_j = \gamma_j/(2\pi D_0 K)$, aquifer width $L = 1250$ m). $R_r = 200$ m, $R_{r\ up} = 100$ m, $R_{r\ dw} = 1000$ m, $L_R = 250$ m, $R = 250$ m, $a = c\_av = 1$ m, $b = 4$ m, $D_0 = 10$ m and $i = 0.006$, $K = 10^{-3}$ m/s, $n = 1800$, $d' = 5$ m, $s = 6.00$ m. a) no pumping well. b) $x_w = 0$; $y_w = 150$ m, $Q = 150$ m³/h. c) $x_w = 272.4$, $y_w = 125$ m, $Q = 150$ m³/h. The graphs represent strength coefficients between image wells n°742 (at $x = -799.68$ m) and image wells n°1053 (at $x = 799.68$ m).</figcaption></figure>

for different widths of the aquifer (ratios $L/L_R$, $1.5 \le L/L_R \le \infty$). The length of the reservoir along the y-axis is constant ($L_R = 250$ m), and the reservoir leakance parameter is constant ($R_r = 200$ m). The hydraulic head variations increase when the width of the alluvial valley is shorter (low $L/L_R$ ratios), because the lateral extension of the aquifer prevents the piezometric mound from propagating in the direction perpendicular to stream (i.e. beyond the limit). However, because in these cases the stream-aquifer interactions are quite significant ($R_{r\ up} = R_{r\ dw} = 100$ m), for ratios $L/L_R \ge 10$ the variations converge towards the case where the width is infinite ($L/L_R = \infty$), whether parallel (Fig. 9a) or perpendicular (Fig. 9b) to the stream. If the stream-aquifer interactions had been weaker (increase of $R_{r\ up}$ and $R_{r\ dw}$), the hydraulic head variations would have converged at greater $L/L_R$ ratios (i.e. for a larger valley).

Fig. 9c shows the flow calculations on the different segments of the watercourse vs. $L/L_R$ ratios. The inflows ($\Delta Flow > 0$) and outflows ($\Delta Flux < 0$) of the aquifer increase little for the smallest $L/L_R$ ratios and converge rapidly from ratios greater than 10 due to the stabilization of the piezometric mound. These increases in outflows from the aquifer towards the stream (i.e. they become more negative, $L/L_R<10$), up- and downstream of the reservoir, are due to a greater propagation of the difference in hydraulic head parallel to the watercourse when the width of the aquifer decreases, which increases the gradient between the aquifer and the stream, and therefore the outflows towards the stream. Downstream outflows are always the largest because, as mentioned previously, the variation in hydraulic head propagates further than upstream.

# 5. Application to field cases

## 5.1. The removal of the Milltown dam

The Milltown Dam (western Montana, USA) construction was completed in 1907 at the confluence of the Clark Fork and Blackfoot rivers (Fig. 10a). Over the next 100 years, the Milltown reservoir filled with mining and smelter wastes from areas located upstream. The Milltown Dam was removed during the period of 2006 to 2009, with the objectives to restore groundwater quality, provide fish passage, and return the two rivers to a natural and free-flowing state (Berthelote, 2013). The associated reservoir covered about 90 ha, with a length of

about 1.5 km ($2R = 1,500$ m). According to Berthelote (2013), the water height at the dam in the reservoir was 8.5 m ($b = 8.5$ m) and was 2.4 m in the upstream part of the reservoir ($a = 2.4$ m). The reservoir-river systems are located in an unconfined to semi-confined valley aquifer mainly formed by fluvial sands, gravels and blocks. Its width is about 2 km in average, and its thickness ranges from 6 to 60 m ($D_0$ was set at 30 m) with hydraulic conductivity ranging from $10^{-3}$ to $3 \times 10^{-1}$ m/s. The hydraulic gradient (i) is 0.004.

Fig. 10a illustrates the decline ($\Delta h$ in meters) of the water table from March 2006 to March 2010. It shows that declines were not limited to the immediate reservoir area but extended at least 6 km downstream and 2 to 3 km upstream of the reservoir (Berthelote, 2013). Fig. 10b shows a cross section of the water level variation at 300 m from the stream (contour map and some available observation wells) and modelling tests with different reservoir and stream leakance parameters. The computed hydraulic head variations ($\Delta h$) are deduced from the comparison between Eq. (14) and the situation without a reservoir (i.e. Eq. (14) without reservoir). For the modelling, the width $L$ of the aquifer was set at 1,200 m (the Clark Fork river is almost in the middle of the valley), and the length of the reservoir along the ordinate axis was set at 0 ($L_R = 0$), as, on the right bank, the reservoir and the river are almost aligned. The water height in the stream $c\_av$ was set at 0 (no information on this parameter); the dependence on this parameter is discussed further.

The models that best explain the observed data, with RMSE values between 0.15 and 0.2 m, were obtained with $2,000 \le R_r \le 2,500$ m for the reservoir, $800 \le R_{r\ up} \le 1,000$ m for the stream upstream of the reservoir, and $R_{r\ dw}$ between $10^5$ and $10^6$ m downstream. High values of $R_{r\ dw}$ were necessary to reproduce the very asymmetric observed $\Delta h$ curve, showing that in this case the increase in water height towards the downstream of the reservoir dam, as well as the inclination of the aquifer, was not sufficient to explain the hydraulic head variations. In this case, this shows that downstream of the dam, the stream is almost not connected to the aquifer. This result is consistent with previous studies, which showed that the aquifer is more confined downstream of the dam (Berthelote, 2013). In the reservoir, surface water-aquifer exchanges are lower than in the river upstream of the reservoir, which can be explained by the filling of the reservoir with poorly permeable materials over the last 100 years.

## 5.2. A dam and a well field on the Gardon d'Anduze river (Massillargues-Attuech site)

Massillargues-Attuech town is located about 40 km northwest of the city of Nîmes in France (Fig. 11a). Near the city, a series of five dams were built on the Gardon d'Anduze river in the 1970s until the 1980s to raise the groundwater level of the alluvial aquifer that became depleted due to gravel quarrying. Here, investigations were carried out to evaluate the groundwater mound created by the reservoir associated with one of these dams (Tour de Barre Dam) where a well field is located 200 m downstream of the dam. The well field (3 pumping wells less than 90 m apart) withdraws approximately $250$ m³/h from the groundwater for the supply of drinking water (Fig. 11a). The aquifer mainly consists of gravel with a hydraulic conductivity ($K$) of about $1 \times 5 \times 10^{-3}$ m/s and a thickness of about 5 m near the dam and the well field ($D_0 = 5$ m); Frissant and Ladouche (2022). The length of the reservoir is 830 m ($2R$), and because the reservoir is very elongated its width was set 0 ($L_R = 0$). Its water heights are 3.5 m (b) at the dam and 1.0 m in its upstream part (a), which means that the difference in water height of the watercourse between up- and downstream of the dam is 2.5 m. The water height in the stream was set identical to $a$ ($c\_av = 1.0$ m). Fig. 11a shows the water table map established in July 2023. The map was drawn from local field observations and was constrained by the elevation of the watercourse to have a correct characterization of the groundwater mound created, particularly near the dam. The hydraulic gradient (i) was about $3 \times 10^{-3}$. The width of the alluvial valley (on the right bank) at the dam
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Dependence on the reservoir leakance factor ($0 \le R_r \le \infty$ m) in the case of an infinite wide aquifer ($L = \infty$). $R_{r\ up} = R_{r\ dw} = 100$ m, $L_R = 250$ m, $R = 250$ m, $a = c_{av} = 1$ m, $b = 4$ m, $D_0 = 10$ m and $i = 0.006$. a) Variations of hydraulic head ($\Delta h$), parallel to the x-axis ($y = 300$ m). b) Variations of hydraulic head ($\Delta h$), parallel to the y-axis ($y = 150$ m); colour scale is identical to a). c) Variations of flow in and out of the aquifer compared to the situation without reservoir ($\Delta Flow$), $K = 10^{-3}$ m/s. $\Delta h$ are deduced from the comparison of Eq. (14) to the situation without a reservoir.</figcaption></figure>

location is approximately 610 m ($L = 610$ m). At about 850 m downstream of the dam, there is an underground dam (impervious structure; Fig. 11a), which locally increases the water level upstream of it. However, its influence is estimated to be less than 0.1 m at the location of the well field. The influence on the water level of such an underground structure was studied in Dewandel et al. (2023). The objectives of this field example are to explain the observed water level data with the developed semi-analytical model that considers both a reservoir-stream system and a pumping well, and to assess what the depletion of the groundwater level at the well field location might be if the dam is removed. As only the cumulative pumping rate of the three wells was

known, a single pumping well located in the middle was assumed for the modelling (i.e. Eq. (15)). To match the observed water levels, in this case a constant corresponding to the elevation of the impervious basement at the centre of the system was added to computed values (108.7 m a.s.L.; $x = y = 0$).

Fig. 11b-e compare the results of the modelling of hydraulic head (h) with Eq. (15), with elevation data from the water table map (grey dots) along four cross sections. Two are parallel to the river (n°1 at 95 m from the riverbank – Fig. 11b and n°2 at 260 m – Fig. 11c) and two are perpendicular (n°3 at 70 m upstream of the dam – Fig. 11d, and n°4 at 200 m downstream – Fig. 11e). The models that best explain the
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Dependence on dam width, $0 \le L_R/R \le 2$ in the case of an infinite wide aquifer ($L = \infty$). $R_r = 200$ m, $R_{r \text{ up}} = R_{r \text{ dw}} = 100$ m, $R = 250$ m, $a = c \cdot av = 1$ m, $b = 4$ m, $D_0 = 10$ m and $i = 0.006$. a) Variations of hydraulic head ($\Delta h$), parallel to the x-axis ($y = 300$ m). b) Variations of hydraulic head ($\Delta h$), parallel to the y-axis ($y = 150$ m); colour scale is identical to a). c) Variations of flow in and out of the aquifer compared to the situation without reservoir ($\Delta Flow$), $K = 10^{-3}$ m/s. $\Delta h$ are deduced from the comparison of Eq. (14) to the situation without a reservoir.</figcaption></figure>

observed data (RMSE values <0.5 m) were obtained assuming a perfect exchange between the stream and the aquifer upstream of the reservoir, and the reservoir and the aquifer ($R_r = R_{r \text{ up}} = 0$). Such high surface water – groundwater exchange were necessary to describe the mound near the dam (e.g. Fig. 11b and d). However, downstream exchanges appear to be slightly reduced ($0 \le R_r dw \le 1,000$ m), which is consistent with field observations. Indeed, downstream of the dam the bedrock almost outcrops in the river. The hydraulic conductivity of the aquifer was set at $2.5 \times 10^{-3}$ m/s, a value consistent to existing data, to describe the groundwater mound and, in particular, the depletion due to pumping observed on the profiles (e.g. Fig. 11c). The modelling results satisfactorily explain the observed data. However, in the most downstream area, the modelling underestimates field observations, which can be due to the rise of water level induced by the underground dam (Fig. 11b and c). In the most upstream area, modelled water levels overestimate field observations, however the water table map is poorly constrained in this area (Fig. 11b and c).

**Fig. 12** compares the computed hydraulic head to the situation without the reservoir, i.e. after dam removal ($\Delta h$), for different stream-aquifer exchange parameters downstream of the dam ($0 \le R_r dw \le 1,000$ m). According to this modelling, if the dam is removed, the water level should decrease by 1 to 1.5 m at the location of the well field, which means that the thickness of the saturated aquifer should decrease by 20 to 30 %, and that the pumping rate of the well field would probably decrease by the same amount.

## 6. Discussion

Above, steady-state semi-analytical solutions were developed for evaluating the rise or decline of hydraulic head in an aquifer caused by the implementation of a dam reservoir or its removal, with or without a well pumping in the aquifer (Eqs. (14) and (15)).

It assumed that the rise (or decline) of the water table remains small in comparison to the initial saturated thickness of the aquifer (i.e. $|D - D_0|/D_0$ is small). This method of linearization, known as the first method of linearization (or *h*-linearization), is commonly used for solving groundwater hydraulic problems (e.g. Theis, 1935) and was chosen to solve the standard flow equation in steady-state condition (Eq. (1)). According to previous works (de Marsily, 1986; Moutsoopoulos et al., 2022; Upadhyaya and Chauhan, 1998), it can be assumed that if the degree of saturation/desaturation $\beta = |D - D_0|/D_0$ is <10 %, the results can be considered as exact. Zlotnik et al. (2017) based on various works for the groundwater mounding induced by artificial recharge (Bear, 1972; Korkmaz, 2013 and others), assumed that relative differences are around a few percent if $\beta \le 0.5$ and on about 10 % if $0.5 < \beta \le 1$. To evaluate the accuracy of the solution for higher ratios, the developed solution was compared to a solution that uses the second linearization method-or *h*²-linearization (see Appendix B, Eq. (B-5)). This other solution is a closer approximation of the exact solution, and is applicable when the decrease or increase of the saturated thickness does not exceed one-half of the initial saturated thickness, i.e. when $\beta < 50$ %. However, its computation time is multiplied by a factor between 10 and 100, and sometimes more, which does not make it useful for most practical cases. Fig. 13 presents examples of calculating the difference in hydraulic head ($\Delta h$) induced by a triangular reservoir with the two solutions where only the initial saturated thickness of the aquifer differs, $2.5 \le D_0 \le 100$ m. Results show that the solution with the first linearization method gives very satisfactory results when the maximum degree of saturation/desaturation $\beta$ is less than 20 % ($D_0 \ge 5$ m in the figure), which is in agreement with previous studies (Bear 1972, Korkmaz 2013, Moutsoopoulos et al., 2022, etc.). For higher $\beta$ values, results are approximate, but they can be used as a preliminary diagnosis.

The dependence on reservoir and stream leakage factors ($R_r, R_{r \text{ up}}, R_r dw$), on the width of the reservoir ($L_R$) and on the aquifer width ($L$) were discussed above. In particular, both theoretical and field examples demonstrate that leakage factors are probably the main driving parameters, thus reinforcing the fact that the hydraulic head in the aquifer and, the inflow and outflow of the aquifer are mainly controlled by these river-aquifer exchange parameters. The example on Fig. 3 (with $n \ge 1.500$) demonstrates the importance of these factors on flows, but also on hydraulic head, such as: approximately $4.1 \times 10^{-2}$ m³/s flow from the reservoir to the aquifer, $1.1 \times 10^{-2}$ m³/s return to the river upstream of the reservoir, and $2.3 \times 10^{-2}$ m³/s return to the river downstream of it. The balance of these flows is non-zero (excess of $0.7 \times 10^{-2}$ m³/s), which is expected since this part of the flow contributes to create the piezometric dome. However, in most real cases these flows should be very low compared to the discharge of the river. Therefore, a very small portion of the river discharge should recharge the aquifer; this should occur at a higher rate in configurations with the reservoir than without. Note, that on a portion in the upstream part of the reservoir (Fig. 3a), there are some negative flows (flow from the aquifer to the reservoir). This is explained by the fact that the reservoir creates a groundwater mound over a surface much larger than it, which induces in the
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Dependence on aquifer width, $1.5 \le L/L_R \le \infty$. $R_r = 200$ m, $R_{r \text{ up}} = R_{r \text{ dw}} = 100$ m, $R = 250$ m, $L_R = 250$ m, $a = c \cdot av = 1$ m, $b = 4$ m, $D_0 = 10$ m and $i = 0.006$. a) Variations of hydraulic head ($\Delta h$), parallel to the x-axis ($y = 300$ m). b) Variations of hydraulic head ($\Delta h$), parallel to the y-axis ($y = 150$ m); colour scale is identical to a). c) Variations of flow in and out of the aquifer compared to the situation without reservoir ($\Delta$Flow), $K = 10^{-3}$ m/s. $\Delta h$ are deduced from the comparison of Eq. (14) to the situation without a reservoir.</figcaption></figure>

upstream part a hydraulic head in the aquifer greater than that of the reservoir, thus producing these flows.

These riverbanks properties appear to be a common feature of many systems, as often a zone of low permeability characterizes the part of the aquifer immediately adjacent to the stream or the reservoir (e.g. Butler et al., 2001). Stream-aquifer exchanges up- and downstream of the reservoir play an important role in the extent of the groundwater mound

induced by the reservoir, their increase (decrease of leakage parameters $R_{r \text{ up}}$ and $R_{r \text{ dw}}$) limit its extent because a large part of the water which has infiltrated from the reservoir is diverted towards the stream. On the contrary, when the stream-aquifer exchanges are reduced (increase of $R_{r \text{ up}}$ and $R_{r \text{ dw}}$), the groundwater mound can propagate very far along the stream and perpendicular to it, up to several km in the example of the Milltown dam, because less flow is captured by the river
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** The removal of the Milltown Dam (western Montana, USA). a) Location and map of the variation in hydraulic head (March 2006 to March 2010), redrawn from Berthelote (2013). b) Cross section of the water level variation (Δh) parallel to the x-axis at 300 m from the stream (contour map and some available observation wells; see EE' red dotted line on a)) and modelling tests (Eq. (14); curves) with different reservoir and stream leakance parameters. Δh are deduced from the comparison of Eq. (14) to the situation without a reservoir. R = 750 m, L<sub>R</sub> = 0, L = 1200 m, b = 8.5 m, a = 2.4 m, c, av = 0, i = 0.004, D<sub>0</sub> = 30 m, 1,500 ≤ R<sub>r</sub> ≤ 2,500 m, 800 ≤ R<sub>r up</sub> ≤ 1,000 m and 800 ≤ R<sub>r dw</sub> ≤ 10<sup>6</sup> m. RMSE: residual mean square error. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

downstream. This is also demonstrated by the theoretical example in Fig. 14 (same parameters as Fig. 3), which shows, even for a small dam reservoir (R = 500 m and L<sub>R</sub> = 250 m) where stream-aquifer exchanges are relatively high (R<sub>r up</sub> = R<sub>r dw</sub> = 100 m), that the reservoir-induced groundwater mound propagates far in all directions, much further than the immediate reservoir area.

When the reservoir is silted with low permeability sediments, ex-
changes are reduced and in certain cases, their accumulation can cause
the cancellation of exchanges and creates a water level depletion zone in
the aquifer compared to the situation without a reservoir (Fig. 7c). In the
given example, the depletion zone is restricted downstream of the dam
and is the consequence of the barrier to groundwater flow from the low-

permeability banks, creating upstream a rise of water levels and
downstream a depletion. This situation of low permeability is compa-
rable to the behaviour of impermeable underground dams, which induce
a rise in water levels upstream and a depletion downstream (Dewandel
et al., 2023). However, if the width of the dam reservoir is small
compared to its length (small L<sub>R</sub>/R ratio), the phenomenon of obstruc-
tion of the banks to groundwater flow disappears and due to the absence
of flow from the reservoir towards the aquifer, the depletion zone ex-
tends and covers an area greater than that of the reservoir. This is what
happened in the Var River (south-east of France, Sogreah, 2015). The
emplacement of several dams along the river in the 1970 s and 1980 s
led to a rapid increase in the water table during the first months (up to 4
<figure><img src="image_13.png" /></figure>

**Fig. 11.** A dam and a well field on the Gardon d'Anduze river (Tour de Barre Dam, Massillargues-Attuech site, France). a) Location and water-table map in July 2023 (metres above sea level). b)–e) cross sections of water-table elevation data (grey dots) and modelling tests with different reservoir-aquifer and stream-aquifer lea-kance parameters (Eq. (15)); see orange dotted lines on a), colour scale is given on d). b) n°1 parallel to the x-axis at 95 m from the riverbank. c) n°2 parallel to the x-axis at 260 m from the riverbank. d) n°3 parallel to the y-axis, x = 345 m. e) n°4 parallel to the y-axis, x = 610 m. R = 415 m, L<sub>R</sub> = 0, L = 610 m, b = 3.5 m, a = c<sub>av</sub> = 1.0 m, i = 0.003, D<sub>0</sub> = 5 m, R<sub>r</sub> = R<sub>r up</sub> = 0 m and 0 ≤ R<sub>r dw</sub> ≤ 10<sup>4</sup> m. Pumping well: x<sub>w</sub> = 470 m; y<sub>w</sub> = 200 m and Q<sub>p</sub>: 250 m<sup>3</sup>/h. Note: to match the observed water levels, a constant corresponding to the elevation of the impervious basement at the centre of the system (108.7 m a.s.L.; x = y = 0) was added to the computed values. RMSE: residual mean square error, computed for water-table elevation data (grey dots). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

<figure><img src="image_14.png" /></figure>

**Fig. 12.** Variations of hydraulic head (Δh), parallel to the x-axis at the location of the well field (y = 220 m). Δh are deduced from the comparison of Eq. (15) to the situation without a reservoir. 0 ≤ R<sub>r dw</sub> ≤ 10<sup>3</sup> m, other parameters are given in Fig. 11.

m), followed by a decrease to reach, after less than ten years, the level before the construction of dams and even in some places deeper levels. This shows that in a particular context of high accumulation of low permeability sediments, the emplacement of a reservoir dam can lead, after a certain time, to a depletion of groundwater levels instead of increase them.

Increase of the water heights in the reservoir (a and b) or the length of the reservoir (2·R) increases the hydraulic head and the size of the associated groundwater mound. Logically, flows entering the aquifer towards the reservoir increase, as wells as those from the aquifer to the stream, up- and downstream of the reservoir. Decreasing a, b or R has the opposite effect. The water height in the stream (c<sub>av</sub>) is also an important parameter, which induces modifications in the hydraulic head and flows. However, in terms of variations of hydraulic head (Δh) and flow variations (ΔFlow) compared to the situation without reservoir this parameter does not influence the result, because these are the variations in water heights in the watercourse (e.g. the quantity b-c<sub>av</sub>) that control Δh and ΔFlow. That is why c<sub>av</sub> was set to 0 for the Milltown dam example. If c<sub>av</sub> had been set to 1 m for example and, a and b increased by the same amount, result would have been identical.

The geometry of the reservoir-stream system in the developed
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Comparison between the two methods of linearization (dots: 1st method Eq. (14); lines: 2nd method Eq. (B-5) in Appendix). R = 375 m, i = 0.006, 2.5 ≤ D₀ ≤ 100 m, L = 875 m, L_R = 250 m, a = c Av = 1 m and b = 4 m, R_f = 200 m and R_f up = R_f dw = 100 m. β: maximum degree of saturation (β = (D - D₀)/D₀, D computed with Eq. (B-5)). a) Variations of hydraulic head (Δh) parallel to the x-axis (y = 300.0 m). b) Variations of hydraulic head (Δh) parallel to the y-axis (x = 300.0 m). RMSE_D: standardized square error. Δh are deduced from the comparison of Eq. (14) or Eq. (B-5) to the situation without a reservoir.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** 3D (half-space) graph of hydraulic-head disturbance created by reservoir dam (Eq. (14)). $R = 250$ m, $i = 0.006$, $D_0 = 10$ m, $L = 1,250$ m, $L_R = 250$ m, $a = c_{av} = 1$ m and $b = 4$ m, $R_r = 200$ m and $R_{r\ up} = R_{r\ dw} = 100$ m. Black contour lines: hydraulic head (in m) and, red contour lines ($\Delta h$): difference in hydraulic head compared to condition without reservoir (in m). The arrows indicate the direction of groundwater flow. This theoretical example is the one of Fig. 3 ($n = 3,000$). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

solution is assumed to be triangular or rectilinear ($L_R = 0$). Even if such geometries can be characteristic of a large number of field situations, the model is probably not suitable for many large reservoirs because of their often dendritic shapes. A new model based on the same theoretical development can be implemented to address these new situations.

The use of the developed solution is not limited to the case where the watercourse (stream-reservoir) system feeds the aquifer. In the case where the aquifer feeds the watercourse, the solutions can also be used but in this case it will be necessary to express the results in terms of variations in hydraulic head or flow variations (i.e., $\Delta h$, $\Delta Flow$) to assess what impact the emplacement, removal or leveling of a reservoir dam will have on the aquifer.

A fundamental assumption of this model is that the watercourse system, i.e. the dam reservoir and the stream, fully penetrates the aquifer and is separated from the aquifer by a low permeable zone. This configuration was used by Hantush (1965) to evaluate the streamflow depletion caused by pumping from a well located near a stream with low-permeability streambank in the case of a horizontal aquifer. However, in many situations, the stream may partially penetrate the aquifer or be superficial (Butler et al., 2001; Hunt, 2014, etc.). Compared to the hydraulic theory of partially penetrating wells, it is reasonable to think that in the case of partial penetration of the watercourse, the results of the solution developed will be correct for distances to the banks approximately five times greater than the thickness of the aquifer. This approximation is valid if the ratio between the horizontal and the vertical hydraulic conductivity of the aquifer does not exceed 10. However, if the stream is very superficial (negligible penetration), the match of the observed data with the developed model should, *a priori*, tend to overestimate the reservoir-stream leakance parameters ($R_r, R_{r\ up}$ and $R_{r\ dw}$) to counterbalance the low degree of penetration. Future models should consider this aspect.

## 7. Conclusion

This study proposes steady-state semi-analytical solutions for assessing the long-term spatial distribution of the rise or decline of the hydraulic head in an aquifer caused by a dam reservoir on a watercourse or its removal. The aquifer is assumed unconfined, homogeneous, slanted and finite. The solutions take into account both the reservoir-aquifer interaction and those between the stream and the aquifer up- and downstream of the reservoir. The reservoir can have a triangular or rectilinear (i.e. aligned with the stream) geometry. The reservoir-stream system fully penetrates the aquifer, and is separated from the aquifer by semi-pervious banks (Robin condition). Solutions also include the influence of a pumping well, particularly to assess what the depletion of the water level at the well location might be if the dam is removed. This has been tested on the well field located near the Tour de Barre Dam (Massillargues-Attuech site, France).

The solutions developed are based on the method of fundamental solutions (MFS) and this work highlights that this analytical technique is useful to overcome the difficulty of finding an analytical solution in complex systems. Here, MFS was particularly beneficial in accounting variable groundwater-river exchanges along the banks of the reservoir and the stream up- and downstream the reservoir. The accuracy of the method depends upon the distance between image wells ($\sigma$) and their distance to the reservoir-stream banks ($d'$). Sensitivity tests, with and without a pumping well, show that the results are reliable for $d'/\sigma$ ratios ranging between 0.7 and 1.0. It is also demonstrated that, without the influence of a pumping well, the hydraulic head does not depend strictly on aquifer hydraulic conductivity, but on reservoir leakance parameter ($R_r$) and stream leakance parameters up- and downstream of the reservoir ($R_{r\ up}$ and $R_{r\ dw}$).

The solutions assumed that the rise (or decline) of the water table is small in comparison to the initial condition (first linearization method of
the partial differential equation); i.e. $\beta = |D - D_0|/D_0$ is small. Tests were performed with another formulation of the solution, more complex but closer to the exact solution (second method of linearization or $h^2$-linerization), and results showed that the solution with the first linearization method is satisfactory when $\beta < 20\%$, which agrees with previous works. For higher values, results are approximate, but they can be used as a preliminary diagnosis. Otherwise, the solution with the second method of linearization can be used, nevertheless its computation time is multiplied by a factor between 10 and 100, and sometimes more, which does not make it useful for most practical cases.

The solution was used to assess the influence of a variety of factors on hydraulic head and on in- and outflows to the aquifer at the location of the reservoir, and along the stream up- and downstream of the reservoir. It is shown that the reservoir and river leakage parameters up- and downstream of the reservoir, as well as the geometry of the reservoir dam, are the main determining parameters of the extent of the induced groundwater mound. Stream-aquifer exchanges up- and downstream of the reservoir ($R_r$ up and $R_r$ dw) largely controls the geometry of the groundwater mound induced by the infiltration of water from the reservoir towards the aquifer. An increase of these exchanges (i.e. decrease of $R_r$ up and $R_r$ dw) limits the spatial extension of the groundwater mound, because more of the water infiltrated from the reservoir can return to the stream up- and downstream of the reservoir. Conversely, their decrease favors propagating the groundwater mound up- and downstream of the reservoir but also laterally, perpendicular to the stream, because less flow in the aquifer is captured by the stream. For example, for the Milltown Dam (U.S.A.), the variation in hydraulic head caused by the removal of the dam propagated several kilometers up- and downstream of the dam. The influence of the reservoir leakage parameter on hydraulic head has been also tested, and results show that in a particular context of high accumulation in the reservoir of low permeability sediments (i.e. high $R_r$ values), the emplacement of a reservoir dam can lead, after a certain time, to a depletion of groundwater levels instead of its rise. This is what happened in the Var River (France).

The influence of the width of the reservoir ($L_R$) and of the aquifer width ($L$) were also tested. Increasing $L_R$, logically increases the flow from the reservoir to the aquifer, and therefore the size of the groundwater mound particularly downstream of the reservoir. It is also shown that for ratios $L_R/R \le 0.1$, hydraulic head are almost identical to the case where $L_R = 0$ (i.e. reservoir banks aligned with the stream). Therefore, for practical applications, an elongated reservoir can be modelled by a rectilinear model without introducing significant errors. Increasing $L$ allows the flow to propagate further perpendicular to the stream and less far along the stream, which causes a decrease of flow from the aquifer to the stream, particularly downstream of the reservoir. Therefore, large aquifer width ($L$) caused a less pronounced groundwater mound parallel

to the stream, but a larger mound perpendicular to it. The influence of the other parameters (reservoir length, water heights in the reservoir and in the stream) were also discussed.

A fundamental assumption of this model is that the watercourse system, i.e. the dam reservoir and the stream, fully penetrates the aquifer, which may not be the case in many natural conditions. In this case, the solutions developed are foreseen to give correct estimates for distances to the banks at least five times greater than the thickness of the aquifer. When the watercourse is very superficial (negligible penetration), the use of the developed model on field data should, a priori, tend to increase the reservoir-stream leakage parameters to counterbalance the low degree of penetration. Partial penetration of the watercourse could be the subject of future analytical works.

A tool with the proposed semi-analytical solutions will soon be available on the BRGM website.

## Funding information

This study was conducted by BRGM, and was funded by 'Hydraulic impact on groundwater level of dam removal - development of analytical solutions' funded by the French Office for Biodiversity (OFB) and the 'Attuech Underground Dam' project funded by Etablissement Public Territorial de Bassin (EPTB) des Gardons, the French Water Agency Rhône Mediterranean & Corsica (AERMC), Occitanie Region and BRGM.

## CRediT authorship contribution statement

**Benoît Dewandel:** Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Conceptualization. **Sandra Lanini:** Writing – review & editing, Validation, Methodology, Formal analysis. **Nicolas Frissant:** Writing – review & editing, Methodology, Data curation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgements

The anonymous *Journal of Hydrology* reviewer is thanked for his useful remarks and comments, which improved the quality of the paper.

# Appendix A

Expanded forms of the boundary condition (Eq. (10)).

Along the stream, upstream the reservoir ($\varphi = 0, x < -R, y = 0$):

$$ \frac{1}{2\pi D_0 K} \sum_{j=1}^{n} \gamma_j \left( \frac{\partial H_j}{\partial Y} \Big|_{y=0} - \frac{1}{R_{rup}} H_j \Big|_{y=0} \right) = - \frac{c \cdot av}{R_{rup}} \quad (A1) $$

Along the oblique side of the reservoir ($\varphi = \vartheta, -R \le x < R, 0 < y < L_R$):

$$ \frac{1}{2\pi D_0 K} \sum_{j=1}^{n} \gamma_j \left( \sin\theta \frac{\partial H_j}{\partial X} - \cos\theta \frac{\partial H_j}{\partial Y} + \frac{1}{R_r} H_j \right) = \frac{s_0}{Rr} + i \cdot \sin\theta \left( 1 - \frac{b_R}{R_r} \right) \quad (A2) $$

Along the downstream part of the reservoir ($\varphi = 3\pi/2, x = R, 0 < y \le L_R$):

$$ \frac{1}{2\pi D_0 K} \sum_{j=1}^{n} \gamma_j \left( \frac{\partial H_j}{\partial X} \Big|_{x=R} - \frac{1}{R_r} H_j \Big|_{x=R} \right) = - \frac{b}{Rr} + i \left( 1 - \frac{b_R}{R_r} \right) \quad (A3) $$

And along the stream downstream of the reservoir ($\varphi = 0, x > R$):
$$ \frac{1}{2\pi D_0 K} \sum_{j=1}^{n} \gamma_j \left( \frac{\partial H_j}{\partial Y} \Big|_{Y=0} - \frac{1}{R_{rdw}} H_j \Big|_{Y=0} \right) = - \frac{c \cdot av}{R_{rdw}} \quad (A5) $$

# Appendix B

Second method of linearization (h²-linearization) to solve Eq. (1).

Still assuming that the aquifer slope is small ($i = \tan\theta \approx \theta$; $\theta \le 0.02$), defining $Z(x,y) = D^2(x,y) - D_0^2$, and assuming that the quantity $\bar{b} = 0.5(D+D_0)$ is a constant of linearization, the Equation (1) becomes linear in Z and is identical to the solution given by Hantush (1964b) (its Eq. (1b) adapted to steady-state condition with the appropriate coordinate system):

$$ K \left( \frac{\partial^2 Z}{\partial x^2} + \frac{\partial^2 Z}{\partial y^2} \right) - \frac{i}{\bar{b}} K \frac{\partial Z}{\partial x} = 0 \quad (B1) $$

where the constant $\bar{b}$ is estimated after several successive iterations (Hantush, 1964a, b, 1965), see below.

The fundamental solution to this equation is of the form:

$$ Z(x, y) = \frac{Q}{\pi K} e^{\bar{\alpha} x} [K_0(R_1 \bar{\alpha}) + K_0(R_2 \bar{\alpha})] \quad (B2) $$

With $\bar{\alpha} = i/(2\bar{b})$.

The hydraulic head, which replaces Eq. (9), is:

$$ h(x, y) = \sqrt{\frac{1}{\pi K} \sum_{j=1}^{n} \gamma_j e^{\bar{\alpha}(x-x_j)} [K_0(R_{j,1}\bar{\alpha}) + K_0(R_{j,2}\bar{\alpha})] + D_0^2} - ix = \sqrt{\frac{1}{\pi K} \sum_{j=1}^{n} \gamma_j M_j(x, y, x_j, y_j) + D_0^2} - ix \quad (B3) $$

And the boundary condition along the stream and the reservoir, which replaces Eq. (10), is:

$$ \frac{1}{\pi K} \sum_{j=1}^{n} \gamma_j \left( \sin\varphi \frac{\partial M_j}{\partial x} - \cos\varphi \frac{\partial M_j}{\partial y} + \frac{1}{R_{r-k}} M_j \right) = 2\bar{b} \left( \frac{s_0}{R_{r-k}} + i \cdot \sin\varphi \left( 1 - \frac{b_k}{R_{r-k}} \right) \right) \quad (B4) $$

Eq. (B-4) is also solved using matrix calculation. It is of the form $[\tilde{\gamma}_j] = [\mathbf{B}] \cdot [\mathbf{A}]^{-1}$ with $[\tilde{\gamma}_j] = \frac{[\gamma_j]}{\pi K}$.

Finally, the hydraulic head, which replaces Eq. (14), becomes:

$$ h(x, y) = \sqrt{\sum_{j=1}^{n} \tilde{\gamma}_j e^{\bar{\alpha}(x-x_j)} [K_0(R_{j,1}\bar{\alpha}) + K_0(R_{j,2}\bar{\alpha})] + D_0^2} - ix \quad (B5) $$

Evaluations of $\bar{b}$ and $[\tilde{\gamma}_j]$ require several successive iterations (Eq. (B-4)). The first evaluation starts, for example, with $\bar{b} = D_0$ along the stream-reservoir banks, which gives first values of $[\tilde{\gamma}_j]$, then the saturated thickness D along the stream-reservoir banks is computed with these strengthen coefficients by applying a few iterations (Eq. (B-5)). Taking these new values of D for $\bar{b}$ allows evaluating, after a second evaluation, new values of $[\tilde{\gamma}_j]$, and soon on. This process is repeated until the values $[\tilde{\gamma}_j]$ converge. In practice, 5 to 10 iterations are enough, but compared to the solution using the first method of linearization (Eq. (14)) its calculation time is multiplied by a factor of between 10 and 100, and sometimes more.

# Appendix C. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jhydrol.2024.131801.

# References

Åberg, S.C., Korkka-Niemi, K., Rautio, A., Åberg, A.K., 2022. The effect of river regulation on groundwater flow patterns and the hydrological conditions of an aapa mire in northern Finland. J. Hydrol-Reg. Stud. 40, 101044 https://doi.org/10.1016/j.ejrh.2022.101044.

Ashraf, M., Kahlown, M.A., Ashfaq, A., 2007. Impact of small dams on agriculture and Groundwater development: a case study from Pakistan. Agric. Water Manag. 92, 90–98. https://doi.org/10.1016/j.agwat.2007.05.007.

Barlow, P.M., Moench, A.F., 1998. Analytical solutions and computer programs for hydraulic interaction of stream-aquifer systems. USGS Report 98-415A. 98 p. https://pubs.usgs.gov/of/1998/ofr98-415A/.

Baxter, R.M., 1977. Environmental effects of dams and impoundments. Annu. Rev. Ecol. Syst. 8, 255–283. https://doi.org/10.1146/annurev.es.08.110177.001351.

Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier Publishing Company, New York, p. 764.

Bednarek, A.T., 2001. Undamming rivers: a review of the ecological impacts of dam removal. Environ. Manag. 27 (6), 803–814. https://doi.org/10.1007/s002670010189.

Berthelote A.R., 2013. Forecasting groundwater responses to dam removal. Graduate Student Theses, Dissertations, & Professional Papers. 1402. https://scholarworks.umt.edu/etd/1402.

Beyer P.J., 2005. Introduction to the special issue: dams and geomorphology. Geomorphology, v. 71, p. 1-2. Proc. Of the 33rd Binghamton Symp. in Geomorphology held 12-12 Oct.2002. Beyer, P. J. (Ed.). (2005). Dams and Geomorphology. Elsevier.

Bruggeman G.A., 1999. Analytical Solutions of Geohydrological Problems. first ed., vol. 46. G.A. Bruggeman, (Ed.). Elsevier. 959 p. eBook ISBN: 9780080527130.

Butler, J.J., Zlotnik, V.A., Tsou, M.S., 2001. Drawdown and stream depletion produced by pumping in the vicinity of a partially penetrating stream. Ground Water 39 (5), 651–659. https://doi.org/10.1111/j.1745-6584.2001.tb02354.x.

Capart, H., 2013. Analytical solutions for gradual dam breaching and downstream river flooding. Water Resour. Res. 49, 1968–1987. https://doi.org/10.1002/wrcr.20167.
Carslaw, H.S., Jaeger, J.C., 1959. In: Conduction of Heat in Solids, second ed. Oxford Science Publications. Clarendon Press, p. 510.

Chaiyo, K., Rattanadecho, P., Chantasirawan, S., 2011. The method of fundamental solutions for solving free boundary saturated seepage problem. Int. Commun. Heat Mass Transfer 38, 249–254. https://doi.org/10.1016/j.icheatmasstransfer.2011.022.

Chang, Q., Zheng, T., Zheng, X., Zhang, B., Sun, Q., Walther, M., 2019. Effect of subsurface dams on saltwater intrusion and fresh groundwater discharge. J. Hydrol. 576, 508–519. https://doi.org/10.1016/j.jhydrol.2019.06.060.

Collier M., Webb, R.H. Schmidt J.C., 1996. Dams and rivers: primer on the downstream effects of dams, U.S. Geological Survey Circular 1126. https://doi.org/10.3133/cir1126.

Collins, M., Lucey, K., Lambert, B., Kachmar, J., Turek, J., Hutchins, E., Purinton, T., Neils, D., 2007. Stream barrier removal monitoring guide. Gulf Main Council Mar. Environ. 2008. www.gulfofmaine.org/streambarrierremoval.

Constantz, J., Essaid, H., 2007. Influence of groundwater pumping on streamflow restoration following upstream dam removal. Hydrol. Process. 21, 2823–2834. https://doi.org/10.1002/hyp.6520.

de Marsily, G., 1986. In: Quantitative Hydrology for Engineers. Academic Press Inc. p. 215.

Dewandel, B., Lanini, S., Frissant, N., 2023. Steady state semianalytical solutions for assessing the two dimensional hydraulic head distribution induced by an underground dam in a sloping aquifer with artificial aquifer recharge and pumping. Hydrolge. J. 1–16 https://doi.org/10.1007/s10040-023-02734-2.

Fang, Y., Zheng, T., Wang, H., Zheng, X., Walther, M., 2022. Nitrate transport behaviour behind subsurface dams under varying hydrological conditions. Sci. Total Environ. 838, 155903 https://doi.org/10.1016/j.scitotenv.2022.155903.

Farinacci A.J., 2009. Surface water and groundwater exchanges in fine and coarse grained river bed systems and responses to initial stages of dam removal, Milltown, Montana. Graduate Student Theses, Dissertations, & Professional Papers. 10828.

Frissant N., Ladouche B., 2022. Etat des connaissances sur le fonctionnement de deux systèmes aquifères du territoire du SAGE des Gardons et propositions d'amélioration de la connaissance : les alluvions du Gardons et les molasses miocènes du bassin d'Uzes et les calcaires urgoniens de la Fontaine d'Eure (State of knowledge on the functioning of two aquifer systems in the territory of the 'SAGE des Gardons' and proposals for improving knowledge: the alluvial deposits of the Gardon and the Miocènes molasses des les Uzes basin and the Urgonian limestones of the Fontaine d'Eure). Report BRGM/RP-71482-FR; p. 93 http://infoterre.brgm.fr/. In French.

Girard, P., da Silva, C.J., Abdo, M., 2003. River-groundwater interactions in the Brazilian pantanal. The case of the Cuiabá river. J. Hydrol. 283, 57–66. https://doi.org/10.1016/S0022-1694(03)00235-X.

Golberg, M.S., 1995. The method of fundamental solutions for Poisson's equation. Eng. Anal. Bound. Elem. 16, 205–213. https://doi.org/10.1016/0955-7997(95)00062-3.

Graf, W.L., 2005. Geomorphology and American dams: the scientific, social, and economic context. Geomorphology 71, 3–26. https://doi.org/10.1016/j.geomorph.2004.05.005.

Hantush, M.S., 1964a. Hydraulics of wells. In: Chow, V.T. (Ed.), Advances in Hydroscience, vol. 1. Academic Press.

Hantush, M.S., 1964b. Depletion of storage, leakage, and river flow by gravity wells in sloping sands. J. Geophys. Res. 69, 2551–2560. https://doi.org/10.1029/JZ069i012p02551.

Hantush, M.S., 1965. Wells near streams with semi-pervious beds. J. Geophys. Res. 70 (12), 2829–2838. https://doi.org/10.1029/JZ070i012p02829.

Hart, D.D., Poff, N.L., 2002. A special section on dam removal and river restoration. Bioscience 52 (8), 653–655. https://doi.org/10.1641/0006-3568(2002)052[0653:ASSODR]2.0.CO;2.

Hayek, M., 2019. Accurate approximate semi-analytical solutions to the Boussinesq groundwater flow equation for recharging and discharging of horizontal unconfined aquifers. J. Hydrol. 570, 411–422. https://doi.org/10.1016/j.jhydrol.2018.12.057.

Huang, C.S., Yeh, H.D., 2015. Estimating stream filtration from a meandering stream under the Robin condition. Water Resour. Res. 51, 4848–4857. https://doi.org/10.1002/2015WR016975.

Hunt, B., 2014. Review of stream depletion solutions, behavior, and calculations. J. Hydrol. Eng. https://doi.org/10.1061/(ASCE)HE.1943-5584.000076.

Johnson, S.E., Graber, B.E., 2002. Enlisting the social sciences in decisions about dam removal. Bioscience 52, 731–738. https://doi.org/10.1641/0006-3568(2002)052 [0731:ETSSID]2.0.CO;2.

Korkmaz, S., 2013. Transient solutions to groundwater mounding in bounded and unbounded aquifers. Groundwater 51 (3), 432–441.

Kuo, M.C.T., 1990. Predicting pressure distribution in irregularly shaped reservoirs. Soc. Petrol. Eng. Reservoir Eng. 5 (1), 87–94. https://doi.org/10.2118/14315-PA.

Kuo, M.C.T., Wang, W.L., Lin, D.S., Lin, C.C., Chiang, C.J., 1994. An image well method for predicting drawdown distribution in aquifers with irregularly shaped boundaries. Ground Water 32 (5), 794–804. https://doi.org/10.1111/j.1745-6584.1994.tb00921.x.

Kupradze, V.D., Aleksidze, M.A., 1964. The method of functional equations for the approximate solution of certain boundary value problems. USSR Comput. Math. Math. Phys. 4, 82–126. https://doi.org/10.1016/0041-5553(64)90006-0.

Learn S., 2011. Well owners want Pacificorp to pay more for damage done by removal of Condit Dam, The Oregonian. https://www.oregonlive.com/environment/2011/12/well_owners_want_pacificorp_s.html.

Lehner, B., Reidy, L.C., Revenga, C., Vörösmarty, C., Fekete, B., Crouzet, Ph., Döll, P., Endelan, J., Frenken, K., Magome, J., Nilsson, C., Robertson, J.C., Rödel, R., Sindorf, N., Wisser, D., 2011. High-resolution mapping of the world's reservoirs and dams for sustainable river-flow management. Front. Ecol. Environ. 9, 494–502. https://doi.org/10.1890/100125.

Li, Y., Cao, S., Yu, L., Yao, J., Lu, J., 2023. Quantifying the impacts of a proposed hydraulic dam on groundwater flow behaviors and its eco-environmental implications in the large Poyang Lake-floodplain system. J. Environ. Manage. 336, 117654 https://doi.org/10.1016/j.jenvman.2023.117654.

Moutsopoulos, K.N., Papaspyros, J.N., Fahs, M., 2022. Approximate solutions for flows in unconfined double porosity aquifers. J. Hydrol. 615, 128679 https://doi.org/10.1016/j.jhydrol.2022.128679.

Nikoletos, I.A., Katsifarakis, K.L., 2024. The method of images revisited: approximate solutions in wedge-shaped aquifers of arbitrary angle. e2022WR034347 Water Resour. Res. 60. https://doi.org/10.1029/2022WR034347.

Pejchar, L., Warner, K., 2001. A river might run through it again: criteria for consideration of dam removal and interim lessons from California. Environ. Manag. 28, 561–575. https://doi.org/10.1007/s002670010244.

Pinder, J.F., Bredehoeft, J.D., Cooper, H.H., 1964. Determination of aquifer diffusivity from aquifer response to fluctuation in river stage. Wat. Res. Reash. 5 (4), 850–855. https://doi.org/10.1029/WR005i004p00850.

Pizzuto, J.E., 2002. Effects of dam removal on river form and process. Bioscience 52, 683–691. https://doi.org/10.1641/0006-3568(2002)052[0683:EODROR]2.0.CO;2.

Poff, N.L., Hart, D.D., 2002. How dams vary and why it matters for the emerging science of dam removal. Bioscience 52, 659–668. https://doi.org/10.1641/0006-3568 (2002)052[0659:HDVAWI]2.0.CO;2.

Polubarinova-Kochina, P.Y., 1977. Theory of Ground Water Movement. Translated from the Russian by J.M. Roger De Wiest. Princeton University Press, New Jersey Russian, p. 613.

Polubarinova-Kochina P.Y., 1962. Theory of Ground Water Movement. Moscow: Nauka, in Russian. 665p.

Roberts, S.J., Gottgens, J.F., Spongberg, A.L., Evans, J.E., Levine, N.S., 2006. Integrated approach to assessing potential dam removals: an example from the Ottawa River. J. Environ. Manage. Northwestern Ohio. https://doi.org/10.1007/s00267-005-0091-8.

Servière M., 2021. Impact de l'effacement des seuils en rivière sur les eaux souterraines. (Impact of dam removal on groundwater). BRGM/RP-71235-FR, 56 p. http://infote rre.brgm.fr/. In French.

Shafroh, P., Friedman, J.M., Auble, G.T., Scott, M.L., Braatne, J.H., 2002. Potential responses of riparian vegetation to dam removal. Bioscience 52, 703–712. https:// doi.org/10.1641/0006-3568(2002)052[0703:PRORVT]2.0.CO;2.

Shuman, J.R., 1995. Environmental considerations for dam removal alternatives for river restoration: regulated rivers. Res. Manage. 11, 249–261. https://doi.org/10.1002/rrr.3450110302.

SOGREAH, 2015. Etude des incidences de l'abaissement des seuils sur la nappe du Var, entre le seuil 7 et le Bec de l'Esteron (Study of the impacts of lowering of dams on the aquifer on the Var River between the dam 7 and the Bec of l'Esteron'). SOGREAH report N°4240189. 139 p. In French.

Stanley, E.H., Doyle, M.W., 2002. A geomorphic perspective on nutrient retention following dam removal. Bioscience 52, 693–701. https://doi.org/10.1641/0006- 3568(2002)052[0693:AGPONR]2.0.CO;2.

Theis, C.V., 1935. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage. Trans. Am. Geophys. Union 16, 519–524. https://doi.org/10.1029/TR016i002p00519.

Thomas, H.H., 1976. The Engineering of Large Dams. Wiley, London, p. 818.

Upadhyaya, A., Chauhan, H.S., 1998. Solutions of Boussinesq equation in semi infinite flow region. J. Irrig. Drain. Eng. 124 (5), 265–270. https://doi.org/10.1061/(ASCE) 0733-9437(1998)124:5(265).

USACE, 2008. National inventory of dams NID, in United States Army Topographic Engineering Center, ed., http://crunch.tec.mil.nidpublic/webpages/nid.cfm.

USDA, 2010. Environmental Assessment Lower Musconецong River Restoration Project Finesville Dam Vicinity: Somerset. USDA Natural Resources Conservation Service, New Jersey.

Wang, F.Z., Zheng, K.H., 2015. The method of fundamental solutions for steady-state groundwater flow problems. J. Chin. Inst. Eng. https://doi.org/10.1080/ 02533839.2015.1082936.

Whitelaw, E., MacMullan, E., 2002. A framework for estimating the costs and benefits of dam removal. Bioscience 52, 724–730. https://doi.org/10.1641/0006-3568(2002) 052[0724:AFFETC]2.0.CO;2.

Williams, G.P., Wolman, M.G., 1984. Downstream effects of dams in alluvial rivers. U.S. Geol. Surv. Prof. Pap. 1286 https://doi.org/10.3133/pp1286.

World Commission on Dams, 2000. Dams and Development: A New Framework For Decision-Making: London and Sterling. Earthscan Publications Ltd, VA, p. 404.

World Register of Dams of the International Commission on Large Dams (ICOLD): https://www.icold-cigb.org/GB/world_register/world_register_of_dams.asp.

Wyrick, J.R., Rischman, B.A., Burke, C.A., McGee, C., Williams, C., 2009. Using hydraulic modeling to address social impacts of small dam removals in southern New Jersey. J. Environ. Manage. 90 (Supplement 3), S270–S278. https://doi.org/10.1016/j. jenvman.2008.07.027.

Xin Y., Zhou Z., Li M., Zhuang C., 2020. Analytical solutions for unsteady groundwater flow in an unconfined aquifer under complex boundary conditions. Water 12, 75. doi:10.3390/w12010075.
Yeung, K., Chakrabarty, C., 1993. An algorithm for transient pressure analysis in arbitrarily shaped reservoirs. Comput. Geosci. 19 (3), 391–397. https://doi.org/10.1016/0098-3004(93)90093-K.

Young, D.L., Chiu, C.L., Fana, C.M., Tsai, C.C., Lin, Y.C., 2006. Method of fundamental solutions for multidimensional Stokes equations by the dual-potential formulation.

Eur. J. Mech. B/fluids 25, 877–893. https://doi.org/10.1016/j.euromechflu.2006.02.004.

Zlotnik, V.A., Kacimov, A., Al-Maktoumi, A., 2017. Estimating groundwater mounding in sloping aquifers for managed aquifer recharge. Ground Water 55, 797–810. https://doi.org/10.1111/gwat.12530.