2.1. Set-up

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. (Reference Guérin, Devauchelle and Lajeunesse2014), we reproduce this configuration in a narrow tank formed by two vertical Plexiglas® plates ($90\ \mathrm {cm} \times 92\ \mathrm {cm}$), separated by a gap of width $W = 5\ \mathrm {cm}$ (figure 1). This tank is filled with glass beads of size either $1$ or $3\ \textrm {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\ \mathrm {mm}$, length $50\ \mathrm {cm}$) uniformly distribute water over the aquifer's surface, with a standard deviation of approximately $3\,\%$.

Figure 1.

Experimental set-up 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\ \mathrm {cm}$. Its length and depth are $L=90$ cm and $H=62\ \mathrm {cm}$, respectively. A series of $24$ capillary tubes (diameter $1.2\ \mathrm {mm}$, length $50\ \mathrm {cm}$) placed above the tank simulate rainfall, by distributing water uniformly over the aquifer's surface.

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\ \mathrm {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\ \mathrm {cm}$ is the length of the aquifer. The thinness of the set-up ($W/L\approx 0.05$), and its cross-wise 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.3. Transient

Guérin et al. (Reference Guérin, Devauchelle and Lajeunesse2014) used an experiment similar to the present one, but for the position of the outlet which, in their set-up, 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. (Reference Guérin, Devauchelle and Lajeunesse2014), we first use the same glass beads (diameter $d_s=1\ \textrm {mm}$) and rainfall rate ($Q_{in} \approx 8\ \textrm {mL}\ \textrm {s}^{-1}$). Following the procedure of the previous section, we then record the time series of the discharge, often called hydrograph in this context. Less than five 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. (Reference Guérin, Devauchelle and Lajeunesse2014), 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 millimetres) than in the shallow one (10 to 20 centimetres).

Figure 2.

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

The difference between our observations and those of Guérin et al. (Reference Guérin, Devauchelle and Lajeunesse2014) 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\ \textrm {mm}$ glass beads for the $1\ \textrm {mm}$ ones, as they prove more convenient to track the propagation of the dye.

4.1. Flow equations

Our experimental aquifer is made of a homogeneous and isotropic porous material. The Reynolds number of the flow is ${Re} = U d_p/\nu$, where $\nu$ 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}\ \textrm {cm}\ \textrm {s}^{-1}$ to $4\times 10^{-1}\ \textrm {cm}\ \textrm {s}^{-1}$, so that ${Re}$ varies between $10^{-2}$ and $0.4$. Groundwater thus flows according to the simplest version of Darcy's law:

(4.1)\begin{equation} {\boldsymbol q} ={-} K {\boldsymbol\nabla} \left ( \frac{p}{\rho g} + y \right), \end{equation}

where $\boldsymbol {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,

(4.2)\begin{equation} \nabla^2\phi = 0, \end{equation}

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

To solve it, we must complement (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:

(4.3)\begin{equation} \frac{\partial\phi}{\partial n} = 0, \end{equation}

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

(4.4)\begin{equation} \phi = y \end{equation}

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 just the water balance across the water table (Van de Giesen, Parlange & Steenhuis Reference Van de Giesen, Parlange and Steenhuis1994). In steady state, the latter reads

(4.5)\begin{equation} \frac{\partial \phi}{\partial y} - \frac{\partial \phi}{\partial x} \frac{\partial h}{\partial x} = \frac{R}{K} \end{equation}

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 (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\times 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 Reference Haitjema and Mitchell-Bruker2005; Guérin et al. Reference Guérin, Devauchelle, Robert, Kitou, Dessert, Quiquerez, Allemand and Lajeunesse2019). The ratio of these two quantities is therefore a small parameter, with respect to which we may linearize (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 Reference Dagan1964):

(4.6)\begin{equation} \frac{\partial \phi}{\partial y} = \frac{R}{K}\quad \textrm{in } y= 0. \end{equation}

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 Reference Polubarinova-Kochina1962). An advantage of doing so is that it provides simultaneously the hydraulic head $\phi$ and the stream function $\psi$. Accordingly, we introduce the complex coordinate $z = x + \textrm {i}y$ (with $i^2=-1$), and look for an analytical function $\varPhi (z) = \phi (z) + \textrm {i} \psi (z)$ that satisfies the boundary conditions. The stream function $\psi$, like the hydraulic head $\phi$, are rescaled with the hydraulic conductivity $K$ (they both have units of length).

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

(4.7)\begin{equation} \psi = 0, \end{equation}

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

(4.8)\begin{equation} \psi ={-} \frac{RL}{K} \left(1+\frac{x}{L}\right) \quad \textrm{in } y= 0. \end{equation}

This condition, together with (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 (§ 4.3). As expected for the stream function, the amplitude $[\psi ]$ of this discontinuity is the discharge of the aquifer, which matches the rainfall input: $KW[\psi ] = Q=RLW$.

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

(4.9)\begin{equation} \varPhi_k(z) = \frac{RL}{K} \cosh\left[k {\rm \pi}\left(\frac{\textrm{i}z}{L} - a\right) \right], \end{equation}

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

(4.10)\begin{equation} \varPhi(z) = \sum_{k=1}^{\infty}C_k \varPhi_k(z), \end{equation}

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.

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

Noting that, for each mode $k$, the stream function $\psi _k$ reduces to a cosine at the surface, we use a 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 Reference Farlow1993). We thus find the expression for the complex velocity potential in the form of a Fourier series:

(4.11)\begin{equation} \varPhi(z) ={-} 2\frac{RL}{K}\sum_{k=1}^{\infty} \frac{\cosh\left(k {\rm \pi}(\textrm{i}z/L - a) \right)}{k{\rm \pi} \sinh( k{\rm \pi} a)}. \end{equation}

The imaginary part of this velocity potential yields the stream function $\psi (x,y)$, the contours of which are streamlines. To draw them, we need to (1) define a numerical grid over the aquifer's extent, (2) evaluate the stream function on this grid, (3) extract the contours of this function. Accordingly, we start by building a triangular mesh of the aquifer, using the internal mesher of FreeFem++ (Hecht (Reference Hecht2012), Python routines available online: https://github.com/odevauchelle/pyFreeFem). Based on (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 \ell$, 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 Reference Hamming1973, p. 534). By trial and error, we opt for $\delta \ell /L = 1.25\times 10^{-2}$ and $100$ modes, a good trade-off between accuracy and computation speed.

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 (4.11). (a–c) Dimensionless travel time of rainwater, $\tau ^*$, as a function of relative landing position, $x_{start}/L$. Solid blue line: numerical approximation (§ 5.2). Black dashed line : expansion near the outlet, (5.2). (d–f) Streamlines (black dashed lines) and hydraulic head (blue shading).

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

(4.12)\begin{equation} V_e = W \int_{{-}L}^{0} h(x)\,{\textrm{d}x}. \end{equation}

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 with the exit condition takes over the flow (§ 4.3). Although unrealistic, this mathematical behaviour makes it possible for (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 set-up.

4.3. Flow near the outlet

The series expansion of § 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\rightarrow \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 Reference Polubarinova-Kochina1962), we find that the velocity potential then reads (figure 7a)

(4.13)\begin{equation} \varPhi \underset{H \to +\infty}{\longrightarrow}\varPhi_{deep} = \frac{2}{\rm \pi}\frac{RL}{K} \left[ \log\left( 1 - \exp({\textrm{i}{\rm \pi} z/L})\right) - \textrm{i} {\rm \pi}\left(\frac{z}{L}+1\right) \right]. \end{equation}

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

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$ is the inclination of the streamlines with respect to horizontal.

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

(4.14)\begin{equation} \varPhi \underset{z \to 0}{\longrightarrow} \varPhi_{outlet} = \frac{RL}{K} \left[ \frac{2}{\rm \pi} \log\left( \frac{z}{L} \right) +\textrm{i} \right]. \end{equation}

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 (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, means 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_{start}$ (figure 6):

(4.15)\begin{equation} x_{start} ={-}\frac{2 \theta}{\rm \pi} L. \end{equation}

One sometimes uses a series of boreholes to measure how the chemical composition of groundwater changes with depth (Haria & Shand Reference Haria and Shand2004). 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.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, ${\boldsymbol q}/s$ (Guyon, Hulin & Petit Reference Guyon, Hulin and Petit2001). Accordingly, the time that a particle needs to move between two points of complex coordinates $z_0$ and $z_1$ reads:

(5.1)\begin{equation} t = \frac{s}{K}\int_{z_0}^{z_1} \frac{\mathrm{d}l}{|\boldsymbol{\nabla} \phi|}, \end{equation}

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 (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 Reference Saffman1959). 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 § 6. First, however, we ignore its influence, and compute the travel time associated with advection only.

Using (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 (§ 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_\textrm {start}$, of the streamline along which it will travel underground. We may therefore calculate the travel time $\tau$ of a raindrop by integrating numerically (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.

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

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

(5.2)\begin{equation} \tau_{outlet} = \frac{\rm \pi}{4} \frac{sL}{R}\left(\frac{\epsilon}{L}\right)^2. \end{equation}

The radial symmetry of the logarithmic singularity makes $\tau _\textrm {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 (5.1), we first note that $|{\boldsymbol \nabla } \phi | = |\varPhi '|$ and $\textrm {d} l = |\textrm {d} z| = |\textrm {d} \varPhi | / |\varPhi '|$, where $'$ denotes the complex derivative. Along a streamline, the stream function is constant, while the hydraulic head decreases, so $|\textrm {d}\varPhi | = -\textrm {d}\phi$. Combining all these expressions, we rewrite (5.1) as

(5.3)\begin{equation} \tau_s ={-} \frac{s}{K}\int_{\phi_{in}}^{\phi_{out}} \frac{\textrm{d}\phi}{|\varPhi^{\prime}|^2}, \end{equation}

where $\phi _{in}$ and $\phi _{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 (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:

(5.4)\begin{equation} \varPhi_s = \phi_c \left[ 1- \frac{{\rm \pi}^2}{2} \left(\frac{z-z_s}{L}\right)^2 \right], \end{equation}

where $z_s$ is the position of the stagnation point ($z_s = -L-\textrm {i}H$ in the lower-left corner and $z_s =-\textrm {i}H$ in the lower-right one). The hydraulic head at the stagnation point, $\phi _c$, is that of the far field, namely $\phi _c = \Re (\varPhi (z_s))$.

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

(5.5)\begin{equation} \varPhi_s^{\prime} ={-} \sqrt{2}{\rm \pi} \frac{\phi_c}{L} ~\sqrt{1 - \frac{\varPhi_s}{\phi_c}}, \end{equation}

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

(5.6)\begin{equation} \tau_s = \frac{sL}{R} \frac{\sinh({\rm \pi} a)}{2{\rm \pi}} \left| \textrm{asinh}\left( \frac{ \phi_{in} - \phi_c}{\psi_0}\right) \right|, \end{equation}

where $\psi _0$ defines the streamline along which the water parcel travels.

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\times 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 Reference Zwietering1959; Nauman & Buffham Reference Nauman and Buffham1983). 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 § 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 (§ 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 colour 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, ${\rm \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 detail in § 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 behaviour of the travel time near the singularity, and near the outlet.

Figure 9.

Dimensionless travel time of a parcel of dye, $\tau ^*$, as a function of the distance to the outlet $-x_\textrm {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, ${\rm \Delta} \tau /\langle \tau \rangle$; horizontal ones indicate initial width of the parcel. Coloured lines: theoretical travel times after § 5.2 for aspect ratios $a=0.35$ (green), $0.7$ (blue, our laboratory experiment) and $1.4$ (orange).

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 (5.2), the duration of this travel decreases as the square of the distance that separates its landing point from the outlet:

(5.7)\begin{equation} \tau^* = \frac{\rm \pi}{4a}\left(\frac{x_{start}}{L}\right)^2. \end{equation}

Figure 10(a) shows that our experiment conforms to this asymptotic behaviour.

Figure 10.

Travel times through a rectangle aquifer. (a,b) Dimensionless travel time, $\tau ^*$, as a function of distance to outlet, $-x_{start}/L$, and distance to divide, $(L+x_{start})/L$. Coloured lines: theoretical travel times, after § 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, ${\rm \Delta} \tau /\langle \tau \rangle$, horizontal ones indicate the initial width of the parcel. (c,d) Probability density functions (PDF) of travel times, $f(\tau ^*)$, for several aspect ratios. Black dashed line: asymptotic regimes (5.10) and (5.11).

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 behaviour, 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 (4.11). Taking again advantage of this observation, we approximate the velocity potential by its first mode, $\varPhi _1$, and use (5.3) to calculate the corresponding travel time. Following a procedure similar to the one applied to the stagnation points in § 5.2, we find the following expression for the travel time near the drainage divide (appendix C):

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

which is a good approximation of the actual travel time (figure 10b). As expected, (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:

(5.9)\begin{equation} f(\tau^*) ={-} \frac{1}{L}\frac{\textrm{d}x_{start}}{\textrm{d}\tau^*}. \end{equation}

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 (5.9) numerically, based on the results of § 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 (5.9), yields the short-time behaviour of the PDF:

(5.10)\begin{equation} f(\tau^*) \sim \left( \frac{a}{\rm \pi}\right)^{1/2} \tau^{*-1/2}. \end{equation}

This expansion consistently represents the short-time behaviour of the numerical PDF (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 (5.9). We find that the tail of each PDF decreases exponentially, according to

(5.11)\begin{equation} f(\tau^*) \sim \frac{4{\rm \pi} a}{\sinh({\rm \pi} a) (1+\cosh({\rm \pi} a))}\exp({-\tau^*/\tau^*_c}), \end{equation}

with a characteristic dimensionless time

(5.12)\begin{equation} \tau^*_c = \frac{\sinh({\rm \pi} a)}{{\rm \pi} a}. \end{equation}

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 (figures 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.