# Vadose Zone Journal

- Copyright © by the Soil Science Society of America, Inc.

## Abstract

Surface runoff is commonly described in numerical models using either the diffusion wave or kinematic wave equations, which assume that surface runoff occurs as sheet flow with a uniform depth and velocity across the slope. In reality, overland water flow and transport processes are rarely uniform. Local soil topography, vegetation, and spatial soil heterogeneity control directions and magnitudes of water fluxes. These spatially varying surface characteristics can generate deviations from sheet flow such as physical nonequilibrium flow and transport processes that occur only on a limited fraction of the soil surface. In this study, we first adapted the HYDRUS-1D model to solve the diffusion wave equation for overland flow at the soil surface. The numerical results obtained by the new model produced an excellent agreement with an analytical solution for the kinematic wave equation. Additional model tests further demonstrated the applicability of the adapted model to simulate the transport and fate of many different solutes (non-adsorbing tracers, nutrients, pesticides, and microbes) that undergo equilibrium and/or kinetic sorption and desorption and first- or zero-order reactions. HYDRUS-1D includes a hierarchical series of models of increasing complexity to account for both uniform and physical nonequilibrium flow and transport, e.g., dual-porosity and dual-permeability models, up to a dual-permeability model with immobile water. This same conceptualization was adapted to simulate physical nonequilibrium overland flow and transport at the soil surface. The developed model improves our ability to describe nonequilibrium overland flow and transport processes and our understanding of factors that cause this behavior.

- ADE, advection–dispersion equation
- APR, active–passive regions
- APR-H, combined active–passive regions and horizontal mobile–immobile regions
- APR-V, combined active–passive regions and vertical mobile–immobile regions
- BTC, breakthrough curve
- GUI, graphical user interface
- HMIM, horizontal mobile–immobile regions
- UFT, uniform flow and transport
- VMIM, vertical mobile–immobile regions

Contaminants at the soil surface (e.g., pesticides, heavy metals, and pathogenic microbes) can be rapidly transported to streams or locations of surface water storage by overland flow. This has been reported to be the primary transport route for contaminant dissemination in agricultural settings (Carpenter et al.,1998; USGS, 1999; Tyrrel and Quinton, 2003). An understanding of and ability to predict processes that influence the transport and fate of contaminants in runoff water is therefore needed to assess and mitigate the risks of contamination of surface water supplies on human health (Furman, 2008; Vereecken et al., 2016).

The diffusion wave equation describes surface runoff as sheet flow with a uniform depth and velocity across the slope. In reality, overland water flow and pollutant transport are rarely uniform. Local soil microtopography, surface roughness, vegetation, and spatial soil heterogeneity vary across distances of centimeters to meters, and they control the directions and magnitudes of water fluxes and the concentrations of pollutants (Zhang and Cundy, 1989). There is increasing evidence that overland flow and transport processes often cannot be described using classical overland flow and transport models that assume uniform flow and transport (e.g., Smith et al., 2011; Cea et al., 2014; Bradford et al., 2015). At the local scale, spatially varying roughness, vegetation, and microtopography influence the distribution of shear stress and create hydrologically active and passive or relatively immobile flow regions. Surface active regions for overland flow can be formed in the field by rills or connected networks of microdepressions, which route overland water and pollutant fluxes on the soil surface where the velocity can be two to seven times higher than the average flow (Dunkerley, 2003; Chen et al., 2013). Great risk of contamination to water resources occurs from hydraulically active regions because of the presence of large water and contaminant fluxes and reduced time for infiltration and contaminant decay, transformation, and/or interphase mass transfer (e.g., interactions with soil particles). Surface passive or relatively immobile regions are formed by regions of depression and obstruction storage that retain water and pollutants and inhibit overland flow or cause shallow, slow-moving flow and exchange with surface active zones. These spatially varying surface characteristics generate nonequilibrium overland flow and transport processes.

Physically based, spatially distributed models have the potential to be an efficient tool to examine and optimize the removal of contaminants from overland flow through land-use changes and best management practices. Many publicly available overland pollutant transport models have been developed (Beven and Kirkby, 1979; Haith and Shoemaker, 1987; Bathurst and Connell, 1992; Flanagan et al., 1995; Arnold et al., 1998). Lumped parameter conceptual models have commonly been used for contaminant transport in overland flow but may be inaccurate because they ignore processes acting on water and chemicals at the soil surface (e.g., Soil and Water Assessment Tool [SWAT; Sadeghi and Arnold, 2002], Hydrologic Simulation Program Fortran [HSPF; Donigian et al., 1995], Integrated Nitrogen Catchment model [INCA; Whitehead et al., 1998; Wade et al., 2002], and COLI [Walker et al., 1990]). Other mechanistic models consider local-scale parameters (e.g., Parallel Flow [ParFlow; Kollet and Maxwell, 2006], CATchment Hydrology [CATHY; Camporese et al., 2010], HydroGeoSphere [HGS; Therrien et al., 2012], and OpenGeoSys [OGS; Kolditz et al., 2012]). For example, models that use the actual varying microtopography based on digital elevation maps have become popular (Gómez and Nearing, 2005; Chen et al., 2013; Zhao and Wu, 2015), but they are computationally demanding and extra effort is needed for collecting the required data. In an attempt to overcome these limitations, some mechanistic models neglect local-scale variations in parameters by considering simplified smooth surfaces with a global roughness coefficient and a constant average slope; e.g., KINematic runoff and EROSion (KINEROS2; Miller et al., 2007; Goodrich et al., 2012; Kennedy et al., 2013) and Precipitation–Runoff Modeling System (PRMS) (Leavesley et al., 1983). These models may provide an excellent approximation of overland flow when calibrated to real data but may not be adequate to describe overland transport that is affected by local-scale parameters. Hence, accurate, physically based modeling of nonequilibrium water flow and solute transport still remains a challenge in the field of surface hydrology, and there is currently no widely accepted physically based model available to simulate overland transport processes of all kinds of reactive chemicals.

Recent studies have demonstrated that equations for overland water flow (the diffusion wave equation) and transport (the advective–dispersion equation, ADE) are equivalent to those for the subsurface (the Richards and ADE equations) when using different functional forms for water content, water capacity, and hydraulic conductivity (Weill et al., 2009; Bittelli et al., 2010). The HYDRUS-1D software (Šimůnek et al., 2016), a popular numerical computer code for solving water flow and reactive transport in variably saturated porous media, can thus be modified to simulate overland flow rather than subsurface flow. Existing subsurface flow and transport models in HYDRUS-1D include a hierarchical series of physical nonequilibrium models of increasing complexity, from dual-porosity and dual-permeability models up to a dual-permeability model with immobile water (Šimůnek and van Genuchten, 2008). All these nonequilibrium models were derived from the Richards and convection–dispersion equations.

The objective of this study was to adapt the existing subsurface version of HYDRUS-1D to simulate similar uniform or physical nonequilibrium flow and reactive solute transport processes during runoff at the soil surface. We first present a detailed description of the implemented equations for overland flow and transport, followed by some examples of model results with a special focus on modeling sorption–desorption processes and physical nonequilibrium flow and transport. A limited comparison of model results and experimental and numerical data that exhibit deviations from sheet flow is also provided. The developed models improve our ability to quantify natural flow and transport processes at the soil surface, and the simulation results enhance our understanding of the factors that cause deviations from sheet flow.

## Numerical Models

A variety of physical and chemical equilibrium and nonequilibrium models is available in HYDRUS-1D to describe water flow and solute transport processes in the subsurface (Šimůnek and van Genuchten, 2008). These models form a hierarchical system of models of increasing complexity to account for both physical equilibrium and nonequilibrium flow, including (i) a uniform flow model, (ii) a dual-porosity (mobile–immobile water) model, (iii) a dual-permeability model, and (iv) a dual-permeability model with immobile water. These conceptual models for subsurface flow were adapted in this study to account for both equilibrium and nonequilibrium overland flow and reactive solute transport, resulting in (i) a uniform flow and transport (UFT) model, (ii) a horizontal mobile–immobile regions (HMIM) model, (iii) a vertical mobile–immobile regions (VMIM) model, (iv) an active–passive regions (APR) model, (v) a combined APR and HMIM (APR-H) model, and (vi) a combined APR and VMIM (APR-V) model. While the mobile–immobile models for overland flow correspond to the dual-porosity models for subsurface flow, the active–passive region models correspond to the dual-permeability subsurface models. Both mobile–immobile and active–passive regions models have two additional subsets with either horizontal or vertical subregions. While the focus in this study was on physical nonequilibrium models that describe conditions when flow and transport occur only on a limited fraction of the soil surface that bypasses regions with little or no flow, the implemented models can also account for chemical nonequilibrium when sorbed concentrations are not in equilibrium with liquid concentrations. The schematics of newly developed overland flow models discussed below are shown in Fig. 1. Below we also describe the governing equations for each of these models. Table 1 provides a comparison of the different definitions of the variables for subsurface and surface models.

### Uniform Flow and Transport Model

Overland flow and solute transport are commonly described using diffusion wave and advection–dispersion equations, respectively. This approach lumps irregular land surface characteristics into effective parameters that are used to simulate uniform sheet flow and transport (Fig. 1a). The diffusion wave equation for overland flow may be written in a similar form to the Richards equation:where *x* is a space coordinate in the direction of flow [L, where L denotes units of length], *t* is time [T, where T denotes units of time], *h* is the surface water depth [L], *q* is the source–sink term [L T^{−1}] accounting for precipitation, evaporation, and infiltration, *z* is the land surface elevation [L], *S* is the mean local slope (dimensionless), *n* is a Manning’s roughness coefficient for overland flow (dimensionless), and *k* is a unit conversion factor [L^{1/3} T^{−1}]. Detailed derivation of this equation can be found, for example, in Weill et al. (2009), Hromadka and Lai (1985), Panday and Huyakorn (2004), and Šimůnek (2015). The *n* parameter is either dimensionless when the units conversion factor *k* is used or has units of [T L^{−1/3}] when it is not. The parameter *k* can also be used to convert Eq. [1] between SI and US customary units and can be left out when consistent units are used. However, it is a standard practice to use *k* = 1 m^{1/3} s^{−1} for SI units and *k* = 1.49 ft^{1/3} s^{−1} for US customary units.

When rainfall (*R* [L T^{−1}]) (or irrigation), evaporation (*E* [L T^{−1}]), and infiltration (*I* [L T^{−1}]) rates are the only source–sink terms, then *q* = *R* – *E* − *I*. Infiltration is usually determined using various empirical, semi-empirical, or physical models. In this work, the infiltration rate *I* is described using Horton’s equation (Horton, 1939):where *f*_{0} is the initial infiltration rate [L T^{−1}], *f*_{c} is the final equilibrium infiltration rate [L T^{−1}], and *k*_{i} is a constant representing the rate of decrease in infiltration [T^{−1}]. According to Eq. [2], infiltration starts at a rate *f*_{0} and then decreases exponentially with time until it reaches an equilibrium infiltration rate *f*_{c}. The use of Horton’s equation is only the first attempt in our model to consider the infiltration process, since the focus of this study is on overland flow and transport processes. This infiltration model will be replaced in the future by other empirical and/or process-based infiltration models.

Solute transport in overland flow is usually described using the ADE of the formwhere *c* is the solute concentration in the aqueous phase [M L^{−3}, where M denotes units of mass], *c*_{r} is the concentration in rainfall water [M L^{−3}], *s* is the sorbed solute concentration at the soil surface area [M L^{−2}], *D* is the effective dispersion coefficient accounting for both molecular diffusion and hydrodynamic dispersion [L^{2} T^{−1}], ϕ is a sink–source term that accounts for various zero- and first-order or other reactions [M L^{−3} T^{−1}], and *Q* is the runoff flow rate [L^{2} T^{−1}]. The parameter *Q* is given aswhere *U* is a depth-averaged velocity [L T^{−1}] calculated using the Manning–Strickler uniform flow formula (Hromadka and Lai, 1985). It should be mentioned that the subsurface pore-water velocity and *U* have the same units, but the subsurface Darcy’s velocity and *Q* do not have the same units. The effect of diffusion on the dispersion coefficient can often be ignored and in this case, *D* can be defined as the product of the dispersivity (λ [L]) and *U*.

A variety of sorption processes can be considered in HYDRUS-1D, including chemical equilibrium described using both linear and nonlinear adsorption isotherms and chemical nonequilibrium described using first-order kinetic models. These same sorption models may also be used with overland transport. However, the units for *s* are different for subsurface [M M^{−1}] and overland [M L^{−2}] transport. Furthermore, the soil bulk density is not required in the overland transport equation, contrary to the subsurface transport equation. In the examples given below, the value of *s* in the equilibrium model is given aswhere *K*_{d} is the distribution coefficient [L]. The one-site kinetic sorption model is written aswhere ω is the first-order rate coefficient [T^{−1}] representing kinetic sorption. The combined equilibrium and kinetic sorption model is written aswhere *f*_{eq} (dimensionless) is the fraction of equilibrium sorption sites, and the subscripts *eq* and *k* on *s* denote equilibrium and kinetic sorption sites given by expressions similar to Eq. [5a] and [5b], respectively. The two-site kinetic sorption model is given aswhere *k*_{att1} [T^{−1}] and *k*_{att2} [T^{−1}] are kinetic sorption rate coefficients, *k*_{det1} [T^{−1}] and *k*_{det2} [T^{−1}] are kinetic desorption rate coefficients, and the subscripts 1 and 2 on the parameters denote the kinetic sorption sites 1 and 2, respectively. Note that the latter model is commonly also used to describe attachment and detachment processes when modeling transport of particular substances, such as viruses, colloids, or pathogens. It should be mentioned that Eq. [5b] can be recast in terms in *k*_{att1} = ω*K*_{d}/*h* and *k*_{det1} = ω.

The upper boundary condition (BC) for water can be either a prescribed water head or flux:where *h*_{0} [L] and *q*_{0} [L^{2} T^{−1}] are the water head and flux at the upper boundary, respectively. The water depth gradient at the lower BC is assumed to be zero:

Two types of BCs can be applied for transport at the upper or lower boundaries. The first-type boundary conditions prescribe the concentration at the boundary:whereas third-type (Cauchy type) boundary conditions may be used to prescribe the concentration flux at the boundary:where *c*_{0} is the concentration of the incoming fluid [M L^{−3}].

### Horizontal Mobile–Immobile Regions Model

A schematic of the HMIM model is shown in Fig. 1b. This model assumes that the soil surface is horizontally divided, parallel to the direction of water flow, into regions with mobile and immobile water. Water and solutes in the immobile region may be stored, retained, and exchanged with the mobile domain. Similar to Eq. [1], the movement of water in the mobile region and moisture dynamics in the immobile region are given as:where *E*_{w}^{MIM} is the water transfer rate between mobile and immobile regions [L T^{−1}], *w*_{m} is the ratio of the width of the mobile region to the total width of the soil surface, and subscripts *m* and *im* denote parameters associated with the mobile and immobile regions, respectively. The parameter *E*_{w}^{MIM} is assumed to be proportional to the difference in water depths between the two regions:where α [T^{−1}] is a first-order mass transfer coefficient (Šimůnek et al., 2003). The average surface water depth of the entire domain is

The governing solute transport equations for the mobile and immobile regions arewhere *E*_{s}^{MIM} is the solute transfer rate between mobile and immobile regions [M L^{−2} T^{−1}]. The parameter *E*_{s}^{MIM} is given aswhere ω_{m} is the solute mass transfer coefficient [T^{−1}] and *c** is a concentration that is equal to *c*_{m} for *E*_{w}^{MIM} > 0 and *c*_{im} for *E*_{w}^{MIM} < 0. Note that solute exchange between the two liquid regions is modeled in Eq. [14] as the sum of an apparent first-order diffusion process and advective transport.

It should be mentioned that Eq. [10] and [13] are written in terms of local-scale mass balances in mobile and immobile regions. To formulate them in terms of the total region, the mass balance equations for mobile and immobile regions need to be multiplied by *w*_{m} and (1 − *w*_{m}), respectively.

### Vertical Mobile–Immobile Regions Model

A schematic of the VMIM model is presented in Fig. 1c. This model assumes that the surface water depth is vertically divided, parallel to the mean surface slope, into an immobile region adjacent to the soil surface and an overlying mobile domain. This conceptual picture is consistent with an immobile depression and/or obstruction storage zone that needs to fill before the initiation of overland flow (Panday and Huyakorn, 2004; Guber et al., 2009). An immobile region is supposed to account for the effects of small surface roughness and/or vegetation that prevents an immediate surface runoff once ponding is reached. It represents a small water layer that needs to be formed before surface runoff can occur. Water flow is described as

In contrast to the HMIM model, the parameter *E*_{w}^{MIM} is now given aswhere *h*_{im}^{max} is the maximum water depth for the immobile region [L], *h*_{T} is the total water depth [L] that is equal to *h*_{m} + *h*_{im}, *Q*_{up} is flow from upstream [L^{2} T^{−1}], and Δ*L* is the length of the surface [L] across which inflow from upstream flows into the immobile region. The Heaviside function (*H*_{o}) in Eq. [16] is equal to 1 when *h*_{im}^{max} > *h*_{T} and is 0 when *h*_{im}^{max} ≤ *h*_{T}. Consequently, water flow in the mobile domain occurs only when the immobile region is filled (*h*_{im}^{max} ≤ *h*_{T}) at a particular location.

Solute transport equations for the mobile and immobile regions are given as

Note that the transport equation for the mobile region does not include the sorption term because this region is not in direct contact with the soil. The parameter *E*_{s}^{MIM} is given as

### Active–Passive Regions Model

A schematic of the APR model is given in Fig. 1d. This model assumes that the soil surface is divided, parallel to the direction of water flow, into hydraulically active (fast flow) and passive (slow flow) domains. This model allows water and solute transport in both active and passive regions and exchange between these regions. Surface water flow in each region is, in analogy of the dual-permeability subsurface flow model of Gerke and van Genuchten (1993), described using separate diffusion wave equations:where *E*_{w}^{APR} is the water transfer rate between active and passive regions [L T^{−1}], *w*_{A} is the ratio of the width of the surface active region and the total surface width (dimensionless), and the subscripts 1 and 2 refer to hydrologically active and passive regions, respectively. Values of *E*_{w}^{APR} and *h* are quantified in a similar manner to Eq. [11] and [12]:

Solute transport is described using separate ADEs for each region:where *E*_{w}^{APR} is the solute transfer rate between active and passive regions [M L^{−2} T^{−1}]. The *E*_{s}^{APR} term is quantified in a manner similar to Eq. [14]:where ω_{12} is the solute mass transfer coefficient for transfer between active and passive regions [T^{−1}].

The total flux is obtained as

and the average flux concentration is

It should be mentioned that Eq. [19] and [22] are written in terms of local-scale mass balances in Regions 1 and 2. To formulate them in terms of the entire region, the mass balance equations for Regions 1 and 2 need to be multiplied by *w*_{A} and (1 − *w*_{A}), respectively.

### Combined Active–Passive Regions and Horizontal Mobile–Immobile Models

A schematic for the combined APR and HMIM (APR-H) models is shown in Fig. 1e. This model assumes that the soil surface is horizontally divided, in the direction parallel to water flow, into four regions: (i) an active (fast flow) region; (ii) a passive (slow flow) region; (iii) an immobile region that is in contact with and exchanges mass (water and/or solute) horizontally with the active region; and (iv) an immobile region that is in contact with and exchanges mass horizontally with the passive region.

Several different surface water depths are considered in the combined APR-H model. Values of *h*_{1} and *h*_{2} are divided into mobile (*h*_{1m} and *h*_{2m}) and immobile (*h*_{1im} and *h*_{2im}) regions using Eq. [12] such that the average water depth for the entire domain (*h*) is given as

The movement of water in active, passive, and immobile zones is given aswhere *E*_{1w}^{MIM} and *E*_{2w}^{MIM} denote the water transfer rate between mobile and immobile domains in active and passive regions [L T^{−1}], respectively. The parameters *E*_{1w}^{MIM} and *E*_{2w}^{MIM} are obtained in an analogous fashion to Eq. [11] using separate values of α, *h*_{m}, and *h*_{im} for active (α_{1}, *h*_{1m}, and *h*_{1im}) and passive (α_{2}, *h*_{2m}, and *h*_{2im}) regions, whereas *E*_{w}^{APR} is determined using Eq. [20].

The following governing equations describe solute transport in the surface active, passive, and immobile regions:where *E*_{1s}^{MIM} and *E*_{2s}^{MIM} denote the solute transfer rate between mobile and immobile domains in active and passive regions [M L^{−2} T^{−1}], respectively. The parameters *E*_{1s}^{MIM} and *E*_{2s}^{MIM} are obtained in an analogous fashion to Eq. [14] using separate values of ω_{m}, *w*_{m}, *h*_{im}, *c*_{m}, *c*_{im}, and *E*_{w}^{MIM} for active (ω_{1}, *w*_{1m}, *h*_{1im}, *c*_{1m}, *c*_{1im}, and *E*_{1w}^{MIM}) and passive (ω_{2}, *w*_{2m}, *h*_{2im}, *c*_{2m}, *c*_{2im}, and *E*_{2w}^{MIM}) regions, whereas *E*_{s}^{APR} is determined using Eq. [23].

### Combined Active–Passive Regions and Horizontal Mobile–Immobile Models

A schematic for the combined APR and VMIM (APR-V) models is shown in Fig. 1f. Similar to the APR model, the combined model assumes that the soil surface is horizontally divided, in the direction parallel to water flow, into active (fast flow) and passive (slow flow) regions. Each of these regions is further divided vertically into an immobile region adjacent to the soil surface and an overlying mobile domain.

The values of *h*_{1} and *h*_{2} are again divided into mobile (*h*_{1m} and *h*_{2m}) and immobile (*h*_{1im} and *h*_{2im}) regions for the combined APR-V model. In this case, the average water depth for the entire domain (*h*) is given as

The movement of water in active, passive, and immobile zones is given as

Values of *E*_{1w}^{MIM} and *E*_{2w}^{MIM} are given in this model similar to Eq. [16] using separate values of *R*, *E*, *h*_{T}, and *h*_{im}^{max} for active (*h*_{1T} and *h*_{1im}^{max}) and passive (*h*_{2T} and *h*_{2im}^{max}) regions.

Solute transport in the surface active, passive, and immobile regions are described using the following model:

Parameters *E*_{1s}^{MIM}, *E*_{2s}^{MIM}, and *E*_{s}^{APR} are determined in a similar manner to Eq. [18], [18], and [23], respectively. Note that the APR-V model assumes that there are no interactions between the two immobile regions under the active and passive flow regions. The additional assumption in the current implementation of the APR-V model is that the maximum water depths *h*_{im}^{max} in the two immobile regions are the same.

## Numerical Implementation

The HYDRUS software uses the Galerkin-type linear finite element method (FEM) for spatial discretization of the governing partial differential equations and the finite difference method to approximate temporal derivatives. The FEM is commonly used in hydrological models, such as in HydroGeoSphere and OpenGeoSys. An important advantage of the FEM over the finite difference method to solve the one-dimensional diffusion wave equation is the ease of increasing the mesh resolution, so that shock fronts can be simulated without any numerical oscillations. Because the governing overland flow equations were formulated using the subsurface flow equations with different definitions of the water content, hydraulic capacity, and conductivity coefficients (see Table 1), the original fully implicit finite difference scheme with Picard linearization was used to solve the overland flow equation as well, while a Crank–Nicholson finite difference scheme was used to solve the ADEs. Details about the numerical solutions of the corresponding subsurface flow and transport models are provided in the HYDRUS-1D manual (Šimůnek et al., 2016). The same graphical user interface (GUI) in HYDRUS-1D is used to select and execute subsurface and overland flow and transport models.

The numerical implementation of the overland flow equation was verified by comparing the numerical results with the analytical solution of the kinematic wave equation (Woolhiser and Liggett, 1967; Singh, 1996). There was excellent agreement between the numerical and analytical solutions at different times and steady-state conditions. A mass balance was also calculated for all of the overland flow and transport models to verify the correctness of the implementation of more complex models, for which analytical solutions either do not exist or cannot be derived. Results demonstrated that the numerical implementation conserved mass.

## Applications

We developed a large number of numerical examples of overland flow and solute transport using the various model formulations. The flow domain is 100 m long and discretized into 101 finite elements. The slope is 1%, and unless otherwise noted, the surface roughness *n* is 0.01. The dispersivity (λ) is set equal to 10 m. A summary of parameter values used in the following numerical examples is provided in Table 2.

### Uniform Flow and Transport Model

Figure 2 presents an example of uniform overland flow over an impervious (Fig. 2a) and permeable (Fig. 2b) soil surface. In this example, inflow was induced by specifying the pressure head equal to 1 cm at the top boundary for 6 min and then simulating overland flow for a total of 20 min. Water depths are shown as a function of time at different locations (0, 20, 40, 60, and 80 m) from the top of the slope. The inlet water pulse moves along the impervious soil surface in Fig. 2a and reaches the bottom of the profile after about 200 s. A maximum water depth of 1 cm that is consistent with the inlet boundary condition is reached and then gradually reduced when the inflow is stopped. In contrast to Fig. 2a, the water depths are not uniform across the permeable soil surface in Fig. 2b due to infiltration. Infiltration reduces the water mass that reaches the bottom boundary and delays its arrival. Mass balance calculations indicated that the sum of cumulative discharge, infiltration, and water at the soil surface equaled the cumulative inflow. This indicates that mass was conserved during infiltration.

The transport of reactive contaminants requires considering advective–dispersive transport with a range of biogeochemical processes such as adsorption, desorption, volatilization (i.e., mass transfer between different phases), degradation, precipitation, dissolution, etc. To demonstrate the applicability of our reactive transport code, the numerical examples of solute transport with different sorption processes are presented in Fig. 3 for the same water flow conditions as in Fig. 2a. Solute transport was induced using a third-type boundary condition at the inlet with a constant unit concentration for 6 min. Solute concentrations are shown as a function of time at different locations along the slope. Figure 3a shows the simulated transport of a conservative tracer. Similar to water flow (Fig. 2a), the tracer reaches the bottom outlet after around 200 s. The tracer was never eluted with solute-free water. Consequently, the solute concentration remains constant after the water and tracer inflow stopped due to a decreased volume of runoff water and solute input. Figure 3b shows the transport of solute that undergoes linear equilibrium sorption. In this case, the retardation coefficient equals 1 + *K*_{D}/*h*, where *K*_{D} was set equal to 1 cm. In comparison to the conservative tracer (Fig. 3a), the transport of the solute undergoing linear equilibrium sorption is delayed. The retardation factor is equal to 2 at the maximum water depth of 1 cm, but drastically increases as *h* goes to zero. Consequently, the delay in solute transport (Fig. 3b) becomes more pronounced during the receding limb of the hydrograph (Fig. 2a). Figure 3c shows the transport of a solute that undergoes kinetic sorption (sorption rate of 0.01 s^{−1}) and desorption (desorption rate of 0.001 s^{−1}). The solute concentration continuously decreases as it is transported over the soil surface during inflow due to kinetic sorption. The solute concentration rapidly increases after inflow ceases due to a decrease in the volume of surface water and continued solute desorption from the land surface. Eventually, the solute concentration reaches an equilibrium concentration level after inflow ceases. The equilibrium aqueous-phase concentration is *k*_{d}/*h*_{min}*k*_{a}, where *h*_{min} [L] is a minimum water depth.

### Horizontal Mobile–Immobile Regions Model

The HMIM model has two additional input parameters for overland flow compared with the UFT model: (i) the water mass transfer coefficient, α; and (ii) the ratio of the width of the mobile domain to the total width of the surface domain, *w*_{m}. The influence of these two parameters on the outflow rate at the bottom boundary as a function of time is demonstrated in Fig. 4. Similar to Fig. 2a, the inlet water depth was equal to 1 cm for 6 min and the soil surface was impervious. A constant value of *w*_{m} = 0.5 and different values of α = 0, 0.001, 0.005, and 0.01 s^{−1} are considered in Fig. 4a. No exchange of water occurs between the mobile and immobile domains when α = 0, and the resulting outflow rate is equal to the uniform flow model. Conversely, a fraction of inflow water is transferred to and from mobile and immobile regions when α > 0. This process slows down overland flow and delays the arrival of the water front relative to the uniform model. It also produces prolonged tailing in the outflow rate compared with the uniform model. Increasing α produces a faster equilibration of water depths between the mobile and immobile zones and causes a greater delay in outflow. It should be mentioned that the relative importance of α on outflow will increase with smaller values of *w*_{m}. A constant value of α = 0.001 s^{−1} and different values of *w*_{m} = 0.3, 0.5, and 0.7 were considered in Fig. 4b. Increasing *w*_{m} produces a slightly earlier outflow arrival time, a higher outflow rate, and less outflow tailing because of the larger mobile domain area. The relative importance of *w*_{m} increases with smaller values of α.

In contrast to the uniform transport model, the HMIM model also depends on the solute mass transfer coefficient (ω_{m}). Figure 5 presents simulated breakthrough curves (BTCs) at the bottom boundary for a 30-min conservative tracer pulse when using the HMIM model with *w*_{m} = 0.5 and ω_{m} = 0, 0.001, 0.005, and 0.01 s^{−1}. Steady-state overland flow conditions were considered in these simulations to isolate the influence of ω_{m} from water exchange (e.g., α). The BTC is equivalent to the uniform transport model when ω_{m} = 0. Increasing ω_{m} produces greater amounts of diffusive solute exchange and a faster equilibration between mobile and immobile domains. Consequently, increasing ω_{m} produces a greater delay in the initial tracer breakthrough but also less long-term concentration tailing for the considered parameters. It should be mentioned that the relative importance of ω_{m} on BTCs will increase with smaller values of *w*_{m}.

### Vertical Mobile–Immobile Regions Model

The VMIM model allows the consideration of depression storage by assigning different values of the maximum water depth for the immobile region (*h*_{im}^{max}). Figure 6a shows plots of the water depth as a function of time at the bottom boundary when using the VMIM model with values of *h*_{im}^{max} = 0, 0.3, and 0.5 cm. In this case, rainfall was uniformly applied over the soil surface for 20 min at a rate of 0.4 cm min^{−1} and infiltration was also considered (Table 2). Figure 6a demonstrates that the water level rises from the land surface up to *h*_{im}^{max} and that overland flow is initiated when *h*_{T} > *h*_{im}^{max}. Overland flow occurs immediately when *h*_{im}^{max} = 0 cm, otherwise overland flow is delayed until the immobile, depression storage, zone is filled. As we expected, higher values of *h*_{im}^{max} resulted in a greater delay in overland flow. The model also simulated the dynamics of overland flow and infiltration for 40 min after rainfall ceased. In particular, overland flow continued when *h*_{T} > *h*_{im}^{max}, whereas infiltration occurred as long as *h*_{T} > 0.

### Active–Passive Regions Model

Overland flow and transport occurs in two parallel surface regions in the APR model. Water exchange between the active and passive regions is determined by the mass transfer coefficient α_{12}. Separate Manning’s roughness coefficients are used in the surface active (*n*_{1}) and passive (*n*_{2}) regions to obtain different velocities. Figure 7a presents simulated active, passive, and total outflow rates as a function of time when using the APR model with *w*_{A} = 0.5, *n*_{1} = 0.01, *n*_{2} = 0.05, and α_{12} = 0, 0.001, and 0.01 s^{−1}. Overland flow across an impervious soil surface was initiated by setting the inlet boundary to a water depth of 1 cm for 20 min. Water moves independently in active and passive domains when α_{12} = 0 (red lines). In this case, two outflow peaks are observed, with outflow starting after about 200 and 1000 s in the active and passive domains, respectively, due to differences in the roughness coefficient. An increase in α_{12} produces faster water exchange and equilibration that cause the outflow peaks for the active and passive regions to decrease and increase, respectively. The APR model approaches the behavior of the uniform flow model at the average velocity when α_{12} is very high. It should be mentioned that an increase in *n*_{2} decreases the water velocity in the passive region, and the APR model, therefore, approaches the HMIM model when *n*_{2} is very large. Figure 7a illustrates that classic nonequilibrium flow behavior (an early arrival, multiple peaks, and long-term tailing in the outflow rate) can be obtained with the APR model.

Figure 8a presents simulated active (dash-dotted), passive (dashed), and total (solid) BTCs at the bottom boundary for a 40-min conservative tracer pulse when using the APR model with *w*_{A} = 0.5, *n*_{1} = 0.01, *n*_{2} = 0.05, and ω_{12} = 0, 0.001, and 0.01 s^{−1}. Steady-state overland flow conditions across an impervious soil surface (an inlet water depth of 1 cm) were considered in these simulations to isolate the influence of ω_{12} from water exchange (e.g., α_{12}). When the exchange coefficient is equal to zero (red lines), solute moves independently through each of the two surface domains. In this case, the total BTC is the weighted superposition of BTCs from the active and passive regions. Increasing ω_{12} causes greater amounts of diffusive solute exchange and a faster equilibration between the active and passive domains. This produces a decrease and an increase in effluent concentrations in the active and passive regions, respectively. When ω_{12} is very high, the BTCs from the two domains converge and resemble the results of the uniform transport model at the average water velocity (purple line in Fig. 4a).

### Combined Active–Passive Regions and Horizontal Mobile–Immobile Model

The simulations associated with Fig. 7a and 8a (the APR model) were repeated for the APR-H model. Additional HMIM parameters for these simulations included: *w*_{A} = 0.5, *w*_{1m} = 0.5, *w*_{2m} = 0.5, α_{1} = 0.01 s^{−1}, α_{2} = 0.01 s^{−1}, ω_{1} = 0.001 s^{−1}, and ω_{2} = 0.001 s^{−1}. The influence of variations in α_{12} (0, 0.001, and 0.01 s^{−1}) and ω_{12} (0, 0.001, and 0.01 s^{−1}) are shown in Fig. 7b and 8b, respectively. The individual hydrographs and BTCs for active and passive regions of the APR model are equivalent to the UFT model when α_{12} = 0 and ω_{12} = 0, respectively. In contrast, the individual hydrographs and BTCs for the active and passive regions of the APR-H model are equivalent to the HMIM model results under these same conditions. The HMIM model can produce an earlier breakthrough and longer tailing in the hydrograph and BTC due to the presence of immobile zones in comparison with the uniform model. Consequently, the APR-H model can exhibit greater physical nonequilibrium characteristics than the APR model. An increase in α_{12} and ω_{12} influences the hydrographs and BTCs, respectively, in a similar manner for both APR and APR-H models. However, the APR-H does not always produce hydrographs and BTCs that are consistent with the uniform model (at the average flow velocity) when α_{12} and ω_{12} are very large due to the presence of immobile zones.

### Combined Active–Passive Regions and Vertical Mobile–Immobile Model

A simulation associated with Fig. 6a (the VMIM model) was repeated for the APR-V model in Fig. 6b. Additional APR parameters for this simulation included: *w*_{A} = 0.5, α_{12} = 0.1 s^{−1}, *n*_{1} = 0.01, *n*_{2} = 0.05, and *h*_{1im}^{max} = *h*_{2im}^{max} = 0.5 cm. Similar to the VMIM model, overland flow in the APR-V model does not occur until the depression storage zone is filled. Overland flow may then occur at faster and slower rates in the active and passive regions, respectively. The water depth will, therefore, be higher in the passive than the active domain for the same rainfall and infiltration rates in both regions. Consequently, the steady-state water depth at the outlet is higher in the APR-V model than in the VMIM model when α_{12} = 0. However, increases in α_{12} cause exchange of water between the passive and active regions such that the steady-state water depths approach each other when α_{12} is very high.

### Solute Wash-off

Figure 9a and 9b shows plots of a runoff hydrograph and solute wash-off from the impervious soil surface, respectively, when using the several model formulations (UFT, HMIM, APR, and APR-H). In this case, the same inflow (1-cm constant head for 6 min), runoff flow rate (45 cm^{2} s^{−1}), initial soil solute concentration (unit concentration), and solute desorption rate (0.01 s^{−1}) were considered. Other model parameters included: (i) the UFT model with *n* = 0.01; (ii) the HMIM model with *n* = 0.01, *w*_{m} = 0.5, α = 0.01, and ω_{m} = 0.001 s^{−1}; (iii) the APR model with *w*_{A} = 0.5, *n*_{1} = 0.005, *n*_{2} = 0.02, α_{12} = 0.001 s^{−1}, and ω_{12} = 0.001 s^{−1}; and (iv) the APR-H model with *w*_{A} = 0.5, *w*_{1m} = 0.5, *w*_{2m} = 0.5, *n*_{1} = 0.005, *n*_{2} = 0.02, α_{1} = 0.01 s^{−1}, α_{2} = 0.01 s^{−1}, ω_{1} = 0.01 s^{−1}, ω_{2} = 0.01 s^{−1}, α_{12} = 0.001 s^{−1}, and ω_{12} 0.001 s^{−1}. Drastically different hydrographs and solute wash-off behavior are observed when using the various physical nonequilibrium model formulations in comparison to the uniform model for the same initial and runoff flow rate. Clearly, the physical nonequilibrium models can produce hydrographs and/or solute wash-off with earlier or delayed arrivals, multiple peaks, and prolonged tailing.

### Validation of Physical Nonequilibrium Models

We next present visual evidence for physical nonequilibrium phenomena during overland flow and transport. The UFT and HMIM models are used to describe published tracer transport data from this system (Bradford et al., 2015), and a comparison is made between simulation results from several one-dimensional models (UFT, HMIM, and APR) and a two-dimensional model domain having spatial variations in the Manning roughness coefficient.

A dye tracer experiment was performed in a laboratory runoff chamber. The chamber was 2.25 m long, 0.15 m wide, and 0.16 m high. Autoclaved ultrapure quartz sand was uniformly packed into the chamber to a depth of 0.1 m. The chamber slope was then set to 11.8%. Steady-state water flow at a rate of 124 mL min^{−1} was achieved in the chamber before initiating a dye tracer experiment using a peristaltic pump connected to a rain simulator at the upslope portion of the chamber. A step pulse (30 min) of 10 mg L^{−1} yellow-colored fluorescent dye tracer was subsequently pumped through the rain simulator at the same flow rate. Figure 10 shows water and dye movement at the soil surface during this experiment. The water flow and dye were clearly not uniformly distributed (sheet flow) across the soil surface. The concentrated flow mostly occurred at the edge of the chamber, whereas other areas had little or no overland flow. These observations provide visual evidence that physical nonequilibrium processes were significantly contributing to overland water flow and tracer transport in this experiment.

Bradford et al. (2015) presented salt tracer (100 mM NaCl solution) data from this same runoff chamber setup when the slope was 5.6, 8.6, and 11.8%. In this case, the eluted solution was collected at the toe of the slope using a fraction collector every 3 min for 30 min, then the chamber was eluted with deionized water for another 90 min. The UFT and HMIM models were used to simulate the tracer BTCs when the Manning’s coefficient was 0.02, the dispersivity was taken as 0.1 of the chamber length (Gelhar et al., 1985, 1992), and other model parameters (*w*_{m} and ω_{m}) were determined by inverse optimization. The fitted parameters were *w*_{m} = 0.39 × 10^{−2}, 0.22 × 10^{−2}, and 0.12 × 10^{−2} and ω_{m} = 0.37 × 10^{−3}, 0.44 × 10^{−3}, and 0.59 × 10^{−3} s^{−1} for the slopes of 5.6, 8.6, and 11.8%, respectively. Figure 11a presents the observed and simulated BTCs. The values of the Pearson’s correlation coefficient (*R*^{2}) for the HMIM models at the slopes of 5.6, 8.6, and 11.8% were 0.94, 0.92, and 0.90, respectively. The simulation results indicated that the HMIM model provided a reasonable fit to the experimental data, whereas the UFT model was unable to accurately describe the initial pulse and the prolonged tailing behavior. The value of *w*_{m} (the ratio of the width of the mobile region and the total width of the soil surface) were generally very small and tended to decrease with increasing chamber slope. The higher slope generates more dynamic water flow and apparently creates a wider range of immobile regions to store the surface water. A greater rate of diffusive mass transfer (ω_{m}) also occurred with increasing slope. Figure 11b shows an enlargement of the simulated effluent tracer concentrations during the first 60 s. The HMIM model generates earlier tracer arrival than the UFT model. It should be mentioned that these differences in the arrival time between HMIM and UFT models will increase with the domain scale.

Most process-based overland flow and transport models neglect local-scale parameter heterogeneity by using average values of the slope and Manning roughness coefficient across the entire surface (Wallach and van Genuchten, 1990; Deng et al., 2005). This approach significantly reduces the input parameter requirements and computational time but may sometimes result in a poor description of the experimental data. The developed HMIM and APR models allowed us to simply represent spatial variability in flow and transport parameters in a simple one-dimensional domain with a limited number of parameters. As an illustration, simulation results from several one-dimensional models (UFT, HMIM, and APR) and a two-dimensional model domain having spatial variations in the Manning roughness coefficient were compared. In this case, HYDRUS-2D was similarly adapted as HYDRUS-1D to simulate uniform (sheet) overland flow and transport over a two-dimensional domain. A hypothetical, spatially varying roughness coefficient was generated using the HYDRUS-2D GUI in a 10- by 100-m simulation domain. The roughness coefficient at the edge of the domain was set to 0.005 and rest of the domain was set to 0.1 (Fig. 12a) to achieve spatial variations in the overland flow field. As expected, solute transport was much faster in the domain with the smaller roughness coefficient (Fig. 12b). The integrated BTCs at the bottom boundary of the two-dimensional domain were subsequently calculated and analyzed using various one-dimensional models (UFT, HMIM, and APR). Figure 13 shows the observed (average effluent concentration from the two-dimensional simulation) and simulated BTCs from the one-dimensional models. The one-dimensional UFT model with an averaged Manning’s roughness coefficient (*n =* 0.1) provided a poor description of the average two-dimensional data. In contrast, the HMIM and especially the APR models provided a much better description of the average two-dimensional data. The Manning’s coefficient set in the APR models were *n*_{1} *=* 0.005 and *n*_{2} *=* 0.1 for active domain and passive domain, respectively. The values of *w*_{A} and ω_{12} were determined by inverse optimization as 0.05 and 0.53 × 10^{−3} s^{−1}, respectively. For the HMIM model, Manning’s coefficient in the mobile domain set equal to 0.005, the inverse *w*_{m} and ω_{m} were equal to 0.09 and 3.50 × 10^{−3} s^{−1}, respectively. In contrast, the HMIM and especially the APR models provided a much better description of the average two-dimensional data. It should be mentioned that the distribution of Manning’s roughness coefficients in the two-dimensional simulation was selected to achieve physical nonequilibrium flow behavior. Additional research is needed to systematically study and assess the ability of one-dimensional physical nonequilibrium models to describe overland flow and transport on two- and/or three-dimensional domains with heterogeneous soil surface properties.

## Summary and Conclusions

The popular HYDRUS-1D code was extended to simulate uniform and physically nonequilibrium overland flow and reactive solute transport (such as salts, nutrients, pesticides, and microbes). This code provides information on the temporal and spatial distribution of water depths and solute concentrations in different phases (e.g., liquid and solid) and regions (e.g., mobile and immobile zones, active and passive regions), as well as on the mass balance of the soil surface, inflow, rainfall, infiltration, and evapotranspiration. These models provide a comprehensive set of tools to numerically investigate many important research problems involving overland flow and reactive transport processes. Physical nonequilibrium models may be better suited for studying hydrological processes at the plot and field scales than equilibrium models when spatial patterns of land surface characteristics are poorly characterized. However, the physical nonequilibrium models also involve a relatively large number of parameters, which may need to be obtained by calibrating against laboratory or field measurements. The HYDRUS-1D model includes provisions to estimate these parameters by inverse optimization as demonstrated above, and a variety of objective functions can be considered based on different measurements (e.g., water fluxes, water depths, resident concentrations, flux concentrations, etc.). Ideally, the model selection should be based on calibration with experimental data. If one model does not provide a satisfactory description of the data, then another model with an additional level of complexity may be considered. The simplest model (with the fewest number of fitting parameters) that accurately describes the experimental data is generally preferred.

Additional modifications to the code are needed to consider other infiltration equations and the full coupling between runoff water and the subsurface.

- Received November 18, 2016.
- Accepted May 6, 2017.

This is an open access article distributed under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).