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 a schematic diagram of soil deformation and soil stress. Soil deformation often occurs in the horizontal and vertical direction (Bronswijk 1990), while it is 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.
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).
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.
Schematic diagram of soil water balance in calculation unit under deformation condition (X and Z are coordinate axes and u is infinitesimal).
Schematic diagram of soil water balance in calculation unit under deformation condition (X and Z are coordinate axes and u is infinitesimal).













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




Lower boundary
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



Change caused by soil self-weight force


Change caused by resultant force
Introduced model parameters
Deformable soil saturated water conductivity








Deformable soil saturated water content
Deformable soil–water characteristic curve










Deformable soil unsaturated hydraulic conductivity
Deformable soil–water capacity
MODEL VERIFICATIONS
Method and materials
Experimental materials
Loessial soil is used to verify the model. Loessial soil, originating after long-term cultivation, 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).
Physical properties of loessial soil
Soil particle/% . | Dry bulk density of soil/(g/cm3) . | density/(g/cm−3) . | Saturated hydraulic conductivity/(cm/min) . | ||
---|---|---|---|---|---|
>50 . | 2–50 . | <2 . | |||
μm . | μm . | μm . | |||
3.50 | 52.00 | 44.5 | 1.40 | 2.628 | 0.014 |
Soil particle/% . | Dry bulk density of soil/(g/cm3) . | density/(g/cm−3) . | Saturated hydraulic conductivity/(cm/min) . | ||
---|---|---|---|---|---|
>50 . | 2–50 . | <2 . | |||
μm . | μm . | μm . | |||
3.50 | 52.00 | 44.5 | 1.40 | 2.628 | 0.014 |
Parameters of loessial soil simulated by three-line model
a . | α1 . | R2 . | b . | α2 . | R2 . | . | c . | α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 |
a . | α1 . | R2 . | b . | α2 . | R2 . | . | c . | α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 |
Curves of loessial soil stress and strain
A . | B . | R2 . |
---|---|---|
1.093 | 0.1044 | 0.99 |
A . | B . | R2 . |
---|---|---|
1.093 | 0.1044 | 0.99 |
Parameters of loessial soil simulated by Lambe model
Deformable soil saturated hydraulic conductivity/(cm·min−1) . | m’ . | R2 . |
---|---|---|
0.014 | 6.18 | 0.97 |
Deformable soil saturated hydraulic conductivity/(cm·min−1) . | m’ . | R2 . |
---|---|---|
0.014 | 6.18 | 0.97 |
Schematic diagram of measuring device of soil stress-strain relationship curve.
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.
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).
Input parameters in the RMSD model
a1 . | α1 . | a2 . | α2 . | a3 . | α3 . | UA . | UB . | US . | θr . | α . | m . | n . | A . | B . | 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 . | n . | A . | B . | 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).
Input parameters in the TRID model
θr . | θs . | α . | m . | n . | Ks . |
---|---|---|---|---|---|
0.06 | 0.47 | 0.02 | 0.422 | 1.731 | 0.014 |
θr . | θs . | α . | m . | n . | 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.
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.
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.
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.
The average relative error (ARE) and Nash efficiency coefficient (NSE) 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 Laboratory of Simulation and Regulation of Water Cycle in River Basin (SKL-WAC) (2016ZY03).