On the development and analysis of coupled surface–subsurface models of catchments. Part 2. A three-dimensional benchmark model and its properties

Abstract

The objective of this three-part work is to formulate and rigorously analyse a number of reduced mathematical models that are nevertheless capable of describing the hydrology at the scale of a river basin (i.e. catchment). Coupled surface and subsurface flows are considered. In this second part, we construct a benchmark catchment scenario and investigate the effects of parameters within their typical ranges. Previous research on coupled surface–subsurface models have focused on numerical simulations of site-specific catchments. Here, our focus is broad, emphasising the study of general solutions to the mathematical models, and their dependencies on dimensionless parameters. This study provides a foundation based on the examination of a geometrically simple three-dimensional benchmark scenario. We develop a non-dimensional coupled surface–subsurface model and extract the key dimensionless parameters. Asymptotic methods demonstrate under what conditions the model can be reduced to a two-dimensional form, where the principal groundwater and overland flows occur in the hillslope direction. Numerical solutions provide guidance on the validity of such reductions, and demonstrate the parametric dependencies corresponding to a strong rainfall event.

1. Introduction

Since the publication of the Stanford Watershed Model by Crawford & Linsley (1966), a wide range of computational models of catchment-scale hydrology have been developed (Singh & Frevert 2003). Indeed, over two hundred models have been identified in the extensive review by Peel & McMahon (2020).

Such computational models are primarily designed in order to predict the evolution of surface and subsurface flow in a particular river basin given the input precipitation via rainfall or snowfall. These so-called rainfall-runoff models are often divided into three classes: empirical, conceptual and physical (Sitterson et al. 2018); this last category of physical models involves those that are developed from the known physical principles of hydrodynamics. For instance, the Richards equation is commonly used to model the subsurface flow through the saturated or unsaturated soil, while the Saint Venant equation is used to model the overland and the channel flow. For a detailed introduction, see Shaw et al. (2010). Such governing equations form the foundation of many currently used computational integrated catchment models, e.g. MIKE SHE (Abbott et al. 1986a,b), HydroGeoSphere (Brunner & Simmons 2012), ParFlow (Kollet & Maxwell 2006) and OpenGeoSys (Kolditz et al. 2012).

However, in contrast to computational studies, there seems to have been more limited work on the systematic mathematical analysis of the fundamental principles of coupled surface–subsurface catchment-scale models. A proper mathematical formulation can allow us to better understand the importance of parameters, establish the limits of simplifications used in computational models and develop analytical or semi-analytical solutions in certain scenarios.

1.1. On the development and benchmarking of computational models

The Stanford Watershed Model IV is a conceptual model, which is considered to be amongst the earliest attempts to computationally model the entire hydrological cycle. Its publication resulted in the subsequent development of an enormous number of independent computational models (Donigian & Imhoff 2006). However, further computational power was needed before the first physically based models were implemented. Notable early examples include TOPMODEL (Kirkby & Beven 1979), MIKE SHE (Abbott et al. 1986a,b) and IHDM (Institute of Hydrology Distributed Model, cf. Beven, Calver & Morris 1987).

The abundance of independent catchment models results in a need to better understand their accuracy and differences. Within the industry, such models are typically assessed by comparing model predictions (usually after earlier calibration) with available data, such as river flow or groundwater depth measurements (see a detailed introduction to rainfall-runoff modelling by Beven 2011). However, there is criticism, e.g. by Hutton et al. (2016), that the models in hydrology are often not reproducible. Beven (2018, p. 6) highlighted some fundamental issues that continue to exist in the state-of-the-art of catchment modelling. He noticed that:

Where model intercomparisons have been done, different models give different results, and it is often the case that the rankings of models in terms of performance will vary with the period of data used, site or type of application. This would seem to be a very unsatisfactory situation for the advancement of the science, especially when we expect that when true predictions are made, they will turn out to be at best highly uncertain and at worst quite wrong.

In response to this problem, many numerical methodologies for calibration, cross-validation and uncertainty estimation have been developed (see e.g. Beven & Binley 1992; Gupta, Beven & Wagener 2006). These methods allow us to assess, in a more unbiased way, the accuracy of the models. However, they do not necessarily point out the reason for potential inaccuracies. As Kirchner (2006) argued, advancing the science of hydrology requires developing not only models that match the available data, but models that are theoretically justified.

Independently, there has been an effort to develop simple (idealised) catchment geometries that can be used as benchmarks to assess the accuracy of integrated catchment models in fully controlled conditions. Kollet & Maxwell (2006) used a tilted V-shaped catchment geometry (figure 1a) to compare predictions for overland flow given by four different hydrological catchment models with an analytical one-dimensional solution. Then, they introduced a simple two-dimensional hillslope (figure 1b), which they used to explore the sensitivity of an integrated ParFlow model for geometry settings (e.g. water table depth, hydraulic conductivity and soil heterogeneities). The same benchmark scenarios were used by Sulis et al. (2010) to compare ParFlow and CATHY models (Bixio et al. 2000). This study was followed by far more extensive intercomparison studies by Maxwell et al. (2014) and Kollet et al. (2017), which used these and other benchmark scenarios to compare the results obtained using a wide range of integrated catchment models.

Figure 1. Illustration of the idealised catchment geometries developed in the works of Kollet & Maxwell (2006) and Gilbert et al. (2016). Geometries (a) and (c) represent a tilted V-shaped river valley with two hillslopes and a river in the middle, with the latter geometry introducing subsurface flow in the third dimension. Geometry (b) represents a single two-dimensional hillslope with a river channel located at the right boundary. (a) Two-dimensional tilted V-shaped catchment, (b) hillslope cross-section and (c) three-dimensional tilted V-shaped catchment.

In the meanwhile, simple catchment/hillslope scenarios have also been used to assess coupled surface and subsurface flow with other models – this includes examination of evapotranspiration (Kollet et al. 2009), atmosphere (Sulis et al. 2017), biochemistry (Cui, Welty & Maxwell 2014), the impact of climate change (Markovich, Maxwell & Fogg 2016) and the effects of different types of heterogeneities, e.g. the heterogeneity of the land surface (Rihani, Chow & Maxwell 2015), soil properties (Meyerhoff & Maxwell 2011) and even flow through fractures (Sweetenham, Maxwell & Santi 2017). The two studies by Jefferson et al. (2015) and Gilbert et al. (2016) introduced a three-dimensional tilted V-shaped catchment with a constant soil depth (figure 1c). The authors used this geometry to perform a sensitivity analysis of integrated catchment models – the first study by Jefferson et al. (2015) focused on the energy flux terms, while the second by Gilbert et al. (2016) studied the heterogeneity of soil permeability. In both studies, the sensitivity analysis results were used to obtain a certain level of dimensionality reduction by applying the active subspace method (Constantine 2015).

An open question remains, however, whether one can simplify the model and its parameter space based on the analysis of the governing equations (even in a simplified catchment scenario), rather than based on the numerical results; this could provide more rigorous insight into the limits of applicability of the above computational reductions.

Another aspect we shall investigate in this work concerns the study of key non-dimensional parameters characterising surface–subsurface hydrological processes. We highlight some prior works that have used non-dimensionalisation in order to analyse governing equations describing individual flow components: for example, this has been applied by Akan (1985) in the Saint Venant equations to study the water infiltration into the ground. It has also been used by e.g. Warrick, Lomen & Islas (1990), Warrick & Hussen (1993) and Haverkamp et al. (1998) for the study of the one-dimensional Richards equation, describing water vertical infiltration through the unsaturated soil.

A notable work, in which non-dimensionalisation plays a prominent role for the case of coupled surface–subsurface models, was performed by Sivapalan, Beven & Wood (1987), and focuses on the TOPMODEL scheme of Kirkby & Beven (1979). A similar study was performed by Calver & Wood (1991) for the IHDM model (Beven et al. 1987). In particular, Calver & Wood (1991) define a list of ten dimensionless parameters, study the dependencies between selected parameters and discuss the properties of the hydrographs. However, the relevant scale of dimensionless parameters is not assessed in this latter work.

1.2. On the development of a simple benchmark model

The modern-day catchment hydrology is studied based on the simulation of complex integrated catchment models. So far, however, the authors have not found many comprehensive studies on the design and analysis of simple benchmark scenarios for coupled surface–subsurface catchment models. Our work in Part 1 (Morawiecki & Trinh 2024) has initiated this task via a thorough examination of the typical parameter sizes. In this part, we focus on the design of a three-dimensional benchmark, study its typical dynamics and discuss its reduction to lower-dimensional models.

Compared with the existing literature, there are three novel elements in our study:

1. (i) Our benchmark scenario is posed on a simple geometry, but the surface/subsurface governing equations are posed in a general three-dimensional dimensionless form.

2. (ii) We use the dimensionless model to provide a rigorous argument behind the simplifications commonly used in computational hydrology. We discuss the reduction of a problem geometry to two dimensions in detail, and comment on the kinematic/dynamic wave approximation. We achieve this by setting clear conditions on the size of dimensionless parameters, and justify them based on the typical values of model parameters obtained in the previous part of our work (see table 1 from Part 1).

3. (iii) We use the benchmark model to numerically explore the impact of the remaining parameters on the system in response to intensive rainfall. Because we attempt to do this in a systematic and analytical way, this work also serves to set a more rigorous benchmark standard for future studies. For example, scaling laws are derived that may serve as a benchmark for other model schemes.

Note that our study is restricted to modelling the formation of storm flow during an intensive rainfall (Guérin et al. 2019); however, similar benchmark scenarios can be considered in order to study other flow regimes. This may include, for instance, drought flow observed during a period without any rainfall (Brutsaert & Nieber 1977), or a sudden drawdown drainage when a rapid change of water level occurs at the outlet (Sanford, Parlange & Steenhuis 1993).

We start by formulating a three-dimensional benchmark scenario in § 2, which is non-dimensionalised in § 3. In § 5 we show that this model can be reduced to a two-dimensional form by neglecting the subsurface and overland flow component in the $y$-direction. Following the numerical methodology from § 6, this model simplification is numerically assessed in § 7. The impact of each parameter in the resulting two-dimensional model is summarised in § 8, which is followed by the discussion in § 9.

Symbols. There are many symbols in this work. For ease of reference, we provide a list of symbols in tables 2 and 3 in Appendix A.

2. Formulation of a simplified three-dimensional catchment model

In this section, we formulate a simplified catchment model, inspired by the infiltration-excess, saturation-excess and tilted V-shaped catchment scenarios from the benchmark study by Maxwell et al. (2014).

We introduce the following three scenarios, as depicted in figure 2.

1. (a) The V-shaped catchment. This scenario, shown in figure 2(a), represents a V-shaped catchment with a thick aquifer, where subsurface water is transferred both through the soil and through the underlying bedrock. The aquifer dimensions are $L_x\times L_y\times L_z$, where $L_z$ is the thickness of the permeable layer of the aquifer. The elevation gradient along the hillslope is denoted as $S_x$, and along the direction of the river as $S_y$. Similar to the V-shaped scenario studied by Maxwell et al. (2014), we shall assume that the channel has a constant width, $w$, and zero depth, $d=0$. Later in § 5, we demonstrate that under certain conditions, the scenario reduces to largely two-dimensional dynamics along the hillslope.

2. (b) The deep aquifer. This scenario, shown in figure 2(b), represents a two-dimensional hillslope with a thick aquifer, where the subsurface water is transferred through both the soil and the underlying bedrock. Following the infiltration- and saturation-excess scenarios discussed in Maxwell et al. (2014), the channel is assumed to have a rectangular $xz$ cross-section with width, $w$, and depth, $d$.

3. (c) The shallow aquifer. This scenario, shown in figure 2(c), represents a catchment with a low-productive aquifer, in which the subsurface water is transferred only through a thin soil layer. Mathematically, the geometry of the problem is equivalent to the deep aquifer scenario with $L_z\ll L_x$. We analyse this scenario in Part 3.

Figure 2. Simplified catchment geometry in the considered scenarios (not to scale). (a) V-shaped catchment scenario, (b) deep aquifer scenario and (c) shallow aquifer scenario.

The focus of work in this Part 2 is the study of the V-shaped catchment scenario and its reduction to a two-dimensional deep aquifer scenario. In Part 3, we shall demonstrate that under the additional restrictions of the shallow aquifer scenario, further analysis can be performed through a long wavelength reduction. In the V-shaped catchment scenario, an orthogonal coordinate system $(x, y, z)$ is chosen such that $z$ is vertical and $y$ is directed along the channel. Using the reflection symmetry of the catchment, we can describe the catchment behaviour by only considering a hillslopes only on one side of the river.

When formulating the governing equations for overland and subsurface flow, we are going to use a more convenient non-orthogonal coordinate system, where the axes $(\hat {x}, \hat {y}, \hat {z})$ are directed along the hillslope edges. Hence, $\hat {x}$ is directed along the hillslope ($\hat {x}=0$ representing the location of the channel), $\hat {y}$ along the channel ($\hat {y}=0$ representing the location of the outlet) and $\hat {z}$ vertically ($\hat {z}=0$ representing the bottom of the aquifer). After the coordinate transformation, the entire catchment can be represented as a cuboid of dimensions $L_{\hat x} \times L_y \times L_z$. The following coordinate transformation is used:

(2.1a–c)\begin{equation} x = \hat{x} \sqrt{1-\left(\frac{S_y}{S_x}\right)^2}, \quad y = \hat{x} \frac{S_y}{S_x} + \hat{y}, \quad z = S_x \hat{x} + S_y \hat{y} + \hat{z}. \end{equation}

We introduced $L_{\hat x}$ to represent the catchment width along the $\hat x$ direction given as

(2.2)\begin{equation} L_{\hat x}=\frac{L_x}{\sqrt{1-\left(\dfrac{S_y}{S_x}\right)^2}}. \end{equation}

The land surface in this geometry corresponds to

(2.3)\begin{equation} H_{surf}(\hat x,\hat y)=z(\hat{x},\hat{y},\hat{z}=L_z)=S_x \hat{x} + S_y \hat{y} + L_z. \end{equation}

Note that real-world systems are characterised by different levels of heterogeneity of the surface, soil, and parent material properties. Here, in order to construct a minimal model, we consider properties to be homogeneous; this is similar to the assumptions made by Maxwell et al. (2014). Thus, the surface is assumed to have uniform roughness, and the properties of soil and rock layer are assumed to be homogeneous, i.e. have a uniform hydraulic conductivity and water-retention curve. Also, we assume that the soil and bedrock do not include the presence of macropores and fractures, which would lead to the formation of preferential flow – see more in the reviews by Bouma (1981) and Neuzil & Tracy (1981). Because of the last assumption, the model may not properly represent the infiltration through the unsaturated zone in many of the real-world systems. As noted e.g. by Beven & Germann (2013), including these effects in the model may significantly affect the time scale of infiltration.

2.1. Asymptotic limits of geometrical parameters

It is convenient to discuss the asymptotic limits of the key non-dimensional parameters that characterise the geometry. First, we have the slope ratios between the channel and hillslope directions

(2.4)\begin{equation} \epsilon = \frac{S_y}{S_x}, \end{equation}

which for a typical UK catchment is $\epsilon \in [0.13,0.25]$ (the estimates represent the interquartile range based on parameters characterising over $1200$ UK catchments, values of which were estimated in Part 1 (here the first quartile is $0.13$, and the third quartile is $0.25$)). We also have the aspect ratio between the catchment height and the catchment dimension along the river

(2.5)\begin{equation} \beta_{zy} = \frac{L_z}{L_y}, \end{equation}

which for a typical UK catchment is $\beta _{zy}\in [0.0007, 0.025]$. Finally, we have the aspect ratio between the catchment height and the catchment length along the hillslope

(2.6)\begin{equation} \beta_{zx} = \frac{L_z}{L_{\hat{x}}}, \end{equation}

which for a typical UK catchment is $\beta _{zx}\in [0.1,2.1]$. Note that as $\beta _{zy}/\beta _{zx} \to 0$, we get long catchments with a width much shorter than their length, while for $\beta _{zy}/\beta _{zx} \to \infty$, we get short catchments with a width much longer than their length.

The impact of these two parameters on the catchment geometry is schematically presented in figure 3. Here, we draw lines of constant topographic elevation on a projection of the catchment onto $z = 0$. Note that, for example in figure 3(a) for $S_y = 0$, surface and subsurface flow will typically occur in the $x$ direction, perpendicular to the river direction. In contrast, for figure 3(c), we may expect to observe a significant flow component parallel to the river direction.

Figure 3. These illustrations provide a guide to understand the impact of changing values of the slope, $S_y$ (a–c), and length, $L_y$ (d–f) on our V-shaped catchment geometry (the channel is shaded). Lines of constant elevation of the topography are represented with dashed lines, drawn on top of a projection of the catchment onto $z = 0$. By the definition of a catchment, the top and bottom boundaries are perpendicular to lines of constant elevation (since an unperturbed flow will follow lines of the steepest descent). These dashed lines help to visualise the geometry of the later contour plots.

It is important to remember that, since our interest is in the study of the benchmark model, we are not necessarily limited to studying only physical regimes. That is, it is still interesting to study the asymptotic limits so that we can establish the qualitative trends.

2.2. Relationship to Maxwell et al. (2014)

Here, we briefly outline how the scenarios introduced above relate to the scenarios presented in the benchmark analysis of Maxwell et al. (2014).

In §§ 4.1 and 4.2, Maxwell et al. (2014) introduce two scenarios called the infiltration excess and saturation excess, respectively. In the infiltration scenario, precipitation exceeds the saturated soil conductivity ($r>K_s$). Only part of the precipitation infiltrates through the soil, while the remaining part accumulates at the surface to form an overland flow (the so-called Horton overland flow). In the saturation-excess scenario ($r< K_s$), overland flow is not generated unless the entire soil becomes fully saturated.

Both scenarios are posed on a single hillslope, which represents a thin layer of soil ($L_z=5$ m) with a slope following the $x$ direction, while the river is assumed to have a fixed surface water height. Thus, this geometry represents the shallow aquifer scenario shown in figure 2, where the flow takes place only in a thin layer of the soil. Note that this geometry does not include water infiltration to the deeper permeable layers of the parent material (as in the deep aquifer scenario in figure 2), which is an effect that characterises the majority of the real-world aquifers (note a small area of aquifers without the groundwater on the UK map in figure 4 from Part 1).

Figure 4. Boundaries defined for the V-shaped tilted catchment. (a) Boundaries for the subsurface flow and (b) boundaries for the surface flow.

A second limitation of the geometries considered by Maxwell is that there is no slope along the river, which drives the flow down the river valley. Although the authors included the slope perpendicular to the hillslope in a separate scenario introduced in their § 4.3 (V-shaped catchment), this benchmark scenario does not include subsurface modelling; therefore, the water infiltration into the soil was not studied.

Our scenarios in this work combine the above two elements, i.e. groundwater flow through deep aquifers and slope both perpendicular and along the river. Therefore, we consider a V-shaped catchment with an additional $z$-dimension allowing the saturation to vary with depth, as in the hillslope scenario. The need to introduce a tilted coordinate system comes from the fact that the elevation gradient (determining the direction of surface flow) is not perpendicular to the river, since it must have a small component along the $y$-axis. In order to satisfy the no-surface flow boundary condition at the catchment boundary, the bottom and top boundaries of the hillslope are thus inclined by a small angle, $\phi =\mathrm {asin}(S_y/S_x)$, relative to the rectangular domain in the infiltration- and saturation-excess scenarios.

Last but not least, we use the typical catchment parameters as estimated in Part 1; note that these values can be significantly different from those numerical values used in the work of Maxwell et al. (2014). Based on our simulations, we observed that if one were to use the parameter values given by Maxwell et al. (2014), this would lead to unrealistic steady states, where the seepage covers almost the entire catchment (even for relatively low levels of mean precipitation).

3. Governing equations (dimensional)

We begin with the dimensional model. As introduced in § 2 of Part 1, we consider three types of flow: the subsurface flow (the three-dimensional (3-D) Richards equation), the overland flow (the 2-D Saint Venant equations) and the channel flow (1-D Saint Venant equation). In this section, we present governing equations for each of the flow components in our benchmark scenario, together with the corresponding boundary conditions. General reviews of these governing equations can be found in the works of Farthing & Ogden (2017), Schaake (1975) and references therein.

3.1. Three-dimensional Richards equation for the subsurface flow

The subsurface flow $\boldsymbol {q_g}(x,y,z,t)$ depends on the pressure head $h_g(x,y,z,t)$. Its evolution in time $t$ is commonly modelled using a 3-D Richards equation (see e.g. Dogan & Motz 2005; Weill, Mouche & Patin 2009), which is given by

(3.1)\begin{equation} \frac{\mathrm{d}\theta}{\mathrm{d} h_g}\frac{\partial h_g}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q_g}, \quad \text{where } \boldsymbol{q_g} = K_s K_r(h_g)\boldsymbol{\nabla}(h_g+z). \end{equation}

Here, $\boldsymbol {\nabla }=({\partial }/{\partial x},{\partial }/{\partial y}, {\partial }/{\partial z})$ is a standard nabla operator, $K_s > 0$ is the saturated soil conductivity and ${\mathrm {d}\theta (h_g)}/{\mathrm {d} h_g}$ is the so-called specific moisture capacity. We assume that the volumetric water content $\theta (h_g)$, and relative hydraulic conductivity $K_r(h_g)$ are functions of the pressure head given by the Mualem–van Genuchten (MvG) model (Van Genuchten 1980)

(3.2a)$$\begin{gather} \theta(h_g) = \begin{cases} \theta_r+\dfrac{\theta_s-\theta_r}{(1+(\alpha_{MvG}h_g)^n)^m} & h_g<0 \\ \theta_s & h_g\ge0 \end{cases}, \end{gather}$$

(3.2b)$$\begin{gather}K_r(h_g) = \begin{cases} \dfrac{(1-(\alpha_{MvG}h_g)^{n-1}(1 +(\alpha_{MvG}h_g)^n)^{-m})^2}{(1+ (\alpha_{MvG}h_g)^n)^{m/2}} & h_g<0 \\ 1 & h_g\ge0 \end{cases}. \end{gather}$$

Here, the value of $h_g=0$ corresponds to the pressure head at the groundwater table surface, which separates the fully saturated zone ($h_g>0$) from the partially saturated zone ($h_g<0$) (see the later figure 6(a) for a reference image). In essence, the MvG model describes the key hydraulic properties of the soil, hydraulic conductivity and saturation as nonlinear functions of the pressure head $h_g$. The model introduces further parameters $\alpha _{MvG}$, $\theta _r$, $\theta _s$, $n$ and $m=1-{1}/{n}$, which depend on the soil properties. The residual water content $\theta _r$ and saturated water content $\theta _s$ represent the lowest and the highest water content, respectively. The $\alpha _{MvG}$ parameter in $\mathrm {m}^{-1}$ represents the scaling factor for the pressure head $h_g$ (m). The $n$ coefficient describes the pore sizes distribution.

3.2. Two-dimensional Saint Venant equations for the overland flow

If the precipitation exceeds the inflow into the soil, water can accumulate on the surface and form overland flow. Typically, and following e.g. Tayfur & Kavvas (1994) and Liu et al. (2004) this flow is described using the 2-D Saint Venant equations that govern the overland water height, $h_s(x, y, t)$. Following the discussion in § 2.2 of Part 1, we consider mass and momentum conservation. Firstly, the continuity equation is given by

(3.3)\begin{equation} \frac{\partial h_s}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q_s}(h_s)+R_{eff}-I, \end{equation}

where $I = I(x,y,t)$ is the infiltration rate, and $R_{eff} = R(x,y,t) - ET(x,y,t)$ is the effective precipitation rate, which we define as the difference between the precipitation rate, $R$, and the evapotranspiration rate, $ET$.

The flux, $\boldsymbol {q_s}$, that appears in the Saint Venant equation (3.3), is commonly obtained in hydrology using an empirical relationship known as Manning's law. Written in vector form, it is given by

(3.4)\begin{equation} \boldsymbol{q_s} = \frac{1}{n_s}h_s^{5/3}\frac{\boldsymbol{S}_f}{\sqrt{|\boldsymbol{S}_f|}}, \end{equation}

where $n_s$ is an empirically determined value known as Manning's coefficient, and describes the overland surface roughness; $\boldsymbol {S_f}$ is a dimensionless friction slope defined as gradient of energy of water per unit weight.

When Manning's law in (3.4) is substituted into the continuity equation (3.3), this yields a single equation for the two unknowns, $h_s$ and $\boldsymbol {S}_f$. In general, the friction slope, $\boldsymbol {S_f}$, is given by momentum conservation (cf. (2.7) in Part 1). However, in computational integrated catchment models, a kinematic approximation is often used which neglects all effects on $\boldsymbol {S}_f$ other than gravity. This approximation is used in e.g. Parflow (Maxwell et al. 2009), although there are others such as e.g. MIKE SHE that implement a more complete, diffusive approximation (MIKE SHE 2017). In the case of the kinematic approximation

(3.5)\begin{equation} \boldsymbol{S}_f\sim\boldsymbol{S}_0, \end{equation}

where $\boldsymbol {S}_0=-\nabla H_{surf}$ is the elevation gradient. In this paper, we shall adopt the above kinematic approximation. This reduction significantly simplifies the problem since, under this approximation, the overland flows only down the hillslope (the $\hat {x}$-direction). As Daluz Vieira (1983) argues, this approximation may give inaccurate predictions when the system is close to reaching a steady state.

3.3. One-dimensional Saint Venant equation for the channel flow

Finally, we need to formulate the governing equation for the surface flow in a rectangular channel of width $w$. The channel is directed along the $\hat {y}$-axis. Following Daluz Vieira (1983) and Chaudhry (2007), the channel flow is modelled as a 1-D Saint Venant equation that governs the channel water height, $z = h_c(\hat {y}, t)$, and is given by

(3.6)\begin{equation} w\frac{\partial h_c}{\partial t} = q_{in} - \frac{\partial q_c}{\partial \hat{y}}, \end{equation}

where $w(h_c, \hat {x})$ is the channel width (constant in the case of a rectangular channel), and $q_{in}$ is a source term governing the total surface and subsurface inflow into the river. As for the overland equations, the flux, $q_c$, is assumed to be given by the empirical Manning's law, which takes the form

(3.7)\begin{equation} q_c = A\frac{\sqrt{S_f^{river}}}{n_c}\left(\frac{A}{P}\right)^{2/3}, \end{equation}

where $A$ is the channel cross-section, $P$ is the channel wetted perimeter, $n_c$ is Manning's coefficient dependent on banks and channel bed roughness and $S_f^{river}$ is the friction slope, which under kinematic approximation is equal to the elevation gradient along the river

(3.8)\begin{equation} S_f^{river}=S_y. \end{equation}

In summary, the solution of the channel flow involves the substitution of Manning's equation (3.7) and the friction slope (3.8) into the Saint Venant equation (3.6). For the case of the V-shaped catchment illustrated in figure 2, where there is a rectangular channel, this involves setting the area $A = wh_c$ and $P = w + 2 h_c$.

The above channel flow model, when coupled to the hillslope forms a challenging numerical computation due to the nonlinearity. Instead, for the purpose of numerical computation, we apply a model simplification of the channel flow similar to what is considered by Maxwell et al. (2014). In this simplification, we set $A=wh_c$ and approximate $P \approx w$, and hence ignore the friction effects of the channel sidewalls. In this case, Manning's equation (3.7) becomes

(3.9)\begin{equation} q_c = w\frac{\sqrt{S_y}}{n_c}h_c^{5/3}. \end{equation}

The later simulations will thus involve the solution of the Saint Venant equation (3.6) with the shallow Manning equation (3.9) and (3.8). The advantage of the above approximation is that the channel flow problem satisfies a similar partial differential equation to the surface flow, but with adjusted coefficient values.

3.4. Boundary conditions

The domain consists of four types of boundaries: (i) the catchment boundary, $\varGamma _B$, both for the surface and subsurface part of the domain, including the bedrock constraining the aquifer from the bottom; (ii) the land surface $\varGamma _s$; (iii) the river bank, $\varGamma _R$; (iv) the river outlet, $\varGamma _O$; and (v) the river inlet $\varGamma _I$ (see figure 4).

1. (i) Firstly, there is no surface flow through the catchment boundary. Also, for simplicity, we will assume that there is no groundwater flow through this boundary – the rainwater can only leave the catchment via the channel flow. Hence, on $\varGamma _B$, we set no-flow conditions for both subsurface and surface flow

   (3.10a)\begin{equation} \boldsymbol{q_g}\boldsymbol{\cdot}\boldsymbol{n} = 0, \quad \boldsymbol{q_s}\boldsymbol{\cdot}\boldsymbol{n} = 0\quad \text{on $\varGamma_B$}. \end{equation}
   Alternatively, one could introduce a free-flow condition for the groundwater flow, $\boldsymbol {q_g}\boldsymbol{\cdot}\boldsymbol {n} = 0$, to allow for the outflow of the groundwater flow through the catchment boundary. In this work, we have chosen no-flow conditions to guarantee that the entirety of the rainfall eventually reaches the channel, which simplifies the resultant water balance.

2. (ii) Next, on the land surface, $\varGamma _s$, continuity of pressure and flow between the groundwater and surface water yields

   (3.10b)\begin{equation} h_s = \begin{cases} 0 & \text{if } h_g<0\\ h_g & \text{if } h_g>0 \end{cases} \quad \text{and} \quad \boldsymbol{q_g}\boldsymbol{\cdot}\boldsymbol{n} = I \quad \text{on $\varGamma_s$}. \end{equation}
   This first condition imposes continuity of pressure only if the groundwater reaches $\varGamma _s$, while the second imposes the condition of rain infiltration, $I$.

3. (iii) On the river bank, $\varGamma _R$, we also impose continuity of pressure between the channel water, which is characterised by a hydrostatic profile, $h(z)=h_c-z$, and the subsurface pressure head, $h_g$

   (3.10c)\begin{equation} h_g = h_c - z \quad \text{on $\varGamma_R$}. \end{equation}

4. (iv) At the inlet, located at the upstream end of the river, $\varGamma _I$, we can impose an inflow from the upstream part of the catchment, which is located outside of the modelled domain. In general, it can change over time, and so

   (3.10d)\begin{equation} q_c = q_{input}(t), \quad \text{on $\varGamma_I$}. \end{equation}
   In our benchmark scenario, we assume for simplicity that $q_{input}=0$, as if the top boundary represents the start of the stream. Such a stream is referred to as a first-order stream (see Strahler 1957), however, in real-world situations the first-order stream does not reach the catchment divide. The presented model can be also generalised to represent higher-order streams by including a non-zero upstream inflow $q_{input}(t)$.

Note that the kinematic approximation (3.5) that we follow in our work reduces the overland and channel equations to advective equations, rather than advective–diffusion equations. Thus, in this approximation, the downstream boundary conditions – at the river bank $\varGamma _R$ (for overland flow) and at the catchment outlet $\varGamma _O$ (for channel flow) – do not have to be imposed.

This means that, effectively, the channel flow does not impact the overland flow. However, overland flow impacts the channel flow thought the inflow term $q_{in}$ in (3.6). According to flow continuity, the input to the channel flow is the sum of the overland flow and the total groundwater flow, integrated over the entire channel perimeter at the given cross-section. Thus

(3.11)\begin{equation} q_{in} =\boldsymbol{q_s}|_{\varGamma_R}\boldsymbol{\cdot}\boldsymbol{n} + \int_{\varGamma_R} \boldsymbol{q_g}\boldsymbol{\cdot}\boldsymbol{n}\,\mathrm{d} l. \end{equation}

Two-way coupling between channel flow and subsurface flow is maintained via boundary condition (3.10c), and two-way coupling between the overland flow and subsurface flow is maintained via (3.10b).

3.5. Initial conditions of the benchmark

The choice of the initial condition is more arbitrary. In contrast to the benchmark scenarios by Maxwell et al. (2014), which assumed a constant groundwater depth, we select a more realistic setting, where the groundwater profile is given by its typical shape for a given catchment. Thus, we find a steady state of $h_g(x,y,z)$, $h_s(x,y)$ and $h_c(x,y)$ given by the time-independent versions of the governing equations (3.1), (3.3) and (3.6), solved for a given mean precipitation rate $R_{eff}=R_0$

(3.12a–c)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q_g} = 0, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q_s} + R_0 - I = 0, \quad\text{and}\quad {q_{in} - \frac{\partial q_c}{\partial \hat{y}}=0}. \end{equation}

Once this initial state is found by solving the above system of equations, we then explore the evolution of $h_g(x,y,z,t)$ and $h_s(x,y,t)$ caused by intensive rainfall, $R_{eff}>R_0$, which moves the system away from the initial state.

4. Governing equations (non-dimensional)

4.1. Non-dimensionalisation

The governing equations for subsurface, surface and channel flow presented in § 2 are now written in tilted coordinates $(\hat x, \hat y, \hat z)$ (cf. (2.1a–c)) and given in dimensional form in Appendix B.1. In order to understand the relative size of the terms appearing in the governing equations, we non-dimensionalise these equations. The following scalings are used:

(4.1)\begin{equation} \left.\begin{array}{c@{}} \hat{x} = L_{\hat x} \hat{x}', \quad h_g = L_z h_g', \quad t = t_0 t', \quad R_{eff} = r R_{eff}', \\ \hat{y} = L_y \hat{y}', \quad h_s = L_s h_s', \quad \theta(h) = \theta'(h'),\quad I = r I',\\ \hat{z} = L_z \hat{z}',\quad h_c = L_c h_c',\quad K_r(h) = K_r'(h'),\quad q_{in} = r L_{\hat x} q_{in}'. \end{array}\right\} \end{equation}

Here, $r$ is an average value of $R_{eff}$. We shall choose the characteristic time, $t_0$, overland water height, $L_s$, and channel water height, $L_c$, according to

(4.2a–c)\begin{equation} t_0 = \frac{L_z}{K_s}, \quad L_s = \left(\frac{L_{\hat x} n_s r}{S_x^{1/2}}\right)^{3/5}, \quad L_c = \left(\frac{n_c r L_{\hat x} L_y}{w S_y^{1/2}}\right)^{3/5}. \end{equation}

The choice of the above quantities comes from balancing the leading terms in the governing equations for subsurface, overland and channel flow, respectively. Their formulation, in terms of tilted coordinates, is presented in (B1), (B2) and (B3) in Appendix B.1.

Additionally, the non-dimensional terms in (4.2a–c) have straightforward physical interpretations. The time scale, $t_0$, describes a characteristic time that rainwater needs to penetrate the aquifer of thickness $L_z$, infiltrating with a characteristic speed $K_s$ (such flow occurs due to gravity if there is no hydraulic gradient, e.g. during uniform rainfall). The quantity $L_s$ represents the height of the overland flow at the river bank in a steady state with rainfall $r$ (assuming that the entire rainfall forms an overland flow, i.e. no infiltration appears). Similarly, $L_c$ is an approximate height of the flow in a wide channel at the river outlet in a steady state. Crucially, we note that the choice of the above scaling seems to be correct for our chosen benchmark, with all relevant dimensionless quantities of typical order unity in the numerical simulations of § 6.2.

It should be noted that, even though $t_0$ is a characteristic time of the vertical flow through the soil, other time scales are present. For example, we shall observe typically shorter time scales for the overland flow, and much longer time scales for the horizontal flow through the soil. Further discussion of the separation of time scales appears in Part 3 of our work.

4.2. Summary of governing equations and parameters

We collect the non-dimensional governing equations from Appendix B.2. To review, our hydrological problems in the 3-D geometry consist of solving three time-dependent partial differential equations for three unknowns: (i) a 3-D Richards equation for the subsurface flow (4.3a); (ii) a 2-D Saint Venant equation for the overland flow (4.3b); and (iii) a 1-D Saint Venant equation for the channel flow (4.3c). In the tilted frame, these are, respectively,

(4.3a)$$\begin{gather} \text{(Subsurface)} \quad \left.\frac{\mathrm{d}\theta}{\mathrm{d} h}\right\vert_{h=h_g}\frac{\partial h_g}{\partial t} = \mathcal{N}_1(h_g) + \beta_{zy}^2 \mathcal{N}_2(h_g) + \epsilon \beta_{zy} \mathcal{N}_3(h_g), \end{gather}$$

(4.3b)$$\begin{gather}\text{(Overland)} \quad \tau_s\frac{\partial h_s}{\partial t} = \frac{\partial }{\partial \hat{x}}(h_s^{5/3}) + R_{eff} - I, \end{gather}$$

(4.3c)$$\begin{gather}\text{(Channel)} \quad \tau_c\frac{\partial h_c}{\partial t} = q_{in} - \frac{\partial }{\partial \hat{y}}(h_c^{5/3}), \end{gather}$$

where the subsurface equations involve operator definitions

(4.4a)\begin{gather} \mathcal{N}_1(h_g) = \frac{\partial }{\partial \hat{z}} \left[K_r(h_g) \left(\frac{\partial h_g}{\partial \hat{z}}+1\right)\right] + \beta_{zx} S_x\frac{\partial }{\partial \hat{x}} \left[K_r(h_g) \left(2\frac{\partial h_g}{\partial \hat{z}}+1\right)\right] \nonumber\\ \quad + \,\beta_{zx}^2 \left(1 + S_x^2\right)\frac{\partial }{\partial \hat{x}} \left[K_r(h_g) \frac{\partial h_g}{\partial \hat{x}}\right] - \left.\frac{\mathrm{d}\theta}{\mathrm{d} h}\right\vert_{h=h_g}\frac{\partial h_g}{\partial t}, \end{gather}

(4.4b)$$\begin{gather} \mathcal{N}_2(h_g) = (1 + S_y^2) \frac{\partial }{\partial \hat{y}} \left[K_r(h_g)\frac{\partial h_g}{\partial \hat{y}}\right]\!, \end{gather}$$

(4.4c)$$\begin{gather}\mathcal{N}_3(h_g) = 2 \beta_{zx} (1+S_x^2) \frac{\partial }{\partial \hat{x}} \left[K_r(h_g)\frac{\partial h_g}{\partial \hat{y}}\right] + S_x\frac{\partial }{\partial \hat{y}} \left[K_r(h_g) \left(2\frac{\partial h_g}{\partial \hat{z}}+1\right)\right]\!. \end{gather}$$

Expressions for $\theta (h_g)$ and $K_r(h_g)$ are provided in Appendix B.2. Each partial differential equation in (4.3) is solved subject to boundary conditions posed on the domain boundaries given by (B8a)–(B8d).

Finally, these equations are characterised by nine independent dimensionless parameters, {$\beta _{zx}$, $\beta _{zy}$, $\sigma _x$, $\sigma _y$, $\tau _s$, $\tau _c$, $\gamma$, $\alpha$, $\rho$}, with definitions provided in Appendix C.

5. Model reduction to a two-dimensional model

In § 2, we formulated a general 3-D catchment model. The purpose of this section is to discuss the non-dimensionalisation of the model, which subsequently allows for the determination of the key dimensionless parameters governing the system. Once these are known, we may use the typical dimensional values established in Part 1 in order to compare the relative strengths of the various physical effects of the system.

We highlight two approximations:

1. (i) Considering either small river slope ($S_y\ll S_x$), short ($L_y\ll L_z$), or long catchment ($L_y\gg L_z$) approximations together with the necessary low channel limit ($L_c\ll L_z$), we may reduce the general 3-D governing equations for $h_g$ and $h_s$ to a 2-D form neglecting the flow along the $y$-axis.

2. (ii) In addition, in the case of the shallow aquifer scenario ($L_z\ll L_x$), we may apply a shallow-water approximation to further reduce the 2-D hillslope model to a 1-D model.

In this section, the approximations given in (i) are discussed. The regime of (ii) and its consequences are explored in Part 3 of our work.

5.1. Discussion of the low channel height limit, $L_c\ll L_z$

The 3-D model in $(x,y,z)$ can be formally approximated by a 2-D model in $(x,z)$ if the subsurface profile, $h_g(x,y,z,t) \sim h_{g,0}(x,z,t)$, and the surface profile, $h_s(x,y,t) \sim h_{s,0}(x,t)$, and both asymptotic approximations are consistent with the initial and boundary conditions (at the leading order).

We observe that the boundary condition (3.10c), along the channel, $\varGamma _R$, depends in general on the channel water height, $h_c(y,t)$, which can vary along the catchment. For example, the dimensional height can vary from $h_c=0$ at $y=L_y$ (if there is no inflow to the river from the upstream point) to a dimensional height $h_c=L_c$ at $y = 0$ (in the case of the steady-state outflow). Returning to non-dimensional values for $h_g$, $h_c$ and $z$, (3.10c) yields

(5.1)\begin{equation} h_g = \frac{L_c}{L_z}h_c - z, \quad \text{on $\varGamma_R$}. \end{equation}

Typically, the values of $L_c/L_z$ are very small: based on the UK catchment data from Part 1, we can extract the interquartile range for $L_c/L_z$, namely $[0.0011, 0.0147]$ (i.e. the middle half of UK catchments have $L_c/L_z$ within this interval). Thus, even though the channel water height may vary along the channel, it is negligibly small comparing with the typical variation of the pressure head. In the limit of $L_c/L_z \to 0$, we see that the subsurface boundary condition is

(5.2)\begin{equation} h_g \sim-z, \quad \text{on $\varGamma_R$}, \end{equation}

which is no longer $y$-dependent.

5.2. An asymptotic expansion for small river slopes, in $\epsilon = S_y/S_x$

Although the remaining boundary conditions ((3.10a)–(3.10d) without (3.10c)) are not explicitly $y$-dependent, the solution $h(x,z,t)$ may still exhibit leading $y$-dependent effects due to e.g. the topography. However, there are certain approximations in which these effects are very small – for example, when the slope along the channel $S_y$ is much lower than the slope along the hillslope $S_x$. Note that the aspect ratio introduced in § 2.1, $\epsilon = S_y/S_x$, typically has small values (half of UK catchments have $\epsilon$ between $0.13$ and $0.25$). Here, we shall demonstrate that when $\epsilon \ll 1$ (equivalent to $S_y\ll S_x$), the solution is expected to be predominantly two-dimensional.

Firstly, we rewrite the set of dimensionless governing equations for the subsurface and overland flows, (B4) and (B5), in a simpler form highlighting its structure

(5.3a)$$\begin{gather} \text{(Subsurface)} \quad \left.\frac{\mathrm{d}\theta}{\mathrm{d} h}\right\vert_{h=h_g}\frac{\partial h_g}{\partial t} = \mathcal{N}_1(h_g) + \beta_{zy}^2 \mathcal{N}_2(h_g) + \epsilon \beta_{zy} \mathcal{N}_3(h_g), \end{gather}$$

(5.3b)$$\begin{gather}\text{(Overland)} \quad \tau_s\frac{\partial h_s}{\partial t} = \frac{\partial }{\partial \hat{x}}(h_s^{5/3}) + R_{eff} - I, \end{gather}$$

where the nonlinear operators, $\mathcal {N}_i$, for $i = 1, 2, 3$ are defined in (B13) in Appendix B.2. Note that these operators are dependent on $h_g$ and $h_s$, and independent of $\epsilon$ and $\beta _{zy}$, which are the only dimensionless parameters involving $S_y$ and $L_y$.

When $\epsilon = 0$, we can verify that the solutions are independent of $\hat {y}$, i.e. they can be written as $h_g(\hat x,\hat y,\hat z)=h_{g,0}(\hat x,\hat z)$ and $h_s(\hat x,\hat y)=h_{s,0}(\hat x)$. This is caused by the combination of three facts:

1. (i) term $\mathcal {N}_1(h_g)$ is independent of $\hat {y}$;

2. (ii) operators $\mathcal {N}_2$ and $\mathcal {N}_3$ applied to a function independent of $\hat {y}$ become $0$; and

3. (iii) the no-flux boundary condition at $\hat {y}=0, \, 1$ in (B8a) is then

   (5.4)\begin{equation} \frac{\partial h_g}{\partial \hat{y}}-\frac{\epsilon}{1-\epsilon^2} \left(\frac{\beta_{zx}}{\beta_{zy}}\frac{\partial h_g}{\partial \hat{x}}-\epsilon\frac{\partial h_g}{\partial \hat{y}} -\frac{2S_x}{\beta_{zy}}\frac{\partial h_g}{\partial \hat{z}}\right)= 0. \end{equation}
   Hence, for $\epsilon =0$, the above boundary condition is satisfied by $h_g=h_{g,0}(\hat x,\hat z)$.

4. (iv) From (5.3b) only $R_{eff}$ and $I$ can be $\hat {y}$-dependent terms, but in the considered scenario $R_{eff}$ is constant, and $I(\hat {x},\hat {y})$ is $\hat {y}$ independent as long as $h_g$ is.

Essentially, $\epsilon =0$ is associated with a zero gradient along the river, i.e. there is no forcing flow in the $\hat {y}$-direction, and the domain becomes transitionally symmetric in that direction.

There is an important consideration in the formal limit as $S_y\propto \epsilon \to 0$. In this limit, holding other parameters fixed, the $L_c$ defined in (4.2a–c) tends to $\infty$, and so the $L_c\ll L_z$ condition from (5.1) is no longer satisfied. This is due to the fact that when reducing the gradient along the channel, $S_y$, the channel water height must increase in order to maintain a significant channel flow (cf. Manning's law (3.9)). Therefore, we would expect for a 2-D dynamics to dominate when

(5.5)\begin{equation} \left(\frac{n_c r L_{\hat x} L_y}{w L_z^{5/3}}\right)^2\ll S_y\ll S_x. \end{equation}

Above, the expression on the left-hand side is obtained from the definition of $L_y$ from (4.2a–c), and represents the value of $S_y$, for which $L_c=L_z$. In real-world situations, $S_y$ (with a median value of $0.014$ based on the data collected in Part 1) is higher by a few orders of magnitude over this threshold (with a median value of $6.1\times 10^{-11}$), and so only the second approximation in (5.5) needs to be considered.

5.3. Asymptotic expansions for short ($\beta _{zy} \gg 1$) and long ($\beta _{zy} \ll 1$) catchments

There are additional limits that allow us to reduce the 3-D problem into simpler 2-D formulations at the leading order, and these involve the non-dimensional geometrical parameter

(5.6)\begin{equation} \beta_{zy} \equiv \frac{L_z}{L_y}. \end{equation}

For instance, in the limit as $\beta _{zy} \to \infty$, the 3-D catchment reduces to an infinitely thin hillslope profile with a negligible flow in the perpendicular direction to the hillslope (since we imposed no-flow conditions at $\hat {y}=0$ and $\hat {y}=1$). Equivalently, this corresponds to an asymptotically short section of a river. From (5.4), we see that the leading-order profile should satisfy the $\operatorname {d\!}{}{h_{g,0}}/\operatorname {d\!}{}{\hat {y}}=0$ condition at $\hat {y} = 0, 1$, which is automatically satisfied for a $\hat {y}$-independent solution. As argued in the previous section, we also conclude that such a $\hat {y}$-independent solution will also satisfy the governing equations (5.3).

Using a similar analysis to the one presented in the previous section, by balancing leading terms in the boundary conditions (5.4), we can show that the full three-dimensional solution can be expanded in terms of $\beta _{zy}^{-1}$

(5.7a)$$\begin{gather} h_g(\hat x,\hat y,\hat z) =h_{g,0}(\hat x,\hat z)+\beta_{zy}^{-1} h_{g,1}(\hat x,\hat y,\hat z) + O(\beta_{zy}^{-2}), \end{gather}$$

(5.7b)$$\begin{gather}h_s(\hat x,\hat y) = h_{s,0}(\hat x)+\beta_{zy}^{-1} h_{s,1}(\hat x,\hat y) + O(\beta_{zy}^{-2}). \end{gather}$$

The last interesting limit we discuss is $\beta _{zy}\rightarrow 0$, which corresponds to the situation of an asymptotically long river. Similarly, the 2-D solution satisfies the governing equations (5.3). However, this time, it does not satisfy the no-flow boundary condition (5.4). Therefore, we expect to observe a boundary layer around $\hat {y}=0$ and $\hat {y}=1$ (see figure 5). Consequently, $h_{g,0}$ and $h_{s,0}$ are understood to represent the ‘outer’ asymptotic solutions, valid for $0 < \hat {y} < 1$.

Figure 5. This graphic shows the steady-state depth of the groundwater table, according to the 3-D model. We note that the solution is mostly $\hat {y}$-independent, except for two apparent boundary layers around $\hat {y}=0$ and $\hat {y}=1$; near these points, the groundwater table aligns with the lines of constant elevation. The figure is generated using the solver described in § 6 using the parameter values given in table 1, except for two values: we use $L_y=13680$ m ($\beta _{zy}=0.05$) and $S_y=0.0075$ ($\epsilon =0.1$); as a result, this graphic matches an inset in figure 9. The boundary layer thicknesses, $\delta _1$ and $\delta _2$, tend to zero as $\beta _{zy} \to \infty$.

Table 1. Columns on right present the default values and ranges of the parameters used to perform a sensitivity analysis. Columns on the left present parameters which were not varied during the sensitivity analysis. The length of the catchment was selected to be $L_y^{trib}$ estimated from Part 1, which represents the typical distance between the main tributaries of the river.

Without loss of generality, let us consider the boundary layer near $\hat {y} = 0$. We rescale $\hat {y}=\delta \hat {y}'$ where $\delta (\beta _{zy})$ is a characteristic size of the boundary layer. After applying this transformation, the governing equations (5.3) become

(5.8a)$$\begin{gather} \text{(Subsurface)} \quad \left.\frac{\mathrm{d}\theta}{\mathrm{d} h}\right\vert_{h=h_g}\frac{\partial h_g}{\partial t} = \mathcal{N}_1(h_g) + \frac{\beta_{zy}^2}{\delta^2} \mathcal{N}_2(h_g) + \epsilon \frac{\beta_{zy}}{\delta} \mathcal{N}_3(h_g), \end{gather}$$

(5.8b)$$\begin{gather}\text{(Overland)} \quad \tau_s\frac{\partial h_s}{\partial t} = \frac{\partial }{\partial \hat{x}}(h_s^{5/3}) + R_{eff} - I. \end{gather}$$

The 3-D effects given by $\mathcal {N}_2$ and $\mathcal {N}_3$ become significant when the characteristic size of the boundary layer is of the order $O(\beta _{zy})$. Hence, we conclude that the thickness of the boundary layer is $O(1/\beta _{zy})$. Therefore, if we consider the ${\beta _{zy}\rightarrow 0}$ limit, the solution for the problem becomes two-dimensional except for an infinitely thin boundary layer around the boundaries.

5.4. Summary of the two-dimensional model

To summarise, we considered three approximations for small river slope ($\epsilon \ll 1$), short catchments ($\beta _{zy}\gg 1$), and long catchments ($\beta _{zy}\ll 1$). We showed that, under any of these approximations, the model can be approximately represented in the following 2-D form:

(5.9a)$$\begin{gather} \text{(Subsurface)} \quad \left.\frac{\mathrm{d}\theta}{\mathrm{d} h}\right\vert_{h=h_{g,0}}\frac{\partial h_{g,0}}{\partial t} = \mathcal{N}_1(h_{g,0}), \end{gather}$$

(5.9b)$$\begin{gather}\text{(Overland)} \quad \tau_s\frac{\partial h_{s,0}}{\partial t} = \frac{\partial }{\partial \hat{x}}(h_{s,0}^{5/3})+R_{eff}-I, \end{gather}$$

with

(5.10)\begin{align} \mathcal{N}_1(h_{g,0}) &= \frac{\partial }{\partial \hat{z}} \left[K_r(h_{g,0}) \left(\frac{\partial h_{g,0}}{\partial \hat{z}}+1\right)\right] + \beta_{zx} S_x\frac{\partial }{\partial \hat{x}} \left[K_r(h_{g,0}) \left(2\frac{\partial h_{g,0}}{\partial \hat{z}}+1\right)\right] \nonumber\\ &\quad + \beta_{zx}^2 (1 + S_x^2)\frac{\partial }{\partial \hat{x}} \left[K_r(h_{g,0}) \frac{\partial h_{g,0}}{\partial \hat{x}}\right]\!. \end{align}

Both the groundwater and overland flows reaching the channel contribute to the river flow in the $\hat {y}$-direction governed by (B6). Thus

(5.11)\begin{equation} \text{(Channel)} \quad \tau_c\frac{\partial h_c}{\partial t} = q_{in} - \frac{\partial }{\partial \hat{y}}(h_c^{5/3}), \end{equation}

where the inflow $q_{in}=q_{in}(t)$ is given by the total groundwater and overland inflow to the channel. This equation allows us to find the hydrograph at the outlet of the catchment $Q(t)=q_c(\hat {y}=0,t)$. However, note that after the reduction to a 2-D model, (5.9a) and (5.9b) are uncoupled from (5.11), and so $q_{in}$ can be found without solving the last equation. Therefore, we can study the properties of the 2-D model by solving (5.9a) and (5.9b) alone, and exploring the properties of the river inflow term $q_{in}(t)$. We follow this approach in the study of the 2-D model in § 8, and in Part 3 of our work. The above system of partial differential equations forms a model describing surface and subsurface flow in a 2-D hillslope cross-section, as presented in figure 6. The boundaries are now one-dimensional, but the boundary conditions are the same as in the 3-D model, as given by (B8a)–(B8d). As before, for the initial condition, we consider the steady state of the above system for $R_{eff}=R_0$. In § 7, we will explore the accuracy of this approximation numerically, investigating the size of 3-D features of the full solution and their behaviour in limits formulated above.

Figure 6. An illustration of a 2-D hillslope geometry in Cartesian coordinates (a) and in the tilted coordinate system (b). These represent a cross-section (in the $xz$-plane) of the original V-shaped catchment.

A two remarks are in order:

1. (i) Firstly, there are some remaining terms in the dimensionless governing equations, which are small; however, they should not be neglected. For example, coefficient $\tau_s = L_s/(t_0 r)$ is small, which means that the characteristic time scale of accumulation of surface water $(L_s/r)$ is much shorter than the characteristic time scale of vertical subsurface flow $(t_0)$. However, temporal term of (5.9b) becomes significant for small values of $t$. Since the typical rainfall times are much lower than the characteristic time of groundwater transfer (estimated as $t_0\approx 6.8\times 10^7~\mathrm {s}\approx 2$ years), we are often interested in the short-time behaviour, and therefore, this term should not be neglected.

2. (ii) Secondly, the $\beta _{zx}$ term is also small in the shallow aquifer scenario compared with the leading term representing the flow in the $\hat z$-direction, and therefore it can be neglected in regions with significant temporal effects (e.g. in partially saturated zones impacted by the rainfall). However, in the fully saturated zone, where $h_g>0$, we have ${\mathrm {d}\theta }/{\mathrm {d} h}=0$. In this zone, the balance between the two remaining terms needs to be maintained – the horizontal flow becomes high enough to balance the vertical flow. Therefore, the $\beta _{zx}$ term cannot be neglected in the fully saturated zone; however, another simplification based on the shallow-water approximation can be considered. This will be explored in the Part 3 of our work.

6. Numerical methodology

In order to validate the reduction of the 3-D model to the 2-D approximation and quantify the impact of model parameters on the observed peak flow, we follow a numerical approach.

Here, we present the numerical method used to implement the coupled surface– subsurface model based on the governing equations introduced in § 3. To summarise, these are

(6.1a)$$\begin{gather} \text{(Subsurface)} \quad \frac{\mathrm{d}\theta}{\mathrm{d} h_g}\frac{\partial h_g}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\cdot} (K_s K_r(h_g)\boldsymbol{\nabla}(h_g+z)), \end{gather}$$

(6.1b)$$\begin{gather}\text{(Overland)} \quad \frac{\partial h_s}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{1}{n_s}h_s^{5/3} \frac{\boldsymbol{S}_0}{\sqrt{|\boldsymbol{S}_0|}}\right)+R_{eff}-I, \end{gather}$$

(6.1c)$$\begin{gather}\text{(Channel)} \quad w\frac{\partial h_c}{\partial t} = q_{in} - w\frac{\sqrt{S_y}}{n_c}h_c^{5/3}, \end{gather}$$

subject to boundary conditions (3.10a)–(3.10d).

Our numerical implementation has two applications in this study. Firstly, in § 7 we use the numerical scheme based on a discrete version of this equations to verify our reductions to the 2-D problem introduced in § 5. Secondly, in § 8.2 we numerically analyse features of a benchmark scenario in a reduced 2-D analysis. We use the same equations as above, excluding $y$-dependent terms and channel flow. The source codes were written in MATLAB and are available in our GitHub repository (Morawiecki 2022).

6.1. Model discretisation

The implementation of the 3-D and 2-D models is based on the finite volume method. The entire hillslope is divided into $N_x\times N_y\times N_z$ cells. Here, $N_z$ is additionally split into $N_{z,1}$ cells representing the layer of soil with the same depth as the channel, and $N_{z,2}=N_z-N_{z,1}$ cells representing deeper layers of the aquifer, as illustrated in figure 7(a). In the case of the shallow aquifer scenario, we set $N_{z,2}=0$. The implementation allows for a mesh refinement by varying the cells extent, according to a geometric series (see figure 7b).

Figure 7. (a) Discretisation of the 3-D catchment, representing the V-shaped catchment scenario. In the case of the 2-D deep aquifer and shallow aquifer model, we set $N_y=0$. (b) Example of mesh refinement. The size of edges is given by the geometric series with ratios $\mu$ and $\nu$.

Depending on the type of simulation, we handle the channel differently. In the case of 2-D simulations, we will assume the water level in the channel to be equal to the channel depth (unless stated otherwise). In the case of 3-D simulations, $h_c(y,t)$ will be iteratively computed. However, following the V-shaped catchment scenario by Maxwell et al. (2014), we consider the channel flow in the same way as the overland flow, just with a different Manning coefficient. This way, we neglect the interactions with the river banks; however, the simulations are significantly more stable.

A value of $h_g$ is assigned to each cell to represent $h_g$ at the centroid of the cell. Due to the pressure continuity condition at the surface, $h_s=h_g$, so there is no need to consider $h_s$ at the surface as an independent variable (the same applies to $h_c$ in 3-D simulations). One challenge is the significant difference in the time scales between the overland and groundwater flow (the ratio of which is quantified with the dimensionless parameter $\tau _s\approx 2.8\times 10^{-4}$). A stable numerical scheme for overland flow requires a shorter time step than groundwater flow. Therefore, for each groundwater flow step, we compute several steps of the surface flow, while simultaneously satisfying the continuity boundary condition at the surface.

The groundwater in each time step is found using an implicit scheme. The following discretised version of Richards equation (B4) for each cell is used:

(6.2a)\begin{equation} V_i\frac{\theta(h_i^{\prime t+1})+ \left.\dfrac{\mathrm{d}\theta}{\mathrm{d} h}\right\vert_{h_i^{\prime t+1}} (h_i^{t+1}-h_i^{\prime t+1})-\theta(h_i^t)}{\Delta t} = \sum_{j\in\,{neighbours}} \mathcal{G}_{ij}, \end{equation}

where

(6.2b)\begin{align} \mathcal{G}_{ij} &=\frac{S_{i,j}}{\cos(\eta_{i,j})} \left(K_{i,j}'+\frac{\mathrm{d} K_{i,j}'}{\mathrm{d} h}(h_{u(i,j)}^{t+1} -h_{u(i,j)}^{\prime t+1})\right) \nonumber\\ &\quad \times \frac{h_j^{t+1}+z_j-h_i^{t+1}-z_i}{\|\boldsymbol{r}_{i\rightarrow j}\|}\frac{\boldsymbol{\beta}\boldsymbol{\cdot}\boldsymbol{r}_{i\rightarrow j}}{\|\boldsymbol{r}_{i\rightarrow j}\|}. \end{align}

Few remarks should be made here:

1. (i) The left-hand side represents the estimated rate of change of the $i$th cell's water content. Here, $V_i$ is the cell's volume, $\Delta t$ is the time step duration, $h_i^t$ and $h_i^{t+1}$ are the pressure head ($h_g$) in cell $i$ at time step $t$ (previous one) and $t+1$ (current one), respectively, $h_i^{\prime t}$ is the pressure head computed in the previous iteration of the implicit scheme and $\theta$ is the saturation given by the MvG model (B7).

2. (ii) The right-hand side represents the sum of all flows between cell $i$ and its neighbours. Here, $S_{i,j}$ is the area of the face between cell $i$ and $j$, $\eta _{i,j}$ is the angle between this face and the line joining these cells’ centroids, $\boldsymbol {r}_{i\rightarrow j}$ is the vector from the centroid of cell $i$ to the centroid of cell $j$ and $z_i$ is the $z$ coordinate of the $i$th cell centroid. Additionally, $\boldsymbol {\beta }=(\beta _{zx}^2,\beta _{zy}^2,1)$ is a vector describing the anisotropy coming from different scaling factors in the dimensionless model. Also, $K_{i,j}'$ is the hydraulic conductivity of the face between cell $i$ and $j$. It is computed using the upwind scheme, i.e. its value is computed using MvG model (B7b) for $h=h_{u(i,j)}^{\prime t+1}$, where $u_{i,j}=i$ if the flow is going from cell $i$ to $j$, and $u_{i,j}=j$ otherwise.

3. (iii) The change of both $\theta$ and $K$ is estimated using the first two terms of the Taylor series. If the time step is short enough, the algorithm converges to $h_i^{t+1}$ satisfying the continuity condition. Equation (6.2) is linear in $h_i^{\prime t+1}$ for all $i$, and therefore, it can be solved using standard methods for linear algebraic equations.

4. (iv) After each iteration of groundwater flow, a number of iterations of overland flow is performed. In order to guarantee numeric stability, the Courant number defined as

   (6.3)\begin{equation} C=\frac{u_{i,j}\Delta t}{\|\boldsymbol{r}_{i\rightarrow j}\|}\quad\text{with } u_{i,j}=K_{i,j}'\frac{h_j^{t+1}+z_j-h_i^{t+1}-z_i}{\|\boldsymbol{r}_{i\rightarrow j}\|}, \end{equation}
   where $u_{i,j}$ represents the flow speed between cell $i$ and $j$, should be lower than $1$ for all pairs of computational cells. In order to achieve this, an adaptive time stepping is used to keep the Courant number at a given threshold value; however, additionally, a minimum number of iterations is also preset to maintain high accuracy.

After each groundwater solver step for each cell with a face on the land surface, we compute the total volume of the water (surface and subsurface) divided by the total area of the top/bottom face ($\Delta x \Delta y$). Let us denote this ratio as $f_{i,j}$, where $i$ and $j$ are the given cell's indices. In each iteration of the surface solver, $h_s$ (and $h_c$ in 3-D simulations) is computed for each cell as $h_{i,j}=f_{i,j}-f^{max}_{i,j}$ for $f_{i,j}>f^{max}_{i,j}$ and $h_{i,j}=0$ otherwise. Here, $f^{max}_{i,j}$ is the maximum possible value of $f_{i,j}$, corresponding to a saturated cell with $h_g=0$. Then $f_{i,j}$ is updated using the following explicit scheme for flow given by the discretised form of the 2-D Saint Venant equation (B5)

(6.4)\begin{equation} f_{i,j}^{t+1} = f_{i,j}^t + \frac{1}{\Delta x}\big(h_{i+1,j}^{5/3}-h_{i,j}^{5/3}\big)\sqrt{S_x} \frac{\Delta t}{n_s} + R_{eff}. \end{equation}

After the last iteration of the overland flow scheme, $f_{i,j}$ values are used to calculate the pressure head $h_g$ in the subsurface cells bordering the land surface, after which the next time step for subsurface flow is computed. After each step we evaluate the output flow. In the case of the 3-D model, we calculate the river flow at the outlet $Q(t)=q_c(\hat y=0, t)$, and in the case of the 2-D model, we calculate the river inflow $q_{in}(t)$. These functions can then be presented in the form of a hydrograph.

In addition to the above time-dependent solver, a steady-state solver was also implemented. It is based on the discretisation in (6.2), where the left-hand side (temporal term) is equal to $0$. The overland flow is included as an additional flow component between the surface cells and is solved simultaneously with the Richards equation.

The implementation described in this section was verified by successfully replicating the benchmark results by Sulis et al. (2010) obtained for a hillslope using the ParFlow integrated catchment model, and by Maxwell et al. (2014) for the V-shaped catchment using the process-based adaptive watershed simulator (PAWS) model (and other coupled surface–subsurface models). The results of this comparison are presented in Appendix D.

6.2. Example three-dimensional solution

Before proceeding to the quantitative analysis, we dedicate this section to a qualitative discussion of the general properties of a typical solution for the presented model. Let us consider a scenario of an intensive rainfall over the V-shaped catchment that initially remains in equilibrium with the mean precipitation characterising the given region. In the experiments, we find a steady-state solution for a mean precipitation $r_0$, which is set as the initial condition. We then simulate the reaction of the system to a constant precipitation $r>r_0$, and analyse the resulting river flow.

In the simulations, unless stated otherwise, all catchment parameters will be set to the typical values characterising UK catchments extracted in Part 1, which are summarised in table 1. For the catchment length, $L_y$, we took the average length of the river between consecutive tributaries, $L_y^{trib}=945$ m, since at this scale, the drainage network no longer exhibits its fractal-like finger pattern. This way, we can treat our benchmark model as a representation of a single first-order catchment (Dietrich & Dunne 1993), or a first-/higher-order stream, forming the base element of more complex drainage networks (Strahler 1957). The simulation results, covering a 24 hour rainfall, are collected and presented in figure 8.

Figure 8. The illustration summarises the key properties of the 3-D solution obtained for the default values of parameters from table 1. Panel (a) depicts the initial condition (steady state for mean rainfall $r_0=2.95 \times 10^{-8}\ \textrm {m}\ \textrm {s}^{-1}$). During the subsequent rainfall $r=2.36 \times 10^{-7}\ \textrm {m}\ \textrm {s}^{-1}$, the water level in the channel rises, as presented in (b), causing the flow at the river outlet to increase, as depicted by the hydrograph (c). We observe slightly different solution profiles depending on the location along the $\hat y$ axis – their main features are outlined in cross-sections (d,e). The surface water height $h_s$ was magnified 5000 times, to make it visible.

Initially, the system remains in a steady state, in which the pressure head $h_g$ increases with depth following an approximately hydrostatic profile (figure 8a). The interesting dynamics responsible for generating the flow occurs near the surface ($\hat {z}=1$).

The surface water is present near the channel and extends further away from it for lower $\hat {y}$ values (figure 8d–f). However, around $\hat y=1$, we do not observe surface water at all. This is caused by a non-zero gradient $S_y$ along the $\hat {y}$ direction, forcing the groundwater flow in that direction. We will refer to the zone in which the groundwater reaches the surface and forms an overland flow as the seepage zone (see figure 6).

Two distinct effects are observed when the intensive rainfall starts. Firstly, the rainfall over the seepage zone starts to accumulate, causing the surface water to quickly rise (as highlighted in figure 8e). Increased overland flow, leads to a rapid rise of the channel water and resulting outflow from the catchment within the first hour (figure 8b,c).

Secondly, the rainfall outside the seepage zone starts to infiltrate the unsaturated section of the soil, forming a characteristic wetting front (as highlighted in figure 8d). After the infiltrating water reaches the groundwater, its level starts to rise. The rising groundwater eventually reaches the surface, causing the growth of the seepage zone (as in figure 8f), increasing the area from which overland flow reaches the river.

An essential observation for this time-dependent solution is that the characteristic time scale of overland flow dynamics is much shorter than the characteristic time of groundwater flow (their ratio is given by the dimensionless parameter $\tau _s\approx 2.8\times 10^{-4}$). This time scale separation is reflected in the shape of the hydrograph in figure 8(c), which shows the dependence between river flow at the outlet $Q(t)=q_c(\hat y,t)$ and time $t$.

A multiscale behaviour can be observed, with an early-time fast rise of total flow dominated by a rising overland flow fed by the rainfall over the seepage zone, and a late-time slow rise of total flow caused by rising groundwater and the resulting slow expansion of the seepage zone over time. This observation allows us both to understand the importance of model parameters (§ 8.2) and to further simplify the problem in Part 3. Note that, for the typical parameters studied in this work, the outlet flow will continue rising, eventually reaching a new steady state with $\lim _{t\rightarrow \infty }Q(t)=rA$. However, in case of the default simulation settings discussed above, it requires months of constant high rainfall for the solution to approach the new steady state. Thus, this effect is not observable over the typical day-long scales of the presented simulations.

7. Verification of 3-D to 2-D reduction

The time-dependent simulations presented in the previous section demonstrate some of the 3-D features that are visible in the simulations. In this section, we further investigate these features, and show how they depend on two model parameters characterising the catchment geometry along the $\hat {y}$ direction, namely the catchment length $L_y$ and the slope parallel to the channel $S_y$. Alternatively, in terms of the non-dimensional quantities, this corresponds to $\beta _{zy}$ and $\epsilon$.

7.1. Three-dimensional features of steady-state solution

In order to develop a better understanding of the 3-D effects, we performed a series of numerical experiments in which we found a steady state for varying values of $L_y$ and $S_y$, while keeping other parameters constant with the values provided in table 1. The groundwater table shape corresponding to the selected steady states is presented in figure 9.

Figure 9. Groundwater table depth in steady states obtained for varying catchment length $L_y=\beta _{zy}^{-1}L_z$ and slope $S_y=\epsilon S_x$. Dashed lines represent lines of constant elevation. The entry with $\epsilon = 0.1$ and $\beta _{zy} = 0$ is the same figure as presented in figure 5.

In all cases we observe that the solution becomes less $\hat {y}$-dependent as $\epsilon \rightarrow 0$, as expected from § 5. However, the dependence on $\beta _{zy}=L_z/L_y$ is more complex. The phase space can be divided into three regions:

1. (i) When $L_y\ll {L_{\hat x}}/{\epsilon }$, the lines of constant elevation are approximately perpendicular to the hillslope (e.g. $\epsilon =0.1$, $L_y=100$). As shown in § 5, in such a case, the leading-order (2-D) solution of the governing equations for small $\epsilon$ also satisfies the boundary conditions in the leading order. In this case, we observe that the leading-order solution follows the lines of constant elevation over the entire domain.

2. (ii) When $L_y \gg {L_{\hat x}}/{\epsilon }$, the lines of constant elevation are approximately parallel to the hillslope (e.g. $\epsilon =0.1$, $L_y=10^6$). In such a case, the leading-order (2-D) solution for the governing equations does not satisfy the flow boundary condition at $\hat y=0$ and $\hat y=1$. As a consequence, a boundary layer is developed near these two boundaries, in which the lines of constant groundwater table depth become parallel to lines of constant elevation, while in the outer solution they become perpendicular to the hillslope. The thickness of these boundary layers $\delta$ decreases inversely proportional to $L_y$ (see figure 10), which is consistent with the theoretical scaling derived in § 5.3.

3. (iii) In the intermediate region, when $L_y = O({L_{\hat x}}/{\epsilon })$ and $\delta =O(1)$ (e.g. $\epsilon =0.1$ and $L_y=3000$), the leading-order solution does not satisfy the boundary conditions, and the ‘boundary layer’ thickness becomes large enough to impact the solution over a major part or even effectively the entire domain. In such cases, the solution does not seem to satisfy the 2-D approximation unless $\epsilon$ is small enough.

Figure 10. Boundary layer thickness at $\hat {y}=0$ (a) and $\hat {y}=1$ (b) as a function of $\beta _{zy}$. The boundary thickness was measured based on the groundwater depth profile along $\hat {x}\approx 0.46$. The boundary was defined as $\hat {y}$, for which the groundwater depth $H(0.46,\hat {y})$ is further than $\pm 5\,\%$ from the groundwater depth evaluated in the middle of the domain, $H(0.46, \hat {y}=0.5)$. For small $\beta _{zy}$ values, the boundary width follows $\delta \propto \beta _{zy}$ scaling. As $\beta _{zy}$ increases, $\delta$ reaches $0.5$, for which the boundary condition affects effectively the entire domain.

7.2. Analysis of 3-D to 2-D reduction error

Here, we follow the qualitative discussion from the previous section by quantifying the difference between the full solution and its 2-D approximation.

In § 5, we argued that the solution for the 3-D problem can be represented as $h_g(\hat x,\hat y,\hat z, t)=h_{g,0}(\hat x,\hat y, t)+\epsilon h_{g,1}(\hat x,\hat y,\hat z, t)+O(\epsilon ^2)$, where $h_{g,0}$ is a 2-D solution for $\epsilon =0$. In this section, we verify this theoretical result numerically in the case of a steady-state solution, to which the same argument also applies.

Finally, we estimated the mean absolute difference between $h_g$ and $h_{g,0}$ by averaging its values for all computational cells weighted by their volume. The dependence of this mean error on $L_y$ and $\epsilon$ is presented in figure 11. It confirms the asymptotic analysis from § 5 and the qualitative observations from § 7.1. Firstly, it confirms that the error of the 2-D approximation increases proportionally to $\epsilon$ for large values of $\epsilon$. However, for small values of $\epsilon$, the error increases, because the effect of $y$-dependent water height at the channel is no longer negligible (see § 5.2). Secondly, it shows that the error is small for very small values of $L_y$ (2D solution is satisfied everywhere) and very large $L_y$ values (2-D solution is satisfied everywhere apart from a thin boundary layer at $\hat {y}=0$ and $\hat {y}=1$), but the error is the highest for intermediate values (here around $L_y=3000$).

Figure 11. The mean absolute difference between the full 3-D solution and its 2-D approximation for small $\epsilon$ as a function of $\epsilon$ and $\beta _{zy}$.

8. Impact of physical parameters on the 2-D model

Following § 5, the inflow to the river in our benchmark scenario for $S_y\ll S_x$ can be approximated by a 2-D model. In this section, we use the numerical procedure described in § 6 to quantify the impact of model parameters on the peak flows observed after an intensive rainfall and link them to the key physical processes accounting for the flow generation.

8.1. Structure of typical hydrographs

In this section, we examine, in more detail, the hydrographs that correspond to two representative simulations. Under many sensible parameter choices, we have observed that many flow experiments can be roughly described into these two prototypical classes. Although this is only described qualitatively in this work, we shall justify it more rigorously using the reduced model of Part 3.

For these simulations, we use the same parameter values as in § 6.2, with a catchment initially remaining in a steady state with rainfall $r_0$. The rainfall then rises to $r>r_0$ at $t=0$. Additionally, we set $S_y=0$ to reduce the problem dimension.

The numerical simulations are based on two experiments differentiated by their $r_0$ values: experiment (A) with $r_0 = 2.95 \times 10^{-8}\ \textrm {m}\ \textrm {s}^{-1}$; and experiment (B) with $r_0 = 2 \times 10^{-9}\ \textrm {m}\ \textrm {s}^{-1}$. The two corresponding hydrographs, $Q(t)$ vs $t$, are presented in the top-left insets of figures 12 and 13. For each hydrograph, solutions are presented at four times and given in insets (a) through (d). In the insets, areas shaded blue represent the saturated groundwater zone (with $h_g>0$), while shaded green areas represent the unsaturated zone. Surface water height was magnified $2000$ times, and its initial height was highlighted in a darker blue. Only a small part of the catchment near the river is presented.

Figure 12. Numerical solution of 2-D model for $r_0=2.95 \times 10^{-8}$ m s$^{-1}$ with an initial seepage zone. All other parameters were set to the default values presented in table 1.

Figure 13. Numerical solution of 2-D model for $r_0=2 \times 10^{-9}$ m s$^{-1}$ without an initial seepage zone. All other parameters were set to the default values presented in table 1.

The main difference between these two hydrographs is the existence of surface water in the initial condition. We observe that in the case presented in figure 12, the groundwater flow is not sufficient to transfer rainwater to the channel, and a fraction of the catchment area (namely the seepage zone) is initially covered with surface water. In contrast, in figure 13, initially there is no overland flow, i.e. the groundwater never reaches the surface (except for the channel boundary). There is a significant qualitative difference between these two cases.

8.1.1. Experiment (A): a case with an initial seepage zone

In experiment (A), we observe that the hydrograph can be roughly divided into two phases. We propose that in these two phases, the flow increase is determined by different physical mechanisms (similar to the 3-D case presented in figure 8).

During an early time phase (roughly first $12$ minutes), we observe a significant rise in total flow reaching the river. This corresponds to evolution between states (a) and (b) in figure 12. We may interpret this early-time rise as a result of rainfall accumulating over the seepage zone, enhancing the already existing overland flow. This causes the river flow to rise by the rainfall excess over the initial seepage zone $(r-r_0)A_s$, where $A_s$ is the initial area of the seepage zone (which can be measured in the simulation).

As a result, in a short time, we define the flow

(8.1)\begin{equation} Q_{crit} \equiv r_0A+(r-r_0)A_s. \end{equation}

The above quantity we shall refer to as the critical flow. Here, $r_0A$ represents the initial flow, and $A=L_xL_y$ is the catchment's area. In fact, the dashed line plotted in the hydrograph of figure 12 is calculated via (8.1) and seems to coincide with the change in gradient of the hydrograph.

At later times, we observe a slow growth of the total flow as a result of rising groundwater. Within this regime, the groundwater flow increases and the seepage zone slowly grows; consequently, there is an increased area over which an overland flow is generated. These effects cause the river flow, $Q(t)$, to exceed the critical flow $Q_{crit}$. If $Q(t) \ll Q_{crit}$, then the river flow is mostly caused by the early-time mechanism, while if $Q(t)\gg Q_{crit}$, then the late-time mechanism dominates.

It should be noted that, here, we have introduced the intuition of the critical flow in (8.1) as a way to better interpret the numerical results. Shortly in § 8.2, we will justify based on sensitivity analysis of the model that many numerical solutions in the phase space do exhibit this behaviour (saturation to the critical flow). Moreover, in Part 3 of our work, we will derive $Q_{crit}$ in a more rigorous way based on the asymptotic analysis based on a shallow-water approximation. For this case, the analogue to (8.1) will be developed asymptotically.

8.1.2. Experiment (B): case with no initial seepage zone

In experiment (B), there is no initial seepage zone. If the rainfall, $r$, is smaller than a certain value (dependent on soil geometry and properties around the channel), we may observe a slow rise in the groundwater table gradient around $x=0$, leading to an increase in the groundwater flow. If the rainfall is higher than this threshold value (as in the case presented in figure 12), then the gradient of the groundwater table eventually reaches the elevation gradient.

For typical values of rainfall much higher than the threshold value, this initial phase is very short and, in practice, not noticeable in the presented hydrograph. After that moment, a seepage zone starts to grow, giving rise to the overland flow, which slowly increases as the saturation front propagates. This is similar to the late-time behaviour of the first hydrograph. Additionally, we observe a rise in the groundwater flow as a result of the growing pressure head in the groundwater around the stream forced by the rising groundwater table. In the first case, the rise of the groundwater table was taking place far from the channel (relatively to its dimension), and so its effect on the groundwater recharge to the channel is not observable.

8.2. Sensitivity analysis

In order to understand the relations between the described dynamics and model parameters, we conducted a sensitivity analysis. We chose eight physical parameters: catchment width $L_x$, aquifer depth $L_z$, elevation gradient along the hillslope $S_x$, hydraulic conductivity $K_s$, precipitation rates $r$ and $r_0$, Manning's constant $n_s$ and the $\alpha _{MvG}$ parameter. Each parameter was varied within the range of its typical values presented in table 1 following Morawiecki & Trinh (2024), while keeping other parameters constant. In figure 14, we present the peak flow and its components after 24 hours, each as a function of the different parameter values. The critical flow calculated using (8.1) is also shown on the graphs as a dashed line.

Figure 14. Results of the sensitivity analysis, showing the dependence of model parameters on the peak flow (light blue) and initial flow (dark blue). The $y$-axis on each figure represents the flow expressed in $\textrm {m}^2\ \textrm {s}^{-1}$. The peak flow is measured for a rainfall of duration of $t=24$ hours. The critical flow, represented with a dashed line, is defined by (8.1).

Based on this analysis and the investigation of the numerical solutions, the following conclusions can be drawn:

1. (i) The critical flow generated by the precipitation accumulating over the initial seepage zone is a significant component of the peak flow; this description is consistent over the different model parameters.

2. (ii) The size of this seepage zone depends on the difference between (a) the total precipitation in the initial condition, and (b) the total groundwater flow. The former, (a), is a product of the precipitation rate, $r_0$, and the catchment area, $A=L_xL_y$, both of which are positively correlated with the seepage zone size. The latter, (b), following Darcy's law, depends on hydraulic conductivity $K_s$, pressure gradient (dependent on slope $S_x$) and the aquifer depth $L_z$, all of which are negatively correlated with the size of the seepage zone.

3. (iii) The precipitation rate, $r$, has a significant impact on both the critical flow as given by (8.1), and on the further growth of the overland flow. This is because it is responsible for the speed of groundwater rising and for surface water accumulation in the growing seepage zone.

4. (iv) The speed of the seepage zone growth is slower for higher slope, $S_x$, values, since $S_x$ determines how deeply the groundwater table is located beneath the surface and how much rainwater it can absorb before reaching the surface. Also, $\alpha _{MvG}$ has a small effect on the seepage zone growth, since it determines the soil saturation above the groundwater table. A higher $\alpha _{MvG}$ causes the soil saturation to drop faster with height, allowing it to absorb more rainwater before it saturates. A similar effect is observed when varying other MvG model parameters ($\theta _S$, $\theta _R$, $n$). The impact of other model parameters on the hydrograph shape after reaching critical flow is very small.

5. (v) The soil depth, $L_z$, has a significant impact on the groundwater flow only for small values (comparable to the depth of the channel). Increasing $L_z$ above 30 m has little impact on the solution, since the flow at such depths is insignificant.

6. (vi) Manning's constant, $n_s$, seems to have almost no impact on the hydrograph. Its main contribution is in affecting the overland flow speed via Manning's law (3.4), and so it affects the characteristic time scale given by the $\tau _s$ parameter. This time scale, however, is shorter than the duration of the simulated rainfall. The effect of the $n_s$ parameter can be significant if the rainfall duration is shorter than the time required to reach the critical flow (which is dependent on $n_s$). We will derive an analytical expression for this time in Part 3 of our work.

9. Discussion

The central question presented in our work is quite simple: What is the simplest 3-D model of coupled surface–subsurface flow on a hillslope?

Despite the fundamental nature of the above question, we have been surprised at the lack of mathematical and fluid dynamical research on issues of this nature in the literature. As mentioned throughout, we have been strongly motivated by the recent work of Maxwell et al. (2014), who designed benchmark scenarios for the purpose of comparing computational catchment models. Here, our philosophy has been more comprehensive in nature, and we are interested in the analytical and computational properties of the model rather than using it as a means to an end. Our benchmark involves several improvements over those proposed previously, allowing us to replicate hydrographs similar to the ones observed in real-world systems.

This work provides deeper insight into the mathematical structure of coupled surface–subsurface models. We extract and interpret nine key dimensionless parameters. As we show using asymptotic methods, under certain initial and geometric conditions ($S_y\ll S_x$, $L_y\ll L_x$ or $L_y\gg L_x$), the original formulation of the 3-D model can be reduced to a 2-D form. We then numerically investigated the shape and scale of the 3-D features, which subsequently allows us to quantify the error in the 3-D-to-2-D reduction.

Our sensitivity analysis of the key physical parameters reveals several interesting dependencies. As we demonstrate, the peak flows observed during sufficiently long rainfalls are usually caused by two mechanisms. First, there is an early-time rise due to surface water accumulating in the part of the catchment already saturated before the rainfall was initiated. Second, there is a late-time effect due to the slow propagation of the seepage zone. This two-scale behaviour can be rigorously justified based on asymptotic analysis of the governing equations. In the accompanying Part 3 of our work, we study the situation of aquifers with a depth much smaller than the catchment width (the shallow aquifer scenario in figure 2). There, we shall demonstrate that a shallow-water approximation allows us to derive analytical scaling laws for the hydrograph, and hence precise quantification of the peak flows mentioned above.

We note some potential consequences of our benchmark model for future research. The (relative) simplicity of our benchmark, and the clear isolation of properties such as peak flow formation and their parametric dependencies, means that the benchmark can be used in future studies for intermodel comparison. For example, data-based methods, such as conceptual and statistical models, may exhibit a different dependence on catchment properties. Then, by isolating the reasons for such discrepancies, we may better understand the limitations of different classes of models. This potentially leads to the development of more theoretically justified models in the future, which may offer improvements in accuracy over a wider range of scenarios.