Simulating water flow in variably saturated soils: a comparison of a 3D model with approximation-based formulations

An increasing number of hydrological models simulate variably saturated flow solving the governing partial differential equations in a three-dimensional domain (Panday & Huyakorn 2004; Loague et al. 2005; Maxwell & Miller 2005; Rigon et al. 2006; Qu & Duffy 2007; Simunek et al. 2008; Kumar et al. 2009; Camporese et al. 2010; Anagnostopoulos & Burlando 2012; Mirus & Loague 2013; Anagnostopoulos et al. 2015). While such an approach is theoretically consistent, its practical advantages when compared with simpler one-dimensional models and simplified subsurface flow routing algorithms are still under debate. Specifically, fully explicit three-dimensional (3D) models based on the Richards equation are computationally demanding but allow a theoretically sounding solution for spatially complex domains, heterogeneous soils, and solute transport, among other advantages (Kollet & Maxwell 2006; Simunek & van Genuchten 2008; Maxwell et al. 2014). Conversely, simplified subsurface flow algorithms provide faster run times and are generally versatile and less prone to numerical instability in the treatment of extreme conditions, such as very dry soils or heterogeneous soil properties (Tocci et al. 1997). However, they suffer from assumptions adopted to simplify the underlying physical problem (Downer & Ogden 2004; Hilberts et al. 2007; Ivanov et al. 2010; Shen & Phanikumar 2010; Fatichi et al. 2012a). Advantages and disadvantages of models with different complexity have led to numerous approaches reported in scientific literature (Singh & Woolhiser 2002; Loague et al. 2005; Kampf & Burges 2007; Ebel et al. 2009; Ogden et al. 2015), providing reasonable estimates of flow and soil moisture in real hillslopes and catchments (Ivanov et al. 2004a, 2004b; Ebel et al. 2008; Jones et al. 2008; Sciuto & Diekkruger 2010; Fatichi et al. 2012b; Cordano & Rigon 2013). The existence of a wide range of approaches calls for studies that elucidate model assumption hierarchy in terms of impacts on practical applications.

Models represent spatial connections, hydraulic head gradients, and mass flows within hydrologic systems with different degrees of detail. It is not clear a priori in which configuration or under which boundary conditions one can simplify the simulation problem and still produce adequate estimates of subsurface flow, runoff generation, and soil moisture distribution. In this context, it is important to ask at which scale solving the nonlinear dynamics of subsurface flow with a fully 3D model becomes essential. Undoubtedly, hydrological models need to represent subsurface flow processes from the plot and hillslope to the catchment scale accurately enough to reproduce the physical behavior of the system. At the same time, however, excessive computational demands or issues of numerical instabilities need to be avoided. As exact analytical solutions involving both unsaturated and saturated zones are unavailable, relevant insights can be obtained from comparisons of models that simulate the same quantities but exhibit different degrees of complexity. Specifically, we target models that maintain representations of theoretically measureable states and quantities; in layman terms of the community, these are ‘spatially-distributed, physically-based/process-based’ models. The approach adopted here aims at shedding light on a hierarchy of model assumptions in order to establish the relative importance of process representation in the hydrological response at the hillslope scale. This is different from analyzing differences induced by the adopted numerical schemes and types of coupling between surface and subsurface processes, which were addressed in previous intercomparison studies of models with similar assumptions in the process-representation (Sulis et al. 2010; Maxwell et al. 2014).

In this study, we use three different models with a hierarchy of decreasing complexity to describe variably saturated flows: (1) the HYDRUS-3D model solves the complete three-dimensional formulation of the Richards equation (Simunek et al. 2006); (2) Real-time Integrated Basin Simulator + VEGetation Generator for Interactive Evolution (tRIBS + VEGGIE) implements a quasi-3D framework based on the one-dimensional vertical Richards model and the Boussinesq equation under the Dupuit–Forchheimer assumptions for lateral flows in the saturated zone and the kinematic approximation for lateral flows in the unsaturated zone (Ivanov et al. 2008, 2010); and (3) Tethys-Chloris (T&C) also uses a quasi-3D framework and 1D vertical Richards equation but relies on the kinematic approximation to resolve both saturated and unsaturated lateral flow exchanges (Fatichi et al. 2012a, 2012b). These two flow approximations were chosen as they represent commonly used approaches for describing water flow in hydrological models. As a default condition, we assume that solving the 3D Richards equation capable of describing flows for a range of geometric and soil heterogeneities represents currently the best possible numerical approximation of the physical system and thus HYDRUS-3D is considered to be the virtual reality. This implies that we neglect preferential flow paths and non-equilibrium flow dynamics, the treatment of which would require an additional layer of complexity (Gerke & van Genuchten 1993; Ross & Smettem 2000; Simunek et al. 2003; Beven & Germann 2013).

A thorough model comparison demands the identification of hillslope geometries, soil properties, and initial conditions that could reveal specific differences among methods resolving subsurface flow. For instance, Beven (1981) showed that the kinematic wave formulation for subsurface flow is a reasonable approximation of the extended Dupuit–Forchheimer assumptions for values of the dimensionless index less than 0.75, where β is the slope angle, K_s is the saturated hydraulic conductivity, and p is the steady flow input per unit area. Steep slopes, high saturated hydraulic conductivities, and low input rates all tend to reduce the value of λ and support the kinematic wave approximation (Beven 1981, 1982a, 1982b). It is therefore reasonable to expect that the kinematic method (i.e. the case of the T&C model) of subsurface exchange would fail for relatively shallow slopes and for wet conditions, where the transient effects of developing water table may be more pronounced and flow will not follow topographic gradients. While theoretical limitations of the Dupuit–Forchheimer assumptions have been investigated in detail for a set of specific conditions (Polubarinova-Kochina 1962; Troch et al. 2003; Hilberts et al. 2004, 2005), most of these were cases with a simple geometry. It is less clear whether the Dupuit–Forchheimer assumption can sufficiently approximate 3D variations of fluxes in domains of complex topography (i.e. the case of the tRIBS + VEGGIE model). Generally, one can expect that flow approximation-based models will agree with 3D numerical solutions in the case of simple geometries and high-to-moderate slopes but will increasingly disagree as the flow domain slope becomes either very shallow or very steep. One can also expect a better agreement for isotropic medium as opposed to cases with spatial variations of soil properties (e.g. soil layering). The outcome of increasing the spatial dimension of the flow domain is less certain, although one could expect a possible accumulation of errors associated with flow velocity approximations and the neglect of the vertical flow component, especially for steeper aquifer conditions.

A detailed comparison of the consequences of different assumptions for solving subsurface flow in a wide range of conditions has not been performed so far and this numerical study offers a first set of assessments. Specifically, this research explores a range of conditions using numerical experiments corresponding to event-based simulations at the scale of a hillslope–zero-order catchment that exhibits a certain degree of complexity in the organization of surface slopes. We limit the simulations to the event timescale in order to be able to focus solely on the representation of the water flow problem since long-term simulations would require including evaporation and transpiration processes. However, the way evapotranspiration is implemented and parameterized in the analyzed models varies greatly and would undoubtedly influence simulated moisture dynamics. Therefore, it would be impossible to separate the effects on water flow caused by model-specific evapotranspiration fluxes from those caused by the differences in flow descriptions.

We examine the effects of slope angle, domain size, rainfall magnitude, and initial moisture conditions on runoff generation and simulated soil moisture dynamics. We explore whether the different model strategies are competing approaches and whether the hierarchy in model simplifying assumptions outlined above is in fact reflected in the results. Using a pure modeling study we rely on previous knowledge of hillslope hydrologic response processes. However, the study advances the understanding of how different formulations of flow processes respond in complex hillslope domains and provides guidance for model development addressing the following specific research questions: (i) How do different subsurface flow routing methods compare? (ii) Under which conditions do simpler models agree with the 3D formulation? (iii) Are there sets of conditions in which a 3D model offers absolute superiority?

In this study, we used three different numerical models to describe the dynamics of variably saturated flow at the hillslope scale: HYDRUS-3D solves the 3D Richards equation; T&C and tRIBS + VEGGIE use approximations of varying complexity.

In contrast to models such as HYDRUS-3D that simultaneously compute fluxes and flow directions for all locations as a result of solving the governing system of equations, simplified models, such as tRIBS + VEGGIE and T&C, approximate the three-dimensional nature of a domain with 1D elements (cells) and determine water flow directions among the cells, often using the topographic gradients inferred from digital elevation data. The delineation of flow directions is an important step in the characterization of hydraulic connections among computational elements. Numerous methods have been proposed and a substantial body of literature has dealt with the issue of flow directions in the past decades (O'Callaghan & Mark 1984; Quinn et al. 1991; Tarboton et al. 1991; Costa-Cabral & Burges 1994; Tarboton 1997; Orlandini et al. 2003; Seibert & McGlynn 2007; Orlandini & Moretti 2009).

The flow directions are commonly determined using single and multiple flow direction algorithms (Nardi et al. 2008; Orlandini & Moretti 2009). The difference is in the partition of the flow from the upslope cell to the neighboring cells with lower elevations. In single flow direction methods, all of the flow is concentrated toward a single downstream cell, whereas multiple flow direction algorithms subdivide the flow among several cells, at least two.

The simplest flow direction method is called D-8. It identifies a single downstream direction among eight adjacent cells. The receiving cell is the one for which the steepest slope is found (O'Callaghan & Mark 1984). The most common multiple-direction method, D-mult (Quinn et al. 1991), attempts to solve the major limitation of the D-8 algorithm, the one-dimensional representation of flow, by weighting the mass flow among all lower neighboring cells in proportion to the actual slope gradient to these cells.

From an extensive review of the literature, Nardi et al. (2008) conclude that single flow direction methods are incapable to efficiently simulate flow on hillslopes and that multiple flow direction methods cannot accommodate concentrated channel flow; the latter often produce an excessive dispersion of flow, which may be inconsistent with the physical drainage (Orlandini et al. 2003).

A reasonable compromise between the simplicity of the single flow method and the sophistication introduced in multi-flow formulations was proposed by Tarboton (1997). The method, usually referred to as D-infinity (D-inf), selects a flow direction as the steepest downward slope among eight triangular facets formed by the elevation field of a cell under consideration and its surrounding cells, in both cardinal and diagonal directions. The partition of mass flow among the downstream cells (one or two, as maximum) is then calculated according to how close the selected flow direction is to the nearest cardinal and diagonal flow directions. A further evolution that reduces D-inf to a single flow direction has been proposed by Orlandini et al. (2003) and a combination of the single and multiple flow directions in a morphologically significant manner has been developed by Orlandini & Moretti (2009).

A coupled dynamic model of vegetation-hydrology interactions tRIBS + VEGGIE, the TIN (Triangulated Irregular Networks) based Real-time Integrated Basin Simulator + VEGetation Generator for Interactive Evolution (Ivanov et al. 2004a, 2004b, 2008), is used in this study. The model mimics principal water and energy processes over complex topography of a river basin. Computational elements are Voronoi cells. The model includes a detailed description of rainfall interception, evapotranspiration, infiltration with continuous soil moisture accounting, lateral moisture transfer in the unsaturated and saturated zones, and runoff routing. Generation of runoff occurs through different mechanisms. Energy and mass budgets are computed for a single canopy layer and soil surface.

Only the soil-hydrology module of tRIBS + VEGGIE is used in this study. An adaptation of the quasi-3D framework of the subsurface flow module (Ivanov et al. 2008, 2010) to the mixed formulation (i.e. containing both water content and water potential) of the 1D Richards equation (Celia et al. 1990) permits the computation of dynamics of an unconfined aquifer. Simple averages are used to average unsaturated hydraulic conductivity between layers. A model based on the non-linear Boussinesq equation, accounting explicitly for slope, under the Dupuit–Forchheimer assumptions for lateral flows in the saturated zone and the kinematic approximation for flows in the unsaturated zone are used (Freeze & Cherry 1979). Net lateral flow is applied as a sink/source term in the 1D Richards equation in each vertical layer. The D-inf method is used to compute flow directions for subsurface flow. For the saturated flow, the directions are computed at each computational time step to account for the effective slopes of the hydraulic head surface; the topographic surface is used to compute flow directions for unsaturated flows. This is computed at the beginning of simulations.

The ecohydrological model T&C is a physically-based tool that has been developed to account for the coupled interactions of energy-water-vegetation in a variety of environments and climates where water is the key component (Fatichi 2010; Fatichi et al. 2012a, 2012b, 2015). All essential components of the hydrological cycle are included. A quasi-3D representation of saturated and unsaturated soil water dynamics is used to obtain the subsurface mass exchanges.

The governing equations describing hydraulics are the one-dimensional Richards equation for vertical subsurface flow in variably saturated soils, and the kinematic wave equation for lateral subsurface, overland, and channel flow. The 1D Richards equation is solved in T&C using a finite volume approach with the method of lines, which discretizes the spatial domain and allows reducing the partial differential equation to a system of ordinary differential equations in time (Lee et al. 2004; Fatichi et al. 2012a). Simple averages are used to obtain the unsaturated hydraulic conductivity between layers. The subsurface kinematic wave is solved in a discretized space domain (a lattice of square cells) for each soil layer. An artificial time lag is used in the subsurface routing, since the routing is made at the end of each time step: for a given cell the inflow is the outflow of the upstream cells in the preceding time step. This artificial time lag allows the model to be easily parallelized and run on multiple processors. Incoming lateral flow is then applied as a source term in the 1D Richards equation in each vertical layer, while outgoing lateral flow is computed concurrently to the Richards equation with the kinematic equation and represents a sink term.

Only the soil hydrology module of T&C is used in this study. Flow directions in T&C can be estimated with multiple methods, two methods are used here: D-mult (Quinn et al. 1991; Schwanghart & Kuhn 2010) and D-inf (Tarboton 1997). The flow direction matrix computed using the surface topography obtained with one of the two methods is successively used to route the subsurface water flow. Another approximation is that the flow window width in subsurface exchanges is assigned to be the length of the grid square in the cardinal direction. The same length is also used to compute the actual distance covered by the flow regardless of whether the movement occurs in a diagonal or cardinal direction.

In HYDRUS-3D, the flow domain was discretized into 25,110 nodes (45,866 triangular prisms) for the small domain and 29,106 nodes (52,394 triangular prisms) for the large domain. In both domains, nodes were organized in 18 mesh layers, with distance between the layers increasing from 2 cm at the surface to 8 cm at the base of the domain.

In tRIBS + VEGGIE and T&C, the flow domains were spatially represented by a regular grid of cells in the horizontal plane, with a vertical discretization into 17 layers (‘grid cell layers’) that had the same resolution as layers represented in HYDRUS-3D. In total, 450 grid cells and 7,650 finite volumes, for the small hillslope domain, and 1,250 grid cells and 21,250 finite volumes, for the large domain, were used by both of the models.

In addition to the domain size, two different slope magnitudes were considered: 10 and 30 °. Homogeneous soil with loamy sand texture was adopted, assuming the van Genuchten–Mualem soil hydraulic model (van Genuchten 1980) (see Table 1 for soil hydraulic parameters based on the Carsel and Parrish soil catalog; Carsel & Parrish 1988). Hysteresis in the soil hydraulic functions was not considered. Simulations were initialized with uniformly distributed soil moisture at either θ_i = 0.07 or θ_i = 0.26. With loamy sand texture, these soil water contents correspond to pressure head h = –1 m and h = –0.1 m, respectively.

The applied rainfall input was assumed to represent typical storm events during spring season over the period of March–June in humid, temperate regions, such as the US Pacific Northwest. Rainfall forcing was distributed over the first 9 hours of the simulations, followed by 60 days (or 1,440 h) of drainage. Two total rainfall depths were applied: 50 mm (at a temporally uniform input rate of 5.6 mm h^–1 over 9 hours) and 85 mm (at a temporally uniform input rate of 9.4 mm h^–1 over 9 hours). The storm duration of 9 hours was derived from mean rainfall statistics of 29 years of data from the H. J. Andrews Experimental Forest, Oregon, PRIMET weather station (publicly available at http://andrewsforest.oregonstate.edu/lter/). Input rates typical for a 9-hour-event in the Western Cascades of Oregon were derived from intensity-duration-frequency curves (ODOT 2005). The lower rate corresponds to a return interval of 2 years and the higher rate to a return interval of 50 years.

Combining all permutations led to 16 scenarios: two domain slopes × two rainfall rates × two initial conditions × two domain sizes.

Depending on the model, the formulation of boundary conditions slightly varies. In all models, no-flux boundary condition is applied at the soil bottom and all sides of the domain, except the downslope face of the hillslope (Figure 1) where subsurface flow will leave the hillslope domain. At this downslope face a seepage face boundary condition is implemented, which allows water to leave the flow domain through nodes where local pressure head is simulated as h ≥ 0 m, i.e. where nodes experience saturation.

In HYDRUS-3D, a prescribed flux is specified at the surface for the first 9 hours of the simulation, i.e. during the period of rainfall input. After that, the boundary condition at the surface is also switched to a seepage face boundary condition. This is to accommodate saturation excess overland flow and to permit runoff through the surface during the drainage of the hillslope. Water that leaves the domain through the seepage face at the surface is removed immediately and cannot re-infiltrate further downhill.

In tRIBS + VEGGIE and T&C, the flux boundary condition is specified at the soil surface, which is either positive (the first 9 hours of simulation) or zero. The seepage face is approximated for pre-defined cells located at the downslope face as moisture sinks in nodes where pressure head is larger than zero. The sink strength is taken as flow to the soil medium with zero pressure head component. In both models surface runoff is also immediately removed from the domain; no surface-channel routing and re-infiltration processes are simulated to match HYDRUS-3D treatment.

Furthermore, maps of spatially distributed soil moisture at key time steps (time of peak discharge for the respective simulations) and integrated over time were also produced for comparison.

Overall, tRIBS + VEGGIE results agreed quite well with HYDRUS-3D with NSE > 0.75 in most of the cases, indicating a very similar model response (Figure 7). RMSE was generally small and exhibited dependence on flow magnitude. More significant deviations were discernible for dry initial conditions, particularly for the large domains (Figures 5–7, NSE <0.5).

The comparison with T&C was more variable. For the small domain with an overall slope angle of 10°, neither the partitioning nor the runoff magnitudes were represented particularly well (Figure 3). In the case of dry antecedent moisture conditions (Figure 3, the top two rows of plots), T&C simulated surface runoff, whereas HYDRUS-3D and tRIBS + VEGGIE did not. With wet initial conditions (Figures 3 and 4, the bottom two rows of plots), T&C simulated less subsurface flow but considerably more surface runoff than HYDRUS-3D and tRIBS + VEGGIE, resulting in an overestimation of the total peak flows by a factor of almost 2. Accordingly, the NSE was low or even negative in many cases (Figure 7). The agreement between T&C and HYDRUS-3D was much better in the case of the small hillslope domain with a 30° slope angle, with respect to simulated total runoff as well as the partition between surface and subsurface runoff (Figure 4). For the large domains, T&C in general reproduced HYDRUS-3D results quite well for wet antecedent moisture conditions (ism 26), as also indicated by the high NSE and low RSME (Figures 5 and 6). As in the case with the small domains, the partition between surface and subsurface runoff was quite different for the large domains. This was especially true for dry initial conditions.

Comparing the two flow routing algorithms that were implemented in T&C suggests that differences in timing and magnitude of simulated runoff are most pronounced in the scenarios of the hillslope domains with a low slope angle (Figures 3 and 5). For 30° slopes, both flow routing algorithms produce similar hydrologic responses (Figures 4 and 6). The disagreement between the two routing schemes was most pronounced for very small runoff values.

The dynamics of soil moisture were often very similar, particularly between HYDRUS-3D and tRIBS + VEGGIE (Figures 3–6). The same applied to T&C except for the cases with the small 10° hillslope and wet initial conditions (Figure 3). The spatial variation of depth-averaged soil moisture showed few differences between HYDRUS-3D and tRIBS + VEGGIE on one hand and T&C on the other hand, particularly for the small domains (Figures 3 and 4).

The results for total runoff and soil moisture dynamics (the mean and the coefficient of spatial variation) were often very similar among all of the three models despite their different assumptions and complexity (Figures 3–6). The larger relative differences in runoff occur for dry initial conditions and for the large hillslope (Figure 2) when the total runoff from the hillslope is very small. Since evapotranspiration is neglected in the current simulations, differences in total runoff may occur only as a consequence of a different storage (spatial distribution of soil moisture) of event water within the hillslope domain, which apparently is only the case for the dry conditions. Since differences in total runoff occur only when runoff amounts are small, event runoff coefficients are almost identical in the three models for all the analyzed cases. The results for tRIBS + VEGGIE based on the Dupuit–Forchheimer assumption for the saturated zone dynamics are almost indistinguishable from the HYDRUS-3D solution in all of the cases with wet initial condition. The results start to diverge for the dry initial conditions, when the dynamics of flow through the seepage face likely become sensitive to representation of gradients near the boundary. The simplest model approximation, the kinematic approximation of T&C is able to produce consistent total runoff hydrographs (in terms of shape, timing, and amount) in most of the situations with the exception of the small hillslope with 10° slope for which the hydrograph peaks are overestimated by a factor of two. Despite similar total runoff hydrographs and mean soil moisture dynamics, the partition between subsurface and surface runoff production is often different in T&C, as compared to tRIBS + VEGGIE/HYDRUS-3D simulations, particularly for the domains with 10° slope (Figures 3 and 5). T&C tends to simulate larger rates of surface runoff as a consequence of larger saturated areas induced by more concentrated, topographically governed flow (see map in Figure 8 for an example). The fact that different runoff generation mechanisms lead to similar total runoff hydrographs can be related to the relatively small size of the hillslope, and the relative fraction of topographically convergent areas with respect to the total domain area. Differences in the types of runoff production decrease for the larger domains where the fraction of topographically convergent areas is smaller. This suggests that the simulation of runoff production mechanisms can be related to the adopted spatial resolution and that flow-approximation models can be similar to 3D models in terms of runoff generation mechanisms at spatial scales where runoff generation is commonly simulated.

Overall, for the dry initial condition, the disagreement among the simulated hydrographs becomes larger. However, the flow magnitudes are very small and runoff is mostly produced as seepage flow through the open face of the downstream boundary. Hydrograph shapes are also likely to be affected by the numerical features of the code implementations, e.g. by the type of averaging the unsaturated conductivities (Szymkiewicz & Helmig 2011) and the tolerance magnitudes used in the solution of differential equations. This result highlights that detailed modeling of low flows is problematic even in idealized and homogeneous cases and is likely to be a nearly impossible task in real-world situations that are characterized by soil heterogeneity and a high degree of uncertainty of the boundary conditions. Given the magnitudes of low flows, this is unlikely to represent a significant practical problem in most applications, since larger temporal scales are usually of interest for low flows and, as shown here, the total runoff production in the three models is very similar.

In agreement with theoretical expectations, the kinematic wave approximation for the subsurface flow (T&C) is adequate for steeper and larger hillslopes. A more satisfactory performance is achieved for wet rather than dry initial conditions, as also found by Hilberts et al. (2007). Given the differences in the formulations of T&C and tRIBS + VEGGIE, the poor performance of the former for the 10° hillslope can be attributed to the predefined, topographically controlled flow directions and the lack of accounting for the actual water table dynamics (i.e. what is accounted for by the Dupuit–Forchheimer approximation). Note that both models are based on the 1D Richards equation applied in a quasi-3D framework and thus the observed differences can be attributed to resolving (or not) the hydraulic head gradient in the saturated zone. When soil water is routed according to the local hydraulic head gradients (i.e. the tRIBS + VEGGIE approximation), the one-dimensional solution applied in a quasi-3D fashion is nearly identical to the fully 3D solution. Such a result contrasts the initial expectation of large differences among the three models and raises the question of necessity of incorporating the third dimension in the formulation of subsurface flow dynamics for shallow homogeneous soils. The effect of accounting for the third dimension appears to be more significant for the initially dry cases with very low flows (Figures 5 and 6). However, it is impossible to state with a high level of certainty whether this is due to the sensitivity of gradient representation near the seepage face or due to the effects of vadose zone dynamics, in which the 3D process of suction may play a non-negligible role.

Based on the presented cases, one can argue that concerning longer-term hydrological dynamics, such as water budget or features of spatio-temporal variability of depth-integrated soil moisture (Fatichi et al. 2015), all three models perform quite similarly and are, basically, identical for steeper slopes and larger domains. Therefore kinematic or Dupuit–Forchheimer are valid approximations. At the event scale, when the actual flow time series (or sub-hillslope soil moisture variations) are concerned, the hydrograph shape and peak flow cannot be well captured by the kinematic wave approximation, but Dupuit–Forchheimer is very efficient in approximating 3D flows, as also shown by Hilberts et al. (2007).

The two different subsurface flow direction algorithms implemented in the model T&C demonstrate that for steeper domains, the two methods tend to converge and the differences are minimal. However, the D-mult approximation yields results that are in better agreement with HYDRUS-3D, as compared to the results obtained with the D-inf type of flow routing. The D-mult algorithm produces a higher flow dispersion that better mimics the hydraulic head gradients impacting lateral redistribution of moisture, thereby reducing the concentration of flow and preventing the formation of strictly topography-controlled saturated areas. Note that these inferences are applicable when the topographic surface is used to represent hydraulic head gradients. In cases when flow direction methods are used to approximate the 2D hydraulic head surface to derive the distribution of gradients (i.e. the Dupuit–Forchheimer approximation), the D-inf algorithm is the most suitable method and produces results that are more consistent with HYDRUS-3D output, when compared to other flow directions schemes (not shown in this study).

This study tested if flow-approximation models that implement a quasi-3D framework based on the 1D vertical Richards equation can produce adequate estimates of the dynamics of both unsaturated and saturated zones in hillslope-scale domains. The flow-approximation models tRIBS + VEGGIE (using the Dupuit–Forchheimer assumption) and T&C (using the kinematic approximation) were compared to results produced by the model HYDRUS-3D that solves the complete 3D formulation of the Richards equation. Event-based simulations at the scale of a hillslope exhibiting real-world complexity in the organization of surface topography were carried out to isolate the water flow problem from other land-surface hydrological dynamics. The effects of slope angle, domain size, rainfall magnitude, and initial moisture conditions on runoff generation and simulated soil moisture dynamics were examined.

Results show that total runoff and soil moisture spatial variability for most of the cases are very similar between the two flow-approximation model formulations and HYDRUS-3D. The results from the flow-approximation models compare very well with HYDRUS-3D results for steeper and larger hillslopes. In line with the results of Hilberts et al. (2007), a more satisfactory performance is achieved for wet rather than dry initial conditions. However, in the dry cases the total runoff is very small and reasons for a less favorable performance can stem both from the details of numerical implementation (all models), particular sensitivities of flow approximation near the seepage face of the downstream boundary (tRIBS + VEGGIE and T&C), and an exaggeration of hydraulic gradients (specific to T&C that demonstrated larger flows).

At the scale of a storm event, when the actual flow time series or details of sub-hillslope soil moisture variations are concerned, the Dupuit–Forchheimer approximation in combination with the D-inf flow direction algorithm is very efficient in approximating 3D flows and can capture both the hydrograph shape and peak flow reasonably well, as compared to the explicit 3D model. The performance of the kinematic wave approximation with the two types of flow direction definition is similar for large events, steep slopes, and high initial water storage. For other cases, using topographic gradients instead of gradients of hydraulic head is not justified when hydrologic response is required to be resolved at fine temporal resolution (i.e. sub-event).

For longer-term hydrologic dynamics, such as water budget partition or features of spatio-temporal variability of depth-integrated soil moisture, all three models perform quite similarly and thus flow-approximation models are valid tools when hydrologic investigation is concerned at larger than storm event time scales. Note, however, that this pilot study did not include cases with multiple storms, soil heterogeneities and other processes of the land-surface hydrologic cycle and therefore such a conjecture should be interpreted with caution. Soil heterogeneity, in particular, could add a new dimension to the problem with the possibility to channel subsurface flow preferentially through permeable layers or create pools of stagnant water or saturation in regions or layers with poorly permeable soils. Soil heterogeneity can also increase the importance of numerical choices such as the averaging of conductivity values between soil layers or numerical convergence criteria. Preliminary analyses with heterogeneous soils (not shown) are suggesting that the similarity in total runoff and soil moisture dynamics among models remain valid also in the occurrence of heterogeneous soils but additional testing is required to fully explore the role of spatial variability in soil properties.

The initial expectation of considerable differences between the two-flow approximation models and the explicit three-dimensional model has not been fully confirmed by this study. The presented analysis was limited to homogeneous soil conditions and uniform soil depth (1 m). We hypothesize, however, that the third dimension in the flow formulation can become relevant for (i) flatter domains and deeper soils, as compared to the domains used in this study, (ii) long, sloped thick soils with large saturated reservoir, and (iii) topographically more complex domains with heterogeneous soil conditions, non-uniform soil thicknesses and permeable bedrock. These hypotheses can be tested in future studies.

V. Y. Ivanov was supported by NSF Grant EAR 1151443 and the Visiting Faculty grant at the Institute of Environmental Engineering, ETH Zürich. S. Fatichi would like to thank Enrica Caporali and the University of Firenze (Italy) for support at the beginning of this research. L. Hopp would like to thank J. J. McDonnell and M. Goeckede for supporting this work.