Infiltration-runoff model for layered soils considering air resistance and unsteady rainfall

Infiltration and runoff are the dynamic processes of water movement at the soil–air interface during rainfall periods. The accurate performance of the infiltration and runoff is of critical importance in hydrological modeling. As one of the most usable models for practical field problems, the Green–Ampt model can reasonably estimate the depth of the wetting front, infiltration capacity, and cumulative infiltration. It has been subjected to many modifications (Mein & Larson 1973; Jia & Tamai 1997; Chu & Marino 2005) applied in many semi-distributed hydrological models, such as the SWAT (Neitsch et al. 2002a, 2002b) and WEP (Jia et al. 2006).

The Green–Ampt model (Green & Ampt 1911) was originally used as an infiltration model for dry and uniform soil with surface ponding. Due to its simplicity, researchers have developed many modifications, such as Mein & Larson (1973), Ali et al. (2013, 2016), and Vatankhah (2015), etc. All these modified versions have been found to be suitable for simulating the infiltration in homogeneous soil under ponding or steady rainfall conditions. Moreover, it has been extended to model the infiltration in non-uniform soil (Bouwer 1969), during unsteady rainfall periods (Chu 1978), in two-layered soil profiles (Moore & Eigel 1981), and in multi-layered soil profiles with unsteady rainfall (Jia & Tamai 1997).

The soil profile above the wetting front, in the Green–Ampt model, is assumed to be fully saturated, indicating that the air in the soil can escape freely from the soil pores and the air resistance can be abstractly neglected in the water infiltration process (Hammecker et al. 2003). However, many previous studies showed that the air may be trapped and compressed, and the infiltration rate was apparently reduced in the assumed saturated soil profile (Morel-Seytoux & Khanji 1975; Grismer et al. 1994; Wang et al. 1997, 1998; Hammecker et al. 2003). Wang et al. (1998) found that the volume of the residual entrapped air, on the condition of air-confining, increased by 7% on average. The infiltration rate was also decreased between three-fold and ten-fold compared to the air-draining condition. Ma et al. (2011) proposed that the ratios of the actual water content to the saturated water content in the soil profile above the wetting front were 0.80 and 0.86 for the laboratory and field infiltration experiments, respectively. The actual water pressure at the wetting front was determined by the capillary and soil air pressures (Morel-Seytoux 1973). Peck (1965a, 1965b) speculated that the required initiation of the escaping air was expected to be equal to the water pressure at the bottom depth of the saturated zone plus the air entry pressure of the material. The escaping air would cease when the value of pressure infinitely closed to 0 to allow the air escape route to be sealed by effective saturation. At that time, the air pressure would start to increase again with further water penetration. Wang et al. (1997) confirmed Peck's speculation and investigated the relationship of two extreme air pressure relationships in a porous medium, which are referred to as the ‘air-breaking value’ and ‘air-closing value’, respectively. The air pressure was found to be constantly increased with the rising of the wetting front (Wang et al. 1997).

Furthermore, based on the indoor soil column experiments, the changing process of the air pressure in the soil, during water infiltration under air confined conditions, was studied by Grismer et al. (1994), Grismer (2016), Li & Fei (2005), and Li (2007). The results showed that the changing process was divided into three stages: the initial stage, transitional stage, and steady stage. During the initial stage and transitional stage, the soil air pressure automatically changed between the closure and breakthrough pressures. During the steady stage, the soil air pressure was relatively fixed. In addition, certain zones above the wetting front were found to be saturated, while the air pressure under the saturated zone was equal to air pressure of the soil under the wetting front.

To account for the effects of the air resistance, the effective hydraulic conductivity, in place of the saturated value, was applied in some of the Green–Ampt models in many previous studies. Bouwer (1966) suggested that the effective hydraulic conductivity was half of the saturated value. Grismer et al. (1994), Morel-Seytoux & Khanji (1975), and Hammecker et al. (2003) developed a simple flow model based on the Green–Ampt model, taking the air compression and air counter-flow into consideration, while it was only tested to be suitable for a homogeneous soil profile. Ma et al. (2011) developed a simple model for layered soil profiles with constant ponding heads, and introduced a saturated coefficient (S_a) into the model to quantitatively demonstrate the air compression and air counter-flow. Wang et al. (1997) proposed an infiltration model which considered the potential effects of the air compression and air counter-flow during water infiltration in a porous medium, and the air pressure was determined as the sum of the ponding depth, wetting front distance and half of air-bubbling pressure and water-bubbling pressure. Nevertheless, the air pressure tends to infinity as wetting front distance tends to infinity, which is not in compliance with the fact.

In summary, with the development of the basin hydrological models, such as SWAT (Neitsch et al. 2002a, 2002b) and WEP (Jia et al. 2006), the simulations of the soil infiltration runoff process gradually took more factors into consideration in the hydrological models, including soil stratification characteristics, air resistance, and unstable rainfall process, to improve the accuracy of the models and solve problems such as the simulation distortions of the models. However, certain deficits were found in the models currently applied and the applicability of the models was limited, and these are systematically summarized in this paper. Based on these deficits, a stratified soil infiltration runoff model under unsteady rainfall, considering the effects of air resistance, was proposed and proven to be practical via indoor experimental testing.

The objectives of this study were as follows: (1) to develop a modified Green–Ampt model to simulate water infiltration in layered soil profiles under unsteady rainfall conditions, considering the influence of air resistance; (2) to test the applicability of the proposed model with one-dimensional infiltration-runoff experiments with different combinations of soil layered patterns during rainfall events in multi-layered soil columns.

In Table 1, the summaries of the Green–Ampt model and other primary modified models are shown, indicating that the previous models have many deficiencies and need to be further improved. In this study, to make up for the deficiency of the above model research, based on the indoor experimental research, the new stratified soil unsteady rainfall, the Green–Ampt infiltration runoff model considering air resistance effect (MGAM) is proposed.

In Figure 1, the schematic diagram of infiltration into a layered soil profile is illustrated. The soil profile under consideration consists of n layers with thickness Z_i (i= 1, 2, …, n), and an initial water content θ_0,i. The saturated hydraulic conductivity and the saturated water content of each layer of the soil are k_s,i and θ_s,i, respectively. Z_i represents the thickness of each layer. The origin of the coordinate system is set at the surface of soil, and the coordinate system is positively downward.

In view of the deficiency of the above models, the improvement of the modified Green–Ampt model (MGAM) is illustrated in Figure 2. First, GGAM was originally used as an infiltration model for the uniform soil. In this case, it extends to model infiltration in the layered soil profiles. As performed in the traditional Green–Ampt model (TGAM), the water content of each soil layer below or above the wetting front is assumed to be uniform, and the soil profile is separated into an upper wetted zone and a lower unwetted zone by wetting front. In the unwetted zone, the soil water content is assumed to be the initial value. H₀ and H_L are the hydraulic heads at the surface and wetting front of the soil, respectively. Z represents the depth of wetting front. Second, as a result of the air entrapment, the soil pores in the water content profile of the soil above the wetting front cannot be completely filled with water (Bouwer 1966; Hammecker et al. 2003). Therefore, it is necessary to take the actual water content (θ_a,i) and the actual hydraulic conductivity (k_a,i) of the ith layer in the soil water content profile above the wetting front as the water content and hydraulic conductivity at the residual air saturation (Hammecker et al. 2003). Previously, θ_a,i and k_a,i, which are less than θ_s,i and k_s,i, have been used in the TGAM. As in Bouwer Green–Ampt model (BGAM) and GGAM, a saturated coefficient S_aθ,i (0 ≤ S_aθ,i ≤ 1) is also introduced to determine the proportion between θ_a,i and θ_s,i in the ith soil layer (θ_a,i = S_aθ,i · θ_s,i). In regard to the relationship between k_a,i and k_s,i, S_ak,i (k_a,i = S_ak,i · k_s,i), representing the relative water conductivity, is adopted. However, in this case, the parameters' calculation methods are improved to be suitable for layered soil profiles.

In addition, due to the effects of the air and capillary pressures at the wetting front, the actual water pressure SW_wf,i at the wetting front is given by SW_wf,i = SW_cf,i − SW_af,i (Morel-Seytoux 1973), where SW_af,i represents the instantaneous air pressure below the wetting front, and SW_cf,i represents the average capillary pressure at the wetting front.

In order for an air bubble to escape from a ponded soil surface, the entrapped air phase needs a sufficiently high air pressure to overcome the static water pressure (Wang et al. 1998). During the infiltration process, the soil water content profile above the wetting front can be divided into four zones (Bodman & Colman 1944, 1945; Shao et al. 2006), namely, the saturated, transitional, transmitting, and wetting zones (Figure 3). In the saturated zone, the soil pores are completely filled with water, and no pores are interconnected except when the air eruption occurs. Therefore, the air in the soil cannot move freely, and SW_af = 0. The remaining zones cannot be completely filled with water, and the macro pores are interconnected. Hence, the air in the soil can move freely, and SW_af >0 (Grismer et al. 1994; Li 2007).

Wang et al. (1997) suggested that the h_ab and h_wb can be evaluated by and (⁠ is equal to 0 to approximately 2 cm for sandy soil; 2 to 5 cm for loamy soil; and 8 to 10 cm for clay soil), where is a parameter of the main drainage curve.

In the research conducted by Whisler & Bouwer (1970), it was proposed that SW_cf, the average capillary pressure at the wetting front, could be replaced by the air-bubbling capillary pressure h_ab in the model presented by Brooks & Corey (1964, 1966).

In regard to the infiltration during unsteady rainfall events, the rainfall consists of a series of time intervals of equal length, within which the rainfall rate remains constant. The infiltration process at time (t_x-1 ∼ t_x) can be determined simultaneously by the ponding depth h₀ of the initial time (ponding water depth in the early stages of the period), rainfall rate P, and ponding infiltration rate f_pt, where t_x is the time at x (⁠⁠), and n represents the number of times the discrete time periods occurred in the rainfall record. Since the ponding water depth h₀ has a different unit from the rainfall intensity P, it is necessary to first divide the ponding water depth in the early stage of this period by time during the calculation process (⁠⁠), and then substitute this into the model for the calculation.

As illustrated in Figure 4, according to the ponding depth h₀ in the early stage of the period, along with the rainfall intensity P’, rainfall rate P, and ponding infiltration rate f_pt, the infiltration process can be divided as follows.

Case a: h₀ = 0, P > f_pt ≥ ⁠, where is the average actual hydraulic conductivity in the transmission zone behind the wetting front. In this case, ponding will occur, and the infiltration process can be divided into non-ponding infiltration f_npt (before t_p) and ponding infiltration f_pt (after t_p), where f_npt is the non-ponding infiltration rate and t_p is the time to surface ponding.

Case b: h₀ > 0, P’ + P<≤f_pt. In this case, no ponding will occur, and the infiltration process can be divided into ponding infiltration f_pt (before t_p) and non-ponding infiltration f_npt (after t_p).

Case c: h₀> 0, P’ + P ≥ f_pt ≥ ⁠. In this case, the infiltration process is under the ponding condition f_pt.

Case d: h₀ = 0, P < ≤f_pt. In this case, the infiltration process is under the non-ponding condition f_npt.

During the non-ponding infiltration process under rainfall conditions, all of the rainfall infiltrates into the soil layers, and no surface runoff occurs. During the ponding infiltration process, the surface runoff can be calculated by the mass balance equation (Equation (23)). The infiltration can then be calculated by Equations (12) and (13) before the wetting front enters the second soil layer (Z≤), and by using Equations (30) and (31) after it enters the mth layer.

To the authors' knowledge, infiltration-runoff experiments regarding three-layered soil under unsteady rainfall conditions have not yet been reported. The combinations of the layered patterns of the different soil types are known to have great effects on the infiltration process when the wetting front enters the interface between the different soil types. Therefore, further investigations via infiltration experiments are required. In addition, during the process of rainfall infiltration, the air in the soil is difficult to discharge, as it is sealed within the soil body. The air pressure produced by the enclosed air hinders the soil infiltration process. The influence of the air resistance on the rainfall infiltration-runoff process also requires further experimental researches.

Figure 7 shows the schematic of this study's experimental setup. The experiments were performed on a uniformly packed column. The length and inner diameter of this column were 50 cm and 29 cm, respectively. The soil column bottom was sealed to prevent the air in the soil from emitting from the bottom. A Marriott tube (with an inner diameter of 10 cm, and height of 50 cm) was used to supply water, and the rainfall intensity was controlled by a medical needle with a No. 5 pinhead simulating 0.10 to 0.32 mm/min; a No. 7 pinhead simulating 0.32 to 1.9 mm/min; and a No.12 pinhead simulating larger than 1.9 mm/min. The rainfall simulator was kept at a certain distance (approximately 2.0 m from the surface of the soil).

In this study, three types of air drying soil, which had been obtained from China's Hunan Province, were contained in the aforementioned column, and formed the clay, loam, and sand layers. In order to obtain sand, clay, and loam, three sampling fields were selected, and three sampling sites were set up. The soil samples (three replications for each type of soil) were collected with a soil wreath knife (100 cm³) at sampling sites to determine soil dry bulk density.

The loading patterns are shown in Table 2. We chose these three-layer soil loading patterns rather than other kinds of soil loading patterns based on two considerations: (1) to save the experimental work task because there is a total of six loading patterns for the three kinds of soil; (2) the three soil loading patterns in Table 2 are believed to be more frequently encountered in the field than other kinds, for example, the sand/loam/clay loading pattern and the sand/loam/clay loading pattern in farmlands because of land reclamation and tillage activities.

The air-dried soil was filled in layers at a certain bulk density (the same as sampling sites), with a number of small increments in the vertical cylinder. Each increment was set at 5 cm, and each soil layer was 15 cm thick. The dry soil was divided into layers according to volume weight (5 cm/layer), as shown in Table 2, and placed into the soil column. During the loading process, in order to avoid artificial effects on soil layer interface, when the upper soil samples were being loaded, the lower compacted soil samples were first loosened, and then filled by the upper soil samples. Then, in accordance with the test settings, the filling thickness of each type of soil was 15 cm, and the total filling thickness 45 cm. During the compacting process, five TDR probes were installed into the column to measure the water content of the soil. Considering the soil layered characteristic and soil water content change with increase of soil depths, the TDR probes were installed at the depths of 5 cm, 15 cm, 25 cm, 35 cm, and 45 cm. The infiltration in the layered soil experiments was designed with three replications of three treatments, and each treatment included a three-layered soil profile.

The particle-size distribution of the soil was measured by the sieve-pipette method. The saturated soil water content was determined by a soak test (Klute 1986). Then, the tested samples were prepared using the same procedure relative to packing, and had the same initial water content and bulk density as those in the column. The saturated soil conductivity was determined by a constant-head method, as illustrated in Figure 8, with a head difference of 10 cm (Huo et al. 2010). All of the soil samples were slowly saturated from the bottom for approximately 1 day, in order to ensure that they became fully saturated before the saturated water content and hydraulic conductivity were measured. The initial water content and bulk density of the soil were determined by an oven drying method with three replications. The field moisture capacity was determined by the Colman method (Colman 1947). The physical properties of the soil are shown in Table 3. A pressure plate method was used to determine the water characteristic curve (US-made 1500-type pressure membrane instrument) of the soil, and the van Genuchten (1980) model was adopted for fitting (Figure 9). The fitting parameters are shown in Table 4.

The infiltration-runoff experiments in this study were conducted with a rainfall process. During the experiment, the water table was a Mariotte bottle, and the wetting front distance and runoff were measured to calculate the cumulative infiltration, infiltration rate, and wetting front distance movement. The experiment was terminated when the wetting front reached the bottom of the last soil layer. Then, after the infiltration experiment was completed, the soil cores at different depths were sampled to measure the water content using an oven-dry method. The temperature remained at 25 ± 1 °C and the evaporation relevant to the infiltration was minimal, and thus neglected.

The successful performance of the MGAM is dependent on the determination of the appropriate values of actual hydraulic conductivity (k_a,i = S_ak,i · k_s,i), actual soil water content (θ_a,i = S_aθ,i · θ_s,i), average capillary pressure at the wetting front (SW_cf), air-bubbling capillary pressure (h_ab), and water-bubbling pressure (h_wb).

The values of θ_ra and θ_rw are difficult to determine, and complex laboratory infiltration experiments and precise setups are required, such as the use of a tension-pressure infiltrometer and an air flowmeter (Wang et al. 1998). Chong et al. (1981) proposed that the saturation coefficient of a soil layer could be readily calculated with the measured water content in the soil water profile above the wetting front at the termination of the infiltration experiment (θ_a) by S_aθ = θ_a/θ_s. However, it was found to be more effective to develop an independent approach which only required the physical properties and hydraulic parameters of the soil to determine the saturation coefficient.

During the infiltration process, the soil water content profile above the wetting front can be divided into four zones (Bodman & Colman 1944, 1945; Shao et al. 2006), including the saturated, transitional, transmitting, and wetting zones. The water content of the soil in the saturated zone was found to be equal to the saturated water content. In the transitional zone, the water content was determined to be close to the saturated water content, and the transmitting and wetting zones were close to the field moisture capacity (Gan et al. 2012). The water content of the soil above the wetting front was also divided into a saturated zone and a transitional zone. The saturated zone was approximately half the area of the soil water profile above the wetting front (Wang et al. 2002). In this study, for the sake of simplicity, the soil water content profile above the wetting front was also divided into saturated and transitional zones, and the saturated zone was half of the soil water content profile above the wetting front: θ_a = θ_s (S_aθ = 1) in the saturated zone, and θ_a = θ_h (S_aθ = θ_h/θ_s) in the transitional zone, where θ_s is the saturated water content and θ_h is the field moisture capacity.

The air pressure (SW_af) at the wetting front varied dynamically from the air-closing value (H_c) to the air-breaking value (H_b). Also, the corresponding value of the average capillary pressure (SW_cf) at the wetting front automatically varied following a scanning drainage curve towards a scanning wetting curve. Whisler & Bouwer (1970) proposed that the average capillary pressure SW_cf could be replaced by the air-bubbling capillary pressure in the wetting retention model of Brooks & Corey (1964, 1966). For example, SW_cf = h_ab = 1/α_d, where the air pressure SW_af = h₀ + 0.5(h_ab + h_wb), and h₀ is the ponding depth. The value of h_ab = 1/α_d, where α_d is the parameter of the main drainage curve and h_wb = h_ab/2-δ (δ equals approximately 0 to 2 cm for sandy soil; 2 to 5 cm for loamy soil, and 8 to 10 cm for clay soil).

The MLGAM was found to be suitable for simulating the infiltration-runoff process in a homogeneous soil with uniform initial moisture content. However, when the MLGAM model was applied for the simulation of the infiltration-runoff process in a layered soil profile, it was found that the uniform soil-water movement parameters were difficult to obtain. Therefore, in order to obtain the uniform soil-water movement parameters of a layered soil profile, a harmonic mean was applied to calculate the average soil saturated hydraulic conductivity along with the average suction head of the wetting front. Then, a mathematical mean was applied to calculate the saturated water and soil initial water contents of the soil. The effective hydraulic conductivity was then calculated with Equation (38). The total rainfall P and runoff volume R were the measured values. The values of P and R were different between different loading patterns. Thus, CN and k_e have different values for different loading patterns.

In all of the models (with the exception of the MLGAM), the soil saturated water content θ_s, field moisture capacity θ_h, and saturated hydraulic conductivity k_s were measured. Also, the relative water conductivity (S_ak), saturated coefficient (S_aθ), soil suction head SW_cf  at the wetting front, air-bubbling capillary pressure h_ab, and water-bubbling pressure h_wb were determined by the soil properties. The models did not conduct any parameter adjustments, with the exception of the MLGAM. The soil hydraulic parameters, and the model parameters for the laboratory experiments, are summarized in Tables 5–7, respectively.

The successful performance of the MGAM was dependent on determining the appropriate input parameters. For the MGAM, the following factors needed to be considered: the determination of which parameter was dominant; the reduction of the calibration uncertainty; and the enhancement of the optimization efficiency. However, the means by which to achieve these still remained unclear. Therefore, the MGAM required further examination by a parameter sensitivity analysis. An LH-OAT parameter sensitivity analysis was adopted in the MGAM model. The LH-OAT is a type of new method which integrates the LH (Latin-Hypercube) sampling, and the OAT (one-factor-at-a-time) sensitivity analysis. Therefore, it has the advantages of both the LH sampling and the OAT sensitivity analysis.

A final effect (sensitivity degree SA) was calculated by averaging the partial effects of each loop for all of the Latin Hypercube points. The final effects could then be ranked with the calculated sensitivity degree, as shown in Table 8.

In order to analyze the sensitivity of the input parameters of the MGAM model, it was applied to the infiltrations in a three-layered soil profile during rainfall conditions. Then, 12 key parameters of the MGAM were selected for the sensitivity analysis. The observed cumulative infiltration and runoff rate were selected as the targeted output variables. Table 9 lists the name, range, and description of each input parameter which was investigated in this study. The MGAM model parameters were the relative water conductivity (S_ak,1, S_ak,2, S_ak,3), air-bubbling pressure values (h_ab1, h_ab2, h_ab3), parameters (δ₁, δ₂, δ₃), and initial water content deficit (Δθ₁, Δθ₂, Δθ₃).

The sampling process divided the entire parameter space into 20 layers, and 12 parameter sets were obtained through the sampling. An OAT sampling of each parameter set was made in each parameter space LH sampling. Since each OAT sampling only changed one parameter value, and there were 12 parameters in each group, then 12 samplings were made of all the OAT samplings of the LH parameter set in the group. Therefore, the LH-OAT sampling number was 260 = 20*(1 + 12). Also, the increment of each input parameter was set at 5%.

Since the MGAM model is suitable to simulate the infiltration-runoff process during unsteady rainfall, the parameters' values change according to soil properties, initial water content, etc. Therefore, the default values of these parameters are changed in a relative way over a certain range. The soil types can be divided into three major types of soil (clay, loam, and sand) according to soil particle-size distribution, thus the default parameters' values change among three major soil types of soil physical properties, for instance over the range 5.0 cm (sand) ∼130 cm (clay) for air-bubbling pressure.

S_ak is a dimensionless coefficient that reflects the effluence of saturated coefficient (S_aθ) on soil hydraulic conductivity. Moreover, the saturated coefficient reflects the saturation degree of the soil water profile above the wetting front. h_ab is the air-bubbling capillary pressure value of the material, cm; δ is a parameter for calculating the water bubbling pressure (h_wb), which equals approximately 0 to 2 cm for sandy soil, 2 to 5 cm for loamy soil, and 8 to 10 cm for clay soil (Wang et al. 1997); Δθ is the initial water content deficit, cm³/cm³.

The results of the sensitivity analysis for the infiltration-runoff process are shown in Table 10. It can be seen that the parameter sensitivity displayed clear differences between the different layers. For the soil's cumulative infiltration and runoff rate, the sensitivity of relative water conductivity values in the first layer were found to be very high (3.42 and 3.18, respectively). Those of the air-bubbling pressure in the first and second layer were high, and all the values ranged from 0.26 to 0.40. The parameters δ in the first and second layers were determined to be medium, and all the values ranged from 0.09 to 0.18. The initial moisture deficits in the first and second layers were found to be high or medium, and all of the values ranged from 0.13 to 0.28. The sensitivities of the other parameters were determined to be small for the infiltration-runoff process, and all sensitivity values were less than 0.05.

For the cumulative infiltration and runoff rate, the parameter sensitivities were also reflected in the class index, and the sensitivity values of relative water conductivity in the first layer were in class IV, and for the second and third layers, the values were in class I. The sensitivity values of the air-bubbling capillary pressure of the first and second layers were in class III, and were in class I for the third layer. The sensitivity values of the parameters δ in the first and second layers were in class II, and in class I for the third layer. The sensitivity values of the initial moisture deficit in the first and second layers were in class III or II, and for the third layer were in class I.

In summary, the parameter sensitivity analysis showed that the sensitivity of relative water conductivity was very high in the first layer, and those of the air bubbling pressure and initial moisture deficit were high or medium in the first and second layers, whereas those of the other parameters were found to be small.

The comparison of the measured and simulated soil moisture distribution, by the MGAM and BGAM (Bouwer 1966); TGAM (Jia & Tamai 1997), and MLGAM (Mein & Larson 1973), at the times when the wetting front entered 15, 30, and 45 cm, are shown in Figure 10. With the exception of the results from the MLGAM, the observed and the simulated soil water distribution were determined by the combination of the soil layered patterns.

The TGAM significantly overestimated the soil water content, and the BGAM underestimated those values for the fine over coarse layered patterns. The MGAM results were in better agreement with those of the measured values. The simulated values by the MLGAM were determined to be larger than those of the observed values in transitional zones, less than those of the observed values in saturated zones for the loading pattern of the clay-loam-sand, and larger than those of the observed values for loading pattern of the sand-clay-loam. Therefore, when compared to the TGAM, BGAM, and MLGAM, the simulated results of the MGAM were found to be in better agreement with the measured soil water content.

The performances of the TGAM, BGAM, MLGAM, and MGAM in the simulations of the soil water content were also reflected in their associated MAE, ARE, and RMSE values, as shown in Table 11. In regard to the loading patterns of the clay-loam-sand, loam-sand-clay, and sand-clay-loam, most of the MAE values for the BGAM ranged from approximately 0.003 to 0.067 cm³/cm³; for the TGAM, they ranged from approximately 0.026 to 0.043 cm³/cm³; and for the MLGAM, they ranged from approximately 0.065 to 0.073 cm³/cm³, whereas for the MGAM the range was from approximately 0.002 to 0.006 cm³/cm³. In regards to loading patterns of the clay-loam-sand, loam-sand-clay, and sand-clay-loam, it was found that most of the ARE values for the BGAM ranged from approximately −9.72% to 0.46%; for the TGAM, the range was approximately from 3.95% to 10.56%; and for the MLGAM, the range was from approximately 32.29% to 72.95%, whereas for the MGAM the range was determined to be from approximately −0.46% to 3.50%. In regards to the loading patterns of the clay-loam-sand, loam-sand-clay, and sand-clay-loam, the NSE values for the BGAM were found to range from approximately −3.06 to 0.75; for the TGAM, the range was from −5.11 to −0.82; and for the MLGAM, the range was found to be from −21.25 to −8.44, whereas for the MGAM the range was determined to be from approximately 0.75 to 0.89.

The overestimation of the soil water content provided by TGAM assumed that the soil water content profile above the wetting front was fully saturated. However, the soil water content profile above the actual wetting front could not be saturated, due to the entrapment of air during the infiltration process, and the actual water content in the soil water content profile above the wetting front was less than the saturated water content. The underestimation of the soil water content by the BGAM for the fine over coarse layered soil was due to the assumption that hydraulic conductivity of the soil water content profile above the wetting front was only half of the saturated hydraulic conductivity. However, the actual hydraulic conductivity in the first layer above the wetting front was larger than the half of the saturated hydraulic conductivity. Therefore, the BGAM underestimated the soil cumulative infiltration, and provided a smaller depth of wetting front than that of the measured results. The poor performance of the MLGAM for all of the loading patterns was due to the fact that the MLGAM was suitable for simulating the infiltration-runoff process in homogeneous soil, but not in the layered soil. In addition, to all treatments, the values of saturated water content in MLGAM are obtained by mathematical mean method, so all treatments share the same values above the wetting front.

The comparison of the measured and calculated results of the soil cumulative infiltration by the TGAM, BGAM, MGAM, and MLGAM are shown in Figure 11. With the exception of the results from the MLGAM, the observed and simulated soil cumulative infiltrations were determined by the combination of the soil layered patterns.

Also, with the exception of the MLGAM, in the loading patterns of the clay-loam-sand, loam-sand-clay, and sand-clay-loam, the calculated and measured values were found to undergo sudden changes when the wetting front entered 15 cm and 30 cm. The reasons for these changes may have been that the capillary suction and hydraulic conductivity underwent major changes when the wetting front entered the interfaces between the layers.

The simulated soil cumulative infiltration by the TGAM was found to be larger than the measured results, and the differences between them were determined to be related to a combination of soil layered patterns. In regard to the loading patterns of the clay-loam-sand, loam-sand-clay, and sand-clay-loam, the measured total cumulative infiltration values at the termination of the experiments were 14.62, 16.50, and 15.73 cm, respectively. The simulated values by the MGAM were 15.12, 16.6, and 15.53 cm, respectively. However, the simulated values by the BGAM were 10.26, 15.72, and 16.33 cm, respectively; 15.05, 20.28, and 25.21 cm, respectively, by the TGAM; and 58.21, 26.66, and 49.66 cm, respectively, by the MLGAM. Therefore, when compared with the BGAM, TGAM, and MLGAM, the simulated values for the MGAM were found to be in better agreement with the observed values.

The performances of the TGAM, BGAM, MLGAM, and MGAM in the simulations of the cumulative infiltration were also reflected in their associated MAE, ARE, RSME, and NSE values, as shown in Table 12. The majority of the MAE for the BGAM ranged from 0.39 to 3.39 cm; TGAM ranged from 0.64 to 1.37 cm; and the MLGAM ranged from 10.83 to 52.85 cm, whereas the MGAM was less than 0.40 cm. The majority of the ARE for the BGAM were in the range of −48.0% to 10.4%; in the TGAM they were larger than 25.0%; and in the MLGAM they were larger than 25.63%, whereas for the MGAM the range was from −7.9% to 6.8%. The majority of the RSME for the BGAM were in the range of 0.46 to 3.78 cm; for the TGAM, the range was from 0.73 to 2.84 cm; and for the MLGAM, the RSME was larger than 12.83 cm, whereas the RSME for the MGAM was less 0.60 cm. Also, the majority of the NSE values for the BGAM were 0.69, 0.90, and 0.91, respectively; for the TGAM, they were 0.98, 0.91, and 0.70; and for the MLGAM, the values were −16.02, −11.22, and −0.63. However, the NSE values for the MGAM were all determined to be larger than 0.98.

The overestimation of the soil cumulative infiltration by the TGAM was due to the fact that it ignored the entrapment of air in the upper soil water profile above the wetting front, and also the air compression below the wetting front. In reality, the upper soil water profile above the wetting front was unsaturated. The hydraulic conductivity of the soil water profile above the wetting front was found to be smaller than that of the saturated zone. Therefore, the TGAM overestimated the cumulative infiltration. In regard to the loading pattern of the clay-loam-sand and loam-sand-clay, the underestimation of the cumulative infiltration by the BGAM was due to the fact that it underestimated the actual hydraulic conductivity of the first layer in the soil water profile above the wetting front. For the loading pattern of the sand-clay-loam, the overestimation of the cumulative infiltration by the BGAM was due to the fact that it overestimated the actual hydraulic conductivity of the first layer in the soil water profile above the wetting front. The overestimation of the cumulative infiltration by the MLGAM was due to the fact that all of the MLGAM parameters experienced no changes during the infiltration process. However, when the wetting front entered the interfaces between two soil layers, the changes of the soil's physics resulted in changes in the actual hydraulic conductivity, capillary pressure head, actual water content, and initial water content of the soil, which caused changes in the infiltration process.

The comparison of the measured and calculated results of the runoff rate are shown in Figure 12. The numbers I, II, and III in Figure 12 represent the repetitions of each treatment. With the exception of the results of the MLGAM simulation, the observed and simulated runoff characteristics were determined by a combination of the soil layered patterns.

The mean calculated runoff-yielding time was also found to be related to a combination of the soil layered patterns. In regard to the loading pattern of the clay-loam-sand, the mean simulation runoff-yielding times by the MGAM (10 minutes), BGAM (4 minutes), and TGAM (10 minutes) were less than that of the observed times (29 minutes), and the mean time of the MLGAM (60 minutes) was found to be markedly larger than that of the observed results. For the loading pattern of the loam-sand-clay, the mean simulation runoff-yielding times by the MGAM (26 minutes) and the TGAM (45 minutes) were determined to be larger than that of the observed results (25 minutes). Also, the BGAM (6 minutes) and MLGAM (6 minutes) were found to be less than that of the observed results. In regards to the loading pattern of the sand-clay-loam, the mean simulated runoff-yielding times by the MGAM (5.0 minutes), BGAM (5.0 minutes), and MLGAM (0 minutes) were determined to be less than that of the observed results (6.0 minutes). However, those of the TGAM (11 minutes) were found to be markedly larger than that of the observed results.

The runoff process is known to fluctuate with the fluctuations in rainfall. The simulated values of the TGAM were found to be less than the observed values. The simulated values of the BGAM were larger than those of the observed values for the loading pattern of the clay-loam-sand and loam-sand-clay, and less than those of the observed values for the loading pattern of the sand-clay-loam. The simulated values by the MLGAM were determined to be less than that of the observed values for the loading pattern of the clay-loam-sand, larger than that of the observed values before the wetting front entered the third layer, and less than that of the observed values after the wetting front entered the third layer for the loading pattern of the loam-sand-clay. They were larger than that of the observed values before the wetting front entered the second layer, and less than that of the observed values after the wetting front entered the second layer for the loading pattern of the sand-clay-loam. Therefore, when compared to the TGAM, BGAM, and MLGAM, the simulated results of the MGAM were found to be in better agreement with the measured runoff rate. The performances of the TGAM, BGAM, MLGAM, and MGAM in the simulations of the runoff rate were reflected in their MAE, ARE, RSME, and NSE values, as shown in Table 13. The majority of the MAE for the BGAM were found to be in the range of 0.006 to 0.083 cm/min; the TGAM range was from 0.007 to 0.048 cm/min; and the MLGAM ranged from 0.036 to 0.152 cm/min, whereas the MAE for the MGAM was less than 0.021 cm/min. The majority of the ARE for the BGAM range from −25.5% to 34.3%; the TGAM were found to be less than −28.12%; and that of the MLGAM ranged from −14.00% to 107.23%, whereas the range for the MGAM was from −10.1% to 9.4%. The majority of the RSME values for the BGAM were in the range of 0.009 to 0.24 cm/min; the TGAM were larger than 0.014 cm/min; and the MLGAM were found to be larger than 0.04 cm/min, whereas the MGAM were less than 0.028 cm/min. Also, the majority of the NSE values for the MGAM were found to be larger than 0.75; the BGAM were less than 0.77; the TGAM less than 0.73; and the MLGAM were less than −0.54.

The underestimation of the runoff rate by the TGAM was due to the fact that it ignored the air entrapment and air compression. The TGAM assumed that the soil water content profile above the wetting front was fully saturated with water (Green & Ampt 1911). The water content and hydraulic conductivity of the soil water content profile above the wetting front were referred to as the saturated soil water content, and hydraulic conductivity, respectively. In reality, during the rainfall-infiltration-runoff process, the air in the soil was displaced and compressed ahead of the wetting front. The compression of the air in the soil led to a substantial decrease in the infiltration rate. In addition, having been affected by the entrapped air, the water content profile of the soil above the wetting front was unsaturated. The hydraulic conductivity of the soil water content profile above the wetting front was smaller than the saturated one. Therefore, the TGAM overestimated the infiltration rate, and underestimated the runoff rate.

The overestimation of the BGAM for the loading pattern of the clay-loam-sand and loam-sand-clay was due to the fact that the hydraulic conductivity of the soil water content profile above the wetting front was assumed to be only half of the saturated hydraulic conductivity. In reality, the actual hydraulic conductivity of the first layer in soil water content profile above the wetting front was larger than half of the saturated hydraulic conductivity. In regard to the loading pattern of the sand-clay-loam, the BGAM overestimated the actual hydraulic conductivity of the first layer in the soil water content profile above the wetting front. Therefore, for the loading pattern of the clay-loam-sand and loam-sand-clay, the BGAM underestimated the infiltration rate, and overestimated the runoff rate, and vice versa for the loading pattern of the sand-clay-loam.

In this study, the MLGAM was found to be suitable for the simulations of homogeneous soil with uniform initial moisture content under steady rainfall. For all of the loading patterns, the underestimations or overestimations of runoff rates by the MLGAM were due to the fact that all of the MLGAM parameters experienced no changes during the infiltration process. In fact, when the wetting front entered the interface between two soil layers, the changes in the soil physics resulted in changes in the infiltration and runoff rates, due to the different water movement parameters.

This study presented a modified Green–Ampt model, and then evaluated the model's performance by the laboratory observed soil water content, runoff rate, and soil cumulative infiltration under different soil loading patterns (different combinations of sand, loam, and clay) during rainfall events. In particular, the model was compared with the BGAM (Bouwer 1966), TGAM (Jia & Tamai 1997), and MLGAM (Mein & Larson 1973). The influence of the mechanism of the soil properties, and the order of the soil layers during the infiltration-runoff process, were very complicated. MGAM was able to achieve a better match with the observation than the TGAM, BGAM, and MLGAM, by considering the air resistance in the layered soil. In addition, all of the parameters in MGAM had physical meaning, and could be determined from the physical and hydraulic properties of the soil. This indicated that the modified Green–Ampt model had the ability to effectively depict the soil type and water content of the layered soil infiltration-runoff process during unsteady rainfall events.

Although MLGAM (Mein & Larson 1973), GGAM (Grismer 2016) did not consider the soil stratification, the model parameters' correction based on field rainfall simulation results may compensate for the lack of consideration of stratification. These results may be mainly on soils with high water permeability. The soil stratification may not be obvious in these watersheds (such as Tahoe Lake soil, volcanic ash soil, etc.), so their parameters' modification results are still relatively good. However, for clearly stratified soils (including natural geological stratification, and human-affected strata, such as farmland and urban soils, etc.), it is very significant to consider soils' stratification.

MGAM considered many influencing factors, such as the soil layers, air resistance, unsteady rainfall, and so on, which potentially could have great effects on the infiltration-runoff process in a field or real catchment scale, to easily eliminate interference occurring in a field or a real catchment scale. These may include effects as such unbalanced amounts of water, large observational errors, and so on. Therefore, in order to better validate the MGAM, corresponding laboratory experiments were designed, and then the model was evaluated by the laboratory data.

When the MGAM model is applied to real catchment scale, in order to obtain model parameters, based on the principle of geostatistics, the soil types and stratified features in the real catchment scale were investigated, stratified sampling was carried out, and then brought back into the room to measure the water characteristic parameters, including the water characteristic curve, saturated hydraulic conductivity, and field water capacity, etc.

This study was sponsored by the following projects: the National Basic Research Program of China (973 Program) (2015CB452701); the Projects of the National Science Foundation of China (51779272); the National Key Research and Development Program of China (2017YFC0404603); the Natural Science Foundation of Qinghai Province of China (2018-ZJ-965Q); Autonomous Research of the State Key Laboratory of Simulations and Regulations of Water Cycles in River Basins (SKL-WAC) (2016ZY03). The authors would like to express their gratitude to the reviewers of this study. Their suggestions contributed significantly to the improvement of the manuscript.