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. (Reference Maxwell2014)

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

In §§ 4.1 and 4.2, Maxwell et al. (Reference Maxwell2014) 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. (Reference Maxwell2014). Based on our simulations, we observed that if one were to use the parameter values given by Maxwell et al. (Reference Maxwell2014), this would lead to unrealistic steady states, where the seepage covers almost the entire catchment (even for relatively low levels of mean precipitation).

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 Reference Dogan and Motz2005; Weill, Mouche & Patin Reference Weill, Mouche and Patin2009), 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 Reference Van Genuchten1980)

(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 (Reference Tayfur and Kavvas1994) and Liu et al. (Reference Liu, Chen, Li and Singh2004) 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. Reference Maxwell, Kollet, Smith, Woodward, Falgout, Ferguson, Baldwin, Bosl, Hornung and Ashby2009), 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 (Reference Daluz Vieira1983) 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 (Reference Daluz Vieira1983) and Chaudhry (Reference Chaudhry2007), 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. (Reference Maxwell2014). 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 Reference Strahler1957), 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. (Reference Maxwell2014), 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.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.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.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. (Reference Maxwell2014), 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. (Reference Sulis, Meyerhoff, Paniconi, Maxwell, Putti and Kollet2010) obtained for a hillslope using the ParFlow integrated catchment model, and by Maxwell et al. (Reference Maxwell2014) 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 Reference Dietrich and Dunne1993), or a first-/higher-order stream, forming the base element of more complex drainage networks (Strahler Reference Strahler1957). 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.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.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.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 (Reference Morawiecki and Trinh2024), 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.