Effective Saturated Hydraulic Conductivity for Representing Field-Scale Infiltration and Surface Soil Moisture in Heterogeneous Unsaturated Soils Subjected to Rainfall Events

Abstract

Spatial heterogeneity in soil properties has been a challenge for providing field-scale estimates of infiltration rates and surface soil moisture content over natural fields. In this study, we develop analytical expressions for effective saturated hydraulic conductivity for use with the Green-Ampt model to describe field-scale infiltration rates and evolution of surface soil moisture over unsaturated fields subjected to a rainfall event. The heterogeneity in soil properties is described by a log-normal distribution for surface saturated hydraulic conductivity. Comparisons between field-scale numerical and analytical simulation results for water movement in heterogeneous unsaturated soils show that the proposed expressions reproduce the evolution of surface soil moisture and infiltration rate with time. The analytical expressions hold promise for describing mean field infiltration rates and surface soil moisture evolution at field-scale over sandy loam and loamy sand soils.

1. Introduction

Determination of infiltration rates and surface soil moisture variation at a range of scales is useful in hydrologic and climatic modeling, and prediction studies [1,2]. The Richards equation is a physically-based model that is utilized to discern infiltration rates and soil moisture variation in unsaturated soils at field-scale. This equation is highly non-linear and imposes a heavy computational burden, especially when describing water movement in initially dry soils [3]. Moreover, detailed data on soil properties in the study area are needed for implementing the Richards equation at its full potential. Such data are rarely available, and studies often resort to calibration exercises to estimate unknown quantities.

The description of infiltration rates and surface soil moisture at different spatial and temporal scales is complicated further by the spatial heterogeneity exhibited by the hydraulic properties of the field-soils [4,5,6]. To estimate the spatial and temporal variation of surface soil moisture content and infiltration rates, ground-based measurements are most commonly used. These measurements provide information on soil moisture data and infiltration rates at a local-scale [7], even though conflicting results from in situ-measurements of saturated hydraulic conductivity by different devices have been obtained. Moreover, the measurements are expensive, time consuming, provide data at few selected points, and require accurate calibration of the measuring instrument. Further, reliable interpretation of spatial and temporal patterns of infiltration rates and surface soil moisture at different scales from measurements available at few selected points is very challenging [8].

An alternate approach for predicting mean flow behavior at the field-scale has been upscaling i.e., to represent the heterogeneous domain by an equivalent homogeneous domain with upscaled (effective or macroscopic) flow properties [9,10]. The upscaling procedures presently available can be broadly categorized as (1) forward upscaling and (2) inverse upscaling.

Forward upscaling relies on deriving the effective parameters and upscaled flow equations by considering information about the spatial structure and variability of soil hydraulic properties measured at the local-scale [10,11]. Stochastic perturbation-based methods [12,13,14,15,16] and stochastic streamtube-based approaches [17,18,19,20,21] have been widely popular in this regard. Perturbation-based methods are constrained to cases of small variations in soil properties, and the resulting equations for effective soil hydraulic parameters and field-scale flow are computationally intensive because of high dimensionality, intercoupling, and nonlinearity [10]. Stochastic streamtube models are better suited for highly heterogeneous soils. However, obtaining field-scale solutions from a local-scale solution in a stochastic streamtube framework often requires computationally intensive Monte-Carlo simulations [22].

Inverse scaling methods utilize information about the measurement and observations of spatial variability of soil properties, state variables of soil water processes (e.g., soil moisture content, or soil water tension) and boundary fluxes at various scales to derive effective parameters of the heterogeneous system at a large scale [9,10]. These methods mainly assume that the effective flow equations are known a-priori at the large-scale. Few studies have also derived field-scale effective properties using inverse scaling methods [23,24,25]. Yeh et al. [25] and Ye et al. [24] used spatial moment analysis to derive expressions for three-dimensional effective unsaturated hydraulic conductivity tensor from soil moisture observations under transient flow-conditions. Ward et al. [23] applied effective medium approximation method to estimate upscaled water retention characteristics and effective unsaturated hydraulic conductivity from the observed moisture distribution for a field-scale water injection test. The proposed inverse scaling methods for estimating effective parameters at field-scale require extensive soil moisture measurements and are suited for low values of soil moisture.

Many watershed models, e.g., SWAT [26], GSSHA [27], and KINEROS2 [28] partition the study area into smaller “field-scale” units and assume homogeneous properties within each unit. When using such models, a fundamental question of incorporating the effects of subgrid heterogeneity needs to be addressed. A viable approach is to replace the heterogeneous soil system with equivalent soil properties so that a deterministic homogeneous solution mimics field-scale behavior with sufficient accuracy. Since a straightforward statistical average does not meet this need when the governing processes and transformations are nonlinear, procedures for estimation of effective saturated soil hydraulic conductivity in such contexts is useful, and is the subject of this paper.

Additionally, it is noteworthy that most of the field-scale studies have focused on estimating analytical expressions for effective soil hydraulic parameters for steady or gravity dominated flow conditions [29]. However, for field conditions where the soil is subjected to rainfall or irrigation events, the boundary conditions change from being rainfall-controlled to soil-controlled in different parts of the field, and transient flow conditions are arguably just as prevalent [30]. Obtaining analytical expressions for effective soil hydraulic properties at field-scale for non-linear Richards equation under such a boundary condition is more complex, and has not received enough attention in the literature. The KINEROS2 model does have a feature for representing the spatial variability of $K_s$, but this feature needs to be tested.

As upscaling of Richards equation under such dynamic boundary conditions is complex, approximate local models are sought for enabling mathematically-tractable field-scale solutions [31]. When seeking field-scale solutions, researchers note that errors incurred by using simplified local models tend to cancel out due to large-scale variability in soil parameters [17,18,32,33,34]. Under this hypothesis, the evolution of surface soil moisture and infiltration rates at a local-scale under a single rainfall event are described using a sharp front model for water movement. Subsequently, utilizing the field-scale spatial variability in saturated hydraulic conductivity, analytical solutions for effective saturated hydraulic conductivity are developed.

We assume that the 3-D soil field is represented as a collection of independent streamtubes with only vertical flow, and the saturated hydraulic conductivity, $K_s$, follows a lognormal distribution over the field. Other soil parameters of Brooks and Corey model, $\mathsf{\lambda },{\mathsf{\theta }}_s,{\mathsf{\theta }}_i,{\mathsf{\psi }}_b$ (see below Equation (2)) are considered as constants. The streamtube model is a good approximation at shallow depth under high recharge rates, and when the horizontal correlation length scale is much larger than the vertical length scale [17,18,20,31,35,36,37,38]. The surface saturated hydraulic conductivity has been found to exhibit a variation that is typically represented by a lognormal distribution, e.g., [17,18,22,35]. Specifically, Russo and Bresler [39] found that the variation in parameters, $\mathsf{\lambda },{\mathsf{\theta }}_s,{\mathsf{\theta }}_i,{\mathsf{\psi }}_b$, to be smaller than that associated with $K_s$, a fact that is supported by experimental studies as well, e.g., [40,41].

The primary goal of the study is to develop analytical expressions for effective saturated hydraulic conductivity to describe the temporal evolution of surface soil moisture and infiltration rates under a rainfall event. These expressions can be used to mimic field-scale variability in current watershed models (such as SWAT, KINEROS, GSSHA and others) that utilize the Green-Ampt method to describe infiltration rate in homogeneous soils. First, local-scale models of infiltration rate and soil moisture are described from a sharp-front model for water movement. Then, assuming a lognormally distributed horizontal random field for saturated hydraulic conductivity, analytical expressions for the effective saturated hydraulic conductivity are developed. Comparisons of the results obtained from this model with extensive sets of Monte Carlo simulations, 3-D numerical solutions and field-scale analytical results are performed to evaluate the applicability of the expressions for effective hydraulic conductivity.

2. Materials and Methods

The governing differential equations and notation in this section are similar to those of previous studies [30,31]. The one-dimensional vertical Richards equation for water movement in unsaturated soils is,

$$\frac{\partial \mathsf{\theta }}{\partial t}-\frac{\partial }{\partial z}\left[K(\mathsf{\psi })\left(\frac{\partial \mathsf{\psi }}{\partial z}-1\right)\right]=0$$

(1)

where $\mathsf{\theta }$ is the volumetric water content, $t$ is time and $z$ is the vertical coordinate pointing downwards. The Brooks and Corey model [42] is often used to describe unsaturated soil properties in Equation (1). The soil-water saturation, pressure head, and conductivity are related as,

$$K_r(\mathsf{\psi })=\frac{K(\mathsf{\psi })}{K_s}={\left(\frac{{\mathsf{\psi }}_b}{\mathsf{\psi }}\right)}^n;S=\frac{\mathsf{\theta }-{\mathsf{\theta }}_r}{{\mathsf{\theta }}_s-{\mathsf{\theta }}_r}={\left(\frac{{\mathsf{\psi }}_b}{\mathsf{\psi }}\right)}^{\mathsf{\lambda }}\mathsf{\psi }>{\mathsf{\psi }}_b;n=2+3\mathsf{\lambda }$$

(2)

where the conductivity $K(\mathsf{\psi })$ depends on the capillary pressure head, $\mathsf{\psi }$, saturated hydraulic conductivity, $K_s$, air entry capillary pressure head, ${\mathsf{\psi }}_b$, and the pore size distribution index, $\mathsf{\lambda }$. Here ${\mathsf{\theta }}_s$ and ${\mathsf{\theta }}_r$ represent the saturated and residual water contents, respectively. Note that the present analysis is only applicable when the soil characteristic curves are described by the Brooks and Corey model.

Under the sharp-front model for water movement, the cumulative infiltrated depth of water into the soil for constant rainfall rate $r$ is [43],

$$\begin{array}{ll}F=rt,& t\le t_p\\ F=F_p+\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }\ln \left(\frac{\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }+F}{\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }+F_p}\right)+K_s\left(t-t_p\right),& t>t_p\end{array}$$

(3)

In Equation (3), $\mathrm{\Delta }\mathsf{\psi }$ is the capillary pressure head difference between the soil surface and the front, and can be computed as $\mathrm{\Delta }\mathsf{\psi }={\mathsf{\psi }}_b/\left(1-\mathsf{\eta }\right)$ [18], and $\mathrm{\Delta }\mathsf{\theta }=\left({\mathsf{\theta }}_s-{\mathsf{\theta }}_i\right)$ is the difference between the water content at saturation and the initial water content. The subscript $p$ refers to conditions at ponding. The ponding time $t_p$ and cumulative infiltrated depth of water at ponding are given by

$$F_p=\frac{K_s\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }}{r-K_s};t_p=\frac{F_p}{r}r>K_s$$

(4)

The following dimensionless variables are defined.

$$\begin{array}{l}{\mathsf{\theta }}_t^{\ast }=\frac{{\mathsf{\theta }}_t-{\mathsf{\theta }}_i}{{\mathsf{\theta }}_s-{\mathsf{\theta }}_i};t^{\ast }=\frac{tr}{\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }};y=\frac{K_s}{r};F^{\ast }=\frac{F}{\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }};\\ t_p^{\ast }=\frac{t_{p}r}{\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }}=\frac{y}{1-y}\\ F_p^{\ast }=\frac{F_p}{\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }}=\frac{y}{1-y}\end{array}$$

(5)

where, ${\mathsf{\theta }}_t$ is the surface soil moisture at any time $t$.

Field-Scale Solution

At the field scale, the quantity $y$ is assumed to vary spatially following a log-normal distribution with the density function.

$$f_y\left(y;{\mathsf{\mu }}_Y,{\mathsf{\sigma }}_Y\right)=\frac{1}{\sqrt{2\mathsf{\pi }}{\mathsf{\sigma }}_{Y}y}\exp \left[-\frac{{\left(\ln y-{\mathsf{\mu }}_Y\right)}^2}{2{\mathsf{\sigma }}_Y^2}\right]$$

(6)

where, ${\mathsf{\mu }}_Y$ and ${\mathsf{\sigma }}_Y$ are the mean and standard deviation of $Y=\ln y=\ln \left(K_s/r\right)$.

The truncated moments of the log-normal distribution are [4]:

$$\begin{array}{ll}{G}_Y\left(y_c,\mathsf{\zeta }\right)& ={\displaystyle \underset{0}{\overset{y_c}{\int }}{y}^{\mathsf{\zeta }}{f}_y\left(y\right)dy}\\ & =\exp \left(\mathsf{\zeta }{\mathsf{\mu }}_Y+\frac{{\mathsf{\sigma }}_Y^2{\mathsf{\zeta }}^2}{2}\right)\left[1-\frac{1}{2}erfc\left(\frac{\ln y_c-{\mathsf{\mu }}_Y}{\sqrt{2}{\mathsf{\sigma }}_Y}-\frac{{\mathsf{\sigma }}_Y\mathsf{\zeta }}{\sqrt{2}}\right)\right]\end{array}$$

(7)

where $erfc\left(\right)$ is the complimentary error function.

Effective Value for Replicating Field-Scale Infiltration Rates

Often, field-scale rainfall and runoff data are available from which field-scale infiltration rates can be computed. The problem of estimating the effective value of $K_s$ that will honor these observations is challenging because local-scale infiltration, $I$ is known to be highly nonlinear, with infiltration rates governed by rainfall rate until ponding time, and then controlled by soil properties. When developing an expression for effective saturated hydraulic conductivity $K_s$ (or y), we are governed by preserving the correct limiting properties so that the infiltration rates obtained from effective properties using the deterministic Green-Ampt model replicate field-scale infiltration rates at small and large times. By utilizing the theories developed in earlier studies of Govindaraju et al. [4] the following expression for effective value of $y_{eff}$ is proposed to incorporate spatial variability of soil hydraulic properties for predicting infiltrations rates at the field-scale.

$$y_{eff}=\frac{K_{s-eff}}{r}=\left[1-G_Y\left(1,0\right)\right]+G_Y\left(1,1\right)$$

(8)

Effective Value for Replicating Surface Soil Moisture

Because the Green-Ampt equation is an approximate form of the Richards equation, the effective value of y in Equation (8), while suitable for field-scale infiltration rates, does not replicate field-scale surface soil moisture well, and a different tack has to be employed. The expression for effective saturated hydraulic conductivity that can describe mean field-scale surface soil moisture evolution for any duration of rainfall $t_r$, and $t_r^{\ast }=\frac{t_{r}r}{\mathrm{\Delta }\mathsf{\psi }\mathrm{\Delta }\mathsf{\theta }}$ is obtained by equating the water content values at local and field-scale at the cessation of rainfall. Thus, the evolution of soil moisture at local-scale just prior to cessation of rain ($t=t_r$) can be expressed as,

$${\mathsf{\theta }}_{tr}^{\ast }={\left\{\frac{1}{y_{eff}}\left[\frac{t_r^{\ast }}{t_r^{\ast }+1}\right]\right\}}^{2\mathsf{\lambda }/(7\mathsf{\lambda }+3)}$$

(9)

The mean-field soil moisture content just prior to cessation of rain ($t=t_r$) can be evaluated (see Appendix A.1.2) as,

$$E\left[{\mathsf{\theta }}_{tr}^{\ast }\right]={\left(\frac{t_r^{\ast }}{t_r^{\ast }+1}\right)}^{2\mathsf{\lambda }/(7\mathsf{\lambda }+3)}\left[G\left(\infty ,\frac{-2\mathsf{\lambda }}{7\mathsf{\lambda }+3}\right)-G\left(K_{cr}^{\ast },\frac{-2\mathsf{\lambda }}{7\mathsf{\lambda }+3}\right)\right]+G\left(K_{cr}^{\ast },0\right)$$

(10)

where $K_{cr}^{\ast }=\frac{t_r^{\ast }}{1+t_r^{\ast }}$ By setting ${\mathsf{\theta }}_{tr}^{\ast }=E\left[{\mathsf{\theta }}_{tr}^{\ast }\right]$ in Equations (10) and (11), we obtain $y_{eff}$ as,

$$y_{eff=}{\left[G\left(\infty ,\frac{-2\mathsf{\lambda }}{7\mathsf{\lambda }+3}\right)-G\left(K_{cr}^{\ast },\frac{-2\mathsf{\lambda }}{7\mathsf{\lambda }+3}\right)+G\left(K_{cr}^{\ast },0\right){\left(\frac{t_r^{\ast }}{t_r^{\ast }+1}\right)}^{-2\mathsf{\lambda }/(7\mathsf{\lambda }+3)}\right]}^{-\left(\left(7\mathsf{\lambda }+3\right)/2\mathsf{\lambda }\right)}$$

(11)

This value of $y_{eff}$ holds for $t>t_r$ as well and is more suited for initially dry sandy loam soils $\left({\mathsf{\theta }}_i\approx {\mathsf{\theta }}_r\approx 0;\mathsf{\lambda }\approx 1\right)$.

3. Design of Numerical Schemes

The analytical expressions for $y_{eff}$ developed in this study were used in the approximate 1-D analytical Green-Ampt model and 1-D Richards equation for field-scale surface soil moisture and infiltration rates. For comparison the following numerical schemes were also utilized:

(a)

Numerical-FD: A numerical scheme denoted as “Numerical-FD” was implemented for Monte-Carlo simulations. In this scheme, from a prescribed lognormal distribution, a random value of $K_s$ was drawn. Further, the numerical solution of the Richards equation (Equations (1)–(3)) was obtained using an implicit finite difference scheme. The water content, $\mathsf{\theta }\left(z,t\right)$, fields so generated were stored for this single realization of $K_s$, and 2000 such Monte Carlo simulations were utilized to obtain the field-scale surface soil moisture, $E\left[{\mathsf{\theta }}_t\right]$ and infiltration rates, $E\left[I\right]$.

(b)

Numerical 3DSH: A 3-D layered domain of dimensions 12 m × 12 m × 3 m was set up in Hydrus-3D [44], with 80 horizontal layers. Finite element mesh was generated, with fine spacing near the soil surface to capture the sharp fronts. In the horizontal plane the average size of triangular elements was 25 cm, leading to a total of 3691 finite element nodes in one horizontal layer. The layered 3-D model allows for setting up a spatially heterogeneous and vertically homogeneous or heterogeneous soil models. Initially, a 2-D spatially heterogeneous domain with the specified parameters for $\ln K_s$ was generated. To set-up a spatially heterogeneous and vertically homogeneous physical domain, the generated values of $\ln K_s$ were assigned to each element for all horizontal layers of the 3D-domain. A no-flux boundary condition was assigned to the vertical side walls and at the bottom boundary. At the soil surface, before ponding is achieved, the boundary condition is of surface flux equal to the constant rainfall intensity. After ponding, a fixed pressure head of zero is applied at the surface node.

(c)

Numerical-3D: The above mentioned 1-D and 3-D numerical schemes only consider horizontal heterogeneity; whereas, in realistic field-conditions heterogeneity exists in both horizontal and vertical directions. Therefore, a 3-D numerical scheme in Hydrus-3D, with heterogeneity in both horizontal and vertical directions was also considered for comparison. For this scheme, uncorrelated saturated hydraulic conductivity values were generated in both horizontal and vertical directions with the specified parameters of $\ln K_s$. The dimensions of the 3-D domain, mesh-spacing, parameters and boundary conditions are similar to those adopted for “Numerical-3DSH” scheme.

For the 3-D numerical solutions, in Hydrus-3D the governing flow equation is given by the modified form of Richards equation [44]

$$\frac{\partial \mathsf{\theta }}{\partial t}=\frac{\partial }{\partial x_i}\left[K\left(K_{ij}^A\frac{\partial \mathsf{\psi }}{\partial x_i}+K_{iz}^A\right)\right]$$

(12)

where, $\mathsf{\psi }$ is the pressure head, $x_i\left(1,2\right)$ are the spatial coordinates, t is time, $K_{ij}^A$ are components of a dimensionless anisotropy tensor $K^A$, and $K$ is the unsaturated hydraulic conductivity.

Given

$$K\left(h,x,y,z\right)=K_s\left(x,y,z\right)K_r\left(h,x,y,z\right)$$

(13)

$K_r$ is the relative hydraulic conductivity and $K_s$ the saturated hydraulic conductivity.

The effect of hysteresis has been ignored and the soil water retention $\mathsf{\theta }\left(\mathsf{\psi }\right)$, and hydraulic conductivity function $K\left(\mathsf{\psi }\right)$ according to Brooks and Corey [42] are given by Equation (2).

These numerical schemes are the same as those used in Ojha and Govindaraju [31] to test upscaling strategies.

4. Results

The ability of the analytical expressions of $y_{eff}$ (Equations (8) and (11)) proposed in this study for predicting field-scale infiltration rates and surface soil moisture contents were tested with the help of numerical simulations and analytical expressions for field-scale surface soil moisture evolution [31] and infiltration rates [4].

4.1. Surface Soil Moisture

4.1.1. Comparisons with Numerical Results

For the numerical experiment, the following parameter values were considered: ${\mathsf{\psi }}_b=-34.22\mathrm{cm}$, ${\mathsf{\theta }}_r=0.01$, ${\mathsf{\theta }}_s=0.415$, ${\mathsf{\theta }}_i=0.06$, $r=0.5\mathrm{cm}/\mathrm{h}$, $t_r=2\mathrm{h}$, $\mathsf{\lambda }=0.91$, ${\mathsf{\mu }}_{\ln K_s}=2$, and ${\mathsf{\sigma }}_{\ln K_s}=0.2$ in Equations (1) and (2). $K_s$ has units of $\mathrm{cm}/\mathrm{h}$. Soils with these properties are characterized as sandy loams [45].

The ability of the derived analytical expressions to predict field-scale surface soil moisture variation was validated for different parameter sets against numerical solutions. The mean field-scale surface soil moisture evolution obtained using $y_{eff}$ in 1-D Richards equation is denoted as “Analytical”. Figure 1 shows comparisons between mean field-scale surface moisture profiles generated for three different values of $\mathsf{\lambda }$, ${\mathsf{\mu }}_{\ln K_s}=2$, ${\mathsf{\sigma }}_{\ln K_s}=0.2$, and $r=0.5\mathrm{cm}/\mathrm{h}$. A good agreement is observed between the analytical and numerical schemes for all values of $\mathsf{\lambda }$. Though the analytical expressions are more suited for $\mathsf{\lambda }\approx 1$, for smaller values of $\mathsf{\lambda }$, the conformity between analytical and 3D-numerical solutions is relatively higher for redistribution conditions, and the agreement improves with time. This is probably due to dominance of vertical flow conditions for lower values of $\mathsf{\lambda }$ in 3-D soil fields. The analytical model reproduces the general trend for surface soil moisture evolution as obtained from the three numerical schemes.

Earlier studies have suggested that the sharp front approximations perform better with increase in field-scale variability, as the local errors due to sharp front approximation cancel out during the averaging process. Therefore, the theoretical results are also validated with numerical results for changes in ${\mathsf{\sigma }}_{\ln K_s}$, as shown in Figure 2. The mean saturated hydraulic conductivity, ${\mathsf{\mu }}_{\ln K_s}=2$, $\mathsf{\lambda }=0.91$ and other soil parameters are kept constant. A conformity between theoretical model and numerical results is good as in case of ${\mathsf{\sigma }}_{\ln K_s}=0.2$ and ${\mathsf{\sigma }}_{\ln K_s}=0.4$; however, for ${\mathsf{\sigma }}_{\ln K_s}=0.6$ the deviation between the results increases for the rainfall duration. For higher variability in ${\mathsf{\sigma }}_{\ln K_s}$, 3-D effects become more prominent, promoting lateral flows thus causing analytical model assumptions to be less met.

The performance of the results from the analytical expressions was also tested for changes in ${\mathsf{\mu }}_{\ln K_s}$ (with other soil parameters kept constant) by comparisons with numerical results. Figure 3 shows comparison between theoretical model results for varying ${\mathsf{\mu }}_{\ln K_s}$, with ${\mathsf{\sigma }}_{\ln K_s}=0.2$, $\mathsf{\lambda }=0.91$ and $r=0.5\mathrm{cm}/\mathrm{h}$. A very good agreement is observed between the theoretical results and the numerical schemes for the rainfall duration as ${\mathsf{\mu }}_{\ln K_s}$ increases. For higher values of ${\mathsf{\mu }}_{\ln K_s}$, the soil drains faster in independent streamtube columns, thus reducing 3-D lateral effect. The soil moisture content values are higher for low value of ${\mathsf{\mu }}_{\ln K_s}=1.5$, and decrease with increase in ${\mathsf{\mu }}_{\ln K_s}$. With increase in ${\mathsf{\mu }}_{\ln K_s}$, a greater fraction of soil-surface is governed by rainfall-limited conditions and pre-ponding conditions dominate at the field-scale over a longer duration. Overall, the expression developed for $y_{eff}$ in this study reproduces the field-scale surface soil moisture well when compared to numerical results.

4.1.2. Comparisons with Analytical Results

The mean field-scale surface moisture profiles generated using $y_{eff}$ are compared with field-scale moisture profiles obtained from the analytical expression (Equations (A6) and (A7) in the Appendix A) proposed by Ojha and Govindaraju [31] for different set of parameters. The mean field-scale surface soil moisture profile generated using $y_{eff}$ in the analytical expression for local-scale surface soil moisture variation derived from sharp-front approximation (Equations (A1)–(A4)) is denoted as “A-Local”. The mean field-scale surface soil moisture generated using Equations (A6) and (A7) in Appendix A is denoted as “A-Field”.

Figure 4 shows comparisons between the analytical results for varying ${\mathsf{\mu }}_{\ln K_s}$, with ${\mathsf{\sigma }}_{\ln K_s}=0.2$, $\mathsf{\lambda }=0.91$ and $r=0.5\mathrm{cm}/\mathrm{h}$. The other parameters are similar to those used for the numerical experiments. A very good agreement is observed between both the theoretical results (i.e., “A-Local and A-Field”) prior to cessation of rainfall for all values of ${\mathsf{\mu }}_{\ln K_s}$. As the value of ${\mathsf{\mu }}_{\ln K_s}$ increases, the conformity between the two analytical results increases during redistribution. Further, as the value of ${\mathsf{\mu }}_{\ln K_s}$ increases, a relatively faster drop in the values of $E\left[{\mathsf{\theta }}_t\right]$ for redistribution conditions is observed. When $K_s$ is high, the soil drains faster, resulting in low moisture content values at the soil surface after rainfall stops.

The proposed analytical expression for $y_{eff}$ is more suited for sandy-loam soils. However, the efficacy of the developed analytical expression for $y_{eff}$ in predicting mean field-scale surface soil moisture of soils with very low $K_s$ (higher clay content) is also examined. The following parameter values are considered: ${\mathsf{\psi }}_b=-1\mathrm{cm}$, ${\mathsf{\theta }}_r=0.01$, ${\mathsf{\theta }}_s=0.415$, ${\mathsf{\theta }}_i=0.01$, $r=1\mathrm{cm}/\mathrm{h}$, $t_r=1\mathrm{h}$, $\mathsf{\lambda }=1.1$, and ${\mathsf{\mu }}_{\ln K_s}=-0.5$. Figure 5 shows comparison between theoretical model results for varying ${\mathsf{\sigma }}_{\ln K_s}$. For these soil properties, the exact form of the analytical expression (Equation (14)) for mean field-scale surface soil moisture proposed by Ojha and Govindaraju [31] is used. The soil moisture profile is generated by considering the first 50 terms of the series expansion.

$$\begin{array}{l}E\left[{\mathsf{\theta }}_t^{\ast }\right]={\left(\frac{t_r^{\ast }}{t_r^{\ast }+1}\right)}^{2\mathsf{\lambda }/(7\mathsf{\lambda }+3)}{\left[D\right]}^{-\mathsf{\lambda }/\mathsf{\eta }}{\displaystyle \sum _{i=0}^{\infty }\left(\begin{array}{l}-\mathsf{\lambda }/\mathsf{\eta }\\ i\end{array}\right)}{D}^{-i}\left[\begin{array}{l}G\left(\infty ,-\left[i\left(\frac{\mathsf{\lambda }-1}{7\mathsf{\lambda }+3}\right)+\frac{\mathsf{\lambda }}{\mathsf{\eta }}\right]\right)\\ -G\left(y^{\prime }+K_{c{t}_r}^{\ast },-\left[i\left(\frac{\mathsf{\lambda }-1}{7\mathsf{\lambda }+3}\right)+\frac{\mathsf{\lambda }}{\mathsf{\eta }}\right]\right)\end{array}\right]+\\ {\left(\frac{t_r^{\ast }}{t_r^{\ast }+1}\right)}^{2\mathsf{\lambda }/(7\mathsf{\lambda }+3)}{\displaystyle \sum _{i=0}^{\infty }\left(\begin{array}{l}-\mathsf{\lambda }/\mathsf{\eta }\\ i\end{array}\right)}{D}^i\left[\begin{array}{l}G\left(y^{\prime }+K_{c{t}_r}^{\ast },i\left(\frac{\mathsf{\lambda }-1}{7\mathsf{\lambda }+3}\right)-\frac{2\mathsf{\lambda }}{7\mathsf{\lambda }+3}\right)\\ -G\left(K_{c{t}_r}^{\ast },i\left(\frac{\mathsf{\lambda }-1}{7\mathsf{\lambda }+3}\right)-\frac{2\mathsf{\lambda }}{7\mathsf{\lambda }+3}\right)\end{array}\right]+\\ {\left[1+\frac{\mathsf{\eta }}{\mathsf{\lambda }}\left(\frac{t_r-t_r^{\ast }}{1+t_r^{\ast }}\right)\right]}^{-\mathsf{\lambda }/\mathsf{\eta }}G\left(K_{c{t}_r}^{\ast },0\right)\end{array}$$

(14)

where

$$\begin{array}{l}D=\frac{\mathsf{\eta }}{\mathsf{\lambda }}{\left(\frac{t_r^{\ast }}{t_r^{\ast }+1}\right)}^{\frac{6\mathsf{\lambda }+4}{7\mathsf{\lambda }+3}}\left[\frac{t^{\ast }}{t_r^{\ast }}-1\right]\\ y^{\prime }={\left(\frac{1}{D}\right)}^{\frac{7\mathsf{\lambda }+3}{\mathsf{\lambda }-1}}\end{array}$$

(15)

With decrease in ${\mathsf{\sigma }}_{\ln K_s}$, the agreement between “A-Local” and “A-Field” reduces for the rainfall duration. Further, for the same ${\mathsf{\mu }}_{\ln K_s}$ with increase in ${\mathsf{\sigma }}_{\ln K_s}$ post-ponding conditions gain prominence, and the magnitude of soil moisture content is high. After the rainfall stops, for a lower value of ${\mathsf{\mu }}_{\ln K_s}$ a gradual decrease in $E\left[{\mathsf{\theta }}_t\right]$ is observed, with the rate of decrease reducing with increase in ${\mathsf{\sigma }}_{\ln K_s}$.

4.2. Infiltration Rates

4.2.1. Comparisons with Numerical Results

The parameters used for the numerical validation of mean field-scale surface soil moisture were also utilized to validate the ability of the derived theoretical expression for $y_{eff}$ (Equation (8)) in predicting mean field-scale infiltration rates. While many cases with different parameter combinations were considered, Figure 6 shows typical results for the case where $\mathsf{\lambda }=0.91$, ${\mathsf{\sigma }}_{\ln K_s}=0.2$, and $r=0.5\mathrm{cm}/\mathrm{h}$, and three different values of mean (a) ${\mathsf{\mu }}_{\ln K_s}=1.5$, (b) ${\mathsf{\mu }}_{\ln K_s}=2.0$ and (c) ${\mathsf{\mu }}_{\ln K_s}=2.5$. The theoretical model reproduces the mean field-infiltration rates predicted by the numerical models well. For the rainfall duration, pre-ponding conditions dominate at the field-scale and the $E\left[I\right]$ is equal to the rainfall rate for different values of ${\mathsf{\mu }}_{\ln K_s}$.

4.2.2. Comparisons with Analytical Results

Similar to surface soil moisture, the mean field-scale infiltration rates generated using $y_{eff}$ are compared with mean field-scale infiltration rates obtained from the analytical expression (Equation (A8) in the Appendix A) proposed by Govindaraju et al. [4]. The curve for mean field-scale infiltration rates, generated using $y_{eff}$ in the analytical expression for local-scale infiltration rate variation derived from sharp-front approximation (Equation (3)), is denoted as “A-Local”. The mean field-scale infiltration rates generated using Equation (A8) in Appendix A are denoted as “A-Field”.

Figure 7 shows comparisons between the analytical results for varying ${\mathsf{\sigma }}_{\ln K_s}$, while the other parameter values considered are: ${\mathsf{\psi }}_b=-1\mathrm{cm}$, ${\mathsf{\theta }}_r=0.01$, ${\mathsf{\theta }}_s=1$, ${\mathsf{\theta }}_i=0.01$, $r=1\mathrm{cm}/\mathrm{h}$, $t_r=28\mathrm{h}$, $\mathsf{\lambda }=1.1$, and ${\mathsf{\mu }}_{\ln K_s}=-0.5$. Using Equation (8), the values of $y_{eff}$ for constant mean ${\mathsf{\mu }}_{\ln K_s}=-0.5$ and different standard deviations, ${\mathsf{\sigma }}_{\ln K_s}=1$, ${\mathsf{\sigma }}_{\ln K_s}=0.5$, ${\mathsf{\sigma }}_{\ln K_s}=0.1$ are obtained as 0.617, 0.634 and 0.61, respectively. The corresponding time to ponding, $t_p$ values are 1.61, 1.73, and 1.56 h, respectively. The parameters were selected such that the ability of the proposed expression for $y_{eff}$ in predicting the infiltration rates under rainfall at small and large times can be validated. For ${\mathsf{\sigma }}_{\ln K_s}=1$, and $t\le 5\mathrm{h}$ the analytical “A-Local” model over-predicts the mean field-scale infiltration rates (“A-Field”). In case of ${\mathsf{\sigma }}_{\ln K_s}=0.5$, initially the “A-Local” and “A-Field” models show very good agreement, but as time progresses for $t\le 5\mathrm{h}$ the infiltration rate values are higher in magnitude for the “A-Local” model. The mean field-scale infiltration rate values are very well predicted by the “A-Local” model for ${\mathsf{\sigma }}_{\ln K_s}=0.1$. Further, at large times, i.e., $t\to \infty$, the agreement between infiltration rate values for “A-Local” and “A-Field” is very high. The fractions of the field where pre- and post-ponding conditions exist are governed by the distributional properties of $K_s$ and are easily determined. With increase in value of ${\mathsf{\sigma }}_{\ln K_s}$, the control shifts from pre-ponding conditions to post-ponding conditions resulting in a decrease in value of $E\left[I\right]$ at the field-scale at short times (as observed for “A-Field”). However, as time progresses the entire field experiences post-ponding conditions irrespective of the field-scale variability, and the infiltration rate values match well for both “A-Local” and “A-Field”.

5. Conclusions

In this study, analytical expressions for effective saturated hydraulic conductivity were developed for a single rainfall event. The proposed effective saturated hydraulic conductivity expressions, when combined with 1-D Richards equation or local-scale simplified models based on sharp front approximation can be utilized to determine field-scale mean-surface moisture content and infiltration rates. The dimensionless equations for effective saturated hydraulic conductivity were derived from a sharp front model. Comparisons with 1-D Monte Carlo-based simulations, 3-D numerical simulations, and analytical results proposed by Govindaraju et al. [4] and Ojha and Govindaraju [31] showed that the proposed theory provides a simple alternative to determining the surface soil moisture and infiltration rates in the field for a single rainfall event, with insights into role of rainfall- and soil-controlled behavior for pre- and post-ponded conditions. The theoretical model readily describes field-scale surface soil moisture and infiltration rates where otherwise extensive amount of numerical computations would be needed.

Given the assumptions, the proposed analytical model for infiltration rates is better suited for recharge rates that are significantly greater than the mean saturated hydraulic conductivity. The field size should be much larger than the horizontal correlation length scale for meaningful spatial averages. The analytical model over-predicts the infiltration rate values at short time as field-scale variability increases. The numerical computation of the right hand side of Equation (7) also needs special care as pointed out by van Genuchten [46]. Nevertheless, the theoretical model proposed in this study provides a useful alternative for describing mean field-scale surface soil moisture and infiltration rates in heterogeneous unsaturated soils.

As noted, expressions for effective properties (Equations (8) and (11)) are better suited for sandy loam soils where the approximation of a sharp wetting front is better met than clayey soils where the wetting front tends to extend over greater depths. The assumption of initially dry soils allows us to neglect the very small conductivity values associated with initial conditions. Our experience with numerical simulations shows that these are not strict constraints, as the local deviations tend to cancel out in the face of variability resulting from field scale heterogeneity.

There are important implications when using effective saturated hydraulic conductivity to simulate field-scale behavior. First, the effective value is dependent on the rainfall intensity (see also KINEROS2) and not solely a function of $K_s$ variability as one would expect from using equivalent media theories. This is a consequence of the infiltration process. Thus, if no rainfall excess and subsequent runoff are generated over a field under a rainfall event, the effective value of $K_s$ is not identifiable as any value greater than rainfall intensity will yield the same result. If the rainfall intensity is varying in time, then the effective conductivity to be used in Green-Ampt model for field-scale infiltration rate will also need to be varying with time. This is a more complex problem and deserves future attention.

A second important consideration revealed here is the need for one effective $K_s$ value to replicate field-scale runoff data, and another one for reproducing field-scale surface soil moisture. This is an artifact of the use of the Green-Ampt model, an approximation that is being used to mimic field-scale behavior as determined by Richards equation. Specifically, the effective $K_s$ from runoff data does not replicate field-scale surface moisture and vice versa. Conversely, if runoff data are used to calibrate for an effective $K_s$ value, this estimate will be different from the value obtained from calibrating an effective $K_s$ value based on estimates of surface soil moisture at numerous locations. If runoff and surface soil moisture data are used together to estimate effective $K_s$ values, then it is likely that such values will not yield good calibration and validation statistics. This is a common limitation in hydrology, and reaffirms that calibration yields effective values, but these values may not be intrinsic medium properties—this problem remains as a challenge.