On the computation of the Basal Envelope Surface of Talwegs using the Analytic Element Method

Nicolas Le Moine

doi:10.5802/crgeos.164

On the computation of the Basal Envelope Surface of Talwegs using the Analytic Element Method

Nicolas Le Moine ¹

¹ UMR 7619 Metis, Sorbonne Université/CNRS/EPHE, 4 Place Jussieu, 75252 Paris cedex 05, France

Comptes Rendus. Géoscience, Volume 355 (2023) no. S1, pp. 79-97.

Résumé

The stream network is a major feature of a landscape, conveying water, sediment, and solute from hillslopes to the ocean. Noticeably, from a large-scale point of view, the elevation of the talwegs of perennial streams is an important head boundary condition for both surface and groundwater flow originating from hillslopes. Assuming a wireframe (1D) representation of talweg lines, the problem of interpolating elevation between talwegs has received attention for applications such as flood mapping using Height Above Nearest Drainage (HAND, [Nobre et al., 2011]), or groundwater level interpolation in low-conductivity aquifer systems. In this study we propose an alternate definition of this large-scale base level concept introduced by [Wyns et al., 2004], namely the Basal Envelope Surface of Talwegs (BEST) and the associated Height Above the Basal Envelope Surface of Talwegs (HABEST), along with a procedure to compute it using the Analytic Element Method (AEM). It can be defined as the head distribution satisfying Laplace equation (Darcy flow with vanishing divergence), with stream segments set as Dirichlet boundary conditions. The BEST is thus the real part of a complex analytic (holomorphic) function which can be modeled using analytic slit elements, with very low computational requirements and without the need for kriging, as it is often seen in the literature. This analytic model is extended to the case of a non-zero, uniform divergence flow (head distribution satisfying a Poisson equation) which can be useful to analyse groundwater levels at catchment scale.

Métadonnées

Reçu le : 2022-04-01
Révisé le : 2022-08-01
Accepté le : 2022-09-14
Première publication : 2022-10-21
Publié le : 2024-02-23

DOI : 10.5802/crgeos.164

Mots clés : Analytic Element Method, Envelope Surface, Talwegs, Geomorphology, Groundwater, Hillslope, Diffusive processes

Affiliations des auteurs :

Nicolas Le Moine ¹

¹ UMR 7619 Metis, Sorbonne Université/CNRS/EPHE, 4 Place Jussieu, 75252 Paris cedex 05, France

Licence :

CC-BY 4.0

Droits d'auteur : Les auteurs conservent leurs droits

@article{CRGEOS_2023__355_S1_79_0,
     author = {Nicolas Le Moine},
     title = {On the computation of the {Basal} {Envelope} {Surface} of {Talwegs} using the {Analytic} {Element} {Method}},
     journal = {Comptes Rendus. G\'eoscience},
     pages = {79--97},
     publisher = {Acad\'emie des sciences, Paris},
     volume = {355},
     number = {S1},
     year = {2023},
     doi = {10.5802/crgeos.164},
     language = {en},
}

TY  - JOUR
AU  - Nicolas Le Moine
TI  - On the computation of the Basal Envelope Surface of Talwegs using the Analytic Element Method
JO  - Comptes Rendus. Géoscience
PY  - 2023
SP  - 79
EP  - 97
VL  - 355
IS  - S1
PB  - Académie des sciences, Paris
DO  - 10.5802/crgeos.164
LA  - en
ID  - CRGEOS_2023__355_S1_79_0
ER  -

%0 Journal Article
%A Nicolas Le Moine
%T On the computation of the Basal Envelope Surface of Talwegs using the Analytic Element Method
%J Comptes Rendus. Géoscience
%D 2023
%P 79-97
%V 355
%N S1
%I Académie des sciences, Paris
%R 10.5802/crgeos.164
%G en
%F CRGEOS_2023__355_S1_79_0

Nicolas Le Moine. On the computation of the Basal Envelope Surface of Talwegs using the Analytic Element Method. Comptes Rendus. Géoscience, Volume 355 (2023) no. S1, pp. 79-97. doi : 10.5802/crgeos.164. https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.5802/crgeos.164/

Version originale du texte intégral (Proposez une traduction )

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.

Figure 1.
Illustration of different potentiometric levels on a topographic/groundwater transect, adapted from Wyns et al. [2004] in order to show both the Basal Surface of Talwegs and the Elevation of Nearest Drainage of Nobre et al. [2011].

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

\frac{b}{a} ≲ 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} + \vec{𝛻} \cdot \vec{q} = r with \vec{q} = - T (h) \vec{𝛻} h

(1)

where h is the groundwater head, S the storage coefficient,

\vec{q}

the 2D (vertically-integrated) Darcy flux density, T the transmissivity (m²⋅s⁻¹), and r the (positively defined) recharge rate (m⋅s⁻¹). Hence,

\vec{𝛻} \cdot (T (h) \vec{𝛻} h) = S \frac{\partial h}{\partial t} - r (x, y, t) .

(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₀, hence:

T_{0} 𝛥 h = S \frac{\partial h}{\partial t} - r (x, y, t) .

(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:

\{\begin{matrix} 𝛥 h = \frac{\partial^{2} h}{\partial x^{2}} + \frac{\partial^{2} h}{\partial y^{2}} = 0 \\ h (x, y) = z_{topo} (x, y) & if (x, y) is in the stream network \end{matrix}

(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 $\vec{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.

Figure 2.
(Left) Example of an envelope surface, in the form of the canopy of an umbrella; (right) definition by Häsing [1964].

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.

Figure 3.
Vector representation of the talweg or stream network. It is a collection of elementary segments/slits, each defined by two endpoints which can be positioned in the complex $C$ -plane with their complex affix.

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:

\vec{q} = - \vec{𝛻} 𝛷 .

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:

\vec{q} = - K \cdot (h - B) \vec{𝛻} h .

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

\begin{matrix} \vec{𝛻} 𝛷 = - \vec{q} = + K \cdot (h - B) \vec{𝛻} h \\ \vec{𝛻} 𝛷 = K \cdot (h - B) \vec{𝛻} (h - B) = \frac{K}{2} \vec{𝛻} {(h - B)}^{2} \\ \Rightarrow 𝛷 = \frac{K}{2} {(h - B)}^{2} . \end{matrix}

This expression hence leads to a nonlinear relation between potential 𝛷 (in m³∕s) and head h (in m). In the following, we will linearize the problem and consider a uniform transmissivity T₀ 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₀, the potential 𝛷 is indeed linearly related to groundwater head h:

\vec{𝛻} 𝛷 = - \vec{q} = + T_{0} \vec{𝛻} h

𝛷 = T_{0} h + constant .

(5)

In order to model the envelope surface, we will work in the complex plane

C

. The start- and end-point of each segment j in the vector talweg network, A_j and B_j, are first given a complex affix in the

C

-plane (see Figure 3):

\{\begin{matrix} a_{j} = x_{A_{j}} + i y_{A_{j}} & affix of start point A_{j} of \\ element j in the complex plane \\ b_{j} = x_{B_{j}} + i y_{B_{j}} & affix of end point B_{j} of \\ element j in the complex plane \end{matrix}

(6)

with i, a number such that i² = −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 $\vec{q}$ , the real-valued potential 𝛷 is the real part of a complex-valued potential 𝛺:

𝛺 (z) = 𝛷 (x, y) + i 𝛹 (x, y)

where 𝛹 is the associated, real-valued stream function. The complex potential 𝛺(z), which is defined in the complex plane

C

and has values also in

C

, has a distinct property: it is analytic or holomorphic, which means that the derivative of 𝛺 with respect to the complex variable z is defined (this is not true in general for a function from

C

C

\begin{matrix} {\frac{d 𝛺}{d z}|}_{z_{0} \in C} & = & {lim}_{z \to z_{0}} \frac{𝛺 (z) - 𝛺 (z_{0})}{z - z_{0}} exists: does not vary \\ depending on the direction from which \\ z_{0} is approached . \end{matrix}

Synonymously, it means that the real potential 𝛷 and the stream function 𝛹 satisfy the Cauchy–Riemann conditions:

\{\begin{matrix} \frac{\partial 𝛷}{\partial x} = \frac{\partial 𝛹}{\partial y} & (i) \\ \frac{\partial 𝛷}{\partial y} = - \frac{\partial 𝛹}{\partial x} & (ii) \\ \frac{\partial^{2} 𝛷}{\partial x^{2}} + \frac{\partial^{2} 𝛷}{\partial y^{2}} = 𝛥 𝛷 = 0 & (iii), consequence of (i) and (ii); \\ the actual property we want \\ for 𝛷 \\ \frac{\partial^{2} 𝛹}{\partial x^{2}} + \frac{\partial^{2} 𝛹}{\partial y^{2}} = 𝛥 𝛹 = 0 & (iv), consequence of (i) and (ii) . \end{matrix}

Contours of the stream function define streamlines; the area between two streamlines defines a stream tube, in which the flux is conserved.

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:

𝛺_{j} (z) = Q_{j} \times (\frac{ln Z_{j} (z)}{2 π} - 1) + \sum_{k = 1}^{n_{coef}} c_{j, k} Z_{j}^{- k} (z) .

(7)

Q_j is the net algebraic discharge (in m³∕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{matrix} 𝜁_{j} (z; a_{j}, b_{j}) & = & \frac{z - \frac{1}{2} (a_{j} + b_{j})}{\frac{1}{2} (b_{j} - a_{j})} \end{matrix}

(8)

\begin{matrix} Z_{j} (z) & = & 𝜁_{j} (z) + \sqrt{𝜁_{j} (z) + 1} \sqrt{𝜁_{j} (z) - 1} . \end{matrix}

(9)

The function $Z (𝜁) = 𝜁 + \sqrt{𝜁 + 1} \sqrt{𝜁 - 1}$ is the Joukowsky transformation [Joukowsky 1910]. Figure 4 illustrates the two steps of the transformation $Z_{j} (z) = 𝜁_{j} (z) + \sqrt{𝜁_{j} (z) + 1} \sqrt{𝜁_{j} (z) - 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 $R \times 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.

Figure 4.
Illustration of the transformation Z_j(z) associated to slit element j, with endpoints A_j and B_j. It is obtained by the composition of a similarity 𝜁_j(z), and the Joukowsky transformation $Z_{j} (z) = 𝜁_{j} (z) + \sqrt{𝜁_{j} (z) + 1} \sqrt{𝜁_{j} (z) - 1}$ . Both transformations are 2D conformal mappings or holomorphisms: they preserve angles locally (see red triangle), hence orthogonality (see grid).

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:

𝛺_{tot} (z) = 𝛷_{tot} (z) + i 𝛹_{tot} (z) = 𝛷_{0} - \bar{q_{0}} z + \sum_{j = 1}^{n_{slit}} 𝛺_{j} (z)

(10)

where q₀ = 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).

\bar{q_{0}} = q_{0, x} - i q_{0, y}

is its complex conjugate, and 𝛷₀ 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 𝛷₀, 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:

Figure 5.
General form of the potential created by a slit element. In each figure, the gray levels show the value of 𝛷 = Re{𝛺(z)} while the dashed white lines are contours of the stream function 𝛹 = Im{𝛺(z)}. Black arrows show the flux vector $q = - \bar{d 𝛺 ∕ d z}$ . (Left) Potential created by a pure sink. (Center) Potential created by a single, real-valued coefficient in the Laurent series; note that the real potential 𝛷 is continuous between the “top” and “bottom” of the slit, while the stream function 𝛹 is discontinuous. (Right) Potential created by a single, pure imaginary coefficient in the Laurent series; there is now a jump in the real potential 𝛷 between the two sides of the slit while the stream function, in turn, is continuous.

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₀ = (Q∕2πT₀)lnr + constant [e.g., De Marsily 1986],
the potential 𝛺(z) = c₁ Z⁻¹(z) with a real-valued c₁ coefficient. In this case, the real potential 𝛷 (gray levels) is continuous,
the potential 𝛺(z) = c₁ Z⁻¹(z) with a pure imaginary c₁ 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₀:

\begin{matrix} 𝛺_{tot}^{'} (z) & = & \frac{1}{T_{0}} 𝛷_{tot} (z) + i \frac{1}{T_{0}} 𝛹_{tot} (z) = 𝛷_{0}^{'} - \bar{q_{0}^{'}} z \\ + \sum_{j = 1}^{n_{slit}} [Q_{j}^{'} \times (\frac{ln Z_{j} (z)}{2 π} - 1) + \sum_{k = 1}^{n_{coef}} c_{j, k}^{'} Z_{j}^{- k} (z)] \end{matrix}

(11)

with the normalized quantities:

\{\begin{matrix} Q_{j}^{'} = \frac{Q_{j}}{T_{0}} & normalized discharge extracted by \\ slit j (depth or drawdown in m) \\ 𝛷_{0}^{'} = \frac{𝛷_{0}}{T_{0}} \\ q_{0}^{'} = \frac{1}{T_{0}} q_{0} \\ c_{j, k}^{'} = \frac{1}{T_{0}} c_{j, k} . \end{matrix}

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_{tot} (z_{p}) = Re {𝛺_{tot}^{'} (z_{p})}$ 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} = Re \{\frac{\partial 𝛺_{tot}^{'} (z_{p})}{\partial 𝛽_{l}}\}

where 𝛽_l is the l-th real degree of freedom. The entry of the matrix on row p and column l reads as:

\{\begin{matrix} J_{p, l} = 1 & if 𝛽_{l} = 𝛷_{0}^{'} \\ J_{p, l} = - Re {z_{p}} & if 𝛽_{l} = q_{0, x}^{'} \\ J_{p, l} = + Im {z_{p}} & if 𝛽_{l} = q_{0, y}^{'} \\ J_{p, l} = Re \{\frac{ln Z_{j} (z_{p})}{2 π}\} - 1 & if 𝛽_{l} = Q_{j}^{'} \\ J_{p, l} = + Re {Z_{j}^{- k} (z_{p})} & if 𝛽_{l} = Re {c_{j, k}^{'}} \\ J_{p, l} = - Im {Z_{j}^{- k} (z_{p})} & if 𝛽_{l} = Im {c_{j, k}^{'}} . \end{matrix}

(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:

J {[𝛽_{1}, 𝛽_{2}, \dots, 𝛽_{N}]}^{T} = {[𝛷_{1}^{'}, 𝛷_{2}^{'}, \dots, 𝛷_{M}^{'}]}^{T} .

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

q = q_{x} + i q_{y} = - [\frac{\partial 𝛷}{\partial x} + i \frac{\partial 𝛷}{\partial y}] = - \bar{\frac{d 𝛺}{d z}} .

Using the chain rule

\frac{d 𝛺_{j}}{d z} = \frac{d 𝛺_{j}}{d Z_{j}} \frac{d Z_{j}}{d 𝜁_{j}} \frac{d 𝜁_{j}}{d z},

the contribution of each slit to the total flux hence reads:

\begin{matrix} q_{j} (z) & = & - \bar{\frac{d 𝛺_{j}}{d z}} \\ = & \bar{\frac{2}{(b_{j} - a_{j}) \sqrt{𝜁_{j} (z) + 1} \sqrt{𝜁_{j} (z) - 1}} [- \frac{Q_{j}}{2 π} + \sum_{k = 1}^{n_{coef}} c_{j, k} k Z_{j}^{- k} (z)]} \\ q_{j} (z) & = & \bar{\frac{2}{(b_{j} - a_{j}) \sqrt{𝜁_{j} (z) + 1} \sqrt{𝜁_{j} (z) - 1}}} \\ \times [- \frac{Q_{j}}{2 π} + \sum_{k = 1}^{n_{coef}} (Re {c_{j, k}} - i Im {c_{j, k}}) k \bar{Z_{j}^{- k} (z)}] \end{matrix}

or, in the normalized form:

\begin{matrix} q_{j}^{'} (z) & = & \bar{\frac{2}{(b_{j} - a_{j}) \sqrt{𝜁_{j} (z) + 1} \sqrt{𝜁_{j} (z) - 1}}} \\ \times [- \frac{Q_{j}^{'}}{2 π} + \sum_{k = 1}^{n_{coef}} (Re {c_{j, k}^{'}} - i Im {c_{j, k}^{'}}) k \bar{Z_{j}^{- k} (z)}] . \end{matrix}

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² [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.

Figure 6.
Map of the Upper Evel catchment and the nested Coët-Dan subcatchment. The vector stream network is extracted with a support area of 0.1 km².

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.

Type of elements	Number of elements	Number of sink terms $Q_{j}^{'}$	Number of $Re {c_{j, k}^{'}}$	Number of $Im {c_{j, k}^{'}}$	Number of degrees of freedom
Stream segm.	1562	0	2 × 1562	0	3124
Boundary segm.	176	0	2 × 176	2 × 176	704
				Total incl. $𝛷_{0}^{'}$	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 𝛷₀ 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:

𝜁_{m} = cos (π \frac{m - \frac{1}{2}}{M}) (m = 1, 2, \dots, M) .

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 $𝛷_{BEST}^{'}$ 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.

Figure 7.
Elevation map of the Basal Envelope Surface of Talwegs (BEST) resulting from the least squares regression against talweg elevations, with a Neumann condition on the boundary.

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 $𝛷_{HABEST}^{'} = (𝛷_{topo}^{'} - 𝛷_{BEST}^{'}) = (z_{topo} - 𝛷_{BEST}^{'})$ . Contours of this surface can be used for flood mapping in the same way as HAND contours.

Figure 8.
Map of the Height Above the Basal Envelope Surface of Talwegs (HABEST) index, defined as the difference between topographic elevation $z_{topo} = 𝛷_{topo}^{'}$ and the elevation $𝛷_{BEST}^{'}$ ; this index is closely related to the HAND index.

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:

𝛷_{GW}^{'} - 𝛷_{BEST}^{'} = a (𝛷_{topo}^{'} - 𝛷_{BEST}^{'}) + b

(13)

with

𝛷_{GW}^{'} = h

the groundwater head (or real potential normalized by T₀), and

𝛷_{topo}^{'} = z_{topo}

the topographic elevation.

Figure 9.
Close-up on the elevation map of the Basal Envelope Surface of Talwegs in the Coët Dan subcatchment (dashed white line). Black squares indicate the location of the 74 boreholes in which the groundwater level was surveyed during the month of February, 1996.

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

\{\begin{matrix} a = 0.912 \\ b = - 1.74 m \end{matrix}

with R² = 0.962 and a rooted mean squared error (RMSE) of 1.67 m.

Figure 10.
Result of the regression of relative groundwater above the Basal Envelope Surface of Talwegs, $(𝛷_{GW}^{'} - 𝛷_{BEST}^{'})$ , against the Height Above the Basal Envelope Surface of Talwegs $(𝛷_{topo}^{'} - 𝛷_{BEST}^{'})$ .

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

Figure 11.
Rooted-Mean-Squared-Error (RMSE) of the regression of $(𝛷_{GW}^{'} - 𝛷_{BEST}^{'})$ against $(𝛷_{topo}^{'} - 𝛷_{BEST}^{'})$ , as a function of the support area used for extracting the talweg network from the DEM. We see a slight dependence, with an optimal value of the support area about A_c = 0.4 km².

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} 𝛥 h = T_{0} 𝛥 𝛷^{'} = - r_{0}

where r₀ > 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:

𝛥 𝛷^{'} = - \frac{r_{0}}{T_{0}} = + 𝛾_{0}

where the sign of 𝛾₀ is set according to the convention used in Strack [2017]. Just as we normalized

Q_{j}^{'} = Q_{j} ∕ T_{0}

for a slit sink term, 𝛾₀ = −r₀∕T₀ is a normalized areal extraction rate term which has the dimension of the inverse of a depth (e.g., m⁻¹).

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 𝛾₀. 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} d A_{j}

or, equivalently,

Q_{j}^{'} = \frac{r_{0}}{T_{0}} d A_{j} = - 𝛾_{0} d A_{j} .

Figure 12.
For each stream slit j, we compute dA_j the area of the subcatchment drained by the slit, that is, the difference between the value of the flow accumulation map at the downstream endpoint of the slit minus the value at the upstream endpoint. For “leaf” slits, dA_j is simply set to the drained areal at the downstream end of the slit. Each slit is then allowed to extract a normalized discharge $Q_{j}^{'} = (r_{0} ∕ T_{0}) d A_{j} = - 𝛾_{0} 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 𝛾₀. 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 𝛾₀ = 1.

Figure 13.
Contour map of real potential created by an areal sink of unit normalized extraction rate 𝛾₀ = −r₀∕T₀ = 1. The Laplacian of the potential uniformly equals 𝛾₀ inside the polygon (white dashed line), and it is harmonic outside ( $𝛥 𝛷_{areal}^{'} = 0$ ). Note the very large magnitude of the values (in the order of 10⁸): of course, it has to be scaled by a much smaller factor |𝛾₀|≪1 in order to model topographic or groundwater elevation.

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

\begin{matrix} 𝛷_{GW}^{'} (z) = 𝛷_{0}^{'} + Re {\bar{v_{0}^{'}} z} + 𝛾_{0} 𝛷_{areal}^{'} (z; 𝛾_{0} = 1) \\ + \sum_{j = 1}^{n_{sink}} [- 𝛾_{0} d A_{j} Re \{\frac{ln Z_{j} (z)}{2 π} - 1\} + \sum_{k = 1}^{n_{coef}} c_{j, k} Re {Z_{j}^{- k} (z)}] \\ (stream slit elements) \\ + \sum_{j^{'} = 1}^{n_{bound}} \sum_{k = 1}^{n_{coef}} Re {c_{j^{'}, k} Z_{j^{'}}^{- k} (z)} \\ (boundary slit elements) . \end{matrix}

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 𝛾₀ = −1.19 × 10⁻⁴ m⁻¹.

Figure 14.
Simulated map of the groundwater potential $𝛷_{GW}^{'}$ in the Coët-Dan subcatchment, with the additional degree of freedom 𝛾₀ (normalized areal exfiltration rate). The regression yields 𝛾₀ = −1.19 × 10⁻⁴ m⁻¹ (net infiltration). The error at the 74 observation borehole is represented with graduated triangles.

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, $𝛷_{topo}^{'} - {\hat{𝛷}}_{GW}^{'}$ as a function of observed depth $𝛷_{topo}^{'} - 𝛷_{GW}^{'}$ . The RMSE is substantially improved compared to the previous method (0.70 m versus 1.63 m).

Figure 15.
Scatter plot of observed versus estimated depth to groundwater level with the additional degree of freedom 𝛾₀ (normalized areal exfiltration rate).

Interestingly, if we have an independent estimate of recharge rate r₀, 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₀ ≃ 4.5 × 10⁻⁹ m⋅s⁻¹. We can then estimate the equivalent, uniform transmissivity ${\hat{T}}_{0}$ with:

{\hat{T}}_{0} = \frac{r_{0}}{| 𝛾_{0} |} \sim 3.8 \times 1 0^{- 5} m^{2} \cdot 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⁻³ m²⋅s⁻¹ for the superficial schist alterite, and 1.5–3.0 × 10⁻⁵ m²⋅s⁻¹ 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:

\{\begin{matrix} z_{topo} = 𝛷_{topo}^{'} = 𝛷_{BEST}^{'} + \underset{HABEST}{\underset{︸}{(𝛷_{topo}^{'} - 𝛷_{BEST}^{'})}} \\ h = 𝛷_{GW}^{'} = 𝛷_{BEST}^{'} + \underset{GW level relative to BEST}{\underset{︸}{(𝛷_{GW}^{'} - 𝛷_{BEST}^{'})}} \end{matrix}

Figure 16.
Decomposition of topographic elevation ( $𝛷_{topo}^{'}$ , top row) and groundwater level elevation ( $𝛷_{GW}^{'}$ , bottom row) using the ancillary level $𝛷_{BEST}^{'}$ . We see that in both cases this level plays the role of a base level for diffusive processes in the hillslope domain (diffusive sediment transport on the topographic surface, and diffusive Darcy flow for groundwater).

As can be seen on the figure, the topographic height above the envelope surface of talwegs ( $𝛷_{topo}^{'} - 𝛷_{BEST}^{'} = HABEST$ ) and the relative groundwater level above this envelope surface ( $𝛷_{GW}^{'} - 𝛷_{BEST}^{'}$ ) 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.

Bibliographie

[Anderson and Anderson, 2010] R. Anderson; S. Anderson Geomorphology: The Mechanics and Chemistry of Landscapes, Cambridge University Press, Cambridge, 2010 | DOI

[Beven and Kirkby, 1979] K. J. Beven; M. J. Kirkby A physically based, variable contributing area model of basin hydrology, Hydrol. Sci. Bull., Volume 24 (1979) no. 1, pp. 43-69 | DOI

[BRGM, 2019] BRGM 195AC01 – Socle métamorphique dans le bassin versant de l’Evel de sa source au Blavet (nc), 2019 https://sigesbre.brgm.fr/files/fiches/BDLISA/LISA_Bretagne_195AC01.pdf (BD LISA database, accessed: 2022-03-30)

[Cohen et al., 2015] Y. Cohen; O. Devauchelle; H. F. Seybold; R. S. Yi; P. Szymczak; D. H. Rothman Path selection in the growth of rivers, Proc. Natl. Acad. Sci. USA, Volume 112 (2015) no. 46, pp. 14132-14137 https://www.pnas.org/content/112/46/14132

[Deffontaines, 1990] B. Deffontaines Développement d’une méthodologie d’analyse morphostructurale et morphonéotectonique. Analyse des surfaces enveloppes du reseau hydrographique et des modeles numeriques de terrain; application au nord-est de la France, Ph. D. Thesis, Université Pierre et Marie Curie (1990)

[De Marsily, 1986] G. De Marsily Quantitative Hydrogeology: Groundwater Hydrology for Engineers, Elsevier Science, Amsterdam, 1986

[Douglas and Peucker, 1973] D. H. Douglas; T. K. Peucker Algorithms for the reduction of number of points required to represent a digitized line or its caricature, Cartogr. Int. J. Geogr. Inform. Geovisuali., Volume 10 (1973), pp. 112-122

[Dunne and Black, 1970] T. Dunne; R. Black Partial area contribution to storm runoff in a small New England watershed, Water Resour. Res., Volume 6 (1970), pp. 1296-1311

[Dunne, 1969] T. Dunne Runoff production in a humid area, Ph. D. Thesis, Johns Hopkins University, Baltimore, MD (1969)

[Dunne, 1980] T. Dunne Formation and controls of channel networks, Prog. Phys. Geogr., Volume 4 (1980), pp. 211-239

[ESRI, 2022] ESRI ArcGIS Pro geoprocessing tool reference, 2022 https://pro.arcgis.com/en/pro-app/2.8/tool-reference/spatial-analyst/stream-link.htm (Accessed: 2022-02-21)

[Fovet et al., 2018] O. Fovet; L. Ruiz; G. Gruau; N. Akkal; L. Aquilina; S. Busnot; R. Dupas; P. Durand; M. Faucheux; Y. Fauvel; C. Fléchard; N. Gilliet; C. Grimaldi; Y. Hamon; A. Jaffrezic; L. Jeanneau; T. Labasque; G. Le Henaff; P. Mérot; J. Molénat; P. Petitjean; A.-C. Pierson-Wickmann; H. Squividant; V. Viaud; C. Walter; C. Gascuel-Odoux Agrhys: An observatory of response times in agro-hydro systems, Vadose Zone J., Volume 17 (2018) no. 1, 180066 | DOI

[Häsing, 1964] J. Häsing Bestimmung der Glättungstiefe rauher Flächen, PTB-Mitt., Volume 4 (1964), pp. 339-340

[Institut Géographique National, 2018] Institut Géographique National RGE ALTI $^{®}$ version 2.0 – Descriptif de contenu, 2018 https://geoservices.ign.fr/sites/default/files/2021-07/DC_RGEALTI_2-0.pdf

[Janković and Barnes, 1999] I. Janković; R. Barnes High-order line elements in modeling two-dimensional groundwater flow, J. Hydrol., Volume 226 (1999) no. 3, pp. 211-223 https://www.sciencedirect.com/science/article/pii/S0022169499001407

[Joukowsky, 1910] N. Joukowsky Über die Konturen der Tragflächen der Drachenflieger, Z. Flugtechnik Motorluftschiffahrt, Volume 1 (1910) no. 22, pp. 281-285

[Martelat and Lachassagne, 1995] A. Martelat; P. Lachassagne Bassin Versant Représentatif Expérimental du Coët Dan (Naizin, Morbihan). Hydrogéologie : détermination des caractéristiques hydrodynamiques du système aquifère au lieu-dit Le Stimoës (1995) no. BRGM/RR-38474-FR (Technical report)

[Mougin et al., 2004] B. Mougin; E. Thomas; R. Wyns; R. Blanchin; F. Mathieu Qualité des eaux en Bretagne – Ruissellement – Infiltration – Temps de réponse – Bassin versants du Yar (22), de l’Horn (29), et du Coët Dan (56) – Rapport final (2004) no. BRGM/RP-52731-FR (Technical report)

[Nobre et al., 2011] A. Nobre; L. Cuartas; M. Hodnett; C. Rennó; G. Rodrigues; A. Silveira; M. Waterloo; S. Saleska Height above the nearest drainage – a hydrologically relevant new terrain model, J. Hydrol., Volume 404 (2011) no. 1, pp. 13-29 https://www.sciencedirect.com/science/article/pii/S0022169411002599

[Nobre et al., 2016] A. D. Nobre; L. A. Cuartas; M. R. Momo; D. L. Severo; A. Pinheiro; C. A. Nobre HAND contour: a new proxy predictor of inundation extent, Hydrol. Process., Volume 30 (2016) no. 2, pp. 320-333 | DOI

[Pauwels et al., 1996] H. Pauwels; A. Martelat; J. C. Foucher; P. Lachassagne Dénitrification dans les eaux souterraines du bassin du Coët-Dan : Suivi géochimique et hydrogéologique du processus (1996) no. BRGM/RR-39055-FR (Technical report)

[Petroff et al., 2012] A. P. Petroff; O. Devauchelle; A. Kudrolli; D. H. Rothman Four remarks on the growth of channel networks, C. R. Géosci., Volume 344 (2012) no. 1, pp. 33-40 https://www.sciencedirect.com/science/article/pii/S1631071312000028

[Steward, 2015] D. R. Steward Analysis of discontinuities across thin inhomogeneities, groundwater/surface water interactions in river networks, and circulation about slender bodies using slit elements in the Analytic Element Method, Water Resour. Res., Volume 51 (2015) no. 11, pp. 8684-8703 | DOI

[Steward, 2020] D. R. Steward Analytic Element Method: Complex Interactions of Boundaries and Interfaces, Oxford University Press, Oxford, 2020 | MR

[Strack, 1989] O. D. Strack Groundwater Mechanics, Prentice Hall, Hoboker, NJ, 1989

[Strack, 2017] O. D. Strack Analytical Groundwater Mechanics, Cambridge University Press, Cambridge, 2017

[Système d’Information sur l’Eau/Sandre, 2020] Système d’Information sur l’Eau/Sandre Document de présentation – Description du référentiel hydrographique (BD TOPAGE $^{®}$ ) version 1, 2020 http://www.sandre.eaufrance.fr/ftp/documents/fr/pre/eth/1/sandre_pres_eth_1.pdf

[Wyns et al., 2004] R. Wyns; J.-M. Baltassat; P. Lachassagne; A. Legchenko; J. Vairon; F. Mathieu Application of proton magnetic resonance soundings to groundwater reserve mapping in weathered basement rocks (Brittany, France), Bull. Soc. Géol. Fr., Volume 175 (2004) no. 1, pp. 21-34 | DOI

[Yao et al., 2017] T. K. Yao; O. Fouché; K. E. Kouadio; M.-S. Oga Discontinuous nature of phreatic aquifers in granitic rocks at watershed scale: A stratiform model from perennial streams and well data, Asian Rev. Environ. Earth Sci., Volume 4 (2017) no. 1, pp. 20-27

Commentaires - Politique

Ces articles pourraient vous intéresser

A hydrogeological acrostic: in honour of Ghislain de Marsily

Vazken Andréassian; Valérie Plagnes; Craig Simmons; ...

C. R. Géos (2023)

Astrometry of directly imaged exoplanets with optical interferometry

Sylvestre Lacour

C. R. Phys (2023)

Caractérisation des profils de formations superficielles par pénétrométrie dynamique à énergie variable : application aux marnes noires de Draix (Alpes-de-Haute-Provence, France)

Olivier Maquaire; Arnaud Ritzenthaler; Denis Fabre; ...

C. R. Géos (2002)