2.1. Subsurface flow

The subsurface flow that governs the pressure head, $h_g(\boldsymbol {x}, t)$, is commonly modelled using a 3-D Richards equation (Dogan & Motz Reference Dogan and Motz2005; Weill, Mouche & Patin Reference Weill, Mouche and Patin2009),

(2.1)\begin{equation} C(h_g)\frac{\partial{h_g}}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\cdot} [K_s K_r(h_g)\nabla(h_g+z)]. \end{equation}

Here, $K_s > 0$ is the saturated soil conductivity and $C(h)={\mathrm {d}\theta (h)}/{\mathrm {d}h}$ is the so-called specific moisture capacity. The pressure head is equal to zero ($h_g=0$) at the surface of the groundwater table, which separates the fully and partially saturated region of the soil. We assume that the volumetric water content, $\theta (h)$, and relative hydraulic conductivity, $K_r(h)$, are given by the Mualem–van Genuchten (MvG) model (Van Genuchten Reference Van Genuchten1980),

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

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

In essence, the MvG model describes the key hydraulic properties of the soil, such as hydraulic conductivity and saturation, as nonlinear functions of the pressure head, $h$, as well as other parameters $\theta _r$, $\theta _s$, $\alpha _{MvG}$, $n$ and $m=1-{1}/{n}$. The residual water content, $\theta _r$, and saturated water content, $\theta _s$, represent the lowest and the highest water content, respectively. The parameter $\alpha _{MvG}$ $[\mathrm {m}^{-1}]$ is a scaling factor for pressure head $h$ [m]. Finally, the coefficient, $n$, characterises the distribution of pore sizes. Since the MvG model parameters describe the soil/rock properties at a given location, in general, they are functions of the spatial coordinates $(x,y,z)$.

2.2. Overland flow

If rainfall exceeds the infiltration capacity, water can accumulate on the surface and form an overland flow. Typically, following e.g. Chow & Ben-Zvi (Reference Chow and Ben-Zvi1973), Tayfur & Kavvas (Reference Tayfur and Kavvas1994) and Liu et al. (Reference Liu, Chen, Li and Singh2004), this flow is described by the 2-D Saint-Venant equations, which govern the overland water height, $z = h_s(x, y, t)$. The first equation is the continuity equation for the overland flow, written as

(2.4)\begin{equation} \frac{\partial{h_s}}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}+R-I-E, \end{equation}

where the source term includes the precipitation rate $R = R(x,y,t)$, the infiltration rate $I = I(x,y,t)$ and the evapotranspiration rate $E = E(x,y,t)$. The flow $\boldsymbol {q}$ can be expressed as $\boldsymbol {q}=h_s\boldsymbol {u}_s$, where $\boldsymbol {u}_s = \boldsymbol {u}_s(x,y,t)$ is the mean flow speed.

Therefore, (2.4) is a continuity equation with two unknowns, $h_s$ and $u_s$ (or $q_s$). The second equation required to close this system is provided by momentum conservation. In the general form of the Saint-Venant equations, the following equation is used:

(2.5)\begin{equation} \frac{\partial{\boldsymbol{u}_s}}{\partial t}+\boldsymbol{u}_s\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_s+g\boldsymbol{\nabla} h+g(\boldsymbol{S}_{\boldsymbol{f}}-\boldsymbol{S}_{\boldsymbol{0}})=0, \end{equation}

where $g$ is the gravitational acceleration, $\boldsymbol {S}_{\boldsymbol {0}}$ is the elevation gradient (bed slope) and $\boldsymbol {S}_{\boldsymbol {f}}$ is the friction slope, i.e. the rate at which energy is lost along the $x$ and $y$ directions. One can model $\boldsymbol {S}_{\boldsymbol {f}}$ by specifying the shear stress between the overland flow and the surface. However, in the computational hydrological models, rather than solving (2.5) directly, the flux, $\boldsymbol {q}=h\boldsymbol {u}$ is often obtained using an empirical relationship known as Manning's equation (see Shaw et al. (Reference Shaw, Beven, Chappell and Lamb2010, chap. 14.3) for more details). This equation, originally formulated to describe a turbulent channel flow, is also commonly applied to the two-dimensional (2-D) turbulent overland flow over a rough terrain. In the vector form, this equation can be written as

(2.6)\begin{equation} \boldsymbol{q}(\boldsymbol{x}, h_s) = \frac{1}{n_s}\frac{\boldsymbol{S}_{\boldsymbol{f}}}{\sqrt{\|S_f\|}}h_s^{5/3}, \end{equation}

where $n_s$ is an empirically determined value known as Manning's coefficient and describes the land surface roughness. Following (2.5), the friction slope $\mathrm {S_f}$ is expressed as

(2.7)\begin{equation} \boldsymbol{S}_{\boldsymbol{f}}=\boldsymbol{S}_{\boldsymbol{0}}-\boldsymbol{\nabla}h_s -\frac{1}{2g}\boldsymbol{\nabla}{v^2}-\frac{1}{g}\frac{\partial{\boldsymbol{v}}}{\partial t}. \end{equation}

Often in practice, the last two terms of (2.7) can be neglected; this forms the so-called diffusion approximation. An additional common simplification is the kinematic approximation, in which the second term in (2.7) is also neglected, giving

(2.8)\begin{equation} \boldsymbol{S}_{\boldsymbol{f}}=\boldsymbol{S}_{\boldsymbol{0}}. \end{equation}

Note that in general, the gradient varies spatially, therefore $\boldsymbol {S}_{\boldsymbol {0}}$ is a function of spatial coordinates $(x,y)$. Similarly, $n_s$ varies not only spatially, but also depends on time through seasonal variations in vegetation (Song et al. Reference Song, Schmalz, Xu and Fohrer2017).

There is an ongoing discussion whether using the above 2-D model (2.4)–(2.7) is an appropriate model for the overland flow, since the actual overland flow reaches the channel through a series of channels (natural or artificial ones), rather than being a uniform sheet of surface water $h_s(x,y)$ (see an overview by Leibowitz et al. (Reference Leibowitz, Wigington, Schofield, Alexander, Vanderhoof and Golden2018)). Nevertheless, we use this formulation as the current standard used in computational physical catchment models.

2.3. Channel flow

Finally, consider the channel flow illustrated in figure 1. We ignore the flow dynamics in the transverse directions of the channel and consider a channel located along $(x(s), y(s))$, for a parameterisation parameter, $s$, defined as the distance from the catchment outlet measured along the channel. Then, the channel water height is described by the Saint-Venant equation applied to the height $z = h_c(s, t)$. There are three main differences between the channel flow equations and the Saint-Venant equations used in the 2-D formulation of overland flow in § 2.2: (i) the direction of flow is given by the direction of the channel; (ii) precipitation adds a negligible contribution to the channel flow – instead, the channel flow is primarily affected by both the surface flow passing through the channel perimeter and the overland flow passing over the river banks; and (iii) the roughness is introduced over the entire perimeter of the channel (e.g. in the case of a rectangular channel, its bottom and the submerged part of its walls).

Mathematically, these assumptions lead to the following governing equation (Vieira Reference Vieira1983; Chaudhry Reference Chaudhry2007):

(2.9)\begin{equation} \frac{\partial{A}}{\partial t} = w\frac{\partial{h_c}}{\partial t} = q_{in} - \frac{\partial q}{\partial s}, \end{equation}

where $A = A(s, h_c)$ is the channel cross-section, $w = w(h_c, s) = {{\rm d} A}/{{\rm d} h_c}$ is the channel width (constant in the case of a rectangular channel) and $q=q(s,t)$ is the mean flow in the channel. The source term, $q_{in}$, is given by considering the total surface and subsurface inflow into the river; hence $q_{in}$ is linked to the boundary values of the solutions of the Richards (2.1) and 2-D Saint-Venant equations (2.4). As for the overland equations, the flux, $q = \mathrm {area} \times \mathrm {velocity} = Av$, rather than being solved using the full momentum equation, is again assumed to be given by the empirical Manning's law, which in the case of a channel takes the form

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

where $P = P(s, h_c)$ is the channel wetted perimeter, and $n_c$ is Manning's coefficient dependent on the banks and channel bed roughness. The quantity $S_f$ is the friction slope, as defined by (2.7) (or its simplified forms), with $S_0$ representing the elevation gradient along the river, denoted here as $S_y$. In the case of a rectangular channel, $A=wh_c$ and $P=w+2h_c$. Additionally, we denote the surface water height ($h_c$) at the outlet as $d$. As before, the kinematic approximation (2.8) is often assumed.

We use the above Manning's law since it is a standard approach in computational hydrological models. However, there may be some scenarios in which this assumption may not be a valid approach. In supercritical flows (such as in the case of flash floods), the flow may rapidly vary along the $y$ direction, and other approaches should be considered instead (see e.g. Mujumdar Reference Mujumdar2001).

This completes our formulation of the system of three coupled partial differential equations (PDEs) used to describe the key water flow components at the catchment scale, namely Richards equation (2.1), the 2-D Saint-Venant equations for overland flow (2.4) and the 1-D Saint-Venant equation for channel flow (2.4).

2.4. Boundary and initial conditions

The solution of the governing PDEs (2.1), (2.4), (2.9) for the respective subsurface $h_g$, surface $h_s$ and channel $h_c$ flows must be accompanied by appropriate boundary and initial conditions. For instance, boundary conditions are required to specify the mass exchange between flow components, and these conditions may introduce additional catchment-dependent parameters such as the channel depth and the particular geometry at the river outlet. Computational models also require the specification of initial conditions, which are generally unknown. Typically, the simulations can be run for a burn-in period to allow the results to be independent of the initial conditions. Example formulations of boundary and initial conditions will be studied in Part 2 of our work, where we focus on the mathematical and numerical analysis of model catchments.

2.5. Typical values of model parameters reported in the literature

All parameters appearing in the equations are hence summarised in . Firstly, let us provide a brief review of the typical values of these parameters, known from the literature.

1. (i) The range of saturated soil conductivity $K_s$ can vary from $10^0$ to $10^{-3}\ \mathrm {m\ s}^{-1}$ for very productive aquifers (for well-sorted sand and gravel, and highly fractured rocks) to below $10^{-9}\ \mathrm {m\ s}^{-1}$ for impervious rocks (Bear Reference Bear1972).

2. (ii) According to Chow (Reference Chow1959), Manning's roughness coefficient for channels, $n_c$, can vary from $0.01\ \mathrm {s\ m}^{-1/3}$ for artificial (e.g. cement) channels to over $0.1\ \mathrm {s\ m}^{-1/3}$ for channels with dense vegetation. Its value for flood plains, $n_s$, can vary from $0.03\ \mathrm {s\ m}^{-1/3}$ for pastures and cultivated areas without crops to $0.1\ \mathrm {s\ m}^{-1/3}$ in densely forested areas. In addition, the value of $n_c$ for a specific area can also vary with time due to seasonal variation in vegetation.

3. (iii) The evapotranspiration rate, $E$, in computational physics-based catchment models is usually estimated using models such as the evapotranspiration model by Kristensen & Jensen (Reference Kristensen and Jensen1975) and the two-layer unsaturated zone/evapotranspiration (UZ/ET) model by Yan & Smith (Reference Yan and Smith1994), both of which are for example used in the MIKE SHE (Système Hydrologique Européen) integrated model. However, the formulation of these models is beyond the scope of this report. Cole et al. (Reference Cole, Slade, Jones and Gregory1991) reports the mean monthly precipitation and evaporation values for different regions of the UK. Precipitation $R$ highly depends on the UK region, varying from $645.5$ mm year$^{-1}$ in central and East England to $1601.9$ mm year$^{-1}$ in Northwest and North Scotland, with the highest precipitation levels observed between October and January. According to Faulkner & Prudhomme (Reference Faulkner and Prudhomme1998), the highest daily precipitation measured throughout the year varies from $25$ mm (during a single day) in the east of England to over $80$ mm in some sites in mountainous regions of Wales and Scotland. The evapotranspiration rate $E$ is similar for all regions, but highly varies in time, from 6–12 mm month$^{-1}$ in January to 63–78 mm month$^{-1}$ in July.

4. (iv) Apart from specifying the precipitation, one also needs to specify catchment geometry – its size, terrain topography and aquifer depth. The publications on integrated catchment models mostly focus on one or a few real-world catchments. There are also works aiming to understand the characteristic properties of catchments and channel networks. Many catchment characteristics were described by Horton (Reference Horton1932), including drainage density, $D_D$, defined as the ratio of total stream length, $L$, and catchment area, $A$ (the authors found typical values for investigated catchments to be $A = 0.64\unicode{x2013} 1.367\ \mathrm {km}^{-1}$), average distance between streams (0.73–1.56 km), average overland flow distance (0.38–0.80 km), slope along the streams (0.014–0.038) and slope along the land (0.072–0.177). The cited values were estimated manually based on topographic maps. Together with the development of computing power and the collection of digital terrain models (DTMs), the methods for estimation of the above quantities were automated (Tarboton & Ames Reference Tarboton and Ames2001).

5. (v) Grieve, Mudd & Hurst (Reference Grieve, Mudd and Hurst2016) compared three different estimation methods for the hillslope width, $L_x$, by estimating its distribution for a number of catchments in the USA. The first estimate was obtained by dividing the catchment area, $A$, by twice the total stream length, $2L$, which is equivalent to the drainage density $(2D_D)^{-1}$. Note that the factor of two accounts for each stream having a hillslope on either side. The second method used a DTM model to find streamlines following the direction of the steepest descent. The lengths of the streamlines are interpreted as a hillslope width. The last method applied by the authors used a slope–area plot to separate areas dominated by channels from hillslopes. The hillslope width, $L_x$, defined using the streamline (flow routing) method is higher than estimates using the drainage density method, however, both have a similar order of magnitude ranging from $30$ to $130$ m.

6. (vi) A systematic study of the typical parameter ranges and their effect on model predictions was done by Doummar, Sauter & Geyer (Reference Doummar, Sauter and Geyer2012). The authors used the MIKE SHE software for a karst system (the Gallusquelle spring in the Southwest Germany). The model parameters were calibrated to minimise the root-mean-square error of the daily observed discharge, and the relative importance of each parameter was numerically investigated using sensitivity analysis. The most significant parameters include the hydraulic conductivity, the moisture content of the unsaturated rock matrix and van Genuchten parameters.

Table 1. List of parameters appearing in the formulation of the integrated catchment model, and where they are discussed in this work. As discussed in the text, these parameters vary spatially, and in some cases also temporally.

3.1. Catchment dimensions ($L_x$, $L_y$)

Capturing catchment dimensions and topography is challenging, as real-world river networks have a fractal-like structure (Rodriguez-Iturbe & Rinaldo Reference Rodriguez-Iturbe and Rinaldo2001). An ambiguous quantity concerns the estimation of the river length characterising the catchment length, $L_y$, since there are multiple ways in which it can be defined, and, furthermore, its value depends on the precision of the dataset used for estimation. Similarly, the estimate of the catchment/hillslope length, $L_x$, is challenging on account of the high spatial variation and ambiguous definition. The aquifer thickness, $L_z$, which depends on the properties of the soil, will be discussed in § 3.3.

In our work, we use three OS datasets, providing different data formats and levels of detail:

1. (i) the OS Open Rivers only consists of major rivers represented as spatial lines;

2. (ii) the OS VectorMap District includes all surface water bodies – wide rivers and standing bodies of water (e.g. lakes and ponds) are represented using spatial polygons, while narrow streams, artificial channels and ditches are represented using spatial lines;

3. (iii) the OS Water Network includes all rivers/streams forming the drainage network, but data can be downloaded only for small user-specified regions, not for the entire UK as in the case of the other two datasets.

Figure 2 demonstrates a comparison of the typical information provided by the three types of datasets.

Figure 2. Comparison of three spatial datasets of river locations (thick black lines): (b) and (c) offer a similar level of accuracy. The background colours represent groundwater depth from the BGS Groundwater Levels dataset in metres, indicating the areas of missing streams in (a), where groundwater reaches the surface in a characteristic finger-like pattern. The northings and eastings on the axes are expressed in kilometres.

3.1.1. Estimating the catchment length, $L_y$

In order to estimate the characteristic scale of the catchment length ($L_y$), we propose the following four measures, graphically represented in figure 3. Note that depending on their size, the OS Open Rivers dataset represent rivers in the form of polygons (we refer to them as major rivers) or lines (minor rivers).

1. (i) Length of main rivers ($L_y^{main}$) – the total length of rivers as defined in the OS Open Rivers dataset; we can use this measure to estimate the total flow into the drainage network.

2. (ii) Length of all rivers ($L_y^{all}$) – the total length of major rivers and minor streams as defined in the OS VectorMap District dataset; it serves the same function as the previous measure, but at a higher spatial resolution; we also use this measure later as one way of measuring the catchment/hillslope width.

3. (iii) Length of the longest river ($L_y^{long}$) – the length of the longest river measured from the spring to the catchment's outlet, extracted from the OS Open Rivers dataset (this measure is approximately the same regardless of whether the minor rivers are included or not); it may be used to investigate the characteristic length of the channel when studying the channel flow, given by (2.9).

4. (iv) Distance between river tributaries ($L_y^{trib}$) – the average distance between river tributaries estimated from the OS Open Rivers dataset, including the first-order streams (see figure 3d); it is the only intensive quantity in this list (i.e. it does not scale with the catchment size) and therefore can be used as the characteristic length of a river in which hillslope flow is not disturbed by the presence of tributaries.

Figure 3. Illustration of the different length scales characterising the drainage network, presented for the River Frome catchment, with the outlet located at Bishop's Frome. Each $L_y$ estimate is equal to the total length of all highlighted streams. Note this total length gradually decreases from (a) to (d), therefore $L_y^{all}\leq L_y^{main} \leq L_y^{long} \leq L_y^{trib}$.

Calculating each of the above measures poses some challenges. Firstly, the river banks and other natural boundaries have a fractal structure, which means that their length can be very sensitive to the chosen spatial resolution. Here, we have used the data in the original resolution to maintain high accuracy in the case of estimates (i)–(iii). In the case of estimate (iv), we are interested in the typical lengths of hillslopes rather than river streams. In this case, the length of the former should not be affected by small-scale local meanders. Therefore, in case (iv), before measuring the river length, we simplify the geometry using the Douglas–Peucker algorithm (Douglas & Peucker Reference Douglas and Peucker1973) with tolerances of $1000$ m. This algorithm effectively smooths river meanders shorter than the chosen tolerance and follows a similar methodology as in Stuetzle, Franklin & Cutler (Reference Stuetzle, Franklin and Cutler2009).

Secondly, when using the OS VectorMap District dataset for calculating the total river/stream length, data corresponding to standing bodies of water (lakes and ponds) should not be taken into account. This cannot be easily implemented since such standing bodies are represented in the same way as for wide rivers within the dataset. Nevertheless, since lakes have much shorter total boundary length than rivers for the majority of UK regions, including such data does not significantly impact the estimates. Therefore, we include all surface water bodies when estimating (ii), which we obtain by adding the sum of the lengths of all spatial lines and the sum of the perimeters of all spatial polygons divided by two (since each of the two banks is counted separately).

For each UK catchment specified in the NRFA, we calculated the total length of all major rivers ($L_y^{main}$), and all major and minor rivers ($L_y^{all}$), as well as the length of the longest stream ($L_y^{long}$) and the mean distance between major tributaries ($L_y^{trib}$). Note that from their definition for each catchment we have, $L_y^{all}\leq L_y^{main} \leq L_y^{long} \leq L_y^{trib}$ (see figure 3). Their median values taken over all considered catchments are: $\mathrm {med}(L_y^{main})=71.6$ km; $\mathrm {med}(L_y^{all})=130$ km; $\mathrm {med}(L_y^{long})=22.8$ km; $\mathrm {med}(L_y^{trib})=945$ m.

3.2. Catchment topography ($S_x$, $S_y$)

The gradient of the terrain is important in order to estimate the size of the surface flow, as given by Manning's law (2.6). In order to estimate the typical values of the elevation gradient perpendicular to the river ($S_x$) and along the river ($S_y$), we use the DTM from the OS Terrain 50 dataset.

3.3. Soil and rock properties ($L_z$, $K_s$, $\alpha _{MvG}$, $\theta _s$, $\theta _r$, $n$)

In this section, we focus on determining the hydraulic properties of soil and rock, which are important to estimate the saturated and relative hydraulic conductivities, $K_s$ and $K_r$, appearing in (2.1).

The geological structure and hydrological properties of soils differ significantly across the UK. In some areas, such as highly productive chalk aquifers in Southeast England, the soil has a very high conductivity, and almost all the rainwater reaches the groundwater table. On the other end of the spectrum, there are aquifers with essentially no groundwater, where the entire rainfall reaches the river either as subsurface flow through the soil or forms an overland flow. Based on the 625 K digital hydrogeological map of the UK developed by the BGS, 15 % of the UK area is classified as a highly productive aquifer, 26 % as moderately productive, 47 % as low productive and the remaining 12 % as rock with essentially no groundwater. The geographical distribution of aquifers of different productivity levels is presented in figure 4(a).

Figure 4. The hydrodynamic properties of aquifers in the UK as provided by the 625 K digital hydrogeological map by the BGS. (a) The map presents the productivity of the aquifer, while (b) presents the dominating mechanisms for the groundwater flow.

3.3.1. Estimating the soil depth, $L_z$

The typical depth of the conductive layer of soil can be estimated based on the UK3D dataset, which was constructed by the British Geological Survey based on measurements from 372 deep boreholes (Waters et al. Reference Waters, Terrington, Cooper, Raine and Thorpe2016). The authors of this dataset interpolated vertical profiles, which were then interpolated to form a national network, or ‘fence diagram model’, of bedrock geology cross-sections. Each segment of a cross-section is assigned to one of the following qualitative aquifer designations (see EA 2017).

1. (i) Principal aquifers: layers with high intergranular and/or fracture permeability that can support river base flow on a strategic scale.

2. (ii) Secondary aquifers A: permeable layers capable of forming an important source of base flow at a local rather than strategic scale.

3. (iii) Secondary aquifers B: predominantly lower permeability layers.

4. (iv) Secondary undifferentiated: assigned in cases where it has not been possible to attribute either category A or B to a rock type.

5. (v) Unproductive strata: layers with low permeability and negligible significance for river base flow.

We estimate the typical depth of each layer by calculating the total area of each type of rock layer and dividing it by the total length of the cross-sections, which are available in the dataset (approximately $20\,000$ km). More details about the data format and our processing methods are presented in Appendix A.4. Our analysis yields the following mean thicknesses, rounded to the nearest metre, for the different types of aquifers presented in the classification (i)–(v) above:

(3.1)\begin{equation} \left.\begin{array}{c@{}} \textrm{Principal} = 405\ \mathrm{m}, \quad \textrm{Secondary A} = 946\ \mathrm{m}, \quad \textrm{Secondary B} = 946\ \mathrm{m}, \\ \textrm{Secondary Undifferentiated} = 132\ \mathrm{m}, \quad \textrm{Unproductive} = 75\ \mathrm{m}. \end{array}\right\} \end{equation}

Thus, the mean principal aquifer thickness is $L_z = 405$ m. However, in some catchments in the UK, this layer may be much smaller – even practically reaching zero thickness, when the groundwater flow is insignificant (notice the classification of an aquifer with no groundwater sites in figure 4b). Therefore, such catchments need to be studied separately. In these low-productive aquifers, the subsurface flow takes place dominantly through the thin layer of soil confined from the bottom by impenetrable bedrock. Detailed quantitative data on soil thickness is not available; however, according to the maps shared by the UK Soil Observatory (Lawley Reference Lawley2009), over a half of the UK consists of deep soils with thickness exceeding 1 m, and nearly half of the soils are less than 1 m thick (see figure 5). Therefore, in the limiting case of catchments with no groundwater, we propose $L_z=1$ m as a reasonable choice for the typical soil thickness in regions with shallow impenetrable bedrock.

Figure 5. Thickness of soil according to the BGS Soil Parent Material Model, which provides a wide range of physical and chemical characteristics of the top layer of soil over the UK.

3.3.2. Estimating the MvG parameters

The necessary quantitative data on rock hydraulic properties (hydraulic conductivity and MvG parameters) would ideally be estimated directly through detailed experiments. However, they can also be calculated indirectly based on known soil composition. The ESDAC shares the 3-D Soil Hydraulic Database Europe, which uses European pedotransfer functions to estimate all MvG parameters (Tóth et al. Reference Tóth, Weynants, Pásztor and Hengl2017). The data is available at $1$ km and $250$ m resolutions for seven different depths (0, 5, 15, 30, 60, 100 and 200 cm). However, these estimates are not very precise since each property in the original database takes only one of a few discrete values, corresponding to specific soil types (see figure 6). In order to provide a single estimate of these parameters, we use data corresponding to $30$ cm and extract their mean value for individual catchments.

Figure 6. The MvG model parameters (from left to right: $K_S$; $\theta _R$; $\theta _S$; $\alpha _{MvG}$; and $n$) at a depth of 15 cm (a) and 1 m (b) according to the 3-D Soil Hydraulic Database Europe. As the legend indicates, each parameter in this dataset can only have one of a few discrete values.

3.3.3. Estimating hydraulic conductivity $K_s$

The values of the conductivity from the 3-D Soil Hydraulic Database, which range from around $10^{-8}$ to $10^{-6}$ m s$^{-1}$ (or approx. 0.001–0.1 m day$^{-1}$). However, as we shall discuss below, the direct measurements conducted in several highly productive chalk aquifers in England give higher estimates, likely because of the flow through macropores and fractures, which dominates over the intergranular flow in many areas in the UK (see figure 4b). Many studies, therefore, focus not only on analysing the hydraulic conductivity of the rock matrix, but also the bulk hydraulic conductivity, which includes the effect of both micropores forming the matrix and naturally occurring macropores/fractures.

In table 2, we cite several prior studies that have shown that the ratio of field-to-laboratory hydraulic conductivity values can be even greater than $1:1000$. For instance, measurements conducted by Robins & Buckley (Reference Robins and Buckley1988) on several boreholes in the Permian and Triassic aquifers in Southwest Scotland showed that hydraulic conductivity measured in the field varies from $0.1$ m day$^{-1}$ (Solway basin) to $20$ m day$^{-1}$ (Dumfries basin). In the study of chalk aquifers in the South Downs, Jones & Robins (Reference Jones and Robins1999) measured hydraulic conductivity values ranging from $0.15$ to $0.67$ m s$^{-1}$. Gardner et al. (Reference Gardner, Cooper, Wellings, Bell and Hodnett1990) observed that the hydraulic conductivity of English chalk aquifers increases from 1–6 mm day$^{-1}$ for low hydraulic potential, up to 100–1000 mm day$^{-1}$, as the potential is increased above $5$ kPa. This large increase in conductivity is caused by the fractures becoming saturated after increasing the potential.

Table 2. Comparison of laboratory (matrix) and field (fractures and matrix) hydraulic conductivity [m day$^{-1}$] for three different types of aquifers discussed by Robins & Ball (Reference Robins and Ball2006). Field and laboratory values can greatly differ.

Therefore, for highly conductive aquifers, we take typical values of conductivity ranging from $K_s=10^{-6}$ m s$^{-1}$ ($0.09$ m day$^{-1}$) to $K_s=10^{-4}$ m s$^{-1}$ ($9$ m day$^{-1}$), but these values can be significantly lower in low-productive aquifers. In the limiting scenario of catchments with no groundwater flow, the typical conductivity is given by the conductivity of the top layer of soil.

According to Kirkham (Reference Kirkham2014), the hydraulic conductivity of natural soils can vary from $0.05$ m day$^{-1}$ for a clay to $30$ m day$^{-1}$ for a silty clay loam. The variation is higher for disturbed soil materials, ranging from $0.02$ m day$^{-1}$ for silt and clay to $600$ m day$^{-1}$ for gravel. Therefore, for the soil, we may assume the same range. Even though the variation is higher than in the aforementioned studies for rock bulk conductivity, the typical values remain similar – between $10^{-6}$ and $10^{-4}$ m s$^{-1}$.

3.4. Water balance terms ($R$, $Q$, $E$)

The source term of the Saint-Venant equation for the overland flow (2.4) can be found by estimating the mean value of terms included in the water balance. It is given by

(3.2)\begin{equation} R=Q-E-\Delta S, \end{equation}

where $R$ is precipitation over a given catchment, $Q$ is river flow (runoff) at the outlet, $E$ is total evapotranspiration and $\Delta S$ is the change in the water volume stored within the catchment (e.g. in groundwater and reservoirs). Here, we not only estimate the values of $R$, $Q$ and $E$ for the UK, but also their annual peak values, which can provide a good parametrisation to model extreme precipitation events (and are used in many flood estimation models, e.g. by Kjeldsen, Jones & Bayliss (Reference Kjeldsen, Jones and Bayliss2008)). We define these quantities in terms of the mean flow per unit area of the catchment (measured in metres per second).

The precipitation and river flow data for all UK catchments are available in the NRFA. We use the standardised annual average rainfall (1961–90) to estimate $P$, and the mean of gauged daily flow to estimate $Q$ (divided by the catchment's area, also provided by NRFA). By averaging all terms over a long period of time (several years), we should expect $\Delta S$ to tend to $0$, which allows us to estimate the mean actual evapotranspiration as $E=Q-R$.

The relationship between $Q$ and $R$ is presented in figure 7(a). The mean runoff is smaller than the mean rainfall, as expected (especially for catchments with low mean rainfall). However, there are some outliers, especially among dry catchments. They correspond to situations where the mean gauged river flow exceeds the mean rainfall rate in a given catchment. Several possible reasons for this observation include the inflow of water from neighbouring catchments, a long-term trend of decreasing groundwater storage ($\Delta S$), or the mismatch between the time for which the precipitation and runoff data were available (in some catchments, the runoff data is available only for a limited period).

Figure 7. Dependencies between water balance terms obtained from the NRFA. Panel (a) shows the relationship between the mean rainfall, $R$, and the mean runoff, $Q$. The difference, $R-Q$, which is a result of evapotranspiration, is approximately constant, $E \approx 1.4\times 10^{-8}\ \textrm {m s}^{-1}$ for the UK catchments. Panel (b) shows the relationship between the median of the annual maximum rainfall (RMED) and the mean rainfall, $R$. The thick solid lines in both plots correspond to the line of best fit; in (b) the fitted line confirms that RMED scales approximately linearly with $R$.

Based on this graph, we can deduce that the evapotranspiration is approximately independent of the other parameters $R$ and $Q$. Rainfall $R$ values vary from $10^{-8}$ to $10^{-7}$ m s$^{-1}$, runoff $Q$ from $10^{-9}$ to $10^{-7}$, while mean evapotranspiration $E$ is usually around 1–1.5$\ \times \ 10^{-8}$ m s$^{-1}$.

Finally, we estimate annual peak values of precipitation, as a measure of how intensive rainfall should be considered in rainfall–runoff models. Following Faulkner & Prudhomme (Reference Faulkner and Prudhomme1998), the standard quantity we can use is the RMED. As shown in figure 7(b), it is closely correlated with mean precipitation $R$, and is approximately 8–10 times higher.

3.5. Manning's coefficients ($n_s$, $n_c$)

Manning's roughness coefficient, which appears in the overland flow equation of (2.6), is not measured directly at the catchment scale, but is either obtained by physical model calibration or in controlled small-scale experiments. Therefore, we take typical values from engineering tables.

The Manning's coefficient of the surface, $n_s$, depends on the type of the surface (Chow Reference Chow1959; Brunner Reference Brunner2016). The NRFA includes information about the area of each catchment covered by arable lands and grasslands ($n_s\approx 0.035\ \mathrm {s\ m}^{-1/3}$), mountains ($n_s\approx 0.025\ \mathrm {s\ m}^{-1/3}$), urban areas ($n_s\approx 0.015\ \mathrm {s\ m}^{-1/3}$) and woodlands ($n_s\approx 0.16\ \mathrm {s\ m}^{-1/3}$). Following these estimates, for each investigated catchment, we calculated a mean value of $n_s$ weighted by the area covered by each of the above terrain types.

Similarly, Manning's coefficient for the channel, $n_c$, depends on the type of the channel. Since there is no UK-wide database on channels types, we took the typical range of $n_c$ values from Chow (Reference Chow1959). They vary from $n_c=0.01\ \mathrm {s\ m}^{-1/3}$ (for smooth cement channels) to $n_c=0.1\ \mathrm {s\ m}^{-1/3}$ (for natural channels with very weedy reaches), with typical values around $n_c=0.035\ \mathrm {s\ m}^{-1/3}$.

3.6. Channel dimensions ($w$, $d$)

Interestingly, estimating the typical channel dimensions (width $w$ and depth $d$), which appear in Manning's law for the channel flow (2.10), turns out to be quite challenging. Firstly, there is no database that provides a complete set of data on channel dimensions for UK rivers. Secondly, channel dimensions vary greatly depending on the river discharge, which can itself vary greatly depending on the catchment (though mainly on the catchment area and average rainfall). Instead, we perform the analysis via an indirect route.

Our central reference is the work by Nixon (Reference Nixon1959), who measured the width and depth of channels of several major rivers at 29 gauging stations in England and Wales. Values ranged from a width of $w=35$ ft ($11$ m) and a depth of $d=3.1$ ft ($0.94$ m) for the River Blackwater at the Shallowfield station to a width of $w=257$ ft ($78$ m) for the Thames at the Kingston Day's Weir station and a depth of $d=16.37$ ft ($4.99$ m) for the Wye River at the Cadora station. The data, together with similar measurements done for American rivers, were used to fit power laws describing the dependence between the channel dimensions and the full-bank discharge $Q$. From Nixon (Reference Nixon1959), the derived power laws are

(3.3a)$$\begin{gather} d=C_1Q^{1/3}\quad\text{with}\ C_1=0.545\ \mathrm{s}^{1/3}, \end{gather}$$

(3.3b)$$\begin{gather}w=C_2Q^{1/2}\quad\text{with}\ C_2=1.65\ \mathrm{s}^{1/2}\ \mathrm{ft}^{{-}1/2}=3.00\ \mathrm{s}^{1/2}\ \mathrm{m}^{{-}1/2}, \end{gather}$$

where $Q$ is the full-bank discharge, while $C_1$ and $C_2$ are fitted constants. These scaling laws can also be justified by analysis of mechanical equilibrium of the sediment bed, see e.g. Henderson (Reference Henderson1963) and Devauchelle et al. (Reference Devauchelle, Petroff, Lobkovsky and Rothman2011).

Because data on channel dimensions is not available for most UK catchments, we use the above power laws to estimate the channel dimensions at the outlet of each studied UK catchment. We note that this is a reasonable approach: Nixon (Reference Nixon1959) derived such laws based on empirical relations for a representative sample of UK catchments. Equally, Wharton (Reference Wharton1992) used similar scaling laws to estimate the peak annual flows in 70 UK catchments based on the channel dimensions.

3.7. A summary of results on characteristic values

Table 3 summarises the typical values of catchment model parameters for all gauged catchments in the UK, with boundaries defined as in the NRFA. A sample of raw data and their spatial distributions is included in Appendix B. In cases where no spatial data was available, an estimate from the literature was obtained (these entries are marked with a ${\dagger}$). Where available estimates obtained with more than one method are compared. The median value of the parameters presented in table 3 will be used in Part 2 to formulate a simple benchmark scenario characterised by similar physical properties to the UK catchments.

Table 3. Summary of catchment model parameters values for UK catchments.

$^{\dagger}$Numbers represent a range of values obtained from the literature; these quantities are not estimated individually on the catchment level.

$^\textrm {a}$Hydraulic conductivity of both bedrock and soil estimated based on various studies. $^\textrm {b}$Hydraulic conductivity of the soil as given in Hydraulic Database of Europe by ESDAC.

Based on table 3, we can make the following notes.

1. (i) The total stream length in the catchment is approximately twice as high when including all surface water streams and bodies from the OS VectorMap ($L_y^{all}$) than when looking only at the main rivers, which are included in OS Open Rivers dataset ($L_y^{main}$), cf. figure 2. The longest stream ($L_y^{long}$) typically contributes to from 20 % to 50 % of the total length of the rivers in OS Open River, with a higher percentage corresponding to shorter rivers (with fewer and shorter tributaries). The estimated lengths greatly vary as a result of the wide range of gauged catchments areas used in the study. The exception is the distance between tributaries ($L_t^{trib}$) ranging from $725$ to $1212$ m. As an intensive quantity, this quantity does not scale with the catchment area.

2. (ii) Hillslope width estimates in the majority of catchments are consistent between the two proposed methods. A difference can be observed in catchments with the lowest density of streams, for which streamlines turn out to be much shorter than the width between streams. The reason is that the data includes numerous streamlines that start far from the river, but terminate at a local elevation minimum, and hence do not reconnect with the river. In our study, such streamlines cover over $18\,\%$ of the UK's area. We do not account for such streamlines when estimating the hillslope width.

3. (iii) Table 3 provides an indication that the slope along the hillslope $S_x$ is typically much larger than the slope along the channel $S_y$, which is a general feature of the topography of river valleys. By the 2-D Manning's equation (2.6), the flow follows the direction of the steepest descent. Therefore, we would then expect that the rainfall water is transferred towards the main river through streams mainly directed along the hillslopes (with $S_x$ slope). Only after it reaches the channel, it is transferred down the river with $S_y$ slope. The consequences of this observation will be explored in detail in Part 2.

4. (iv) The highest disagreement can be observed between the soil conductivity as reported in the ESDAC dataset and the literature. The latter measured the bulk soil conductivity of the soil (and bedrock), which is a better macroscopic estimate than the one provided in the ESDAC dataset.

5. (v) The last two columns of table 3 include example values from the literature. The first column includes parameter ranges suggested by Doummar et al. (Reference Doummar, Sauter and Geyer2012) to model a real-world Karst catchment in Southwest Germany. The catchment width is higher than typical values for the UK. Notable differences between the suggested parameter ranges and estimations for UK catchments are in soil conductivity, where higher values come from measurements of observed flow velocities using artificial tracers in highly conductive areas of the investigated system. These velocities can be very high in well-fractured Karst systems, as is also observed in some measurements in the UK (cf. the chalk aquifer of Yorkshire studied by Ward, Williams & Chadha (Reference Ward, Williams and Chadha1997)).

6. (vi) The last column includes parameter values suggested in two simple benchmark scenarios by Maxwell et al. (Reference Maxwell2014) for integrated catchment model intercomparison. The first scenario includes a single hillslope tilted in the $x$-direction, ending in a flat river of constant depth. Only subsurface flow in a thin 5 m-deep soil layer is modelled. The second scenario includes a 2-D tilted V-shaped catchment with a different elevation gradient along and perpendicular to the river. However, this latter one is only used to model surface flow – subsurface flow is not involved.

   One should note a few differences between parameter values used in these scenarios. The rainfall values are higher, since they correspond to typical intensive rainfall, not to an annual mean estimated for UK catchments. The next difference is a much lower aquifer thickness, since in benchmark scenarios, it represents only the soil layer ($L_z^*=5$ m), ignoring the subsurface water percolation to the usually much deeper permeable bedrock layer (on average $L_z=684$ m based on the UK3D dataset). Another important difference is the low soil conductivity in the scenarios: the saturated hydraulic conductivity, $K_s$, between $1.16\times 10^{-7}$ to $1.16\times 10^{-5}$ m s$^{-1}$ is lower than the reported typical values ranging from $10^{-4}$ to $10^{-6}$ m s$^{-1}$ (estimated based on large-scale studies in the UK, see § 3.3.3). Note that the combination of the two previous factors (thin aquifer with low conductivity) typically leads to model results, in simple benchmarks, where the majority of land is covered with surface water even for small rainfalls (cf. modelling and numerical simulations in Part 2).

4.1. Regional variability and correlation of parameters

Note that the estimates presented in table 3 refer to mean values across the UK. However, regional averages can significantly differ, as is represented by the interquartile range. As an example, consider the distribution of certain parameter mean values for UK catchments, as listed in the NRFA and presented in Appendix B. One can observe that highland areas, characterised by a higher elevation gradient (both along rivers and hillslopes), are also accompanied by higher rainfalls and lower soil conductivity (and usually low aquifer thickness). These observations are confirmed by the correlation matrix presented in the figure 8. The higher precipitation in mountainous regions was already observed by Duckstein, Fogel & Thames (Reference Duckstein, Fogel and Thames1973). Its impact on the aquifer properties (such as soil capacity and infiltration rate) was used by Bell & Moore (Reference Bell and Moore1998) to formulate the grid model and its successor, the grid-to-grid model (Bell et al. Reference Bell, Kay, Jones and Moore2007). The high correlation of catchment area and river length comes from the fractal structure of drainage networks (Rodriguez-Iturbe & Rinaldo Reference Rodriguez-Iturbe and Rinaldo2001), which develop to uniformly cover available space. Finally, there is a strong correlation between the $\alpha _{MvG}$, $n$ and $\theta _S$ parameters from the MvG model parameters, since all were estimated by Tóth et al. (Reference Tóth, Weynants, Pásztor and Hengl2017) based on the same soil composition dataset.

Figure 8. Correlogram of physical parameters summarised in table 3. It presents the Pearson correlation coefficient calculated based on values of the given parameters, which were estimated for UK catchments specified in the NRFA. Incomplete records were omitted.

4.2. Classification of catchments using cluster analysis

In order to better understand the similarities between different types of catchments, characterised by different combinations of physical parameters, we perform a cluster analysis. The goal is to construct a specified number of typical parameter combinations that reflect the spatial variation of UK catchments. Here, we use k-means clustering on the parameters summarised in table 3 to find a set of three clusters (see e.g. Kaufman & Rousseeuw Reference Kaufman and Rousseeuw2009). Parameters not directly associated or extracted with specific individual catchments in mind are removed from consideration ($w$, $d$, $n_c$, $L_z^*$). Also ignored are quantities that scale with the catchment areas, since these parameters are dependent on the location of the catchment outlet, rather than the physical properties characterising a given region ($L_y^{all}$, $L_y^{main}$ and $L_y^{long}$). Finally, the residual water content, $\theta _r$, was removed from consideration since it is almost always zero. This leaves us with the remaining 15 variables. Only catchments for which all these parameters are available are used for clustering.

The distribution of parameter values for each cluster is presented in figure 9. The figure highlights some key features of each cluster.

1. (i) From figure 9(f,g), we note that the first cluster is characterised by the highest elevation gradients, which mostly represent catchments located in highlands, while the third cluster is associated with the lowest elevation gradients, representing relatively flat lowlands. Intermediate values are represented by the largest second cluster.

2. (ii) Figure 9(i–l) shows that the third cluster, representing the lowlands, has significantly different soil properties from the other two clusters. It has the highest hydraulic conductivity, $K_s$, and the MvG $\alpha _{MvG}$ parameter, as well as the lowest saturated water content, $\theta _S$, and MvG coefficient, $n$.

3. (iii) Finally, in figure 9(b,c), we note that the first cluster, representing the highlands, has significantly higher mean precipitation, $R$, and runoff, $Q$, values, which is related to the fact that precipitation is higher in eastern parts of the UK, coinciding with many of the highland and mountainous areas.

Figure 9. Number of catchments belonging to each cluster and box plots showing the distribution of parameters among the catchments belonging to each graph.

4.3. Catchment characterisation using PCA

In order to better understand the separation between the class properties, one can use PCA, which allows us to find the linear combination of the original parameters with the highest variation in the dataset. Table 5 summarises the linear coefficients (eigenvectors) representing the first four principal components. Based on the parameters with the highest eigenvalues, our intuitive interpretation of the components is as follows. One can interpret the first principal component as strongly related to the topography of the terrain. The highland terrain with higher slopes and higher rainfall/runoff values is characterised by high values of the first principal components, while the lowlands correspond to low values. The second principal component is high for catchments with high aquifer thickness, high soil saturation, low $\alpha _{MvG}$ and/or high Manning's $n$ roughness coefficient. Typically, such a configuration of parameters represents catchments, in which groundwater flow has a relatively high importance. On the contrary, catchments with a low value of the second component will be characterised by a higher contribution of the overland flow component.

Figure 10 presents the clusters in a 2-D space spanned by the first two principal components. Here the separation between classes is apparent – as observed before, they vary in terrain topography. Additionally, the second cluster represents the most productive aquifers with the highest importance of groundwater. In the other two classes, the lower importance of groundwater is balanced by higher overland flows. In the case of cluster 1, the groundwater flow is limited by low hydraulic conductivity of the soil (cf. figure 9i), while in cluster 3, it is limited by the very flat terrain (cf. figure 9f,g), which can easily become saturated.

Figure 10. The illustration of the first two principal components. The graph on the left shows the values of principal components for catchments belonging to each of the three clusters. The map on the right shows the geographic distribution of these clusters.

Figure 11 also illustrates the spatial distribution of the first two principal components, as well as the geographical distribution of catchments belonging to different classes. Based on visual inspection of the distribution of regions associated with different values of the first principal component in figure 11(a), the measure seems to be highly spatially autocorrelated. This is due to the UK's distinct lowland, midland and highland regions, characterised by low, medium and high values of the first principal component, respectively. The spatial autocorrelation of the second principal component can also be observed in figure 11(b); however, this component varies much more locally – we expect it to be dependent on the local geological structure of the aquifer and the local soil composition.

Figure 11. Spatial distribution of values of the first (a) and second (b) principal components of the catchment parameters. Note that visually the value of the first component seems to coincide with the main lowland, midland and highland regions in the UK.