DRYP 1.0: a parsimonious hydrological model of DRYland Partitioning of the water balance

Dryland regions are characterised by water scarcity and are facing major challenges under climate change. One difficulty is anticipating how rainfall will be partitioned into evaporative losses, groundwater, soil moisture, and runoff (the water balance) in the future, which has important implications for water resources and dryland ecosystems. However, in order to effectively estimate the water balance, hydrological models in drylands need to capture the key processes at the appropriate spatio-temporal scales. These include spatially restricted and temporally brief rainfall, high evaporation rates, transmission losses, and focused groundwater recharge. Lack of available input and evaluation data and the high computational costs of explicit representation of ephemeral surface–groundwater interactions restrict the usefulness of most hydrological models in these environments. Therefore, here we have developed a parsimonious distributed hydrological model for DRYland Partitioning (DRYP). The DRYP model incorporates the key processes of water partitioning in dryland regions with limited data requirements, and we tested it in the data-rich Walnut Gulch Experimental Watershed against measurements of streamflow, soil moisture, and evapotranspiration. Overall, DRYP showed skill in quantifying the main components of the dryland water balance including monthly observations of streamflow (Nash–Sutcliffe efficiency, NSE, ∼ 0.7), evapotranspiration (NSE > 0.6), and soil moisture (NSE ∼ 0.7). The model showed that evapotranspiration consumes > 90 % of the total precipitation input to the catchment and that < 1 % leaves the catchment as streamflow. Greater than 90 % of the overland flow generated in the catchment is lost through ephemeral channels as transmission losses. However, only ∼ 35 % of the total transmission losses percolate to the groundwater aquifer as focused groundwater recharge, whereas the rest is lost to the atmosphere as riparian evapotranspiration. Overall, DRYP is a modular, versatile, and parsimonious Python-based model which can be used to anticipate and plan for climatic and anthropogenic changes to water fluxes and storage in dryland regions.

1 Introduction

Drylands are regions where potential evapotranspiration far exceeds precipitation and where water is scarce. Consequently, the water balance in such areas is highly sensitive to climatic forcing in terms of the delivery of precipitation and the evaporative demand from the atmosphere (Pilgrim et al., 1988; Goodrich et al., 1997; Zoccatelli et al., 2019; Kipkemoi et al., 2021). A key challenge is anticipating how rainfall partitioning into evaporative losses, groundwater, soil moisture, and runoff is likely to change under a future climate. Hydrological models provide important insights into the translation of climate information to water partitioning at or below the land surface. However, drylands exhibit several key hydrological processes that are distinct from humid regions and which are typically omitted from most current hydrological models (Huang et al., 2017). The lack of simple, computationally efficient hydrological models for drylands undermines efforts to anticipate and plan for climatic and anthropogenic changes to water storage and fluxes in catchments, with implications for water resources for ecosystems and society (Huang et al., 2017). Drylands cover around 40 % of the global land surface (Cherlet et al., 2018) and support a population of around 2 billion people (White and Nackoney, 2003), yet there are no widely available, parsimonious models for simulating the dryland water balance. Climatically, dryland regions are characterised by high rates of evapotranspiration and low annual precipitation delivered with high spatial and temporal variability (Wheater et al., 2007; Zoccatelli et al., 2019; Aryal et al., 2020). Precipitation events are characterised by high-intensity and low-duration rainfall over restricted spatial areas (Pilgrim et al., 1988). This results in a highly dynamic hydrological system prone to flash flooding, and also to water scarcity and food insecurity, societal risks that are exacerbated by climate change, population growth, and dryland expansion (Reynolds et al., 2007; Giordano, 2009; Siebert et al., 2010; Taylor et al., 2012; Huang et al., 2015, 2017; Wang et al., 2017; Cuthbert et al., 2019a).

In drylands, runoff occurs mainly as infiltration excess (Hortonian) overland flow due to high-intensity precipitation events, and it leads to the development of short-lived streamflow in ephemeral streams. These ephemeral streams play an important role in the water balance because high transmission losses of water through porous streambeds are the main source of aquifer recharge in such environments – a mechanism called focused recharge (Abdulrazzak, 1995; Coes and Pool, 2007; Goodrich et al., 2013; Shanafield and Cook, 2014; Cuthbert et al., 2016; Goodrich et al., 2018; Schreiner-McGraw et al., 2019). In contrast, diffuse recharge, which is the result of local infiltration of water below the evaporation zone within the soil, is typically limited in drylands due to low precipitation and high rates of evapotranspiration (Taylor et al., 2013; Schreiner-McGraw et al., 2019). These conditions result in dryland environments having no significant long-term storage of water within soils (Pilgrim et al., 1988; Huang et al., 2017).

The complexity of rainfall regimes, runoff generation processes, and subsurface flow paths in drylands create challenges for data collection, resulting in a paucity of data and consequent restrictions on the use of numerical models to enhance understanding of the water balance (Abbott et al., 1986; Woolhiser, 1990; Ewen et al., 2000; Woodward and Dawson, 2000; Michaelides and Wainwright, 2002; Ivanov et al., 2004; Šimunek et al., 2006; Wheater et al., 2007; Noorduijn et al., 2014; Schreiner-McGraw et al., 2019; Cuthbert et al., 2019b). Existing hydrological models, operating at catchment to regional scales, are challenged in drylands due to their inherent assumptions about key flow processes and due to the hard-coded parameterisations required to ensure convergence and numerical stability (e.g. physically based models; Kampf and Burges, 2007). Models also generally lack the ability to represent the development of ephemeral streams and their potential hydraulic interactions with groundwater systems (Quichimbo et al., 2020; Zimmer et al., 2020). Despite the recent improvement in models to include transmission losses (e.g. Hughes et al., 2006; Hughes, 2019; Lahmers et al., 2019; Mudd, 2006), the availability of appropriate numerical tools that allow for a better description of surface–groundwater interactions is still limited at catchment, regional, and global scales. Ephemeral flow in streams and water losses to the subsurface are currently underrepresented in medium- to large-scale models, despite representing half of the global stream network length (Datry et al., 2017; Messager et al., 2021). Additionally, the degree of complexity of existing models and their inherently high computational cost does not allow for comprehensive sensitivity and uncertainty analysis, which would support the evaluation and interpretation of model results.

In this context, it is important for models to capture the linkages between the spatially and temporally variable climate, nonuniform runoff generation, soil moisture, and focused groundwater recharge to support predictive capability of how the dryland water balance may shift with changes in climate. Models also need to include groundwater processes in drylands where the low regional hydraulic gradient governs the redistribution of groundwater resources, such as water availability for evapotranspiration in riparian areas (Maxwell and Condon, 2016; Mayes et al., 2020). Only a few large-scale hydrological models include gradient-based (diffuse) groundwater flow processes (Vergnes et al., 2012; de Graaf et al., 2015; Reinecke et al., 2019). These processes should be kept simple to make them transferable to different catchments regardless of their scale. Useful dryland models should also be able to employ the limited information available, while being numerically efficient enough to allow for evaluation of the model performance and uncertainty.

Here, we present the development of a parsimonious model which considers the main processes and spatio-temporal timescales that control the water partitioning, fluxes, and changes in water storage in dryland regions for estimation of runoff, soil moisture, actual evapotranspiration, and groundwater recharge. We do not intend for this model to accurately simulate event-based flood hydrographs, for example, for flood hazard analysis. Instead, we aimed to develop a model that captures the long-term behaviour of the water balance in dryland regions. Here, we apply and test our new model in the Walnut Gulch Experimental Watershed (WGEW), southeastern Arizona, USA, where availability of high-resolution data enabled us to evaluate the model performance.

2 DRYP: a parsimonious model for DRYland region water Partitioning

2.1 Model overview

The main hydrological processes that control fluxes and storage of water in dryland regions are shown in Fig. 1a. The movement of water through the different storage components within the catchment is characterised as follows: spatially distributed rainfall falling during individual events over the surface is partitioned into infiltration and runoff, depending on the temporal and spatial characteristics of the rainfall and the antecedent soil moisture conditions prior to the rainfall event (Goodrich et al., 1997; Zoccatelli et al., 2019). Water infiltrated into the soil can be extracted by plant evapotranspiration and/or soil evaporation, or it can percolate to the water table as diffuse recharge. Runoff is routed to the nearest stream based on topographic gradient. In each stream reach, water may be added through groundwater discharge as base flow or water may be lost through the porous boundaries by transmission losses as it moves downstream. The volumes of both base flow and transmission losses are dependent on the water table depth (Quichimbo et al., 2020). Transmission losses into the near-channel alluvial sediments increase the water available for plant evapotranspiration in the riparian zone and also generate focused recharge when the water holding capacity of the sediments in the riparian zone is exceeded (Schreiner-McGraw et al., 2019). Groundwater discharge into streams depends on the hydraulic gradient and occurs when the water table elevation is higher than stream stage elevation. Additionally, when the water table is close to the surface, capillary rise increases the root zone water availability for riparian plant evapotranspiration. Finally, anthropogenic activities, such as localised stream and groundwater abstraction as well as irrigation, may affect the storage and fluxes of the water balance.

Figure 1 Schematic representation of DRYP showing (a) the main hydrological processes controlling water partitioning in dryland regions; (b) distributed data sets needed to derive input parameters; (c) vertical and horizontal discretisation and representation of topographically driven surface runoff, vertical flow in the unsaturated zone, and hydraulic gradient-driven groundwater flow in the saturated component; and (d) model structure and potential processes within a single grid cell for the surface component (see Sect. 2.2), unsaturated zone (see Sect. 2.3), and saturated zone (see Sect. 2.4). Arrows represent flow directions, and red lines represent anthropogenic fluxes.

The only forcing variables in DRYP are spatially explicit fields of precipitation and potential evapotranspiration. The partitioning of the water balance then depends on the combination of this forcing and its interactions with spatially distributed parameters representing topography, land cover, soil hydraulic properties, hydrogeological characteristics of the aquifer, and anthropogenic activities (Fig. 1b). Hydrological processes in DRYP are structured into three main components: (i) a surface water component (SW) where precipitation is partitioned into infiltration and overland flow, which is then routed through the model domain based on the topographic gradient; (ii) an unsaturated zone (UZ) component that represents the soil and a riparian area parallel to streams; and (iii) a saturated zone (SZ) component that represents groundwater flow (Fig. 1c). All three components in DRYP are discretised as square grid cells, and all components are vertically integrated into a computational one-way sequential scheme (Fig. 1c). However, all components are hydraulically interconnected, allowing for gradient-driven and potentially bidirectional water exchange (Fig. 1c, d).

DRYP is written in Python and uses the Python-based “Landlab” package, which has the versatility to handle gridded data sets and model domains (Hobley et al., 2017; Barnhart et al., 2020). DRYP is structured in a modular way to allow user flexibility to control the desired level of process and parameter complexity, as well as the grid size and time-stepping choices appropriate for the desired application of the model. The grid size is the same for all layers, but the time step for different components may vary flexibly as described below. All grid cells potentially consist of all the process elements shown in Fig. 1d. However, the stream and riparian components can be excluded if stream channel characteristics are not provided, in which case all generated runoff in a cell will simply be routed to the next downstream cell with no additional losses or interactions. The scale of the stream and riparian zone is only limited by the grid size.

For all cells, at the beginning of every time step, the input rainfall (P) is partitioned into surface runoff (RO) and infiltration (I) depending on the available water content of the unsaturated zone (UZ). Water in the UZ can be extracted as actual evapotranspiration (AET), a combination of soil evaporation and plant transpiration, and/or percolate (R) to the saturated zone (SZ), depending on the water content and hydraulic properties of the unsaturated zone. If a cell is defined as a stream, transmission losses (TL) or groundwater discharge contributing to base flow (BF) and a riparian unsaturated zone (RUZ) are included in the local partitioning. The riparian zone is defined as an area parallel to the stream with a specified width. The riparian zone receives contributions from TL and a volume of infiltrated water proportional to the riparian area. Water within the riparian zone can either percolate, becoming focused recharge, or it can be extracted by plants as riparian evapotranspiration. Focused and diffuse recharge are combined as the main inputs to the SZ, which may also interact with the UZ depending on the water table elevation as it rises and falls through the simulation. The movement of water in the SZ is driven by the lateral hydraulic gradient. Additionally, anthropogenic interactions in the model are implemented as localised fluxes from the saturated zone (ASZ) and streams (AOF), whereas water abstraction for irrigation (AUZ) is delivered to the surface where it then contributes to infiltration into the unsaturated zone.

2.2 Model input files and parameter settings

DRYP requires spatial characterisation of key input parameters and data including a digital elevation model (DEM), channel properties in cells where streams are explicitly defined (length, width, and saturated hydraulic conductivity), land cover (plant rooting depth), various soil hydraulic properties, and aquifer properties (specific yield, aquifer thickness, and saturated hydraulic conductivity) (Fig. 1). Where local-scale information is not available for parameterisation, publicly available data at high spatial resolution at the regional and global scales can be considered for model parameterisation (e.g. Pelletier et al., 2016; Leenaars et al., 2018; Dai et al., 2019, etc.). A summary of model parameters for the different model components and structures is presented in Table 1. If parameters are not provided, “global” default values are used as defined in Table 1. The following sections describe the implementation of each process included in DRYP in detail. Precipitation and potential evapotranspiration are the only forcing variables and can be supplied as either spatially variable gridded data sets in netCDF (network common data form) format or as spatially uniform values for each time step. Gridded data sets must be interpolated or aggregated to match the model grid resolution.

Table 1 Model parameters for different processes considered in the model; some required parameters depend on the infiltration approach (“Inf. method”). For soil hydraulic properties, default values correspond to a sandy loam soil texture (Clapp and Hornberger, 1978; Rawls et al., 1982).

^* Default values correspond to a flow velocity of ∼ 1 m s⁻¹ over a 300 m straight path.

2.3 Surface component

Two main processes are considered in the surface component: (i) the partitioning of precipitation into infiltration and runoff, and (ii) runoff routing and the partitioning of runoff into streamflow and transmission losses in stream cells. These are described below.

2.3.1 Infiltration and runoff

The partitioning of precipitation into infiltration and runoff at the land surface is a key process in drylands and a potentially major source of uncertainty in the overall water partitioning for these regions. Hence, four different infiltration approaches have been included in DRYP, which can be toggled on or off within the main control file (prior to simulation) to allow the user to experiment with different infiltration model structures. These approaches include two point-scale methods: the Philip infiltration approach and the modified Green–Ampt method; and two upscaled methods for summarising infiltration over larger areas: the upscaled Green–Ampt and the multi-scale Schaake approach.

Method 1: infiltration based on Philip's equation

In this option, infiltration, f [L T⁻¹] during a rainfall event is based on the explicit solution of the infiltrability depth of Philip's equation (Philip, 1957):

$$\begin{array}{}\text{(1)}& f\left(t_c\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{S}_p{t}_c^{-\frac{\mathrm{1}}{\mathrm{2}}}+K_{\mathrm{sat}},\end{array}$$

where K_sat is saturated hydraulic conductivity [L T⁻¹], S_p is sorptivity [L² T${}^{\mathrm{1}/\mathrm{2}}$], and t_c is time since the beginning of the precipitation event [T]. The sorptivity term is estimated using the following equation (Rawls et al., 1982):

$$\begin{array}{}\text{(2)}& S_p={\left[\mathrm{2}{K}_s\left({\mathit{\theta }}_{\mathrm{sat}}-\mathit{\theta }\right)\left|{\mathit{\psi }}_f\right|\right]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

where θ is volumetric water content [L³ L⁻³], θ_sat is volumetric water content under saturated conditions [L³ L⁻³], and ψ_f is suction head [L]. ψ_f is estimated as follows (Clapp and Hornberger, 1978):

$$\begin{array}{}\text{(3)}& \left|{\mathit{\psi }}_f\right|={\mathit{\psi }}_a{\displaystyle \frac{\mathrm{2}\mathit{\lambda }+\mathrm{2.5}}{\mathit{\lambda }+\mathrm{2.5}}},\end{array}$$

where ψ_a is maximum suction head [L], and λ is a parameter that represents the pore size distribution of the soil [–] (Clapp and Hornberger, 1978).

The total infiltration depth in any given cell, I [L], during a precipitation event is estimated by solving the integral of Eq. (1) over the event duration. The integral of Eq. (1) is solved using the time compression approach (TCA) (Sherman, 1943; Holtan, 1945; Mein and Larson, 1973; Sivapalan and Milly, 1989), assuming that infiltration after ponding depends on the cumulative infiltrated volume. Therefore, to match the initial infiltration rate at the beginning of each time step with the infiltration at the end of the previous time step, the start time of infiltration is shifted to match the total cumulative infiltration. A more detailed description and the analytical solution of the approach can be found in Assouline (2013) and Chow et al. (1988).

Method 2: infiltration based on a modified Green–Ampt method

We have implemented a modified version of the Green–Ampt approach defined by the following equation (Scoging and Thornes, 1979; Michaelides and Wilson, 2007):

$$\begin{array}{}\text{(4)}& f\left(t_c\right)=K_{\mathrm{sat}}+{\displaystyle \frac{B}{t_c}},\end{array}$$

where B represents the initial suction head [L], and t_c is the same as Eq. (1); here, we use sorptivity (Eq. 2) as a proxy for the initial head owing to the non-linear dependency of sorptivity on the water content of the soil.

The integral of Eq. (4) was also solved using the time compression approach (Sherman, 1943; Holtan, 1945; Mein and Larson, 1973; Sivapalan and Milly, 1989). However, as there is no explicit solution for Eq. (4), we used an implicit solution.

Method 3: infiltration based on an upscaled Green–Ampt method

This method is based on the semi-analytical solution of the Green–Ampt equation for spatially heterogeneous hydraulic conductivity developed by Craig et al. (2010):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{I}\left(t_c\right)& ={\displaystyle \frac{p}{\mathrm{2}}}\mathrm{erfc}\left({\displaystyle \frac{\ln \left(pX\right)-{\mathit{\mu }}_Y}{{\mathit{\sigma }}_{Y\sqrt{\mathrm{2}}}}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}X}}\log \left|K_{\mathrm{sat}}\right|\mathrm{erfc}\left({\displaystyle \frac{{\mathit{\sigma }}_Y}{\sqrt{\mathrm{2}}}}-{\displaystyle \frac{\ln \left(pX\right)-{\mathit{\mu }}_Y}{{\mathit{\sigma }}_{Y\sqrt{\mathrm{2}}}}}\right)\\ & +p{\int }_{\mathrm{0}}^{X\left(t_c\right)}\mathit{\epsilon }\left(X\left(t\right),K_{\mathrm{sat}}\right)f_k\left(K_{\mathrm{sat}}\right)\mathrm{d}{K}_{\mathrm{sat}},\end{array}\end{array}$$

where $\widehat{I}$ is the mean infiltration rate [L T⁻¹], p is the precipitation rate [L T⁻¹], t_c is the same as in Eq. (1), f_k is the probability density function of K_sat, μ_Y and σ_Y are mean and standard deviation of the log saturated hydraulic conductivity, ${\mathit{\mu }}_Y=\ln \left|K_{\mathrm{sat}}\right|-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\sigma }}_Y$, and X is a dimensionless time. X is estimated as follows:

$$\begin{array}{}\text{(6)}& X={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{\mathit{\alpha }}{P{t}_c}}},\end{array}$$

where $\mathit{\alpha }=\left|{\mathit{\psi }}_f\right|({\mathit{\theta }}_{\mathrm{sat}}-\mathit{\theta })$, with ψ_f representing the suction head.

The ε(X,K_s) in Eq. (5) is an error function that can be estimated by the following approximation (Craig et al., 2010):

$$\begin{array}{}\text{(7)}& \mathit{\epsilon }\approx \mathrm{0.3632}\cdot {\left(\mathrm{1}-X\right)}^{\mathrm{0.484}}\cdot {\left(\mathrm{1}-{\displaystyle \frac{K_{\mathrm{sat}}}{pX}}\right)}^{\mathrm{1.74}}{\left({\displaystyle \frac{K_{\mathrm{sat}}}{pX}}\right)}^{\mathrm{0.38}}.\end{array}$$

The f_k(K_s) is assumed as a lognormal distribution following Craig et al. (2010):

$$\begin{array}{}\text{(8)}& f_K\left(K_{\mathrm{sat}}\right)={\displaystyle \frac{\mathrm{1}}{K_{\mathrm{sat}}{\mathit{\sigma }}_Y\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{{\left(\ln \left(K_{\mathrm{sat}}\right)-{\mathit{\mu }}_Y\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_Y^{\mathrm{2}}}}\right).\end{array}$$

As suggested by Craig et al. (2010), we solve the integral of the Eq. (5) efficiently using a two-point Gauss–Lagrange numerical integration method.

Method 4: infiltration based on the multi-scale Schaake method

The Schaake et al. (1996) approach is based on the assumption that rainfall and infiltration rates follow an exponential distribution to approximate the spatial heterogeneity of soil properties. Therefore, the spatially averaged infiltration I [L] is estimated as

$$\begin{array}{}\text{(9)}& I={\displaystyle \frac{P{I}_c}{P+I_{\mathrm{c}}}},\end{array}$$

where P is total rainfall [L], and I_c is cumulative infiltration capacity [L].

Infiltration capacity is estimated as follows (Schaake et al., 1996):

$$\begin{array}{}\text{(10)}& I_{\mathrm{c}}=\left({\mathit{\theta }}_{\mathrm{sat}}-\mathit{\theta }\right)\left(\mathrm{1}-\exp \left(-k_{\mathrm{d}t}\right)\right),\end{array}$$

where k_dt is a constant that depends on soil hydraulic properties.

Following Chen and Dudhia (2001), we define k_dt as

$$\begin{array}{}\text{(11)}& k_{\mathrm{d}t}=k_{\mathrm{d}{t}_{\mathrm{ref}}}{\displaystyle \frac{K_{\mathrm{sat}}}{K_{\mathrm{ref}}}},\end{array}$$

where K_ref [L T⁻¹] is a reference hydraulic saturated conductivity equal to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ (Wood et al., 1998; Chen and Dudhia, 2001), and the parameter $K_{\mathrm{d}{t}_{\mathrm{ref}}}$ is specified as a scale calibration parameter.

2.3.2 Runoff routing and transmission losses

Rainfall that does not infiltrate (i.e. precipitation, P, minus infiltration, I) into the unsaturated component is routed over the model domain based on topography. The flow routing scheme varies depending on whether a cell is defined as a stream. A simple flow accumulation approach is used in cells without a defined stream, whereas for defined stream cells, an additional flux term is added to the flow accumulation approach to account for groundwater interactions via the riparian zone. This flux will either be a transmission loss or a base flow contribution from the saturated component.

Flow routing in cells without streams

Runoff produced in any given cell is instantaneously routed to the next downstream cells using the flow accumulation approach implemented in Landlab (Braun and Willett, 2013; Hobley et al., 2017).

The next runoff downstream cell is estimated using a D8 flow direction approach (eight potential directions based on adjacent cells). The flow accumulation method adds the amount of runoff from the upstream cells:

$$\begin{array}{}\text{(12)}& Q_i={\sum }_{i=\mathrm{1}}^N{Q}_{{\mathrm{in}}_i},\end{array}$$

where Q_in [L³] is the volume of water that discharges from upstream cells into the current cell i, N is the number of upstream cells discharging into the current cell, and Q_i [L³] is the volume of water in the cell.

Flow routing in stream cells

In defined stream cells, the amount of water entering the cell, q_in [L³ T⁻¹], is instantly reduced by any transmission losses, i_ch [L³ T⁻¹], and any remaining water, q_out [L³ T⁻¹], is moved to the next downstream cell:

$$\begin{array}{}\text{(13)}& q_{\mathrm{out}}=q_{\mathrm{in}}-i_{\mathrm{ch}}.\end{array}$$

Water from the upstream cell, q_in, is assumed to be released to the next cell following a linear reservoir approach:

$$\begin{array}{}\text{(14)}& q_{\mathrm{in}}=q_{\mathrm{0}}{e}^{-k_T{t}^{*}},\end{array}$$

where k_T [T⁻¹] is a recession term that is equal to the inverse of the residence time of the streamflow at each cell, t^* represents time [T], and q₀ is the initial flow rate of water entering the channel. q₀ is estimated as

$$\begin{array}{}\text{(15)}& q_{\mathrm{0}}=\left(Q_{\mathrm{in}}+S_{\mathrm{SW}}-Q_{\mathrm{ASW}}\right)k_T,\end{array}$$

where Q_ASW [L³] is the volume of water abstracted from the stream, and S_SW [L³] is water stored in the channel.

It is assumed that the sediments in the streambed are homogenous. In order to use an explicit approach while also maintaining the simplicity of the model, the channel cross section is assumed to be rectangular. Consequently, the rate of infiltration depends on the wetted perimeter of the channel, and the infiltration rate, i_ch, at the stream cell is estimated assuming a unit gradient Darcian flow across the wetter perimeter:

$$\begin{array}{}\text{(16)}& i_{\mathrm{ch}}=K_{\mathrm{ch}}\left(\mathrm{2}y+W\right)L_{\mathrm{ch}},\end{array}$$

where K_ch [L T⁻¹] is saturated hydraulic conductivity of the streambed, L_ch [L] is channel length for a given cell, W is channel width [L], and y is streamflow stage [L]. If the rate of water entering the stream cell is less than the potential channel infiltration rate, flow to the next downstream cell is set to zero (all water is lost via infiltration) and i_ch=q_in.

Stream stage, y, is estimated by assuming that flow velocity does not change along the channel in any given cell (no flow acceleration). Therefore, the streamflow stage and the volume at any time along the channel are kept constant in any given stream cell. A constant velocity approach assumes that there are no backward effects on the streamflow routing approach. Thus, the stream stage is estimated as the height of the rectangular prism with area A=WL_ch and volume at time t as follows:

$$\begin{array}{}\text{(17)}& y={\displaystyle \frac{q_{\mathrm{in}}}{A}}.\end{array}$$

After substituting Eq. (17) into Eq. (16) and then into Eq. (13), the time integral of Eq. (13) represents the total amount of water, Q_out [L³], that moves to the next downstream channel cell (becoming Q_in):

$$\begin{array}{}\text{(18)}& Q_{\mathrm{out}}={\int }_{\mathrm{0}}^{\min [t_{q=\mathrm{0}},\mathrm{\Delta }t]}\left[q_{\mathrm{0}}{e}^{-k_{T}t}-L_{\mathrm{ch}}{K}_{\mathrm{ch}}\left(\mathrm{2}{\displaystyle \frac{q_{\mathrm{0}}{e}^{-k_{T}t}}{W{L}_{\mathrm{ch}}}}+W\right)\right]\mathrm{d}t.\end{array}$$

Note that the time step choice is important to bear in mind with respect to the size of the catchment modelled, as it represents the minimum travel time for flow to reach the catchment outlet.

The amount of water stored in the channel is estimated by applying a mass balance of all inputs and outputs of the channel:

$$\begin{array}{}\text{(19)}& S_{\mathrm{RO}}^t=Q_{\mathrm{in}}+S_{\mathrm{SW}}^{t-\mathrm{1}}-Q_{\mathrm{ASW}}-Q_{\mathrm{TL}}-Q_{\mathrm{out}},\end{array}$$

where t represents the current time step, and Q_TL [L³] is transmission losses estimated as the integral of the second term of Eq. (18). The total of Q_TL is restricted to the storage available in the aquifer:

$$\begin{array}{}\text{(20)}& Q_{\mathrm{TL}}=\min [Q_{\mathrm{TL}},\max [\left(z-h\right)A{S}_{\mathrm{y}},\mathrm{0}],\end{array}$$

where z is the surface elevation [L], h is water table elevation [L], A is the area of cell [L²], and S_y is aquifer specific yield [–].

2.4 Unsaturated component

Water infiltrated into the soil or through the stream channel becomes a flux input to the UZ (Fig. 1d). The unsaturated component comprises the soil and the riparian zone, both of which are simulated using a linear `bucket' soil moisture balance model (Fig. 2a), following an approach similar to the water balance model from the Food and Agriculture Organization (FAO) of the United Nations (Allen et al., 1998):

$$\begin{array}{}\text{(21)}& \mathrm{\Delta }{S}_{\mathrm{UZ}}=I+Q_{\mathrm{TL}}-\text{AET}-R,\end{array}$$

where ΔS represents storage change [L], AET represents actual evapotranspiration rate [L T⁻¹], and R represents the potential recharge rate [L T⁻¹]. The term Q_TL is only defined for stream cells. Diffuse potential recharge results from the local vertical percolation of the unsaturated zone, whereas focused potential recharge is produced in the riparian unsaturated zone (see Fig. 1).

Figure 2 Schematic illustration of the unsaturated component. The right panel represents the variation in the ratio of potential to actual evapotranspiration in relation to the water content of the soil. The reader is referred to Sect. 2.3 and 2.4 for a detailed explanation of the terms shown here.

The amount of water available for plant evapotranspiration in the UZ, L [L], is estimated as the product of the rooting depth, D_root [L], and the volumetric water content, θ [L³ L⁻³]. The maximum amount of water that the soil can store is limited by the field capacity of the soil (L_fc), whereas the minimum amount is constrained by the wilting point (L_wp). Thus, the total available water, L_TAW, for plant transpiration is estimated by the difference between L_fc and L_wp (see Fig. 2).

The potential amount of water that plants can remove from the UZ as transpiration, PET [L T⁻¹], is the result of the product between a crop coefficient, k [–], and the reference potential evapotranspiration, ET₀ [L T⁻¹] (Allen et al., 1998). When there is enough water to supply plant energy demands, water can be extracted from the UZ at a rate equal to the PET. However, when there is not enough water in the UZ to supply the PET, plants are considered to be under stressed conditions and the actual evapotranspiration (AET) is constrained as

$$\begin{array}{}\text{(22)}& \mathrm{AET}=I+\mathit{\beta }\left(\mathrm{PET}-I\right),\end{array}$$

where β is a dimensionless parameter that depends on the water content. β is estimated by

$$\begin{array}{}\text{(23)}& \mathit{\beta }={\displaystyle \frac{L-L_{\mathrm{TAW}}}{L_{\mathrm{TAW}}\left(\mathrm{1}-c\right)}},\end{array}$$

where c is the fraction of L_TAW [–] at which plants can extract water from the UZ without suffering water stress; c is set to 0.5, as recommended by the FAO guidelines (Allen et al., 1998), although this can be varied in DRYP.

If, after accounting for infiltration and AET, there is a surplus of water in the soil that exceeds the field capacity, diffuse recharge (R) to the groundwater system occurs. If the model is run at daily time steps, we assume that all water content above field capacity will percolate and produce R. However, for sub-daily time steps, it is more realistic to assume that the soil slowly releases water as R when it is above the field capacity, depending on the soil water retention curve. Hence, in this case, we assume that percolation to the water table depends on the water content and occurs only under the influence of gravity as follows:

$$\begin{array}{}\text{(24)}& D_{\mathrm{UZ}}{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}=-K\left(\mathit{\theta }\right).\end{array}$$

where D_UZ is the depth of the unsaturated zone [L] (see also Fig. 3a). K(θ) is estimated using the Brooks and Corey (1964) relations and the Clapp and Hornberger (1978) parameters (see Eq. 3):

$$\begin{array}{}\text{(25)}& K\left(\mathit{\theta }\right)=K_{\mathrm{sat}}{\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{sat}}}}\right)}^{\left(\mathrm{2}\mathit{\lambda }+\mathrm{2.5}\right)}.\end{array}$$

We then substitute Eq. (25) into Eq. (24) and assume that the soil drains immediately into the groundwater component after evapotranspiration loss. Hence, an analytical solution based only on drainage without considering other inputs or outputs is specified by

$$\begin{array}{}\text{(26)}& \mathit{\theta }=\exp \left(-\left(\mathrm{2}\mathit{\lambda }-\mathrm{1.5}\right)\log \left|{\mathit{\theta }}^{-\mathrm{2}\mathit{\lambda }-\mathrm{1.5}}-{\displaystyle \frac{\mathrm{\Delta }t\left(\mathrm{2}\mathit{\lambda }+\mathrm{1.5}\right)K_{\mathrm{sat}}}{D_{\mathrm{UZ}}{\mathit{\theta }}_{\mathrm{sat}}^{\mathrm{2}\mathit{\lambda }+\mathrm{2.5}}}}\right|\right).\end{array}$$

The UZ model component in DRYP can also change its behaviour when the head in the SZ component beneath restricts the downward movement of water. This case is described below in Sect. 2.5.1 (unsaturated–saturated zone interactions).

By default, the riparian zone uses the same hydraulic properties of the soil unsaturated zone except for the saturated hydraulic conductivity, which is assumed to be the same as the channel streambed K_ch; however, these parameters are also user-defined. The size of the riparian zone has a user-defined width (default is 20 m), and the length is the same as the stream.

2.5 Saturated component

Lateral saturated flow underneath the unsaturated zone assumes the Dupuit–Forchheimer conditions for the Boussinesq equation and Darcian conditions for flow in/out of each model cell:

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial h}{\partial t}}+q_{\mathrm{s}}+q_{\mathrm{riv}}={\displaystyle \frac{\mathrm{1}}{S_{\mathrm{y}}}}\mathrm{\nabla }\cdot \left(-K_{\mathrm{sat}}h\mathrm{\nabla }h\right)+R-Q_{\mathrm{ASZ}},\end{array}$$

where K_aq is the saturated hydraulic conductivity of the aquifer [L T⁻¹], S_y is the specific yield [–], q_s is saturation excess [L T⁻¹] (see Sect. 2.5.1), q_riv is discharge into stream [L T⁻¹] (see Sect. 2.5.2), Q_ASZ [L T⁻¹] is groundwater abstraction, ∇ represents the gradient operator, and ∇⋅ represents the divergence operator. Where the saturated thickness of the aquifer is relatively constant over the simulation period, transmissivity, T [L² T⁻¹] (the product of the aquifer thickness and the saturated hydraulic conductivity of the aquifer), may be held constant, thereby linearising Eq. (25). Additionally, an exponential function based on Fan et al. (2013) has been added to represent the reduction of transmissivity in relation to depth:

$$\begin{array}{}\text{(28)}& T=K_{\mathrm{sat}}{f}_{\mathrm{D}}\exp \left(-{\displaystyle \frac{z-h}{f_{\mathrm{D}}}}\right),\end{array}$$

where f_D is effective aquifer depth [L]. These different transmissivity parameterisation options can be toggled on or off in the main model control file.

Equation (27) is solved using a forward time central space (FTCS) finite difference approach. FTCS is an explicit finite difference approximation whose solution is sensitive to grid size and time step. Thus, in order to obtain a stable convergence of Eq. (27), a time-variable approach was adopted. The maximum allowable time step for the saturated component is estimated based on the Courant number criteria (we use 0.25 as a default value but this may be changed by the user):

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{T\mathrm{\Delta }t}{S_{\mathrm{y}}\mathrm{\Delta }{x}^{\mathrm{2}}}}\le \mathrm{0.25}.\end{array}$$

If the maximum time step of the SZ component is greater than the minimum time step of any other component of the model, the time step of the SZ component is reduced to the time step of the minimum time step of the model (see Sect. 2.6 for more details on the model time step options).

2.5.1 Unsaturated–saturated zone interactions

Unsaturated–saturated zone interactions are implemented using a variable-depth unsaturated zone as follows (Fig. 3a). Unsaturated zone thickness (D_uz) is equal to the rooting depth when the water table elevation (h) is below the rooting depth, but when the water table is above the rooting depth, the thickness of the unsaturated zone is reduced to the depth of the water table:

$$\begin{array}{}\text{(30)}& D_{\mathrm{uz}}=min\left[D_{\mathrm{root}},z-h\right].\end{array}$$

When the water table is below the rooting elevation, z_root, there is no two-way interaction between the soil and the groundwater compartment (only one-way, as recharge); thus, no updates to the water table elevation are required (see Fig. 3a1). However, when the water table crosses the z_root threshold, either via recharge or lateral groundwater flow, the water table is updated depending on the change in groundwater storage:

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{\Delta }{S}_{\mathrm{SZ}}}{\mathrm{\Delta }t}}=\mathrm{\nabla }\cdot \left(-K_{\mathrm{sat}}h\mathrm{\nabla }h\right)+R-Q_{\mathrm{ASZ}},\end{array}$$

where ΔS_SZ is the change in groundwater storage per unit area [L³ L⁻²]. Specifically, if an SZ cell is being recharged and the water table rises past the rooting depth in a given time step, the water table is updated according to

$$\begin{array}{}\text{(32)}& h_t={\displaystyle \frac{\mathrm{1}}{{\mathit{\theta }}_{\mathrm{sat}}-{\mathit{\theta }}_t}}\left[\mathrm{\Delta }{S}_{\mathrm{SZ}}-\left(z_{\mathrm{uz}}-h_{t-\mathrm{1}}\right)S_{\mathrm{y}}\right]+z_{\mathrm{uz}},\end{array}$$

whereas when the water table is draining and passes the rooting depth in a given time step,

$$\begin{array}{}\text{(33)}& h_t=-{\displaystyle \frac{\mathrm{1}}{S_{\mathrm{y}}}}\left[\mathrm{\Delta }{S}_{\mathrm{SZ}}-\left(h_{t-\mathrm{1}}-z_{\mathrm{root}}\right)\left({\mathit{\theta }}_{\mathrm{sat}}-{\mathit{\theta }}_{\mathrm{fc}}\right)\right]+z_{\mathrm{root}}.\end{array}$$

When the water table is above the rooting depth elevation, the water table elevation will be updated according to

$$\begin{array}{}\text{(34)}& h_t={\displaystyle \frac{\mathrm{\Delta }{S}_{\mathrm{SZ}}}{{\mathit{\theta }}_{\mathrm{sat}}-{\mathit{\theta }}_{\mathrm{fc}}}}+h_{t-\mathrm{1}},\end{array}$$

whereas if it is below the rooting depth elevation, the water table elevation is simply

$$\begin{array}{}\text{(35)}& h_t={\displaystyle \frac{\mathrm{\Delta }{S}_{\mathrm{SZ}}}{S_{\mathrm{y}}}}+h_{t-\mathrm{1}}.\end{array}$$

When the water table is above z_root, there is more water potentially available for evapotranspiration, as it can be taken from the groundwater reservoir via capillary rise or direct root water uptake. Thus, the potential maximum amount of water taken up from the groundwater reservoir, PAET_SZ [L T⁻¹], is computed as the remaining PET after AET from the unsaturated component as follows:

$$\begin{array}{}\text{(36)}& {\mathrm{PAET}}_{\mathrm{SZ}}=\mathrm{PET}-\mathrm{AET}.\end{array}$$

For a shallow water table, upward capillary fluxes may also be taken from the groundwater reservoir. Thus, the rate of actual evapotranspiration from the SZ (AET_SZ), including both plant water uptake and capillary rise, is estimated as a linear function of the water table depth as follows:

$$\begin{array}{}\text{(37)}& {\mathrm{AET}}_{\mathrm{SZ}}=\max \left[{\mathrm{PAET}}_{\mathrm{SZ}}\left({\displaystyle \frac{h-z_{\mathrm{root}}}{D_{\mathrm{root}}}}\right)\mathrm{\Delta }t,\mathrm{0}\right].\end{array}$$

Figure 3 Schematic representation of (a) UZ–SZ interactions: (a1) indicates no UZ–SZ interaction, whereas (a2) indicates UZ–SZ interaction (soil depth, D_root, is reduced to D_uz). (b) surface water–groundwater (SW–GW) interactions in stream cells: boundary conditions change from no-flow to head-dependent flux conditions once the streambed or ground surface is intersected by the water table. The upper part of panel (b) shows the numerical implementation of SW–GW interactions in a stream cell.

2.5.2 Surface–groundwater interactions

Surface–groundwater interactions are characterised in DRYP through transmission losses, as described in Sect. 2.3.2. In addition, when the water table intersects a cell's defined streambed elevation it produces discharge into the stream, q_riv [L T⁻¹], and when the water table reaches the ground surface it produces saturation excess, q_s [L T⁻¹] (Fig. 3b) (Eq. 27).

Discharge into streams, q_riv, is quantified using a head-dependent flux boundary condition (similar to that used in MODFLOW Harbaugh, 2005):

$$\begin{array}{}\text{(38)}& q_{\mathrm{riv}}=C\left(h-h_{\mathrm{riv}}\right),\end{array}$$

where C is a conductance term [L² T⁻¹] estimated as

$$\begin{array}{}\text{(39)}& C={\displaystyle \frac{K_{\mathrm{ch}}{L}_{\mathrm{ch}}W}{\mathrm{0.25}\mathrm{\Delta }x}}.\end{array}$$

To avoid numerical instabilities, we use a regularisation approach implemented via a smooth switch between the flux boundary condition and a constant head boundary (and vice versa) using a convex function (Marçais et al., 2017):

$$\begin{array}{}\text{(40)}& q_{\mathrm{s}}=f_u\left({\displaystyle \frac{h-h_{\mathrm{b}}}{z-h_{\mathrm{b}}}}\right)f_g\left(\mathrm{\nabla }\cdot \left(-K_{\mathrm{sat}}h\mathrm{\nabla }h\right)+R-q_{\mathrm{riv}}\right),\end{array}$$

where h_b is the aquifer bottom elevation [L], and f_u is the continuous function between [0, 1]. f_u is specified as follows (Marçais et al., 2017):

$$\begin{array}{}\text{(41)}& f_u=\exp \left(-{\displaystyle \frac{\mathrm{1}-u}{r}}\right),\end{array}$$

where r is a dimensionless regularisation factor r > 0, which has been specified as 0.001 following Marçais et al. (2017). f_g is the Heaviside step function:

$$\begin{array}{}\text{(42)}& f_g=\left\{\begin{array}{l}\mathrm{0},u<\mathrm{0}\\ u,u\ge \mathrm{0}.\end{array}\right.\end{array}$$

After both q_s and q_riv are estimated, their corresponding volumes are estimated by multiplying the flow rate, the time step, and the corresponding surface area (cell or stream). The volume is then added as additional runoff in the surface component (Sect. 2.3.2). The water table is updated to its topographical elevation and kept as a constant head boundary condition. The boundary switches back to a flux condition if the water table drops back below the water table.

2.6 Numerical implementation and time step

DRYP is a fully open-source, grid-based model with a layer-based structure, developed using the Landlab architecture (Hobley et al., 2017) and its Python library. Landlab was chosen due to its versatility and modular design, the latter of which allows the user to plug in multiple modules for different levels of complexity and processes using grid-based objects (Hobley et al., 2017; Barnhart et al., 2020).

As most hydrological processes in DRYP, except the SZ component and the modified Green–Ampt infiltration, are described according to explicit analytical solutions, it is possible to run DRYP at hourly or sub-hourly time steps at a low computational cost.

The three main DRYP components (i.e. surface, unsaturated, and saturated components) can run at different time steps, from sub-hourly to daily. The riparian zone of the unsaturated component can also be run at a different time step to that of the unsaturated component. Where different time steps are used between components, the fluxes and state variables are temporally aggregated in DRYP by accumulating and/or averaging them over the specified time step as appropriate and then transferring them to the next component. In addition, and as described above, for the saturated component, an internal time step is also automatically considered to ensure the stability of the numerical solution.

3 Model evaluation methods

3.1 Evaluation using synthetic experiments

The use of synthetic experiments is an important aspect of model development in hydrology which is welcome but not often used (Clark et al., 2015). The objective of synthetic experiments is to better understand the structural controls on the physical processes represented in the model, for example, on groundwater–soil interactions (Rahman et al., 2019; Batelis et al., 2020). Here, we perform a set of numerical experiments to evaluate the stability and convergence of DRYP components, particularly the coupling of both the surface and unsaturated zone with the groundwater component. Convergence and stability of the numerical solution of the groundwater component using the FTCS finite difference approach and the regularisation have been well documented in different studies (e.g. Wang and Anderson, 1982; Anderson et al., 2015; Marçais et al., 2017). Hence, here we have considered two sets of model evaluations: (i) a quantitative evaluation of the model performance in relation to the well-known numerical model, MODFLOW, for a simple surface–groundwater interaction test represented as a draining condition (see Sect. 3.1.1), and (ii) a qualitative evaluation of the model performance with respect to the desired skill of the model to seamlessly allow interactions between groundwater and the land surface and surface water components (see Sect. 3.1.2).

3.1.1 Comparing DRYP and MODFLOW

For the quantitative evaluation, a 1-D synthetic experiment considering an inclined plane aquifer was set up using DRYP (see Fig. 4a). The length and width of the model domain were specified as 10 and 1 km respectively. Hydraulic saturated conductivity and aquifer specific yield were specified as 1.2 m d⁻¹ and 0.01 respectively. Boundary conditions were specified as no-flow for both the right and left side as well as the bottom of the model domain. The model grid size was set to 1 km × 1 km.

A model with identical geometry, grid size, and hydraulic properties was built in MODFLOW using the “FloPy” Python package (Bakker et al., 2016b, a). Boundary conditions for the MODFLOW model were the same as DRYP except for the top boundary condition, which was specified using the “drain” package (Harbaugh et al., 2000). The elevation at which the water starts to drain was specified as the top surface elevation of the model domain. A high value of the conductivity term (500 m² d⁻¹) was used in order to capture the seepage process and to ensure convergence as well as minimal water balance errors (Batelaan and Smedt, 2004).

Figure 4 Model domain for the synthetic experiments: (a) 1-D model, and (b) tilted-V catchment and flow boundary conditions specified for model simulations.

The synthetic test consisted of a free-draining condition for an unconfined aquifer with a water table depth equal to zero (at the surface level). The time step used for evaluation was 1 d. The evaluation considered the temporal variation in the water table for both DRYP and MODFLOW models, as well as the water balance errors. Errors were evaluated at all locations along the aquifer. Mass balance errors were estimated by the algebraic sum of inputs, outputs, and the storage change.

3.1.2 Qualitative analysis of surface–groundwater interactions

The geometry of the model domain for the qualitative tests consisted of a tilted-V catchment (Fig. 4) with a size of 7×10 square cells on a 1 km resolution grid. Land use and soil hydraulic characteristics were specified as uniform over the entire model domain, and the saturated zone was considered as a homogenous and unconfined aquifer. Boundary conditions were specified as no-flow boundaries for all sides as well as at the bottom of the model domain. The initial water table was set as a horizontal plane at the level of the catchment outlet (100 m) for all simulations (Fig. 4). For experimental purposes, hydraulic characteristics of both the unsaturated and saturated zone were arbitrarily chosen. Thus, a loamy sand soil texture with K_sat= 29.9 cm h⁻¹, θ_sat= 0.40, θ_fc= 0.175, and θ_wp= 0.075 was chosen for the unsaturated zone, whereas, for the saturated zone, the hydraulic conductivity of the aquifer (K_aq) was specified as 6 m d⁻¹, and the specific yield (S_y) was set as 0.01. The high value of K_aq combined with S_y and boundary conditions of the aquifer were applied in order to allow a fast increase/decrease in the water table and the observation of surface–groundwater interaction over a short period of time.

Three main scenarios were analysed using synthetic time series of precipitation and evapotranspiration and changing hydraulic parameters of the UZ:

1. an “infiltration–discharge” scenario, where all precipitation was allowed to infiltrate into the catchment and no infiltration excess was produced over the model domain;

2. an “infiltration–evapotranspiration–discharge” scenario was simulated by adding a time-variable potential evapotranspiration as input into the model;

3. an “infiltration–runoff–evapotranspiration–discharge” scenario was designed to evaluate the production of runoff and focused groundwater recharge, as well as groundwater discharge. For this last scenario, the saturated hydraulic conductivity of the soil was decreased by 1 order of magnitude to produce infiltration excess and, consequently, runoff.

For all three scenarios, precipitation events were specified at a constant value of 0.25 [mm h⁻¹] over 10 d followed by a 20 d dry period. Potential evapotranspiration was specified as a sinusoidal function with a 24 h period and a maximum rate of 0.10 [mm h⁻¹]. These experimental values of precipitation and evapotranspiration combined with the hydraulic properties of the unsaturated and saturated zone allowed a visual evaluation of surface–groundwater interactions under different conditions, such as increasing and decreasing water table through the model run and its interaction with the unsaturated zone.

3.2 Model evaluation based on observed catchment data at Walnut Gulch, USA

In addition to evaluating the DRYP model with synthetic experiments, the model was also evaluated at the Walnut Gulch Experimental Watershed (WGEW), a 149 km² basin near Tombstone, Arizona, USA (31^∘43^′ N, 110^∘41^′ W) (Fig. 5). The climate of the region is semi-arid with low annual rainfall: long-term average of 312 mm yr⁻¹ (Goodrich et al., 2008). The ephemeral channels of WGEW are comprised of mixed sedimentary beds that promote high transmission losses, leading to downstream declining discharge in all but the largest streamflow events (Singer and Michaelides, 2014; Michaelides et al., 2018). WGEW was chosen because it has a long and spatially explicit record of runoff (Stone et al., 2008) for multiple flumes as well as high-density event-based rainfall data for 95 operational gauging stations (Goodrich et al., 2008) which were used to analyse trends in rainfall characteristics (Singer and Michaelides, 2017) and from which the STORM model was created (Singer et al., 2018). In addition, water content from a cosmic-ray neutron sensor as well as latent heat flux from an eddy-covariance flux tower are also available in the basin (Emmerich and Verdugo, 2008; Zreda et al., 2012). Together, these data provide an ideal opportunity to assess many components of a model of dryland water balance and partitioning (Emmerich and Verdugo, 2008; Goodrich et al., 2008; Keefer et al., 2008; Stone et al., 2008; Scott et al., 2015). The water table at WGEW is deep (∼ 50 m at the catchment outlet), so the potential interaction between surface and groundwater is generally unidirectional (Quichimbo et al., 2020).

Figure 5 Geographic location of the Walnut Gulch Experimental Watershed and the location of monitoring stations.

3.2.1 Model setting, inputs, and parameters

For WGEW, model simulations were performed using the modified Green–Ampt infiltration approach because of its ability to describe the high potential infiltration rates at the beginning of the precipitation event, which are particularly important in this setting. The time step of both the surface component and unsaturated component was specified as 1 h, whereas we use a time step of 1 d for the riparian zone to reduce computational time. The high temporal resolution for the unsaturated component was used to capture the observed high-intensity, low-duration rainfall at WGEW, as well as the influence of diurnal fluctuations in evapotranspiration. As the water table is deep below the ground surface and surface water–groundwater interactions are known to be limited, the groundwater component was not included in model simulations for WGEW.

Spatial and temporal information required as inputs and for model parameters were obtained for WGEW from https://www.tucson.ars.ag.gov/dap/ (last access: 15 October 2021). A model domain of 104×41 square cells on a 300 m resolution grid was developed. A digital elevation model with a spatial resolution of 30 m×30 m was obtained from the Shuttle Radar Topography Mission (SRTM) 1 Arc-Second Global map (available at https://earthexplorer.usgs.gov, last access: 15 May 2021). The DEM was aggregated by averaging cells to the 300 m grid size. Textural characteristics of soil and land cover, obtained as polygon files from https://www.tucson.ars.ag.gov/dap/ (last access: 10 October 2021), were converted into model gridded inputs by considering the feature with the biggest area as the raster value. Based on the soil texture, baseline hydraulic properties required for modelling were obtained from Rawls et al. (1982) and Clapp and Hornberger (1978). Values of field capacity and wilting point required to estimate L_TAW (Fig. 2a) were obtained assuming a matric potential of −33 and −1500 kPa following the FAO guidelines (Walker, 1989).

Stream positions were estimated from the 30 m×30 m DEM. The routing network at the 30 m×30 m grid resolution was specified by defining a minimal upstream drainage area threshold of 65 ha, which corresponds to the medium stream network resolution specified in Heilman et al. (2008). Stream cells were then aggregated to the model grid size, 300 m × 300 m, to obtain the stream length at any given cell. Stream width was assumed to be 10 m for the whole model domain based on average values observed across the whole catchment (Miller et al., 2000). Point measurements of rainfall were obtained from 95 rainfall stations that are well distributed within the basin (Goodrich et al., 2008). Rainfall data, at every location, were temporally aggregated to 1 h and then spatially interpolated using a natural neighbour algorithm to a 30 m ×30 m grid size to preserve the high spatial and temporal variability of stations located at distances smaller than the model grid size. Finally, rainfall was spatially aggregated to the grid size of the model domain. Potential evapotranspiration (PET) was calculated using hourly data from the ERA5-Land reanalysis (Hersbach et al., 2020) because this data set enabled high temporal resolution (1 h) and has the potential to drive hydrological and land surface models (Albergel et al., 2018; Alfieri et al., 2020; Tarek et al., 2020; Singer et al., 2021). Data from ERA5 have a spatial resolution of ∼ 9 km at the Equator. The Penman–Monteith approach was chosen to estimate hourly PET due to its high accuracy with respect to producing evapotranspiration values under different climates and locations, and also because it is considered a standard method by the FAO (Allen et al., 1998).

High-resolution temporal measurements of runoff at three flumes (F01, F02, and F06) along the main Walnut Gulch channel were used in the evaluation of runoff generation (Fig. 5). To evaluate modelled soil moisture, we used data from a cosmic-ray neutron sensor station from the COSMOS network (Zreda et al., 2012), located within the Kendall subcatchment of WGEW (Fig. 5). The raw data (publicly available at http://cosmos.hwr.arizona.edu/Probes/StationDat/010/index.php, last access: 20 June 2021) were corrected for atmospheric pressure (Desilets et al., 2010), atmospheric vapour pressure (Rosolem et al., 2013), above-ground biomass, and variation in background intensity using the standardised data processing Cosmic-Ray Sensor PYthon tool (Power et al., 2021) for the period between mid-2010 and 2018.

Finally, data from the AmeriFlux site Kendall Grassland (US-Wkg; available at https://ameriflux.lbl.gov/sites/siteinfo/US-Wkg, last access: 20 June 2021), were used for the evaluation of simulated AET (Fig. 5). Uncertainty in flux tower data is mainly attributed to instrumental and random errors, and it increases with flux magnitude (Richardson et al., 2006; Schmidt et al., 2012). Mean relative errors for AmeriFlux sites are around −5 % with deviations of ±16 % (Schmidt et al., 2012). Historical records from mid-2006 to 2018 were available for model evaluation.

3.2.2 Model sensitivity analysis and calibration

An initial trial-and-error calibration of the model was performed to explore the parameter sensitivities of DRYP and to reduce the a priori parameter ranges used in the second step. This first trial-and-error calibration considered only the performance of the model to represent streamflow at the catchment outlet (flume F01). The calibration was performed by applying spatially constant multiplicative factors kW, kK_sat, kD_root, kK_ch, and kk_T, to model parameters W, K_sat, D_root, K_ch, and k_T respectively. These parameters were used because they control the storage and the water partitioning of components (surface and subsurface components) in the DRYP model for WGEW. Parameters W, K_ch, and k_T were assumed to be uniform over the entire catchment due to the lack of spatial information, whereas the rest of parameters listed in Table 1 vary depending on their mapped spatial distribution. The initial manual calibration enabled a set of parameter ranges to be defined for a Monte Carlo experiment to analyse the multi-parameter uncertainty of the model results. Then, a set of 1000 realisations was implemented for the analysis with parameters randomly generated using a uniform distribution for each parameter.

The Generalized Likelihood Uncertainty Estimation (GLUE) framework (Beven and Binley, 1992) was used as the uncertainty analysis framework. The GLUE framework considers that, owing to the uncertainty in the input data, model structure, and limitations of boundary condition, there are multiple set of parameters that can produce acceptable simulations. To determine which simulations were considered acceptable (i.e. behavioural), we used a combination of two different “goodness of fit” indices: the Nash–Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) and the per cent bias (PBIAS), which are defined as follows:

$$\begin{array}{}\text{(43)}& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(O_i-S_i\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}}};\text{(44)}& {\displaystyle }\mathrm{PBIAS}\left(\mathit{\%}\right)=\mathrm{100}\cdot {\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{O}_i-{\sum }_{i=\mathrm{1}}^n{S}_i}{{\sum }_{i=\mathrm{1}}^n{O}_i}}.\end{array}$$

Here, O represents the observation, $\stackrel{\mathrm{‾}}{O}$ is the arithmetic mean of observations, S represents the model simulations, and n is the number of observations.

In order to define behavioural models, a set of thresholds was specified for the three indices. For streamflow, NSE values higher than 0.50 and PBIAS values less than 20 % (i.e. less than 1 % of the total water budget of the study area) were considered as acceptable simulations. For soil moisture and actual evapotranspiration, only NSE values greater than 0.5 were accepted.

In order to combine these measures into a single performance metric, models that did not meet these conditions were assigned a value of 0, whereas the indexes were linearly scaled between 0 and 1 for the rest of the models. Scaling of NSE values was performed according to the following range: 0 for the minimum value (NSE  = 0.5) and 1 for the maximum value of NSE which is also 1. For PBIAS, absolute values were scaled by considering the maximum value (PBIAS = 20 %) equal to 0 and the minimum (PBIAS  = 0) equal to 1. The combined performance measure was calculated as the product of all indexes considered in the analysis:

$$\begin{array}{}\text{(45)}& p_i={\prod }_{k=\mathrm{1},\mathrm{2},\mathrm{6}}N{S}_k^{\ast }\cdot {\mathrm{PBIAS}}_k^{\ast },\end{array}$$

where p is the combined performance measure for the ith parameter set, the ^* signifies scaled values, and k represents the variable considered in the analysis.

For soil moisture, a direct comparison between observations and simulations was not possible due to differences between the representative soil depths of measurements and simulations. Modelled soil moisture represents the water content of the entire soil column specified by the rooting depth, whereas the observed soil moisture represents the water content over a depth-averaged value, which can be characterised by an effective soil depth that depends on the soil moisture itself (Franz et al., 2012). A direct comparison would result in the misrepresentation of high values of observed soil water content by the model due to the attenuation of peak values over larger soil depths. This problem has been solved by using exponential models that need to be calibrated using measurements at different soil depths (e.g. Wagner et al., 1999; Albergel et al., 2008). Therefore, to enable model–data comparisons that capture the variation in both high and low values of soil moisture observations, we scaled observed soil moisture using the following expression:

$$\begin{array}{}\text{(46)}& O^{*}=O^{\mathit{\alpha }}+S_{\min },\end{array}$$

where ^* refers to the scaled value, and α is estimated by

$$\begin{array}{}\text{(47)}& \mathit{\alpha }={\displaystyle \frac{\log \left|S_{\max }-S_{\min }\right|}{\log \left|O_{\max }-O_{\min }\right|}}.\end{array}$$

The period between 1 January 2007 and 1 January 2018 was the temporal domain for model simulations at WGEW, with a warm-up period of 1 year prior to this period. This period matches the overlapping period of streamflow observations and flux tower observations. Soil moisture was evaluated for a shorter period of available data between 1 October 2010 and 1 January 2018. Additionally, modelled soil moisture for COSMOS site was obtained by spatial averaging the nine cells located around the COSMOS station to match the typical COSMOS footprint diameter (∼ 700 m) (Desilets and Zreda, 2013; Schrön et al., 2017).

4 Results

4.1 Comparing DRYP to MODFLOW

The modelling results show a good agreement between DRYP and MODFLOW, and both models ran with negligible mass balance errors ($\mathrm{1.79}\times {\mathrm{10}}^{-\mathrm{15}}$ m³ for DRYP, and $\mathrm{6.95}\times {\mathrm{10}}^{-\mathrm{8}}$ m³ for MODFLOW) (see Fig. 6). Both the temporal and spatial variations in the MODFLOW model are well captured by DRYP. Differences in water table elevations are in the range of just 0.022 m at the beginning of the simulation when the aquifer starts to drain, and these differences only decrease as the water table decreases (see Fig. 6a). Temporal variation in the groundwater storage for both models shows consistency, with higher values due to high gradients at the beginning of the simulation, decreasing as the water table declines (see Fig. 6b). More fluctuations are observed in the MODFLOW simulations, which can be attributed to the time step used for the simulation, which needs to be reduced in order to smoothly capture the variation in water table depth as the model switches boundary conditions (Eq. 41). DRYP captures this variation smoothly due to the exponential function and by automatically reducing the time step to assure numerical convergence.

Figure 6 (a) Simulated head along the aquifer for different time steps (in months, represented using “M” in panel a) estimated by DRYP (solid lines) and MODFLOW (dashed lines), and (b) temporal variation in the mass balance error for DRYP.

Figure 7 Temporal variation in (a) precipitation (black line) and evapotranspiration (grey line), (b) water content of the unsaturated zone, (c) groundwater recharge, (d) runoff/discharge, and (e) water table elevations. The right panels represent zoomed-in sections of the shaded areas in the left panels. Solid lines represent the variation at the catchment outlet, whereas dashed lines represent the temporal variation in the stream at a distance of 4 km from the catchment outlet. For panels (b)–(e), blue lines represent the “infiltration–discharge” scenario, green lines represent the “infiltration–evapotranspiration–discharge” scenario, and red lines represent “infiltration–runoff–evapotranspiration–discharge” scenario.

4.2 Evaluation of qualitative synthetic experiments

Figure 7 shows the temporal variation in fluxes and state variables for the three simulated scenarios at two evaluation points located along the channel, one at the catchment outlet and the second 4 km from the catchment outlet. The respective scenarios' results are described below:

1. In the “infiltration–discharge” scenario (blue lines), precipitation immediately infiltrates into the unsaturated zone when it falls over the catchment (Fig. 6a), increasing the water content of the soil. As there are no losses due to evapotranspiration, the water content steadily increases until it reaches field capacity (Fig. 7b). At field capacity, given that the soil cannot hold any excess water, it starts to release water as diffuse recharge. The soil remains at field capacity for the rest of the simulation, allowing the water from the next rainfall event to move directly to the saturated zone producing recharge (Fig. 7c). Recharge produces an increase in groundwater storage and, consequently, increases the discharge at the outlet of the catchment (Fig. 7d). In the early precipitation events, the contribution of groundwater discharge is minimal. However, this contribution keeps increasing until a dynamic steady state is eventually reached (by 14 600 h, not shown in Fig. 7). Discharge closely follows the temporal variation in the precipitation, due to the high transmissivity of the aquifer and the saturation of the soil; a sharp increase in discharge means that precipitation has become the main contributor to discharge changes because the water table is at the surface (Fig. 7d).

2. In the “infiltration–evapotranspiration–discharge” scenario (green lines), the addition of evapotranspiration in this experiment produces a reduction in soil water content (Fig. 7b). As precipitation is much higher than evapotranspiration, soil moisture quickly reaches a dynamic steady state at the end of the second precipitation event. At cells located close to the catchment outlet, the rise of the water table to ground level reduces the thickness of the unsaturated zone to zero, and as no water can be infiltrated, the soil water content is kept at its highest value during the precipitation event. After the precipitation event, the rate of evapotranspiration, which is greater than the rate of lateral groundwater inflow, gradually reduces the amount of water in the cell. However, as the storage of the SZ keeps increasing, the thickness of the UZ decreases and the rate of lateral groundwater flow becomes greater than the rate of evapotranspiration; this also results in quick changes in the water content of the soil (Fig. 7b). Recharge is also reduced and, as expected, only occurs when the soil moisture reaches field capacity (Fig. 7c). Discharge is also reduced as a result of decreased aquifer recharge due to upward losses by AET. For cells close to saturation, the storage in the groundwater reservoir is affected by evapotranspiration losses (not observed in right panel due to the y axis scale), which in turn results in daily fluctuations in discharge that are the inverse of evapotranspiration fluctuations.

3. In the “infiltration–runoff–evapotranspiration–discharge” scenario (red lines), a reduced K_sat results in the development of infiltration–excess overland flow. The rate of infiltration at the beginning of the precipitation event is high enough to provide water for evapotranspiration without reducing the soil water storage (Fig. 7b), which explains the similarity in soil moisture behaviour with the second scenario. When cells start to produce runoff as a result of infiltration excess, discharge also starts to rise. At stream cells with a deep water table, the increase in streamflow is the result of flow accumulation along the channel during the precipitation event (e.g. Fig. 7d, left panel at 6600 to 6700 h). At cells where the water table interacts with the surface, groundwater discharge gradually increases the streamflow at the catchment outlet at much longer temporal scales (Fig. 7d). At the catchment outlet, streamflow is also affected by the fluctuation of the water table due to the daily variation in evapotranspiration losses (Fig. 7d).

Figure 8 Cumulative volume of the main components of the water balance for the simulated scenarios: (a) infiltration–discharge, (b) infiltration–evapotranspiration–discharge, and (c) infiltration–runoff–evapotranspiration–discharge. P is the precipitation, R is recharge, Q is discharge at the catchment outlet, AET is actual evapotranspiration, GWS is the change in groundwater storage, and Err is the water balance error of the simulation.

Figure 8 shows the cumulative volumes of different components of the water balance as well as the cumulative mass balance error of the model. Mass balance errors are low in comparison to the total amount of water entering the catchment, with values of less than 0.12 % for the first case (only precipitation). For the other two cases where evapotranspiration is included, errors are less than 0.02 %. The higher error for the first scenario is attributed to the concentration of flow at the catchment outlet, which leads to an increase in the number of cells discharging into the surface and the channel as well as the resulting minor numerical artefacts.

Coupling of surface and groundwater processes often results in numerical instabilities and in convergence problems (Batelaan and Smedt, 2004; Marçais et al., 2017). However, the results of these synthetic experiments illustrate DRYP's ability to produce realistic hydrological process behaviours by providing a stable solution for representing surface–groundwater interactions without producing numerical artefacts. DRYP is effective at handling the complex coupling and dynamic switching of different types of hydraulic boundary conditions, producing acceptable results with negligible mass balance errors.

4.3 Model performance at WGEW

4.3.1 Spatio-temporal visualisation of the model process simulation at WGEW

The ability of the model to capture the dynamics of dryland hydrological processes is illustrated for WGEW in Fig. 9. The best model (see following section) captures the emergence of ephemeral flow conditions for specific storms, as well as the spatio-temporal changes in soil moisture. It can be seen how, for a given initial soil moisture condition, the production of runoff due to a rainfall event falling over only the central part of catchment results in the concentration of flow along the stream. As water moves downstream, the stream loses water due to transmission losses, which ultimately consumes almost all of the available water by the time runoff reaches the catchment outlet (flume F01 in Fig. 9c).

Figure 9 Spatio-temporal visualisation of the model process simulation at WGEW: (a) the rainfall event, (b) soil moisture prior to the rainfall event, (c) the ephemeral stream for the rainfall event, and (c) soil moisture after the rainfall event. The x and y axis distance units are in metres.

4.3.2 Characterisation of the temporal variation in simulated variables

Calibration using the trial-and-error method, showed that streamflow displayed particular sensitivity to the parameters K_sat, D_root, K_ch, and k_T. This informed a set of parameters ranges that were used in the Monte Carlo analysis as follows: for hydraulic conductivity at the channel, kK_ch, the range varies between 0.10 and 0.30; for kK_sat, the range varies between 0.20 and 0.50; for kk_T, the range varies between 3 and 10; and for kD_root, the range varies between 0.80 and 1.20. This resulted in 21 behavioural models with values of p above zero. The calibrated parameters for the best simulation were kK_ch= 0.21, D_root = 1.02, kK_sat= 0.30, and kk_T= 7.7. A factor of kk_T= 7.7 applied to a default value of k_T (0.083) represents a flow velocity of 0.41 m s⁻¹ in the channel.

Soil moisture

The DRYP model demonstrates skill at capturing the dynamics of the soil moisture (Fig. 10b) with NSE values of around 0.69. Discrepancies in the magnitude of peak values are likely the result of scaling, so simulations are not able to account for the variation in the effective measurement depth of COSMOS water content estimates (Franz et al., 2012, 2013). The effective COSMOS measurement depth is greater for low values of soil water content (around 33 cm), whereas the effective measurement depth is shallow for higher values of water content (around 16 cm). However, discrepancies may also reflect the limited ability of the soil moisture model to represent high variations occurring at shallow depths of the soil layer, due to the use of a single store.

Evapotranspiration

The DRYP model also skilfully captures (NSE ∼ 0.7) the seasonality and the overall temporal variation in evapotranspiration, a dominant component of the water budget in drylands (Fig. 10c), although peak values are generally overpredicted after long periods of dry conditions. Nevertheless, discrepancies between flux tower data and simulated AET up to 15 % for 1 year have been reported for grassland vegetation in previous studies (Twine et al., 2000; Scott, 2010), and such errors are mainly attributed to the inherent uncertainty in rainfall and latent heat flux measurements (Scott, 2010).

Figure 10 Comparison between observed and simulated values of monthly temporal variation (a–f) and monthly distribution (b–f) of (a) monthly precipitation (left axes) and yearly precipitation (right axes), (b) soil moisture at the COSMOS Kendall location, (c) actual evapotranspiration at Kendall, (d) streamflow at flume F06, (e) streamflow at flume F02, and (f) streamflow at flume F01. See Fig. 4 for station locations.

Streamflow

DRYP is also able to reproduce the seasonality and the monthly production of runoff at the outlet of the catchment (F01, NSE ∼ 0.9) (Fig. 10f), as well as at the two upstream flumes (F02, F06) considered in the analysis (NSE > 0.60) (Fig. 10d, e). However, monthly values at flumes F02 and F06 are overpredicted in 2012, perhaps reflecting the development of a crusting layer in previous dry years (e.g. 2009, 2011), a process not included in the model. On the other hand, low production of runoff during wet years (e.g. 2015) may be attributed to the energy of high-intensity rainfall events removing such a crusting layer from the top of the soil, which would result in the increase in infiltration rates (Becker et al., 2018). Additionally, the spatial aggregation of the DEM causes slight inaccuracies in the estimated contributing areas for different streams. This affects not only the volume but also the timing of streamflow events, which may result in over-/underprediction of streamflow events and may ultimately affect the overall water budget.

Water balance

Precipitation shows high annual variability for the evaluated period, with the lowest value of 200 mm yr⁻¹ in 2011 up to 400 mm yr⁻¹ in 2015 (Fig. 10a), which translates into variability in the annual water partitioning for WGEW. For the evaluation period, 1 January 2007 to 1 January 2018, water balance estimates from the best model show that ∼ 92 % of the total precipitation infiltrates into the soil (see Fig. 11). However, almost all infiltrated water returns to the atmosphere as evapotranspiration losses, representing 89 % of the total precipitation. A small proportion, ∼ 3 % of the total precipitation, remains in the soil and the riparian zone, and this stored water corresponds mainly to wetter years of the simulation period (2014 and 2015). Only a small percentage, less than 0.03 %, percolates as diffuse recharge contributing to groundwater storage. Water that does not infiltrate into the soil (8 % of the precipitation) is routed downstream. However, this amount of water is consistently reduced by transmission losses, representing ∼ 7 % of the precipitation. Water entering the riparian zone via transmission losses is partitioned into evapotranspiration and focused recharge. Evapotranspiration consumes up to 60 % of these transmission losses, representing ∼ 4.5 % of the total precipitation. This is broadly consistent with previous studies showing values of 20 mm yr⁻¹ or 5.5 % to 7 % of the total precipitation (Renard, 1970; Renard et al., 2008). The amount of surface water leaving the catchment represents less than 1.0 % of the total amount of precipitation falling over the catchment. These values highlight the impact of transmission losses on the streamflow and aquifer recharge. The main contributor to the total amount of groundwater recharge is focused recharge (∼ 2.5 % of precipitation).

Figure 11 Average fluxes of different components of the water budget of WGEW for the simulated period, between 1 January 2007 and 1 January 2018. Water stored in the soil and the riparian zone corresponds to ∼ 3 %. Blue arrows show input fluxes, green arrows represent water leaving the catchment, orange arrows represent internal surface and unsaturated zone fluxes, and yellow arrows represent water moving to the saturated zone (groundwater flow not modelled).

5 Conclusions

We have developed and presented a parsimonious model to estimate water partitioning in dryland regions (“DRYP”). DRYP is parsimonious in that the model structure has few tunable parameters but still captures the essential elements of dryland hydrology and skilfully represents hydrological stores and fluxes in drylands. It is designed to allow the user to identify mechanisms and factors that affect the water balance in dryland regions, where most existing hydrological models may not be able to capture them appropriately (e.g. the role of focused versus diffuse discharge). We have provided a technical description of all components of DRYP and evaluated it under different scenarios. We first evaluated the ability of DRYP to provide stable numerical simulations of the interaction of surface and subsurface components through synthetic model experiments. We then evaluated DRYP using streamflow, soil moisture, and evapotranspiration data from the semi-arid Walnut Gulch Experimental Watershed (Arizona, USA). We tested the ability of the model to produce behavioural simulations based on multi-parameter Monte Carlo experiments evaluated against a range of objective performance metrics. A comparison between DRYP and MODFLOW for a simple draining-case model showed excellent agreement with an error in hydraulic head of < 0.022 m and a mass balance error of $\mathrm{6.95}\times {\mathrm{10}}^{-\mathrm{8}}$ m³. An evaluation of surface–groundwater interactions using numerical experiments over a synthetic model domain indicated that DRYP shows skill at producing stable simulations for the main components of the water balance with low mass balance errors (< 0.12 %). Thus, we conclude that DRYP has the potential to be robustly applied in environments where surface–subsurface interactions play an important role in the overall mass balance of the catchment.

For Walnut Gulch, DRYP effectively captures the spatio-temporal variability of the main components of the dryland water balance at monthly timescales. We find that focused recharge represents ∼ 2.5 % of the total amount of rainfall, whereas diffuse recharge is below 0.03 %. Evapotranspiration is the dominant process representing 90 % of the water leaving the catchment. Evapotranspiration from riparian areas also plays an important role in groundwater recharge, as the amount of water becoming focused recharge is only around ∼ 40 % of the transmission losses.

Considering the combination of explicit solutions of surface and subsurface components, the parsimonious structure, and the low computational cost, it is possible for DRYP to perform long runs using hourly or sub-hourly time steps. These characteristics enable DRYP to test long-term and seasonal changes in water availability to plants and humans in limited-water environments under different scenarios and future climatic conditions such as anthropogenic activities or during droughts. Additionally, given the minimal data requirements, DRYP has the potential to be applied in areas where only information at large scales is available.

Furthermore, improving the soil–vegetation interaction in the unsaturated zone to capture the temporal variation in plant water demand will likely enhance the performance of the model. A more complex representation of the highly dynamic behaviour of ephemeral streamflow will be considered in future developments in order to enhance the ability of the model to represent flooding conditions. Additionally, we also plan to use DRYP in conjunction with stochastic rainfall simulation tools (such as STORM; Singer et al., 2018) to explore the impact of the variability of precipitation on the water balance.

Appendix A

Table A1 List of parameters and model variables.