# Flow and residence time in a two-dimensional aquifer recharged by rainfall

Valentin Jules, Eric Lajeunesse, Olivier Devauchelle, Adrien Guérin, Claude

Jaupart, Pierre-Yves Lagrée

► To cite this version:

Valentin Jules, Eric Lajeunesse, Olivier Devauchelle, Adrien Guérin, Claude Jaupart, et al.. Flow and residence time in a two-dimensional aquifer recharged by rainfall. Journal of Fluid Mechanics, 2021, 917, pp.A13. 10.1017/jfm.2021.221. hal-03207646

HAL Id: hal-03207646

https://hal.science/hal-03207646

Submitted on 26 Apr 2021

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
# Flow and residence time in a two-dimensional aquifer recharged by rainfall

V. Jules<fnref n="1" />, E. Lajeunesse<fnref n="1" />†, O. Devauchelle<fnref n="1" />, A. Guérin<fnref n="2" />, C. Jaupart<fnref n="1" /> and P-Y. Lagrée<fnref n="3" />

(Received xx; revised xx; accepted xx)

We investigate the flow of water in a two-dimensional laboratory aquifer recharged by artificial rainfall. As rainwater infiltrates, it forms a body of groundwater which can exit the aquifer only through one of its sides. The outlet, located high above the aquifer bottom, drives the flow upwards. Noting that the water table barely departs from the horizontal, we linearize the boundary condition at the free surface, and solve the flow equations in steady state. We find an approximate expression for the velocity potential, which accounts for the shape of the streamlines, and for the propagation of dye through the aquifer. Based on this theory, we calculate the travel time of water through the experiment, and find that its distribution decays exponentially, with a characteristic time that depends on the shape of the aquifer. We find that the hydrodynamic dispersion that occurs at the pore scale does not matter much for this distribution, which depends essentially on the geometry of the groundwater flow.

## 1. Introduction

An aquifer is a porous rock or deposit saturated with water. In most geological settings, groundwater is in mechanical contact with the atmosphere through the pores of the overlying soil. The aquifer is then said to be “unconfined”, and the free surface of the groundwater body is called the “water table” (Ingebritsen & Sanford 1999). Unconfined aquifers are recharged by precipitations. Rainwater infiltrates into the soil, and percolates through the unsaturated zone until it reaches the water table (Kirchner et al. 2000). This groundwater then flows through the aquifer, and eventually emerges into streams to form surface runoff.

The time rainfall needs to travel through the aquifer controls the chemical weathering of the porous rock, and the propagation of dissolved contaminants towards the drainage network (Kirchner et al. 2000; Maher 2010; Godsey et al. 2010; Rempe & Dietrich 2014). As it travels through the aquifer, a parcel of groundwater gradually spreads under the combined action of three processes: (1) molecular diffusion, which results from thermal fluctuations at the molecular scale; (2) hydrodynamic dispersion, which results from velocity fluctuations at the pore scale; and (3) advection by the average velocity of groundwater, which varies at the scale of the aquifer (Saffman 1959; De Marsily 1986; Charlaix et al. 1987). The travel time of water through the aquifer results from the combination of these processes, which, in its turn, depends on the geometry and scale of the groundwater flow (Le Borgne et al. 2007; Dentz et al. 2011).

Field measurements, based on isotopic ratios or conservative tracers, reveal that the mean travel time of water typically ranges from a few months to millennia (Gleeson et al. 2016;

<footnote n="1">¹Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France</footnote>

<footnote n="2">²MSC, Université Paris Diderot, CNRS (UMR 7057), 75013 Paris, France</footnote>

<footnote n="3">³Sorbonne Université, CNRS - UMR 7190, Institut Jean Le Rond d'Alembert, F-75005 Paris, France</footnote>

<footnote>† Email address for correspondence: lajeunes@ipgp.fr</footnote>
Cartwright *et al.* 2017). Over such a long time, runoff balances recharge, and one may consider that the flow is in steady state (Cardenas 2007).

In nature, most aquifers extend much further horizontally than vertically, which makes the groundwater flow amenable to the shallow-water approximation. Combined with Darcy's law (Darcy 1856), the latter provides a simplified description of the flow, in the form of the Dupuit-Boussinesq equation (Dupuit 1863; Boussinesq 1903). This theory, which neglects the vertical component of the flow, consistently accounts for the fluctuations of the water table caused by intermittent rainfall (Brutsaert & Nieber 1977; Troch *et al.* 2013; Guérin *et al.* 2014, 2019). When averaged over a sufficiently long time, however, these fluctuations vanish, and the average water table conforms to the steady-state Dupuit-Boussinesq equation (Guérin *et al.* 2019). On average, one then finds that the water table mirrors the land's surface, with summits far from rivers, and valleys that line the drainage network (Petroff *et al.* 2011; Devauchelle *et al.* 2012).

In reality, however, the stream that drains a catchment often lies above the base of the aquifer. Groundwater must therefore rise to reach the outlet (Lehr 1963), and the vertical component of the resulting flow violates the Dupuit-Boussinesq approximation (Haria & Shand 2004; Powers *et al.* 1967; Read 1993; Bresciani *et al.* 2014). In this configuration, the velocity potential that drives the flow obeys the Laplace equation, the solution of which crucially depends on the boundary conditions (Polubarinova-Kochina 1962). One often assumes that the flow is bounded by a horizontal impervious layer at its base, while impervious boundaries represent the water divides on its sides (Toth 1963). The groundwater pressure must also match the atmospheric one along the water table. This last boundary condition, however, is more challenging than the others: the position of the water table is not specified a priori, and needs to be solved for.

Turning the similarity between the water table and the land surface into a boundary condition, Toth (1963) assumed that the water-table elevation is proportional to that of the land surface, and set the groundwater pressure accordingly. Mathematically, this is a Dirichlet condition: the value of the potential is specified at the boundary (Farlow 1993). The resulting flow is then entirely determined by the land's topography: groundwater sinks below hills, and rises towards valleys (Toth 1963; Cardenas 2007).

This hypothesis, although leading to realistic-looking flows, is nevertheless unjustified. In reality, the water table freely adjusts its shape to accommodate the rainfall input, and there is no reason to believe that, in doing so, it follows the shape of the land's surface exactly. Instead of specifying the water table a priori, one must solve the flow equations to get it (Polubarinova-Kochina 1962), and this free-boundary problem requires not one, but two boundary conditions at the water table. In addition to setting the groundwater pressure, one must also account for the flux of rainwater that reaches the water table (Powers *et al.* 1967). Mathematically, this is a Stefan problem: the values of both the potential and its derivative are specified at the boundary, the position of which must be solved for (Rubinstein 2000). One thus needs to solve the Laplace equation inside a domain that is a priori unknown.

How does the aquifer's geometry influence such a flow? Is it possible to infer this geometry based on the residence time of water in the aquifer? Inspired by Guérin *et al.* (2014), we address these questions in a laboratory experiment. We begin with the description of the experimental setup (§2). We then focus on the steady-state regime (§3). We derive an analytical solution for this regime and compare it to experimental streamlines (§4). We then compute the travel time of water, and show how the aquifer's aspect ratio affects it (§5). Finally, we estimate the magnitude of hydrodynamic dispersion in our experiment (§6).
<figure><img src="image_2.png" /></figure>

Figure 1: Experimental setup and notations. The aquifer is made of glass beads (shaded area) piled between two vertical glass plates (black solid lines) separated by a gap of width $W = 5$ cm. Its length and depth are $L = 90$ cm and $H = 62$ cm, respectively. A series of 24 capillary tubes (diameter 1.2 mm, length 50 cm) placed above the tank simulate rainfall, by distributing water uniformly over the aquifer's surface.

## 2. Laboratory aquifer

### 2.1. *Setup*

At its simplest, an unconfined aquifer consists of a homogeneous porous medium, the upper surface of which is in contact with the atmosphere. Inspired by Guérin *et al.* (2014), we reproduce this configuration in a narrow tank formed by two vertical plexiglas plates (90 × 92 cm), separated by a gap of width $W = 5$ cm (figure 1). This tank is filled with glass beads of size either 1 or 3 mm, which form an artificial aquifer bounded by an impermeable bottom. An impervious vertical wall closes the left side of the tank, mimicking the drainage divide of a catchment. We generate an artificial rain with an array of 24 capillary tubes evenly spaced above the tank. These tubes (diameter 1.2 mm, length 50 cm) uniformly distribute water over the aquifer's surface, with a standard deviation of about 3%.

We want to investigate the flow in a deep aquifer. To generate it, we bound the right-hand side of the aquifer with an impervious wall extending up to $H = 62$ cm above the aquifer's bottom (figure 1). The remaining height (30 cm) is covered with a permeable vertical grid, which retains the glass beads, while allowing water to seep out. The aquifer's outlet thus lies high above its bottom. The corresponding aspect ratio is $H/L \approx 0.7$, where $L = 90$ cm is the length of the aquifer. The thinness of the setup ($W/L \approx 0.05$), and its crosswise invariance, ensure that the flow remains virtually two-dimensional.

When the grid that maintains the beads in the aquifer is directly in contact with the atmosphere, surface tension induces a pressure jump, which might affect the pressure field in the aquifer. To mitigate this inconvenience, we spread a soft plastic sheet over the outer side of the grid. This device maintains a thin film of water at atmospheric pressure over the grid surface, thus eliminating the effect of surface tension.

### 2.2. *Method*

Before each experimental run, we fill the aquifer with water until it overflows. After a few tens of minutes, the discharge at the aquifer's outlet has virtually vanished; the water table
<figure><img src="image_3.png" /></figure>

Figure 2: Discharge at the outlet of our deep experimental aquifer (blue line) and of a shallow one (black dashed line, Guérin, 2015). Light color: raw signal, dark color: average over 10 successive data points. Vertical dotted line: beginning of rainfall. Both aquifers are 5 cm wide and 90 cm long; they are made of the same 1 mm glass beads, and are submitted to similar rainfall inputs: $Q_{in} = 7.9 \text{ mL s}^{-1}$ for the deep aquifer and $8.3 \text{ mL s}^{-1}$ for the shallow one.

is then horizontal, and intersects the outlet. Starting from this initial condition, we switch on the artificial rain. An electromagnetic flowmeter (Kobold MIK 0.05-1 L/min, accuracy 2%) measures the rainfall input, $Q_{in}$, from which we deduce the rainfall rate $R = Q_{in}/(LW)$. The latter then remains constant during the entire experiment.

As the rain falls on the aquifer's surface, it infiltrates into the porous medium, and seeps down until it reaches the water table, which deforms to accommodate this input. The asymmetry of the boundaries causes it to bend down towards the outlet (figure 1). As the resulting pressure field drives water towards the outlet, the aquifer's discharge increases, until it eventually reaches steady state, and balances the rainfall input ($Q = RLW$).

A bucket collects the water that exits the aquifer. It lies on a scale which measures its weight every 0.1 s with a precision of 1 mg. The water discharge $Q$ is the derivative of the bucket weight with respect to time. As water drips from the outlet, the discharge delivered to the bucket fluctuates around its average. As a result, the standard deviation of the steady-state discharge is about $0.3 \text{ mL s}^{-1}$, a value that is small compared to the range of discharge explored in our experiments (table 1).

Once the flow is steady, we inject small volumes of blue dye in the aquifer, and track their motion to trace the groundwater streamlines. We shall discuss the shape of these streamlines in section 3. Before we do so, however, we assess how the depth of our aquifer affects its response to rainfall, as compared to the shallow aquifer of Guérin *et al.* (2014).

## 2.3. *Transient*

Guérin *et al.* (2014) used an experiment similar to the present one, but for the position of the outlet which, in their setup, lies right on the bottom of the aquifer ($H = 0$). In the following, we refer to this configuration as the “shallow aquifer” because the water table remains close to the bottom.

For comparison with Guérin *et al.* (2014), we first use the same glass beads (diameter $d_s = 1$ mm) and rainfall rate ($Q_{in} \approx 8 \text{ mL s}^{-1}$). Following the procedure of the previous
<table><caption>Table 1: Experimental parameters. The aquifer, made of 3 mm glass beads, has a hydraulic conductivity $K = 9.1 \pm 0.3$ cm s$^{-1}$, and a permeability $k \approx 9.3 \times 10^{-9}$ m$^2$. The pore size is $d_p \sim \sqrt{k} \approx 0.1$ mm.</caption><thead><tr><th>Run no.</th><th>Rainfall input, Q<sub>in</sub> [mL s<sup>-1</sup>]</th><th>Rainfall rate, R [cm s<sup>-1</sup>]</th><th>Ratio R/K [-]</th></tr></thead><tbody><tr><td>1</td><td>6.3</td><td>1.4 &middot; 10<sup>-2</sup></td><td>1.5 &middot; 10<sup>-3</sup></td></tr><tr><td>2</td><td>12.5</td><td>2.8 &middot; 10<sup>-2</sup></td><td>3.1 &middot; 10<sup>-3</sup></td></tr><tr><td>3</td><td>5.0</td><td>1.1 &middot; 10<sup>-2</sup></td><td>1.2 &middot; 10<sup>-3</sup></td></tr><tr><td>4</td><td>16.7</td><td>3.7 &middot; 10<sup>-2</sup></td><td>4.1 &middot; 10<sup>-3</sup></td></tr><tr><td>5</td><td>10.3</td><td>2.3 &middot; 10<sup>-2</sup></td><td>2.5 &middot; 10<sup>-3</sup></td></tr></tbody></table>

section, we then record the time series of the discharge, often called hydrograph in this context. Less than 5 seconds after the beginning of rainfall, the discharge of our aquifer increases abruptly (figure 2). After a characteristic time of 10 to 15 seconds, it saturates, and reaches steady state. In the shallow aquifer of Guérin *et al.* (2014), the transient is surprisingly longer: it lasts several minutes (figure 2). The geometry of the flow also changes: the deformation of the water table is much smaller in the deep aquifer (a few millimeters) than in the shallow one (ten to twenty centimeters).

The difference between our observations and those of Guérin *et al.* (2014) results from the boundary condition at the outlet. In the shallow aquifer, the flow is virtually parallel to the aquifer's bottom. According to the Dupuit-Boussinesq approximation, the discharge is then proportional to the water depth. The water level must therefore rise significantly before the flow can transport an amount of water that matches the input: it is this swelling that takes time. By contrast, the flow through a deep aquifer is free to develop a vertical component instantly, which allows water to flow through the whole aquifer right from the beginning of the rainfall event. The large volume of groundwater that an increase of pressure can mobilize thus allows the flow to respond instantly.

Unfortunately, the structure of this transient, two-dimensional flow is out of our reach in our experiment. In the next section, we focus instead on the steady state, which we visualize with dye. To this end, we substitute 3 mm glass beads for the 1 mm ones, as they prove more convenient to track the propagation of the dye.

## 3. Streamlines

Once the flow is in steady state, we inject small parcels of blue dye into the aquifer, at different locations along its surface, and monitor their propagation. The dye consists of a mixture of deionised water (90%) and food coloring (Matfer, 10%). Its density (997 kg m$^{-3}$) is close to that of pure water; we therefore assume that it behaves as a passive tracer. To inject it in the aquifer, we use a set of six hypodermic needles, each of them connected to a syringe through a flexible pipe. The syringes are placed on a press (SPLab06, Schenchen Baoding), which allows us to simultaneously inject small volumes of blue dye (3 mL) at six different locations, regularly spaced along the surface of the aquifer (figure 3).

Entrained by the flow, the parcels of dye then travel through the aquifer. A camera facing
<figure><img src="image_4.png" /></figure>

Figure 3: Propagation of parcels of dye through the aquifer during run no 1 (table 1, time indicated on each snapshot). Dye parcels are injected through the surface of the aquifer at distance $l = 87.3, 72.7, 57.6, 42.5, 27.4, \text{ and } 11.4 \text{ cm}$ from the outlet. The aquifer, made of 3 mm glass beads, receives a rainfall input $Q_{in} = 6.3 \text{ mL s}^{-1}$. Black dashed lines: streamlines after equation (4.11). White dots: position of the dye parcels after equation (5.1). Theory has no adjustable parameters.
the setup records an image of their position every 15 seconds (figure 3). A light panel, placed behind the aquifer, increases the contrast of the image.

We estimate the porosity of the aquifer by comparing the weights of a reference volume of beads successively filled with air and water. We find a porosity $s \sim 0.4$, as expected for the random close packing of spheres (Andreotti et al. 2013). Following Lobkovsky et al. (2004), we also measure the hydraulic conductivity of the aquifer using two communicating vessels connected by a Darcy column filled with 3mm glass beads, and saturated with water. The release of a finite volume of water in one of the vessel induces a pressure gradient in the Darcy column, which drives a flow that tends to equilibrate the water levels in the vessels. The discharge then relaxes exponentially. The characteristic time of this relaxation depends on the hydraulic conductivity $K$ of the pack of beads. Using this method, we find $K = 9.1 \pm 0.3 \text{ cm s}^{-1}$, from which we deduce the permeability of our aquifer, $k = vK/g \approx 9.3 \cdot 10^{-9} \text{ m}^2$, where $g$ is the acceleration of gravity, and $v$ the kinematic viscosity of water. Finally, we estimate the pore size from the permeability: $d_p \sim \sqrt{k} \approx 0.1 \text{ mm}$. This value is about ten time smaller than the bead diameter, as observed in previous experiments (Souzy et al. 2020).

In total, we performed five experimental runs, varying the rainfall input between 5 and $16.7 \text{ mL s}^{-1}$ (table 1). We observed the same behavior during each run: at first, each parcel of dye is slowly entrained downwards, until it eventually reaches the deepest point of its trajectory (figure 3 and supplementary material, movie 1). It then accelerates as it rises towards the surface and exits the aquifer through the outlet. As it moves through the aquifer, the parcel barely diffuses in the direction perpendicular to the flow. On the contrary, it rapidly spreads along the streamline, a manifestation of the well-known anisotropy of hydrodynamic dispersion in porous media (Saffman 1959).

The travel time of a dye parcel depends on the position of its injection point (figure 3 and supplementary material, movie 1). In an aquifer made of 3 mm glass beads and receiving a rainfall input $Q = 6.3 \text{ mL s}^{-1}$, a parcel injected 11.4 cm away from the outlet leaves the aquifer within 3 minutes. The same volume of dye, injected 87.3 cm away from the outlet, needs about 3 hours to reach the outlet. The course of a rain drop through the aquifer thus strongly depends on its landing point.

To better illustrate this phenomenon, we merge all the pictures that record the propagation of the dye (about 700 images) into a single one by keeping, for each pixel, its darkest value over the course of the experiment. The resulting image reveals the trajectories of each dye parcel in the form of six blue stripes, the width of which is that of the parcel. Each of the two edges of a stripe corresponds to a streamline (figure 4). The curvature of the streamlines confirms that the flow is two-dimensional. The flow even becomes essentially vertical near the upper surface of the aquifer, where the streamlines head downwards. As groundwater enters deeper in the aquifer, its trajectory gradually curves inwards, and reaches its maximum depth. Groundwater then rises towards the aquifer's outlet, where all the streamlines converge. The vertical component of the flow is specially strong for the streamlines that originate far from the outlet. As a result, the further a streamline starts from the outlet, the deeper its excursion below the surface.

In the next section, we derive an approximate analytical expression for the streamlines, and compare it to our observations.
<figure><img src="image_5.png" /><figcaption>Figure 4: Trajectories of the parcels of dye revealed by merging 700 pictures (figure 3) into one, keeping at each point the darkest pixel. Experimental conditions are those of figure 3. Black dashed lines: streamlines after equation (4.11).</figcaption></figure>

# 4. Hydraulic head and stream function

## 4.1. Flow equations

Our experimental aquifer is made of a homogeneous and isotropic porous material. The Reynolds number of the flow is $Re = Ud_p/v$, where $v$ is the kinematic viscosity of water, and $U$ its typical velocity in the aquifer. Based on the propagation of the dye in our experiment, we find that the latter varies in the range $10^{-2} \text{cm s}^{-1}$ to $4 \cdot 10^{-1} \text{cm s}^{-1}$, so that $Re$ varies between $10^{-2}$ and 0.4. Groundwater thus flows according to the simplest version of Darcy's law:

$$ q = -K \nabla \left( \frac{p}{\rho g} + y \right) \qquad (4.1) $$

where $q$, $\rho$, and $p$, are the flux of water, its density, and pressure, respectively. The horizontal and vertical coordinates, $x$ and $y$, are measured with respect to the position of the aquifer's outlet (figure 1). In steady state, combining Darcy's law and the incompressibility condition yields the Laplace equation,

$$ \nabla^2 \phi = 0 \qquad (4.2) $$

where $\phi = p/\rho g + y$ is the hydraulic head.

To solve it, we must complement equation (4.2) with boundary conditions. We first assume that the flow does not vary across the experiment; the two plexiglas plates thus confine it to a vertical plane. The three walls bounding the aquifer are impervious, and the normal velocity of groundwater vanishes along them:

$$ \frac{\partial \phi}{\partial n} = 0, \qquad (4.3) $$

where $n$ denotes the direction of the normal to the wall.

We now turn our attention to the water table, $y = h(x)$. It is the surface where the pressure is atmospheric. Arbitrarily setting the latter to zero, this definition translates into

$$ \phi = y \qquad (4.4) $$
<figure><img src="image_6.png" /><figcaption>Figure 5: First three modes of the complex velocity potential (4.11) in an aquifer of aspect ratio a = 1. (a-c) Dimensionless stream function $K\psi/(RL)$ at the free surface ($y = 0$). (d-f) Streamlines (black dashed line) and hydraulic head (blue shading, normalized to 1).</figcaption></figure>

at the free surface ($y = h(x)$).

When rain falls, the water table deforms to accommodate the input of water. One must therefore deduce its shape from the flow equations: it is part of the solution. This free-boundary problem thus requires an additional boundary condition, which is but the water balance across the water table (Van de Giesen et al. 1994). In steady state, the latter reads:

$$ \frac{\partial \phi}{\partial y} - \frac{\partial \phi}{\partial x} \frac{\partial h}{\partial x} = \frac{R}{K} \qquad (4.5) $$

at the free surface ($y = h(x)$), where the rainfall rate $R$ and the hydraulic conductivity $K$ have units of velocity.

The Laplace equation (4.2), subjected to the boundary conditions (4.3), (4.4), and (4.5), describes how rainfall drives the flow of groundwater through the aquifer in steady state. In the absence of rainfall ($R = 0$), the water table does not deform ($h = 0$). In this trivial case, the pressure is hydrostatic ($\phi = 0$) and the groundwater remains still.

## 4.2. Weak rainfall

We could not find any closed-form solution of the system of equations (4.2), (4.3), (4.4), and (4.5). Fortunately, in our experiments, the rainfall rate is much smaller than the hydraulic conductivity, and the ratio of the two, $R/K$, never exceeds $4.1 \cdot 10^{-3}$ (table 1). The same observation often applies in nature, where the average recharge rate is typically much smaller than the conductivity of unconfined aquifers (Haitjema & Mitchell-Bruker 2005; Guérin et al. 2019). The ratio of these two quantities is therefore a small parameter, with respect to which we may linearize equation (4.5). When this ratio is vanishingly small ($R/K \ll 1$), the water table is virtually horizontal, and its deformation barely visible (1 cm at most in our experiment). Accordingly, we now linearize the boundary condition at the free surface (4.5) around the hydrostatic solution ($h = 0, \phi = 0$, Dagan, 1964):

$$ \frac{\partial \phi}{\partial y} = \frac{R}{K} \quad \text{in } y = 0. \qquad (4.6) $$
In the weak-rainfall regime, the flow obeys the Laplace equation (4.2) with boundary conditions (4.3) and (4.6). As anticipated, the elevation of the water table, h, drops out of the equations at this order: we no longer need to solve a Stefan problem. We solve this linear system using complex analysis (Polubarinova-Kochina 1962). An advantage of doing so is that it provides simultaneously the hydraulic head φ and the stream function ψ. Accordingly, we introduce the complex coordinate z = x + iy (with i² = -1), and look for an analytical function Φ(z) = φ(z) + iψ(z) that satisfies the boundary conditions. The stream function ψ, like the hydraulic head φ, are rescaled with the hydraulic conductivity K (they both have units of length).

Along the impervious walls (x = 0, x = -L, and y = -H), equation (4.3) translates into

$$ \psi = 0, \qquad (4.7) $$

where we have set the arbitrary constant of the stream function to zero. Using the Cauchy-Riemann relations, we rewrite equation (4.6) in terms of the stream function, and integrate the result along the x-axis:

$$ \psi = -\frac{RL}{K} \left(1 + \frac{x}{L}\right) \quad \text{in } y = 0. \qquad (4.8) $$

This condition, together with equation (4.7), means that the stream function is discontinuous at the outlet. This discontinuity is an inevitable consequence of the linearization of the free surface, which reduces the aquifer outlet to a point (section 4.3). As expected for the stream function, the amplitude [ψ] of this discontinuity is the discharge of the aquifer, which matches the rainfall input: KW[ψ] = Q = RLW.

The flow of groundwater in our aquifer thus obeys equation (4.2) with boundary conditions (4.7) and (4.8). To solve this system, we first note that any periodic function of the form

$$ \Phi_k(z) = \frac{RL}{K} \cosh \left[ k\pi \left( \frac{iz}{L} - a \right) \right], \qquad (4.9) $$

where k is an integer, satisfies boundary condition (4.7) in an aquifer of aspect ratio a = H/L (figure 5). Equation (4.9) is not the solution of our problem, however, as it does not fulfill the kinematic condition (4.8) at the water table. Fortunately, our problem is now linear, and we may look for its solution in the form of a series of periodic modes:

$$ \Phi(z) = \sum_{k=1}^{\infty} C_k \Phi_k(z). \qquad (4.10) $$

where C_k are real constants. We now need to find the C_k coefficients for which the series converges to boundary condition (4.8) at the surface.

Noting that, for each mode k, the stream function ψ_k reduces to a cosine at the surface, we use Fourier transform to compute the C_k coefficients. This is tantamount to finding the Fourier series of the saw-tooth function — a textbook problem (Farlow 1993). We thus find the expression for the complex velocity potential in the form of a Fourier series:

$$ \Phi(z) = -2 \frac{RL}{K} \sum_{k=1}^{\infty} \frac{\cosh(k\pi(iz/L - a))}{k\pi \sinh(k\pi a)}. \qquad (4.11) $$

The imaginary part of this velocity potential yields the stream function ψ(x, y), the contours of which are streamlines. To draw them, we need to (i) define a numerical grid over the aquifer's extent, (ii) evaluate the stream function on this grid, (iii) extract the contours of this function. Accordingly, we start by building a triangular mesh of the aquifer, using the internal mesher of FreeFem++ (Hecht, 2012, Python routines available online:
<figure><img src="image_7.png" /><figcaption>Figure 6: Steady-state flow in aquifers of aspect ratio *a* = 0.35, 0.7 (as in our laboratory experiment), and 1.4, respectively, after equation (4.11). (a-c) Dimensionless travel time of rainwater, τ*, as a function of relative landing position, *x*<sub>start</sub>/L. Solid blue line: numerical approximation (section 5.2). Black dashed line: expansion near the outlet, equation (5.2). (d-f) Streamlines (black dashed lines) and hydraulic head (blue shading).</figcaption></figure>

https://github.com/odevauchelle/pyFreeFem). Based on equation (4.11), we then evaluate the stream function on this mesh. Finally, we use the Matplotlib library to extract its contours, which are the streamlines (figure 6). The accuracy of this procedure depends on the mesh size $\delta l$, and on the truncation of the series (4.11). This truncation generates large, high-frequency oscillations near the singularity, an instance of the Gibbs phenomenon, which we attenuate with the Lanczos $\sigma$ factors (Hamming, 1973, p. 534). By trial and error, we opt for $\delta l/L = 1.25 \cdot 10^{-2}$ and 100 modes, a good trade-off between accuracy and computation speed.

The streamlines corresponding to each mode of the velocity potential (4.11), when plotted individually, divide the water table between regions where water enters the aquifer, and regions where water seeps out of it (figure 5); the incoming flux balances the outgoing one exactly. The summation of a large number of modes preserves this mass balance, but distributes the flux differently. Rainfall uniformly enters the aquifer through the water table, and all the streamlines converge towards a point sink, through which the entire discharge exits the aquifer (figure 6). The mass balance then requires that the flow velocity diverges at the outlet, thus creating a mathematical singularity. We shall characterize this singularity in the next section.

The streamlines that correspond to equation (4.11) are consistent with our experimental observations (figure 6). Oriented downward near the water table, they gradually curve as they plunge into the aquifer. Eventually, they rise towards the aquifer's outlet, where they
all converge. Streamlines are especially bent near the two stagnation points located in the bottom corners of the aquifer.

A more quantitative comparison with the experiments reveals that the streamlines drawn on figure 4 after equation (4.11) consistently reproduce the travel of the dye through our laboratory aquifer. This linearized theory even accounts for the thinning of the dye parcel near the bottom of the aquifer (figure 4). These observations support the use of the weak-rainfall approximation in our experiments.

Under this approximation, the intensity of rainfall affects only the norm of the groundwater velocity, leaving the geometry of the flow unchanged. The shape of the streamlines depends solely on the aspect ratio of the aquifer, $a = H/L$. Not surprisingly, we find that the streamlines approach the horizontal as we decrease the latter (figure 6). Yet, even for a very shallow aquifer, the vertical component of the flow remains predominant near the outlet and near the divide, in contradiction with the Dupuit-Boussinesq approximation (figure 6d).

As the water table bulges to accommodate the rainfall input, the volume of water stored in the aquifer increases. In steady-state, the intensity of rainfall thus controls the volume of water the aquifer stores above the outlet. One can easily measure this quantity by collecting all the water that seeps out of the experiment after the rain has stopped. It is therefore tempting to calculate this volume by integrating the elevation of the water table above the outlet along the aquifer's length:

$$V_e = W \int_{-L}^{0} h(x) dx. \qquad (4.12)$$

In the present theory, we use the boundary condition (4.4) to find the water table elevation at the free surface, $h(x) \approx \phi(x, y = 0)$. As one might expect, near the drainage divide, the water table lies above the outlet ($h > 0$). It decreases towards the outlet, however, and eventually plunges below it ($h < 0$), as the singularity associated to the exit condition takes over the flow (section 4.3). Although unrealistic, this mathematical behavior makes it possible for equation (4.12) to yield a vanishing or even negative volume. In fact, integrating this equation, we find that the positive and negative contributions of $h$ compensate for each other exactly (Appendix A). Regardless of the rainfall rate or the aquifer's aspect ratio, the volume stored during rainfall always vanishes. Thus, although the weak rainfall approximation faithfully represents the flow streamlines, it cannot account for the water stored in the aquifer during rainfall. The estimation of this volume requires to solve the original, non-linear Stefan problem.

Our experiments, like the linear theory, are limited to the weak-rainfall regime. At large rainfall rates, we expect that the deformation of the water table will eventually break the linear approximation, thus changing the flow geometry. However, this non-linear regime is beyond the scope of the present paper, and the capacity of our experimental setup.

### 4.3. Flow near the outlet

The series expansion of section 4.2 becomes singular near the outlet. To characterize this singularity mathematically, we could expand the modes of the series as $z$ approaches zero, but we would then find that nested series expansions obscure the issue. Instead, we propose here a longer route that, hopefully, offers a clearer view of the singularity.

The flow near the aquifer's outlet should be insensitive to the position of the bottom wall. Taking advantage of this, we consider an infinitely deep aquifer, by sending the bottom wall to infinity ($H \to \infty$). The problem then reduces to the Laplace equation (4.2) with boundary condition (4.6) at the water table, while boundary condition (4.7) applies only along the right and left walls. Using conformal mapping (Polubarinova-Kochina 1962), we find that
<figure><img src="image_8.png" /><figcaption>Figure 7: (a) Flow in an infinitely deep aquifer. (b) Asymptotic regime near the singularity at the aquifer outlet. Black dashed line: streamlines (linearly spaced), blue shading: hydraulic head. $\theta$: inclination of the streamlines with respect to horizontal.</figcaption></figure>

the velocity potential then reads (figure 7a):

$$ \Phi \underset{H \to +\infty}{\longrightarrow} \Phi_{\text{deep}} = \frac{2}{\pi} \frac{RL}{K} \left[ \log \left( 1 - e^{i\pi z/L} \right) - i\pi \left( \frac{z}{L} + 1 \right) \right]. \quad (4.13) $$

The corresponding streamlines resemble that of an aquifer of finite depth, except near the drainage divide where they go down to infinity (figure 7a).

Near the outlet, the flow is also insensitive to the left wall. Accordingly, sending the latter to infinity, or, equivalently, expending (4.13) about $z = 0$, yields the expression of the velocity potential near the outlet:

$$ \Phi \underset{z \to 0}{\longrightarrow} \Phi_{\text{outlet}} = \frac{RL}{K} \left[ \frac{2}{\pi} \log \left( \frac{z}{L} \right) + i \right]. \quad (4.14) $$

As expected for a point sink in a Laplacian field, the velocity potential diverges logarithmically, with a prefactor proportional to the ratio of the discharge per unit width, $Q/W = RL$, to the conductivity $K$ (figure 7b). The corresponding streamlines form spokes converging towards the outlet, in agreement with the trajectory of the dye parcels in our experiment (figure 4). This flow pattern is independent of the prefactor in equation (4.14). In addition, not too close to the exit grid, the flow is also insensitive to the details of the outlet's geometry. Transposed into a natural setting, this observation means that the river's shape does not affect the far field of the groundwater flow, although it may add to the pressure head a constant, which depends on the water level in the river.

The spoke pattern formed by the streamlines, and their even spacing, mean that the flux of water is uniformly distributed around the outlet. The same observation holds at the surface of the aquifer: since the rainfall is uniform, streamlines are evenly spaced there too. The water mass balance then imposes a linear relation between the inclination of a streamline near the outlet, $\theta$ (figure 7b), and its entry point into the aquifer, $x_{\text{start}}$ (figure 6):

$$ x_{\text{start}} = -\frac{2\theta}{\pi} L. \quad (4.15) $$

One sometimes uses a series of boreholes to measure how the chemical composition of
groundwater changes with depth (Haria & Shand 2004). Equation (4.15), if it holds in the field, would allow us to relate this depth to the path of groundwater through the aquifer.

Now that we have established a mathematical description of the streamlines, we turn our attention to the time it takes for a fluid parcel to travel along them. Most likely, the outlet singularity, which concerns only the end of the travel, will barely affect this travel time.

# 5. Travel time

## 5.1. Travel of a parcel of dye

On average, fluid parcels travel along the streamlines at a velocity equal to the flux of water divided by the aquifer porosity, $q/s$ (Guyon et al. 2001). Accordingly, the time that a particle needs to move between two points of complex coordinates $z_0$ and $z_1$ reads:

$$t = \frac{s}{K} \int_{z_0}^{z_1} \frac{dl}{|\nabla \phi|}, \qquad (5.1)$$

where the integral follows the streamline that connects $z_0$ and $z_1$ ($l$ is the distance along this streamline).

The travel time computed from equation (5.1) only takes into account the advection by the flow. In reality, the velocity of water varies significantly within a single pore, and between pores (Saffman 1959). The dispersion thus induced is visible in our experiment: the parcels of dye spread out along their path (figure 3). We shall discuss how this dispersion influences the distribution of travel times in section 6. First, however, we ignore its influence, and compute the travel time associated to advection only.

Using equation (5.1), we compute numerically the positions of the six parcels of dye injected in our experimental aquifer as a function of time. To do so, we extract the shape of the six streamlines that originate from the injection points (section 4.2). We then evaluate the hydraulic head along their course, compute the magnitude of its gradient, and evaluate integral (5.1) numerically. This procedure yields the travel time as a function of the position along the streamline. It accounts for the motion of the parcels, without any free parameter, but cannot represent their spreading (figure 3). Encouraged by this observation, we now investigate how this travel time depends on where a raindrop lands.

## 5.2. Travel times

Where a raindrop lands sets the origin, $x_{start}$, of the streamline along which it will travel underground. We may therefore calculate the travel time $\tau$ of a raindrop by integrating numerically equation (5.1) along this streamline, as in the previous section. In practice, however, the existence of three singular points — two stagnation points on the left and right corners of the aquifer bottom, and a singularity at the outlet — complicates this calculation. To deal with this problem, we build a mesh of the aquifer that skirts the two stagnation points and the aquifer outlet (figure 8). To do so, we trim three corners off the mesh along circles of radius $\epsilon$. Some of the streamlines cut through these circles. Equation (5.1) provides the contribution to the travel time of the segments located outside the corresponding disks. To complete the calculation, we need to work out the travel time through the disk that surrounds each singularity.

Replacing the velocity potential with its expression near the outlet (4.14) in equation (5.1), we compute the time water takes to cross the disk surrounding the outlet (figure 8):

$$\tau_{\text{outlet}} = \frac{\pi sL}{4R} \left(\frac{\epsilon}{L}\right)^2. \qquad (5.2)$$
<figure><img src="image_9.png" /><figcaption>Figure 8: Mesh (grey triangles) of an aquifer skirting the two stagnation points and the aquifer outlet (grey disks). Black dotted lines : streamlines (logarithmically spaced). Blue shading : hydraulic head. The aquifer aspect ratio is a = 0.5. For illustration, we set the radius of the quarter circles to $\epsilon/L = 0.1$ and the maximum mesh size to $\delta\ell/L = 0.025$. In practice, travel times are calculated with $\epsilon/L = 0.05$ and $\delta\ell/L = 0.0125$.</figcaption></figure>

The radial symmetry of the logarithmic singularity makes $\tau_{\text{outlet}}$ independent of the direction along which the streamline approaches the outlet.

We now turn our attention to the two stagnation points (figure 8). To compute the integral in equation (5.1), we first note that $|\nabla\phi| = |\Phi'|$ and $dl = |dz| = |d\Phi|/|\Phi'|$, where $'$ denotes the complex derivative. Along a streamline, the stream function is constant, while the hydraulic head decreases, so $|d\Phi| = -d\phi$. Combining all these expressions, we rewrite equation (5.1) as

$$ \tau_s = -\frac{s}{K} \int_{\phi_{\text{in}}}^{\phi_{\text{out}}} \frac{d\phi}{|\Phi'|^2}, \quad (5.3) $$

where $\phi_{\text{in}}$ and $\phi_{\text{out}}$ are the values of the hydraulic head at the entry and exit points of the disk, respectively. Near a stagnation point, the velocity potential is dominated by the first mode of equation (4.11). Taking advantage of this, we approximate the velocity potential by this mode only. Expanding the latter near the stagnation point, we obtain the following expression for the velocity potential:

$$ \Phi_s = \phi_c \left[ 1 - \frac{\pi^2}{2} \left( \frac{z-z_s}{L} \right)^2 \right]. \quad (5.4) $$

where $z_s$ is the position of the stagnation point ($z_s = -L - iH$ in the lower-left corner and $z_s = -iH$ in the lower-right one). The hydraulic head at the stagnation point, $\phi_c$, is that of the far field, namely $\phi_c = \mathcal{R}(\Phi(z_s))$.

Writing the derivative of the potential (5.4) as a function of the potential $\Phi_s$ itself, namely

$$ \Phi'_s = -\sqrt{2\pi} \frac{\phi_c}{L} \sqrt{1 - \frac{\Phi_s}{\phi_c}}, \quad (5.5) $$

allows us to turn the curvilinear integral (5.3) into a classical one. Integrating the latter over $\phi$, we find that the time spent by a fluid parcel inside the disk surrounding a stagnation point reads:

$$ \tau_s = \frac{sL}{R} \frac{\sinh(\pi a)}{2\pi} \left| \operatorname{asinh} \left( \frac{\phi_{\text{in}} - \phi_c}{\psi_0} \right) \right|, \quad (5.6) $$

where $\psi_0$ defines the streamline along which the water parcel travels.
<figure><img src="image_10.png" /><figcaption>Figure 9: Dimensionless travel time of a parcel of dye, $\tau^*$, as a function of the distance to the outlet $-x_{start}/L$ (lower axis). Upper axis: streamline inclination at the outlet, $\theta$. Blue markers: experimental runs (table 1). Vertical error bars indicate the exit duration, $\Delta\tau/\langle\tau\rangle$; horizontal ones indicate initial width of the parcel. Colored lines: theoretical travel times after section 5.2 for aspect ratios $a = 0.35$ (green), $0.7$ (blue, our laboratory experiment), and $1.4$ (orange).</figcaption></figure>

We can now evaluate the travel time of water through the aquifer by adding to integral (5.1) the contributions of the outlet (5.2) and, if need be, that of the two stagnation points (5.6). The accuracy of this calculation depends on the radius of the quarter-circles, $\epsilon$, and on the maximum mesh size $\delta\ell$. By trial and error, we choose $\delta\ell/L = 1.25 \cdot 10^{-2}$ and $\epsilon = 0.05$.

The mean travel time of water in the aquifer is the ratio of the volume of water it contains to its discharge: $\langle\tau\rangle = saL/R$ (Zwietering 1959; Nauman & Buffham 1983). This expression, which resembles a simple definition, is actually an exact consequence of the water-mass balance, combined with Green's theorem (appendix B). In the following, we rescale all travel times with this average, thus introducing the dimensionless travel time $\tau^* = \tau/\langle\tau\rangle$.

Figure 6 a-c shows the dimensionless travel time $\tau^*$ of a raindrop as a function of its landing point, $x_{start}$, for several values of the aquifer's aspect ratio. We find that $\tau^*$ increases continuously with the distance $|x_{start}|$ that separates the landing point from the aquifer outlet. A careful examination of its variation reveals two singularities: the travel time vanishes near the outlet, and diverges near the divide. We shall derive the asymptotic expressions that account for these regimes in section 5.3. Before that, however, we compare our theory to the travel time of dye in our experiment.

To estimate this travel time, we use the images that record the propagation of the dye (section 3), and decompose each of them into their hue, saturation, and value in the range [0, 255]. We select pixels with a saturation above 40 to ensure that their hue is well-defined, and retain the ones that fall into the hue range that corresponds to the color of the dye ([80, 130]). This procedure yields a series of well-contrasted pictures, where the parcels of dye are highlighted. Using these images, we measure the time at which the front of each parcel reaches the aquifer outlet, and the time at which its rear end finally leaves it. We then estimate the average travel time of the parcel as the mean of these two times, $\tau$. Their difference corresponds to the exit duration, $\Delta\tau$, which we use as an estimate of the
measurement dispersion. In fact, the latter results from pore-scale dispersion, a phenomenon we shall address in more details in section 6. Figure 9 compares our experiments to the theory. Although the latter only accounts for the advection by the flow, we find that it matches our observations. Encouraged by this result, we now examine the two asymptotic regimes that represent the behavior of the travel time near the singularity, and near the outlet.

### 5.3. Asymptotic regimes

When a rain drop lands near the outlet, its travel is dominated by the singularity that accelerates the flow towards the exit point. According to equation (5.2), the duration of this travel decreases as the square of the distance that separates its landing point from the outlet:

$$ \tau^* = \frac{\pi}{4a} \left( \frac{x_{\text{start}}}{L} \right)^2. \qquad (5.7) $$

Figure 10 a shows that our experiment conforms to this asymptotic behavior.

To the contrary, the travel time diverges near the left-hand wall of the aquifer, which simulates a drainage divide (figure 6a-c). To account for this behavior, we note that a streamline originating near the divide remains close to the walls along its entire course (figure 6d-f). In this region, the velocity potential is dominated by the first mode of equation (4.11). Taking again advantage of this observation, we approximate the velocity potential by its first mode, $\Phi_1$, and use equation (5.3) to calculate the corresponding travel time. Following a procedure similar to the one applied to the stagnation points in section 5.2, we find the following expression for the travel time near the drainage divide (Appendix C):

$$ \tau^* \approx -\frac{\sinh(\pi a)}{\pi a} \left[ \log \left( 1 + \frac{x_{\text{start}}}{L} \right) - \log \left( \frac{4}{1 + \cosh(\pi a)} \right) \right], \qquad (5.8) $$

which is a good approximation of the actual travel time (figure 10b). As expected, equation (5.8) features, at leading order, the logarithmic divergence that is the signature of the stagnation points. It also bears the mark of the aspect ratio of the aquifer, which affects both the prefactor of the logarithm, and the additive constant.

### 5.4. Distribution of travel times

In our experiment, rain falls uniformly onto the aquifer's surface. We may thus write the probability density function (PDF) of travel times as:

$$ f(\tau^*) = -\frac{1}{L} \frac{dx_{\text{start}}}{d\tau^*}. \qquad (5.9) $$

This expression is but a recasting of the function that relates the travel time of a drop of water to its landing spot.

We approximate equation (5.9) numerically, based on the results of section 5.2, and plot it for different values of the aspect ratio (figure 10c). The resulting distributions diverge at the origin. From there, they gradually decrease at a rate which depends on the aspect ratio of the aquifer.

Short travel times correspond to raindrops entering the aquifer near the outlet. Accordingly, the expression of the travel time near the outlet (5.7), combined with equation (5.9), yields the short-time behavior of the PDF:

$$ f(\tau^*) \sim \left(\frac{a}{\pi}\right)^{1/2} \tau^{*-1/2}. \qquad (5.10) $$

This expansion consistently represents the short time behavior of the numerical PDF (figure
<figure><img src="image_11.png" /><figcaption>Figure 10: Travel times through a rectangle aquifer. (a) and (b) Dimensionless travel time, $\tau^*$, as a function of distance to outlet, $-x_{start}/L$, and distance to divide, $(L + x_{start})/L$. Colored lines: theoretical travel times, after section 5.2. Black dashed lines: asymptotic regimes (5.7) and (5.8). Blue markers: experimental data (same symbols as in figure 9), vertical error bars indicate the exit duration, $\Delta\tau/\langle\tau\rangle$, horizontal ones indicate the initial width of the parcel. (c) and (d) Probability density functions (PDF) of travel times, $f(\tau^*)$, for several aspect ratios. Black dashed line: asymptotic regimes (5.10) and (5.11).</figcaption></figure>

10c). The divergence of the latter at short times is a consequence of the acceleration of the flow near the outlet.

Long travels, on the contrary, correspond to raindrops entering the aquifer near the drainage divide. Streamlines originating from this area entrain their water near the two stagnation points, where it lingers for a logarithmically long time. Accordingly, we now combine the travel time near the drainage divide (5.8) with equation (5.9). We find that the tail of each PDF decreases exponentially, according to

$$f(\tau^*) \sim \frac{4\pi a}{\sinh(\pi a)(1 + \cosh(\pi a))} e^{-\tau^*/\tau_c^*}, \quad (5.11)$$
<figure><img src="image_12.png" /><figcaption>Figure 11: Relative dispersion, $\Delta\tau/\tau$, as a function of the distance to outlet rescaled with pore size $-x_{start}/d_p$ (lower axis). Upper axis: streamline inclination near the outlet, $\theta$. Blue markers: experimental data (table 1). Black dashed line: equation (6.1) fitted to data.</figcaption></figure>

with a characteristic dimensionless time

$$ \tau_c^* = \frac{\sinh(\pi a)}{\pi a} \qquad (5.12) $$

Again, this expansion accords with the numerical PDF (figure 10d). As anticipated, the characteristic time $\tau_c^*$ increases with the aquifer's aspect ratio: deep aquifers induce long travels. In particular, $\tau_c^*$ tends to infinity for an infinitely deep aquifer, whereas it tends towards 1 for a shallow one ($a = 0$). In short, the tail of the travel time distribution depends primarily on the aspect ratio of the aquifer.

Although it only accounts for the advection by the flow, the above theory proves a fair representation of our experiment: this approximation suffices to reproduce the travel of individual parcels of dye (figure 3 and 9). In reality, however, the pore-scale hydrodynamic dispersion affects the travel time of groundwater. In the next section, we evaluate the magnitude of its contribution.

## 6. Hydrodynamic dispersion

Hydrodynamic dispersion gradually spreads a parcel of dye as it travels along its streamline (figure 3). In terms of travel time, this spreading becomes a dispersion time, $\Delta\tau$, defined as the time separating the moment when the front of the parcel reaches the aquifer's outlet, from the moment it has left it entirely. In section 5.2, we estimated this dispersion time, together with the average travel time, for each parcel of dye in our experiment. We now compare these two quantities, and find that the ratio of the dispersion time to the average travel time varies in the range $\Delta\tau/\tau \approx 0.4 - 1.3$: hydrodynamic dispersion indeed affects the propagation of dye through the aquifer (figure 11).

Our measurements also reveal that the ratio of the dispersion time to the travel time, $\Delta\tau/\tau$, decreases with $x_{start}$ (figure 11). The relative influence of hydrodynamic dispersion therefore decreases with the length of the travel. We now interpret this observation with scaling arguments. The time necessary for a parcel of dye to cross the aquifer's outlet scales
like $\Delta\tau \sim \ell_p/U$, where $\ell_p$ is the length of the parcel when it reaches the outlet, and $U$ is the characteristic velocity at which the parcel travels along the streamline. Assuming that the length of a streamline is proportional to the position of its origin, $-x_{start}$, this velocity scales like $U \sim -x_{start}/\tau$. As the pore-scale dispersion spreads the parcel, the length of the latter scales like $\ell_p \sim \sqrt{D\tau}$, where $D$ is the dispersion coefficient. Neglecting molecular diffusion, this coefficient is the product of the fluid velocity $U$ and a dispersion length $\lambda$ ($D = \lambda U$, Hillel, 2003). In a homogeneous aquifer, this dispersion length is proportional to the pore size, $d_p$. Combining these four equations, we find that $\Delta\tau/\tau$ decreases with the distance $x_{start}$ according to:

$$ \frac{\Delta\tau}{\tau} \approx \left( \frac{\lambda}{-x_{start}} \right)^{1/2}, \qquad (6.1) $$

A fit of equation (6.1) to our data consistently reproduces the evolution of the ratio $\Delta\tau/\tau$ in our experiment, with $\lambda \approx 1288$ $d_p \approx 12.9$ cm, a value close to that measured in sand columns (De Marsily 1986).

Equation (6.1) indicates that the relative importance of hydrodynamic dispersion at the pore scale with respect to advection at the aquifer scale is essentially governed by the ratio of the dispersion length, $\lambda$, to the distance, $x_{start}$, that separates the landing point of the raindrop from the aquifer's outlet.

## 7. Conclusion

The experiments presented in this paper reveal the flow pattern in a deep, unconfined aquifer recharged by rainfall. Since the outlet stands high above the aquifer's bottom, groundwater needs to go up to reach it, which breaks down the Dupuit-Boussinesq approximation. We find an analytical expression of the velocity potential, in the form of an infinite series of modes, which accounts for the shape of the experimental streamlines. This linear theory represents correctly the travel of water through the aquifer, provided that the intensity of rainfall is much less than the aquifer's conductivity, a condition often met in nature (Haitjema & Mitchell-Bruker 2005).

The underground fate of a raindrop depends on where it enters the aquifer. A raindrop that lands near the outlet flows out almost immediately. On the contrary, the trajectory of a drop fallen near the drainage divide reaches down to the bottom of the aquifer, before it rises up towards the outlet. The duration of this long travel is further increased by the stagnation points that drainage divides generate. These stagnation points control the tail of the travel-time distribution, which we find to decrease exponentially, with a characteristic time that depends on the aquifer's aspect ratio.

Our theory bears similarity to the model of Toth (1963). Like us, Toth restricts his analysis to small deformations of the water table, and linearizes the boundary condition at the free surface accordingly. However, whereas our theory explicitly accounts for the rainfall rate that recharges the aquifer, Toth (1963) imposes the elevation of the water table. Assuming that the latter inherits its shape from the land above it, he represents it with an uniform slope modulated by a sinusoidal wave (Toth, 1963, equation 5d). The sine mimics the hills and depressions of the land surface, while the slope represents the regional gradient. If the latter vanishes, the solution found by Toth (his equation 6) coincides with the first mode of equation (4.11) of the present paper. As this mode dominates the flow far below the water table, both models yield the same trajectory for a droplet fallen near the drainage divide. Toth's model can therefore account for the tail of the travel time distribution. Closer to the surface, however, Toth's model departs from the present theory, and, therefore, cannot represent the faster paths we observe in our experiment.
Our estimate of the travel time only takes into account advection by the flow. In reality, hydrodynamic dispersion affects the distribution of travel times. Our experiments indicate that the magnitude of this effect depends on the ratio of a dispersion length, proportional to the pore size, to the length of the streamline. In the field, these lengths most likely differ by several orders of magnitude, and we thus speculate that pore-scale dispersion is negligible.

Near the aquifer's outlet, the streamlines flow radially, with an inclination that depends on the location of their origin. This result, if it holds in nature, provides a simple method to reconstruct the distribution of travel times through an aquifer. Indeed, a series of boreholes, drilled at different depths near the edge of a river, would collect groundwater traveling along different streamlines. Sampling water in each borehole and measuring the time it spent underground, based on isotopic ratios or conservative tracers, could thus yield the distribution of travel times through the aquifer. One might then infer the aquifer's effective depth from the characteristic time of the distribution's tail.

Our experiments are based on the measurement of the position and size of dye parcels travelling through the aquifer. Such measurements are hardly possible in the field. Instead, one usually injects a chemical tracer into a well, and records the evolution of its concentration in another downstream one. The breakthrough curve, i.e. the time series of the concentration in the second well, is then used to infer the flow within the aquifer (Gelhar et al. 1992). In our experiment, the absorption of light by the dye might, in principle, allow us to follow its distribution in the aquifer. This method would allow us to produce breakthrough curves, and to compare them with field measurements.

Finally, our experiments reveal that the reaction of the water table to rainfall is considerably faster in a deep aquifer than in a shallow one. This observation calls for an investigation of the stormflow regime, which is crucial to flood forecasting (Guérin et al. 2019). Our experiments suggest that this regime mobilizes the entire volume of the aquifer instantly, but a theory of this dynamical regime remains to be established.

We are grateful to S. Deboeuf, F. James, P. Gondret, Y. Meheust, F. Métivier, J. Neufeld, C. Soulaine, and P. Szymczak for their insightful comments. We thank M. Vulliet and M. Colledge, who performed some of the experiments as part of a tutored project. We also express our gratitude to O. Bour for his guidance through the travel-time literature. Finally, we thank an anonymous referee for suggesting the use of the Lanczos sigma factors.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCID.** V. Jules, https://orcid.org/0000-0003-1594-8569; E. Lajeunesse, https://orcid.org/0000-0002-0950-6054; O. Devauchelle, https://orcid.org/0000-0002-7295-4896; A. Guérin, https://orcid.org/0000-0001-6197-4587; C. Jaupart, https://orcid.org/0000-0001-7961-0538; P-Y. Lagrée, https://orcid.org/0000-0002-3931-6622.

## Appendix A. Volume of groundwater stored in the aquifer

During rainfall, the water table bulges, and the volume of water stored in the aquifer increases. In steady-state, the volume of water that the aquifer stores above the outlet reads :

$$V_e = W \int_{-L}^{0} h(x) dx, \qquad (A 1)$$

where $h(x)$ is the elevation of the water table with respect to the outlet. At leading order, it is the hydraulic head at $y = 0$: $h(x) = \phi(x, y = 0)$. The stored volume is therefore:

$$V_e = W \int_{-L}^{0} \phi(x, y = 0) dx. \qquad (A 2)$$
To calculate this integral, we first rewrite it in terms of the complex potential,

$$V_e = W\mathfrak{R}\left\{\int_{-L}^{0} \Phi(z)dz\right\}. \qquad (A3)$$

Taking advantage of the Cauchy integral theorem, we evaluate this integral along a contour that follows the aquifer's walls,

$$V_e = \mathfrak{R} \left\{ \int_{-L}^{-L(1+ia)} \Phi(z)dz + \int_{-L(1+ia)}^{-Lia} \Phi(z)dz + \int_{-Lia}^{0} \Phi(z)dz \right\}. \quad (A4)$$

As only the second integral has a finite real part, the stored volume simplifies into:

$$V_e = \mathfrak{R} \left\{ \int_{-L(1+ia)}^{-Lia} \Phi(z) dz \right\}. \qquad (A5)$$

Using equation (4.11), we rewrite this integral as:

$$V_e = -2 \frac{RL}{K} \sum_{k=1}^{\infty} \frac{1}{k\pi \sinh(k\pi a)} \int_{-L}^{0} \cos\left(\frac{k\pi x}{L}\right) dx, \qquad (A6)$$

which yields:

$$V_e = 0. \qquad (A7)$$

In the weak rainfall approximation, the volume of water stored in the aquifer during rainfall is zero. This, of course, is only a first-order result. To find how much water a deep aquifer can accommodate in steady state, we need to return to the original, non-linear Stefan problem.

## Appendix B. Mean travel time of water in the aquifer

As the flow is in steady state, we can attribute an age $a$ to each point in our aquifer, that corresponds to the time it takes for water to travel from the surface to this point. The mean travel time of groundwater is the flux of age through the contour $\Gamma$ that bounds the aquifer, normalized by the outlet discharge, $Q = RLW$:

$$\langle \tau \rangle = \frac{W}{Q} \oint_{\Gamma} a \mathbf{q} \cdot \mathbf{n} \, dl, \qquad (B1)$$

where $l$ is the distance along the contour of the aquifer. Using Gauss's theorem, we rewrite equation (B1) as:

$$\langle \tau \rangle = \frac{W}{Q} \left( \iint_{\Sigma} \mathbf{q} \cdot \nabla a \, d\Sigma - \iint_{\Sigma} a \nabla \cdot \mathbf{q} \, d\Sigma \right), \qquad (B2)$$

where $\Sigma$ is the flow domain, that is, a vertical section of the aquifer of area $aL^2$ (figure 1).

As a parcel of groundwater travels through the aquifer, its age $a$ increases with time according to:

$$\frac{Da}{Dt} = \frac{\partial a}{\partial t} + \mathbf{v} \cdot \nabla a = 1, \qquad (B3)$$

where $\mathbf{v} = \mathbf{q}/s$ is the average velocity of the parcel (Davis & Bentley 1982; Varni & Carrera 1998). In steady state, the above equation becomes:

$$\mathbf{q} \cdot \nabla a = s. \qquad (B4)$$

Combining equations (B2), (B4), and the incompressibility condition $(\nabla \cdot \mathbf{q}) = 0$, we
finally get the expression of the mean travel time :

$$ \langle \tau \rangle = \frac{saL}{R}. \qquad (B 5) $$

The mean travel time is therefore equal to the volume of water in the aquifer divided by its discharge (Zwietering 1959; Nauman & Buffham 1983). It is independent of the aquifer's shape.

## Appendix C. Expansion of the travel time near the drainage divide

A streamline originating near the divide remains close to the walls, where the velocity potential is dominated by the first mode of equation (4.11) (figure 6d-f). We therefore approximate the velocity potential by $\Phi_1$, and use equation (5.3) to express the dimensionless travel time :

$$ \tau^* = -\frac{2R}{aLK} \int_{\phi_{\text{in}}}^{\phi_{\text{out}}} \frac{d\phi_1}{|\Phi'_1|^2}, \qquad (C 1) $$

where we integrate along the streamline. Here, we invoked the symmetry of the first mode, and integrated between the origin of the streamline, $z_{\text{in}} = -L + \delta$, and its mid-course point, $z_{\text{out}} = -L/2 + i(\delta - H)$, where $\delta \ll L$ is the distance separating the origin of the streamline from the divide.

Expressing $\Phi'_1$ as a function of $\Phi_1$, we turn the above curvilinear integral into a classical one,

$$ \tau^* = -\frac{2 \sinh(\pi a)}{\pi a} \int_{\phi_{r, \text{in}}}^{\phi_{r, \text{out}}} \frac{d\phi_r}{|\Phi_r| |2 + \Phi_r|}. \qquad (C 2) $$

where we introduced the rescaled potential $\Phi_r = \phi_r + i\psi_r = \cosh(\pi(iz/L - a)) - 1$. The bounds of this integral are the values of $\phi_r$ in $z_{\text{in}}$ and $z_{\text{out}}$.

Near the stagnation point, both $\phi_r$ and $\psi_r$ at small. At leading order, the integrand $I(\phi_r, \psi_r)$ of equation (C 2) becomes:

$$ I \approx I_0 = \frac{1}{2|\Phi_r|} = \frac{1}{2\sqrt{\phi_r^2 + \psi_r^2}}. \qquad (C 3) $$

As expected, this expression is not integrable when $\psi_r$ vanishes: the travel time diverges as the flow line approaches the stagnation point. Using the above expression only, the integral of equation (C 2) has a closed-form expression, namely

$$ \tau_0^* = -\frac{\sinh(\pi a)}{\pi a} \left( \operatorname{asinh} \left( \frac{\phi_{r, \text{in}}}{\psi_r} \right) - \operatorname{asinh} \left( \frac{\phi_{r, \text{out}}}{\psi_r} \right) \right), \qquad (C 4) $$

which diverges logarithmically as $\psi_r$ vanishes, or equivalently, as the landing point of a raindrop approaches the divide:

$$ \tau_0^* \approx -\frac{\sinh(\pi a)}{\pi a} \log \left( 1 + \frac{x_{\text{start}}}{L} \right). \qquad (C 5) $$

The logarithmic divergence (C 5) accounts for the behavior of the travel time near the divide, but only within an additive constant, $\tau_1^*$. To calculate this constant, we subtract the singularity from the leading order integrand:

$$ I - I_0 = \frac{2 - |2 + \Phi_r|}{2|\Phi_r||2 + \Phi_r|}. \qquad (C 6) $$
By design, this quantity does not diverge as $\psi_r$ vanishes, and we thus approximate it with

$$I - I_0 \approx -\frac{\operatorname{sign}(\phi_r)}{2(2 + \phi_r)}. \qquad (C 7)$$

Integrating the latter between $\phi_{r,in} = \cosh(\pi a) - 1$ and $\phi_{r,out} = -1$ (the former is positive, and the latter negative) yields the constant we are looking for:

$$\tau_1^* = \frac{\sinh(\pi a)}{\pi a} \log \left( \frac{4}{1 + \cosh(\pi a)} \right). \qquad (C 8)$$

Adding equations (C 5) and (C 8) yields the expansion of the travel time near the drainage divide, equation (5.8).

## REFERENCES

* ANDREOTTI, BRUNO, FORTERRE, YOËL & POULIQUEN, OLIVIER 2013 *Granular media: between fluid and solid*. Cambridge University Press.

* BOUSSINESQ, JOSEPH 1903 Sur un mode simple d'écoulement des nappes d'eau d'infiltration à lit horizontal, avec rebord vertical tout autour lorsqu'une partie de ce rebord est enlevée depuis la surface jusqu'au fond. *CR Acad. Sci* **137** (5), 11.

* BRESCIANI, ETIENNE, DAVY, PHILIPPE & DE DREUZY, J-R 2014 Is the dupuit assumption suitable for predicting the groundwater seepage area in hillslopes? *Water Resources Research* **50** (3), 2394–2406.

* BRUTSAERT, W. & NIEBER, J.L. 1977 Regionalized drought flow hydrographs from a mature glaciated plateau. *Water Resources Research* **13** (3), 637–643.

* CARDENAS, M BAYANI 2007 Potential contribution of topography-driven regional groundwater flow to fractal stream chemistry: Residence time distribution analysis of tóth flow. *Geophysical Research Letters* **34** (5).

* CARTWRIGHT, IAN, CENDÓN, DIONI, CURRELL, MATTHEW & MEREDITH, KARINA 2017 A review of radioactive isotopes and other residence time tracers in understanding groundwater recharge: Possibilities, challenges, and limitations. *Journal of Hydrology* **555**, 797–811.

* CHARLAIX, E, HULIN, JP & PLONA, TJ 1987 Experimental study of tracer dispersion in sintered glass porous materials of variable compaction. *The Physics of fluids* **30** (6), 1690–1698.

* DAGAN, G 1964 Second order linearized theory of free-surface flow in porous media. *La Houille Blanche* (8), 901–910.

* DARCY, H. 1856 *Les fontaines publiques de la ville de Dijon*. Dalmont.

* DAVIS, STANLEY N & BENTLEY, HAROLD W 1982 Dating groundwater: A short review. ACS Publications.

* DE MARSILY, GHISLAIN 1986 Quantitative hydrogeology. *Tech. Rep.* Paris School of Mines, Fontainebleau.

* DENTZ, MARCO, LE BORGNE, TANGUY, ENGLERT, ANDREAS & BIJELJIC, BRANKO 2011 Mixing, spreading and reaction in heterogeneous media: A brief review. *Journal of contaminant hydrology* **120**, 1–17.

* DEVAUCHELLE, OLIVIER, PETROFF, ALEXANDER P, SEYBOLD, HANSJÖRG F & ROTHMAN, DANIEL H 2012 Ramification of stream networks. *Proceedings of the National Academy of Sciences* **109** (51), 20832–20836.

* DUPUIT, J. 1863 *Études théoriques et pratiques sur le mouvement des eaux dans les canaux découverts et à travers les terrains perméables*. Dunod.

* FARLOW, STANLEY J 1993 *Partial differential equations for scientists and engineers*. Courier Corporation.

* GELHAR, LYNN W, WELTY, CLAIRE & REHFELDT, KENNETH R 1992 A critical review of data on field-scale dispersion in aquifers. *Water resources research* **28** (7), 1955–1974.

* VAN DE GIESEN, NC, PARLANGE, J-Y & STEENHUIS, TS 1994 Transient flow to open drains: Comparison of linearized solutions with and without the dupuit assumption. *Water Resources Research* **30** (11), 3033–3039.

* GLEESON, TOM, BEFUS, KEVIN M, JASECHKO, SCOTT, LUIJENDIJK, ELCO & CARDENAS, M BAYANI 2016 The global volume and distribution of modern groundwater. *Nature Geoscience* **9** (2), 161–167.

* GODSEY, SARAH E, AAS, WENCHE, CLAIR, THOMAS A, DE WIT, HELEEN A, FERNANDEZ, IVAN J, KAHL, J STEVE, MALCOLM, IAIN A, NEAL, COLIN, NEAL, MARGARET, NELSON, SARAH J & OTHERS 2010 Generality of fractal 1/f scaling in catchment tracer time series, and its implications for catchment travel time distributions. *Hydrological Processes* **24** (12), 1660–1671.
GUÉRIN, ADRIEN 2015 Dynamics of groundwater flow in an unconfined aquifer. Theses, Université Paris Diderot (Paris 7), Sorbonne Paris Cité.

GUÉRIN, ADRIEN, DEVAUCHELLE, OLIVIER & LAJEUNESSE, ERIC 2014 Response of a laboratory aquifer to rainfall. *Journal of Fluid Mechanics* **759**, R1.

GUÉRIN, ADRIEN, DEVAUCHELLE, OLIVIER, ROBERT, VINCENT, KITOU, THIERRY, DESSERT, CÉLINE, QUIQUEREZ, AMÉLIE, ALLEMAND, PASCAL & LAJEUNESSE, ÉRIC 2019 Stream-discharge surges generated by groundwater flow. *Geophysical Research Letters* **46** (13), 7447–7455.

GUYON, E., HULIN, J.P. & PETIT, L. 2001 *Hydrodynamique physique*. EDP Sciences/CNRS Editions.

HAITJEMA, HENK M & MITCHELL-BRUKER, SHERRY 2005 Are water tables a subdued replica of the topography? *Groundwater* **43** (6), 781–786.

HAMMING, RW 1973 Numerical methods for scientists and engineers.

HARIA, ATUL H & SHAND, PAUL 2004 Evidence for deep sub-surface flow routing in forested upland wales: implications for contaminant transport and stream flow generation. *Hydrology and Earth System Sciences Discussions* **8** (3), 334–344.

HECHT, F. 2012 New development in freefem++. *J. Numer. Math.* **20** (3-4), 251–265.

HILLEL, DANIEL 2003 *Introduction to environmental soil physics*. Elsevier.

INGEBRITSEN, STEVEN E & SANFORD, WARD E 1999 *Groundwater in geologic processes*. Cambridge University Press.

KIRCHNER, J.W., FENG, X., NEAL, C. & OTHERS 2000 Fractal stream chemistry and its implications for contaminant transport in catchments. *Nature* **403** (6769), 524–527.

LE BORGNE, TANGUY, DE DREUZY, J-R, DAVY, PHILIPPE & BOUR, OLIVIER 2007 Characterization of the velocity field organization in heterogeneous media by conditional correlation. *Water resources research* **43** (2).

LEHR, JAY H 1963 Groundwater: Flow toward an effluent stream. *Science* **140** (3573), 1318–1320.

LOBKOVSKY, ALEXANDER E, JENSEN, BILL, KUDROLLI, ARSHAD & ROTHMAN, DANIEL H 2004 Threshold phenomena in erosion driven by subsurface flow. *Journal of Geophysical Research: Earth Surface* **109** (F4).

MAHER, K 2010 The dependence of chemical weathering rates on fluid residence time. *Earth and Planetary Science Letters* **294** (1-2), 101–110.

NAUMAN, E BRUCE & BUFFHAM, BA 1983 *Mixing in continuous flow systems*. John Wiley & Sons Inc.

PETROFF, AP, DEVAUCHELLE, O, ABRAMS, DM, LOBKOVSKY, AE, KUDROLLI, A & ROTHMAN, DH 2011 Geometry of valley growth. *Journal of Fluid Mechanics* **673** (245-254), 6.

POLUBARINOVA-KOCHINA, P.YA. 1962 Theory of ground water movement princeton university press. Princeton, NJ.

POWERS, WL, KIRKHAM, DON & SNOWDEN, G 1967 Orthonormal function tables and the seepage of steady rain through soil bedding. *Journal of Geophysical Research* **72** (24), 6225–6237.

READ, WW 1993 Series solutions for laplace's equation with nonhomogeneous mixed boundary conditions and irregular boundaries. *Mathematical and Computer Modelling* **17** (12), 9–19.

REMPE, DANIELLA M & DIETRICH, WILLIAM E 2014 A bottom-up control on fresh-bedrock topography under landscapes. *Proceedings of the National Academy of Sciences* **111** (18), 6576–6581.

RUBINSTEIN, L.I. 2000 *The stefan problem*, vol. 8.

SAFFMAN, PG 1959 A theory of dispersion in a porous medium. *Journal of Fluid Mechanics* **6** (3), 321–349.

SOUZY, MATHIEU, LHIASSIER, HENRI, MÉHEUST, YVES, LE BORGNE, TANGUY & METZGER, BLOEN 2020 Velocity distributions, dispersion and stretching in three-dimensional porous media. *Journal of Fluid Mechanics* **891**, A16.

TOTH, JOZSEF 1963 A theoretical analysis of groundwater flow in small drainage basins. *Journal of geophysical research* **68** (16), 4795–4812.

TROCH, PETER A, BERNE, ALEXIS, BOGAART, PATRICK, HARMAN, CIARAN, HILBERTS, ARNO GJ, LYON, STEVE W, PANICONI, CLAUDIO, PAUWELS, VALENTIJN, RUPP, DAVID E, SELKER, JOHN S & OTHERS 2013 The importance of hydraulic groundwater theory in catchment hydrology: The legacy of wilfried brutsaert and jean-yves parlange. *Water Resources Research* **49** (9), 5099–5116.

VARNI, MARCELO & CARRERA, JESÚS 1998 Simulation of groundwater age distributions. *Water Resources Research* **34** (12), 3271–3281.

ZWIETERING, THN 1959 The degree of mixing in continuous flow systems. *Chemical Engineering Science* **11** (1), 1–15.