## 1. Introduction

It is well known since the pioneering work of Dunne [1969, 1980] and Dunne and Black [1970] that saturated groundwater flow (both subsurface and deep) plays an important role in forming and sustaining channelized flow, especially in humid environments. From this perspective, neither can streams be considered as merely fixed (Dirichlet) boundary conditions for the groundwater system, nor can the groundwater flow feeding the streams be considered a mere “baseflow contribution”: the energy level in the aquifer system and in the streams are strongly related for all flow conditions.

Wyns et al. [2004] showed that in low conductivity, shallow aquifer systems (such as weathered crystalline basement settings), there is a strong relation between the groundwater level at a given location, and the drop in topographic elevation from this location to the nearest wet talweg. More precisely, there is a correlation between the height above a basal envelope surface of talwegs (height labeled a on Figure 1), and the height of the groundwater level above this same envelope surface (height labeled b). In such systems, the low hydraulic conductivity implies a steepening of groundwater head gradients in order to convey groundwater from hillslopes to channels.

Obviously, the steepest gradients are achieved when the groundwater potentiometric surface is close to the topographic surface, i.e.

$$\frac{b}{a}\lesssim 1.$$ |

Interestingly, a parallel concept has been introduced in the field of flood mapping, namely the Height Above Nearest Drainage (HAND) of Nobre et al. [2011]. There are a couple of differences between the two relative elevations defined respectively by the Height Above the Basal Envelope Surface of Talwegs (HABEST) and the HAND. These differences are illustrated in Figure 1:

- the HAND computation procedure relies on a hydrologically-conditioned Digital Elevation Model (DEM) in which a reference drainage network is identified using a minimum support area. It then requires a downstream drainage allocation of each pixel, so that it creates a jump at the crossing of a water divide when two pixels on either side of the divide are “snapped” to two distinct drainages. Put differently, the reference potentiometric surface from which the HAND is computed is not continuous (we should call it the Elevation of Nearest Drainage, though it is not explicitly introduced by Nobre et al. [2011]).
- In contrast, the surface constructed following Wyns et al. [2004] is continuous. It is produced by interpolating point-scale elevation data in the talweg network.

From the point of view of free surface (channel) flow, a surface of minimum potential energy joining all the talwegs is of course, an important reference level: the transient rise of water during a flood event corresponds to a relative elevation above the “normal” drainage level, i.e. a local HAND value [Nobre et al. 2016].

Since the computation of such a base level has several applications, there is a need to give a formal definition and an unambiguous procedure to compute it.

## 2. A formal definition of the Basal Envelope Surface of Talwegs (BEST)

Let us consider the vertically-integrated, transient diffusivity equation for groundwater flow with the Dupuit–Forchheimer assumption [e.g., De Marsily 1986]:

$$S\frac{\partial h}{\partial t}+\overrightarrow{\mathit{\Nabla}}\cdot \overrightarrow{q}=r\phantom{\rule{1em}{0ex}}\text{with}\phantom{\rule{0.3em}{0ex}}\overrightarrow{q}=-T\left(h\right)\overrightarrow{\mathit{\Nabla}}h$$ | (1) |

^{2}⋅s

^{−1}), and r the (positively defined) recharge rate (m⋅s

^{−1}). Hence,

$$\overrightarrow{\mathit{\Nabla}}\cdot \left(T\left(h\right)\overrightarrow{\mathit{\Nabla}}h\right)=S\frac{\partial h}{\partial t}-r\left(x,y,t\right).$$ | (2) |

Assuming a uniform or slowly-varying transmissivity across the domain, we can take T out of the divergence operator. In the following, we will assume a strictly uniform transmissivity T_{0}, hence:

$${T}_{0}\mathit{\Delta}h=S\frac{\partial h}{\partial t}-r\left(x,y,t\right).$$ | (3) |

We can consider a hypothetical situation where the groundwater level reaches a minimum (∂h∕∂t ≈ 0) in conjunction with a vanishing recharge rate (r(x,y,t) ≈ 0), with the perennial streams still connected with the shallow aquifer. The head distribution would then satisfy:

$$\left\{\begin{array}{c}{\displaystyle}\mathit{\Delta}h=\frac{{\partial}^{2}h}{\partial {x}^{2}}+\frac{{\partial}^{2}h}{\partial {y}^{2}}=0\phantom{\rule{0ex}{2.0pt}}\phantom{\rule{1em}{0ex}}\hfill \\ h\left(x,y\right)={z}_{\text{topo}}\left(x,y\right)\phantom{\rule{1em}{0ex}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}\left(x,y\right)\phantom{\rule{0.3em}{0ex}}\text{isinthestreamnetwork}\hfill \end{array}\right.$$ | (4) |

where z_{topo} denotes the elevation of the topographic surface. Petroff et al. [2012] as well as Cohen et al. [2015] have studied in great detail the interaction of a stream network with a Poisson diffusive field in the more general case of a nonlinear transmissivity T(h) = K ⋅ (h − B) for an unconfined aquifer, where B(x,y) denotes the elevation of the underlying impermeable layer and K the saturated hydraulic conductivity of the aquifer. In the following, we will restrict ourselves to the linear case, as it allows simple analytic developments based on the superposition principle. Here, the condition 𝛥h = 0 means that the flux density $\overrightarrow{q}$ is assumed to have a vanishing divergence under hillslopes (it is a solenoidal vector field). The resulting piezometric surface is a minimum-curvature surface joining the 3-dimensional talweg lines, in a way similar to the gores of an umbrella’s canopy stretched over the ribs (Figure 2). Such a surface is an envelope surface in the sense of Häsing [1964], which can be defined in several ways.

Usually [see e.g., Deffontaines 1990; Wyns et al. 2004; Yao et al. 2017], this surface is constructed by kriging of scattered talweg elevation data: a more straightforward way of constructing it is to solve (4). In fact, under the assumption of uniform transmissivity we can directly model the analytic function using the Analytic Element Method [AEM, see e.g. Strack 1989; Steward 2020], instead of solving the steady-state equilibrium equation numerically with stream elevations as Dirichlet boundary conditions. As will be shown further, the procedure boils down to a simple linear least-square optimization, without the need for variography.

## 3. Mathematical modeling with analytic slit elements

### 3.1. Vector representation of the talweg graph

The approach proposed in this study relies on a vector representation of the talweg or stream network. Such representations are usually provided by national geographic agencies, as it is the case for the BD Topage in France [Système d’Information sur l’Eau/Sandre 2020], but can also be derived from the analysis of Digital Elevation Models (DEM). The network is a set of links, i.e., sections of a stream channel connecting two successive junctions, a junction and the outlet, or a junction and the drainage divide [ESRI 2022]. We assume that the topology of the links is also available in the form of an oriented graph. This general setting is illustrated in Figure 3. Each link is a polyline that can be decomposed into segments, each defined by a start point A_{j} and an end point B_{j}.

### 3.2. General form of the potential

In AEM, the head distribution is related to a real-valued potential 𝛷 [e.g., De Marsily 1986], such that:

$$\overrightarrow{q}=-\overrightarrow{\mathit{\Nabla}}\mathit{\Phi}.$$ |

For unconfined aquifers such as the shallow groundwater systems developed in weathered bedrock, our focus in this study, the general form of the flux is:

$$\overrightarrow{q}=-K\cdot \left(h-B\right)\overrightarrow{\mathit{\Nabla}}h.$$ |

In the case of a piecewise uniform elevation B of the substratum (i.e. $\overrightarrow{\mathit{\Nabla}}B=\overrightarrow{0}$ by subdomains), we have the integral form for the potential [e.g., Steward 2020]:

$$\begin{array}{ccc}{\displaystyle}& {\displaystyle}\overrightarrow{\mathit{\Nabla}}\mathit{\Phi}=-\overrightarrow{q}=+K\cdot \left(h-B\right)\overrightarrow{\mathit{\Nabla}}h& {\displaystyle}\\ {\displaystyle}& {\displaystyle}\overrightarrow{\mathit{\Nabla}}\mathit{\Phi}=K\cdot \left(h-B\right)\overrightarrow{\mathit{\Nabla}}\left(h-B\right)=\frac{K}{2}\overrightarrow{\mathit{\Nabla}}{\left(h-B\right)}^{2}& {\displaystyle}\\ {\displaystyle}& {\displaystyle}\Rightarrow \mathit{\Phi}=\frac{K}{2}{\left(h-B\right)}^{2}.& {\displaystyle}\end{array}$$ |

This expression hence leads to a nonlinear relation between potential 𝛷 (in m^{3}∕s) and head h (in m). In the following, we will linearize the problem and consider a uniform transmissivity T_{0} of the groundwater system, independent of h. This linearization is motivated by several considerations:

- first, as weathering develops from the surface down, it is no stronger an assumption to consider a rather uniform depth for the lower limit of the transmissive layer, than to consider a piecewise uniform elevation for this limit; this amounts to assume a rather uniform saturated thickness at first order, if the groundwater table is always close to the surface;
- second, our objective is to use observed groundwater levels in boreholes as well as topographic elevation of streams in order to infer large-scale hydrodynamic properties of the system such as, precisely, an “effective”, large-scale estimate of transmissivity. This is much easier to do when using groundwater head h directly as a proxy, linearly related to potential 𝛷: it will be clear when we present the degrees of freedom used to compute the value of potential 𝛷 using AEM in the next paragraphs.

With the assumption of a uniform transmissivity T_{0}, the potential 𝛷 is indeed linearly related to groundwater head h:

$$\overrightarrow{\mathit{\Nabla}}\mathit{\Phi}=-\overrightarrow{q}=+{T}_{0}\overrightarrow{\mathit{\Nabla}}h$$ |

$$\mathit{\Phi}={T}_{0}h+\text{constant}.$$ | (5) |

_{j}and B

_{j}, are first given a complex affix in the $\mathbb{C}$-plane (see Figure 3):

$$\left\{\begin{array}{cc}{a}_{j}={x}_{{A}_{j}}+\text{i}{y}_{{A}_{j}}\phantom{\rule{1em}{0ex}}\hfill & \text{affixofstartpoint}\phantom{\rule{0.3em}{0ex}}{A}_{j}\phantom{\rule{0.3em}{0ex}}\text{of}\hfill \\ \phantom{\rule{1em}{0ex}}\hfill & \text{element}\phantom{\rule{0.3em}{0ex}}j\phantom{\rule{0.3em}{0ex}}\text{inthecomplexplane}\phantom{\rule{0ex}{2.0pt}}\hfill \\ {b}_{j}={x}_{{B}_{j}}+\text{i}{y}_{{B}_{j}}\phantom{\rule{1em}{0ex}}\hfill & \text{affixofendpoint}\phantom{\rule{0.3em}{0ex}}{B}_{j}\phantom{\rule{0.3em}{0ex}}\text{of}\hfill \\ \phantom{\rule{1em}{0ex}}\hfill & \text{element}\phantom{\rule{0.3em}{0ex}}j\phantom{\rule{0.3em}{0ex}}\text{inthecomplexplane}\hfill \end{array}\right.$$ | (6) |

with i, a number such that i^{2} = −1. Hereafter, the letter z will denote an affix z = x + i y in the complex plane (with the exception of z_{topo} which will denote the topographic elevation). In the case of a divergence-free flux density $\overrightarrow{q}$, the real-valued potential 𝛷 is the real part of a complex-valued potential 𝛺:

$$\mathit{\Omega}\left(z\right)=\mathit{\Phi}\left(x,y\right)+\text{i}\phantom{\rule{2.77695pt}{0ex}}\mathit{\Psi}\left(x,y\right)$$ |

$$\begin{array}{ccc}{\displaystyle}{\left.\frac{\text{d}\mathit{\Omega}}{\text{d}z}\right|}_{{z}_{0}\in \mathbb{C}}& =& {\displaystyle}{lim}_{z\to {z}_{0}}\frac{\mathit{\Omega}\left(z\right)-\mathit{\Omega}\left({z}_{0}\right)}{z-{z}_{0}}\phantom{\rule{1em}{0ex}}\text{exists:doesnotvary}\\ {\displaystyle}& & {\displaystyle}\text{dependingonthedirectionfromwhich}\\ {\displaystyle}& & {\displaystyle}{z}_{0}\phantom{\rule{0.3em}{0ex}}\text{isapproached}.\end{array}$$ |

$$\left\{\begin{array}{cc}{\displaystyle}\frac{\partial \mathit{\Phi}}{\partial x}=\frac{\partial \mathit{\Psi}}{\partial y}\phantom{\rule{1em}{0ex}}\hfill & \text{(i)}\phantom{\rule{0ex}{6.0pt}}\hfill \\ {\displaystyle}\frac{\partial \mathit{\Phi}}{\partial y}=-\frac{\partial \mathit{\Psi}}{\partial x}\phantom{\rule{1em}{0ex}}\hfill & \text{(ii)}\phantom{\rule{0ex}{6.0pt}}\hfill \\ {\displaystyle}\frac{{\partial}^{2}\mathit{\Phi}}{\partial {x}^{2}}+\frac{{\partial}^{2}\mathit{\Phi}}{\partial {y}^{2}}=\mathit{\Delta}\mathit{\Phi}=0\phantom{\rule{1em}{0ex}}\hfill & \text{(iii),consequenceof(i)and(ii);}\phantom{\rule{0ex}{-6.0pt}}\hfill \\ \phantom{\rule{1em}{0ex}}\hfill & \text{theactualpropertywewant}\hfill \\ \phantom{\rule{1em}{0ex}}\hfill & \text{for}\phantom{\rule{0.3em}{0ex}}\mathit{\Phi}\hfill \\ {\displaystyle}\frac{{\partial}^{2}\mathit{\Psi}}{\partial {x}^{2}}+\frac{{\partial}^{2}\mathit{\Psi}}{\partial {y}^{2}}=\mathit{\Delta}\mathit{\Psi}=0\phantom{\rule{1em}{0ex}}\hfill & \text{(iv),consequenceof(i)and(ii)}.\hfill \end{array}\right.$$ |

Following the approach of Steward [2015, 2020], we will then consider each segment j of the talweg network as a slit element creating a complex potential 𝛺_{j}(z) at location z, with the form:

$${\mathit{\Omega}}_{j}\left(z\right)={Q}_{j}\times \left(\frac{ln{Z}_{j}\left(z\right)}{2\pi}-1\right)+\sum _{k=1}^{{n}_{\text{coef}}}{c}_{j,k}{Z}_{j}^{-k}\left(z\right).$$ | (7) |

Q_{j} is the net algebraic discharge (in m^{3}∕s) entering the slit, and the c_{j,k} are the complex coefficients of a Laurent series that can be adjusted to match various boundary conditions. If Q_{j} > 0, the slit is a sink to the system (it subtracts water from the aquifer) while if Q_{j} < 0, it is a source (it leaks water to the aquifer). The mapping function Z_{j}(z) for each slit element is given by:

$$\begin{array}{ccc}{\displaystyle}{\mathit{\zeta}}_{j}\left(z;{a}_{j},{b}_{j}\right)& =& {\displaystyle}\frac{z-\frac{1}{2}\left({a}_{j}+{b}_{j}\right)}{\frac{1}{2}\left({b}_{j}-{a}_{j}\right)}\end{array}$$ | (8) |

$$\begin{array}{ccc}{\displaystyle}{Z}_{j}\left(z\right)& =& {\displaystyle}{\mathit{\zeta}}_{j}\left(z\right)+\sqrt{{\mathit{\zeta}}_{j}\left(z\right)+1}\sqrt{{\mathit{\zeta}}_{j}\left(z\right)-1}.\end{array}$$ | (9) |

The function $Z\left(\mathit{\zeta}\right)=\mathit{\zeta}+\sqrt{\mathit{\zeta}+1}\sqrt{\mathit{\zeta}-1}$ is the Joukowsky transformation [Joukowsky 1910]. Figure 4 illustrates the two steps of the transformation ${Z}_{j}\left(z\right)={\mathit{\zeta}}_{j}\left(z\right)+\sqrt{{\mathit{\zeta}}_{j}\left(z\right)+1}\sqrt{{\mathit{\zeta}}_{j}\left(z\right)-1}$. The first step (mapping from z-plane to 𝜁-plane, Equation (8)) is just a similarity, i.e. a combination of translation, rotation and scaling that sends the two endpoints A_{j} and B_{j} of slit j to affixes −1 and +1 respectively. Then, the Joukowsky transformation (from 𝜁-plane to Z-plane, Equation (9)) maps the $\mathbb{R}\times \mathbb{R}$ plane to the outside of the unit circle: it amounts to “pop-opening” the slit segment [A_{j} B_{j}] conformally, i.e. with a transformation that conserves angles.

By superposition, the total (complex) potential at location z in the complex plane is simply the sum of the contributions of all slit elements in the domain:

$${\mathit{\Omega}}_{\text{tot}}\left(z\right)={\mathit{\Phi}}_{\text{tot}}\left(z\right)+\text{i}\phantom{\rule{2.77695pt}{0ex}}{\mathit{\Psi}}_{\text{tot}}\left(z\right)={\mathit{\Phi}}_{0}-\overline{{q}_{0}}\phantom{\rule{0.3em}{0ex}}z+\sum _{j=1}^{{n}_{\text{slit}}}{\mathit{\Omega}}_{j}\left(z\right)$$ | (10) |

_{0}= q

_{0,x}+ i q

_{0,y}is a complex number giving the background, regional direction of the flow (“tilt” in the groundwater table; since it is uniform, it does not add any divergence to the flux). $\overline{{q}_{0}}={q}_{0,x}-\text{i}\phantom{\rule{2.77695pt}{0ex}}{q}_{0,y}$ is its complex conjugate, and 𝛷

_{0}an offset value for the real potential 𝛷 = Re{𝛺}.

For each element j, we have (2n_{c} + 1) real degrees of freedom: the sink term Q_{j}, the n_{c} real parts Re{c_{j,k}} and n_{c} imaginary parts Im{c_{j,k}} with c_{j,k} = Re{c_{j,k}} + i Im{c_{j,k}}. Together with 𝛷_{0}, v_{0,x}, and v_{0,y}, the total number of real degrees of freedom is thus N = 3 + n_{s} × (2n_{c} + 1).

In Figure 5 we show the general form of the potential for three cases:

- a pure sink element with Q≠0 and all terms set to zero in the Laurent series; far from the slit, the real potential created by such an element is identical to the well known Thiem solution h = 𝛷∕T
_{0}= (Q∕2πT_{0})lnr + constant [e.g., De Marsily 1986], - the potential 𝛺(z) = c
_{1}Z^{−1}(z) with a real-valued c_{1}coefficient. In this case, the real potential 𝛷 (gray levels) is continuous, - the potential 𝛺(z) = c
_{1}Z^{−1}(z) with a pure imaginary c_{1}coefficient. In this case, the real potential 𝛷 is discontinuous across the slit.

The two last situations do no produce any net inflow nor outflow.

Of course, in the general case we have no idea of the transmissivity of the porous medium; in order to work directly with head levels/drops we will normalize the expression of the complex potential by the (unknown) uniform transmissivity T_{0}:

$$\begin{array}{ccc}{\displaystyle}{\mathit{\Omega}}_{\text{tot}}^{\prime}\left(z\right)& =& {\displaystyle}\frac{1}{{T}_{0}}{\mathit{\Phi}}_{\text{tot}}\left(z\right)+\text{i}\phantom{\rule{2.77695pt}{0ex}}\frac{1}{{T}_{0}}{\mathit{\Psi}}_{\text{tot}}\left(z\right)={\mathit{\Phi}}_{0}^{\prime}-\overline{{q}_{0}^{\prime}}\phantom{\rule{2.77695pt}{0ex}}z\\ {\displaystyle}& & {\displaystyle}+\phantom{\rule{2.77695pt}{0ex}}\sum _{j=1}^{{n}_{\text{slit}}}\left[{Q}_{j}^{\prime}\times \left(\frac{ln{Z}_{j}\left(z\right)}{2\pi}-1\right)+\sum _{k=1}^{{n}_{\text{coef}}}{c}_{j,k}^{\prime}{Z}_{j}^{-k}\left(z\right)\right]\\ {\displaystyle}& & {\displaystyle}\end{array}$$ | (11) |

$$\left\{\begin{array}{cc}{Q}_{j}^{\prime}={\displaystyle}\frac{{Q}_{j}}{{T}_{0}}\phantom{\rule{1em}{0ex}}\hfill & \mathit{normalized}\phantom{\rule{0.3em}{0ex}}\text{dischargeextractedby}\phantom{\rule{0ex}{-5.39996pt}}\hfill \\ \phantom{\rule{1em}{0ex}}\hfill & \text{slit}\phantom{\rule{0.3em}{0ex}}j\phantom{\rule{0.3em}{0ex}}\text{(depthordrawdowninm)}\hfill \\ {\mathit{\Phi}}_{0}^{\prime}={\displaystyle}\frac{{\mathit{\Phi}}_{0}}{{T}_{0}}\phantom{\rule{1em}{0ex}}\hfill & \phantom{\rule{0ex}{5.39996pt}}\hfill \\ {q}_{0}^{\prime}={\displaystyle}\frac{1}{{T}_{0}}{q}_{0}\phantom{\rule{1em}{0ex}}\hfill & \phantom{\rule{0ex}{5.39996pt}}\hfill \\ {c}_{j,k}^{\prime}={\displaystyle}\frac{1}{{T}_{0}}{c}_{j,k}.\phantom{\rule{1em}{0ex}}\hfill & \hfill \end{array}\right.$$ |

### 3.3. Optimization procedure

Given a set of control points {z_{p}}_{1⩽p⩽M} where we want to specify a value ${h}_{p}={h}_{\text{tot}}\left({z}_{p}\right)=\text{Re}\left\{{\mathit{\Omega}}_{\text{tot}}^{\prime}\left({z}_{p}\right)\right\}$ of the groundwater head, we can build a linear system to estimate the optimum value of all degrees of freedom. With the assumption of a uniform transmissivity, h is linear with respect to all degrees of freedom (11) so the regression matrix is simply the [M × N] Jacobian matrix:

$${J}_{p,l}=\text{Re}\left\{\frac{\partial {\mathit{\Omega}}_{\text{tot}}^{\prime}\left({z}_{p}\right)}{\partial {\mathit{\beta}}_{l}}\right\}$$ |

_{l}is the l-th real degree of freedom. The entry of the matrix on row p and column l reads as:

$$\left\{\begin{array}{cc}{J}_{p,l}=1\phantom{\rule{1em}{0ex}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}{\mathit{\beta}}_{l}={\mathit{\Phi}}_{0}^{\prime}\hfill \\ {J}_{p,l}=-\text{Re}\left\{{z}_{p}\right\}\phantom{\rule{1em}{0ex}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}{\mathit{\beta}}_{l}={q}_{0,x}^{\prime}\hfill \\ {J}_{p,l}=+\text{Im}\left\{{z}_{p}\right\}\phantom{\rule{1em}{0ex}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}{\mathit{\beta}}_{l}={q}_{0,y}^{\prime}\phantom{\rule{0ex}{1.79993pt}}\hfill \\ {J}_{p,l}=\text{Re}\left\{\frac{ln{Z}_{j}\left({z}_{p}\right)}{2\pi}\right\}-1\phantom{\rule{1em}{0ex}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}{\mathit{\beta}}_{l}={Q}_{j}^{\prime}\phantom{\rule{0ex}{1.79993pt}}\hfill \\ {J}_{p,l}=+\text{Re}\left\{{Z}_{j}^{-k}\left({z}_{p}\right)\right\}\phantom{\rule{1em}{0ex}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}{\mathit{\beta}}_{l}=\text{Re}\left\{{c}_{j,k}^{\prime}\right\}\phantom{\rule{0ex}{1.79993pt}}\hfill \\ {J}_{p,l}=-\text{Im}\left\{{Z}_{j}^{-k}\left({z}_{p}\right)\right\}\phantom{\rule{1em}{0ex}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}{\mathit{\beta}}_{l}=\text{Im}\left\{{c}_{j,k}^{\prime}\right\}.\hfill \end{array}\right.$$ | (12) |

The control points are positioned along the segments of the stream network. Since there are 2n_{c} + 1 degrees of freedom for each slit, the system is overdetermined if we position more control points than this value along each segment. The vector of optimum values for the parameters (the 𝛽_{l}’s) is the least squares solution to:

$$\mathbf{J}\phantom{\rule{0.3em}{0ex}}{\left[{\mathit{\beta}}_{1},{\mathit{\beta}}_{2},\dots ,{\mathit{\beta}}_{N}\right]}^{T}={\left[{\mathit{\Phi}}_{1}^{\prime},{\mathit{\Phi}}_{2}^{\prime},\dots ,{\mathit{\Phi}}_{M}^{\prime}\right]}^{T}.$$ |

Instead of Dirichlet conditions at control points, we can also specify Neumann conditions on the normalized velocity vector

$$q={q}_{x}+\text{i}{q}_{y}=-\left[\frac{\partial \mathit{\Phi}}{\partial x}+\text{i}\frac{\partial \mathit{\Phi}}{\partial y}\right]=-\overline{\frac{\text{d}\mathit{\Omega}}{\text{d}z}}.$$ |

Using the chain rule

$$\frac{\text{d}{\mathit{\Omega}}_{j}}{\text{d}z}=\frac{\text{d}{\mathit{\Omega}}_{j}}{\text{d}{Z}_{j}}\frac{\text{d}{Z}_{j}}{\text{d}{\mathit{\zeta}}_{j}}\frac{\text{d}{\mathit{\zeta}}_{j}}{\text{d}z},$$ |

$$\begin{array}{ccc}{\displaystyle}{q}_{j}\left(z\right)& =& {\displaystyle}-\overline{\frac{\text{d}{\mathit{\Omega}}_{j}}{\text{d}z}}\\ {\displaystyle}& =& {\displaystyle}\overline{\frac{2}{\left({b}_{j}-{a}_{j}\right)\sqrt{{\mathit{\zeta}}_{j}\left(z\right)+1}\sqrt{{\mathit{\zeta}}_{j}\left(z\right)-1}}\left[-\frac{{Q}_{j}}{2\pi}+\sum _{k=1}^{{n}_{\text{coef}}}{c}_{j,k}k\phantom{\rule{2.77695pt}{0ex}}{Z}_{j}^{-k}\left(z\right)\right]}\\ {\displaystyle}{q}_{j}\left(z\right)& =& {\displaystyle}\overline{\frac{2}{\left({b}_{j}-{a}_{j}\right)\sqrt{{\mathit{\zeta}}_{j}\left(z\right)+1}\sqrt{{\mathit{\zeta}}_{j}\left(z\right)-1}}}\\ {\displaystyle}& & {\displaystyle}\times \left[-\frac{{Q}_{j}}{2\pi}+\sum _{k=1}^{{n}_{\text{coef}}}\left(\text{Re}\left\{{c}_{j,k}\right\}-\text{i}\phantom{\rule{2.77695pt}{0ex}}\text{Im}\left\{{c}_{j,k}\right\}\right)k\phantom{\rule{2.77695pt}{0ex}}\overline{{Z}_{j}^{-k}\left(z\right)}\right]\end{array}$$ |

$$\begin{array}{ccc}{\displaystyle}{q}_{j}^{\prime}\left(z\right)& =& {\displaystyle}\overline{\frac{2}{\left({b}_{j}-{a}_{j}\right)\sqrt{{\mathit{\zeta}}_{j}\left(z\right)+1}\sqrt{{\mathit{\zeta}}_{j}\left(z\right)-1}}}\\ {\displaystyle}& & {\displaystyle}\times \left[-\frac{{Q}_{j}^{\prime}}{2\pi}+\sum _{k=1}^{{n}_{\text{coef}}}\left(\text{Re}\left\{{c}_{j,k}^{\prime}\right\}-\text{i}\phantom{\rule{2.77695pt}{0ex}}\text{Im}\left\{{c}_{j,k}^{\prime}\right\}\right)k\phantom{\rule{2.77695pt}{0ex}}\overline{{Z}_{j}^{-k}\left(z\right)}\right].\end{array}$$ |

This expression is again linear with respect to all real degrees of freedom.

## 4. Two applications of the Basal Envelope Surface of Talwegs

### 4.1. Study site: the Evel catchment and Coët-Dan subcatchment

We illustrate the concepts of the previous section on a small catchment located in Brittany (Western France), the Coët-Dan catchment, where surveys have been conducted since the early 1970s [Fovet et al. 2018]. This catchment is developed on brioverian shale with low hydraulic conductivity. In order to avoid edge effect, we consider the larger catchment of the Upper Evel River in which the Coët-Dan catchment is completely nested (gauging site J5618322—Evel at Remungol).

#### 4.1.1. Digital elevation model and derived vector datasets

Figure 6 shows the map of topographic elevation. The original data is the 1-m resolution RGE ALTI [Institut Géographique National 2018]; in order to speed up the extraction of the hydrographic network without altering channel elevations, we downscale this dataset to a resolution of 5 m by keeping the minimum elevation on each 5 × 5 pixels block. The resulting DEM is then hydrologically conditioned, and the flow direction and flow accumulation maps are computed. The stream network is extracted in vector form from the DEM with a support area of 0.1 km^{2} [ESRI 2022]. The boundary of the Evel catchment is also extracted in vector form, as the segments making up this boundary will be assigned a zero normal velocity condition in the least squares regression. In order to reduce the geometrical complexity, the stream links and the boundary polyline are simplified using the Douglas–Peucker algorithm [Douglas and Peucker 1973] with a planar tolerance of 40 m.

#### 4.1.2. Estimation of the Basal Envelope Surface of Talwegs

The total number of segments in the geometric setting is given in Table 1, along with the associated degrees of freedom.

**Table 1.**

List of the degrees of freedom in the least squares problem

Type of elements | Number of elements | Number of sink terms ${Q}_{j}^{\prime}$ | Number of $\text{Re}\left\{{c}_{j,k}^{\prime}\right\}$ | Number of $\text{Im}\left\{{c}_{j,k}^{\prime}\right\}$ | Number of degrees of freedom |
---|---|---|---|---|---|

Stream segm. | 1562 | 0 | 2 × 1562 | 0 | 3124 |

Boundary segm. | 176 | 0 | 2 × 176 | 2 × 176 | 704 |

Total incl. ${\mathit{\Phi}}_{0}^{\prime}$ | 3829 |

We use n_{coef} = 2 complex coefficients for the Laurent series on each slit element, that is, up to 4 real degrees of freedom per slit. However, we recall that in order for the real potential 𝛷 to be continuous between the two “sides” of a slit, the imaginary parts of these coefficients must be zero: hence there are only 2 real degrees of freedom for stream slit elements. This condition of continuity is not required for boundary slit elements: here we only want the normal velocity component to be zero on the “inner” side of each slit. We completely disregard the value of the real potential outside of the boundary, so we use both a real part Re{c_{j,k}} and an imaginary Im{c_{j,k}} part for the 2 complex coefficients of each boundary slit element. We set the background velocity to zero (v_{0,x} = v_{0,y} = 0), but the constant 𝛷_{0} is needed in any case and yields an additional degree of freedom, rising this number to 3829.

In order to get a least squares solution for all degrees of freedom, we use the overspecification principle of Janković and Barnes [1999] with M = 8 control points on each slit, whose positions are given by:

$${\mathit{\zeta}}_{m}=cos\left(\pi \frac{m-\frac{1}{2}}{M}\right)\phantom{\rule{1em}{0ex}}\left(m=1,2,\dots ,M\right).$$ |

The total number of control points is then M × n_{slit} = 8 × (1562 + 176) = 13,904, much larger than the number of degrees of freedom (3829). This resulting system is 13,904 × 3829, which is rather large, but still very easily solved by classical linear algebra packages. If we were to solve an even larger system, we could use the iterative solution introduced by Janković and Barnes [1999], sequentially solving for the coefficients of each slit to match its boundary conditions while holding all other coefficients constant for other slits, until convergence is reached [see also, e.g., Steward 2020].

Figure 7 shows the resulting envelope surface, that is, the real part ${\mathit{\Phi}}_{\text{BEST}}^{\prime}$ of the normalized complex potential 𝛺′. Elevation ranges from about 50 m a.s.l. at the outlet to about 140 m in the North of the catchment.

### 4.2. HABEST as a continuous substitute of HAND for flood mapping

Having built a continuous energy level joining all talweg lines, we can now compute the difference between the Basal Envelope Surface of Talwegs (BEST) and the topographic surface. This difference can be coined Height Above the Basal Envelope Surface of Talwegs (HABEST). It is conceptually the same thing as the Height Above Nearest Drainage’ (HAND), with two advantages:

- It is continuous in space. The HAND computation procedure requires a downstream drainage allocation of each pixel, so it creates a jump at the crossing of a water divide when two pixels on either side of the divide are “snapped” to two distinct drainages;
- if a reference hydrographic network is available in vector form, there is no need for hydrological preconditioning of the DEM (cf. Section 4.1.1).

Figure 8 shows the resulting index ${\mathit{\Phi}}_{\text{HABEST}}^{\prime}=\left({\mathit{\Phi}}_{\text{topo}}^{\prime}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}\right)=\left({z}_{\text{topo}}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}\right)$. Contours of this surface can be used for flood mapping in the same way as HAND contours.

### 4.3. Groundwater level interpolation

#### 4.3.1. Regression using HABEST as predictor

In this section, we use the resulting HABEST index as a predictor for groundwater level, following the idea of Wyns et al. [2004]. Figure 9 is a zoom on the map of Figure 7 on the Coët-Dan subcatchment; it also shows the location of 74 boreholes in which the groundwater level was surveyed during the month of February, 1996 [Pauwels et al. 1996]. The model reads:

$${\mathit{\Phi}}_{\text{GW}}^{\prime}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}=a\left({\mathit{\Phi}}_{\text{topo}}^{\prime}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}\right)+b$$ | (13) |

_{0}), and ${\mathit{\Phi}}_{\text{topo}}^{\prime}={z}_{\text{topo}}$ the topographic elevation.

The result of the regression is shown in Figure 10. The parameters are:

$$\left\{\begin{array}{cc}a=0.912\phantom{\rule{1em}{0ex}}\hfill & \phantom{\rule{0ex}{-1.79993pt}}\hfill \\ b=-1.74\phantom{\rule{0.3em}{0ex}}\text{m}\phantom{\rule{1em}{0ex}}\hfill & \hfill \end{array}\right.$$ |

^{2}= 0.962 and a rooted mean squared error (RMSE) of 1.67 m.

It is worth noting that the estimate of groundwater level using (13) is not completely insensitive to the value of the support area chosen for extracting the channel network. The RMSE obtained for different values of the support area is shown in Figure 11; we see that the optimal value is about A_{c} = 0.4 km^{2}.

#### 4.3.2. A non-zero divergence solution with area-sink

The definition of the envelope surface as a harmonic potential (corresponding to a divergence-free velocity field) is straightforward, and the result shown in Figure 7 exactly corresponds to the specifications so far. However, if we look closer at Figure 9, we see that slits corresponding to “leaves” in the network graph often correspond to a maximum of the envelope surface, just as if the surface was locally “hanging” from those elements, with the flux oriented from the slit to the domain (see points denoted with a “L”). It is no big issue if we do not want to give any particular physical meaning to this envelope surface; however, this is a spurious behaviour if we want to interpret it, for instance, as a (dry-season) groundwater level.

We can achieve a slightly more satisfactory interpolating surface by computing a non-zero divergence velocity field, still using AEM. For this we will set a non-zero sink term Q_{j} for each stream slit, in conjunction with an area sink [Strack 2017] with uniform infiltration rate. The Laplacian of the potential now satisfies the following steady-state equation away from the slits:

$${T}_{0}\mathit{\Delta}h={T}_{0}\mathit{\Delta}{\mathit{\Phi}}^{\prime}=-{r}_{0}$$ |

_{0}> 0 is a recharge depth per unit time (positively-defined, source term). As the transmissivity is unknown, it is again more convenient to rewrite this equation as:

$$\mathit{\Delta}{\mathit{\Phi}}^{\prime}=-\frac{{r}_{0}}{{T}_{0}}=+{\mathit{\gamma}}_{0}$$ |

_{0}is set according to the convention used in Strack [2017]. Just as we normalized ${Q}_{j}^{\prime}={Q}_{j}\u2215{T}_{0}$ for a slit sink term, 𝛾

_{0}= −r

_{0}∕T

_{0}is a normalized areal extraction rate term which has the dimension of the inverse of a depth (e.g., m

^{−1}).

Since we add a source term in the system, we need to add one or several sink term(s) somewhere else: obviously, it has to be at stream slits. The idea is to make the sink term Q_{j} for each slit element j a function of 𝛾_{0}. Let consider the portion of the stream slit graph shown in Figure 12: if we assume that we have only very local flow systems, we can consider that each slit cannot extract more than the amount of recharge on the subcatchment, it drains. Denoting dA_{j} the area of this subcatchment, which is a purely geometric quantity that can be easily computed in a Geographic Information System (GIS), we will even consider that the two quantities exactly match:

$${Q}_{j}={r}_{0}\phantom{\rule{2.77695pt}{0ex}}\text{d}{A}_{j}$$ |

$${Q}_{j}^{\prime}=\frac{{r}_{0}}{{T}_{0}}\phantom{\rule{2.77695pt}{0ex}}\text{d}{A}_{j}=-{\mathit{\gamma}}_{0}\phantom{\rule{2.77695pt}{0ex}}\text{d}{A}_{j}.$$ |

It is worth reiterating the fact that even though Q_{j} represents a flux extracted from the groundwater system, it does not necessarily represent a contribution to streamflow. A fraction of this flux might as well represent riparian or wetland evapotranspiration concentrated in the vicinity of the talweg, especially for the smallest streams that often dry up during summertime.

For the sake of brevity, the reader is referred to Strack [2017] for the expression of 𝛷_{areal}(z) and the associated velocity v_{areal}(z). Both are linear with respect to the new degree of freedom 𝛾_{0}. Of course, there is no stream function 𝛹_{areal} associated with 𝛷_{areal}: because of the non-zero divergence there is no holomorphic complex potential. Figure 13 shows the contours of 𝛷_{areal} for a unit extraction rate 𝛾_{0} = 1.

By superposition, the real potential produced by all elements in the domain finally reads:

$$\begin{array}{ccc}{\displaystyle}& & {\displaystyle}{\mathit{\Phi}}_{\text{GW}}^{\prime}\left(z\right)={\mathit{\Phi}}_{0}^{\prime}+\text{Re}\left\{\overline{{v}_{0}^{\prime}}z\right\}+{\mathit{\gamma}}_{0}{\mathit{\Phi}}_{\text{areal}}^{\prime}\left(z;{\mathit{\gamma}}_{0}=1\right)\\ {\displaystyle}& & {\displaystyle}\phantom{\rule{1em}{0ex}}+\phantom{\rule{2.77695pt}{0ex}}\sum _{j=1}^{{n}_{\text{sink}}}\left[-{\mathit{\gamma}}_{0}\phantom{\rule{2.77695pt}{0ex}}\text{d}{A}_{j}\text{Re}\left\{\phantom{\sum _{k=1}^{{n}_{\text{coef}}}}\frac{ln{Z}_{j}\left(z\right)}{2\pi}-1\right\}+\sum _{k=1}^{{n}_{\text{coef}}}{c}_{j,k}\text{Re}\left\{{Z}_{j}^{-k}\left(z\right)\right\}\right]\\ {\displaystyle}& & {\displaystyle}\phantom{\rule{2em}{0ex}}\text{(streamslitelements)}\\ {\displaystyle}& & {\displaystyle}\phantom{\rule{1em}{0ex}}+\phantom{\rule{2.77695pt}{0ex}}\sum _{{j}^{\prime}=1}^{{n}_{\text{bound}}}\sum _{k=1}^{{n}_{\text{coef}}}\text{Re}\left\{{c}_{{j}^{\prime},k}{Z}_{{j}^{\prime}}^{-k}\left(z\right)\right\}\\ {\displaystyle}& & {\displaystyle}\phantom{\rule{2em}{0ex}}\text{(boundaryslitelements)}.\end{array}$$ |

Just as before, the coefficients of the Laurent series c_{j,k} for stream slits are real-valued in order to ensure the continuity of 𝛷′, while the coefficients c_{j ′,k} for boundary slits can remain complex-valued as long as we are not interested in computing a potential outside of the boundary.

The degrees of freedom are better identified with ancillary potential data at groundwater observation boreholes. Figure 14 shows the final map; the regression yields 𝛾_{0} = −1.19 × 10^{−4} m^{−1}.

In Figure 15, we compare the observed versus estimated groundwater level at the 74 observation boreholes. In order to narrow down the range of values we plot the estimated depth, ${\mathit{\Phi}}_{\text{topo}}^{\prime}-{\hat{\mathit{\Phi}}}_{\text{GW}}^{\prime}$ as a function of observed depth ${\mathit{\Phi}}_{\text{topo}}^{\prime}-{\mathit{\Phi}}_{\text{GW}}^{\prime}$. The RMSE is substantially improved compared to the previous method (0.70 m versus 1.63 m).

Interestingly, if we have an independent estimate of recharge rate r_{0}, we can get an estimate of large-scale, single-layer equivalent transmissivity. Various estimates of annual effective rainfall can be found in the literature for the region, ranging from 180 mm [Mougin et al. 2004] to 393 mm [BRGM 2019]. Taking a mean estimate and considering a 50%–50% partition of this amount between surface runoff and recharge, this leads to a rough estimate of r_{0} ≃ 4.5 × 10^{−9} m⋅s^{−1}. We can then estimate the equivalent, uniform transmissivity ${\hat{T}}_{0}$ with:

$${\hat{T}}_{0}=\frac{{r}_{0}}{\left|{\mathit{\gamma}}_{0}\right|}\sim 3.8\times 1{0}^{-5}\phantom{\rule{0.3em}{0ex}}{\text{m}}^{2}\cdot {\text{s}}^{-1}.$$ |

We can compare this conceptual estimate with values given by aquifer tests: Martelat and Lachassagne [1995] give the ranges 0.8–1.5 × 10^{−3} m^{2}⋅s^{−1} for the superficial schist alterite, and 1.5–3.0 × 10^{−5} m^{2}⋅s^{−1} for the deeper, fissured schist. We see that even if it remains conceptual, the value obtained by regression using the AEM method could be an interesting proxy for characterizing the subsurface at large scale.

## 5. Conclusion and perspectives

In this concluding section we sum up our findings by making a parallel between the decomposition of the topographic surface and that of the groundwater level. Figure 16 illustrates the following decomposition:

$$\left\{\begin{array}{cc}{z}_{\text{topo}}={\mathit{\Phi}}_{\text{topo}}^{\prime}={\mathit{\Phi}}_{\text{BEST}}^{\prime}+\underset{\text{HABEST}}{\underbrace{\left({\mathit{\Phi}}_{\text{topo}}^{\prime}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}\right)}}\phantom{\rule{1em}{0ex}}\hfill & \hfill \\ h={\mathit{\Phi}}_{\text{GW}}^{\prime}={\mathit{\Phi}}_{\text{BEST}}^{\prime}+\underset{\text{GWlevelrelativetoBEST}}{\underbrace{\left({\mathit{\Phi}}_{\text{GW}}^{\prime}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}\right)}}\phantom{\rule{1em}{0ex}}\hfill & \hfill \end{array}\right.$$ |

As can be seen on the figure, the topographic height above the envelope surface of talwegs (${\mathit{\Phi}}_{\text{topo}}^{\prime}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}=\text{HABEST}$) and the relative groundwater level above this envelope surface (${\mathit{\Phi}}_{\text{GW}}^{\prime}-{\mathit{\Phi}}_{\text{BEST}}^{\prime}$) have a very similar behaviour. It comes as no surprise for geomorphologists as, in the hillslope domain, both sediment transfer at the surface and groundwater transfer in the porous medium are ruled by diffusive processes [Anderson and Anderson 2010]. It confirms that the Basal Envelope Surface of Talwegs is an interesting energy level for hydrological and geomorphological analysis.

From the point of view of catchment hydrological modeling, the kind of parsimonious and analytical reconstruction of the groundwater surface proposed in this study, as well as the identification of a conceptual, large-scale estimate of transmissivity ${\hat{T}}_{0}$, could provide useful ancillary information for constraining the calibration of lumped models, with a rationale similar to that of, e.g., TOPMODEL [Beven and Kirkby 1979]. Only a handful of observation boreholes are needed: if the observations are available during both dry and wet seasons, the approach could also yield a large-scale estimate of storage coefficient, $\hat{S}$. The two parameters ${\hat{T}}_{0}$ and $\hat{S}$ can yield first-order estimates of lumped groundwater storage and fluxes which can easily be used as state variables in the model.

Codes used in this study are available in the author’s Basilisk sandbox at http://basilisk.fr/sandbox/nlemoine/AEM/envelope/envelope.c.

## Conflicts of interest

The author has no conflict of interest to declare.

## Acknowledgments

I would like to take the opportunity of this special issue to thank Ghislain de Marsily for all his incredible teaching, especially for passing on to students a taste for analytic approaches.