1. Introduction
In this third and final part of our work, we leverage the parametric study (Part 1: Morawiecki & Trinh Reference Morawiecki and Trinh2024a) and two-dimensional (2-D) benchmark models (Part 2: Morawiecki & Trinh Reference Morawiecki and Trinh2024b) 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:
\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 Reference Grayson, Moore and McMahon1992). 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. Reference Moore, Bell, Cole and Jones2007). 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. Reference Parkin, O'donnell, Ewen, Bathurst, O'Connell and Lavabre1996; Bathurst et al. Reference Bathurst, Ewen, Parkin, O'Connell and Cooper2004; Beven Reference Beven2018).
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. (Reference Sulis, Meyerhoff, Paniconi, Maxwell, Putti and Kollet2010) and Maxwell et al. (Reference Maxwell2014), 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 Reference Cunnane1987; Kjeldsen, Jones & Bayliss Reference Kjeldsen, Jones and Bayliss2008) the comparison of asymptotic scaling laws between different catchment models is not a commonly used approach in the modern hydrology (although see Vieira (Reference Vieira1983) for a comparison of Saint Venant approximations, and Cook, Knight & Wooding (Reference Cook, Knight and Wooding2009) 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. Reference Shaw, Beven, Chappell and Lamb2010). 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 Reference Hálek and Švec2011, chap. 2; Troch et al. Reference Troch2013), 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 Reference Polibarinova-Kochina and Wiest1962; Troch et al. Reference Troch2013). 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 Reference Warrick, Lomen and Islas1990).
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 (Reference Parlange, Rose and Sander1981) and Tao, Wang & Lin (Reference Tao, Wang and Lin2018). These analytical solutions and their approximations are important as they provide benchmark models for testing numerical schemes (e.g. benchmarks by MacDonald et al. (Reference MacDonald, Baines, Nichols and Samuels1995) for the overland flow and by Tracy (Reference Tracy2006) for the subsurface flow), and can be used to develop less computationally demanding modelling approaches such as TOPMODEL by Kirkby & Beven (Reference Kirkby and Beven1979). 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. Reference Maxwell2014).
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 (Reference Dupuit1863) and Forchheimer (Reference Forchheimer1914), which states that the groundwater flow is predominantly horizontal, and that the total flow scales proportionally with the saturated aquifer thickness. Boussinesq (Reference Boussinesq1877) 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 Reference Bartlett and Porporato2018), or the Dupuit–Boussinesq equation (see, e.g. Guérin, Devauchelle & Lajeunesse Reference Guérin, Devauchelle and Lajeunesse2014). 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. (Reference Paniconi, Troch, Van Loon and Hilberts2003) and Cook et al. (Reference Cook, Knight and Wooding2009).
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 (Reference Wooding and Chapman1966), Anderson & Brooks (Reference Anderson and Brooks1996), Troch, Paniconi & Emiel van Loon (Reference Troch, Paniconi and Emiel van Loon2003) and Bartlett & Porporato (Reference Bartlett and Porporato2018). 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:
\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 Reference Henderson and Wooding1964; Woolhiser & Liggett Reference Woolhiser and Liggett1967), 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 Reference Morawiecki and Trinh2024a). 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 Reference Kirkby2019) is not considered. Under this assumption, the heights 
$H_g$ and 
$h_s$ can be combined to form a single dependent variable,
\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
\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. Reference Troch2013) for an unconfined aquifer:
\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. Reference Troch, Paniconi and Emiel van Loon2003), 
$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 Reference Caputo and Stepanyants2008). 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
\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:
\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:
\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$:
\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)
\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,
\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.:
\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:
\begin{equation} Q(1, t) = 0, \end{equation}
where the dimensionless total flow, 
$Q$, is defined as
\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$:
\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
\begin{gather} \sigma =\frac{L_z}{L_xS_x} = \frac{\textrm{thickness of the porous layer}}{\textrm{elevation drop along the hillslope}}, \end{gather}
\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 (Reference Horton1936).
In the second case of (2.8), we have introduced an additional dimensionless parameter:
\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:
\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 Reference Morawiecki and Trinh2024a).
3. Numerical methodology and typical dynamics
Recall that the solutions are essentially characterised by the overland and groundwater heights, and a quadruplet of parameters:
\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:
\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 Reference Morawiecki2022). 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$),
\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.
(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.(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.(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.(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
\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:
\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
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
\begin{gather} \partial_x H(0, t) = 0, \quad \partial_x H(1, t) ={-}\sigma^{{-}1}, \end{gather}
\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
\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
\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.
(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$.(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$.(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
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]$,
\begin{gather} \frac{{Pe}^{{-}1}}{a}\frac{\partial g}{\partial X}+g^k=\rho_0(1-aX) - 1, \end{gather}
\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
\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
\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
\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)$:
\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,
\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$.
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$:
\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,
\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
\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
\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
\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
\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:
\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
\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
\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
\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
\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.
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
\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
\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 ^*$:
\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:
\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):
\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 ^*$:
\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:
\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:
\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
\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
\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
\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 Reference Morawiecki and Trinh2024b), 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
\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$:
\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
\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.
(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}}$.(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 Reference Gustard, Bullock and Dixon1992, Section 3.1.2).(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. Reference Bedient, Huber and Vieux2008, 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. (Reference Kjeldsen, Jones and Bayliss2008). 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 (Reference Morawiecki2023).
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. Reference Maxwell2014). 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š Reference Klemeš1986; Beven Reference Beven2019).
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 Reference Morawiecki and Trinh2023a,Reference Morawiecki and Trinhb).
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. Reference Kjeldsen, Jones and Bayliss2008). In a forthcoming paper (Morawiecki & Trinh Reference Morawiecki and Trinh2023b,Reference Morawiecki and Trinhc), 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 Reference Peel and McMahon2020). 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 Reference Morawiecki and Trinh2023a), 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 Reference Moore2007). The second is the grid-to-grid model by Bell et al. (Reference Bell, Kay, Jones and Moore2007), 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. Reference Parlange, Stagnitti, Heilig, Szilagyi, Parlange, Steenhuis, Hogarth, Barry and Li2001; Mizumura Reference Mizumura2002). 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 Reference Lamb1999). Finally, we highlight the importance of multiporosity regions in hydrological modelling; these lead to effects such as a preferential flow (Beven & Germann Reference Beven and Germann2013). We provide further details on potential extensions of this study in Morawiecki (Reference Morawiecki2023, 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 Reference Pauwels and Uijlenhoet2018) 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 (Reference Abdul and Gillham1989), Kollet et al. (Reference Kollet2017) and Scudeler et al. (Reference Scudeler, Paniconi, Pasetto and Putti2017). 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.
Acknowledgements
We thank Sean Longfield (Environmental Agency) for many useful interactions and for motivating this work via the 7th Integrative Think Tank hosted by the Statistical and Applied Mathematics CDT at Bath (SAMBa). We also thank Thomas Kjeldsen (Bath), Tristan Pryer (Bath) and Rob Lamb (Lancaster/JBA Trust) for insightful discussions. We are indebted to the reviewers and the JFM editorial team for their comments and suggestions that were instrumental in the final development of this paper.
Funding
P.M. is supported by a scholarship from the EPSRC Centre for Doctoral Training in Statistical Applied Mathematics at Bath (SAMBa), under the project EP/S022945/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. List of symbols
For ease of reference, we provide a list of symbols in table 4.
Table 4. List of symbols
Appendix B. Relation between the 1-D Boussinesq equation and the 2-D Richards equation
In this appendix, we provide additional details on the derivation of the 1-D model presented in this paper, in connection with the physical models presented in Part 2. We shall explain how the governing equations (2.3) and (2.5) relate to the 2-D Richards equation given in (5.9a) from Part 2 (Morawiecki & Trinh Reference Morawiecki and Trinh2024b). Some parts of the following are classical, and relate to the derivation of the Boussinesq equation (cf. Bear & Verruijt Reference Bear and Verruijt1987 for details); our new contribution is to consider the influence of the overland flow in the seepage zone on the Boussinesq equation and to couple this latter equation with its standard formulation for the remaining part of the domain.
In this study, we assume that hillslope flow is predominantly 2-D in the 
$xz$ cross-section; this reduction from three dimensions to two dimensions is discussed in detail in Part 2. In addition, we consider the small aspect ratio limit, 
$L_z\ll L_x$ and, hence, 
$\beta _{zx} =L_z/L_x \to 0$. In the groundwater region, following (5.9a) from Part 2, the leading-order flow then satisfies
\begin{equation} \left.\frac{\mathrm{d}\theta}{\mathrm{d} h}\right\vert_{h=h_g'}\frac{\partial h_g'}{\partial t} = \frac{\partial }{\partial \hat{z}} \left[K_r(h_g') \left(\frac{\partial h_g'}{\partial \hat{z}}+1\right)\right], \quad 0 < \hat{z} < 1, \end{equation}
with 
$\hat {z} = 0$ corresponding to the bottom of the aquifer and 
$\hat {z} = 1$ corresponding to the top surface. Note that for 
$\hat {z} < H$, the Richards equation becomes much simpler, since for saturated soil we have 
${\mathrm {d}\theta }/{\mathrm {d} h}=0$ and 
$K_r(h)=1$. Solving (B1) and imposing the no-flow boundary condition at the bottom of the aquifer 
$\hat {z}=0$, we obtain
\begin{equation} h_g'(\hat x,\hat z,t)=H'(\hat x,t)-\hat z, \quad 0 < \hat{z} < H'(\hat{x}, t), \end{equation}
which corresponds to a hydrostatic vertical profile of pressure. Note that based on the above solution, for regions of completely saturated soil, the 2-D function 
$h_g'(\hat x,\hat z,t)$ can be replaced by the 1-D indicator, 
$H'(\hat x,t)$. The curve 
$\hat {z} = H'(\hat {x}, t) \in (0, 1]$ corresponds to the groundwater table, which separates between saturated (where 
$h_g' > 0$) and unsaturated (where 
$h_g' < 0$) regions, if both coexist.
We assume that the system is configured as shown in figure 11. Thus, it is divided into a region B (fully saturated case), where 
$0 \leq \hat {x} \leq a(t)$ and the ground is entirely saturated, with 
$H = 1$. Similarly, we have region A (unsaturated case) for 
$\hat {x} > a(t)$, where 
$H < 1$ and there is an unsaturated column where 
$H < \hat {z} < 1$.
Figure 11. Control volume (A) outside the seepage zone and (B) inside the seepage zone. One-way arrows represent flow in and out of the control volumes.
For ease of interpretation, we shall return, for the next few subsections, to dimensional quantities (related to the governing equation (2.7)). In addition, it is more convenient to use Richards equation expressed in the Cartesian coordinates 
$(x, z)$. In this coordinate system, (B2) becomes
\begin{equation} h_g(x,z,t) = S_x x + H(x,t) - z, \quad \text{for } x\in [0,L_x],\ z\in [0,L_z], \end{equation}
where 
$L_z\hat {z}=z-S_x x$, 
$H=L_zH'$ and 
$h_g=L_zh_g'$. Following Darcy's law, we have 
$\boldsymbol {q}=K_s \boldsymbol {\nabla } (h_g+z)$. After substituting (B3), we obtain the leading term of horizontal flow:
\begin{equation} q_x={-}K_s\frac{\partial }{\partial x} (h_g+z)={-}K_s\left(S_x+\frac{\partial H}{\partial x}\right). \end{equation}
Note that the above point flux is independent of 
$z$. Multiplying by 
$H$, we thus see that the total (unsigned) flow along the hillslope is then given by
\begin{equation} Q_g = H |q_x| = K_s H\frac{\partial H}{\partial x} + K_s H S_x. \end{equation}Next, we consider the unsaturated part of the hillslope (region A) and the seepage zone (region B) separately since the latter requires adding an overland flow component.
B.1. Region A: governing equation outside the seepage zone, $H<1$
Let us consider region A, where there is a layer of unsaturated soil. Recall that the pressure 
$h_g$ is zero at the free surface of the groundwater, where 
$z=S_xx+H_g(x,t)$. This allows us to set the integration constant 
$H$ appearing in (B3) to 
$H(x,t)=H_g(x,t)$.
We consider the control volume, 
$V$, within the column horizontally bounded by 
$[x, x + \Delta x]$, as shown in figure 11. By conservation and the divergence theorem, it is argued that the total flux integral around the boundary, 
$\partial V$, is zero and, hence,
\begin{align} 0 &= \oint_{\partial V} \boldsymbol{\nabla} (h+z) \boldsymbol{\cdot} \boldsymbol{n}\,\mathrm{d}l \nonumber\\ &={-} H(x)q_x(x) + H(x+\Delta x)q_x(x+\Delta x) + \left.\int_x^{x+\Delta x} q_z\right|_{z=S_x x + H(x)} \,\mathrm{d}\kern 0.06em x. \end{align}
Richards equation in the fully saturated groundwater region reduces to 
$\boldsymbol {\nabla }\boldsymbol{\cdot}\boldsymbol {q}=K_s\nabla ^2(h+z)=0$ (as discussed following (B1)). Thus, in the 
$\Delta x\rightarrow 0$ limit, we can approximate the vertical flux at the top of the saturated zone as
\begin{equation} q_z^{top}=q_z|_{z=S_x x + H(x)}={-}\frac{\partial }{\partial x}(H q_x)=K_s\frac{\partial }{\partial x}\left(HS_x+H\frac{\partial H}{\partial x}\right), \end{equation}
where we have used the expression for 
$q_x$ derived in (B4).
Now, let us consider the unsaturated zone above the groundwater table (dashed region in figure 11). The inflow from the top boundary is 
$r(t)$, and the outflow to groundwater is 
$q_z^{top}$. If during a rainfall the inflow is greater than the outflow, the total volume of water in the soil column per surface area,
\begin{equation} \mathcal{V}=\int_0^{L_z}\theta(h(z,t))\,\mathrm{d} z, \end{equation}
increases, eventually leading to a rise of the groundwater table 
$H$ over time.
In order to find the exact rate of change of 
$H(x,t)$, one should solve (B1) and find 
$H(x,t)$ for which 
$h_g(x,z=H(x,t),t)=0$. Using the chain rule, we have
\begin{equation} \frac{\partial }{\partial t}h_g(z=H(x,t),t)=\left.\frac{\partial h_g(x,z,t)}{\partial z}\right|_{z=H(x,t)}\frac{\partial H(x,t)}{\partial t}+\left.\frac{\partial h_g}{\partial t}\right|_{z=H(x,t)}=0, \end{equation}from which the groundwater growth rate is
\begin{equation} \frac{\partial H(x,t)}{\partial t}={-}\left.\left(\frac{\partial h_g(x,z,t)}{\partial z}\right)^{{-}1}\frac{\partial h_g}{\partial t}\right|_{z=H(x,t)}. \end{equation}However, rather than solving Richards equation in the unsaturated zone, a standard time-dependent Boussinesq-based approach to groundwater modelling is based on introducing a drainable porosity defined as
\begin{equation} f \equiv \dfrac{{\rm d}\mathcal{V}}{{\rm d}H}. \end{equation}
In essence, the drainable porosity describes, for a given change in the height 
$H$, the corresponding change in the subsurface water volume 
$\mathcal {V}$ (B8). It is often assumed that the drainable porosity 
$f$ is a constant parameter characterising a given soil. However, as we argue in Appendix C, this assumption is not consistent with the Richards-based approach, in which the growth of the volume 
$\mathcal {V}$ does not have to immediately lead to a rise in the groundwater volume. Therefore, in the same appendix, we introduce a mean drainable porosity 
$f(x)$, which does not fully represent the dynamics governed by the Richards equation but still allows us to exactly reproduce the time when the groundwater reaches the surface in the Richards-based approach.
Since the volume of subsurface water changes as a result of precipitation 
$r(x, t)$ and inflow/outflow 
$q_z^{top}$, we have
\begin{equation} f(x)\frac{\partial H}{\partial t}=\dfrac{{\rm d}\mathcal{V}}{{\rm d}t}=q_z^{top}+r(x,t). \end{equation}
Substituting (B7) into the above expression gives us the governing equation for 
$H$:
\begin{equation} f(x)\frac{\partial H}{\partial t}=K_s\frac{\partial }{\partial x}\left(HS_x+H\frac{\partial H}{\partial x}\right)+r(x,t). \end{equation}
This equation is equivalent to (2.3) under the additional assumptions that infiltration is equal to precipitation: as in Part 2, we ignore the effects of evapotranspiration and, in addition, we assume that no overland flow is generated unless the soil becomes fully saturated (i.e. rainfall never exceeds soil infiltration capacity, 
$r(x,t)\leq K_s$).
B.2. Region B: governing equation for seepage zone $H\geq 1$
In the case of the seepage zone, the groundwater height is fixed at 
$H_g=L_z$. However, we need to determine the additional surface water height, 
$h_s(x,t)$, which is given by solving the Saint Venant equation (see (3.3) from Part 2):
\begin{equation} \frac{\partial h_s}{\partial t}=\frac{\partial }{\partial x}[Q_s(h_s)] + r(x, t) - I(x, t), \end{equation}
where the surface flow, 
$Q_s(h_s)$, is given by the Manning's equation (2.6), and 
$I$ is the infiltration rate, equal to the negated vertical flow, 
$-q_z$, at ground level.
To determine the vertical flow, 
$q_z$, we perform a similar conservation argument as in Appendix B.1. First, at the interface between the subsurface and surface flows, we set continuity of pressure and flow. The first condition allows us to specify the pressure head, 
$h_g$, at the surface:
\begin{equation} h_g(x, z, t) = h_s(x, t) \quad \text{at } z = S_x x + L_z. \end{equation}
Consequently, we set the integration constant 
$H$, appearing in (B3), to 
$H(x,t)=L_z+h_s(x,t)$. Substituting it into (B4), we obtain
\begin{equation} q_x={-}K_s\left(S_x+\frac{\partial h_s}{\partial x}\right) \quad \text{for }0 < z < L_z. \end{equation} Now, let us consider the control volume of groundwater contained in the vertical column 
$[x, x + \Delta x]$, as illustrated in figure 11 in region B. In this case, we have 
$H_g=L_z$, and the mass balance equation (B7) yields
\begin{equation} q_z|_{z=S_x x + L_z}={-}L_z\frac{\partial q_x}{\partial x} = \frac{\partial }{\partial x}\left(K_sL_zS_x+K_sL_z\frac{\partial h_s}{\partial x}\right) = \frac{\partial }{\partial x}\left(K_sL_z\frac{\partial h_s}{\partial x}\right), \end{equation}
after substituting the expression (B16) for the groundwater flow, 
$q_x$. Thus, we have the required expression for the infiltration rate:
\begin{equation} I(x, t) ={-}q_z|_{z=S_x x + L_z} ={-}\frac{\partial }{\partial x}\left(K_sL_z\frac{\partial h_s}{\partial x}\right). \end{equation}Note that the infiltration is positive if surface water infiltrates into the soil or negative if groundwater emerges to the surface.
By substituting 
$Q_s$ from (2.6) and infiltration 
$I$ from (B18) into (B14), we obtain
\begin{equation} \frac{\partial h_s}{\partial t}=\frac{\partial }{\partial x}\left(K_sL_z\frac{\partial h_s}{\partial x}+\frac{\sqrt{S_x}}{n_s} h_s^k\right)+r(x,t). \end{equation}
Note that in the seepage zone, 
$H=H_g+h_s$, where 
$H_g=L_z$ is constant. Hence, this equation is equivalent to the governing equation (2.5). Thus, we have derived both cases of the governing equation (2.7), forming the basis of this work.
Appendix C. Mean drainable porosity function, $f(x)$
C.1. On the mean drainable porosity
In Appendix B, we defined the drainable porosity 
$f$ as the volumetric change in the subsurface volume for a given change in the groundwater height, 
$\textrm {d}{\mathcal {V}}/\textrm {d}{H}$. However, this approach is not fully consistent with the solution for the Richards equation for subsurface flow.
Typically, in scenarios where the precipitation increases to a constant value 
$r>1$, a characteristic wetting front is observed. This front can be seen in figure 8(d) in Part 2, and its propagation is discussed in detail by Caputo & Stepanyants (Reference Caputo and Stepanyants2008). The front moves downward towards the groundwater table, changing soil saturation and eventually leading to the rise of the groundwater table. This behaviour is not captured by a standard time-dependent Boussinesq equation or our 1-D model.
However, our analysis in §§ 6.1 and 6.2 shows that the most important mechanism by which groundwater contributes to the peak flow generation is by extending the seepage zone. Moreover, in Appendix F, we show that horizontal groundwater flow is negligibly slow and does not affect the speed at which the groundwater is rising over a short timescale characterising a typical rainfall. This means that the groundwater becomes saturated when rainwater fills the available drainable volume 
$v_h$, i.e. after time 
$t=v_h/r$. The same time is predicted by the 1-D model by setting an 
$x$-dependent mean drainable porosity 
$f(x)$, defined as
\begin{equation} f(x)=\frac{v_H(x)}{D(x)}, \end{equation}
where 
$D(x)$ is the thickness of the unsaturated soil layer. The above choice of mean drainable porosity is the one that we use throughout this work.
To summarise, even though setting a time-independent porosity and constant groundwater recharge does not capture the delay required for the wetting front to reach the groundwater, it allows for the correct prediction of the soil critical time and the resultant peak flows observed in this work. In this model, note that 
$H(x,t)$ does not represent the exact thickness of the saturated zone that forms groundwater but corresponds to the amount of water absorbed by the soil, even if it has not yet reached the saturated zone. However, this difference in interpretation does not seem to affect the main results obtained in this paper. Formulating this complex system in terms of a 1-D partial differential equation facilitates the analysis and computation.
C.2. Computation of mean drainable porosity
Here, we formally derive an expression for the mean porosity given the 1-D model parameters for the considered benchmark scenario. First, following the Richards equation, we find the pressure head 
$h_g$ profile above the groundwater table and use it to evaluate the drainable volume for a column of soil of height 
$D(x)$ above the groundwater table.
We note that in the considered case of 
$\beta _{zx} = L_z/L_x \ll 1$, the leading solution of the 2-D Richards equation involves only vertical flow along the 
$\hat {z}$ axis, as given by (B1). In the scenario considered in this paper, we assume that the system is initially in a steady state for a constant rainfall 
$r_0$. Under these assumptions, we integrate the time-independent version of (B1) to form a first-order nonlinear ODE for the pressure head 
$h_g(\hat z)$:
\begin{equation} K_r(h_g) \left(\dfrac{{\rm d}h_g}{{\rm d}\hat{z}}+1\right) = \frac{r_0}{K_s}, \end{equation}
where the constant of integration on the right-hand side has been chosen to match the dimensionless infiltration 
$r_0/K_s$, passing through the surface. Let us consider a column of the soil above a groundwater table, and we take 
$h_g|_{\hat {z}=0}=0$. Given the pressure head, we can find the saturation 
$\theta (h_g(\hat z))$ following the Mualem–van Genuchten model:
\begin{equation} \theta(h_g) =\begin{cases} \theta_r+\dfrac{\theta_s-\theta_r}{(1+(\alpha h_g)^n)^m} & h_g<0, \\ \theta_s & h_g\ge 0, \end{cases} \end{equation}
where 
$\theta _s$ and 
$\theta _r$ are the saturated and residual water content, whereas 
$\alpha$, 
$n$ and 
$m=1-{1}/{n}$ are other Mualem–van Genuchten parameters characterising a given soil. We can use (C3) to compute the drainable volume and resulting mean drainable porosity (C1):
\begin{equation} f(x) = \frac{v_H(x)}{D(x)} = \frac{1}{D(x)} \int_{0}^{D(x)} [\theta_s-\theta(h_g(\hat z))] \,{\rm d}\hat z. \end{equation}
Solving (C2) and then numerically integrating (C4) allows us to calculate 
$f$.
C.3. Analytical approximations of mean drainable porosity
As an alternative to the integral expression above, let us develop an approximation for 
$v_H$, based on the assumption that the groundwater table is located near the land surface. Since for 
$h\rightarrow 0$, 
$K_r(h)\rightarrow 1$, the leading-order solution for (C2) around 
$\hat z=0$ satisfies
\begin{equation} h_g(\hat z) \sim \left(\frac{r_0}{K_s} - 1\right) \hat z. \end{equation}Then, integrating (C4) with the Mualem–van Genuchten model (C3) gives
\begin{equation} f_1(x) = (\theta_s-\theta_r)\left[{}_2F_1\left(m,\frac{1}{n};1+\frac{1}{n};\left({-}a\left(\frac{r_0}{K_s} - 1\right)D(x)\right)^n\right)-1\right], \end{equation}
where 
${}_2F_1$ is a hypergeometric function. One can also find the leading-order approximation of (C6) for 
$D \ll 1$, which yields
\begin{equation} f_2(x) = \frac{m}{n+1}(\theta_s-\theta_r)\left[-\left(\frac{r_0}{K_s} - 1\right) \alpha D(x)\right]^n. \end{equation}However, the approximation (C6) is more accurate. Functions (C6) and (C7) are compared with the full numerical solution in figure 12.
Figure 12. (a) Initial vertical profile of the pressure head 
$h_g$ and corresponding saturation 
$\theta$ in a column of soil. (b) Dependence of mean porosity 
$f={v_h}/{D}$ on the depth of the groundwater below the surface. As shown, approximations (C5)–(C7) accurately describe soil properties close to the groundwater table. We used Mualem–van Genuchten parameter values from the previous part of our work, i.e. 
$\alpha =3.367$ m
$^{-1}$, 
$\theta _s=0.388$, 
$\theta _R=0.115$ and 
$n=1.282$.
Appendix D. Derivation of steady-state overland flow for $\mu \to \infty$
The steady state for the seepage zone is given by (5.1b):
\begin{equation} 1+\sigma\frac{\partial h_0}{\partial x}+\mu h_0^k=\rho_0(1-x), \end{equation}
where we used 
$h_0(x)=H(0)-1$ to represent the surface water height.
As we have noted previously, our regime of interest is where 
$\mu$ is large. It is convenient to rescale via 
$x = aX$, where 
$a$ was defined as the size of the seepage zone, so that the domain is fixed in 
$X\in [0, 1]$. In the limit 
$\mu \to \infty$, we anticipate from dominant balance that the overland flow is small, 
$h_0 \to 0$, and hence we rescale 
$h_0 = \mu ^{-1/k}g(X)$ in order to balance the third term with the right-hand side of (D1). With the boundary conditions (4.3a,c), this rescaling gives
\begin{gather} \frac{{Pe}^{{-}1}}{a}\frac{\partial g}{\partial X}+g^k=\rho_0(1-aX) - 1, \end{gather}
\begin{gather}g'(0) = 0 \quad \text{and} \quad g(1) = 0. \end{gather}In (D2a), for convenience, we have defined the small parameter
\begin{equation} {Pe}^{{-}1} = \frac{\sigma}{\mu^{1/k}} \ll 1, \end{equation}
which, according to table 1, is approximately 
${Pe}^{-1} \approx 10^{-5}$. As we show, the asymptotic analysis as 
${Pe}^{-1} \to 0$ involves the analysis of an outer solution, where 
$0 \leq X < 1$, and an inner solution with a boundary layer, where 
$X \to 1$.
In the outer region, we expand the solution and contact line as
\begin{equation} g(X) = g_0(X) + {Pe}^{{-}1} g_1(X) + \cdots \quad \text{and} \quad a = a_0 + {Pe}^{-\gamma} a_1 + \cdots, \end{equation}
where 
$\gamma > 0$ and is to be determined later. Then we develop the following leading-order approximation of the outer surface water height via (D2a),
\begin{equation} g(X) \sim g_0(X) = [\rho_0(1-a_0 X)-1]^{1/k}. \end{equation}
The above outer expansion exhibits an infinite gradient, possibly before the contact line is reached, at the point 
$X = (1 - \rho _0^{-1})/a_0$. However, we show in the following that 
$a_0$ is such that this point corresponds to 
$X = 1$.
In order to develop the solution within the boundary layer, we consider rescaling
\begin{equation} g(X) = {Pe}^{-\beta} \hat{G}(\hat{X}) \quad \text{and} \quad X = 1 - {Pe}^{-\alpha} \hat{X}, \end{equation}
with 
$\alpha$, 
$\beta$ and 
$\gamma$ from (D4) to be determined. Applying these transformations to (D2a) gives
\begin{align} &{Pe}^{-(1+\beta-\alpha)} \left[-\frac{1}{a_0}\frac{\partial \hat{G}}{\partial \hat{X}}\right] +{Pe}^{{-}k\beta} [ \hat{G}^k] \nonumber\\ &\quad =[\rho_0(1-a_0) - 1] + {Pe}^{-\alpha} [\rho_0 a_0 \hat{X}] - {Pe}^{-\gamma}[\rho_0 a_1] + o({Pe}^{-\alpha}). \end{align}
The balance of 
${Pe}$ terms is achieved for 
$1+\beta -\alpha =k\beta =\alpha =\gamma$, which corresponds to 
$\alpha = \gamma = k/(2k-1)$ and 
$\beta =1/(2k-1)$. Equation (D7) for the leading 
${O}(1)$ term is 
$\rho _0(1-a_0) - 1 = 0$, and hence this approximation gives us an estimate for the leading-order contact line position,
\begin{equation} a\sim a_0=1-\frac{1}{\rho_0}. \end{equation}
Indeed, this equation confirms that the leading-order outer solution, (D5), predicts an infinite gradient as 
$X \to 1$. In order to derive an approximation for the gradient of the surface water height near the contact line, we must proceed to the next order in the inner region.
Equation (D7) at 
${O}({Pe}^{-{k}/{(2k-1)}})$ yields
\begin{equation} -\frac{1}{a_0}\frac{\partial \hat{G}}{\partial \hat{X}}+\hat{G}^k= \rho_0 a_0 \hat{X} - \rho_0 a_1, \end{equation}along with the two boundary conditions:
\begin{equation} \hat{G}(0) = 0 \quad \text{and} \quad \hat{G}(\hat{X}\rightarrow\infty) \sim (\rho_0-1)^{1/k} \hat{X}^{1/k}. \end{equation}
The second boundary condition above corresponds to the inner limit of the outer solution (D5). The above first-order ODE, together with the two boundary conditions, forms a boundary value problem with eigenvalue 
$a_1$. It cannot be solved analytically. In any case, if desired, the above problem can be computed numerically, and it would yield the correction to the contact line, 
$a_1 = a_1(\rho _0)$.
Note finally that the gradient of 
$h_0$ at 
$x=a$ can now be estimated directly from (D1). Since 
$h_0(a)=0$, we have, using the expansion (D4),
\begin{equation} \left.\frac{\partial h_0}{\partial x}\right\vert_{x=a} = \frac{\rho_0(1-a) - 1}{\sigma} \sim \frac{\rho_0(1-a_0) - 1 - \rho_0 a_1{Pe} ^{-\gamma}}{\sigma} ={-} \frac{\rho_0 a_1}{\sigma}{Pe} ^{-\gamma}. \end{equation} For Manning's law, 
$k= 5/3$, so 
$\gamma = k/(2k-1) = 5/7$. Hence, we expect the gradient to be 
$\textrm {d}{h_0}/\textrm {d}\kern 0.06em {x} = {O}({Pe}^{-5/7})$ while the contact line position equally satisfies 
$a - a_0 = {O}({Pe}^{-5/7})$. These two scaling laws are confirmed by solving the boundary value problem given by (D2a) numerically and transforming it back to the original variables (see figure 13).
Figure 13. (a) The size of the seepage zone relative to its first-order approximation 
$a-a_0$. (b) The gradient of 
$H$ at this point. The results were obtained by solving (D1). The fitted power law is consistent with the theoretical exponent 
$\gamma =5/7$.
Appendix E. Derivation of the initial groundwater table for $\sigma \rightarrow 0$
Here, we find the leading-order solution for the ODE (5.6) under the limit 
$\epsilon =\sigma \rightarrow 0$:
\begin{equation} \epsilon \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$, and where 
$a(t) \sim a_0$ is the leading-order contact line position as 
${Pe} \to \infty$. Expanding 
$H_0$ in powers of 
$\epsilon$, we obtain the approximation:
\begin{equation} H_0(x)= \rho_0(1-x) +\rho_0 \epsilon +{O}(\epsilon^2). \end{equation}
This solution does not satisfy the boundary condition at 
$x=a_0$. For 
$x\to a_0$, we develop an inner expansion by re-scaling 
$x=a_0+\epsilon X$ and 
$H_0(x) = g(X)$. The ODE now becomes
\begin{equation} \dfrac{{\rm d}g}{{\rm d}\kern 0.06em X} = \frac{\rho_0(1-a_0-\epsilon X)}{g(X)} - 1, \quad \text{for } X>0, \end{equation}
with 
$g(0)=1$. Solving now to the first two orders, we have
\begin{equation} g(X)=1+\epsilon \rho_0(1-{\rm e}^{{-}X}-X) +{O}(\epsilon^2), \end{equation}
where we have used 
$a_0 = 1-1/\rho _0$ from (5.5).
The composite solution is obtained by adding the outer and inner asymptotic expansions from (E2) and (E4), respectively, and subtracting their common part. The final result, to two orders of accuracy, is the following approximation:
\begin{equation} H_0(x)=\rho_0(1-x+\epsilon-\epsilon \,{\rm e}^{-{(x-a_0)}/{\epsilon}}) + {O}(\epsilon^2). \end{equation}
After replacing 
$\epsilon$ with 
$\sigma$, we obtain the solution (5.7) stated in the main text.
Appendix F. Derivation of the time-dependent groundwater solution for $t \to 0$ and $\rho \to \infty$
We provide additional details for the early-time analysis of § 6.1, particularly in connection with the boundary-layer asymptotics. Consider (4.3b) outside the seepage zone (
$H<1$):
\begin{equation} f(x) \frac{\partial H}{\partial t} = \frac{\partial }{\partial x}\left(\sigma H\frac{\partial H}{\partial x} + H\right) + \rho. \end{equation}
As noted in §6.1, we consider the initial condition, 
$H(x, t = 0) = H_0(x; \, \rho _0)$, to be the steady-state response of the system to a precipitation rate, 
$\rho _0$. That is, the initial groundwater solution satisfies (5.6):
\begin{equation} (\sigma H_0 H_0' + H_0)' + \rho_0 = 0, \end{equation}subject to the boundary conditions (4.3b,d) in § 5, i.e.
\begin{equation} H(x=a(t))=1 \quad\text{and}\quad \left.\left(\sigma H\frac{\partial H}{\partial x} + H\right) \right\vert_{x=1} = 0. \end{equation}
Let us consider the short-time behaviour of the time-dependent equation. Let 
$t=\epsilon t'$, 
$H(x,t)=H_0(x) + \epsilon H'(x,t)$, where 
$\epsilon \ll 1$. Then, (F1) becomes
\begin{equation} f(x) \frac{\partial H'}{\partial t'} = \rho - \rho_0 + \epsilon \frac{\partial }{\partial x}\left(\sigma H_0\frac{\partial H'}{\partial x} + \sigma H'\frac{\partial H_0}{\partial x} + H'\right) + {O}(\epsilon^2). \end{equation}Therefore, the first two leading terms of the small-time groundwater solution are
\begin{equation} H(x,t) = H_0(x; \, \rho_0) + \left[\frac{\Delta \rho}{f(x)} \right]t + {O}(t^2), \quad \text{as } t \to 0 \text{ with } x > a(t), \end{equation}
where 
$\Delta \rho \equiv \rho -\rho _0$. Thus, we see that (F5) predicts that the groundwater height increases linearly with time by an amount proportional to the sudden impulse of rain (or, rather, its difference 
$\Delta \rho$).
However, the above asymptotic solution assumes that 
$x - a(t) = {O}(1)$, and indeed it fails to account for the fact that 
$H = 1$ when 
$x = a(t)$. We must consider it to be an outer solution, valid away from 
$x = a(t)$. A comparison between the outer approximation (F5) and the full solution is shown in figure 7(a). In the case of the fully coupled surface–subsurface flow, the solution of the PDE diverges from the outer solution (F5) within a small boundary layer around the seepage front.
The size of the boundary layer tends to zero as 
$\rho \to \infty$. This is tested as follows. First, we compute the full PDE model (2.8a) and find the place where 
$H(x_{numeric})=1$. Then, we estimate the boundary layer thickness by finding the difference between 
$x_{numeric}$ and 
$x_{approx}=1-1/\rho$, which is the leading-order (outer) approximation of 
$a(t)$. As figure 14 demonstrates, it depends both on time 
$t$ and precipitation represented by 
$\rho$.
The above is confirmed via a dominant balance. Let 
$x=a+\delta x'$. Then (F4) becomes
\begin{equation} f(x) \frac{\partial H'}{\partial t'} = \rho - \rho_0 + \frac{\epsilon}{\delta^2} \frac{\partial }{\partial x'}\left(\sigma H_0\frac{\partial H'}{\partial x'} + \sigma H'\frac{\partial H_0}{\partial x'} + \delta H'\right) + {O}(\epsilon^2). \end{equation}The diffusion terms become significant and balance the precipitation 
$\rho$ when 
$\epsilon /\delta ^2 \sim \Delta \rho$. Then the thickness of the boundary layer is on the order of 
$\delta =\epsilon ^{1/2}\Delta \rho ^{-1/2}$, i.e. it increases proportionally to 
$t^{1/2}$ and 
$\Delta \rho ^{-1/2}$. These trends are confirmed in figure 14.
Appendix G. Derivation of the explicit solution for $Q(t)$
In this appendix, we show that the implicit solution (6.16) for late time (i.e. for 
$Q_s^* > \rho a_0$),
\begin{equation} t^*(Q_s^*) = \frac{1}{\rho}(Q_s^*)^{1/k}+\mu^{1/k} \mathcal{T}\left(\frac{Q_s^*}{\rho}\right), \end{equation}
can be written in an explicit form given by (6.18a–c), as long as the increase of flow after the critical time is negligibly small compared with the critical flow value 
$Q_{\mathit {crit}}=\rho a_0$. This assumption is motivated by the observation that the dynamics of the flow increase after reaching critical flow slows down significantly.
Let us consider the flow 
$Q_s^*$ near its critical value, 
$Q_s^* = \rho a_0 + \epsilon Q'$, and consider the asymptotic limit 
$\epsilon \rightarrow 0$. Then (G1) becomes
\begin{equation} t^* = \frac{1}{\rho}(\rho a_0 + \epsilon Q')^{1/k}+\mu^{1/k} \mathcal{T}(Q_s^*). \end{equation}Therefore, we have
\begin{equation} t^* = \underbrace{\frac{1}{\rho}(\rho a_0)^{1/k}}_{t_{\mathit{crit}}}+\mu^{1/k} \mathcal{T}(Q_s^*)+{O}(\epsilon). \end{equation}
By neglecting the 
${O}(\epsilon )$ term, we get
\begin{equation} t^* - t_{\mathit{crit}} = \mu^{1/k} \mathcal{T}(Q_s^*), \end{equation}
which allows us to obtain an explicit equation for 
$Q_s^*$:
\begin{equation} Q_s^* = \mathcal{T}^{{-}1}(\mu^{{-}1/k}(t^*-t_{\mathit{crit}})). \end{equation} Now we need to find an inverse function for 
$t=\mathcal {T}(a(t))$, i.e. a function that provides the location of the seepage front for a given time. Function 
$\mathcal {T}$ is defined by (6.2) as
\begin{equation} \mathcal{T}(x=a(t)) \sim \frac{f(x=a(t))}{\rho - \rho_0} (1 - H_0(x=a(t))). \end{equation}
We use the approximation of 
$f(x)$ given by (C7), and the approximation of 
$H(x)$ given by (5.7),
\begin{gather} f(x)=CD(x)^n,\quad \text{where } C=\frac{m}{n+1}(\theta_s-\theta_r)\left[-\left(\frac{r_0}{K_s} - 1\right) \alpha\right]^n, \end{gather}
\begin{gather}D(x)=1-H(x)=1-\rho_0(1-x+\sigma-\sigma \,{\rm e}^{-{(x-a_0)}/{\sigma}}). \end{gather}By substituting (G7a) to (G6) we obtain
\begin{equation} \mathcal{T}(x=a(t)) = \frac{1}{\rho - \rho_0} C D(x=a(t))^{n+1}. \end{equation}
Now, we solve this equation for 
$D(x=a(t))$:
\begin{equation} [C^{{-}1}(\rho - \rho_0)\mathcal{T}(x=a(t))]^{{1}/{(n+1)}} = D(x=a(t)). \end{equation}After substituting (G7b) we get
\begin{equation} \rho s(t) = A\rho t^{{1}/{(n+1)}} = 1-\rho_0(1-x+\sigma-\sigma\, {\rm e}^{-{(x-a_0)}/{\sigma}}), \end{equation}
where we introduced 
$A=[C^{-1}(\rho - \rho _0)]^{{1}/{(n+1)}}$ and 
$t=\mathcal {T}(x=a(t))$ to shorten the notation. The solution of this equation for 
$x$ is
\begin{equation} x = a(t) = s(t)+\sigma+\underbrace{1-\frac{1}{\rho_0}}_{a_0}\,+\,\sigma W_0(-{\rm e}^{{-}1-s(t)/\sigma}), \end{equation}
where 
$W_0({\cdot })$ is the Lambert 
$W$ function. Therefore, the solution can be written as
\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}This concludes the derivation of the explicit solution (6.18). The accuracy of this approximation is demonstrated in figure 9.
Open access
