On the development and analysis of coupled surface-subsurface models of catchments.

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


Introduction
Since the publication of the Stanford Watershed Model by Crawford and Linsley (1966) a wide range of computational models of catchment-scale hydrology have been developed (Singh and Frevert 2003).Indeed, over two hundred models have been identified in the extensive review by Peel and 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 and Simmons 2012), ParFlow (Kollet and Maxwell 2006), and OpenGeoSys (Kolditz et al. 2012).
However, in contrast to computational studies, there seems to have been more limited † Email address for correspondence: piotr.morawiecki@bath.edu‡ Email address for correspondence: p.trinh@bath.ac.uk 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 and Imhoff 2006).However, further computational power was needed before the first physically-based models were implemented.Notable early examples include TOPMODEL (Kirkby and Beven 1979), MIKE SHE (Abbott et al. 1986b,a), and IHDM (Institute of Hydrology Distributed Model, cf.Beven et al. 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's predictions (usually after earlier calibration) to 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 andBinley 1992 andGupta et al. 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 and Maxwell (2006) used a tilted V-shaped catchment geometry (fig.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 (fig.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.
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  Kollet and 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.evapotranspiration (Kollet et al. 2009), atmosphere (Sulis et al. 2017), biochemistry (Cui et al. 2014), the impact of climate change (Markovich et al. 2016), and the effects of different types of heterogeneities, e.g. the heterogeneity of the land surface (Rihani et al. 2015), soil properties (Meyerhoff and Maxwell 2011), and even flow through fractures (Sweetenham et al. 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 (fig.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, though, 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 nondimensional parameters characterising surface-subsurface hydrological processes.We highlight some prior works that have used nondimensionalisation 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 et al. (1990); Warrick and Hussen (1993); 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 nondimensionalisation plays a prominent role for the case of coupled surface-subsurface models, was performed by Sivapalan et al. (1987), and focuses on the TOPMODEL scheme of Kirkby and Beven (1979).A similar study was performed by Calver and Wood (1991) for the IHDM model (Beven et al. 1987).In particular, Calver and 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.

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 and Trinh 2022) 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 to the existing literature, there are three novel elements in our study: (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.(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 2D 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).(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 and Nieber 1977), or a sudden drawdown drainage when a rapid change of water level occurs at the outlet (Sanford et al. 1993).
We start by formulating a three-dimensional benchmark scenario in section 2, which is non-dimensionalised in section 3.In section 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 section 6, this model simplification is numerically assessed in section 7. The impact of each parameter in the resulting twodimensional model is summarised in section 8, which is followed by the discussion in section 9.
Symbols.There are many symbols in this work.For ease of reference, we provide a list of symbols in table 2 and table 3 in appendix A.

Formulation of a simplified three-dimensional catchment model
In this section, we formulate a simplified catchment model, inspired by the infiltrationexcess, 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 fig. 2.
(a) The V-shaped catchment.This scenario, shown in fig.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 × L y × 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 section 5, we demonstrate that under certain conditions, the scenario reduces to largely two-dimensional dynamics along the hillslope.(b) The deep aquifer.This scenario, shown in fig.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 saturationexcess scenarios discussed in Maxwell et al. (2014), the channel is assumed to have a rectangular xz cross-section with width, w, and depth, d.(c) The shallow aquifer.This scenario, shown in fig.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 ≪ L x .We analyse this scenario in Part 3. 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 (x, ŷ, ẑ) are directed along the hillslope edges.Hence, x is directed along the hillslope (x = 0 representing the location of the channel), ŷ along the channel (ŷ = 0 representing the location of the outlet), and ẑ vertically (ẑ = 0 representing the bottom of the aquifer).After the coordinate transformation, the entire catchment can be represented as a cuboid of dimensions L x × L y × L z .The following coordinate transformation is used: We introduced L x to represent the catchment width along the x direction given as The land surface in this geometry corresponds to 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 and 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 and Germann (2013), including these effects in the model may significantly affect the timescale of infiltration.

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, which for a typical UK catchment is ϵ ∈ [0.13, 0.25] †.We also have the aspect ratio between the catchment height and the catchment dimension along the river, which for a typical UK catchment is β zy ∈ [0.0007, 0.025].Finally, we have the aspect ratio between the catchment height and the catchment length along the hillslope: which for a typical UK catchment is β zx ∈ [0.1, 2.1].Note that as β zy /β zx → 0, we get long catchments with a width much shorter than their length, while for β zy /β zx → ∞, 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 fig. 3. Here, we draw lines of constant topographic elevation on a projection of the catchment onto z = 0. Note that, for example in fig.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 fig.3(c), we may expect to observe a significant flow component parallel to the river direction.
Increasing Sy (increasing ) 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.† 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).Maxwell et al. (2004) Here, we briefly outline how the scenarios introduced above relate to the scenarios presented in the benchmark analysis of Maxwell et al. (2014).

Relationship to
In sections 4.1 and 4.2, Maxwell et al. (2014) introduce two scenarios called the infiltration-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 fig.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 fig.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 fig. 4 from Part 1).
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 section 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, ϕ = 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).

Governing equations (dimensional)
We begin with the dimensional model.As introduced in sec. 2 of Part 1, we consider three types of flow: the subsurface flow (the 3D Richards equation), the overland flow (the 2D Saint Venant equations), and the channel flow (1D 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 and Ogden (2017); Schaake Jr (1975) and references therein.

3D Richards equation for the subsurface flow
The subsurface flow 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 three-dimensional Richards equation (see e.g.Dogan and Motz (2005) and Weill et al. (2009)), which is given by: dθ dh g ∂h g ∂t = ∇ • q g , where q g = K s K r (h g )∇(h g + z). (3.1) Here, ∇ = ( ∂ ∂x , ∂ ∂y , ∂ ∂z ) is a standard nabla operator, K s > 0 is the saturated soil conductivity and dθ(hg) dhg is the so-called specific moisture capacity.We assume that the volumetric water content θ(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): 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 fig.6a 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 α MvG , θ r , θ s , n, and m = 1 − 1 n , which depend on the soil properties.The residual water content θ r and saturated water content θ s represent the lowest and the highest water content, respectively.The α MvG parameter in m −1 represents the scaling factor for the pressure head h g [m].The n coefficient describes the pore sizes distribution.

2D 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 and Kavvas (1994) and Liu et al. (2004) this flow is described using the two-dimensional Saint Venant equations that govern the overland water height, h s (x, y, t).Following the discussion in sec.2.2 of Part 1, we consider mass and momentum conservation.Firstly, the continuity equation is given by 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, 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 where n s is an empirically determined value known as Manning's coefficient, and describes the overland surface roughness; 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 S f .In general, the friction slope, S f , is given by momentum conservation [cf.eqn (2.7) in Part 1].However, in computational integrated catchment models, a kinematic approximation is often used which neglects all effects on S f other than gravity.This approximation is used in e.g.Parflow (Maxwell et al. 2009), though there are others such as e.g.MIKE SHE that implement a more complete, diffusive approximation (by DHI 2017).In the case of the kinematic approximation, where S 0 = −∇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 x-direction).As Vieira (1983) argues, this approximation may give inaccurate predictions when the system is close to reaching a steady state.

1D 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 ŷ-axis.Following Vieira (1983) and Chaudhry (2007), the channel flow is modelled as a one-dimensional Saint Venant equation that governs the channel water height, z = h c (ŷ, t), and is given by: where w(h c , 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: 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 river f is the friction slope, which under kinematic approximation is equal to the elevation gradient along the river: 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 fig.2, where there is a rectangular channel, this involves setting the area A = wh c and P = w + 2h 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 ≈ w, and hence ignore the friction effects of the channel side walls.In this case, Manning's equation (3.7) becomes The later simulations will thus involve the solution of the Saint Venant equation (3.6) with the shallow Manning's 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.

Boundary conditions
The domain consists of four types of boundaries: (i) the catchment boundary, Γ B , both for the surface and subsurface part of the domain, including the bedrock constraining the aquifer from the bottom; (ii) the land surface Γ s ; (iii) the river bank, Γ R ; (iv) the river outlet, Γ O ; and (v) the river inlet Γ I (see fig. 4).(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 Γ B , we set no-flow conditions for both subsurface and surface flow: Alternatively, one could introduce a free-flow condition for the groundwater flow, q g • 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.(ii) Next, on the land surface, Γ s , continuity of pressure and flow between the groundwater and surface water yields: This first condition imposes continuity of pressure only if the groundwater reaches Γ s , while the second imposes the condition of rain infiltration, I.

(iii)
On the river bank, Γ 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 : At the inlet, located at the upstream end of the river, Γ 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 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 the 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 inlows 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 Γ R (for overland flow) and at the catchment outlet Γ 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 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).

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

Governing equations (nondimensional)
4.1.Nondimensionalisation The governing equations for subsurface, surface, and channel flow presented in section 2 are now written in tilted coordinates (x, ŷ, ẑ) [cf.(2.1)] 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 nondimensionalise these equations.The following scalings are used: 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: 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 eqns (B 1), (B 2), and (B 3) in appendix B.1.Additionally, the non-dimensional terms in (4.1) have straightforward physical interpretations.The timescale, 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 section 6.2.
It should be noted that even though t 0 is a characteristic time of the vertical flow through the soil, other timescales are present.For example, we shall observe typically shorter timescales for the overland flow, and much longer timescales for the horizontal flow through the soil.Further discussion of the separation of timescales appears in Part 3 of our work.

Summary of governing equations and parameters
We collect the nondimensional governing equations from appendix B.2.To review, our hydrological problems in the three-dimensional geometry consist of solving three timedependent partial differential equations for three unknowns: (i) a 3D Richards equation for the subsurface flow (4.2a); (ii) a 2D Saint Venant equation for the overland flow (4.2b); and (iii) a 1D Saint Venant equation for the channel flow (4.2c).In the tilted frame, these are respectively where the subsurface equations involve operators definitions: Expressions for θ(h g ) and K r (h g ) are provided in appendix B.2.Each partial differential equation in (4.2) is solved subject to boundary conditions posed on the domain boundaries given by (B 8a)-(B 8d).

Model reduction to a two-dimensional model
In section 2, we formulated a general three-dimensional catchment model.The purpose of this section is to discuss the nondimensionalisation 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: (i) Considering either small river slope (S y ≪ S x ), short (L y ≪ L z ), or long catchment (L y ≫ L z ) approximations together with the necessary low channel limit (L c ≪ L z ), we may reduce the general three-dimensional governing equations for h g and h s to a two-dimensional form neglecting the flow along the y-axis.(ii) In addition, in the case of the shallow aquifer scenario (L z ≪ L x ), we may apply a shallow-water approximation to further reduce the 2D hillslope model to a 1D 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 ≪ L z The three-dimensional model in (x, y, z) can be formally approximated by a twodimensional model in (x, z) if the subsurface profile, h g (x, y, z, t) ∼ h g,0 (x, z, t), and the surface profile, h s (x, y, t) ∼ 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, Γ 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 nondimensional values for h g , h c ,, and z, (3.10c) yields 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 to the typical variation of the pressure head.In the limit of L c /L z → 0, we see that the subsurface boundary condition is which is no longer y-dependent.

5.
2. An asymptotic expansion for small river slopes, in ϵ = 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 section 2.1, ϵ = S y /S x , typically has small values (half of UK catchments have ϵ between 0.13 and 0.25).Here we shall demonstrate that when ϵ ≪ 1 (equivalent to S y ≪ 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, (B 4) and (B 5), in a simpler form highlighting its structure: where the nonlinear operators, N i , for i = 1, 2, 3 are defined in (B 13) in appendix B.2.Note that these operators are dependent on h g and h s , and independent of ϵ and β zy , which are the only dimensionless parameters involving S y and L y .
When ϵ = 0, we can verify that the solutions are independent of ŷ, i.e. they can be written as h g (x, ŷ, ẑ) = h g,0 (x, ẑ) and h s (x, ŷ) = h s,0 (x).This is caused by the combination of three facts: operators N 2 and N 3 applied to a function independent of ŷ become 0; and (iii) the no-flux boundary condition at ŷ = 0, 1 in (B 8a) is then: (5.4) Hence, for ϵ = 0, the above boundary condition is satisfied by h g = h g,0 (x, ẑ).(iv) From eqn (5.3b)only R eff and I can be ŷ-dependent terms, but in the considered scenario R eff is constant, and I(x, ŷ) is ŷ-independent as long as h g is.
Essentially, ϵ = 0 is associated with a zero gradient along the river, i.e. there is no forcing flow in the ŷ-direction, and the domain becomes transitionally symmetric in that direction.
There is an important consideration in the formal limit as S y ∝ ϵ → 0. In this limit, holding other parameters fixed, the L c defined in (4.1) tends to ∞, and so the L c ≪ 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 increases in order to maintain a significant channel flow [cf.Manning's law (3.9)].Therefore, we would expect for two-dimensional dynamics to dominate when (5.5) Above, the expression on the left-hand side is obtained from the definition of L y from (4.1), 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 • 10 −11 ), and so only the second approximation in (5.5) needs to be considered.
5.3.Asymptotic expansions for short (β zy ≫ 1) and long (β zy ≪ 1) catchments There are additional limits that allow us to reduce the three-dimensional problem into simpler two-dimensional formulations at the leading order, and these involve the nondimensional geometrical parameter (5.6)For instance, in the limit as β zy → ∞, the three-dimensional 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 ŷ = 0 and ŷ = 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 dh g,0 / dŷ = 0 condition at ŷ = 0, 1, which is automatically satisfied for a ŷ-independent solution.As argued in the previous section, we also conclude that such a ŷ-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 β −1 zy : (5.7b) The last interesting limit we discuss is β zy → 0, which corresponds to the situation of an asymptotically long river.Similarly, the two-dimensional 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 ŷ = 0 and ŷ = 1 (see fig. 5).Consequently, h g,0 and h s,0 are understood to represent the 'outer' asymptotic solutions, valid for 0 < ŷ < 1.
Without loss of generality, let us consider the boundary layer near ŷ = 0. We rescale ŷ = δ ŷ′ where δ(β zy ) is a characteristic size of the boundary layer.After applying this transformation, the governing equations (5.3) become: (5.8b) The three-dimensional effects given by N 2 and N 3 become significant when the characteristic size of the boundary layer is of the order O(β zy ).Hence, we conclude that the thickness of the boundary layer is O(1/β zy ).So, if we consider the β zy → 0 Figure 5: This graphic shows the steady-state depth of the groundwater table, according to the 3D model.We note that the solution is mostly ŷ-independent, except for two apparent boundary layers around ŷ = 0 and ŷ = 1; near these points, the groundwater table aligns with the lines of constant elevation.The figure is generated using the solver described in section 6 using the parameter values given in table 1, with except for two values: we use L y = 13680 m (β zy = 0.05) and S y = 0.0075 (ϵ = 0.1); as a result, this graphic matches an inset in fig.9.The boundary layer thicknesses, δ 1 and δ 2 , tend to zero as β zy → ∞. limit, the solution for the problem becomes two-dimensional except for an infinitely thin boundary layer around the boundaries.

Summary of the two-dimensional model
To summarise, we considered three approximations for small river slope (ϵ ≪ 1), short catchments (β zy ≫ 1), and long catchments (β zy ≪ 1).We showed that under any of these approximations, the model can be approximately represented in the following two-dimensional form: Both the groundwater and overland flows reaching the channel contribute to the river flow in the ŷ-direction governed by (B 6).Thus, 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 (ŷ = 0, t).However, note that after the reduction to a two-dimensional model, equations (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 2D model by solving equations (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 2D model in section 8, and in Part 3 of our work.
The above system of PDEs forms a model describing surface and subsurface flow in a 2D hillslope cross-section, as presented in fig.6.The boundaries are now one-dimensional, but the boundary conditions are the same as in the three-dimensional model, as given by (B 8a)-(B 8d).As before, for the initial condition, we consider the steady state of the above system for R eff = R 0 .In section 7, we will explore the accuracy of this approximation numerically, investigating the size of three-dimensional features of the full solution and their behaviour in limits formulated above.A two remarks are in order: (i) Firstly, there are some remaining terms in the dimensionless governing equations, which are small; however, they should not be neglected.For example, coefficient, is small, which means that the characteristic timescale of accumulation of surface water is much shorter than the characteristic timescale of vertical subsurface flow.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 ≈ 6.8 • 10 7 s ≈ 2 years), we are often interested in the short-time behaviour, and therefore, this term should not be neglected.
(ii) Secondly, the β zx term is also small in the shallow aquifer scenario compared to the leading term representing the flow in the ẑ-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 dθ dh = 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 β 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.

Numerical methodology
In order to validate the reduction of the 3D model to the 2D 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 section 3. To summarise, these are: subject to boundary conditions (3.10a)-(3.10d).
Our numerical implementation has two applications in this study.Firstly, in section 7 we use the numerical scheme based on a discrete version of this equations to verify our reductions to the 2D problem introduced in section 5. Secondly, in section 8.2 we numerically analyse features of a benchmark scenario in a reduced two-dimensional 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).

Model discretisation
The implementation of the 3D and 2D models is based on the finite volume method.The entire hillslope is divided into N x × N y × N z cells.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 fig.7a.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 fig. 7b).
Depending on the type of simulation, we handle the channel differently.In the case of two-dimensional simulations, we will assume the water level in the channel to be equal to the channel depth (unless stated otherwise).In the case of three-dimensional 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's 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 Example of mesh refinement.The size of edges is given by the geometric series with ratios µ and ν.
h s at the surface as an independent variable (the same applies to h c in three-dimensional simulations).One challenge is the significant difference in the timescales between the overland and groundwater flow (the ratio of which is quantified with the dimensionless parameter τ s ≈ 2.8 • 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 (B 4) for each cell is used: where Few remarks should be made here: (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, ∆t is the time step duration, h t i and h t+1 i is the pressure head (h g ) in cell i at time step t (previous one) and t + 1 (current one) respectively, h ′t i is the pressure head computed in the previous iteration of the implicit scheme, and θ is the saturation given by the Mualem-Van Genuchten model (B 7).(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, η i,j is the angle between this face and the line joining these cells' centroids, r i→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, β = (β 2 zx , β 2 zy , 1) is a vector describing the anisotropy coming from different scaling factors in the dimensionless model.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 Mualem-Van Genuchten model (B 7b) for h = h ′t+1 u(i,j) , where u i,j = i if the flow is going from cell i to j, and u i,j = j otherwise.(iii) The change of both θ 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 t+1 i satisfying the continuity condition.Equation (6.2) is linear in h ′t+1 i for all i, and therefore, it can be solved using standard methods for linear algebraic equations.(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 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 (∆x∆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 3D simulations) is computed for each cell as 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 discretized form of the 2D Saint Venant equation (B 5): 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 3D model, we calculate the river flow at the outlet Q(t) = q c (ŷ = 0, t), and in the case of the 2D 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 PAWS model (and other coupled surface-subsurface models).The results of this comparison are presented in appendix D.

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 parameter default value parameter range 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 trib y estimated from Part 1, which represents the typical distance between the main tributaries of the river.
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 trib y = 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 and 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 fig.8. Initially, the system remains in a steady state, in which the pressure head h g increases with depth following an approximately hydrostatic profile (fig.8a).The interesting dynamics responsible for generating the flow occurs near the surface (ẑ = 1).
The surface water is present near the channel and extends further away from it for lower ŷ values (fig.8d-f).However, around ŷ = 1, we do not observe surface water at all.This is caused by a non-zero gradient S y along the ŷ 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 fig. 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 fig.8e).Increased overland flow, leads to a rapid rise of the channel water and resulting outflow from the catchment within the first hour (fig.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 fig.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 fig.8f), increasing the area from which overland flow reaches the river.
An essential observation for this time-dependent solution is that the characteristic timescale of overland flow dynamics is much shorter than the characteristic time of groundwater flow (their ratio is given by the dimensionless parameter τ s ≈ 2.8 • 10 −4 ).This timescale separation is reflected in the shape of the hydrograph in fig.8c, which shows the dependence between river flow at the outlet Q(t) = q c (ŷ, 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 (section 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→∞ 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.

Verification of 3D to 2D reduction
The time-dependent simulations presented in the previous section demonstrate some of the three-dimensional 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 ŷ direction, namely the catchment length L y and the slope parallel to the channel S y .Alternatively, in terms of the nondimensional quantities, this corresponds to β zy and ϵ.

Three-dimensional features of steady state solution
In order to develop a better understanding of the three-dimensional 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 fig.9.In all cases we observe that the solution becomes less ŷ-dependent as ϵ → 0, as expected from section 5.However, the dependence on β zy = L z /L y is more complex.The phase space can be divided into three regions: When L y ≪ L x ϵ , the lines of constant elevation are approximately perpendicular to the hillslope (e.g.ϵ = 0.1, L y = 100).As shown in section 5, in such a case, the leading-order (two-dimensional) solution of the governing equations for small ϵ 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.(ii) When L y ≫ L x ϵ , the lines of constant elevation are approximately parallel to the hillslope (e.g.ϵ = 0.1, L y = 10 6 ).In such a case, the leading-order (2D) solution for the governing equations does not satisfy the flow boundary condition at ŷ = 0 and ŷ = 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 δ decreases inversely proportional to L y (see fig. 10), which is consistent with the theoretical scaling derived in section 5.3.(iii) In the intermediate region, when L y = O L x ϵ and δ = O(1) (e.g.ϵ = 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 two-dimensional approximation unless ϵ is small enough.

Analysis of 3D to 2D reduction error
Here, we follow the qualitative discussion from the previous section by quantifying the difference between the full solution and its two-dimensional approximation.
In section 5, we argued that the solution for the three-dimensional problem can be represented as h g (x, ŷ, ẑ, t) = h g,0 (x, ŷ, t) + ϵh g,1 (x, ŷ, ẑ, t) + O(ϵ 2 ), where h g,0 is a twodimensional solution for ϵ = 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 ϵ is presented in fig.11.It confirms the asymptotic analysis from section 5 and the qualitative observations from section 7.1.Firstly, it confirms that the error of the two-dimensional approximation increases proportionally to ϵ for large values of ϵ.However for small values of ϵ, the error increases, because the effect of y-dependent water height at the channel is no longer negligible (see section 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 (2D solution is satisfied everywhere apart from a thin boundary layer at ŷ = 0 and ŷ = 1), but the error is the highest for intermediate values (here around L y = 3000).

Impact of physical parameters on the 2D model
Following section 5, the inflow to the river in our benchmark scenario for S y ≪ S x can be approximated by a two-dimensional model.In this section, we use the numerical procedure The boundary thickness was measured based on the groundwater depth profile along x ≈ 0.46.The boundary was defined as ŷ, for which the groundwater depth H(0.46, ŷ) is further than ±5% from the groundwater depth evaluated in the middle of the domain, H(0.46, ŷ = 0.5).For small β zy values, the boundary width follows δ ∝ β zy scaling.As β zy increases, δ reaches 0.5, for which the boundary condition affects effectively the entire domain.The mean absolute difference between the full three-dimensional solution and its two-dimensional approximation for small ϵ as a function of ϵ and β zy .described in section 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.

Structure of typical hydrographs
In this section, we examine, in more detail, the hydrographs that correspond to two protypical 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 section 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 • 10 −8 m/s; and Experiment (B) with r 0 = 2 • 10 −9 m/s.The two corresponding hydrographs, Q(t) vs t, are presented in the top-left inserts of fig. 12 and fig.13.For each hydrograph, solutions are presented at four times and given in the 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.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 fig.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 fig.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 three-dimensional case presented in fig.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 fig.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 The above quantity we shall refer to as the critical flow.Here, r 0 A represents the initial flow, and A = L x L y is the catchment's area.In fact, the dashed line plotted in the hydrograph of figure12 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 the flow is mostly caused by the early-time mechanism, while if Q(t) ≫ 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 section 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.

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

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 α MvG parameter.Each parameter was varied within the range of its typical values presented in table 1 following Morawiecki and Trinh (2022), while keeping other parameters constant.In fig.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.
Based on this analysis and the investigation of the numerical solutions, the following conclusions can be drawn: 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.(ii) The size of this seepage zone depends on the difference between (i) the total precipitation in the initial condition, and (ii) the total groundwater flow.The former, (i), is a product of the precipitation rate, r 0 , and the catchment area, A = L x L y , both of which are positively correlated with the seepage zone size.The latter, (ii), 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.(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.(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, α MvG has a small effect on the seepage zone growth, since it determines the soil saturation above the groundwater table.A higher α 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 Mualem-van Genuchten model parameters (θ S , θ R , n).The impact of other model parameters on the hydrograph shape after reaching critical flow is very small.

(v)
The soil depth, L z , has a significant impact on the groundwater flow only for small values (comparable with the depth of the channel).Increasing L z above 30m has little impact on the solution, since the flow at such depths is insignificant.(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 the Manning's law (3.4), and so it affects the characteristic timescale given by the τ s parameter.This timescale, 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.

Discussion
The central question presented in our work is quite simple: What is the simplest threedimensional 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 surfacesubsurface models.We extract and interpret nine key dimensionless parameters.As we show using asymptotic methods, under certain initial and geometric conditions (S y ≪ S x , L y ≪ L x or L y ≫ L x ), the original formulation of the 3D model can be reduced to a 2D form.We then numerically investigated the shape and scale of the three-dimensional features, which subsequently allows us to quantify the error in the 3D-to-2D 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 fig.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.

Appendix A. List of symbols
For convenience, we provide a list of symbols in tables 2 and 3. terms of asymptotic expansion of h s q g , q s , q c groundwater, overland, and channel flow v s , v c velocity of overland and channel flow S f friction slope for the overland flow S river face area between cell i and j r i→j vector from the centroid of cell i to the centroid of cell j β vector of β parameters, β = (β 2 zx , β 2 zy , 1) K ′ i,j hydraulic conductivity of the face between cell i and j u i,j function returning the index of the uplift cell (i or j) f t i,j water volume in surface cell (i, j) divided by its base area ∆x, ∆y extent of the cell in the x and y direction The Saint Venant equation (3.3), together with (3.4) in transformed coordinates x and ŷ becomes: where R eff = R − ET is the effective precipitation.The channel flow is given by equations (3.6) and (3.9) combined, which for our simplified catchment give the last governing equation: All boundary conditions in the dimensional form are listed in section 3.4.

B.2. Dimensionless form
Now, we rewrite eqs (B 1), (B 2), (B 3) using the dimensionless quantities introduced in section 4.1.Here and henceforth, we shall drop the primes, and assume that all subsequent quantities are dimensionless.The dimensionless governing equations are as follows.First, the 3D Richards equation for pressure head, h g (x, ŷ, ẑ): . (B 4) The 2D Saint Venant equation for overland water height h s (x, ŷ): The definition of dimensionless parameters (β zx , β zy , τ s , τ c , γ), their interpretation, and estimated values are presented in appendix C. Numerical values under the equations represent the typical order of magnitude of parameters multiplying the given term.However, note that terms marked with " †" symbol include the ŷ-derivative of the solution, which, as was discussed in section 5, is ŷ-independent in the leading order.The effect of the relative size of these terms is much smaller (by approximately one order of magnitude) than indicated by the provided values of prefactors.
In the above equations, the dimensionless θ(h) and K r (h) functions are given by: dθ(h) dh = mn(θs−θr) h where α = α MvG L z is a dimensionless MvG α parameter.Finally, as divided into the enumerations of section 3.4, the non-dimensional boundary conditions are now as follows: (i) On the catchment boundary, Γ B : On the land surface, Γ s : (iii) At the channel, Γ R : Finally, in the river inlet, Γ I : In the above boundary conditions, we have introduced three new dimensionless parameters: ρ = r Ks , λ s = Ls Lz , and λ c = Lc Lz .The river inflow terms in (3.11) and of the soil.In this case, only a part of the rainwater infiltrates through the ground, while the remaining part forms a so-called Horton overland flow.
All the necessary model parameters are presented in the aforemention publications, so no calibration is required.However, information about initial and boundary conditions is missing from the works.We used the same boundary conditions as presented in section 2, while for the initial condition, we assumed a constant depth of the groundwater table, with pressure head h decreasing linearly with depth z (this corresponds to no initial vertical flow through the soil).
We compared the results obtained using the finite volume solver described in section 6 with the results obtained using ParFlow presented by Sulis et al. (2010).As before, we used an image processing tool to extract the data from the graphs presented in these publications.
Fig. 16 demonstrates that the solver very accurately reproduces the results from the original paper in both scenarios for a dense computational mesh (∆z = 0.0125 m).Fig. 17 additionally shows that the solver also produces almost identical output for lower resolution (∆z = 0.1 m, ∆z = 0.2 m), which demonstrates the similarity of our discretisation and numerical artefacts.Since Maxwell et al. (2014) showed that ParFlow results are consistent with other currently used physical catchment models, we conclude that our solver properly represents all their assumptions within the framework of the considered simple scenario.

Figure 1 :
Figure1: Illustration of the idealised catchment geometries developed in the works ofKollet and Maxwell (2006) andGilbert 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.

Figure 2 :
Figure 2: Simplified catchment geometry in the considered scenarios (not to scale).

Figure 4 :
Figure 4: Boundaries defined for the V-shaped tilted catchment.

Figure 6 :
Figure6: An illustration of a two-dimensional hillslope geometry in Cartesian coordinates (left) and in the tilted coordinate system (right).These represent a cross-section (in the xz-plane) of the original V-shaped catchment.

Figure 7 :
Figure 7: (a) Discretisation of the 3D catchment, representing the V-shaped catchment scenario.In the case of the 2D 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 µ and ν.

Figure 8 :
Figure 8: The illustration summarises the key properties of the 3D solution obtained for the default values of parameters from table 1. Figure (a) depicts the initial condition (steady state for mean rainfall r 0 = 2.95 • 10 −8 m/s).During the subsequent rainfall r = 2.36 • 10 −7 m/s, 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 ŷ 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.

Figure 9 :
Figure 9: Groundwater table depth in steady states obtained for varying catchment length L y = β −1 zy L z and slope S y = ϵS x .Dashed lines represent lines of constant elevation.The entry with ϵ = 0.1 and β zy = 0 is the same figure as presented in fig. 5.

Figure 10 :
Figure 10: Boundary layer thickness at ŷ = 0 (a) and ŷ = 1 as a function of β zy .The boundary thickness was measured based on the groundwater depth profile along x ≈ 0.46.The boundary was defined as ŷ, for which the groundwater depth H(0.46, ŷ) is further than ±5% from the groundwater depth evaluated in the middle of the domain, H(0.46, ŷ = 0.5).For small β zy values, the boundary width follows δ ∝ β zy scaling.As β zy increases, δ reaches 0.5, for which the boundary condition affects effectively the entire domain.

Figure
Figure11: The mean absolute difference between the full three-dimensional solution and its two-dimensional approximation for small ϵ as a function of ϵ and β zy .

Figure 12 :
Figure 12: Numerical solution of 2D model for r 0 = 2.95 • 10 ms −1parameters were set to the default values presented in table 1.

Figure 13 :
Figure 13: Numerical solution of model for r 0 = 2 • 10 −9 ms −1 without an initial seepage zone.All other parameters were set to the default values presented in table 1.

Figure 14 :
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 [m 2 /s].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 equation (8.1).
rate r eff effective rainfall (defined as R − ET ) n s , n cManning's n coefficient for surface and channel g gravitational acceleration q in total surface and subsurface flow to the channel A area of channel cross-section P channel wetted perimeter Catchment L x , L x, L y , L z catchment/hillslope dimension along x, x, y, and z geometry H surf (x, ŷ) elevation of the land surface S 0 = −∇H surf elevation gradient (slope) S x , S y slope measured along x and y direction ϕ angle between the direction of the steepest descent and the x L c characteristic overland and channel water height v x s,0 , v ŷ s,0 , v c,0 characteristic scale of overland and channel flow velocity Table 2: First list of symbols group symbol description Dimensionless α dimensionless α MvG parameter parameters β zx , β zy aspect ratio of cross-section along the hillslope and channel ϵ ratio of S y to S x slope γ aspect ratio of the stream's cross-section λ s , λ c ratio between surface/channel characteristic water height and aquifer thickness L z ρ dimensionless rainfall defined as ρ = r Ks τ s , τ c ratio of overland/channel flow and groundwater timescale Numerical N x , N y , N z number of mesh cells along the x, y, and z axes method ∆t time step duration V i volume of cell i S i,j 1D Saint Venant equation for channel water height h s (ŷ):
The definitions of the dimensionless parameters is provided in appendix C. MvG model parameters quantifying pore size distribution Overland and h s , h c water height on the land surface and in the channel channel flow h s,0 , h s,1

Table 3 :
Second list of symbols Appendix B. Governing equations is tilted coordinates B.1.Dimensional formHere we write down the governing equations introduced in section 2 in (x, ŷ, ẑ) coordinates as given by transformation (2.1).