# A semi-analytical solution for modelling the impact on the water table at the pumping well location due to the influence of a finite-sized underground physical barrier in an inclined, homogenous and infinite aquifer

Benoît Dewandel, Sandra Lanini, Augustin Gouy, Nicolas Frissant

► To cite this version:

Benoît Dewandel, Sandra Lanini, Augustin Gouy, Nicolas Frissant. A semi-analytical solution for modelling the impact on the water table at the pumping well location due to the influence of a finite-sized underground physical barrier in an inclined, homogenous and infinite aquifer. Advances in Water Resources, 2025, 205, pp.105091. 10.1016/j.advwatres.2025.105091. hal-05235519

HAL Id: hal-05235519

https://brgm.hal.science/hal-05235519v1

Submitted on 2 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A semi-analytical solution for modelling the impact on the water table at the pumping well location due to the influence of a finite-sized underground physical barrier in an inclined, homogenous and infinite aquifer

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

## ARTICLE INFO

### Keywords:

- Semi-permeable barriers

- Subsurface dam

- Underground dam

- Method of fundamental solutions

- Groundwater hydraulics

- Partial lateral occultation

## ABSTRACT

In hydrogeology, an underground physical barrier with partial lateral occultation refers to a low-permeability screen constructed across groundwater flow without reaching the lateral limits of the aquifer, and whose objective is to modify the elevation of the water table. Due to the complexity of the hydraulic-head distribution associated, several practical design issues arise if the exploitation of the resource is to be ensured. Transient and steady-state semi-analytical solutions, based on the method of fundamental solutions, are proposed for evaluating the hydraulic-head distribution caused by a semi-permeable rectangular barrier that fully penetrates an infinite, homogeneous and sloping aquifer, with or without a pumping well. Based on hundreds of simulations with different aquifer, pumping and barrier characteristics, a set of type curves was produced for evaluating i) the duration of equilibrium of the water-table, ii) the length of the barrier to optimize the pumping rate of a given well or its location, and iii) the hydraulic conductivity and width of the barrier. Type curves are used to explain the time-variations in water-table induced by the construction of a railway station basement, and to design a hypothetical underground barrier to raise the water-table level at the location of a pumping well. The variability of barrier characteristics introduced by uncertainty in hydrogeological parameters is also addressed. The analysis of developed solutions includes sensitivity tests and satisfactory comparisons with numerical modelling. Solutions and type-curves can be implemented as operational tools for improving the design of such structures and are meant to supplement existing hydrogeological models.

## 1. Introduction

The reduced availability of water resources due to climate change, combined with the increasing pressure on water resources, exerted by the increasing demand for domestic, agricultural and industrial water supplies, has motivated humans since ancient times to develop various water technologies to augment existing water resources (Mays 2010; Bouwer 2002; Dillon, 2005). Underground physical barriers are one such water engineering technologies (e.g. BCEOM 1978; Hanson and Nilsson, 1986; Furukawa and Kuroda, 2011; Bhusari et al., 2016; Nishigaki et al., 2004; Izady et al., 2021; Liu et al., 2023; Choi et al., 2023; Masayuki, 2024; Apaydin, 2024). Such physical barriers are also used to manage fresh groundwater in coastal aquifers and to mitigate seawater intrusion with the aim of securing the groundwater supply (e.g.

Kaleris and Ziogas, 2013; Abdoulhalik and Ahmed, 2017; Wu et al., 2020; Zheng et al., 2020; Zhang et al., 2021; Abdoulhalik et al., 2022, 2024; Chang et al., 2023; Wu and Lu, 2023; Yin et al., 2023; Gao et al., 2024).

Typically, subsurface physical barriers consist of semi-permeable or impermeable structures that completely penetrate the aquifer down to the impermeable bedrock, and whose length completely obstructs the width of the aquifer. Sometimes, they may partially penetrate the aquifer and cover only a limited part of its width, so the barrier may not completely obstruct groundwater flow. In these latter cases, the underground physical barrier modifies groundwater flow, causing the groundwater level to rise upgradient and decline downgradient of the barrier. Due to the complexity of the groundwater fluxes and hydraulic head distribution associated with such structures, several practical design questions arise if the exploitation of the resource is to be ensured:

<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.</footnote>

<footnote>E-mail address: b.dewandel@brgm.fr (B. Dewandel).</footnote>
List of symbols

<table><tr><td>D</td><td>saturated thickness of the aquifer, in m</td></tr><tr><td>D<sub>0</sub></td><td>saturated thickness of the aquifer without barrier or pumping well, in m</td></tr><tr><td>d'</td><td>distance of image wells to the side of the barrier, in m</td></tr><tr><td>e<sup>x</sup></td><td>exponential function</td></tr><tr><td>h</td><td>hydraulic head, in m</td></tr><tr><td>i</td><td>hydraulic gradient, i=tanθ, dimensionless</td></tr><tr><td>K</td><td>hydraulic conductivity of the aquifer, in m/s</td></tr><tr><td>K<sub>0</sub>(·)</td><td>modified Bessel function of the second kind of the zero order; K<sub>0</sub>(w) = 1/2 ∫<sub>0</sub><sup>∞</sup> e<sup>-γ-<sup>w<sup>2</sup></sup> / 4γ</sup> dγ</td></tr></table>

<table><tr><td>K<sub>1</sub>(·)</td><td>modified Bessel function of the second kind of the first order</td></tr><tr><td>k'</td><td>hydraulic conductivity of the barrier, in m/s</td></tr><tr><td>Lx</td><td>width of the barrier along the x-axis, in m</td></tr><tr><td>Ly</td><td>length of the barrier along the y-axis (perpendicular to the hydraulic gradient i), in m</td></tr><tr><td>Ln</td><td>logarithmic natural function</td></tr><tr><td>n</td><td>number of image wells</td></tr><tr><td>n&#x20D7;</td><td>vector normal to the interface barrier/aquifer oriented toward the barrier</td></tr><tr><td>Q<sub>p</sub></td><td>constant pumping (or injecting) rate, in m<sup>3</sup>/s</td></tr><tr><td>r</td><td>radius, r = √x<sup>2</sup> + y<sup>2</sup>, in m</td></tr><tr><td>R<sub>j</sub></td><td>distances to the image wells, in m</td></tr><tr><td>s</td><td>variation of the aquifer saturated thickness, s = D-D<sub>0</sub>, in m</td></tr><tr><td>t</td><td>time, in second</td></tr><tr><td>T</td><td>aquifer transmissivity (T = D<sub>0</sub>K), in m<sup>2</sup>/s</td></tr></table>

<table><tr><td>x, y</td><td>coordinates of a Cartesian system, in m</td></tr><tr><td>x<sub>p</sub>, y<sub>j</sub></td><td>coordinates of image wells, in m</td></tr><tr><td>x<sub>w</sub>, y<sub>w</sub></td><td>coordinates of the well pumping (or injecting) well, in m</td></tr><tr><td>W(u,w)</td><td>Hantush leaky function, W(u,w) = ∫<sub>u</sub><sup>∞</sup> e<sup>-γ-<sup>w<sup>2</sup></sup> / 4γ</sup> dγ</td></tr><tr><td>α</td><td>parameter, α = i/(2D<sub>0</sub>)</td></tr><tr><td>β</td><td>degree of saturation, dimensionless</td></tr><tr><td>γ<sub>j</sub></td><td>strength coefficient of the j<sup>th</sup> image well, m<sup>3</sup>/s</td></tr><tr><td>γ&#x20D7;<sub>j</sub></td><td>strength coefficient of the j<sup>th</sup> image well, m</td></tr><tr><td>φ</td><td>angle between the normal to the side and the x-axis, in radian</td></tr><tr><td>θ</td><td>dip angle of the impermeable basement, in radian</td></tr><tr><td>σ</td><td>distance between image wells, in m</td></tr><tr><td>Δζ</td><td>incremental spacing, in m</td></tr><tr><td>Δh<sub>t</sub></td><td>maximum hydraulic-head variation at the barrier sides (centre of the barrier) in transient, in m</td></tr><tr><td>Δh<sub>t=∞</sub></td><td>maximum hydraulic-head variation at the barrier sides (centre of the barrier) in steady-state condition, in m</td></tr><tr><td>Δh<sub>wd</sub></td><td>variation of the hydraulic head at the location of the pumping well induced by the construction of the barrier, in m</td></tr><tr><td>Δh<sub>wd_max</sub></td><td>maximum increase of the water-table at the well location, in m</td></tr><tr><td>Δh<sub>Lx.K/(L,y.k')</sub></td><td>maximum hydraulic-head variation at the barrier sides (centre of the barrier) for a semi-permeable barrier, in m</td></tr><tr><td>Δh<sub>Lx.K/(L,y.k')=∞</sub></td><td>maximum hydraulic-head variation at the barrier sides (centre of the barrier) for an impermeable barrier, in m</td></tr><tr><td>Lx.K/(L,y.k')</td><td>hydraulic resistance parameter of the barrier</td></tr></table>

i) After how long will the water-table recover an equilibrium state (i.e. steady-state)? ii) What geometry (length) of the underground physical barrier is required to optimize the pumping rate of a well, or to increase the water level in it and by how much? iii) How far from the pumping well should it be built? or iv) What width of structure and with what type of material, i.e. its hydraulic conductivity, should the structure be constructed?

The influence of underground physical barrier on groundwater were evaluated using field -mainly for barriers that obstructs the width of the aquifer, laboratory observations and numerical modelling for quantifying fluxes responses, the transport of contaminants or the effects of such structures on saltwater intrusion (e.g. Attard et al., 2014; Sun et al., 2019; Liu et al., 2023; Wu and Lu 2023; Masayuki, 2024; Abdoulhalik and Ahmed, 2017a, 2017b; Abdoulhalik et al., 2017, 2022, 2024; Abd-Elaty et al., 2024; Dewandel et al., 2024a). However, such studies require a large number of field observations, which are not available in most places, or numerous simulations with complex numerical models. Therefore, to improve our understanding and enable a large number of rapid simulations, a simpler approach is needed to model changes in groundwater levels due to underground physical barriers, which can be based on analytical solutions. The objective of this study is to examine the behaviour of groundwater flux around an underground physical barrier that obstructs a limited part of the width of a homogeneous aquifer, where the rise in water level upgradient can serve to secure the pumping of the resource. Although an analytical model cannot represent aquifer - underground physical barrier interactions with the same degree of detail as a numerical model, an analytical model can be a useful tool to examine the influence of various factors and, for practical problems, to provide estimates commensurate with the level of data available. This study complements previous analytical works and contributes to broadening the applicability of analytical methods to problems that

were previously only possible with numerical modelling.

Few theoretical studies have presented analytical solutions for evaluating the transient or the long-term impact of disturbances to groundwater flow caused by underground physical barriers. Some of this work has resulted in useful semi-empirical equations for steady-state conditions deduced from laboratory experimental and numerical modelling (Anwar, 1983; Pujades et al., 2012), others consider barriers with zero thickness (Fitts, 1997), or barriers of infinite length with transient calculations (El-Hames, 2012). Recently, Wu et al., (2023) proposed analytical steady-state solutions to assess the depth location of the seawater wedge toe and the maximum safe pumping rate from a well in rectangular coastal aquifers, when a fully penetrating barrier of finite length is constructed. Their solutions were based on the Schwartz-Christoffel conformal mapping method and the potential theory. Additionally, Dewandel et al., (2024a, b) presented steady-state semianalytical solutions for evaluating the hydraulic-head distribution created by an impermeable rectangular underground physical barrier that fully penetrates a sloping aquifer, where the length of the barrier concerns a limited part of the aquifer width. Their solutions were based on the first and the second method of linearization (h-linearization, Dewandel et al., 2024a, and h<sup>2</sup>-linearization, Dewandel et al., 2024b), and utilize the method of fundamental solution (MFS) to solve the no-flux condition through the four sides of the barrier. The solution assuming the first method of linearization was also compared to numerical modelling, providing satisfactory results (Dewandel et al., 2024a). The second linearization method was used to provide a closer approximation of the exact solution and to define the conditions for applying solutions using the first, simplest method.

The objective of the present study is to develop an analytical model from which theoretical curves will be established to answer the practical questions of designing an underground physical barrier posed
previously. Particular attention will be paid to the evaluation of the hydraulic conductivity of the barrier, which in some cases can reduce its efficiency and the amount of the stored water in the reservoir created upgradient of the barrier due to non-negligible leaks (e.g. Shirahata et al., 2022; Masayuki, 2024).

This study seeks transient and steady-state semi-analytical solutions based on the method of fundamental solutions (MFS), for computing hydraulic-head disturbance created by a semi-permeable rectangular underground physical barrier of insignificant storage capacity and a pumping well in an infinite sloping and homogeneous aquifer. The solutions are developed for hydrogeological situations where the aquifer can be considered relatively homogeneous, thus without significant spatial variations in hydraulic conductivity and/or storage coefficient. MFS is a mesh-free method, which is attractively simple and easy to implement in comparison to mesh-based methods. It is based on the fundamental solutions of governing equations (Kupradze and Aleksidze, 1964; Golder, 1995), which gives it true advantages over conventional numerical models (no model to build, no mesh generation, no discretization of physical boundaries). 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; Johansson et al., 2011; Young et al., 2006; Tsai et al., 2006; Cheng and Hong, 2020; Williams, 2021). Its basic concept is to obtain the solution of a partial differential equation and its boundary conditions by superposition of the fundamental solution on imaginary wells 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 (e.g. Tsai et al., 2006; Cheng and Hong, 2020; Dewandel et al., 2024a, 2024b, 2024c). However, MFS does not provide solution in the domain where imaginary wells are implemented. 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) and was addressed to the transient problem of pumping wells in irregularly shaped aquifers or estimating flux filtration along meandering streams with semi-permeable banks caused by pumping wells (Kuo, 1990; Yeung and Chakrabarty, 1993; Kuo et al., 1994; Huang and Yeh, 2015). Most recently, Lin et al. (2025) used this technique to assess the time-varying depth of the freshwater-saltwater interface for island aquifers with complex coastlines under various pumping and natural recharge conditions. This technique has been also used to model the two-dimensional groundwater hydraulic head in steady-state condition created by an impermeable rectangular physical barrier (Dewandel et al., 2024a, 2024b) and the groundwater hydraulic head depletion induced by river dam removal (Dewandel et al., 2024c).

Results of the proposed semianalytical solutions for an underground physical barrier alone and including a pumping well are analysed and discussed. The performance of transient MFS solution has been compared with numerical modelling. The numerical accuracy of the semi-analytical solution is also discussed, as well as the linearization method used to solve the governing partial differential equation. Hundreds of simulations allow to establish dimensionless graphs relating to i) the time required for the water-table to recover equilibrium following the construction of the underground physical barrier, ii) the length of the barrier required to optimize the pumping rate of a given well or its location, and iii) the effectiveness of the structure according to its width and its hydraulic conductivity. The time-to-equilibrium type-curve is then applied to an example taken from the literature (Pujades et al., 2012) and is used to explain the difference in water table between two observation wells located on either side of the basement of a railway station for 290 days after construction stopped. Finally, the type curves are used to design a hypothetical underground physical barrier in an alluvial aquifer to raise the water-table level at the location of a pumping well due to levels being too deep during summer periods, sometimes

causing water shortage. The variability of barrier characteristics introduced by uncertainty in hydrogeological parameters is also addressed.

The proposed solutions may seem relatively complex to solve. However, all calculations reported in this paper were performed with an Excel macro on a personal laptop and took less than 2–3 min, which demonstrates that they can be used for practical applications.

## 2. Mathematical statements

### 2.1. Transient solution

The aquifer sketch considers an unconfined, infinite, sloping, homogeneous aquifer with a saturated thickness $D_0$ (in m) with no barrier no pumping. The dip of the aquifer is aligned with the x-axis direction and $\theta$ is the dip angle of the impermeable basement (in radians). The natural flow direction is parallel to the aquifer bottom, the natural (no barrier, no pumping) hydraulic gradient is $i$. It is assumed that the aquifer slope is small ($i=\tan\theta \approx 0$; $\theta \le 0.02$), therefore that groundwater flow is parallel to the substratum (Dupuit–Forchheimer assumption). The aquifer is characterized by constant and uniform hydraulic conductivity $K$ (in m/s) and storage capacity $S$ (-). The underground barrier has a rectangular shape, fully penetrates the aquifer, is semi-permeable (hydraulic conductivity $k'$ in m/s) and is of insignificant storage capacity. It is located at the centre of a Cartesian coordinate system $(x, y)$, and is defined by its length along the x-axis $Lx$ (in m) and its length along the y-axis, i.e. perpendicular to natural hydraulic gradient, $Ly$ (in m). A pumping well with a constant rate $Q_p$ is located anywhere in the aquifer.

The governing transient non-linear groundwater flow equation is defined by the Boussinesq equation:

$$ \frac{\partial}{\partial x} \left( KD \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial y} \left( KD \frac{\partial h}{\partial y} \right) = S \frac{\partial h}{\partial t} \quad (1) $$

where $h(x, y, t) = D(x, y, t) - ix$ is the hydraulic head (in m), $t$ is the time (second), $K$ the hydraulic conductivity of the aquifer (in m/s), and $S$ its storativity (-). Without the barrier and the pumping well, the initial condition is given by: $D(x, y, 0) = D_0$ thus $h(\pm\infty, \pm\infty, t) = \pm\infty$. At the location of the pumping well, $x_w, y_w$, boundary conditions are $\lim_{r_p \to 0} (r_p \frac{\partial D}{\partial r_p}) = -\frac{Q_p}{2\pi D_0 K}$ with $r_p = \sqrt{(x-x_w)^2 + (y-y_w)^2}$, and its pumping rate $Q_p$ is constant (in m³/s, $Q_p>0$).

The boundary conditions are: $D(\pm\infty, \pm\infty, t) = D_0$ thus $h(\pm\infty, \pm\infty, t) = \pm\infty$. At the location of the pumping well, $x_w, y_w$, boundary conditions are $\lim_{r_p \to 0} (r_p \frac{\partial D}{\partial r_p}) = -\frac{Q_p}{2\pi D_0 K}$ with $r_p = \sqrt{(x-x_w)^2 + (y-y_w)^2}$, and its pumping rate $Q_p$ is constant (in m³/s, $Q_p>0$).

Along the structure, which assumes insignificant storage capacity, the barrier-groundwater interactions assume the Robin condition (or Fourier condition), therefore along the four sides of the barrier (see Appendix A), the boundary condition is:

$$ \frac{\partial h}{\partial \vec{n}^{-}} K = K \left( \cos\varphi \frac{\partial h}{\partial x} + \sin\varphi \frac{\partial h}{\partial y} \right) = \frac{k'}{Lx \text{ or } Ly} (h_{\text{Opp(side)}} - h_{\text{Adj(side)}}) \quad (2) $$

where $k'$ is the hydraulic conductivity of the physical barrier (m/s), $Lx$ its width and $Ly$ its length. $h_{\text{Adj(side)}}$ is the hydraulic head on the side of the barrier where the flux has to be computed and $h_{\text{Opp(side)}}$, the hydraulic head at the opposite side. $\vec{n}^{-}$ is the vector normal to the interface barrier/aquifer oriented outward the aquifer, and $\varphi$ is the angle of the vector according to the x-axis. $\varphi=0$ upgradient of the barrier $(x=-Lx/2)$ and $\varphi=\pi$ downgradient $(x=Lx/2)$, and $\varphi = \frac{3\pi}{2}$ and $\varphi = \frac{\pi}{2}$ for the two other sides $(y=Ly/2$ and $y=-Ly/2)$.

Notice that the Robin condition admits two extreme cases, no barrier, when $\frac{k'}{Lx}$ and $\frac{k'}{Ly} \to \infty$, therefore $h$ is identical to initial condition (e.g. $h(x, y, t) = D_0 - ix$), and the no-flow condition (no exchange), when $\frac{k'}{Lx} = 0$ and $\frac{k'}{Ly} = 0$, therefore $(\cos\varphi \frac{\partial h}{\partial x} + \sin\varphi \frac{\partial h}{\partial y}) K = 0$ in Eq. (2).

Assuming that the aquifer slope is small (i.e. $i=\tan\theta \approx 0$; $\theta \le 0.02$),
which implies that groundwater flow is parallel to the substratum and assuming that $(D - D_0)/D_0$ is small (i.e. the first method of linearization - or h-linearization- is used; e.g. Zlotnik et al., 2017) a linearized form of Eq. (1) can be found by defining $s(x, y, t) = D(x, y, t) - D_0$. 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} = S \frac{\partial s}{\partial t} \quad (3)$$

$$\text{with } s(\pm \infty, \pm \infty, t) = 0 \text{ and } \text{Lim}_{r_p \to 0} \left( r_p \frac{\partial s}{\partial r_p} \right) = - \frac{Q_p}{2\pi D_0 K}$$

Without the barrier and the pumping well, the initial condition is given by: $s(x, y, 0) = 0$. The initial condition for the case where the barrier is implemented at $t = 0$ in an aquifer with a pumping well under steady-state condition is discussed later.

Boundary conditions along the four sides of the underground barrier (Eq. (2)) become (see Appendix A):

For the upgradient ($\varphi = 0$, $x=-Lx/2$, $-Ly/2 \le y \le Ly/2$) and down-gradient ($\varphi = \pi$, $x=Lx/2$, $-Ly/2 \le y \le Ly/2$) sides of the barrier:

$$K \frac{\partial s}{\partial x} \Big|_{x=\pm \frac{Lx}{2}} - \frac{k'}{Lx} \left[ s\left(\frac{Lx}{2}, y, t\right) - s\left(-\frac{Lx}{2}, y, t\right) \right] = i(K - k') \quad (4)$$

And for the two other sides of the barrier ($\varphi = \frac{3\pi}{2}$, $y=Ly/2$, $-Lx/2 \le$

$x \le Lx/2$ and $\varphi = \frac{\pi}{2}$, $y=-Ly/2$, $-Lx/2 \le x \le Lx/2$):

$$K \frac{\partial s}{\partial y} \Big|_{y=\pm \frac{Ly}{2}} - \frac{k'}{Ly} \left[ s\left(x, \frac{Ly}{2}, t\right) - s\left(x, -\frac{Ly}{2}, t\right) \right] = 0 \quad (5)$$

Note that in the absence of a pumping well, or if the well is on the x-axis, the second member of the left side of Eq. (5) is zero due to the symmetry of the hydraulic head distribution with respect to the x-axis (i.

e. s($x, \frac{Ly}{2}, t$) = s($x, -\frac{Ly}{2}, t$), therefore in these cases only $K \frac{\partial s}{\partial y} \Big|_{y=\pm \frac{Ly}{2}} = 0$.

The flux conditions through the four sides of the underground physical barrier are represented by applying the method of fundamental solutions (MFS). The MFS assumes that the solution of the equation describing the hydraulic head in the sloping aquifer with an underground barrier and a pumping well can be obtained from a linear combination of strengthened terms produced by the application of the appropriate fundamental solution at a series of image wells (Fig. 1b). Along the barrier sides, the superposition of the perturbations induced by each image well has to respect the time-dependent boundary conditions given by both Eq. (4) and Eq. (5).

The fundamental transient solution of the linearized equation (Eq. (3)) describing the hydraulic head in an inclined aquifer with a pumping well at the centre of the system, can be found in Hantush (1964a, 1964b). With the used coordinate system (Figs. 1a, b), the form of the

<figure><img src="image_3.png" /><figcaption>Fig. 1. Sketch of the conceptual model for a semi-permeable underground physical barrier and a pumping well in an infinite, homogeneous and sloping aquifer. a) Section view (the hydraulic gradient has been exaggerated) and b) Plan view with the location of image wells; the equation on the right side corresponds to Eq. (2) and the one on the left corresponds to Eq. (7).</figcaption></figure>
fundamental solution is given by the saturated thickness D:

$$D(x, y, t) - D_0 = \frac{Q}{4\pi D_0 K} e^{\alpha x} \int_0^t \frac{1}{t} e^{-a^2 vt - \frac{r^2}{4vt}} dt = \frac{Q}{4\pi D_0 K} e^{\alpha x} \cdot W(u, r\alpha) \quad (6)$$

where $W(u,r\alpha)$ is the Hantush leaky function, $e$ is the exponential function, $Q$ is the constant pumping rate ($Q < 0$; in m³/s), $\alpha = i/(2D_0)$, $r = \sqrt{x^2 + y^2}$, $v = \frac{D_0 K}{S}$ is the diffusivity of the aquifer and $u = \frac{r^2}{4vt}$. Note that the same solution can be found in Carslaw and Jaeger (1959; section 10.7) and Polubarinova-Kochina (1977; p.464–465).

To solve the given problem, i.e., designing an underground physical barrier to improve the pumping of a well, a steady-state solution was chosen to describe the pumping (Hantush 1964a, 1964b). Therefore, the initial condition (Eq. (1) at $t = 0$) is given by the steady-state solution of the hydraulic head of a well pumping in an inclined aquifer (Hantush, 1964a, 1964b).

Due to the linearity of Eq. (3), the principle of superposition is used for adding the pumping well. Therefore, according to the method of fundamental solutions, the hydraulic head in the sloping aquifer with an underground physical barrier (implemented at $t = 0$) and a well pumping for a very long time at a constant rate $Q_p$ is of the form:

$$h(x, y, t) = \sum_{j=1}^{n} \frac{\gamma_j^t}{4\pi D_0 K} G_j(x, y, x_j, y_j, t) + D_0 - i \cdot x + \frac{Q_p}{2\pi D_0 K} e^{\alpha x} K_0(r_p \alpha) \quad (7)$$

with $G_j(x, y, x_j, y_j, t) = e^{\alpha(x-x_j)} \cdot W(u_j, R_j\alpha)$, and $u = \frac{R_j^2 S}{4D_0 Kt}$ and, $R_j = \sqrt{(x-x_j)^2 + (y-y_j)^2}$ are distances to the image wells (Fig. 1b), $x, y$ are the coordinates (in m) where the hydraulic head $h$ has to be calculated. $x_j$ and $y_j$ the coordinates of the $j^{th}$ image wells (up to $n$) where the fundamental solution $G_j$, which depends on time, is calculated to represent the underground barrier, and $\gamma_j^t$ are strengths applied to this solution that must be evaluated for each time value $t$. The last term of Eq. (7) corresponds to pumping steady-state solution where $K_0(\cdot)$ is the modified Bessel function of the second kind of the zero order. If there is no-pumping well the last term of Eq. (7) disappears. To simplify, although image wells can be located anywhere inside the barrier, it was chosen to place them regularly within the barrier at a distance $d'$ of the structure sides (Fig. 1b), the distance between them is: $\sigma = 2(Lx + Ly - 4d')/n$. Here, it is considered that image wells are located on the same place regardless of time $t$, although their location could have been different depending on the time value (Tsai et al., 2006).

The evaluation of time-dependent strength coefficients $\gamma_j^t$ is achieved by transposing the flow condition around the barrier sides. Thus, the solution given by Eq. (7) has to respect the boundary condition given by Eq. (4) and Eq. (5), it implies:

For the side upgradient of the barrier ($x=-Lx/2$, $-Ly/2 \le y \le Ly/2$):

$$\frac{1}{4\pi D_0 K} \sum_{j=1}^{n} \gamma_j^t \left\{ K \frac{\partial G_j(x, y, x_j, y_j, t)}{\partial x} \Big|_{x=-\frac{Lx}{2}} - \frac{k'}{Lx} \left[ G_j\left(\frac{Lx}{2}, y, x_j, y_j, t\right) - G_j\left(-\frac{Lx}{2}, y, x_j, y_j, t\right) \right] \right\} = i(K - k') - K \frac{\partial s_p}{\partial x} \Big|_{x=-\frac{Lx}{2}} + \frac{k'}{Lx} \left[ s_p\left(\frac{Lx}{2}, y\right) - s_p\left(-\frac{Lx}{2}, y\right) \right] \quad (8)$$

Where $s_p$ is the term corresponding to the pumping well (steady-state solution, last term of Eq. (7)).

For the side downgradient of the barrier ($x=Lx/2$, $-Ly/2 \le y \le Ly/2$):

$$\frac{1}{4\pi D_0 K} \sum_{j=1}^{n} \gamma_j^t \left\{ K \frac{\partial G_j(x, y, x_j, y_j, t)}{\partial x} \Big|_{x=-\frac{Lx}{2}} - \frac{k'}{Lx} \left[ G_j\left(\frac{Lx}{2}, y, x_j, y_j, t\right) - G_j\left(-\frac{Lx}{2}, y, x_j, y_j, t\right) \right] \right\} = i(K - k') - K \frac{\partial s_p}{\partial x} \Big|_{x=-\frac{Lx}{2}} + \frac{k'}{Lx} \left[ s_p\left(\frac{Lx}{2}, y\right) - s_p\left(-\frac{Lx}{2}, y\right) \right] \quad (9)$$

For the upper side of the barrier ($y=Ly/2$, $-Lx/2 \le x \le Lx/2$):

$$\frac{1}{4\pi D_0 K} \sum_{j=1}^{n} \gamma_j^t \left\{ K \frac{\partial G_j(x, y, x_j, y_j, t)}{\partial x} \Big|_{y=\frac{Ly}{2}} - \frac{k'}{Ly} \left[ G_j\left(x, \frac{Ly}{2}, x_j, y_j, t\right) - G_j\left(x, \frac{Ly}{2}, x_j, y_j, t\right) \right] \right\} = -K \frac{\partial s_p}{\partial y} \Big|_{y=\frac{Ly}{2}} + \frac{k'}{Ly} \left[ s_p\left(x, \frac{Ly}{2}\right) - s_p\left(-x, \frac{Ly}{2}\right) \right] \quad (10)$$

And for the lower side of the barrier ($y=-Ly/2$, $-Lx/2 \le x \le Lx/2$):

$$\frac{1}{4\pi D_0 K} \sum_{j=1}^{n} \gamma_j^t \left\{ K \frac{\partial G_j(x, y, x_j, y_j, t)}{\partial x} \Big|_{y=-\frac{Ly}{2}} - \frac{k'}{Ly} \left[ G_j\left(x, \frac{Ly}{2}, x_j, y_j, t\right) - G_j\left(x, \frac{Ly}{2}, x_j, y_j, t\right) \right] \right\} = -K \frac{\partial s_p}{\partial y} \Big|_{y=-\frac{Ly}{2}} + \frac{k'}{Ly} \left[ s_p\left(x, \frac{Ly}{2}\right) - s_p\left(-x, \frac{Ly}{2}\right) \right] \quad (11)$$

Eqs. (8)–(11) can be decomposed into two terms: ‘F’ on the left side, and ‘b’ on the right side, thus:

$$\frac{1}{4\pi D_0 K} F(\gamma_j^t, x, y, x_j, y_j, t, Lx, Ly, K, k') = b(K, k', i, D_0, Q_p, x_w, y_w) \quad (12)$$

Due to the complexity of finding simple analytical forms of the derivatives of the function $G_j$ according to $x$ and $y$ (e.g. Chaudhry and Zubair, 1994; Hunt, 2005), they are solved numerically using a finite difference approximation (see Supplemental Material).

Strength coefficients $\gamma_j^t$ are evaluated on the sides of the barrier at the closest distance ($d'$) between the image wells and the sides of the barrier (Fig. 1b). As the flow condition (Eqs. (8)–(11)) is computed in $n$ locations, it results in a series of $n$ linear equations that have to be solved simultaneously. They can be represented in matrix form:

$$\frac{1}{4\pi D_0 K} \begin{bmatrix} F_{11} & F_{21} & F_{31} & \dots & F_{n1} \\ F_{12} & F_{22} & F_{32} & \dots & F_{n2} \\ F_{13} & F_{23} & F_{33} & \dots & F_{n3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ F_{1j} & F_{2j} & F_{3j} & \dots & F_{nj} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ F_{1n} & F_{2n} & F_{3n} & \dots & F_{nn} \end{bmatrix} \begin{bmatrix} \gamma_1 \\ \gamma_2 \\ \gamma_3 \\ \vdots \\ \gamma_j \\ \vdots \\ \gamma_n \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \\ \vdots \\ b_j \\ \vdots \\ b_n \end{bmatrix} \quad (13)$$

Where ‘F’ and ‘b’ are left and right terms of Eqs. (8)–(11). Subscripts indicate the location where the calculation is performed and the number of the image well (Fig. 1). Eq. (13) is of the form $\frac{1}{4\pi D_0 K} [\mathbf{A}] [\gamma_j] = [\mathbf{B}]$, where $[\mathbf{A}]$ represents the equations matrix, $[\gamma_j^t]$ is the unknown variable vector and $[\mathbf{B}]$ the known vector. With $[\gamma_j^t] = \frac{[\gamma_j]}{4\pi D_0 K}$ the strength coefficients were evaluated numerically using matrix calculation ($[\gamma_j^t] = [\mathbf{A}]^{-1} \cdot [\mathbf{B}]$). Therefore, Eq. (7) can be rewritten as:
$$h(x, y, t) = \sum_{j=1}^{n} \tilde{\gamma}_j^T G_j(x, y, x_j, y_j, t) + D_0 - i \cdot x + \frac{Q_p}{2\pi D_0 K} e^{\alpha x_p} K_0(r_p \alpha) \quad (14)$$

Therefore, the hydraulic head depends on the geometry of the underground physical barrier (Lx and Ly), the hydraulic gradient (i), the aquifer diffusivity ($D_0K/S$), the time (t), and on the ratio between the hydraulic conductivity of the aquifer (K) and the barrier (k'). If the barrier is impermeable and under steady-state condition, the hydraulic head depends only on the hydraulic gradient and the geometry of the barrier (Pujades et al., 2012; Dewandel et al., 2024a and 2024b).

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. (7)–(14).

If the pumping starts at the time of the building of the structure (at $t > 0$) at a constant pumping rate, an initial hydraulic head parallel to the bottom of the aquifer is necessary (i.e. $h(x, y, 0) = D_0 - i \cdot x$). The transient solution is given in appendix (Appendix B). In this case, boundaries conditions Eq. (8)–(11) have to be solved accordingly. This case has been used to compare the transient MFS solution with a full numerical modelling.

## 2.2. Steady-state solution

The steady-state solution (i.e. $t \to \infty$) of the distribution of the hydraulic head induced by a semi-permeable underground physical barrier and a pumping well, is given by replacing the function $G_j$ in Eq. (7) by its corresponding steady-state solution $H_j$:

$$H_j(x, y, x_j, y_j) = e^{\alpha(x-x_j)} \cdot K_0(R_j \alpha) \quad (15)$$

Therefore, following the procedure describes above, Eq. (14) becomes:

$$\begin{aligned} h(x, y) &= \sum_{j=1}^{n} \frac{\gamma_j}{2\pi D_0 K} H_j(x, y, x_j, y_j) + D_0 - i \cdot x + \frac{Q_p}{2\pi D_0 K} e^{\alpha x_p} K_0(r_p \alpha) \\ &= \sum_{j=1}^{n} \tilde{\gamma}_j H_j(x, y, x_j, y_j) + D_0 - i \cdot x + \frac{Q_p}{2\pi D_0 K} e^{\alpha x_p} K_0(r_p \alpha) \end{aligned} \quad (16)$$

Notice that because of steady-state condition $H_j$, $\gamma_j$ and $\tilde{\gamma}_j$, do not depend on time. The boundary conditions are similar to Eqs. (8)–(11), but rewritten for steady-state conditions (without time t, function $G_j$ is replaced by $2 \cdot H_j$ and, $\frac{\partial \gamma_j}{\partial x}$ and $\frac{\partial \gamma_j}{\partial y}$ by $2 \cdot \frac{\partial H_j}{\partial x}$ and $2 \cdot \frac{\partial H_j}{\partial y}$, the explanation of factor 2 is given by the relationship between functions $W(u, w)$ [with $t \to \infty$] and $K_0(\cdot)$, see list of symbols). In this case, the derivatives can be calculated analytically and have the following forms:

According to x:

$$\frac{\partial H_j}{\partial x} = \alpha e^{\alpha(x-x_j)} \left\{ K_0(R_j \alpha) - (x-x_j) \frac{K_1(R_j \alpha)}{R_j} \right\} \quad (17)$$

And according to y:

$$\frac{\partial H_j}{\partial y} = -\alpha e^{\alpha(x-x_j)} \frac{y-y_j}{R_j} K_1(R_j \alpha) \quad (18)$$

Where $K_1(\cdot)$ is the modified Bessel function of the second kind of the first order. The derivative of the pumping well solution ($\frac{\partial s_p}{\partial x}$, $\frac{\partial s_p}{\partial y}$) follows the same forms, except that x is replaced by $x-x_w$, and y by $y-y_w$. Then, a new system of linear equations is obtained and solved numerically to provide estimates of strength coefficients $\tilde{\gamma}_j$. Notice that in the absence of a pumping well, the terms corresponding to the pumping well (terms with $K_0$) disappear from Eq. (16) and from the equations corresponding to the boundary conditions.

Steady-state solutions are used to compute the hydraulic head distribution on the long term.

# 3. Numerical analysis of the proposed semi-analytical solutions

## 3.1. Comparison with numerical modelling

The performance of the transient MFS solution was evaluated by comparison with a full numerical modelling implemented using MODFLOW 6 (Langevin et al., 2021). The objective is to evaluate the accuracy of the developed solution, and therefore to evaluate the discrepancies between the hydraulic heads computed with time-dependent strength coefficients applied to each image wells ($\gamma_j^t$) and the hydraulic heads resulting from the full numerical modelling. The numerical model simulates confined groundwater flow within a heterogeneous but axisymmetric domain, discretized using a refined logarithmic mesh structure centred on the barrier. The spatial domain spans 40 km × 40 km in the horizontal plane (x, y plane), with cell sizes varying from 1 m at the barrier's centre to 500 m near the model boundaries. Vertically, the aquifer is discretized into 25 4-m-thick layers, representing a total saturated thickness of 100 m. The aquifer exhibits a slope identical to the hydraulic gradient (e.g., i: 0.005). The top and bottom surfaces of the aquifer therefore decrease linearly from west to east across the domain (i.e. similarly to Fig. 1a). The top surface configuration is matched by the boundary conditions, which consist of constant head boundaries on all four sides. For example, with i: 0.005 a constant hydraulic head of 350 m is imposed on the western boundary and 150 m on the eastern boundary, ensuring a linear hydraulic gradient consistent with both topography and regional flow direction. The northern and southern boundaries are assigned constant height values that are linearly interpolated between these two extremes to maintain internal consistency. The model respects the validity of the proposed solution (Eq. (14)).

In this numerical modelling example (Fig. 2a), the hydraulic conductivity of the aquifer (K) is $10^{-5}$ m/s with storage coefficient 0.1. The underground physical barrier is modelled as a semi-permeable and fully-penetrates the aquifer. Its length according to the y-axis (Ly) is 401 m, its width according to the x-axis 9 m (Lx) and its hydraulic conductivity is $4 \cdot 10^{-7}$ m/s (k'). The pumping well is located 21 m from the centre of the barrier along the x-axis and the pumping starts at $t=0$ at a rate of $Q_p=25$ m$^3$/h ($x_w=-21.0$ m, $y_w=0.0$ m).

In this example, the pumping starts at the beginning of the simulation in order to avoid boundary effects on a steady-state solution. The simulation was set up for up to about 19 years (about $10^7$ min) with time steps defined using a logarithmic progression. In this context, the model's boundaries were not reached and the aquifer remained saturated.

Figs. 2b and c present the results of numerical modelling together with analytically computed drawdown values on semi-log scale. For the semi-analytical solution, the distance of image wells to the side of the barrier (d') was 2.6 m, and the number of image wells (n) was set at 700, but the solution starts to give satisfactory from n = 600; the accuracy of the solution is discussed in the next paragraph. As shown in the figures, drawdown values (s) computed with the MFS solution for observation wells near the pumping well, upgradient and downgradient of the barrier, satisfactorily match the values provided by the full numerical modelling. The standardized root mean square error (RMSE D) was calculated for each numerical-analytical comparison of drawdown values (see inserted table on Fig. 2b). RMSE D ranges between 0.007 and 0.014, showing that the proposed transient MFS solution (Eq. (14)) agrees with the numerical modelling. Other examples are given in appendix for an impermeable barrier ($k'=0$ m/s; Fig. B) and, a semi-permeable barrier with a pumping well not located on the x-axis and other aquifer characteristics(Fig. C); others are provided in the Supplementary Materials. All gives similar results.

## 3.2. Accuracy of the semi-analytical solutions

As classically done in studies using MFS, the accuracy of the method
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Comparison with numerical modelling for a 401 m long semi-permeable barrier ($k'=4.10^{-7}$ m/s). See the text for other model parameters. (a) Plan view with the location of the barrier, the pumping well ($x_w=-21.0$ m, $y_w=0.0$ m) and the observation wells (the location of wells are given by their colours (see inserted table), (b) comparison of drawdown evolution (s) vs. time computed with MODFLOW 6 software (triangles) and the MFS solution (Eq. (14); curves), the colours of the curves and of triangles refer to the previous inserted table. Number of image wells (n): 700, distance of image wells to the side of the barrier ($d'$): 2.6 m. RMSE D: standardized root mean square error, (c) hydraulic head cross-sections along the x-axis (at $y=0$) at different times since the start of pumping.</figcaption></figure>

has been explored (e.g. Huang and Yeh 2015; Tsai et al., 2006; Cheng and Hong 2020; Dewandel et al., 2024a, 2024b, 2024c). Examples here are given for the transient condition of an underground physical barrier with and with no pumping well (Eq. (14)), and for the steady-state condition (Eq. (16)).

Figs. 3a, c, e give examples of groundwater flow computations along the barrier according to Eqs. (4)–(5) at different time after the building of the barrier ($t=1$ h, 10 days and 695 days) for the sides parallel to the y-axis ($x=\pm Lx/2$), $F_1 = K_{\text{dx}}^{\text{dh}} D \cdot \Delta\zeta$ (groundwater flow from the aquifer, in m$^3$/s) and $F_2 = \frac{k}{Lx} (h_{\text{Oppside}} - h_{\text{Adjside}}) D \cdot \Delta\zeta$ (groundwater flow from the barrier, in m$^3$/s), with D the saturated thickness in the aquifer near the

underground barrier (m) and $\Delta\zeta$ an incremental spacing along the y-axis (in m). The semi-permeable underground barrier and the aquifer properties are identical for each case, but $d'/\sigma$ ratios differ, with $500 \le n \le 2800$ (Lx=3 m, Ly=650 m, $i=0.006$, $D_0=20$ m, $K=5 \cdot 10^{-3}$ m/s, $S=0.15$, $k'=5 \cdot 10^{-7}$ m/s, $d'=1.3$ m, without a pumping well Eq. (14)).

Figs. 3b, d, f show the difference between $F_1$ and $F_2$. Results are satisfactory when conditions of Eqs. (4)–(5) are satisfied, therefore when $F_1$ is almost identical to $F_2$ (Fig. 2a, c, e) or when their difference becomes negligible (Fig. 3b, d, f). Here, for $d'/\sigma$ ratios lower than 2.2 at $t=10$ and 695 days, results are clearly not satisfactory and become acceptable for higher ratios. At $t=1$ h, $d'/\sigma \ge 1.7$ gives also satisfactory results. So, the question is: what ratio should be used to obtain
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Sensitivity tests performed with Eq. (14) for a semi-permeable underground barrier alone at different time since the building of the structure ($t = 1$ h, 10 days and 695 days). $Lx=3$ m, $Ly=650$ m, $i=0.006$, $D_0=20$ m, $K=5 \cdot 10^{-3}$ m/s, $S = 0.15$, $k'=5 \cdot 10^{-7}$ m/s ($Lx \cdot K/k' = 3000$ m), $x=\pm Lx/2$, $d'=1.3$ m and $500 \le n \le 2800$. $F_1$: groundwater flow from the aquifer and $F_2$: groundwater flow from the barrier. $D$: saturated thickness in the aquifer near the barrier (m), $\Delta\zeta$: incremental spacing along the barrier ($\Delta\zeta=16$ m). a), c), e) [left] Flow computed along the underground barrier for different $d'/\sigma$ ratios. b), d), f) [right] Difference between the flow from the aquifer ($F_1$) and the flow from the barrier ($F_2$).</figcaption></figure>

satisfactory results, does it depend on the time and/or the hydraulic conductivity of the underground barrier? Afterwards, the barrier is defined by its hydraulic resistance parameter, $Lx \cdot K/(Ly \cdot k')$. An infinite value corresponds to an impermeable barrier, and low value to a semi-permeable barrier.

**Fig. 4** shows a series of sensitivity tests of the average groundwater flow balance normalized to the mean ($\sum \frac{F_1-F_2}{(F_1+F_2)/2}$), for different barrier geometries ($1 \le Lx \le 10$ m, $200 \le Ly \le 1000$ m), different aquifer properties ($0.003 \le i \le 0.018$, $5 \le D_0 \le 25$ m, $10^{-4} \le K \le 5 \cdot 10^{-3}$ m/s, $0.05 \le S \le 0.35$),

different hydraulic resistance parameters ($2.0 \le Lx \cdot K/(Ly \cdot k') \le \infty$), different time since the building of the underground structure ($3.5 \le t \le \infty$, in days), different pumping locations and pumping rates (between 70 and 300 m from the barrier, $60 \le Q_p \le 150$ m³/h), different $d'$ ($0.35 \le d' \le 2.3$ m) and different number of image wells ($80 \le n \le 2800$). As expected, 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 barrier. Obviously, increasing $n$ increases the result accuracy, but for reasonable computation times, from a few tens of
<table><thead><tr><th>Symbol</th><th>K (m/s)</th><th>D<sub>0</sub> (m)</th><th>S (-)</th><th>i (-)</th><th>t (d)</th><th>Ly (m)</th><th>Lx (m)</th><th>d' (m)</th><th>Lx·K/(Ly·k') (m)</th><th>x<sub>w</sub> (m)</th><th>y<sub>w</sub> (m)</th><th>Q<sub>w</sub> (m<sup>3</sup>/h)</th></tr></thead><tbody><tr><td>Symbol: X</td><td>2.0E-03</td><td>10.0</td><td>0.35</td><td>0.01</td><td>6.9</td><td>500</td><td></td><td>3</td><td>1.3</td><td>300.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Symbol: ●</td><td>2.0E-03</td><td>10.0</td><td>0.35</td><td>0.01</td><td>69.4</td><td>1000</td><td></td><td>5</td><td>2.3</td><td>5.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Symbol: ◦</td><td>5.0E-04</td><td>20.0</td><td>0.15</td><td>0.006</td><td>694.4</td><td>650</td><td></td><td>3</td><td>1.3</td><td>4.6</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Symbol: ◆</td><td>5.0E-04</td><td>20.0</td><td>0.15</td><td>0.006</td><td>694.4</td><td>650</td><td></td><td>3</td><td>1.3</td><td>4.6</td><td>-300</td><td>0</td><td>150</td></tr><tr><td>Symbol: ■</td><td>5.0E-03</td><td>15.0</td><td>0.05</td><td>0.018</td><td>∞</td><td>450</td><td></td><td>2</td><td>0.8</td><td>2.2</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Symbol: ▲</td><td>1.0E-03</td><td>7.0</td><td>0.1</td><td>0.005</td><td>3.5</td><td>250</td><td></td><td>10</td><td>1</td><td>2.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Symbol: ▲</td><td>1.0E-03</td><td>7.0</td><td>0.1</td><td>0.005</td><td>3.5</td><td>250</td><td></td><td>10</td><td>1</td><td>2.0</td><td>-150</td><td>0</td><td>80</td></tr><tr><td>Symbol: ●</td><td>2.0E-03</td><td>5.0</td><td>0.12</td><td>0.005</td><td>6944.4</td><td>300</td><td></td><td>6</td><td>1.5</td><td>∞</td><td>-150</td><td>0</td><td>50</td></tr><tr><td>Symbol: +</td><td>1.0E-04</td><td>25.0</td><td>0.02</td><td>0.003</td><td>347.2</td><td>200</td><td></td><td>1</td><td>0.35</td><td>10.0</td><td>-70</td><td>0</td><td>60</td></tr><tr><td>Symbol: ■</td><td>1.0E-04</td><td>25</td><td>0.2</td><td>0.003</td><td>347.2</td><td>200</td><td></td><td>1</td><td>0.35</td><td>10.0</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Sensitivity tests performed with Eq. (14) and Eq. (16) (8th case with $t=\infty$; steady-state), $d'/\sigma$ vs. average groundwater flow balance normalized to the mean ($\sum |\frac{F_1-F_2}{(F_1+F_2)/2}|$), for different underground physical barrier geometries, time $t$, different barrier hydraulic resistance parameters (Lx·K/[Ly·k']) and aquifer parameters. $1 \le Lx \le 10$ m, $200 \le Ly \le 1000$ m, $0.003 \le i \le 0.018$, $5 \le D_0 \le 25$ m, $10^{-4} \le K \le 5 \cdot 10^{-3}$ m/s, $0.05 \le S \le 0.35$, $2.0 \le Lx·K/(Ly·k') \le \infty$ (in m), $3.5 \le t \le \infty$ (in days), $0.35 \le d' \le 2.3$ m, $Q_p=50$ to $150$ m$^3$/h and $80 \le n \le 2800$. See the inserted table for pumping well locations ($x_w, y_w$).</figcaption></figure>

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 2.0 and 3.0 (where $\sigma$ is distance between image wells, $\sigma = 2(Lx + Ly - 4d')/n$), the semi-analytical solutions (Eqs. (14) and (16)) give satisfactory results. This result is similar to that obtained for steady-state solutions with impermeable underground barrier with or without pumping well (Dewandel et al., 2024a, 2024b) and shows that the sensitivity of the semi-analytical solutions depends little on the time $t$ and the resistance parameter of the barrier, but mainly on the geometry of the underground barrier. As mentioned previously, for relatively early time, the semi-analytical solution gives satisfactory results for lower $d'/\sigma$ ratio (e. g. the two examples at $t = 3.5$ days on Fig. 4).

## 4. Dependence on parameters

### 4.1. Time required to reach the equilibrium of the water-table

Fig. 5 shows examples of hydraulic head computations with an underground barrier alone (Fig. 5a, c) and, with an underground barrier and a well pumping at a rate of 300 m$^3$/h located at 300 m upgradient of the barrier (Fig. 5b, d). For the four given situations, the geometry of the underground barrier and the hydraulic parameters of the aquifer are identical (Lx=5 m, Ly=800 m, $i = 0.005$, $D_0=25$ m, $K = 8 \cdot 10^{-4}$ m/s, $S = 0.35$), only differs the resistance parameter of the barrier. It is impermeable on Figs. 5a, b (Lx·K/(Ly·k')=∞), and semi-permeable with Lx·K/(Ly·k')=1.0 on Figs. 5c, d (i.e. $k'=4.8 \cdot 10^{-6}$ m/s). The calculations were carried out at different times since the construction of the structure, from 10 min (~7·10$^{-3}$ day) up to $10^6$ min (~695 days). Calculations with steady-state solutions are also presented. In these examples, the well is slightly offset from the x-axis, at 0.1 m (i.e. $y_w$: 0.1 m) to avoid

infinite drawdown at the well location (see the pumping well term in Eqs. (14) and 16). At the early time, the curves (up to 100 min) do not show significant difference according to initial conditions without or with a pumping well. Then, the influence of the underground barrier appears on the curves and show the classical increase of the water-table upgradient and decrease downgradient of the structure. As expected, the variations in hydraulic head compared to initial condition are less pronounced in the semi-permeable barrier situation, that is, it increases less upgradient and decreases less downgradient of the barrier compared to the impermeable barrier situation. This is due to the groundwater flow through the structure for the semi-permeable barrier situation, which reduces the obstruction of the structure and thus the increase and decrease of the hydraulic head up- and downgradient of the latter. Concerning the time to reach the steady-state condition, for these examples about two years were required to stabilize the water-table.

Fig. 6a presents the dimensionless time $t \cdot D_0K/(S \cdot Ly^2)$ vs. the ratio $\Delta h_t/\Delta h_{t=\infty}$, in case of an impermeable underground barrier (Lx·K/(Ly·k')=∞). The graph has been drawn with cases with and without a pumping well (Table 1; 24 cases with different barrier geometries, aquifer properties and pumping). Pumping wells were located close to the x-axis upgradient of the barrier ($0.1 \le y_w \le 0.3$ m), where the structure creates the most increase in hydraulic head. $\Delta h_t$ is the maximum hydraulic-head variation at the barrier sides, i.e. at the centre of the barrier, in transient (Eq. (14)), and $\Delta h_{t=\infty}$ is the maximum hydraulic head variation at the barrier sides (i.e. centre of the barrier) in steady-state condition (Eq. (16)); Fig. 6a. Therefore, when the ratio tends to 1, it means that the time has been long enough to reach water-level equilibrium (i.e. steady-state condition). The curve has the form of a cumulative log-normal distribution function (see Appendix D) and shows that the more the aquifer is diffusive (increase of $D_0K/S$) and the
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Examples of hydraulic head computations with Eq. (14) (transient) and Eq. (16) (steady-state) along the x-axis. a) $Lx=5$ m, $Ly=800$ m, $i = 0.005$, $D_0=25$ m, $K = 8 \cdot 10^{-4}$ m/s, $S = 0.35$, $Lx \cdot K/(Ly \cdot k') = \infty$ (impermeable barrier), $10 \le t \le 10^6$ min (Eq. (14)) and steady-state (Eq. (16)), $d'=2.3$ m and $n = 1300$. b) as a) but with a pumping well $x_w=-300$ m and $y_w=0.1$ m, $Q_p=300$ m³/h. c) as a) but with $Lx \cdot K/(Ly \cdot k')=1.0$ (semi-permeable barrier, $k'=4.8 \cdot 10^{-7}$ m/s). d) as c) but with a pumping well $x_w=-300$ m and $y_w=0.1$ m, $Q_p=300$ m³/h.</figcaption></figure>

less the length of the barrier perpendicular to the regional flow (Ly) is, the quicker the equilibrium is reached.

The analysis has been also performed for different hydraulic resistance parameters of the barrier ($Lx \cdot K/(Ly \cdot k')=0.3, 1.0$ and $10.0$, Fig. 6b, Table 2). The curves for each resistance parameter value are always in the form of a cumulative log-normal distribution function, but the time to reach equilibrium logically decreases as the resistance parameter decreases (i.e. increase in hydraulic conductivity of the barrier). It is also noticed that when the resistance parameter of the barrier approaches a value close to 10 ($Lx \cdot K/(Ly \cdot k')=10.0$; Fig. 6b) the curve is almost identical to the case of an impermeable barrier. This shows that it is not necessary to build a strictly impermeable structure to achieve a similar impact on the water-table, particularly on the rise of the water-table upgradient of the barrier. These curves are therefore useful for evaluating the time required to reach water-table equilibrium after the construction of the underground barrier, and for different resistance parameters of the barrier. The type-curves are provided in Appendix D.

## 4.2. Optimal location of the pumping well and optimal pumping rate according to the geometry and the location of the underground barrier

Another interesting point highlighted by the situations with a pumping well (Fig. 5b and d) is the decrease of the water-table downgradient of the underground barrier compared to the situation with no pumping well, due to the propagation of the drawdown induced by the well further than the barrier. The question then arises of the optimal location of the pumping well in relation to the barrier, i.e. where the increase of water table at the location of the well due to the barrier is maximum. If the well is too close to the barrier, the water level can be lower than before the construction of the structure because of the too strong influence exerted by the no flow or semi-permeable barrier.

Conversely, if the pumping well is located too upgradient of the underground barrier, the well can be outside the groundwater mound created by the barrier, or too far to have a significant, or desired, increase in hydraulic head at the location of the well.

Fig. 7 presents example of the computation of the hydraulic-head in steady-state condition (Eq. (16)) for an impermeable underground barrier of fixed geometry and fixed aquifer hydraulic properties, where the well is placed upgradient of the structure at increasing distance (from 50 to 1200 m of the barrier; $y_w$: 0.15 m). Fig. 7a shows initial condition (pumping without barrier), and Fig. 7b the computed steady-state hydraulic-head with the barrier (Eq. (16)). The situation at 50 m upgradient of the barrier induced an increase of the water-table of about 0.8 m at the well location, while at 1200 m the increase is less than 0.1 m. Fig. 8 analyses the variation of the hydraulic-head at the location of the well ($\Delta h_{wd}$) for the situation described in Fig. 7, and for three other pumping rate cases (100, 200 and 600 m³/h), according to the dimensionless distance of the well to the barrier ($x_w/Ly$). The distance is made dimensionless according to the length of the barrier Ly because it is this length which induces the increase in the water table upgradient of the barrier. The figure shows that if the pumping rate is 100 m³/h, the optimal location of the well is about 25 m upgradient of the barrier with a maximum increase of the water-table $\Delta h_{wd\_max}$ of about 2.3 m, at 200 m³/h it is about 50 m with $\Delta h_{wd\_max}$: 1.76 m, at 400 m³/h it is about 120 m with $\Delta h_{wd\_max}$: 1.13 m and at 600 m³/h it is about 170 m with $\Delta h_{wd\_max}$: 0.79 m. Thus, the optimal location of the well is non-linear, and it depends not only on the barrier geometry but also on the pumping rate of the well, and probably on the aquifer properties. The experience described in Figs. 7 and 8 was carried out for nine other underground barrier (impermeable) and aquifer situations with different pumping rates (more than 400 simulations, Table 3). Figs. 9a and b present the results with dimensionless parameters. Fig. 9a presents
<figure><img src="image_8.png" /><figcaption>Fig. 6. Type curves describing the time (t·D₀K/(S·Ly²), dimensionless) required to reach the equilibrium of the water-table (ratio Δhₜ/Δhₜ=∞). a) impermeable underground barrier. Constructed from 24 cases of hydraulic head modelling from Eq. (14) (transient) and Eq. (16) (steady-state), cases 1 to 24 Table 1. b) semi-permeable underground barrier with Lx·K/(Ly·k')=0.3, 1.0 and 10.0. Constructed from 18 cases of hydraulic head modelling from Eq. (14) (transient) and Eq. (16) (steady-state), cases A to F Table 2. Δhₜ: maximum hydraulic-head variation at the barrier sides (centre of the barrier) in transient, Δhₜ=∞: maximum hydraulic-head variation at the barrier sides (centre of the barrier) in steady-state condition. For the cases with pumping wells, yₗ was set between 0.1 and 0.3 m.</figcaption></figure>

the dimensionless optimal distance of the well to the barrier (x<sub>w_opt</sub>/Ly) as function of the dimensionless pumping rate (Q<sub>p</sub>/[D<sub>0</sub>·K·i·Ly]); x<sub>w_opt</sub> is the optimal distance of the well to the barrier where the increase of water level due to the barrier is maximal. Q<sub>p</sub> is made dimensionless according to the transmissivity of the aquifer (D<sub>0</sub>K) and the product i·Ly, which relates the efficiency of the barrier to increase the water-table upgradient of the barrier. The data follow a linear curve on a log-log scale (see the given equation on the graph). Fig. 9b presents the inverse of the dimensionless efficiency of the barrier on the well (Δh<sub>wd_max</sub>·D<sub>0</sub>·K/Q<sub>p</sub>) as function of the dimensionless optimal distance of the well to the barrier (x<sub>w_opt</sub>/Ly). The dimensionless efficiency relates that the maximum water-table increase (Δh<sub>wd_max</sub>) depends on the pumping rate (Q<sub>p</sub>) and the aquifer transmissivity. Similarly, the data follow a linear curve on a log-log scale (see the given equation on the graph).

The use of the two graphs allows, for example, by knowing the
properties of the aquifer (K, D₀, i) and the pumping rate of the well (Qₚ)
to evaluate the length (Ly) of the impermeable underground barrier

perpendicular to the natural gradient i and, at which distance from the
well it should be constructed (x<sub>w_opt</sub>) to have the maximum increase in
the hydraulic-head at the location of the well. In this case, from the two
empirical relationships defined on Fig.9a, b, the following relationship
for the length of the barrier (Ly) perpendicular to the regional ground-
water flow is:

$$
\text{Ly} = \left( \frac{\Delta h_{\text{max}}}{0.1581 \cdot i^{1.857} \left( \frac{T}{Q_p} \right)^{0.857}} \right)^{0.5385} \quad (19)
$$

Where T is the aquifer transmissivity (T = D₀K, in m²/s).

Then the optimal location of the barrier x<sub>w_opt</sub> is given by (Fig. 9a):

$$
x_{w\_opt} = 0.1892 \cdot \text{Ly}^{-0.241} \cdot \left( \frac{Q_p}{T \cdot i} \right)^{1.241} \quad (20)
$$

It is very important to note that the proposed empirical relationships (Eqs. (19), (20)) are valid under the conditions of which they were established, i.e. 0.03≤x<sub>w_opt</sub>/Ly≤3, 0.2≤Q<sub>p</sub>/L<sub>y·T·i</sub>≤11, 0.001≤Δh<sub>max</sub>T/Q<sub>p</sub>≤1, and for the radius of the pumping well between 0.1 and 0.3 m (0.1≤y<sub>w</sub>≤0.3).

Note that the two curves do not depend on the width of the under-
ground barrier parallel to the natural gradient (Lx) because it is the
length perpendicular to the natural gradient (Ly) which is related to the
increase in the water table upgradient of the barrier. However, Lx in-
tervenes in the hydraulic resistance parameter of the barrier. This is
discussed later.

4.3. Which hydraulic conductivity for the underground barrier?

Fig. 10 shows the dependence of the ratio between the maximum hydraulic-head variation at the underground barrier sides, i.e. at the centre of the barrier, for a semi-permeable barrier (Δh<sub>Lx,K/(Ly,k')</sub>) and, the one with an impermeable barrier (Δh<sub>Lx,K/(Ly,k')=∞</sub>), according to the resistance parameter of the structure (Lx·K/[Ly·k']). The graph was established with simulations in steady-state condition (Eq. (16)) for different barrier geometries, aquifer properties, pumping and barrier hydraulic conductivities (more than 100 simulations, Table 4). The form of the curve follows a cumulative log-normal distribution function (see Appendix D) and precises the efficiency of the underground barrier according to its resistance parameter. It shows that for resistance parameter close to 10.0, the variation in hydraulic head is 10 % less compared to the case of an impermeable barrier. For a resistance parameter for the barrier of about 20.0, the variation in hydraulic head is 5 % less compared to the situation of an impermeable structure. Finally, a difference in hydraulic-head close to 0 (i.e. ratio (Δh<sub>Lx,K/(Ly,k')</sub>/(Δh<sub>Lx,K/(Ly,k')=∞</sub>~1) is reached for resistance parameter values close to 100.0. The drawn type-curve is useful to evaluate the maximum value of hydraulic conductivity of the material that must be used to construct the underground barrier.

5. Application to examples

5.1. Time required to reach water-table equilibrium through a field example

Fig. 11 presents a field example given by Pujades et al. (2012). It corresponds to the basement of the *El Prat de Llobregat* railway station (Spain), where 29-m deep and impermeable diaphragm walls were used to make an enclosure, completely penetrating an aquifer formed by sand and gravel deposits. Its thickness ranges from 12 to 20 m, with hydraulic conductivity ranging from 1 to 5·10⁻⁴ m/s, and storativity ranging between 0.05 and 0.3 (Tubeau, 2004). Therefore, the average diffusivity of the aquifer is about 3·10⁻² m²/s (min. value: 4·10⁻³ m²/s and max. value: 2·10⁻¹ m²/s). The basement of the railway station can be
<caption>**Table 1**</caption>

Impermeable underground barrier, pumping well and aquifer characteristics used for establishing Fig. 6a. For each case hydraulic head was computed for $t = 10, 100, 1000, 10^4, 10^5, 10^6$ min (Eqs.14), and in steady-state condition (Eq.16). For the cases with a pumping well, $y_w$ was set between 0.1 and 0.3 m to represent realistic well radius.

<table><caption>**Table 2**</caption><thead><tr><th>Case</th><th>K (m/s)</th><th>D<sub>0</sub> (m)</th><th>S (-)</th><th>i (-)</th><th>Ly (m)</th><th>Lx (m)</th><th>T/S (m<sup>2</sup>/s)</th><th>Ly/Lx (-)</th><th>Q<sub>p</sub> (m<sup>3</sup>/h)</th><th>x<sub>w</sub> (m)</th></tr></thead><tbody><tr><td>Case 1</td><td>1.0E-03</td><td>15</td><td>0.01</td><td>0.005</td><td>500</td><td>5</td><td>1.50E+00</td><td>100</td><td>-</td><td>-</td></tr><tr><td>Case 2</td><td>1.0E-02</td><td>15</td><td>0.05</td><td>0.01</td><td>800</td><td>4</td><td>3.00E+00</td><td>200</td><td>-</td><td>-</td></tr><tr><td>Case 3</td><td>1.0E-03</td><td>15</td><td>0.1</td><td>0.001</td><td>800</td><td>4</td><td>1.50E-01</td><td>200</td><td>-</td><td>-</td></tr><tr><td>Case 4</td><td>1.0E-03</td><td>15</td><td>0.35</td><td>0.001</td><td>800</td><td>4</td><td>4.29E-02</td><td>200</td><td>-</td><td>-</td></tr><tr><td>Case 5</td><td>1.0E-03</td><td>50</td><td>0.35</td><td>0.01</td><td>800</td><td>4</td><td>1.43E-01</td><td>200</td><td>-</td><td>-</td></tr><tr><td>Case 6</td><td>1.0E-03</td><td>50</td><td>0.35</td><td>0.01</td><td>100</td><td>2</td><td>1.43E-01</td><td>50</td><td>-</td><td>-</td></tr><tr><td>Case 7</td><td>1.0E-03</td><td>15</td><td>0.1</td><td>0.001</td><td>100</td><td>2</td><td>1.50E-01</td><td>50</td><td>-</td><td>-</td></tr><tr><td>Case 8</td><td>1.0E-04</td><td>10</td><td>0.3</td><td>0.008</td><td>600</td><td>2</td><td>3.33E-03</td><td>300</td><td>-</td><td>-</td></tr><tr><td>Case 9</td><td>1.5E-02</td><td>40</td><td>0.03</td><td>0.018</td><td>600</td><td>8</td><td>2.00E+01</td><td>75</td><td>-</td><td>-</td></tr><tr><td>Case 10</td><td>1.5E-02</td><td>40</td><td>0.5</td><td>0.018</td><td>600</td><td>8</td><td>1.20E+00</td><td>75</td><td>-</td><td>-</td></tr><tr><td>Case 11</td><td>1.5E-04</td><td>40</td><td>0.05</td><td>0.002</td><td>600</td><td>8</td><td>1.20E-01</td><td>75</td><td>-</td><td>-</td></tr><tr><td>Case 12</td><td>1.5E-03</td><td>20</td><td>0.25</td><td>0.02</td><td>100</td><td>10</td><td>1.20E-01</td><td>10</td><td>-</td><td>-</td></tr><tr><td>Case 13</td><td>5.0E-04</td><td>20</td><td>0.25</td><td>0.006</td><td>450</td><td>3</td><td>4.00E-02</td><td>150</td><td>-</td><td>-</td></tr><tr><td>Case 14</td><td>1.0E-04</td><td>4</td><td>0.3</td><td>0.02</td><td>15</td><td>5</td><td>1.33E-03</td><td>3</td><td>-</td><td>-</td></tr><tr><td>Case 15</td><td>1.0E-04</td><td>4</td><td>0.3</td><td>0.02</td><td>30</td><td>5</td><td>1.33E-03</td><td>6</td><td>-</td><td>-</td></tr><tr><td>Case 16</td><td>1.0E-04</td><td>50</td><td>0.3</td><td>0.0025</td><td>2500</td><td>10</td><td>1.67E-02</td><td>250</td><td>-</td><td>-</td></tr><tr><td>Case 17</td><td>1.0E-03</td><td>10</td><td>0.3</td><td>0.008</td><td>200</td><td>2</td><td>3.33E-02</td><td>100</td><td>200</td><td>-300</td></tr><tr><td>Case 18</td><td>1.0E-03</td><td>15</td><td>0.25</td><td>0.014</td><td>750</td><td>2</td><td>6.00E-02</td><td>375</td><td>600</td><td>-150</td></tr><tr><td>Case 19</td><td>1.0E-05</td><td>20</td><td>0.35</td><td>0.004</td><td>250</td><td>2</td><td>5.71E-04</td><td>125</td><td>6</td><td>-240</td></tr><tr><td>Case 20</td><td>1.0E-03</td><td>15</td><td>0.06</td><td>0.006</td><td>800</td><td>4</td><td>2.50E-01</td><td>200</td><td>600</td><td>-180</td></tr><tr><td>Case 21</td><td>1.0E-04</td><td>50</td><td>0.3</td><td>0.003</td><td>2500</td><td>10</td><td>1.67E-02</td><td>250</td><td>600</td><td>-180</td></tr><tr><td>Case 22</td><td>5.0E-04</td><td>50</td><td>0.3</td><td>0.006</td><td>2000</td><td>7</td><td>8.33E-02</td><td>285.7</td><td>600</td><td>-180</td></tr><tr><td>Case 23</td><td>2.0E-03</td><td>10</td><td>0.35</td><td>0.01</td><td>800</td><td>5</td><td>5.71E-02</td><td>160.0</td><td>-</td><td>-</td></tr><tr><td>Case 24</td><td>2.0E-03</td><td>10</td><td>0.35</td><td>0.01</td><td>800</td><td>5</td><td>5.71E-02</td><td>160.0</td><td>300</td><td>-300</td></tr></tbody></table>

Semi-permeable underground barrier, pumping well and aquifer characteristics used for establishing Fig. 6b. For each case hydraulic head was computed for $t = 10, 100, 1000, 10,000, 10^5$ and $10^6$ min (Eqs.14), and in steady-state condition (Eq.16), and with $Lx \cdot K/(Ly \cdot k')=0.3, 1.0$ and $10.0$. For the cases with a pumping well, $y_w$ was set between 0.1 and 0.3 m.

<table><thead><tr><th>Case</th><th>K (m/s)</th><th>D<sub>0</sub> (m)</th><th>S (-)</th><th>i (-)</th><th>Ly (m)</th><th>Lx (m)</th><th>T/S (m<sup>2</sup>/s)</th><th>Ly/Lx (-)</th><th>Q<sub>p</sub> (m<sup>3</sup>/h)</th><th>x<sub>w</sub> (m)</th></tr></thead><tbody><tr><td>Case A</td><td>5.0E-04</td><td>50</td><td>0.3</td><td>0.006</td><td>100</td><td>4</td><td>8.3E-02</td><td>25.0</td><td>-</td><td>-</td></tr><tr><td>Case B</td><td>1.0E-03</td><td>50</td><td>0.15</td><td>0.015</td><td>200</td><td>3</td><td>3.3E-01</td><td>66.7</td><td>-</td><td>-</td></tr><tr><td>Case C</td><td>2.0E-04</td><td>30</td><td>0.05</td><td>0.003</td><td>1000</td><td>5</td><td>1.2E-01</td><td>200.0</td><td>-</td><td>-</td></tr><tr><td>Case D</td><td>2.0E-03</td><td>10</td><td>0.35</td><td>0.010</td><td>800</td><td>20</td><td>5.7E-02</td><td>40.0</td><td>-</td><td>-</td></tr><tr><td>Case E</td><td>8.0E-04</td><td>15</td><td>0.2</td><td>0.008</td><td>400</td><td>3</td><td>6.0E-02</td><td>133.3</td><td>200.0</td><td>-150</td></tr><tr><td>Case F</td><td>1.5E-03</td><td>25</td><td>0.15</td><td>0.005</td><td>600</td><td>3</td><td>2.5E-01</td><td>200.0</td><td>-</td><td>-</td></tr></tbody></table>

assimilated to an impermeable underground physical barrier that fully penetrates the aquifer. Its length and width are about 1800 m (Ly) and 200 m (Lx). The hydraulic gradient is $6.6 \cdot 10^{-4}$, and the time-variations of the hydraulic-head from two observation wells, located up- and downgradient near the basement of the railway station, are given in Fig. 11a. The construction of the basement of the railway station stopped in November 2012 and water level increased upgradient and decreased downgradient approximately up to August 2008 with a maximum difference of water level between the two wells of about 1.25 to 1.3 m, thus about 290 days after the stop of the construction. Afterwards, the water levels seem to have been disturbed by other processes (e.g. recharge of the aquifer, pumping, etc.). Fig. 11b presents the type curves given in Figs. 6a, b and the variation in hydraulic head between the two observation wells according to the dimensionless time ($t \cdot D_0K/(S \cdot Ly^2)$). To match the observed data with the type-curve of an impermeable barrier (black curve), a diffusivity value of $4.2 \cdot 10^{-2}$ m$^2$/s, and a maximum difference in hydraulic-head in steady-state condition of 1.36 m ($\Delta h_{t-\infty}$) were necessary. The fit is relatively satisfactory ($R^2=0.80$, RMSE D [standardized root mean square deviation]=0.036). The diffusivity value estimated is very close to field data, and the value for $\Delta h_{t-\infty}$ (steady-state condition) relatively close to that measured 290 days after stop of the construction, which should correspond to about 92–96 % of the maximum expected difference in hydraulic-head. According to the type curve, the time required to reach the equilibrium of the water-table (steady-state condition) must have taken several additional years. Considering $t \cdot D_0K/(S \cdot Ly^2) = 10$, which corresponds to water-table

equilibrium, it would have taken more than 20 years to reach steady-state condition. However, if we consider 99 % of equilibrium ($t \cdot D_0K/(S \cdot Ly^2) = 3$), 7 to 8 years would have been necessary to reach a situation close to steady-state condition.

## 5.2. Design of a hypothetical underground barrier to raise the water-table at the location of a pumping well

An example of the use of the proposed type-curves is given through the following practice problem. In an alluvial aquifer, a well is pumped at a flow rate of 250 m$^3$/h. Because of the decline of the water-table during summer periods, the construction of an underground physical barrier is planned to increase the water level at the location of the pumping well by 1 m. The aquifer has a thickness of 10 m (ranging between 8 and 12 m), its average transmissivity is $2.5 \cdot 10^{-2}$ m$^2$/s (2.2 - $2.8 \cdot 10^{-2}$ m$^2$/s), its average storativity is 5 % (0.04 - 0.06) and the average hydraulic gradient 0.0035 (0.0034 - 0.0036). How long should the underground barrier be (Ly)? How far should it be from the well ($x_{w_opt}$)? After how long will the water-table be in equilibrium? And what width (Lx), what hydraulic conductivity ($k'$) for the barrier? The results will be given with a range of values allowing the uncertainty in the hydrogeological parameters to be considered.

Therefore, with $T = 2.5 \cdot 10^{-2}$ m$^2$/s, $i = 0.0035$ and $Q_p=250$ m$^3$/h, the length of the barrier Ly is evaluated from Eq. (19) to about 1240 m. The uncertainty on Ly due to the variability of hydrogeological parameters (i.e. $T$ and $i$) is about $\pm 100$ m. Then, the optimal location of the barrier
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Steady-state hydraulic head in the case of an impermeable underground barrier and a well pumping at 400 m³/h located at different distances, $x_w$, upgradient of the barrier (Eq. (16)). a) initial condition, with no barrier. b) with the barrier. The geometry of the barrier and the properties of the aquifer are given on the figures. $y_w$: 0.15 m.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Analysis of the variation of the hydraulic head at the location of the pumping well induced by the construction of the underground barrier ($\Delta h_{wd}$, in m) for different dimensionless distance of the well upgradient of the barrier ($x_w/Ly$), steady-state condition, and different pumping rates ($Q_p$). $\Delta h_{wd\_max}$: maximum increase of the water-table at the well location (in m). The geometry of the underground barrier and the properties of the aquifer are given on Fig. 7.</figcaption></figure>
<caption>**Table 3**</caption>

Impermeable underground barrier, pumping well and aquifer characteristics used for establishing Fig. 9a,b. Steady-state condition (Eq. (16)). For the pumping wells, $y_w$ was set between 0.1 and 0.3 m.

<table><thead><tr><th>Case</th><th>K (m/s)</th><th>D<sub>0</sub> (m)</th><th>i (-)</th><th>Ly (m)</th><th>Lx (m)</th><th>Ly/Lx (-)</th><th>Q<sub>p</sub> (m<sup>3</sup>/h)</th></tr></thead><tbody><tr><td>Case p-1</td><td>3.0E-04</td><td>10</td><td>0.005</td><td>600</td><td>4</td><td>150.0</td><td>20, 35, 50, 60</td></tr><tr><td>Case p-2</td><td>1.0E-03</td><td>20</td><td>0.003</td><td>810</td><td>3</td><td>270.0</td><td>125, 250, 500</td></tr><tr><td>Case p-3</td><td>5.0E-04</td><td>25</td><td>0.015</td><td>150</td><td>3</td><td>50.0</td><td>50, 100, 200, 400</td></tr><tr><td>Case p-4</td><td>5.0E-04</td><td>25</td><td>0.003</td><td>510</td><td>3</td><td>170.0</td><td>50, 100, 200, 400</td></tr><tr><td>Case p-5</td><td>5.0E-04</td><td>25</td><td>0.015</td><td>510</td><td>3</td><td>170.0</td><td>100, 200, 400, 600</td></tr><tr><td>Case p-6</td><td>2.5E-03</td><td>10</td><td>0.0035</td><td>470</td><td>3</td><td>156.7</td><td>100, 200, 400, 600</td></tr><tr><td>Case p-7</td><td>1.0E-04</td><td>50</td><td>0.0080</td><td>220</td><td>2</td><td>110.0</td><td>50, 100, 200, 400</td></tr><tr><td>Case p-8</td><td>3.0E-04</td><td>20</td><td>0.0020</td><td>300</td><td>3</td><td>100.0</td><td>25, 40, 60, 80, 110, 130</td></tr><tr><td>Case p-9</td><td>5.0E-06</td><td>100</td><td>0.01</td><td>950</td><td>3</td><td>316.7</td><td>10, 25, 50, 75</td></tr><tr><td>Case p-10</td><td>1.0E-04</td><td>50</td><td>0.0025</td><td>2500</td><td>10</td><td>250.0</td><td>50, 100, 200, 250</td></tr></tbody></table>

$x_{w, opt}$ is given by Eq. (20) and is evaluated approximately 135 m downgradient of the well. The uncertainty of the well location due to the variability of hydrogeological parameters is about ±20 m.

The time necessary to reach the equilibrium of the water-table is given by the type-curve on Fig. 6a (impermeable barrier). With a storativity $S$ of 0.05, the time necessary to reach the equilibrium $(t·T/(S·Ly^2) \sim 10.0)$ is evaluated at 362 days, thus about one year. However, 95 % $(t·T/(S·Ly^2) \sim 0.7)$ of the equilibrium should be reached in less than a month (25 days), and 90 % $(t·T/(S·Ly^2) \sim 0.31)$ in 11 days. Considering the uncertainty due to the variability of hydrogeological parameters $(T, S$ and $i)$ and the evaluation of ranges for the length of the barrier (Ly), the time necessary to reach the equilibrium varies between 330 and 620 days, for 95 % of the equilibrium between 23 and 43 days, and between 10 and 19 days for 90 %. These different time estimates, quite significant, are mainly due to the uncertainty of the aquifer diffusivity $(0.4–0.7 m^2/s)$.

The hydraulic conductivity $k'$ of the material, which should be used to construct the underground barrier, can be estimated with the used of the type-curve in Fig. 10. Fig. 12 presents the range of $k'$ for different degree of efficiency of the barrier and width $Lx$: 90 % $(Lx·K/(Ly·k')=10.0)$, 95 % $(Lx·K/(Ly·k')=20.0)$, 99 % $(Lx·K/(Ly·k')=60.0)$ and 100 % $(Lx·K/(Ly·k')=100.0)$, for the given example, i.e. with $K = 2.5·10^{-3}$ m/s and $Ly=1240$ m. The graph shows that for a barrier width (Lx) between 3 and 4 m, the hydraulic conductivity range of the barrier material must be between 1 and $4·10^{-7}$ m/s to expect a barrier efficiency of 95–99 %. Considering the uncertainty due to the variability of hydrogeological parameters $(T, S$ and $i)$ and the ranges evaluated for the length of the barrier (Ly), the minimum hydraulic conductivity for the barrier with an efficiency of 95–99 % is about $8·10^{-8}$ m/s, and the maximum about $5·10^{-7}$ m/s, thus very close to the previous estimate.

## 6. Discussion

Above, transient and steady-state semi-analytical solutions based on MFS were developed for evaluating the rise and decline of hydraulic head in an aquifer caused by the implementation of a semi-permeable underground physical barrier, with or without a well pumping in the aquifer (Eqs. (14) and (16)). Comparison of the developed transient solution with numerical modelling showed very satisfactory result and conform previous comparison under steady-state condition (Dewandel et al., 2024a). This reinforces the use of the MFS to provide solutions in

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

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Dimensionless curves of the optimal location of a pumping well near an impermeable underground barrier. Constructed from simulations with Eq. (16), with underground barrier and aquifer situations, and pumping rates from Table 3 (cases p-1 to p-10). a) dimensionless optimal distance of the well $(x_{w, opt}/Ly)$ vs. dimensionless pumping rate $(Q_p/[D_0·K·i·Ly])$. b) Inverse of the dimensionless efficiency of the barrier on the well vs. dimensionless optimal distance of the well $(\Delta h_{wd_max}·D_0·K/Q_p)$. R: linear regression coefficient. Note that the proposed empirical relationships are valid under the conditions of which they were established (see text).</figcaption></figure>

complex system, even under transient condition (see also the references given in the Introduction). This work complements previous analytical works aimed at assessing the hydrogeological impact of the construction of hydraulic barriers on aquifers (Anwar, 1983; Pujades et al., 2012; Wu et al., 2023; Dewandel et al., 2024a, 2024b). Since the developed solutions assume a homogeneous aquifer, they should be used in situations where the aquifer can be considered relatively homogeneous, thus without significant spatial variations in hydraulic conductivity and/or storage coefficient, as is often the case in many alluvial aquifers (e.g., the example given in Section 5.1). Also, since the developed solutions neglect the storage of the barrier, its use under transient condition is valid for physical barriers constructed with materials with very low drainable porosity (clay, concrete, etc.), which is assumed to be reasonable for most practical cases.

The proposed solutions 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 (first method of linearization or h-linearization) is commonly
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Efficiency of the underground barrier ($\Delta h_{Lx.K/(Ly.k')}/\Delta h_{Lx.K/(Ly.k')-\infty}$) according to its hydraulic resistance parameter ($Lx.K/[Ly.k']$), steady-state condition. Constructed from simulations with Eq. (16), with barrier and aquifer situations, and pumping rates from Table 4 (cases s-1 to s-11). a) Sketch. b) Dimensionless graph. $\Delta h_{Lx.K/(Ly.k')}$: maximum hydraulic-head variation at the barrier sides (centre of the barrier) for a semi-permeable barrier, $\Delta h_{Lx.K/(Ly.k')=\infty}$: maximum hydraulic-head variation at the barrier sides (centre of the barrier) for an impermeable barrier.</figcaption></figure>

<caption>**Table 4**</caption>

Impermeable and semi-permeable barriers, pumping well and aquifer characteristics used for establishing Fig. 10. nb: number of simulations with different hydraulic conductivity values for the barrier for each case. For the pumping wells, $y_w$ was set between 0.1 and 0.3 m.

<table><thead><tr><th>Case</th><th>K (m/s)</th><th>D<sub>0</sub> (m)</th><th>i (-)</th><th>Ly (m)</th><th>Lx (m)</th><th>Ly/Lx (-)</th><th>k' (m/s)</th><th>Lx.K/(Ly.k') (-)</th><th>Q<sub>p</sub> (m<sup>3</sup>/h)</th><th>x<sub>w</sub> (m)</th><th>nb</th></tr></thead><tbody><tr><td>Case s-1</td><td>1.0E-03</td><td>50.0</td><td>0.006</td><td>660</td><td>3</td><td>220.0</td><td>0-3.3E-5</td><td>0.1-&gt;∞</td><td>-</td><td>-</td><td>8</td></tr><tr><td>Case s-2</td><td>1.0E-04</td><td>20.0</td><td>0.010</td><td>200</td><td>3</td><td>66.7</td><td>0-2.0E-6</td><td>0.8-&gt;∞</td><td>-</td><td>-</td><td>8</td></tr><tr><td>Case s-3</td><td>1.0E-02</td><td>50.0</td><td>0.003</td><td>500</td><td>10</td><td>50.0</td><td>0-1.0E-3</td><td>0.1-&gt;∞</td><td>-</td><td>-</td><td>8</td></tr><tr><td>Case s-4</td><td>5.0E-04</td><td>70.0</td><td>0.014</td><td>750</td><td>2</td><td>375.0</td><td>0-1.2E-5</td><td>0.2-&gt;∞</td><td>-</td><td>-</td><td>8</td></tr><tr><td>Case s-5</td><td>2.0E-03</td><td>15.0</td><td>0.002</td><td>800</td><td>15</td><td>53.3</td><td>0-2.0E-4</td><td>0.1-&gt;∞</td><td>-</td><td>-</td><td>12</td></tr><tr><td>Case s-6</td><td>1.0E-04</td><td>10.0</td><td>0.020</td><td>150</td><td>5</td><td>30.0</td><td>0-2.0E-5</td><td>0.2-&gt;∞</td><td>-</td><td>-</td><td>14</td></tr><tr><td>Case s-7</td><td>2.0E-05</td><td>100.0</td><td>0.008</td><td>1000</td><td>5</td><td>200.0</td><td>0-3.3E-7</td><td>0.3-&gt;∞</td><td>-</td><td>-</td><td>7</td></tr><tr><td>Case s-8</td><td>2.0E-05</td><td>100.0</td><td>0.008</td><td>1000</td><td>5</td><td>200.0</td><td>0-3.3E-7</td><td>0.3-&gt;∞</td><td>200</td><td>-890</td><td>7</td></tr><tr><td>Case s-9</td><td>1.0E-04</td><td>10.0</td><td>0.020</td><td>150</td><td>5</td><td>30.0</td><td>0-2.0E-5</td><td>0.3-&gt;∞</td><td>25</td><td>25</td><td>14</td></tr><tr><td>Case s-10</td><td>2.0E-03</td><td>15.0</td><td>0.002</td><td>800</td><td>15</td><td>53.3</td><td>0-2.0E-4</td><td>0.2-&gt;∞</td><td>400</td><td>-430</td><td>12</td></tr><tr><td>Case s-11</td><td>1.0E-03</td><td>10.0</td><td>0.006</td><td>660</td><td>3</td><td>220.0</td><td>0-3.3E-5</td><td>0.1-&gt;∞</td><td>150</td><td>-130</td><td>8</td></tr></tbody></table>

used for solving groundwater hydraulic problems and was chosen to solve the standard transient flow equation (Eq. (1)). According to previous works (de Marsily, 1986; Moutsopoulos 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), assumed that relative differences are around a few percent if $\beta \le 0.5$. Works on the hydraulic-head distribution induced by the construction of underground barriers and dams on river, with another linearization method which gives a solution closer to the exact solution ($h^2$-linearization Dewandel et al., 2024b, 2024c) showed that the solution with the first linearization method gives very satisfactory results when the maximum degree of saturation/desaturation $\beta$ is less

than 20 %. However, the computation times with the $h^2$-linearization method are extremely long, which makes this solution unsuitable for practical use. In all simulations (about 180) used to build the different type curves (Figs. 6, 9, 10, 13) $\beta$ ranges between $6 \cdot 10^{-4}$ and 0.17 (average value: 0.05). For higher $\beta$ values, results are approximate but can be used as a preliminary diagnosis.

One of the questions that may be asked when considering the construction of an underground barrier is: how long is needed to stabilize the water-table? The type-curve in Fig. 6a allows the duration at different stages to be estimated. The example given in Fig. 11 (basement of a railway station) shows that field data fits satisfactorily the curve of an impermeable underground barrier. However, the curve also shows that in this case, due to the great length of the basement (1800 m) and the diffusivity of the aquifer ($4.2 \cdot 10^{-2} \text{ m}^2/\text{s}$), it probably took twenty
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Field application of the type curves of the time required to reach the equilibrium of the water-table, basement of the *El Prat de Llobregat* railway station (Spain) from Pujades et al., (2012). a) Water level measurements up- and downgradient of the basement (modified from Pujades et al., 2012). b) Analysis of the difference in hydraulic-head between the two observation wells (up- and downgradient) after the stop of the construction compared to theoretical time type curves (Fig. 6), diffusivity of the aquifer: $4.2 \cdot 10^{-2}$ m$^2$/s, $\Delta h_t=1.36$ m. The horizontal bars correspond to the uncertainty on time (±10 days) and, the vertical bars to the uncertainty on aquifer diffusivity (see text) and on $\Delta h_t$ values (±0.02 to 0.05 m). R: linear regression coefficient, RMSE D: standardized root mean square deviation.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Evaluation of the hydraulic conductivity of the underground barrier $k'$ for different barrier efficiencies (compared to impermeable barrier), 90 % ($Lx·K/(Ly·k')=10.0$), 95 % ($Lx·K/(Ly·k')=20.0$), 99 % ($Lx·K/(Ly·k')=60.0$) and 100 % ($Lx·K/(Ly·k')=100.0$). Deduced from type curve of Fig. 10, with $K = 2.5·10^{-3}$ m/s and $Ly=1240$ m.</figcaption></figure>

years to reach the equilibrium of the water-table, despite approximately 95 % of this state was achieved in less than a year. For the theoretical example given, approximately one year is necessary to reach steady-state condition due to the higher diffusivity of the aquifer (0.5 m²/s) and the smaller length of the barrier (1240 m). However, considering the given variability of hydrogeological parameters, and particularly the diffusivity of the aquifer, this time can be doubled. This shows, as expected, that the diffusivity of the aquifer is the main controlling parameter. The optimal length of the barrier (Ly), its width (Lx), its hydraulic conductivity ($k'$) and its optimal distance to the pumping well ($x_{w_opt}$) were also assessed using the type-curves in Figs. 9 and 10. It is shown that, given the realistic ranges of hydrogeological parameters, the ranges evaluated for the technical characteristics of the barrier are relatively low (Ly: 1240±100 m, $x_{w_opt}$: 135±20 m, and $8·10^{-8} ≤ k'≤5·10^{-7}$ m/s with $3≤Lx≤54$ m). Nevertheless, the less precise the knowledge of the aquifer properties, the more imprecise the evaluation of the characteristics of the underground barrier will be.

The hydraulic conductivity $k'$ of the material, which should be used to construct the hypothetical underground barrier, was estimated between 0.8 and $5·10^{-7}$ m/s to expect a barrier efficiency of 95–99 %. This estimate suggests that it is not necessary to build the barrier with very low permeability materials (e.g. concrete, or impermeable diaphragm walls). In addition, thanks to field observations on an existing barrier suspected of leaking, this curve can also be used to evaluate the resistance parameter of the barrier, and the use of proposed solutions will make it possible to estimate the flux passing through it (leakage).

The given examples show that the proposed type-curves allows to estimate rapidly the characteristics of the underground physical barrier and their range of variability. This information can be used further to provide 2-D estimates of the hydraulic-head using the proposed solutions or numerical modelling, particularly if the question of migration/ accumulation of pollutants arises.

An analysis of the barrier resistance parameter was carried out on some barrier situations of Fig. 9 with $Lx·K/(Ly·k')=1.0$ and 10.0 and, aquifer properties, barrier geometries, and pumping rates in Table 3, to evaluate the change of the maximum increase of the water-table at the well location ($\Delta h_{wd,max}$) in case of semi-permeable barriers. The results (Fig. 13a, b) showed that the curve defined in Fig. 9a is unchanged (Fig. 13a), which means that the optimal location of the well ($x_{w_opt}$) is the same whatever the value of the resistance parameter. However, the curve defined in Fig. 9b for impermeable barrier (black dotted curve) is

slightly shifted downwards due to a decrease in $\Delta h_{wd,max}$ by approximately 10 % for barrier resistance parameter value of 10.0 – in this case the barrier is almost impermeable, and significantly shifted downwards due to a decrease in $\Delta h_{wd,max}$ by approximately 50 % for resistance parameter value of 1.0 (semi-permeable barrier). These results are very comparable to the type curve established for evaluating the efficiency of the barrier, for $Lx·K/(Ly·k')=1.0$, $\Delta h_{Lx,K/(Ly·k')}/(\Delta h_{Lx,K/(Ly·k')}=∞)≈0.5$ and, for $Lx·K/(Ly·k')=10.0$, $\Delta h_{Lx,K/(Ly·k')}/(\Delta h_{Lx,K/(Ly·k')}=∞)≈0.9$ (Fig. 10).

With the barrier-aquifer configuration studied, i.e. the barrier only concerns a limited part of the width of the aquifer, the main purpose of the underground barrier is to raise the level of the water-table, especially at the location of a pumping well, and possibly to increase the productivity of the well due to the increase in the saturated thickness. It is therefore a question of securing the exploitation of groundwater during dry periods, when water levels become too deep, to avoid a water shortage provided there is sufficient water coming from the upstream.

This configuration does not constitute a barrier to store groundwater upgradient as underground dams do in case of complete lateral obstruction of a valley (e.g. Sun et al., 2019; Liu et al., 2023; Masayuki, 2024). Furthermore, in the studied configuration, the additional volume of water due to the rise of the water-table may be lower than that extracted annually by pumping, depending on the storativity of the aquifer. Therefore, the use of underground barriers for improving the supply of groundwater is not limited to the case of complete obstruction of the groundwater flow. However, such barriers are not a common situation, which raises questions about their interest, the limits and constraints of their use, and must be examined from different aspects, in particular through field experiments. The limited or partial obstruction of the aquifer of a valley offers the advantage of always allowing groundwater to flow, even at the level of the water-table mound created by the underground barrier, thus limiting the accumulation of pollutants in the groundwater, a phenomenon sometimes encountered in cases of complete obstruction (e.g. Yoshimoto et al., 2010; Zheng et al., 2020; Gao et al., 2024).

Underground barriers cause a rise in water levels upgradient of the structure, which can benefit pumping wells or reduce seawater intrusion, but they can also flood basements or create new wetlands. Conversely, they also induce a decrease of water level downgradient which can decrease the yields of pumping wells -or even cause them to be abandoned, wetlands may fall dry, and seawater can intrude coastal areas. Hydrogeological and predesign studies are thus necessary to
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Dimensionless curves of the optimal location of a pumping well near a semi-permeable underground barrier. Constructed from simulations with Eq. (16), with few barrier and aquifer situations, and pumping rates from Table 3, with barrier hydraulic resistance parameter $L_x \cdot K \cdot (Ly \cdot k')=1.0$ and $10.0$. The black dotted curves are identical to Fig. 9 (impermeable barrier). a) dimensionless optimal distance of the well vs. dimensionless pumping rate. b) Inverse of the dimensionless efficiency of the barrier on the well vs. dimensionless optimal distance of the well.</figcaption></figure>

assess the environmental, economic and social impacts of all proposed underground barriers (e.g. Ishida et al., 2003; Yihdego, 2016; Attard et al., 2016; Li et al., 2018; Sun et al., 2019), and the semi-analytical solutions and type-curves developed here can be used in these studies.

Additionally, the proposed type-curves and semi-analytical solutions can be used in geotechnical studies, e.g. to evaluate the rise and decline of the water table at a given location caused by the construction of the basement of a building, as well as to design a pumping well to avoid flooding of basements located upgradient of the structure.

## 7. Conclusion

Transient and steady-state semi-analytical solutions for evaluating the hydraulic-head distribution caused by a semi-permeable rectangular underground physical barrier of insignificant storage capacity that fully penetrates an infinite, homogeneous and sloping aquifer, with or without a pumping well, are proposed based on the method of fundamental solutions (Eqs. (14), (16)). As the aquifer is assumed to be

infinite, the length of the barrier concerns a limited part of its width. The approach relies on the first method of linearization of the partial differential equation, which according to previous works, give very satisfactory results when the maximum degree of saturation/desaturation $\beta$ is less than 20 %. For higher values, results are approximated but can be used as a preliminary diagnosis.

The semi-analytical transient solution was satisfactorily compared with the full numerical modelling. The accuracy of the semi-analytical solutions was also evaluated, results are satisfactory for a ratio $d'/\sigma$ ranging between 2.0 and 3.0 ($d'$: distance of image wells to the side of the barrier, $\sigma$: distance between image wells), and corroborates previous studies with impermeable underground barriers with or without a pumping well (steady-state solutions, Dewandel et al., 2024a, 2024b).

Based on hundreds of simulations of hydraulic-head in transient and steady-state conditions, with or without a pumping well, different properties of the aquifer, different hydraulic conductivities and geometries of the barrier, a set of type curves was provided for evaluating: i) the duration of equilibrium of the water-table following the construction of the barrier, ii) the length of the underground barrier required to optimize the pumping rate of a given well or its location (or the inverse), and iii) the effectiveness of the barrier according to its width and its hydraulic conductivity. The time type-curve was applied to the construction of a railway station basement provided in the literature and gave satisfactory results. The use of all the type-curves was illustrated through a hypothetical example for the purpose of designing an underground barrier to raise the water-table level at the location of a pumping well. The variability of barrier characteristics introduced by uncertainty in hydrogeological parameters was also addressed. Apart from the time required to reach water-table equilibrium due to its sensitivity to aquifer diffusivity, the length, width and hydraulic conductivity of the barrier were estimated within a reasonable range for the hypothetical case studied. Estimating the hydraulic conductivity range of the underground barrier makes it possible to evaluate what type of material and/or technique should be used for the construction of the barrier. It is shown that it is not necessary to build a strictly impermeable structure to achieve a similar impact on the water-table, particularly on the rise of the water-table upgradient of the barrier.

The semi-analytical solutions developed consider an infinite homogeneous aquifer, future work should consider that the aquifer can be limited in space, either because of the limits of the valley where the aquifer developed (no-flow boundary) and/or because of rivers, a sea or a lake (e.g. constant head boundary), and the spatial variability of aquifer parameters. Results with these schemes should be different from those shown here.

## Funding information

This study was conducted and funded by BRGM and the French Ministry of Higher Education and Research under the ‘Underground barriers to manage groundwater resources: assessment of their potential benefits and limitations for adaptation to climate change’ project.

## CRediT authorship contribution statement

**Benoît Dewandel:** Methodology, Formal analysis, Conceptualization. **Sandra Lanini:** Methodology, Formal analysis. **Augustin Gouy:** Formal analysis, Data curation. **Nicolas Frissant:** Project administration, Methodology, 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.
Acknowledgements

useful remarks and comments, which improved the quality of the paper.

The two anonymous reviewers of the Journal are thanked for their

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.advwatres.2025.105091.

Appendix

A- Robin' condition along the sides of the underground physical barrier
This section explains how Eqs. (2), (4) and (5) were obtained.

<figure><img src="image_17.png" /><figcaption>Fig. A. Sketch of the boundary condition along a semi-permeable underground rectangular physical barrier. a) plan view and b) vertical section view.</figcaption></figure>

With the conventional rule of positive sign when entering the aquifer, the flow through the four sides of a rectangular barrier is given by the Darcy' law.

$$
q = \frac{k'}{L_x \text{ or } L_y} (h_{\text{Opp,side}} - h_{\text{Adj,side}}) \qquad (A-1)
$$

And, still according to Darcy’s law, this flow evaluated from the hydraulic head distribution in the aquifer is given by:

$$
q = \frac{\partial h}{\partial \vec{n}} \cdot K = K \begin{bmatrix} \frac{\partial h}{\partial x} \\ \frac{\partial h}{\partial y} \end{bmatrix} \cdot \begin{bmatrix} \cos \varphi \\ \sin \varphi \end{bmatrix} \qquad (A-2)
$$
Where K is the hydraulic conductivity of the aquifer, k' the hydraulic conductivity of the underground barrier (m/s), Lx and Ly its lengths, the storativity of the barrier is not significant (for transient condition). h<sub>Adj,side</sub> is the hydraulic head on the side of the barrier where the flux has to be computed (i.e. side A, B, C or D on Fig. Aa), Adj.: adjacent) and h<sub>Opp,side</sub>, the hydraulic head at the opposite side (i.e. side A, B, C or D on Fig. Aa). n<sup>-</sup> is the vector normal to the interface barrier/aquifer oriented outward the aquifer, φ is the angle of the vector according to the x-axis.

As the flow is identical at the interface, thus Eq. A-1 equals Eq. A-2, it is deduced the boundary conditions given by Eq. (2) in the text, which is, by definition, a mixed boundary condition, named Robin or Fourier condition. The following gives the boundary conditions on each side of the barrier.

Along the side A (Fig. Ab), upgradient face (φ= 0, n<sup>-</sup> = [1/0], -Ly/2 ≤ y≤Ly/2), the Robin condition is:

$$
K \frac{\partial h}{\partial x} \Big|_{x=\frac{Lx}{2}} = \frac{k'}{Lx} (h_{\text{side C}} - h_{\text{side A}}) = \frac{k'}{Lx} \left[ h\left(\frac{Lx}{2}, y, t\right) - h\left(\frac{Lx}{2}, y, t\right) \right] \quad (A-3)
$$

Along the side C, downgradient face (φ=π, n⁻ = [-1/0], -Ly/2 ≤ y≤Ly/2), the Robin condition is:

$$
-K \frac{\partial h}{\partial x} \Big|_{x=\frac{Lx}{2}} = \frac{k'}{Lx} (h_{\text{side A}} - h_{\text{side C}}) = \frac{k'}{Lx} \left[ h\left(-\frac{Lx}{2}, y, t\right) - h\left(\frac{Lx}{2}, y, t\right) \right] \quad (A-4)
$$

$$
\text{Therefore, } K \frac{\partial h}{\partial x} \Big|_{x=\frac{Lx}{2}} = \frac{k'}{Lx} \left[ h\left(\frac{Lx}{2}, y, t\right) - h\left(-\frac{Lx}{2}, y, t\right) \right] \quad (A-4)
$$

Then, Eqs. A-3 and A-4 follow the following form:

$$
K \frac{\partial h}{\partial x} \Big|_{x=\pm \frac{Lx}{2}} = \frac{k'}{Lx} \left[ h\left(\frac{Lx}{2}, y, t\right) - h\left(-\frac{Lx}{2}, y, t\right) \right] \quad (A-5)
$$

Along the side B, i.e. along the upper side (y=Ly/2) (φ = 3π/2, n⁻ = [0/-1], -Lx/2 ≤ x≤Lx/2), the Robin condition is:

$$
-K \frac{\partial h}{\partial y} \Big|_{y=\frac{Ly}{2}} = \frac{k'}{Ly} (h_{\text{side D}} - h_{\text{side B}}) = \frac{k'}{Ly} \left[ h\left(x, -\frac{Ly}{2}, t\right) - h\left(x, \frac{Ly}{2}, t\right) \right] \quad (A-6)
$$

$$
\text{Therefore } K \frac{\partial h}{\partial y} \Big|_{y=\frac{Ly}{2}} = \frac{k'}{Ly} \left[ h\left(x, \frac{Ly}{2}, t\right) - h\left(x, -\frac{Ly}{2}, t\right) \right] \quad (A-6)
$$

And along the side D, i.e. along the lower side (y=-Ly/2, φ = π/2, n⁻ = [0/1], -Lx/2 ≤ x≤Lx/2), the Robin condition is:

$$
K \frac{\partial s}{\partial y} \Big|_{y=\frac{Ly}{2}} = \frac{k'}{Ly} (h_{\text{side B}} - h_{\text{side D}}) = \frac{k'}{Ly} \left[ h\left(x, \frac{Ly}{2}, t\right) - h\left(x, -\frac{Ly}{2}, t\right) \right] \quad (A-7)
$$

Then, Eqs. A-6 and A-7 follow the following form:

$$
K \frac{\partial h}{\partial y} \Big|_{y=\pm \frac{Ly}{2}} = \frac{k'}{Ly} \left[ h\left(x, \frac{Ly}{2}, t\right) - h\left(x, -\frac{Ly}{2}, t\right) \right] \quad (A-8)
$$

Robin' condition with the first method of linearization

With s(x, y, t) = D(x, y, t) - D₀ and as h(x, y, t) = D(x, y, t) - ix

Eq. A-5 becomes (i.e. Eq. (4)) in text:

$$
K \frac{\partial s}{\partial x} \Big|_{x=\pm \frac{Lx}{2}} - \frac{k'}{Lx} \left[ s\left(\frac{Lx}{2}, y, t\right) - s\left(-\frac{Lx}{2}, y, t\right) \right] = i(K - k') \quad (A-9)
$$

And Eq. A-8 becomes (i.e. Eq. (5)) in text:

$$
K \frac{\partial s}{\partial y} \Big|_{y=\pm \frac{Ly}{2}} - \frac{k'}{Ly} \left[ s \left( x, \frac{Ly}{2}, t \right) - s \left( x, -\frac{Ly}{2}, t \right) \right] = 0 \quad (A-10)
$$

B- Transient semi-analytical solution with the pumping starting at the time of the building of the structure

Initial condition: h(x, y, 0) = D₀ - ix.

Semi-analytical solution:

$$
h(x, y, t) = \sum_{j=1}^{n} \tilde{r}_j^t G_j(x, y, x_j, y_j, t) + D_0 - i \cdot x + \frac{Q_p}{4\pi D_0 K} e^{i\alpha x_p} \cdot W(u, r_p \alpha) \quad (B-1)
$$

where W(u,rα) is the Hantush leaky function, e is the exponential function, Q<sub>p</sub> is the constant pumping rate (Q<sub>p</sub> < 0; in m<sup>3</sup>/s), α = i/(2D<sub>0</sub>), r =
√x<sup>2</sup> + y<sup>2</sup>, u = r<sup>2</sup>/4vt, and v = D<sub>0</sub>K/S is the diffusivity of the aquifer. See the text for the explanation of other parameters (Eq. (14)).
<figure><img src="image_18.png" /><figcaption>**Fig. B.** Comparison with numerical modelling for an impermeable barrier, similar to Fig. 2 but with an impermeable barrier. $K = 10^{-5}$ m/s, $S = 0.1$, $D_0=100$ m, $Ly=401$ m, $Lx=9$ m, $x_w=-21.0$ m, $y_w=0.0$ m, $Q_p=25$ m³/h. (a) Plan view with the location of the barrier, the pumping well and the observation wells (the location of wells are given by their colours (see inserted table)). b) Comparison of drawdown evolution vs. time (s) computed with MODFLOW 6 (triangles) and the MFS solution (Eq. (14); curves). Number of image wells (n): 700, distance of image wells to the side of the barrier (d'): 2.6 m. RMSE D: standardized root mean square error. c) hydraulic head cross-sections along the x-axis (at y = 0) at different times since the start of pumping.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. C.** Comparison with numerical modelling for a semi-permeable barrier, non-axisymmetric pumping well and other barrier and aquifer parameters. $K = 4 \cdot 10^{-6}$ m/s, $S = 0.05$, $i = 0.01$, $D_0=100$ m, $Ly=201$ m, $Lx=5$ m, $k'=10^{-7}$ m/s, $x_w=-21.0$ m, $y_w=20.0$ m $Q_p=10$ m$^3$/h. (a) Plan view with the location of the barrier, the pumping well and the observation wells (the location of wells are given by their colours (see inserted table). b) Comparison of drawdown evolution vs. time (s) computed with MODFLOW 6 (triangles) and the MFS solution (Eq. (14); curves). Number of image wells (n): 430, distance of image wells to the side of the barrier (d'): 2.1 m. RMSE D: standardized root mean square error. c) hydraulic head cross-sections along the x-axis (at y = 0) at different times since the start of pumping.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. D. Type-curves with empirical relationship (cumulative log-normal distribution function). a) Dimensionless time to reach steady-state condition (water-table equilibrium), with Lx·K/(Ly·k')=0.3, 1.0 and ∞ (impermeable barrier). b) Evaluation of the hydraulic conductivity of the underground barrier. R: linear regression coefficient.</figcaption></figure>

<caption>Fig. D</caption>

Data availability

The authors are unable or have chosen not to specify which data has
been used.

References

Abd-Elaty, I., KuriqiA, Pugliese L, Ahmed, A., 2024. Shoreline subsurface dams to protect coastal aquifers from sea level rise and saltwater intrusion. Appl. Water. Sci. 14, 49. https://doi.org/10.1007/s13201-023-02032-y, 2024.

Abdoulhalik, A., Ahmed, A.A., Hamill, G.A., 2017. A new physical barrier system for seawater intrusion control. J. Hydrol. 549, 416–427. https://doi.org/10.1016/j.jhydrol.2017.04.005.

Abdoulhalik, A., Ahmed, A.A., 2017a. How does layered heterogeneity affect the ability of subsurface dams to clean up coastal aquifers contaminated with seawater intrusion? J. Hydrol. 553, 708–721. https://doi.org/10.1016/j.jhydrol.2017.08.044.

Abdoulhalik, A., Ahmed, A.A., 2017b. The effectiveness of cutoff walls to control saltwater intrusion in multi-layered coastal aquifers: experimental and numerical

study. J. Environ. Manag. 199, 62–73. https://doi.org/10.1016/j.jenvman.2017.05.040.

Abdoulhalik, A., Abdelgawad, A.M., Ahmed, A.A., Moutari, S., Hamill, G., 2022. Assessing the protective effect of cutoff walls on groundwater pumping against saltwater upconing in coastal aquifers. J. Environ. Manag. 323, 116200. https://doi.org/10.1016/j.jenvman.2022.116200.

Abdoulhalik, A., Ahmed, A.A., Abd-Elaty, I., 2024. Effects of layered heterogeneity on mixed physical barrier performance to prevent seawater intrusion in coastal aquifers. J. Hydrol. 637, 131343. https://doi.org/10.1016/j.jhydrol.2024.131343.

Anwar, H.O., 1983. The effect of a subsurface barrier on the conservation of freshwater in coastal aquifers. Water Res. 17 (10), 1257–1265. https://doi.org/10.1016/0043-1354(83)90250-6.

Attard, G., Winiarski, Th, Rossier, Y., Eisenlohr, L., 2014. Underground constructions impacts on groundwater flow in urban context. Journées Nationales de Géotechnique et de Géologie de l'Ingénieur JNGG2014 – Beauvais, 8-10 juillet 2014.

Apaydin, A., 2024. Underground dams and Türkiye's experience. Nobel Bilimsel Eserler Sertifika No.: 20779, p. 144. E-ISBN: 978-625-393-982-3.

Attard, G., Winiarski, Th, Rossier, Y., Eisenlohr, L., 2016. Review: impact of underground structures on the flow of urban groundwater. Hydrogeol. J. 24, 5–19. https://doi.org/10.1007/s10040-015-1317-3.

BCEOM, 1978. Les Barrages Souterrains [The Underground Dams]. BCEOM -Paris, FR, Bureau Central d'Études pour les équipements d'Outre-Mer, Ministry of
Cooperation, series on Rural Techniques in Africa, Paris, p. 135. https://www.ircwa.sh.org/sites/default/files/217-78BA.pdf.

Bear, J., 1972. Dynamics of Fluid in Porous Media. Dover, New York.

Bouwer, H., 2002. Artificial recharge of groundwater: hydrogeology and engineering. Hydrogeol. J. 10, 121–142. https://doi.org/10.1007/s10040-001-0182-4.

Bhusari, V., Katpatal, Y.B., Kundal, P., 2016. An innovative artificial recharge system to enhance groundwater storage in basaltic terrain: example from Maharashtra, India. Hydrogeol. J. 24, 1273–1286. https://doi.org/10.1007/s10040-016-1387-x, 2016.

Carslaw, H.S., Jaeger, J.C., 1959. Conduction of heat in solids -2nd edition. 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.2010.11.022.

Chang, Q., Zheng, T., Zheng, X., Gao, C., Song, X., Walther, M., 2023. Repulsion driven by groundwater level difference around cut off walls on seawater intrusion in unconfined aquifers. Sci. Total Environ. 874, 162535. https://doi.org/10.1016/j.scitotenv.2023.162535, 2023.

Chaudhry, S.M., Zubair, M., 1994. Generalized incomplete gamma functions with applications. J. Comput. Appl. Math. 55, 99–124. https://doi.org/10.1016/0377-0427(94)90187-2, 1994.

Cheng, A.H.D., Hong, Y., 2020. An overview of the method of fundamental solutions -Solvability, uniqueness, convergence, and stability. Eng. Anal. Bound. Elem. 120, 118–152. https://doi.org/10.1016/j.enganabound.2020.08.013, 2020.

Choi, J.-R., Sung, J.H., Kang, D-H, Choo, K.-S., Kim, B.-S., 2023. A study on the priority selection method for underground dam installation considering humanities and social factors using fuzzy analytic hierarchy process in Korea. Water. (Basel) 15, 3296. https://doi.org/10.3390/w15183296, 2023.

De Marsily, G., 1986. Quantitative Hydrology for Engineers. Academic Press Inc.

Dewandel, B., Lanini, S., Frissant, N., 2024a. 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. Hydrogeol. J. 32, 801–816. https://doi.org/10.1007/s10040-023-02734-2.

Dewandel, B., Lanini, S., Frissant, N., 2024b. An improved steady-state semianalytical solution for assessing the two-dimensional hydraulic head distribution induced by an underground dam in a sloping unconfined aquifer. Hydrogeol. J. https://doi.org/10.1007/s10040-024-02831-w.

Dewandel, B., Lanini, S., Frissant, N., 2024c. Steady-state semi-analytical solutions for assessing the two-dimensional groundwater hydraulic head depletion induced by river dam removal. J. Hydrol. 642, 131801. https://doi.org/10.1016/j.jhydrol.2024.131801, 2024.

Dillon, P.J., 2005. Future management of aquifer recharge. Hydrogeol. J. 13, 313–316. https://doi.org/10.1007/s10040-004-0413-6.

El-Hames, A.S., 2012. Numerical solution for water table rise estimation behind deep underground dam. Groundwater 50, 621–626. https://doi.org/10.1111/j.1745-6584.2011.00861.x.

Fitts, C., 1997. Analytical modelling of impermeable and resistant barriers. Ground. Water. 35, 312–317. https://doi.org/10.1111/j.1745-6584.1997.tb00088.x.

Furukawa, H., Kuroda, T., 2011. The Underground dam, A New Concept in Water Resource Development to Meet the Challenges of Natural Disaster and Climate Change. Kaichosha Editors, p. 768.

Gao, S., Zheng, T., Zhang, B., Fang, Y., Zheng, X., 2024. Combined effects of aquifer heterogeneity and subsurface dam on nitrate contamination in coastal aquifers. J. Environ. Manage. 351, 119740. https://doi.org/10.1016/j.jenvman.2023.119740.

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

Hantush, M.S., 1964a. 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., 1964b. Hydraulics of Wells. In: Advances in Hydroscience, v.1. V.T. Chow. Academic Press ed.

Hanson, G., Nilsson, A., 1986. Ground-water dams for rural-water supplies in developing countries. Ground. Water. 24, 497–506. https://doi.org/10.1111/j.1745-6584.1986.tb01029.x, 1986.

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., 2005. Flow to vertical and nonvertical wells in leaky aquifers. J. Hydrol. Eng. 10, 477–484. https://doi.org/10.1061/(ASCE)1084-0699(2005).

Ishida, S., Kotoku, M., Abe, E., Fazal, M.A., Tsuchihara, T., Imaizumi, M., 2003. Construction of subsurface dams and these impacts on the environment. RMZ – Mater. Geoenviron. 50 (1), 149–152.

Izady, A., Khorshidi, M.S., Nikoo, M.R., Al-Maktoumi, A., Chen, M., Al-Mamari, H., Gandomi, M.H., 2021. Optimal water allocation from subsurface dams: a risk-based optimization approach. Water Res. Manag. 35, 4275–4290. https://doi.org/10.1007/s11269-021-02946-9 doi.org/10.1111/j.1745-6584.2012.00986.x.

Johansson, B.T., Lesnic, D., Reeve, T.H., 2011. A method of fundamental solutions for two-dimensional heat conduction. Int. J. Comput. Math. 1. https://doi.org/10.1080/00207160.2010.522233, 2011.

Kaleris, V., Ziosig, A.I., 2013. The effect of cutoff walls on saltwater intrusion and groundwater extraction in coastal aquifers. J. Hydrol. 476, 370–383. https://doi.org/10.1016/j.jhydrol.2012.11.007.

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. Reserv. 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.

Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, S., & Provost, A.M. (2021). MODFLOW 6: U.S. Geological Survey Modular Hydrologic Model-User guide. U.S. Geological Survey techniques and Methods, book 6, chap. A60. https://doi.org/10.3133/tm6A60.

Li, F., Chen, X., Liu, C., Lian, Y., He, L., 2018. Laboratory tests and numerical simulations on the impact of subsurface barriers to saltwater intrusion. Natural Hazards 91 (3), 1223–1235. https://doi.org/10.1007/s11069-018-3176-4.

Lin, Y.F., Kurylyk, B.L., Werner, A.D., Liu, C.Y., Solórzano-Rivas, C., Lin, J.H., 2025. Image-well solution for island aquifers with pumping, recharge, and complex coastlines. Adv. Water. Resour. 104977. https://doi.org/10.1016/j.adwatres.2025.104977. ISSN 0309-1708.

Liu, J., Brunner, Ph., Tokunaga, T., 2023. Impacts of subsurface dam construction on downstream groundwater levels and salinity in coastal aquifers. Groundwater 1–14. https://doi.org/10.1111/gwat.13304.

Masayuki, I., 2024. Review of subsurface dam technology based on Japan's experience in the Ryukyu Arc. Water 2024 16, 2282. https://doi.org/10.3390/w16162282.

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.

Mays, L., 2010. Ancient Water Technologies. Springer Dordrecht, p. 280. https://doi.org/10.1007/978-90-481-8632-7 ISBN: 978-90-481-8631-0.

Nishigaki, M., Kankam-Yeboah, K., Komatsu, M., 2004. Underground dam technology in some parts of the world. J. Groundwater Hydrol. 46, 113–130. https://doi.org/10.5917/jagh1987.46.113.

Polubarina-Kochna P.Y., 1977. Theory of ground water movement. Moscow: Nauka, in Russian. 665 p.

Pujades, E., López, A., Carrera, J., Vázquez-Suñé, E., Jurado, A., 2012. Barrier effect of underground structures on aquifers. Eng. Geol. 145–146, 41–49. https://doi.org/10.1016/j.enggeo.2012.07.004, 2012.

Shirahata, C., Yoshimoto, S., Tsuchihara, T., Nakazato, H., Ishida, S., 2022. A method for evaluating coastal underground barrier wall using groundwater tidal response. Groundwater 60 (6), 774–783. https://doi.org/10.1111/gwat.13201. Vol.

Sun, Y., Xu, S.-G., Kang, P.P., Fu, Y.Z., Wang, T.X., 2019. Impacts of artificial underground reservoir on groundwater environment in the reservoir and downstream area. Int. J. Environ. Res. Public Health 16, 1921. https://doi.org/10.3390/ijerph16111921, 2019.

Tsai, C.C., Young, D.L., Fan, C.M., Chen, C.W., 2006. MFS with time-dependent fundamental solutions for unsteady Stokes equations. Eng. Anal. Bound. Elem. 30, 897–908. https://doi.org/10.1016/j.enganabound.2006.04.006, 2006.

Tubau, I., 2004. Estudio Hidrogeológico y propuesta de Proceso Constructivo Para La Excavación De Un Túnel Entre Pantallas En El Delta Del Llobregat. Universitat Politecnica de Catalunya (UPC). Minor thesis.

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

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.

Yihdego, Y., 2016. Evaluation of flow reduction due to hydraulic barrier engineering structure: case of urban area flood, contamination and pollution risk assessment. Geotech. Geol. Eng. 1–14. https://doi.org/10.1007/s10706-016-0071-1.

Yin, J., Wang, N., Lu, C., Tsai, F.T.-C., Chen, H., 2023. Fast desalinization of residual saltwater using subsurface dams combined with saltwater discharge or freshwater recharge. J Hydrol 619, 129282. https://doi.org/10.1016/j.jhydrol.2023.129282, 2023.

Yoshimoto, S., Tsuchihara, T., Ishida, S., Imaizumi, M., 2010. Development of a numerical model for nitrates in groundwater in the reservoir area of the Komese subsurface dam, Okinawa, Japan. Environ. Earth Sci. 70, 2061–2077. https://doi.org/10.1007/s12665-011-1356-6, 2013.

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.

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

Williams, H., 2021. A note on the method of fundamental solutions. Academia Letters, 745. https://doi.org/10.20935/AL745.

Wu, H., Lu, C., Kong, J., Werner, A.D., 2020. Preventing seawater intrusion and enhancing safe extraction using finite-length, impermeable subsurface barriers: 3D analysis. Water. Resour. Res. 56. https://doi.org/10.1029/2020WR027792.

Wu, H., Lu, C., 2023. Seasonal fluctuations in the groundwater level accelerate the removal of residual saltwater upstream of subsurface dams. J. Hydrol. 625, 130026. https://doi.org/10.1016/j.jhydrol.2023.130026, 2023.

Wu, H., Lu, C., Shen, C., Ye, Y., 2023. Using a subsurface barrier to control seawater intrusion and enhance groundwater extraction in coastal aquifers: an analytical
study. J. Hydrol. 621, 129537. https://doi.org/10.1016/j.jhydrol.2023.129537, 2023.

Zhang, J., Lu, C., Shen, C., Zhang, C., Kong, J., Li, L., 2021. Effects of a low-permeability layer on unstable flow pattern and land-sourced solute transport in coastal aquifers. J. Hydrol. 598, 1–13. https://doi.org/10.1016/j.jhydrol.2021.126397.

Zheng, T., Zheng, X., Sun, O., Wang, L., Walther, M., 2020. Insights of variable permeability full-section wall for enhanced control of seawater intrusion and nitrate

contamination in unconfined aquifers. J. Hydrol. 586, 1–10. https://doi.org/10.1016/j.jhydrol.2020.124831.

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.