On the development and analysis of coupled surface–subsurface models of catchments. Part 3. Analytical solutions and scaling laws

Abstract

The objective of this three-part work is to formulate and rigorously analyse a number of reduced mathematical models that are nevertheless capable of describing the hydrology at the scale of a river basin (i.e. catchment). Coupled surface and subsurface flows are considered. In this third part, we focus on the development of analytical solutions and scaling laws for a benchmark catchment model that models the river flow (runoff) generated during a single rainfall. We demonstrate that for catchments characterised by a shallow impenetrable bedrock, the shallow-water approximation allows a reduction of the governing formulation to a coupled system of one-dimensional time-dependent equations for the surface and subsurface flows. Asymptotic analysis is used to derive semi-analytical solutions for the model. We provide simple asymptotic scaling laws describing the peak flow formation, and demonstrate its accuracy through a comparison with the two-dimensional model developed in Part 2. These scaling laws can be used as an analytical benchmark for assessing the validity of other physical, conceptual or statistical models of catchments.

1. Introduction

In this third and final part of our work, we leverage the parametric study (Part 1: Morawiecki & Trinh 2024a) and two-dimensional (2-D) benchmark models (Part 2: Morawiecki & Trinh 2024b) to perform an in-depth asymptotic analysis of a coupled surface–subsurface model of a catchment. We specifically focus on flow within an aquifer characterised by a thin porous layer. The system begins in a steady state, for a constant precipitation $r_0$, which is characterised by an initial seepage zone. Our objective is to understand the response of the catchment, when subjected to intense rainfall $r>r_0$.

One of the main conclusions from the numerical simulations in Part 2 is that in early time, the river inflow rapidly increases as a result of the rainfall accumulating over the initial seepage zone. It eventually reaches a critical flow, followed by a much slower rise in the river inflow caused by the expansion of the seepage zone (see figure 1). One of the primary results we derive in this work, is an analytical solution for both early and late times, including a simple analytical formula for the critical flow:

(1.1)\begin{equation} Q_{\mathit{crit}} = \underbrace{K_s S_x L_z}_{\mathit{groundwater\ flow}} + \underbrace{r L_x \left(1 - \frac{K_s S_x L_z}{r_0 L_x}\right)}_{\mathit{overland\ flow}}. \end{equation}

This formula consists of a contribution from the groundwater flow, and a contribution from the overland flow formed over the seepage zone. The parameters $L_x$, $L_z$ and $S_x$ relate to geometrical features of the hillslope, $K_s$ corresponds to the soil hydraulic conductivity, whereas $r_0$ and $r$ represent the initial and simulated rainfall intensity, respectively.

Figure 1. (a) Studied hillslope geometry; initially groundwater and surface water are in steady state for precipitation rate $r_0$, and therefore river inflow (per unit length) is $Q(0)=r_0L_x$. The simulated rise of river inflow caused by a constant rainfall $r>r_0$ is presented in (b). Note the characteristic fast rise of the river inflow to $Q_{\mathit {crit}}$ at $t\in [0,t_{\mathit {crit}}]$.

We argue that such analytical scaling laws are valuable, both as a tool to diagnose the correctness of other, more complex rainfall-runoff models, and also as a measure for characterising a catchment's propensity to flooding.

The introduction of Part 1 covers the general subject of hydrology and parameter estimation, whereas the introduction of Part 2 covers computational rainfall-runoff models. We begin by discussing the content of this current paper, focusing on asymptotics and scaling laws, in the context of the existing literature.

1.1. On the importance of analytical benchmarks results

Catchment hydrology is one of many areas of engineering where numerical approaches tend to dominate over analytical ones. Due to the complex multiscale nature of hydrology, limited data availability and high computational cost, formulating and solving the correct equations can be enormously challenging (Grayson, Moore & McMahon 1992). Instead of the physical models based on the fundamental laws of hydrodynamics, simpler models such as conceptual and statistical models are often used instead (Moore et al. 2007). They are usually developed following a trial-and-error approach to fit available real-world data. However, these models do not provide any guarantee of model performance when applied to situations that may be underrepresented or missing in the training datasets (Parkin et al. 1996; Bathurst et al. 2004; Beven 2018).

To better understand the theoretical limits of different classes of catchment models, it is crucial to have a solid understanding of the different processes characterising the physical models. This goal can be achieved by developing numerical benchmark scenarios, as done by, e.g. Sulis et al. (2010) and Maxwell et al. (2014), and using them to compare predictions of different models. One perspective is that, in order to ensure consistency between the models over a wide range of input parameters, models should predict the same scaling laws for key features. For example, if one model predicts that the peak river flow, $Q$, is proportional to the catchment area, $A$, so $Q\propto A$, whereas another model predicts that $Q\propto \sqrt {A}$, then regardless of the fitting of these models, they cannot give consistent predictions over the entire range of $A$ values. For example, the second model fitted to a training dataset dominated with measurements conducted in large catchments would tend to overestimate flow in the case of small, often ungauged, catchments.

The above perspective may seem oversimplistic. However, we argue that despite many statistical works demonstrating such scaling laws (cf. Cunnane 1987; Kjeldsen, Jones & Bayliss 2008) the comparison of asymptotic scaling laws between different catchment models is not a commonly used approach in the modern hydrology (although see Vieira (1983) for a comparison of Saint Venant approximations, and Cook, Knight & Wooding (2009) for comparison of Richards and Boussinesq-based models).

The emphasis of our work, here, on deriving of such analytical scaling laws for flow prediction in the case of coupled surface–subsurface flows in catchments. Not only do these scaling laws provide clear guidance on the key dependencies of model parameters, but their analysis often illuminates further model simplifications.

Next, we provide a brief overview of the existing analytical theory, highlighting the lack of research on fully coupled surface–subsurface systems.

1.2. On analytical solutions in catchment hydrology

The typical governing equations used in physical catchment models include Richards equation for the subsurface flow (or the Boussinesq equation for the unconfined groundwater flow) and the Saint Venant equations for the overland and channel flows (see, e.g. the review by Shaw et al. 2010). These equations and their simplifications have been well-studied using analytical methods, although largely in an uncoupled manner. One of the equations studied is the Boussinesq equation, which is commonly used to determine the shape of the groundwater table. In this paper by Boussinesq equation we refer to the Boussinesq equation used to describe unconfined groundwater flow (see, e.g. Hálek & Švec 2011, chap. 2; Troch et al. 2013), in order to distinguish it from other equations and approximations known by the same name. In some cases, analytical solutions for this equation can be derived. Examples include for example steady-state groundwater flow and evolution in one-dimensional (1-D) hillslopes (Polibarinova-Kochina & Wiest 1962; Troch et al. 2013). Similarly, analytical solutions have been developed for the 1-D Richards equation to describe water transfer through the unsaturated soil under constant and time-varying infiltration (Warrick, Lomen & Islas 1990).

For the case of overland flow over a hillslope, analytical solutions have been found for a kinematic approximation of the Saint Venant equations, as done, e.g. by Parlange, Rose & Sander (1981) and Tao, Wang & Lin (2018). These analytical solutions and their approximations are important as they provide benchmark models for testing numerical schemes (e.g. benchmarks by MacDonald et al. (1995) for the overland flow and by Tracy (2006) for the subsurface flow), and can be used to develop less computationally demanding modelling approaches such as TOPMODEL by Kirkby & Beven (1979). However, no analytical solutions have been found so far coupled systems that include governing equations describing both subsurface/groundwater and surface flow. Despite the importance of these models in catchment hydrology, the study of these models has been restricted to numerical solutions only (Maxwell et al. 2014).

1.3. On the shallow-water approximation for subsurface flow

Previously, in Part 2, we introduced deep-aquifer scenario to describe a catchment with a deep aquifer, and shallow-aquifer scenario to describe a catchment where the subsurface flow is predominantly transferred through a thin porous layer near the surface. Mathematically, shallow-aquifer scenario is the limiting case of deep-aquifer scenario in which the aquifer depth is much smaller than the catchment width, $L_z\ll L_x$. We showed that in both cases, under standard initial conditions, the full three-dimensional (3-D) catchment model can be reduced to a simpler 2-D hillslope model. We shall begin with this 2-D assumption as the basis in this paper. We then demonstrate that in the shallow-water scenario, under certain assumptions, the 2-D model can be further reduced to a 1-D model.

The reduction of the 2-D subsurface flow into a 1-D model is not a new concept. The simplification is based on the Dupuit–Forchheimer approximation by Dupuit (1863) and Forchheimer (1914), which states that the groundwater flow is predominantly horizontal, and that the total flow scales proportionally with the saturated aquifer thickness. Boussinesq (1877) used this assumption to develop a 1-D model for the groundwater height; this is now known as the Boussinesq equation for groundwater flow (see, e.g. Bartlett & Porporato 2018), or the Dupuit–Boussinesq equation (see, e.g. Guérin, Devauchelle & Lajeunesse 2014). As we show later, it can be derived from the 2-D Richards model under the aforementioned assumption that $L_z\ll L_x$. The accuracy of this approximation is studied in detail by Paniconi et al. (2003) and Cook et al. (2009).

The Boussinesq equation is commonly used in groundwater modelling, and a wide class of analytical and approximate solutions has been developed. Notable examples are reviewed by Wooding & Chapman (1966), Anderson & Brooks (1996), Troch, Paniconi & Emiel van Loon (2003) and Bartlett & Porporato (2018). These studies, however, concern only groundwater flow, and do not involve the coupling with the overland flow, which is an essential component of a standard physical catchment model.

Here, we extend these studies by including the effect of overland flow in the Boussinesq equation. Our main result in this paper is the derivation and analysis of the following 1-D dimensionless coupled surface/subsurface model:

(1.2) \begin{equation} \frac{\partial H}{\partial t}=\begin{cases} \displaystyle f(x)^{{-}1}\left[\frac{\partial }{\partial x}\left(\sigma H\frac{\partial H}{\partial x} + H\right) + \rho_0 r(x, t)\right] & \text{if } H \leq 1,\\ \displaystyle \frac{\partial }{\partial x}\left(\sigma\frac{\partial H}{\partial x} + \mu\left(H-1\right)^k\right) + \rho_0 r(x, t) & \text{if } H > 1, \end{cases} \end{equation}

where $f(x)$, $\sigma$, $\mu$ and $\rho _0$ are dimensionless parameters explained in detail in § 2, and $H(x,t)$ is the total height of groundwater and surface water, which depends on the distance from the channel $x$ and time $t$. Values $H\leq 1$ represent unsaturated soil without surface water, and $H>1$ represent saturated soil with surface water. The main difference from the classical Boussinesq equation is the second case in the above equation with $H>1$, in which we include an additional term representing the overland flow.

We use the above 1-D coupled surface–subsurface model to develop analytical solutions for the river flow formed by rainfall of a constant intensity (however, the result can be generalised for time-dependent rainfall). Our methodology takes advantage of the negligibly small size of the diffusion terms in most of the seepage zone, which allows us to use the method of characteristics for the study of wave propagation. This approach is similar to previous kinematic treatments of the Saint Venant equations (Henderson & Wooding 1964; Woolhiser & Liggett 1967), but for our problem, the size of the seepage zone increases as a result of rising groundwater, which introduces a secondary dynamics.

The analytical approximations we develop in this work are possible due to a few governing assumptions; they are supported by the analysis of the typical values of catchment parameters described in Part 1 (Morawiecki & Trinh 2024a). The main approximations are: (i) the typical rainfall duration is much shorter than the characteristic timescale of groundwater flow; (ii) the typical timescale of surface flow is much shorter than that of subsurface flow; and (iii) the mean precipitation rate is larger than the maximal groundwater flow passing through the saturated zone.

We start by introducing the above 1-D model in § 2, with its typical dynamics discussed in § 3. For a scenario of single intensive rainfall described in § 4, we find an approximated analytical form of the initial steady state in § 5, followed by a short-time asymptotic analysis in § 6. The accuracy of the developed analytical approximations is assessed in § 7. In § 8, we highlight key hydrograph features predicted by this analytical solution, followed by conclusions in § 9 and further discussion in § 10.

2. Formulation of the 1-D coupled model

In this section, we introduce a 1-D model describing the horizontal groundwater and overland flow along the hillslope, firstly in a dimensional and then in a dimensionless form. Its formal mathematical derivation from the 2-D benchmark model introduced in our previous paper is presented in Appendix B: here we focus on presenting the general structure of the model instead.

2.1. Dimensional model

Let us consider a 2-D hillslope of length $L_x$ with a uniform terrain slope $S_x$, uniform thickness of the porous layer $L_z$, uniform saturated soil hydraulic conductivity $K_s$ and an impenetrable bedrock beneath the hillslope. As shown in the $(x,z)$-plane in figure 2, we denote the thickness of the saturated zone as $H_g(x,t)$ and the height of the surface water as $h_s(x,t)$.

Figure 2. Hillslope geometry used to formulate a 1-D surface–subsurface model.

We shall assume that overland flow can only occur when the soil becomes fully saturated. The overland flow generated by exceeding the soil infiltration capacity (Kirkby 2019) is not considered. Under this assumption, the heights $H_g$ and $h_s$ can be combined to form a single dependent variable,

(2.1)\begin{equation} H(x,t) = H_g(x,t) + h_s(x,t), \end{equation}

defined as the total height of groundwater and surface water. We now review the governing equations for $H$, for which the details are presented in Appendix B.

2.1.1. Groundwater flow

The standard approach to model shallow-water aquifers uses the Dupuit–Forchheimer assumption, which states that groundwater flows horizontally with the pressure head following a hydrostatic profile. Under this approximation, the groundwater flow is given by

(2.2)\begin{equation} Q_g = K_s H\frac{\partial H}{\partial x} + K_s H S_x. \end{equation}

When the soil is not fully saturated, and hence $H< L_z$, the evolution of the groundwater depth is given by the continuity equation, resulting in a standard form of the Boussinesq equation (Troch et al. 2013) for an unconfined aquifer:

(2.3)\begin{equation} f(x, t) \frac{\partial H_g}{\partial t} = \frac{\partial }{\partial x}Q_g + r = \frac{\partial }{\partial x}\left(K_s H_g\frac{\partial H_g}{\partial x} + K_s H_g S_x\right) + r(x,t), \end{equation}

where $r(x,t)$ denotes the groundwater recharge, and $f(x, t)$ is a drainable porosity. In this paper, we assumed that the recharge is equal to the precipitation.

The drainable porosity function is more subtle; as introduced in Appendix C, it is formally defined as the rate of change of groundwater volume, $\mathcal {V}$, given a change in the groundwater level, $H$, i.e. $f={\rm d}\mathcal {V}/{\rm d} H$. Note that $f$ depends on the soil saturation above the groundwater table. For example, higher soil saturation implies that less water is required to raise the groundwater by a given volume and, hence, $f$ is lower. Recall that the soil saturation, $\theta$, is computed as a function of the pressure head, $h_g$, as given by the Mualem–van Genuchten model (C3). In theory, computing $f$ would involve coupling equation (2.3) with a model for $h_g(x, z,t)$.

In the literature (e.g. Troch et al. 2003), $f$ is often assumed to be a parameter with a value specific to the soil type at a given location. In practice, however, $f$ can change over time. For example, during a rainfall, a characteristic wetting front forms, which slowly propagates from the surface towards the groundwater table (Caputo & Stepanyants 2008). In order for this table to rise, the rainwater must first infiltrate through the unsaturated zone. This causes $f$ to significantly change over time.

In this paper, we approximate $f(x, t)$ by a time-independent mean drainable porosity $f_{{mean}}(x)$. Although the model results will not duplicate the full time-dependent behaviour observed in Part 2, the mean value, $f_{{mean}}(x)$, is chosen such that solutions correctly capture the key time when the groundwater reaches the land surface. We show later that the resultant model reproduces the hydrograph during a single precipitation event (see figure 9 later).

In Appendix C, we justify the choice of

(2.4)\begin{equation} f(x, t) \approx f_{{mean}}(x) \equiv \frac{v_H(x)}{D(x)}, \end{equation}

where $v_H(x)$ is the initial drainable volume per unit area at a given location, and $D(x)=L_z-H(x)$ is the depth of the groundwater below the land surface. In other words, the drainable porosity is given by the fraction of the soil volume that can be filled with water. Henceforth, we take the mean approximation (2.4) as the definition of $f$. Further discussion of the drainable porosity function is provided in Appendix C, where we provide formulae for its implementation.

2.1.2. Coupling with the overland flow

Now let us consider the case where the soil is fully saturated, for which $H_g = L_z$, which implies $H>L_z$. In this case, as derived in Appendix B, the surface depth $h_s$ evolves according to the following continuity equation:

(2.5)\begin{equation} \frac{\partial h_s}{\partial t} = \frac{\partial }{\partial x}(Q_g + Q_s) + r = \frac{\partial }{\partial x}\left(K_s L_z \frac{\partial h_s}{\partial x} + \frac{\sqrt{S_x}}{n_s} h_s^k\right) + r(x,t). \end{equation}

In the above equation, we have used Manning's equation to represent the overland flow:

(2.6)\begin{equation} Q_s = \frac{\sqrt{S_f}}{n_s} h_s^k \sim \frac{\sqrt{S_x}}{n_s} h_s^k, \end{equation}

where $k=5/3$ and $n_s$ is the Manning roughness coefficient, which depends on the hillslope surface type and is determined empirically. We use the kinematic approximation ($S_f \sim S_x$), where the friction slope $S_f$ is only dependent on the elevation gradient $S_x$. Alternatively, we could also consider the diffusive approximation $S_f\sim S_x+{\partial h_s}/{\partial x}$; however, we limit this study to the kinematic approximation only for simplicity.

Furthermore, we note that apart from the constant gravitationally-induced groundwater flow, the pressure difference caused by the gradient of surface water height may affect the groundwater flow. Typically, the size of the overland flow is negligibly small compared with the thickness of the porous layer. However, as we discuss in § 5, this approximation fails at the propagating seepage front, which requires us to include the $({\partial }/{\partial x})(K_s L_z ({\partial h_s}/{\partial x}))$ diffusion term.

Now, we can combine (2.3) and (2.5) into a single equation for $H$:

(2.7) \begin{equation} \frac{\partial H}{\partial t} =\begin{cases} \displaystyle f(x)^{{-}1}\left[\frac{\partial }{\partial x}\left(K_s H\frac{\partial H}{\partial x} + K_s H S_x\right) + r\right] & \text{if } H \leq L_z, \\ \displaystyle \frac{\partial }{\partial x}\left[K_s L_z \frac{\partial H}{\partial x} + \frac{\sqrt{S_x}}{n_s} (H-L_z)^k\right] + r & \text{if } H > L_z. \end{cases} \end{equation}

We assume a no-flow boundary condition at the catchment boundary ($x=L_x$). At the location of the river ($x=0$), we assume that the river table is located at the same level as the overland water height (or at the surface if no overland flow is present). Therefore, at $x=0$ we set $H_g=L_z$ and a (flat) free-surface condition for the overland flow, $\partial h_s/\partial x = 0$. In addition, in this work, we study the time evolution of the above system assuming that it is initially in equilibrium for a given mean rainfall $r_0< r$, and then subjected to a rainfall $r$ for $t > 0$. Therefore, for the initial condition, we take the steady-state of equation of (2.7) for a given mean rainfall $r_0$. The boundary and initial conditions are summarised shortly in § 2.3.

2.1.3. Channel flow

As was discussed in § 3.3 of Part 2, the total groundwater and overland flow that is reaching the riverbank form a channel flow. This flow can be described by 1-D Saint Venant equations. However, in this study, our focus is on studying the properties of river inflow from the hillslope, not the subsequent channel flow. Therefore, for the purpose of this study, we assume that the river height is constant, limited to the depth of the channel. Analysing how the river inflow propagates thought the channel (or the entire drainage network) and how the drying of the channel affects the surface and subsurface flows can be interesting extensions of this study.

2.2. Non-dimensional model

The above model can be non-dimensionalised by taking $x = L_x x'$, $t = T_0 t'$, $H=L_z H'$ and $r=r_0r'$. Here, $T_0= L_x/(K_sS_x)$ is a characteristic timescale of the groundwater flow, chosen to balance the temporal term and the $\partial _x(K_s H S_x)$ term.

Once non-dimensionalised, our governing equations (2.7) become (after dropping primes)

(2.8a) \begin{equation} \frac{\partial H}{\partial t} =\begin{cases} \displaystyle f(x)^{{-}1}\left[\frac{\partial }{\partial x}\left(\sigma H\frac{\partial H}{\partial x} + H\right) + \rho_0 r(x, t)\right] & \text{if } H \leq 1, \\ \displaystyle \frac{\partial }{\partial x}\left(\sigma\frac{\partial H}{\partial x} + \mu (H-1)^k\right) + \rho_0 r(x,t) & \text{if } H > 1, \end{cases} \end{equation}

and the dimensionless parameters $\sigma$, $\mu$ and $\rho _0$ are introduced shortly in § 2.3. In this work, we assume that the rainfall is constant and uniform, i.e. $r(x, t) = r = \text {constant}$, except for the initial jump from $r_0$ to $r>r_0$. However, we shall discuss in § 10 that our methodology can be applied in the case of time-dependent rainfall.

In combination with the governing equations (2.8a), we have to specify the dimensionless boundary conditions. First, at $x = 0$, we need to consider two situations. First, if a seepage zone exists for the initial $r_0$, we set a free flow boundary condition,

(2.8b)\begin{equation} H_x(0, t) = 0. \end{equation}

However, as we shall demonstrate in § 3.1, if $r_0$ is low enough, initially the seepage does not exist. Then we assume that $H(x,t)$ representing groundwater is reaching the surface at $x=0$, i.e.:

(2.8c)\begin{equation} H(0, t) = 1. \end{equation}

During a rainfall ($r>r_0$), the groundwater gradient at $x=0$, which is initially negative, increases as the groundwater rises. The seepage starts to grow, when the when it becomes positive, which is when boundary condition (2.8c) is replaced with (2.8b).

At the right-hand edge, by the definition of a catchment, there is zero flow:

(2.8d)\begin{equation} Q(1, t) = 0, \end{equation}

where the dimensionless total flow, $Q$, is defined as

(2.9) \begin{equation} Q(x, t) =\begin{cases} \displaystyle H + \sigma H\frac{\partial H}{\partial x} & \text{if } H \leq 1, \\ \displaystyle 1 + \sigma\frac{\partial H}{\partial x} + \mu(H-1)^k & \text{if } H > 1. \end{cases} \end{equation}

For the initial condition, $H(x,t=0)=H_0(x)$, we take a steady state of (2.8a) for $r=1$:

(2.10) \begin{equation} 0 =\begin{cases} \displaystyle\frac{\partial }{\partial x}\left(\sigma H_0\frac{\partial H_0}{\partial x} + H_0\right) + \rho_0 & \text{if } H_0 \leq 1, \\ \displaystyle\frac{\partial }{\partial x}\left(\sigma\frac{\partial H_0}{\partial x} + \mu \left(H_0-1\right)^k\right) + \rho_0 & \text{if } H_0 > 1. \end{cases} \end{equation}

In this paper, we refer to this model (2.8) as the 1-D model.

2.3. The non-dimensional parameters

In the first case of (2.8), we have introduced two key dimensionless parameters, defined as

(2.11a)\begin{gather} \sigma =\frac{L_z}{L_xS_x} = \frac{\textrm{thickness of the porous layer}}{\textrm{elevation drop along the hillslope}}, \end{gather}

(2.11b)\begin{gather}\rho_0 =\frac{r_0 L_x}{L_z S_x K_s} = \frac{\textrm{precipitation flux}}{\textrm{maximum groundwater flux}}. \end{gather}

Note that $\sigma \to \infty$ as the hillslope becomes increasingly flat. The parameter $\rho _0$ represents the ratio of the total precipitation rate (given by $rL_x$ in $\mathrm {m}^2\ \mathrm {s}^{-1}$) to the maximum possible groundwater flow for fully saturated soil (given by $L_z S_x K_s$ in $\mathrm {m}^2~\mathrm {s}^{-1}$). Introduction of the maximum groundwater discharge is a classic concept in hydrology, see, e.g. Horton (1936).

In the second case of (2.8), we have introduced an additional dimensionless parameter:

(2.12)\begin{equation} \mu = \frac{L_z^{k-1}}{K_sS_x^{1/2}n_s}. \end{equation}

We argue later in § 6.2 that the characteristic size of the overland flow scales as $L_s=\mu ^{-1/k}L_z$. Following that section, we introduce a key dimensionless parameter to describe the dynamics in the seepage zone, namely the Péclet number:

(2.13)\begin{equation} Pe =\frac{\mu^{1/k}}{\sigma}=\frac{\dfrac{\sqrt{S_x}}{n_s}L_s^{k}}{K_sL_z\dfrac{L_s}{L_x}} =\frac{\text{convective overland flow}}{\text{diffusive effect of the groundwater flow}}. \end{equation}

In order to interpret ${Pe}$, we note that the numerator represents the second term on the right-hand side of (2.7) for $H>L_z$, representing convective effects. The denominator represents the size of the first term on the right-hand side of (2.7) for $H>L_z$, representing diffusive effects.

Based on median values of physical parameters used in the above equations provided in table 1, we have $\sigma \approx 10^{-2}$, $\rho _0\approx 1.5$, $\mu \approx 10^7$ and ${Pe}\approx 10^5$. Consequently, our work primarily focuses on the limits of $\mu, \, {Pe} \to \infty$, where convection dominates diffusion.

Table 1. Typical values of physical parameters characterising UK catchments extracted in Part 1 of this paper (Morawiecki & Trinh 2024a).

3. Numerical methodology and typical dynamics

Recall that the solutions are essentially characterised by the overland and groundwater heights, and a quadruplet of parameters:

(3.1)\begin{equation} H = H(x, t; \sigma, \mu, \rho_0, r), \end{equation}

as well as the constant $k = 5/3$ appearing in Manning's formula. These solutions are found by solving the partial differential equation (PDE) (2.8a) subject to the boundary conditions (2.8b)–(2.8d). In the typical configuration, the flow transitions from overland ($H > 1$) to groundwater ($H \leq 1$) at a contact point $x = a(t)$, which is determined as part of the solution.

The model in (2.8a) was implemented in Matlab using the ode15s solver to find its steady state and pdepe to solve the time-dependent problem. We divide the spatial and temporal domain as follows:

(3.2)\begin{equation} x_i = \frac{i}{N_x} \quad \text{and} \quad t_j = \frac{j}{N_t} t_{max}, \quad i, j = 1, 2, 3, \ldots, \end{equation}

where we typically use $N_x = 200$ and $N_t = 300$. We check whether increasing mesh size and time resolution does not significantly affect the obtained solution, and if it does, we refine the mesh. The codes used to generate figures in this work are available in a GitHub repository (Morawiecki 2022). All numerical results in this paper were obtained for the values presented in table 1, unless stated otherwise.

3.1. Typical results for $\rho _0 < 1$ and $\rho _0 > 1$

The existence of the seepage zone in the initial steady state depends on whether the value of $\rho _0$, defined in (2.11b), is higher or lower than $1$, and each case exhibits a different transient behaviour. Typical solutions obtained in these two cases are shown in figure 3.

Figure 3. Schematic presenting early- and late-time dynamics for a catchment with and without an initial seepage zone (corresponding to $\rho _0>1$ and $\rho _0\leq 1$ respectively). In (a), (b) and (c) a seepage zone is observed, and therefore a separate graph is presented to illustrate the evolution of the surface water depth (note that the vertical axis is multiplied by the scaling factor $\mu ^{1/k}\approx 3260$).

First, consider figure 3(a,b). In the case of $\rho _0>1$, we already have a seepage zone in the initial state. In a short timescale, the height of the surface water increases quickly until reaching a seemingly quasi-static state. Afterwards, the flow continues to increase as a result of the saturation front propagating uphill, but this process is characterised by a much longer timescale and a slower rate of flow rise. The difference between the short- and long-time behaviour can also be seen in the produced hydrograph in figure 4. It shows the dependence between the total river inflow, defined as the total flow (2.9) evaluated at the river bank ($x=0$),

(3.3)\begin{equation} Q(t)=Q(x=0,t). \end{equation}

In the presented hydrograph, we have marked the initial fast transition as (A) and the subsequent slow transition as (B).

Figure 4. Schematic representation of the hydrograph obtained for a catchment with and without an initial seepage zone (corresponding to $\rho _0>1$ and $\rho _0\leq 1$, respectively). Points represent times for which the profiles are shown in figure 3, whereas letters A–D refer to corresponding phases from that figure.

For $\rho _0<1$, we do not observe an initial seepage zone, i.e. $H(x,t=0)<1$ for all $x$. For some time the groundwater table is rising, increasing groundwater flow reaching the river, until the groundwater depth gradient at $x=0$ becomes $0$. Then the seepage zone starts to slowly form and propagate away from the channel, increasing the overland flow reaching the river.

However, in practice, in the case of real-world catchments characterised by a thin porous layer (which is a base assumption behind the presented 1-D model), the groundwater flow rate is highly limited. Therefore, in such catchments, we expect the $\rho _0>1$ case to be more prevalent, which is additionally confirmed by low base flow index (BFI) values characterising low-productive catchments. (BFI describes the ratio between the base flow and total flow in the given catchment. High values refer to catchments dominated by the groundwater flow, whereas low values refer to catchments with a significant overland flow component.) Therefore, in this paper, we focus on discussing the mathematical properties of each phase presented in figure 4 only in the case of $\rho _0>1$.

In general, the solution for the PDE model (2.8a) can only be found numerically. However, by taking advantage of the typical sizes of dimensionless parameters, the model can be further simplified. Figure 5 shows the effect of dimensionless parameters, $\rho$, $\rho _0$, $\sigma$ and $\mu$ on the model's solution. The graphs in the left column show how the initial steady state $H_0(x)$ depends on the value of each parameter, whereas the graphs on the right present the effect of each parameter on the hydrograph $Q(x=0,t)$. The conclusions from this numerical experiment are as follows.

1. (i) Parameter $\rho _0$ (typical value $2$), which following (2.11b) characterises the mean precipitation rate in terms of the groundwater flow, has a significant effect on the initial steady state. As discussed in detail previously, $\rho _0<1$ corresponds to a hillslope with no initial seepage zone, and is characterised by different dynamics than the $\rho _0>1$ case, in which we observe a fast rise of flow in the early time. In most of this paper, we consider only the latter case.

2. (ii) Parameter $\rho$ (typical value $\approx$20), which characterises the simulated precipitation rate in terms of the groundwater flow, does not affect the initial steady state, but it does affect how quickly the flow is rising. Higher $\rho$ values lead to both a higher flow over the seepage zone and a faster growth of this zone. Values of $\rho$ can vary significantly depending on the rainfall event considered.

3. (iii) Parameter $\sigma$ (typical value $10^{-2}$), which following (2.11a) characterises the thickness of porous layer compared with the elevation drop along the hillslope, has a significant effect only on the solution outside the seepage zone. In the case of the seepage zone, the term including $\sigma$ is negligibly small compared with the $\mu$ term. However, it affects the speed at which the seepage zone is growing. Note that as $\sigma \rightarrow 0$, the initial groundwater shape becomes a linear function (with a possible small boundary layer at its left border). Even though this limit is not critical in our analysis, it may allow us to approximate the initial groundwater shape without the need to solve the governing equations numerically.

4. (iv) Parameter $\mu$ (typical value $10^6$), which following (2.12) characterises the overland flux, does not have a significant effect on the groundwater table outside the seepage zone, but it has a major effect on the height of the surface water within the seepage zone. Higher $\mu$ values correspond to lower surface water height, which in the limit $\mu \rightarrow \infty$ becomes negligible compared with the variation of the groundwater depth. In addition, in this limit, the seepage zone size reaches a limiting value, $a_0=1-{1}/{\rho _0}$ (see § 5). This limit is strongly supported by real-world data (typical value of $\mu$ for UK catchments is of the order of $10^6$ based on the typical parameter values estimated in Part 1 of this study), and it allows us to derive the formula for a typical hydrograph.

Figure 5. Effect of dimensionless parameters shown for the initial steady state $H_0(x)$ versus $x$ (insets left column) and for the hydrograph $Q(t)$ versus $t$ (insets right column). Inset (a) shows changing $\rho _0$; inset (b) shows changing $\rho$; inset (c) shows changing $\sigma$; and inset (d) shows changing $\mu$. In the case of (a), (c) and (d), solid lines represent solutions for $\rho _0=1.5$ (scenario with an initial seepage zone), and dashed lines represent solutions for $\rho _0=0.6$ (scenario without an initial seepage zone). The surface water ($H>0$) is magnified 1000 times.

4. The formulation of an asymptotic model for intense rain

In the previous section, we presented numerical simulations of the full PDE system (2.8) and showed that under certain parameter choices, the resulting hydrographs could be approximately classified into two behaviours as shown in figure 3. In particular, when the system is initiated with an initial seepage zone, i.e. $\rho _0 > 1$, then in response to an intense rainfall with $r > 1$, the river inflow rapidly increases over time. It is important to note that the duration of a standard intensive rainfall is much shorter compared with the typical timescale of the groundwater flow (i.e. typical travel time along the hillslope), which is approximately

(4.1)\begin{equation} T_0=\frac{L_x}{K_sS_x} \approx \text{1000 days}. \end{equation}

Our focus is to develop the short-time asymptotics to better understand this crucial response. Ultimately, we aim to derive an analytical solution for the river inflow, $Q(x=0,t)$. Based on the physical constraints, we are primarily interested in the following asymptotic limits:

(4.2)\begin{equation} \left.\begin{aligned} \text{small-time} : \quad & t \ll 1, \\ \text{convection-dominated flow in the seepage zone} : \quad & Pe \gg 1, \\ \text{intense rainfall} : \quad & \rho = r \rho_0 \gg 1. \end{aligned}\right\} \end{equation}

To begin, let us reformulate the governing system in terms of a boundary value problem. We assume that there exists a single contact line located at $x = a(t)$ where $H(a(t), t) = 1$. This configuration is illustrated in figure 2. From (2.8a), the evolution of these surfaces is governed by

(4.3a 4.3b)

Here, we have introduced $\rho =\rho _0 r$. Thus, we have a set of two time-dependent equations for $H$ that are second-order in space, along with an additional contact-line position $a(t)$. Consequently, we require five boundary conditions in addition to the initial condition. Two boundary conditions are needed at $x = 0, 1$, and three matching conditions are required at the interface, $x = a(t)$. In total, these conditions are

(4.4a,b)\begin{gather} \partial_x H(0, t) = 0, \quad \partial_x H(1, t) ={-}\sigma^{{-}1}, \end{gather}

(4.4c–e)\begin{gather}H(a^-, t) = 1, \quad H(a^+, t) = 1, \quad \partial_x H(a^-, t) = \partial_x H(a^+, t), \end{gather}

where $a^\pm$ corresponds to the right/left limits as $x \to a$.

The first two boundary conditions are obtained from (2.8b) and (2.8d). The next two boundary conditions arise from defining $a$ as the point where the groundwater table reaches the surface (i.e. where $H=1$). The last interface condition is a consequence of the continuity of flow given by (2.9). Note that a kinematic condition can be derived for the front position. Applying the chain rule to $H(a(t), t) = 1$, we have

(4.5) \begin{equation} \frac{\partial H}{\partial t} + \frac{\partial H}{\partial x} \dfrac{{\rm d}a}{{\rm d}t} = 0 \quad \text{at } x = a(t). \end{equation}

Following (2.10), we set the initial condition given by the steady state of (4.3a)–(4.3b) for $r=1$, which we denote as $H_0(x; \, \rho _0)$. Thus, $H_0$ satisfies

(4.6) \begin{equation} 0 =\begin{cases} (\sigma H_0H_0' + H_0)' + \rho_0 & \text{for } x > a(t=0), \\ (\sigma H_0' + \mu (H_0-1)^k)' + \rho_0 & \text{for } x \leq a(t=0), \end{cases} \end{equation}

where primes $(')$ denote differentiation with respect to $x$. Here, $a(t=0)$ can be regarded as an eigenvalue and determined from this initial condition.

In the next two sections, we use this model to derive an asymptotic solution for the hydrograph $Q(t)$. Our approach involves three main steps.

1. (i) First, in § 5, we study the initial state $H(x, 0) = H_0(x; \, r = 1)$, which is assumed to be the steady-state response to the rain input $r=1$. This is a complicated coupled overland–groundwater problem, but we are able to develop analytical approximations in the limit of $\mu \to \infty$ or equivalently $Pe \to \infty$.

2. (ii) Next, in § 6.1, we study the small-time response of the groundwater configuration and the propagation of the seepage zone relative to this initial steady state. At the time $t = 0$, the rainfall is set to $r > 1$, which causes the groundwater to rise and the seepage zone to shift. Analytical approximations can be developed for the case of ${Pe} \to \infty$ and for large rainfalls, $r \to \infty$.

3. (iii) Finally, in § 6.2, we develop an analytical approach for predicting the evolution of the overland flow, which leverages our analysis of the seepage zone propagation obtained in step (ii). This turns out to be a wave propagation study using the method of characteristics.

5. Asymptotic analysis of the initial condition, $H_0$, with $Pe \to \infty$

In our model, we assume that the system begins at the configuration that corresponds to the particular steady-state solution forced by the ‘typical’ rainfall, $\rho _0$.

By integrating (4.6) and applying the upstream boundary condition $q(1)=0$, we obtain

(5.1a 5.1b)

The fact that the limit $\mu \to \infty$ involves the Péclet number, defined as ${Pe} = \mu ^{1/k}/\sigma$ via (2.13), is not entirely obvious. Note that as $\mu \to \infty$, the dominant balance in the overland equation for $x \leq a$ indicates that $H_0 \sim 1$ in this limit. We rescale $x = aX$ and $H_0 = \mu ^{-1/k} g(X)$, obtaining, for $X \in [0, 1]$,

(5.2)\begin{gather} \frac{{Pe}^{{-}1}}{a}\frac{\partial g}{\partial X}+g^k=\rho_0(1-aX) - 1, \end{gather}

(5.3)\begin{gather}g'(0) = 0 \quad \text{and} \quad g(1) = 0. \end{gather}

In the limit ${Pe}\to \infty$, we note that naively, the diffusion term in (5.2) tends to zero. Then, since $g(1) = 0$, we can approximate $\rho (1 - a) \sim 1$, which gives the front position as $a\sim 1-1/\rho _0$. However, note that in this limit, the leading (outer) solution is given by $g \sim [\rho _0 (1 - aX) - 1]^{1/k}$ and, hence, exhibits an infinite gradient as $X \to 1$. Consequently, it is not obvious that the diffusion term can be neglected a priori as ${Pe}\to \infty$. The gradient of the solution exhibits a boundary layer and thus requires a matched asymptotics approach.

In Appendix D, we show that the contact line, $x = a$, and the gradient at the front can be expanded into an asymptotic expansion. In terms of the original $H_0$, this is

(5.4)\begin{equation} a\sim a_0+a_1{Pe}^{-\beta} \quad \text{and} \quad H_0'(a) \sim \left[-\frac{\rho_0 a_1}{\sigma}\right] {Pe}^{-\beta}, \end{equation}

where $\beta = k/(2k-1)$ and the leading-order contact position is indeed

(5.5)\begin{equation} a_0 = 1-\frac{1}{\rho_0}. \end{equation}

Note that increasing the rainfall rate, $\rho _0 \to \infty$, sends $a \to 1$, and overland water saturates the entire hillslope. In contrast, the limit $\rho _0 \to 1^+$ reduces the seepage zone size to zero, as anticipated in § 3. The correction factor of $a_1$ in (5.4) can be calculated as an eigenvalue via the solution of a boundary-value problem (cf. (D9a)). Finally, note that as ${Pe}^{-1} \to 0$, the gradient at the transition between overland and groundwater flows, $H_0'(a) \to 0$.

As we show in § 6.1, in order to find the speed of the seepage zone growth, we need to find the initial depth of the groundwater outside the seepage zone first. We can find this initial depth by solving (5.1a), which can be rearranged as

(5.6)\begin{equation} \sigma \dfrac{{\rm d}H_0}{{\rm d}\kern 0.06em x} = \frac{\rho_0(1-x)}{H_0(x)} - 1\quad \text{for } x\in[a_0,1], \end{equation}

with a boundary condition $H_0(a_0)=1$. This first-order nonlinear ordinary differential equation (ODE) does not have an explicit analytical solution; it can either be solved numerically, or we can investigate its shape in different limits.

5.1. Analytical solution as $\sigma \to 0$

One quite useful limit is to consider $\sigma \rightarrow 0$, corresponding to the infinitely thin porous layer limit. In Appendix E, we derive the outer asymptotic expansion for $H_0$ in terms of $\sigma$, (E2), its inner expansion around $x=0$ (E4), and finally match them to form the following composite approximation for $H_0(x)$:

(5.7)\begin{equation} H_0(x)=\rho_0(1-x+\sigma-\sigma \,{\rm e}^{-{(x-a_0)}/{\sigma}}), \quad \text{for } x \in [a_0, 1]. \end{equation}

As shown in figure 6, this asymptotic solution provides a good approximation of the groundwater shape both for small $\sigma$ values and, surprisingly, also for large $\sigma$ values. In the latter case, $H_0$ becomes a quadratic function,

(5.8)\begin{equation} 1-H_0\sim \frac{\rho_0}{2\sigma}(x-a_0)^2, \end{equation}

which is also a limiting behaviour of our matched asymptotic solution (5.7) as $x\to a_0$.

Figure 6. Comparison of groundwater depth given by (5.6) (full numerical solution) with the matched asymptotic approximation given by (5.7).

6. Short-time asymptotics

This section relates to the asymptotic limits of $t \to 0$, $\mu \to \infty$ and $\rho = \rho _0 r \to \infty$ in (4.2).

6.1. Groundwater rise and propagation of the seepage zone

Having derived certain analytical properties of the steady-state configuration, $H_0(x; \, \rho _0)$ (used as an initial condition), we can now study the short-time behaviour of the system as the rain input is set to $\rho$. As argued at the start of § 4, this is a good approximation, when the rainfall duration is much shorter than the characteristic time of the groundwater transfer to the channel.

As shown in Appendix F, the outer solution outside the seepage zone (4.3b) can be expanded into a regular series expansion in powers of time $t$:

(6.1)\begin{equation} H_{outer}(x,t) \sim H_0(x; \, \rho_0) + \left[\frac{\rho - \rho_0}{f(x)}\right] t + {O}(t^2). \end{equation}

This approximation assumes that $x - a(t) = {O}(1)$. The approximation (6.1) thus indicates that the groundwater rises in a fashion proportional to time and the difference between current and prior rain input; it correctly describes the shape of the groundwater except for a thin boundary layer at $x=a(t)$ of thickness of ${O}(\sqrt {t/(\rho - \rho _0)})$ (see figure 7a). Therefore, for intense rainfall, $\rho \gg 1$, we can neglect the effect of this boundary layer.

Figure 7. (a) Comparison of the numerical solution of the groundwater shape (solid lines) with the outer solution developed in Appendix F (dashed lines) at different times $t$. The corresponding size of the seepage zone is presented in (b). A small region is magnified to highlight differences between the presented approximations. The lines are not smooth due to the $h(x)$ interpolation error.

We may use the outer groundwater approximation, (6.1), in order to predict the motion of the contact line, $x = a(t)$. Setting $H_{outer} = 1$ gives, in implicit form,

(6.2)\begin{equation} t \sim \frac{f(x=a(t))}{\rho - \rho_0} (1 - H_0(x=a(t))) \equiv \mathcal{T}(x=a(t)). \end{equation}

In order to calculate the above, we must solve two first-order ODEs: (5.6) for the height, $H_0(x)$, and (C2) for the head, $h_g(\hat {z})$, itself used in the calculation of $f(x)$. Alternatively, one can use the analytical approximations for $H_0(x)$ given by (5.7), and $f(x)$ given by (C6) or (C7). Figure 7(b) compares these approximations with the location of the saturation front computed from a full numerical solution of the 1-D model. As we observe based on the difference between the full numerical solution and leading-order approximation (LOA), neglecting the boundary layer around $a(t)$ introduces a small error when estimating the seepage zone size. Replacing the ODEs with analytical approximations for $f(x)$ and $H_0(x)$ in (6.2) also introduces an error, but it is significantly smaller.

6.2. Evolution of the overland flow

Now, knowing how the seepage zone propagates, we can develop a time-dependent solution for the overland flow. Our goal is to extract how the overland flow into the river, $Q_s(x=0, t)$, evolves in time, taking into account the effects of increased rainfall and the seepage zone growth.

6.2.1. Problem reduction under ${Pe}\rightarrow \infty$ limit

The equation for overland flow is given by (4.3a) with the initial condition satisfying steady state (5.1b). We rescale according to

(6.3)\begin{equation} \eta = \mu^{1/k}(H-1) \quad \text{and} \quad T = \mu^{1/k}t. \end{equation}

Here, $\eta = \eta (x, T)$ is the rescaled surface water height $h_s=H-1$. Then (4.3a) can be written as

(6.4)\begin{equation} \frac{\partial \eta}{\partial T} - k \eta^{k-1}\frac{\partial \eta}{\partial x} - {Pe}^{{-}1}\frac{\partial^2 \eta}{\partial x^2}= \rho, \end{equation}

and (5.1b), which provides the initial condition, $H_0$, is

(6.5)\begin{equation} 1 + \eta^k - {Pe}^{{-}1} \dfrac{{\rm d}\eta}{{\rm d}\kern 0.06em x} = \rho_0 (1 - x), \end{equation}

where, as before, ${Pe}^{-1}=\sigma /\mu ^{1/k}$. Following (4.3a) and (4.3c) the boundaries conditions are

(6.6a,b)\begin{equation} \partial_x \eta(0, t) = 0, \quad \eta(a(t), t) = 0. \end{equation}

Note that the characteristic time it takes the overland flow to reach the channel ($\mu ^{-1/k}T_0\approx 0.1\ \textrm {day}$) is much shorter than the characteristic time describing the groundwater flow ($T_0\approx 1000\ \textrm {days}$), and has a similar order of magnitude as a typical rainfall duration. As a result, a short-time approximation is not satisfactory to describe flow variation during a single rainfall event.

The solution for general times can be obtained by considering the ${Pe} \rightarrow \infty$ limit, similarly as we did when analysing the steady state. This limit allows us to neglect the diffusion term everywhere except for a negligibly thin boundary layer around $x=a$.

In the limit ${Pe} \rightarrow \infty$, we expand $\eta = \eta _0 + {Pe}^{-1} \eta _1 + \cdots$, and (6.4) becomes a first-order hyperbolic PDE:

(6.7a)\begin{equation} \frac{\partial \eta_0}{\partial T} - [k \eta_0^{k-1}]\frac{\partial \eta_0}{\partial x}= \rho, \quad x \geq 0. \end{equation}

For $0 \leq x \leq a(t)$, there is an initial condition given by

(6.7b)\begin{equation} \eta_0(x, 0) = \rho_0^{1/k}(a_0-x)^{1/k}, \end{equation}

where we have used the fact shown in Appendix D that $a_0=1-1/\rho _0$ (cf. (5.5)). The above initial condition is defined along the entire initial seepage zone, $x\in [0,a_0]$. Note that neglecting the diffusion term results in a kinematic wave equation, for which the downstream boundary condition (6.6a) is no longer required.

6.2.2. Implicit solution using methods of characteristics

The system (6.7) can be solved using the method of characteristics (Lagrange–Charpit equations). The solution is given by characteristic curves $(T, x, \eta _0)$, now parameterised by $(s, \tau )$, where $\tau$ is the characteristic curve parameter, and $s$ parameterises the initial data. The characteristic equations are

(6.8a) \begin{equation} \dfrac{{\rm d}T}{{\rm d}\tau}=1, \quad \dfrac{{\rm d}\kern 0.06em x}{{\rm d}\tau}={-}k\eta_0^{k-1}, \quad \dfrac{{\rm d}\eta_0}{{\rm d}\tau}=\rho. \end{equation}

The initial conditions are specified along $\tau = 0$ according to two types of characteristics. One set of characteristics emerges from $T = 0$, at the location of the initial water shape, $H_0(x)$, valid for $x\in [0, a(t)]$. Another set of characteristics emerges from the propagating front, $x = a(t)$, representing the groundwater reaching the surface and, hence, initiating surface flow.

Parameterising the initial data by $x = s$, we have

(6.8b) \begin{equation} (T(s, 0), x(s, 0), \eta_0(s, 0)) =\begin{cases} (0, s, H_0(s)), & s \in [0, a(t)], \\ (\mu^{1/k}\mathcal{T}(s), s, 0), & s \in [a(t), \infty). \end{cases} \end{equation}

The first condition will use the initial surface height, $H_0(s) = \rho _0^{1/k}(a_0-s)^{1/k}$ given by (6.7b). The second condition is essentially specified along the moving front, $(T, x, \eta _0) = (T, a(t), 1)$, but we have written it in terms of the $s$-independent variable, and the rescaled function $\mathcal {T}$ in (6.2). In summary, the characteristic solution can be obtained via direct integration of (6.8), giving

(6.9) \begin{equation} \left.\begin{array}{c@{}} T(s,\tau)=T(s,0)+\tau, \\ x(s,\tau)=x(s,0)-\rho^{{-}1}[\eta_0(s,0)+\rho\tau]^k+\rho^{{-}1}[\eta_0(s,0)]^k, \\ \eta_0(s,\tau)=\eta_0(s,0)+\rho\tau. \end{array}\right\} \end{equation}

We show an example of the characteristics and characteristic projections in figure 8.

Figure 8. Characteristic curves given by (6.9) for parameters listed in table 1. Dark blue lines represent curves originating from the initial seepage zone, and light green lines represent curves originating from the propagating front of the seepage zone.

Once the solution is determined, a key quantity of interest is the surface water height at $x = 0$, as it determines the overland flow reaching the river. We denote this critical point along the characteristics as $(T, \eta _0) = (T^*, \eta ^*)$. By setting $x(s,\tau )=0$ in the characteristic equations (6.9) and eliminating $\tau$ from the second equation, we obtain

(6.10)\begin{equation} \begin{pmatrix} T^* \\ \eta^* \end{pmatrix}=\begin{pmatrix} T(s,0)+\dfrac{1}{\rho}(\eta^*-\eta_0(s,0)) \\ (\rho x(s,0)+(\eta_0(s,0))^k)^{1/k} \end{pmatrix}. \end{equation}

We need to consider two cases separately: characteristics starting from the initial seepage zone and characteristics emerging from the propagating seepage front. Each case is given by the different initial conditions, as specified in (6.8b).

In the first case, by substituting the first initial condition from (6.8b) into (6.10), we obtain

(6.11)\begin{equation} \begin{pmatrix} T^* \\ \eta^* \end{pmatrix}=\begin{pmatrix} \dfrac{1}{\rho}(\eta^*-\rho_0^{1/k}(a-s)^{1/k}) \\ (\rho s+\rho_0(a-s))^{1/k} \end{pmatrix}. \end{equation}

By finding $s$ from (6.12a) and substituting into (6.12b) we can express $T^*$ as a function of $\eta ^*$:

(6.12)\begin{equation} T^*(\eta^*)=\frac{1}{\rho}\left[\eta^*-\left(\frac{\rho_0}{\rho-\rho_0}\right)^{1/k}(\rho a_0 - (\eta^*)^k)^{1/k}\right]. \end{equation}

This equation is satisfied for $\eta ^*\in [(\rho _0 a_0)^{1/k}, (\rho a_0)^{1/k}]$. The lower limit corresponds to the initial height, and the upper limit corresponds to the height reached by the characteristic curve starting at $x_0=a_0$. At the upper limit, the characteristic curve reaches the river ($x=0$) at what we refer to as the critical time:

(6.13)\begin{equation} T_{\mathit{crit}}=\frac{1}{\rho}(\rho a_0)^{1/k}. \end{equation}

This critical saturation event is associated with the critical characteristic curve highlighted in figure 8.

For those characteristics starting from the propagating seepage front, $x = a(t)$, we substitute the second initial condition from (6.8b) into (6.10):

(6.14)\begin{equation} \begin{pmatrix} T^* \\ \eta^* \end{pmatrix}=\begin{pmatrix} \mu^{1/k}\mathcal{T}(s)+\dfrac{1}{\rho}\eta^* \\ (\rho s)^{1/k} \end{pmatrix}. \end{equation}

By eliminating $s$, we can express $T^*$ as a function of $\eta ^*$:

(6.15)\begin{equation} T^*(\eta^*) =\mu^{1/k} \mathcal{T}\left(\frac{1}{\rho}(\eta^*)^k\right)+\frac{1}{\rho}\eta^*. \end{equation}

By combining (6.12) and (6.15), we can find the height of the surface water, at the river, $x = 0$, value for all times $T \geq 0$. This is done by solving the implicit equation:

(6.16) \begin{equation} T^*(\eta^*) =\begin{cases} \displaystyle\frac{\eta^*}{\rho}-\frac{1}{\rho}\left(\frac{\rho_0}{\rho-\rho_0}\right)^{1/k}(\rho a_0 - (\eta^*)^k)^{1/k}, & \text{for } \eta^* \leq (\rho a_0)^{1/k},\\ \displaystyle\frac{\eta^*}{\rho}+\mu^{1/k} \mathcal{T}\left(\frac{1}{\rho}(\eta^*)^k\right), & \text{for } \eta^* > (\rho a_0)^{1/k}. \end{cases} \end{equation}

Alternatively, we can express the height of the surface water $\eta ^*$ in terms of the overland component of river inflow, which is represented by the last term in (2.9), $Q_s^*=(\eta ^*)^{{k}}$. This leads to the equation:

(6.17) \begin{equation} t^*(Q_s^*) =\begin{cases} \displaystyle\frac{1}{\rho}(Q_s^*)^{1/k}-\frac{1}{\rho}\left(\frac{\rho_0}{\rho-\rho_0}\right)^{1/k}(\rho a_0 - Q_s^*)^{1/k}, & \text{for } Q_s^* \leq \rho a_0,\\ \displaystyle\frac{1}{\rho}(Q_s^*)^{1/k}+\mu^{1/k} \mathcal{T}\left(\frac{Q_s^*}{\rho}\right), & \text{for } Q_s^* > \rho a_0. \end{cases} \end{equation}

Equation (6.17) represents one of the major results of this work, since it provides an implicit expression for the shape of the hydrograph $Q_s^*(t^*)$.

6.2.3. Approximating the hydrograph in an explicit form

We can obtain an approximated explicit form for the $Q_s^*(t)$ function for $Q_s^* > \rho a_0$. In the limit $\mu \to \infty$ (equivalent to ${Pe} \to \infty$), we may expand around the value of $Q_s^*$ at $t_{\mathit {crit}}$ in (6.13), and write

(6.18a)\begin{equation} Q_s^*(t^*) \sim \rho a(\mu^{{-}1/k}(t^*-t_{\mathit{crit}}))\quad\text{for } t^* \geq t_{\mathit{crit}}, \end{equation}

where we have used $a(t)=\mathcal {T}^{-1}(t)$ from (6.2) to describe the propagation of the wetting front in time. Following approximation (5.7) and (C7), it can be written explicitly as

(6.18b)\begin{equation} a(t)=\underbrace{a_0}_{\mathit{term} 1}+\underbrace{s(t)}_{\mathit{term} 2}+\underbrace{\sigma[1+W_0(-\mathrm{e}^{{-}1-s(t)/\sigma})]}_{\mathit{term} 3}, \end{equation}

where

(6.18c)\begin{equation} s(t)=A t^{{1}/{(n + 1)}} \quad\text{with}\ A = \frac{1}{\rho_0} \left[\frac{n + 1}{m} \frac{\rho-\rho_0}{\theta_s-\theta_r}\left(1-\frac{r_0}{K_s}\right)^{{-}n} \alpha^{{-}n}\right]^{{1}/{(n + 1)}}. \end{equation}

Here, $W_0({\cdot })$ is the Lambert $W$ function, and $\alpha$, $\theta _s$, $\theta _r$, $n$ and $m$ are soil properties used in the Mualem–van Genuchten model (see § C.2).

The first term in (6.18b) corresponds to the flow over the initial seepage zone. The second and third terms represents the growth of flow after reaching the critical point. The third term, in the case of $\sigma \ll 1$, quickly grows from $0$ asymptotically reaching $\sigma$ as $t\rightarrow \infty$, whereas the second term is responsible for further growth of the river flow. Therefore, for thin hillslopes ($\sigma \ll 1$), the growth of river flow after passing the critical point scales proportionally to the $\rho A$ factor.

In the next section, we summarise all approximations derived in this section, and validate them by comparing with each other and full numerical solutions of 1-D and 2-D benchmark models.

7. A numerical comparison between different approximations

7.1. Numerical set-up

In Part 2 (Morawiecki & Trinh 2024b), we performed a numerical verification of the assumption of reducing the 3-D benchmark model to a 2-D model. We also conducted a detailed sensitivity analysis, highlighting the dependencies of model parameters on the resultant peak flows. Here, we continue this analysis by comparing the hydrographs and peak flows between the five different approximations derived and discussed in this paper. The approximations are summarised in table 2 in order from the most complex to the simplest.

Table 2. Summary of the approximations developed in this work.

Similar to the methodology presented in Parts 1 and 2, we assess the performance of the above models in two ways. First, we compare the hydrograph obtained using each model for standard values of parameters characterising UK catchments, as listed in table 3. Second, we run a sensitivity analysis by varying seven model parameters, one at a time, while keeping the others at their default values, and measuring the peak flow in the river after an intensive rainfall. In both numerical experiments, we consider a uniform rainfall over a duration of 24 h.

Table 3. Default values and ranges of parameters used to perform the sensitivity analysis. The part on the right presents parameters not varied during the sensitivity analysis.

7.2. Comparing the hydrographs

A comparison of the hydrographs obtained under the different approximations is presented in figure 9. First, we note that the 1-D models formulated in this paper produce similar hydrographs to the 2-D model from the previous part of our work. However, the 1-D models slightly underestimate the flow, and the solutions are not smooth around the critical point separating early-time and last-time growth.

Figure 9. Hydrograph computed using approximations listed in table 2 for default values of parameters given in table 3. The graph area around the critical point is magnified. Numerical instability are observed for the 1-D model, caused by the finite discretisation of space, which does not allow capturing the exact location of the seepage, and by instabilities related to the governing equation for the seepage zone (2.8a), characterised by a very small diffusive term.

Second, all approximated solutions of the 1-D model produce consistent results for $t\leq t _{\mathit {crit}}$ (except for the explicit solution, which is valid only for $t>t_{\mathit {crit}}$). The results are also similar for $t>t_{\mathit {crit}}$, but inaccuracies related to different approximations start to become noticeable. For example, the implicit solution closely follows the 1-D model solution, but with the flow slightly shifted towards the higher values. This deviation is caused by neglecting the boundary layer characterising the groundwater shape around the critical point. As discussed in Appendix F, this inaccuracy decreases as $\rho$ increases. Replacing numerical solutions for $f(x)$ and $H_0(x)$ with their analytical approximations seems to have a negligible effect on the model for typical sizes of catchment parameters.

Similarly, using approximation (6.18) for the explicit solution leads to the underestimation of the groundwater rise rate, which slows down the growth of flow for higher $t$ values. In addition, the flow around $t=t_{\mathit {crit}}$ is slightly overestimated as a result of neglecting the variation of the $(Q_s^*)^{1/k}$ term appearing in the implicit solution (6.17). Despite these small inaccuracies, the explicit solution still seems to produce excellent qualitative and quantitative agreement. Moreover, due to its simple form, the explicit form allows us to directly understand the effect of various catchment properties on the expected peak flows.

7.3. Sensitivity analysis

We chose seven physical parameters for the sensitivity analysis: catchment width $L_x$, aquifer depth $L_z$, elevation gradient along the hillslope $S_x$, hydraulic conductivity $K_s$, precipitations rates $r$ and $r_0$ and Manning's constant $n_s$. We varied each parameter within the range of its typical values presented in , while keeping the other parameters constant. In each case, we simulated the model's response to a 24-hour-long rainfall event (as shown in figure 9), and then measured the peak river inflow $Q(x=0,t)$ reached at the end of this period. The results of the sensitivity analysis are presented in figure 10.

Figure 10. Sensitivity analysis results showing the dependence between the peak flow reached after 24 h of intensive rainfall and seven different model parameters. The predictions conducted using four models of varying complexity are presented. The dashed region represents the parameter range, for which there is no initial seepage zone ($\rho _0<1$).

We note that the dimensionless parameter determining the existence of the initial seepage zone is given by $\rho _0 = r L_x/(K_s L_z S_x)$. Therefore, as the dimensional parameters are varied, the initial condition may not involve an initial seepage zone if $\rho _0 < 1$. The parameter ranges for which that happens are marked by a dashed region in figure 10. We observe that as long as there is an initial seepage zone, the analytical approximations of the hydrograph are largely accurate over the range of tested parameters.

As expected, the (total) peak flows reached are higher than the maximum levels set by the critical saturation flow, given by (8.2), since the latter only describes the flow reached during the early-time phase. Nevertheless, for most parameter values, the critical flow curve provides a good approximation of the peak flows. In Part 2 of our work, based purely on the 3-D and 2-D simulations, we arrived at the same conclusion.

The cases where the critical flow value highly underestimates the peak flow are around $\rho _0=1$. In these cases, either the seepage zone does not initially exist but the soil is almost fully saturated near the river ($\rho _0$ slightly lower than 1), or it exists but is very small ($\rho _0$ slightly larger). In both cases, rainfall causes the seepage zone to grow significantly relative to its initial size; however, this growth is not captured by the time-independent estimate (8.2).

8. Summary of the key hydrograph features

In figure 9 in the previous section, we showed a typical shape of a hydrograph given by the 1-D model. Here, we summarise the main features of the hydrograph and its importance in benchmarking.

As discussed before, two different phases are visible: phase 1, an early-time fast rise caused by water accumulating over an initial seepage zone; and Phase 2, a late-time slow rise caused by a growing seepage zone. The analytical approximation presented in § 6.2.3 shows that the growth of the overland flow in the second phase can be approximated as $Q_s(t)\approx \rho a(t)$, i.e. it corresponds to the total precipitation rate over the seepage zone $a(t)$ slowly growing in time $t$. Together with the groundwater flow $Q_g=1$, they give a total river inflow $Q(t)=1+\rho a(t)$, which in dimensional units is

(8.1)\begin{equation} Q(t) = K_s S_x L_z + r L_x a(t) = r_0 L_x (1 - a_0) + r L_x a(t), \end{equation}

where $a_0=1-{K_sS_xL_z}/{rL_x}$ is the size of the initial seepage zone.

Following § 6.2, the transition from the first to the second phase corresponds to the moment when the characteristic curve starting from the furthermost point of the initial seepage zone reaches the river ($x=0$). This observation allowed us to estimate the dimensionless critical flow, which in dimensional units correspond to (8.1) with $a(t)=a_0$:

(8.2)\begin{equation} Q_{\mathit{crit}} = \underbrace{K_s S_x L_z}_{\mathit{groundwater\ flow}} + \underbrace{r L_x \left(1 - \frac{K_s S_x L_z}{r_0 L_x}\right)}_{\mathit{overland\ flow}}. \end{equation}

It is reached at the critical time $t_{\mathit {crit}}$ given by (6.13), which in dimensional units is

(8.3)\begin{equation} t_{\mathit{crit}} = \frac{L_z}{r}\left[\frac{S_x^{1/2}K_s n}{L_z^{k-1}} \left(\frac{L_x r}{K_s S_x L_z} - \frac{r}{r_0}\right)\right]^{1/k}. \end{equation}

Following the above event, further growth (8.1) is slow, which is a result of the difference of $\mu ^{1/k}\approx 10^4$ factor between the characteristic timescale of overland flow (responsible for Phase 1) and groundwater flow (responsible for Phase 2). Therefore, $Q_{\mathit {crit}}$ may be a good approximation of the flow even long after the critical time.

We highlight a few additional features of our 1-D benchmark model.

1. (i) When the groundwater component of $Q_{\mathit {crit}}$ (8.2) is much smaller than the overland component (e.g. during intensive rainfalls), the critical flow reached during extreme rainfalls can be approximated by

   (8.4)\begin{equation} Q_{\mathit{crit}} \approx r L_x \left(1 - \frac{K_s S_x L_z}{r_0 L_x}\right). \end{equation}

   Since the consecutive river flow rise is slow, the above estimate can be used as an approximation of the peak flow reached, assuming that the rainfall is long enough to reach the critical time, $t_{\mathit {crit}}$.

2. (ii) Under this approximation, the critical point can be represented as a function of three parameters: rainfall intensity, catchment area (since the flow scales proportionally to both the hillslope width $L_x$ and catchment length $L_y$), and the $K_s S_x L_z/(r_0 L_x)$ factor, which is equal to the fraction of groundwater flow $K_s S_x L_z$ to the mean total flow $r_0 L_x$. This last parameter can be related to what is often referred to as the BFI (Gustard, Bullock & Dixon 1992, Section 3.1.2).

3. (iii) There are some similarities between this expression and other models used in hydrology. Equation (8.4) is a special case of the so-called rational method, which assumes that river flow is proportional to area and precipitation rate (see Bedient et al. 2008, chapter 1.5). The proportionality constant (runoff coefficient) here is identified as $1-\textrm {BFI}$.

   The BFI appears in many statistical methods used in flood estimation, which, unlike our physically-based approach, are based on applying statistical methods such as linear regression to the available catchment data. A notable example is the Flood Estimation Handbook (FEH) flood estimation method by Kjeldsen et al. (2008). It assumes that the median of the annual maximum flow (QMED) scales as $\textrm {QMED}\propto 0.0460^{{BFIHOST}^2}$, where BFIHOST is a soil-based BFI estimator. Note that, similarly to (8.4), the predicted flow decreases with the BFI, but in a nonlinear fashion. Interestingly, other catchment descriptors used in the FEH method include the catchment's area and precipitation, which, like BFI, are also related to the maximum annual flow through nonlinear functions, selected to fit the available data.

Even though our 1-D benchmark model is based on a series of simplifying assumptions (e.g. a thin porous layer and pre-existing seepage zone), which are often not satisfied in real-world catchments, its predictions seem to be reasonable in comparison with data-based catchment models. The connection between our simple scaling laws, shown above and derived analytically from a physical model, and statistical models such as the aforementioned FEH method, which are formulated in a completely different fashion, is intriguing. These results will be presented in a forthcoming work by the present authors, and can also be found in Morawiecki (2023).

9. Conclusions

The primary aim of our work has been to develop and analyse a rigorous benchmark scenario for coupled surface–subsurface flows in a typical catchment. We have achieved this goal by first characterising the typical parameter scales according to the available data on UK catchments (Part 1), formulating and computing the 3-D model and its reduction (Part 2) and finally applying methods in asymptotic analysis to a reduced model valid for catchments dominated by overland dynamics (Part 3).

In this last work, our analysis yields valuable scaling laws for the peak flows (see § 8), which precisely quantify the separation of timescales observed in the hydrographs following an intense period of rain (see a distinct early- and late-time behaviour in figure 9). In particular, we find that the early-time behaviour is governed by rainfall accumulation over a pre-existing seepage zone, followed by a slower flow rise in late time. This latter stage is limited by the speed with which the rising groundwater increases the size of the seepage zone.

All approximations are in good agreement with hydrographs produced by the more complete 1-D and 2-D models, and allow accurate prediction of river peak flows over a wide range of catchment parameters (as long as the underlying assumptions are satisfied). However, different regimes not captured by our model could be studied, including for example behaviour of catchments with no initial seepage ($\rho _0<1$), and late-time catchment behaviour in case of long continuous rainfalls ($t={O}(1)$).

10. Discussion

Our investigations in these three parts have been limited to fairly elementary scenarios and geometries. However, our final results involving the derivation of analytical/asymptotic scaling laws with a clear underlying structure may serve as a valuable benchmark for other hillslope or catchment models. Currently, we observed that the benchmarking of coupled surface–subsurface catchment models has been limited to quantitative numerical comparisons, either with real-world observations or with the numerical output of other schemes (e.g. Maxwell et al. 2014). Although such studies allow practitioners to evaluate the given model's performance in specific conditions, they do not necessarily allow one to draw general conclusions about each model's limit of applicability. As a result, we have no guarantee that a given model will still perform well if applied in situations not captured in the training or validation data as demonstrated in many studies (e.g. Klemeš 1986; Beven 2019).

10.1. Applications of the benchmark to model intercomparisons

Deriving asymptotic estimates for the hydrograph, as done in this work, opens another possibility: models can be compared at a more fundamental level. For instance, we can check if the flows produced by various models scale with the different catchment properties in the expected fashion (and under the conditions we have specified). Detecting situations or limits where the predictions diverge can allow us to better understand the limitations of different models and shed light on how these models can be extended beyond their current limit of applicability. This idea is explored in our two forthcoming works (Morawiecki & Trinh 2023a,b).

As a particular example, statistical models used for flood estimation (such as the FEH method, briefly discussed at the end of § 8) require the selection of empirical catchment descriptors for use in regression formulae to predict flood response (Kjeldsen et al. 2008). In a forthcoming paper (Morawiecki & Trinh 2023b,c), we show that the expression (8.4) can be used to derive a simple expression for predicting peak monthly and annual river flows. This prediction turns out to be highly accurate when applied to real-world data. We then use this result to discuss the limitations of the existing statistical model from the FEH.

Similarly, we can compare the analytical solutions developed within the physical benchmark model with the flow hydrographs generated via conceptual rainfall-runoff models (see the model overview by Peel & McMahon 2020). Although these models use continuous-time rainfall data to generate hydrographs, they are typically not based on the same fluid dynamical models of surface and subsurface flows. Consequently, they can be characterised by different scaling laws than the ones found in this paper. In our forthcoming work (Morawiecki & Trinh 2023a), we demonstrate this discrepancy by studying two models. The first is a probability-distributed model (PDM) used by, e.g. the Environment Agency National Flood Forecasting System (Moore 2007). The second is the grid-to-grid model by Bell et al. (2007), used by the UK Centre for Ecology and Hydrology to provide real-time flow predictions in the UK. In both cases, our simple benchmark scenario allows the identification of key differences between physical and aforementioned conceptual models.

10.2. Extensions and generalisations of the benchmark scenario

Another important line of inquiry is the generalisation of the analysis we have presented to situations that are more representative. Such extensions can involve the analysis of non-uniform rainfall, varying initial conditions to study the response to extended periods of drought, or catchments response to sudden drawdown or outlet water level. These asymptotic regimes have already been studied using the Boussinesq approximation (see, e.g. Parlange et al. 2001; Mizumura 2002). However, our coupled surface–subsurface approach could allow us to better understand the potential role of the seepage dynamics. Then, analytical approximations of the drying process could be used to assess the assumptions of the conceptual rainfall-runoff models, especially since dry periods are sometimes used for a preliminary parameter calibration (Lamb 1999). Finally, we highlight the importance of multiporosity regions in hydrological modelling; these lead to effects such as a preferential flow (Beven & Germann 2013). We provide further details on potential extensions of this study in Morawiecki (2023, chapter 9).

Lastly, it would be interesting to obtain experimental validation of the studied regimes. There have been quite a few lab- and field-scale experimental studies, in which a constant rainfall was artificially generated. However, many of them (e.g. Pauwels & Uijlenhoet 2018) focus on systems limited to the groundwater flow only. There are some experimental studies in which rainfall we observed to yield seepage growth, e.g. Abdul & Gillham (1989), Kollet et al. (2017) and Scudeler et al. (2017). However, in these experiments, the soil was initially dry, i.e. there was no initial seepage, so they correspond only to the $\rho _0=0$ case. It would be interesting to conduct controlled experiments with more realistic settings, i.e. with an already developed groundwater table, and compare the resultant hydrograph with our model predictions.