Abstract

To simulate the soil water movement process of deformable soils, a modified Richards model considering soil deformation (RMSD) was established. In the model, new parameters were introduced, including deformable soil porosity, deformable soil saturated hydraulic conductivity and unsaturated hydraulic conductivity of expansive soils, which varied with soil depth and time under the effect of soil deformation. The newly introduced parameters originated from physical properties of the soils and their calculation formulas were suggested. One-dimensional infiltration–runoff experiments were performed to evaluate the performance of the RMSD. The results showed that average relative errors (ARE) of simulated runoff intensity and cumulative infiltration (by the RMSD) ranged from −10.0% to −1.0% and from −1.0% to 11.0%, respectively, and Nash efficiency coefficients (NSE) of simulated cumulative infiltration (by the RMSD) were larger than 0.90. As the RMSD model is much better than the traditional Richards model (TRID) in fitting the observations of soil cumulative infiltration and runoff intensity, it is believed that the newly suggested model provides a suitable tool to depict the soil water movement in deformable soils.

INTRODUCTION

Soil deformation refers to the phenomenon of strain and displacement of soil under external forces. Deformable soil will swell, shrink or collapse along with the movement of soil water, which will greatly affect hydrological, chemical and biological processes, soil and water conservation, agricultural irrigation and engineering construction (Wang et al. 2006; Lu & Pei 2007; Vereecken et al. 2016). Deformable soil is widely distributed all over the world (Smiles 2000), especially in the Loess Plateau.

The presence of soil deformation introduces difficulties for both description of soil water movement processes and parameter measurements of the deformable soils. Although many hydrological models have been developed to compute the soil water movement in non-deformable soil, the establishment of suitable models for soil deformation is still challenging (Strudley et al. 2008; Alaoui et al. 2011; Vereecken et al. 2016). For instance, Kostiakov and Horton models, the semi-empirical and empirical models, cannot provide detailed information of soil water movement process due to simple equations lacking in physical meaning (Lei et al. 1988; Mishra et al. 2003). The Green–Ampt model and its modified versions have simple formulation and the model parameters can be easily obtained from soil physical properties (Loáiciga & Huang 2007; Ying et al. 2010; Gan et al. 2015). However, the Green–Ampt models are limited in simulating soil infiltration process, soil water redistribution process, and water absorption of the root. As a physically based numerical model, the Richards equation (Richards 1931; Zhang et al. 2003; Li et al. 2009; Wang & Yang 2013) has been extended into many complex conditions (Richards et al. 1997; Dam & Feddes 2000; Peth et al. 2006; Barontini et al. 2007). Consequently, the Richards equation might be the best fitted model to describe the soil and water movement process of soil deformation. However, the Richards equation is strongly non-linear and cannot easily be solved analytically, especially under complex soil deformation conditions. Thus, soil deformation needs to be explicitly and simply interpreted by numerical equations.

Currently, a wealth of analytical and numerical solutions of flow and transport equations of soil deformation have been conducted. The amount of soil deformation depends on applied load (e.g., static and dynamic loads) (Wiermann et al. 1999; Ghezzehei & Or 2001), soil type and moisture status (Lindstrom et al. 1995). However, these studies focus on effects of soil compaction or swelling on soil hydraulic conductivity and properties, lacking systematic study on effects of soil deformation on runoff and infiltration. Garnier et al. (1997) simulated the hydraulic conductivity of deformed soil based on the Eulerian description (Ruy et al. 1999) and the Lagrangian description (Raats & Klute 1968). McGarry & Malafant (1987) proposed the three-line model to describe the relationship between soil water content and soil deformation, and the relationship between soil stress and strain was described by logarithmic function (Zhou et al. 2010; Cui & Miao 2011). Giráldez et al. (1983) used the universal soil volume change equation to simulate the soil shrinkage characteristic curve. Smiles (1997) deduced the equation of mass-energy conservation based on the Darcy law and continuity equation. Horn & Fleige (2003) developed pedotransfer functions to estimate the compaction sensitivity. Assouline (2006a, 2006b) extended models for soil water retention and hydraulic conductivity curves to calculate and predict the hydraulic properties of compacted or tilled soils. However, these studies still had a tedious and limited quantitative understanding and evaluation of deformable soil structure and their effects on hydrological, chemical and biological processes (Logsdon et al. 2013), and most models had limited applications due to many abstract input parameters. In addition, some researchers have studied the effect of soil property on soil deformation of loess in the Loess Plateau, including soil clay content, organic matter content, bulk density, calcium content and dry–wet cycling (Wang et al. 2005; Zhang et al. 2006), while systematic models on the effects of soil deformation on runoff and infiltration were not conducted. Therefore, the effect of soil deformation on runoff and infiltration is still in the exploratory stage, and the related models are not yet mature or not well established.

In this study, a modified deformable soil water movement model, the RMSD (Richards model considering soil deformation) was proposed to simulate the effects of soil deformation on runoff and infiltration, which was based on the Richards equation and soil stress–strain mechanism analysis. Also, the RMSD parameters with physical meanings were introduced to depict the soil deformation, which could be measured by the proposed experimental setup or fitted by the suggested calculation formulas. The main objectives of this study were to: (i) propose a modified soil and water movement model for deformable soil; (ii) formulate the effects of dynamic soil deformation on soil hydraulic properties and simplify the parameters which were measured by experimental setup or empirically fitted by the suggested calculation formulas; and (iii) evaluate the approach where the impacts of the swelling pressure and the self-weight stress on porosity changes and water flow are separated. The proposed model and introduced parameters were verified by comparison with the traditional Richards model (TRID) and simulation of Loess soil under rainfall infiltration runoff test. The research will help to improve the theory of soil deformation and the simulation precision of hydrological models in basins, and it will also be conducive to soil and water conservation and agricultural irrigation.

MATERIAL AND METHODS

Governing equations

Figure 1 shows the schematic diagram of soil deformation and soil stress. The soil deformation often performed in the horizontal and vertical direction (Bronswijk 1990), while it was mainly reflected in the vertical direction in the process of water absorption and soil swelling. Soil deformation is affected by soil swelling force F and self-weight stress G. P is the joint force of soil swelling force and self-weight stress, Z1 is the vertical size of the small soil unit before deformation, Z2 is the vertical size of the small soil unit after deformation, Vs and Vs’ are, respectively, the soil volume of the small soil unit before and after deformation.

Figure 1

Schematic diagram of soil deformation (F is swelling force, N; G is self-weight stress, N; P is the resultant force between soil swelling force and self-weight stress, N; Vs is the volume of soil, cm3).

Figure 1

Schematic diagram of soil deformation (F is swelling force, N; G is self-weight stress, N; P is the resultant force between soil swelling force and self-weight stress, N; Vs is the volume of soil, cm3).

For convenient analysis, it is assumed that: (1) soil deformation only causes soil porosity change; (2) soil deformation is elastic and instantaneous. Hence, for a given unit (, of which = 1 × 1, as shown in Figure 2, u is a microelement), the amount of water content in the soil at any unit time is composed of three parts. They are as follows:

Figure 2

Schematic diagram of soil water balance in calculation unit under deformation condition (X and Z are coordinate axes and u is infinitesimal).

Figure 2

Schematic diagram of soil water balance in calculation unit under deformation condition (X and Z are coordinate axes and u is infinitesimal).

Changes in the water content of the unit  
formula
(1)
Changes in the water content caused by the unit change  
formula
(2)
where is the initial soil porosity before deformation (at the beginning of the period), cm3/cm3; is the deformed soil porosity (at the end of the period of soil moisture), cm3/cm3; is the porosity change of the deformable soil, . According to assumption (1), the amount of porosity change is equal to the volume change of soil, cm3/cm3; is the initial soil moisture content (at the beginning of the period), cm3/cm3; is the deformed soil moisture content (at the end of the period of soil moisture), cm3/cm3.
Since is a high-order small amount, it is ignored here. Therefore, the amount of water content of the whole soil within time can be expressed as: 
formula
(3)
In one dimension, the water flux caused by the water potential in the direction of the Z axis in a unit time can be expressed as: 
formula
(4)
 
formula
(5)
 
formula
(6)
where is the water flow, cm3; q is the water flux, cm/min; is the unsaturated soil hydraulic conductivity, cm/min. is the soil water potential, cm; is the hydrostatic pressure (when the hydrostatic pressure is small, it can be ignored), cm; z is the hydraulic heads, cm; is the soil water suction, cm.
During the unit time , the water flux caused by the water absorption of plant roots in the soil can be expressed as: 
formula
(7)
where is the water absorption of plant roots, cm3; S is the water absorption strength of plant roots, min−1.
According to the mass conservation law, the following can be obtained: 
formula
(8)
where is soil volumetric water content, cm3/cm3; e is soil porosity, cm3/cm3; t is time, min; is the initial soil porosity, cm3/cm3; is initial soil moisture, cm3/cm3.

The boundary conditions

Upper boundary

Based on hydraulic head, h, the upper boundary is divided into the following cases.

When rainfall intensity does not exceed the infiltration capacity of the soil, . Then the upper boundary can be described as: 
formula
(9)
where is either depth of surface ponding or negative hydraulic head on the soil surface, cm; is the matric potential of air-dried soil, cm; is the water flux at time t, cm/min.
When rainfall intensity exceeds the infiltration capacity of the soil and the depth of surface ponding is lower than the maximum value, the upper boundary can be described as: 
formula
(10)
where is maximum value of surface ponding, cm.
When the depth of surface ponding is larger than the maximum value, the upper boundary can be described as: 
formula
(11)

Lower boundary

 
formula
(12)

Variable time step methods and implicit difference scheme are used to solve the one-dimensional moving boundary problem in this paper. Solving methods refer to Hydrus (Kool & Van Genuchten 1991; Leij et al. 1991).

Calculation of soil porosity change

Change caused by soil swelling force

It is assumed that the soil deformation is caused by the change of soil porosity when the water uptake and soil is deformed. Therefore, the porosity change caused by the swelling force of soil can be expressed as: 
formula
(13)
where is the soil porosity change caused by water swelling of soil, cm3/cm3; is the period of soil porosity caused by swelling force, cm3/cm3; is the same period of initial porosity, cm3/cm3.
Assuming the soil is homogeneous and its deformation is free-form, the soil swelling characteristic curve is described using the three-line model (McGarry & Malafant 1987): 
formula
(14)
where v is the specific volume that is the reciprocal of the soil bulk density, cm3/g; U is the mass water content, g/g; α1, α2 and α3 are the slope of soil swelling characteristic curve; UA, UB and US are the critical value of the mass water content of the curve, g/g. a1, a2 and a3 are fitting parameters derived from testing according to the three-line model.
Thus, we can obtain the following equation: 
formula
(15)
where is the bulk density change caused by soil swelling force, cm3/g; the others have the same meaning as described above.

Change caused by soil self-weight force

Similarly, the change of porosity caused by soil self-weight stress can be calculated as follows: 
formula
(16)
where is the soil porosity change that is caused by soil self-weight stress, cm3/cm3.
The bulk density change caused by soil self-weight stress can be calculated by soil stress-strain relation curve (Chen et al. 1994) as follows: 
formula
(17)
where is the bulk density change caused by soil self-weight stress, cm3/g; is the wet bulk density, N/cm3; h is the soil depth, cm; A and B are empirical parameters derived from the testing results of Chen et al. (1994). The others have the same meaning as described above.

Change caused by resultant force

The soil deformation is the result of the joint action of the expansibility, collapsibility and the self-weight stress. Thus, the amount of soil porosity changes caused by the resultant force can be expressed as: 
formula
(18)

Introduced model parameters

Deformable soil saturated water conductivity

Influenced by self-weight stress and soil swelling force, the deformable soil porosity changes with the depth of the soil, which leads to the variation of soil saturated water conductivity with the depth. An improved Lambe model (Lambe & Whitman 1979) can be used to calculate the relationship between deformable soil saturated water conductivity and its depth: 
formula
(19)
 
formula
(20)
where is the deformable soil saturated hydraulic conductivity when soil porosity is e, cm/min; and , measured by constant volume method, are respectively porosity and saturated hydraulic conductivity, cm3/cm3, cm/min; is the saturated soil porosity in depth of h, cm3/cm3; is the fitting parameter related to soil texture; is soil density, g/cm3, which is the fitting parameter; is the slope of soil characteristic curve of deformable soil; U is the mass water content of soil, g/g; , a3, A and B are fitting parameters derived from testing with the improved Lambe model.

Deformable soil saturated water content

When soil is saturated, the saturated water content of soil is equal to the soil porosity. Thus, the saturated water content of soil under any soil depth can be expressed as: 
formula
(21)
where (cm3/cm3), is the saturated soil moisture content where the soil depth is h.

Deformable soil–water characteristic curve

The relationship between soil water suction and soil water content θ is described by the improved van Genuchten model (Van Genuchten et al. 1991): 
formula
(22)
where is the soil saturation coefficient; is the residual water content, cm3/cm3; , m, n are parameters, m = 1–1/n. , , m and n also vary with soil porosity. For the same kind of soil, the variations of , , m and n are small. Therefore, , , m and n at the different depths can be replaced with each other.

Deformable soil unsaturated hydraulic conductivity

The unsaturated hydraulic conductivity is calculated by the improved van Genuchten model. 
formula
(23)
where is the unsaturated hydraulic conductivity, cm/min.

Deformable soil–water capacity

 
formula
(24)
where is the water capacity, min−1.

MODEL VERIFICATIONS

Method and materials

Experimental materials

Loessial soil is used to verify the model. Loessial soil, originating from long-term cultivation of loess, is widely distributed in the Loess Plateau where soil erosion is extremely serious. Loessial soil is loose and porous and has characteristics of wet expansion and dry shrinkage and collapsibility deformation (Wang et al. 2006). Physical properties of loessial soil, shown in Table 1, are measured by the pipette method in the test soil, which is dried by wind and screened by 5 mm sieves. To facilitate the analysis of saturated water movement parameters (shown in Table 2), a specific gravity bottle is used to determine the soil density and Vernier caliper method is used to determine the swelling soil characteristic curve. A constant pressure method (shown in Figure 3) is adopted to determine soil stress-strain relationship curves (shown in Table 3). Table 4 shows the curve of bulk density (soil bulk density is 1.1, 1.2, 1.3, 1.4 and 1.5 g/cm3, respectively) and water saturated conductivity which is measured by the constant volume method (shown in Figure 4).

Table 1

Physical properties of loessial soil

Soil particle/%
 
Dry bulk density of soil/(g/cm3density/(g·cm−3Saturated hydraulic conductivity/(cm/min) 
>50 2–50 <2 
μμμ
3.50 52.00 44.5 1.40 2.628 0.014 
Soil particle/%
 
Dry bulk density of soil/(g/cm3density/(g·cm−3Saturated hydraulic conductivity/(cm/min) 
>50 2–50 <2 
μμμ
3.50 52.00 44.5 1.40 2.628 0.014 
Table 2

Parameters of loessial soil simulated by three-line model

α1 R2 α2 R2   α3 R2 UA UB US 
0.70 0.05 0.95 0.62 0.39 0.96  0.73 0.03 0.89 0.24 0.32 0.34 
α1 R2 α2 R2   α3 R2 UA UB US 
0.70 0.05 0.95 0.62 0.39 0.96  0.73 0.03 0.89 0.24 0.32 0.34 
Table 3

Curves of loessial soil stress and strain

R2 
1.093 0.1044 0.99 
R2 
1.093 0.1044 0.99 
Table 4

Parameters of loessial soil simulated by Lambe model

Deformable soil saturated hydraulic conductivity/(cm·min−1m’ R2 
0.014 6.18 0.97 
Deformable soil saturated hydraulic conductivity/(cm·min−1m’ R2 
0.014 6.18 0.97 
Figure 3

Schematic diagram of measuring device of soil stress-strain relationship curve.

Figure 3

Schematic diagram of measuring device of soil stress-strain relationship curve.

Figure 4

Schematic diagram of measuring device of soil hydraulic conductivity.

Figure 4

Schematic diagram of measuring device of soil hydraulic conductivity.

Experimental design

Before testing, the soil screened by 5 mm was evenly loaded into layers (10 cm depth) into the soil trough. A Trime moisture sensor was installed at the depth of the soil trough at 0, 10, 25 and 40 cm, in order to measure the dynamics of soil moisture. Two outlets were set up in the soil trough, and are the interflow and the surface runoff. When loading the soil, we first laid a 10-cm-thick stone with anti-filter action at the bottom of the soil trough. Then, we laid a gauze on the top of the stone. Finally, the soil sieved through a 5 mm screen was loaded into the soil trough up to the designed density (1.4 g/cm3). Before the test, enough water was applied until the soil water content was above the field water capacity (the bottom steel pipe outflowed continuously). The soil sample was then covered with plastic to prevent evaporation and left undisturbed until fully deformed.

During the experiment, the artificial rainfall intensity was 30 mm/h, and the rainfall duration approximately 350 min. The soil water dynamics was recorded by the Trime system in the test. At the same time, the surface runoff and the interflow were collected and measured. The time interval of runoff observation is determined by the velocity of runoff. Because the interflow is slower than the surface runoff, it is necessary to adopt a longer observation time. The observation will carry on until no more interflow can be measured. The slope of the soil trough is fixed at 5° (Figure 5). We designed four rainfall–infiltration–runoff tests for deformable soils with a thickness of 10, 20, 30 and 40 cm, respectively. The tests for each thickness were repeated twice.

Figure 5

Schematic diagram of the experimental setup.

Figure 5

Schematic diagram of the experimental setup.

The main recorded data are runoff (surface runoff and interflow), dynamic water content of soil profile, runoff duration, rainfall intensity, rainfall duration and sediments.

Model parameters

The RMSD model's parameters have clear physical meaning, which can be measured by the experimental setup or calculated by the equations proposed by the authors. When calculating the soil water movement parameters, several model parameters need to be measured in laboratory experiments. These model parameters mainly include the soil swelling characteristic curve, soil stress–strain relationship curve, soil initial bulk density, soil saturated water content and soil saturated hydraulic conductivity corresponding to the soil initial bulk density, soil water characteristic curve corresponding to the soil initial bulk density, parameter m related to the soil porosity. The soil initial water content at different soil depths are measured. The parameters θr, θsh, α and n are changed with soil porosity change. In this paper, the θsh is replaced by eh, and calculated with Equations (20) and (21); θr, α and n are replaced by these parameters corresponding to the soil initial bulk density (Table 5).

Table 5

Input parameters in the RMSD model

a1 α1 a2 α2 a3 α3 UA UB US θr α m' ρd K0 
0.7 0.05 0.62 0.39 0.73 0.03 0.24 0.32 0.34 0.06 0.02 0.422 1.731 1.093 0.1044 6.18 2.628 0.014 
a1 α1 a2 α2 a3 α3 UA UB US θr α m' ρd K0 
0.7 0.05 0.62 0.39 0.73 0.03 0.24 0.32 0.34 0.06 0.02 0.422 1.731 1.093 0.1044 6.18 2.628 0.014 

However, in the TRID, the soil structure is considered to remain constant during the wet–dry alternation, and the soil water movement parameters are unchanged with soil depth. As the TRID is applied in swelling soils, the soil saturated water movement parameters are difficult to determine. Consequently, these parameters corresponding to the soil initial bulk density are measured and then used in the TRID (Table 6).

Table 6

Input parameters in the TRID model

θr θs α Ks 
0.06 0.47 0.02 0.422 1.731 0.014 
θr θs α Ks 
0.06 0.47 0.02 0.422 1.731 0.014 

RESULTS AND DISCUSSION

Runoff intensity

Figure 6 shows the relationship between measured and simulated runoff intensity varying with time. For all soil depth treatments, the runoff intensity increased with time and then stabilized until steady runoff condition.

Figure 6

Relationship between runoff intensity and time at different depths.

Figure 6

Relationship between runoff intensity and time at different depths.

At the start of the rainfall, both the RMSD model and the TRID model performed less well. Deformation was not instantaneous, requiring time to develop. The RMSD model ignored the lagging time and the TRID model did not consider soil deformation at all. In terms of flow time, the runoff intensity of both the RMSD model and the TRID model lagged behind the observed values.

In the later period, the runoff intensity became stable. The TRID overestimated the runoff intensity when the soil depth was 10 cm, and underestimated the runoff intensity when the soil depth was greater than 10 cm. In comparison, the simulation runoff intensity by the RMSD model was closer to the measured value and better fitted.

To further compare the simulation accuracy of the two models, ARE (average relative error) and NSE (Nash efficiency coefficient) values were adopted to reflect the performances of the RMSD and the TRID in simulating runoff intensity. The results are shown in Table 7.

Table 7

The relative error and Nash efficiency coefficient between measured and simulated runoff intensity

Models Depth ARE/(%) NSE 
TRID 10 cm 0.74 0.58 
20 cm −13.82 0.54 
30 cm −17.01 0.69 
40 cm −27.47 0.49 
RMSD 10 cm −4.04 0.70 
20 cm −9.49 0.60 
30 cm −8.78 0.69 
40 cm −1.38 0.58 
Models Depth ARE/(%) NSE 
TRID 10 cm 0.74 0.58 
20 cm −13.82 0.54 
30 cm −17.01 0.69 
40 cm −27.47 0.49 
RMSD 10 cm −4.04 0.70 
20 cm −9.49 0.60 
30 cm −8.78 0.69 
40 cm −1.38 0.58 

For all soil depth treatments, ARE values between the simulated values of the RMSD and the measured runoff rate were in the range of −10.0% to −1.0%, and the NSE values were all larger than 0.58. While the simulated values of the TRID ARE values were −28.0% to 1.0%, the NSE values were between 0.49 and 0.69. According to NSE values of the TRID model, the TRID model incurred analogue distortion, which illustrated that the TRID model is not fitted for soil deformation. Compared with ARE and NSE values of the RMSD, it could be concluded that the RMSD significantly improved the simulation effectiveness in runoff intensity by considering the effect of soil deformation on the soil water movement process.

Cumulative infiltration amount

Figure 7 shows the results of measured and simulated cumulative infiltration varying with time. The cumulative infiltration amount of soil increased with the time. When the soil was saturated, the measured and simulated soil cumulative infiltration amount changed abruptly. The reason for this phenomenon might be that the water absorption capacity of saturated soil disappeared and the soil water moved under the action of gravity.

Figure 7

Relationship between cumulative soil infiltration and time.

Figure 7

Relationship between cumulative soil infiltration and time.

As time went on, the value simulated by the RMSD model was gradually equal to the measured value. Generally, the simulated value of the RMSD model was closer to the measured value. The TRID without considering soil deformation led to overestimating or underestimating the cumulative infiltration of the soil. In addition, at the initial period of infiltration, the RMSD model overestimated the cumulative infiltration amount of soil. The possible reason was that the soil deformation was not only determined by soil water content and self-weight stress, but also related to the time. During the process of soil infiltration, the soil uptakes water, expands and swells, which is a process that requires time to develop. In other words, the processes of deformation were not instant. However, the RMSD model assumed that soil deformation was an elastic and instantaneous deformation, an assumption which eliminated the model's ability to simulate the effect of deformation lagged on the real infiltration process of soil.

When the deformable soil absorbed water, the soil parameters such as porosity, soil moisture, self-weight and swelling force varied with depth, water content and time. When swelling soils absorbed water, the swelling deformation of soils occurred under the action of the soil swelling pressure and self-weight stress. When the soil layer was thin, the self-weight stress was much smaller than the soil swelling pressure. Under the domination of the soil swelling pressure, the soil volume increased which led to the increasing of soil porosity, the soil saturated water content and soil saturated hydraulic conductivity. With the increase of soil thickness, the soil self-weight stress increased and gradually became the dominant factor. Under the domination of the soil self-weight stress, the soil volume decreased, the soil porosity decreased, and the soil saturated water content and soil saturated hydraulic conductivity decreased. Consequently, when the soil layer was thin (10 cm), the soil infiltration was underestimated by the TRID, which caused the simulated runoff rate to be larger than the measured value. Conversely, when the soil layer was thick (>10 cm), the runoff rate was underestimated by the TRID. When the self-weight was smaller than the swelling force (depth <10 cm), expansion happened. Conversely, when the self-weight was larger than the swelling force (depth >10 cm), the soil collapsed.

Similarly to the analysis of runoff intensity, ARE and NSE values were applied to reflect the performances of the RMSD and the TRID in simulating the cumulative infiltration of loessial soil. As shown in Table 8, for all soil depth treatments, ARE values between the simulated values of the RMSD and the measured cumulative infiltration of loessial soil were in the range of −1.0% to 11.0%, and the NSE values were all larger than 0.90. The simulated values of the TRID ARE values were −10.0 to 32.0% and the NSE values were between 0.58 and 0.94. It could be concluded that the RMSD fitted more the observations of the cumulative infiltration of loessial soil.

Table 8

The relative error and Nash efficiency coefficient between the measured and simulated cumulative infiltration

Model Depth ARE/(%) NSE 
TRID 10 cm −9.02 0.94 
20 cm 10.67 0.94 
30 cm 12.59 0.93 
40 cm 31.67 0.58 
RMSD 10 cm −0.92 0.99 
20 cm 3.25 0.98 
30 cm 5.7 0.98 
40 cm 10.29 0.94 
Model Depth ARE/(%) NSE 
TRID 10 cm −9.02 0.94 
20 cm 10.67 0.94 
30 cm 12.59 0.93 
40 cm 31.67 0.58 
RMSD 10 cm −0.92 0.99 
20 cm 3.25 0.98 
30 cm 5.7 0.98 
40 cm 10.29 0.94 

Conclusively, the RMSD performed better than the TRID in simulating soil and water movement process, owing to the consideration of soil deformation.

Discussion

Some research focus on soil deformation resulted from either dynamic compaction (e.g., tractor traffic, tillage, cattle trampling) (Berli et al. 2003; Cui et al. 2010; Nawaz et al. 2013), hydraulic stress from natural wetting-and-drying cycles or soil swelling force (Garnier et al. 1997; Baveye 1998). In the RMSD model, the dynamic soil deformation induced by dynamic water content, swelling pressure and the self-weight stress were fully studied. Influenced by self-weight stress and soil swelling force, the deformable soil porosity changed with the depth of the soil, which led to the dynamic variation of soil saturated water conductivity with the depth.

By the analysis of stress characteristics for soil deformation, the RMSD model for expansive soil water movement was proposed based on the mass conservation law. Some new parameters including the deformable soil saturated hydraulic conductivity and deformable soil saturated water content in deformed soil were introduced. The newly introduced parameters originated from the physical properties of soil and calculated by the RMSD equations, which were easy to obtain. Therefore, the introduced parameters in the model had clear physical meaning, which would be helpful to model analysis and applications. Meanwhile, the RMSD model can not only simulate the soil infiltration process, but can also simulate soil water redistribution, the root water uptake. Hence, the RMSD model has more extensive applicability. What is more, as for a simplified semi-empirical deformable soil model, we argue that the RMSD would likely be more effective to improve the quantification of soil ecosystem processes and identify the key challenges in upscaling deformable soil properties and processes measured in the laboratory to the field scale.

It should be pointed out that two simplified assumptions are adopted in the section Governing equations to conveniently analyze the process of soil deformation. In fact, the deformation may cause residual soil particle size change and itself includes lagging time, thus the assumptions may need further evaluation. However, even with the two assumptions, the RMSD model still shows improvement in simulating deformable soils.

CONCLUSION

A modified model, the RMSD, was proposed for simulating soil water movement process of deformable soils. The RMSD was based on the Richards model. In the model, new parameters were introduced to quantify the soil deformation, including the deformable soil saturated hydraulic conductivity and the deformable soil saturated water content in deformed soil. The parameters could be determined with physical properties and suggested formulations, which makes the RMSD convenient for application. Although a further evaluation of the model assumptions is desired, the model verifications showed the RMSD is better than the TRID in fitting the observations of soil cumulative infiltration and the runoff intensity, illustrating that the newly suggested model is a suitable tool to depict the soil water movement in deformable soils.

ACKNOWLEDGEMENTS

This study is supported by the National Basic Research Program of China (973 Program) (2015CB452701), the Projects of National Science Foundation of China (51379215, 51779272), and the Project of State Key Lab. of Simulation and Regulation of Water Cycle in River Basin (SKL-WAC) (2016ZY03).

REFERENCES

REFERENCES
Alaoui
A.
,
Lipiec
J.
&
Gerke
H. H.
2011
A review of the changes in the soil pore system due to soil deformation: a hydrodynamic perspective
.
Soil and Tillage Research
115–116
(
5
),
1
15
.
Assouline
S.
2006a
Modeling the relationship between soil bulk density and the water retention curve
.
Vadose Zone Journal
5
,
554
563
.
doi:10.2136/vzj2005.0083
.
Assouline
S.
2006b
Modeling the relationship between soil bulk density and the hydraulic conductivity function
.
Vadose Zone Journal
5
,
697
705
.
doi:10.2136/vzj2005.0084
.
Barontini
S.
,
Ranzi
R.
&
Bacchi
B.
2007
Water dynamics in a gradually nonhomogeneous soil described by the linearized Richards equation
.
Water Resources Research
43
(
8
),
1050
1056
.
Berli
M.
,
Kirby
J.
,
Springman
S.
&
Schulin
R.
2003
Modelling compaction of agricultural subsoils by tracked heavy construction machinery under various moisture conditions in Switzerland
.
Soil and Tillage Research
73
,
57
66
.
doi:10.1016/S0167-1987(03)00099-0
.
Bronswijk
J. J. B.
1990
Shrinkage geometry of a heavy clay soil at various stresses
.
Soil Science Society of America Journal
54
(
5
),
1500
1502
.
Chen
Z.
,
Zhou
J.
&
Wang
H.
1994
Soil Mechanics
.
Tsinghua University Press
,
Beijing
,
China
.
Cui
Y.
&
Miao
L. C.
2011
Experimental study on characteristics of unsaturated compacted expansive soil permeability
.
Chinese Journal of Rock Mechanics and Engineering
32
(
7
),
2007
2012
.
Cui
K.
,
Défossez
P.
,
Cui
Y. J.
&
Richard
G.
2010
Soil compaction by wheeling: changes in soil suction caused by compression
.
European Journal of Soil Sciences
61
,
599
608
.
doi:10.1111/j.1365-2389.2010.01245.x
.
Dam
J. C. V.
&
Feddes
R. A.
2000
Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation
.
Journal of Hydrology
233
(
1–4
),
72
85
.
Gan
Y.
,
Jia
Y.
,
Wang
K.
,
Hao
W.
&
Wei
N.
2015
Model of rainfall infiltration in layered soil considering air resistance
.
Journal of Hydraulic Engineering
46
(
2
),
164
173
.
Garnier
P.
,
Perrier
E.
,
Jaramillo
R. A.
&
Baveye
P.
1997
Numerical model of 3-dimensional anisotropic deformation and 1-dimensional water flow in swelling soils
.
Soil Science
162
(
6
),
410
420
.
Ghezzehei
T. A.
&
Or
D.
2001
Rheological properties of wet soils and clays under steady and oscillatory stresses
.
Soil Science Society of America Journal
65
,
624
637
.
doi:10.2136/sssaj2001.653624x
.
Giráldez
J. V.
,
Delagado
C.
&
Sposito
G.
1983
A general soil volume change equation: I. The two-parameter model
.
Soil Science Society of America Journal
47
(
3
),
419
422
.
Horn
R.
&
Fleige
H.
2003
A method for assessing the impact of load on mechanical stability and on physical properties of soils
.
Soil and Tillage Research
73
,
89
99
.
doi:10.1016/S0167-1987(03)00102-8
.
Kool
J. B.
&
Van Genuchten
M. T.
1991
HYDRUS-One-dimensional Variably Saturated Flow and Transport Model, Including Hysteresis and Root Water Uptake
.
Version 3.3, Research Report No. 124
,
U.S. Salinity Lab, U.S. Depart ment of Agriculture, Agricultural Research Service
,
Riverside, CA
,
USA
.
Lambe
T. W.
&
Whitman
R. V.
1979
Soil Mechanics
.
SI version, Wiley
,
New York
,
USA
.
Lei
Z.
,
Yang
S.
&
Xie
S.
1988
Soil Water Dynamics
.
Tsinghua University Press
,
Beijing
,
China
.
Leij
F. J.
,
Skaggs
T. H.
&
Van Genuchten
M. T.
1991
Analytical solutions for solute transport in three-dimensional semi-infinite porous media
.
Water Resources Research
27
(
10
),
2719
2733
.
Li
X.
,
Wang
H.
,
Sun
H.
,
Zhang
J.
&
Sun
X.
2009
Numerical simulation of soil water movement based on MATLAB
.
Science Technology and Engineering
9
(
20
),
5978
5981
.
Lindstrom
M. J.
,
Voorhees
W. B.
,
Schumacher
T. E.
1995
Soil properties across a landscape continuum as affected by planting wheel traffic
. In:
Site Specific Management for Agricultural Systems
(
Robert
P. C.
,
Rust
R. H.
&
Larson
W. E.
, eds).
ASA-CSSA-SSSA, Madison, WI, USA
, pp.
351
363
.
Loáiciga
H. A.
&
Huang
A.
2007
Ponding analysis with Green-and-Ampt infiltration
.
Journal of Hydrologic Engineering
12
(
1
),
109
112
.
Logsdon
S.
,
Berli
M.
&
Horn
R.
, (eds.)
2013
Quantifying and Modelling Soil Structure Dynamics. Advances in Agricultural Systems Modeling
, Vol.
3
.
SSSA
,
Madison, WI
,
USA
.
Lu
C.
&
Pei
Y.
2007
Numerical simulation of one-dimensional soil water motion for complex upper boundary conditions
.
Chinese Journal of Hydraulic Engineering
2
,
136
142
.
McGarry
D.
&
Malafant
K. W. J.
1987
The analysis of volume changes in unconfined units of soil
.
Soil Science Society of America Journal
51
(
2
),
290
297
.
Mishra
S. K.
,
Tyagi
J. V.
&
Singh
V. P.
2003
Comparison of infiltration models
.
Hydrological Processes
17
,
2629
2652
.
Nawaz
M. F.
,
Bourrié
G.
&
Trolard
F.
2013
Soil compaction impact and modelling. A review
.
Agronomy for Sustainable Development
33
(
2
),
291
309
.
Peth
S.
,
Horn
R.
,
Fazekas
O.
&
Richards
B. G.
2006
Heavy soil loading and its consequence for soil structure, strength, and deformation of arable soils
.
Journal of Plant Nutrition and Soil Science
169
,
775
783
.
doi:10.1002/jpln.200620112
.
Raats
P. A. C.
&
Klute
A.
1968
Transport in soils: the balance of mass
.
Soil Science Society of America Journal
32
(
2
),
161
166
.
Richards
L. A.
1931
Capillary conduction of liquids through porous mediums
.
Journal of Applied Physics
1
(
5
),
318
333
.
Richards
B. G.
,
Baumgartl
T.
,
Horn
R.
&
Gräsle
W.
1997
Modelling the effects of repeated wheel loads on soil profiles
.
International Agrophysics
11
(
3
),
177
187
.
Smiles
D. E.
1997
Water balance in swelling materials: some comments
.
Australian Journal of Soil Research
35
(
5
),
1143
1152
.
Smiles
D. E.
2000
Hydrology of swelling soils: a review
.
Australian Journal of Soil Research
38
(
3
),
501
521
.
Strudley
M. W.
,
Green
T. R.
&
Ii
J. C. A.
2008
Tillage effects on soil hydraulic properties in space and time: state of the science
.
Soil and Tillage Research
99
(
1
),
4
48
.
Van Genuchten
M. T.
,
Leij
F. J.
&
Yates
S. R.
1991
The RETC Code for Quantifying the Hydraulic Functions of Unsaturated Soils
.
U.S. Salinity Laboratory, U.S. Department of Agriculture, Agricultural Research Service
,
Riverside, CA
,
USA
.
Vereecken
H.
,
Schnepf
A.
,
Hopmans
J. W.
,
Javaux
M.
,
Or
D.
,
Roose
T.
&
Allison
S. D.
2016
Modeling soil processes: review, key challenges, and new perspectives
.
Vadose Zone Journal
15
(
5
),
1
57
.
Wang
F.
&
Yang
H.
2013
A finite element mumerical simulation for unsaturated Richards model with complex boundary
.
Journal of Wenzhou University
34
(
01
),
31
35
.
Wang
Y.
,
Wang
Y.
,
Liu
J.
,
Yan
Y.
&
Gao
X.
2005
Study on factors affecting soil swelling in loess region
.
Agricultural Research in Arid Area
23
(
5
),
93
97
.
Wang
H.
,
Yang
G.
,
Jia
Y.
&
Wang
J.
2006
Connotation and evaluation index system of soil and water resources
.
Journal of Hydraulic Engineering
37
(
4
),
389
394
.
Wiermann
C.
,
Way
T. R.
,
Horn
R.
,
Bailey
A. C.
&
Burt
E. C.
1999
Effect of various dynamic loads on stress and strain behavior of a Norfolk sandy loam
.
Soil and Tillage Research
50
,
127
135
.
doi:10.1016/S0167-1987(98)00199-8
.
Ying
M.
,
Feng
S.
,
Su
D.
,
Gao
G.
&
Huo
Z.
2010
Modeling water infiltration in a large layered soil column with a modified Green–Ampt model and HYDRUS-1D
.
Computers & Electronics in Agriculture
71
(
1
),
S40
S47
.
Zhang
H.
,
Chen
S.
&
Chen
S.
2003
Numerical simulation of infiltration of unsaturated soil
.
Rock and Soil Mechanics
24
(
5
),
715
718
.
Zhang
S.
,
Grip
H.
&
Lővdahl
L.
2006
Effect of soil compaction on hydraulic properties of two loess soils in China
.
Soil and Tillage Research
90
,
117
125
.
Zhou
B.
,
Kong
L.
,
Chen
W.
,
Bai
H.
&
Li
X.
2010
Jingmen expansive soil water characteristic curve analysis of characteristic parameters and the unsaturated shear strength prediction
.
Chinese Journal of Rock Mechanics and Engineering
29
(
5
),
1052
1059
.