Analysis of the Hydraulic Properties of Undisturbed Layered Loess in Northwest China

Abstract

Extensive agricultural irrigation in the loess region of Northwest China has seriously damaged the local hydrogeological environment. To properly understand the hydrological processes and the hydraulic properties of the layered soil, the field soil column irrigation test, laboratory soil column infiltration test, and undisturbed soil sample hydraulic experiments were carried out. The results showed that the proposed infiltration model can continuously simulate the infiltration process of the loess–palaeosol sequence well. The layered structure may form a temporary groundwater table at the interface of the two different soils under irrigation conditions. This provides a scientific basis for proposing reasonable irrigation measures.

1. Introduction

The Chinese Loess Plateau, with widespread thick loess–palaeosol sequences, is largely a product of the Quaternary eolian activities. During recent years, South Jingyang Platform, Shaanxi Province, has received significant attention as a loess region in which landslides frequently occur due to extensive agricultural irrigation [1,2,3]. To prevent serious economic losses and casualties, many researchers have reasonably explored the mechanism of loess landslides, which is closely related to the variation in hydrological conditions from the perspective of constitutive relations [4,5,6]. However, few studies have focused on the specific hydrologic processes and hydrogeological properties of the typical loess–palaeosol sequences under the irrigation condition, which has practical significance to guide seasonal agricultural irrigations.

Many researchers have investigated the variation of the water content profile during the infiltration process of layered soil through indoor or field experiments [7,8]. Zohrab et al. [9] showed that when the sublayer is more pervious than the top layer, the wetting front becomes unstable through the second layer. They determined the ratio of the hydraulic conductivities at which the wetting front loses its stability. Li et al. [10] investigated the effects of layered soil on wetting patterns and water distributions from a surface point source. They concluded that an interface existing in the layered soil, whether fine-over-coarse or coarse-over-fine, had a common feature of limiting downward water movement. Yang et al. [11] evaluated the effects of rainfall intensity and duration on finer over coarser layered soils columns. They found that rainfall intensity had a major effect on infiltration in the finer layer but had limited effect on the coarse layer due to the large difference in saturated permeability between the two layers. The variation in cumulative infiltration and infiltration rate with time can also be used to study the characteristics of water transport in layered soil. There are many models that can simulate the relationship between cumulative infiltration and time in homogeneous soils [12,13], such as the Philip model, the Green–Ampt model, and the Kostiakov model. However, infiltration into layered soils is more complex compared to homogeneous soils. Previous research has shown that layered structure has a blocking effect on water movement and has proposed infiltration models applicable to layered soils [14,15,16,17]. The Green–Ampt model has been modified to describe the infiltration characteristics of the layered soils [18,19,20]. However, fewer studies have focused on the infiltration models of undisturbed loess–palaeosol sequences.

Basic soil hydraulic properties, such as the soil water characteristic curve (SWCC) and the hydraulic conductivity, can be used to analyse the mechanisms of the water blocking effect of the layered soils [21,22]. They can be directly tested under laboratory or field conditions [23,24]. Traditional experimental methods include the pressure membrane, tensiometer, and centrifuge methods. However, the lack of consideration of the scale effect in a numerical model may misrepresent the responses and the parameters obtained from the soil samples may not be able to accurately match the parameters of the model with different scales [25]. We can also combine stable numerical solution schemes of the governing flow equations with efficient parameter optimization methods to find the values of the soil hydraulic parameters by inverse modelling [26,27]. During recent years, the Bayesian methods, combining a prior distribution of the soil hydraulic parameters with the soil water state variable observations using Bayes’ theorem, have found widespread application [28]. The Markov chain Monte Carlo (MCMC) simulation has been used recently to estimate the posterior distribution of the parameters [29,30,31]. The DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm has its roots within the Monte Carlo (MC) method, using subspace sampling and outlier chain correction to accelerate convergence to the target distribution [32]. The specific difference of the parameters derived from both the direct methods and the inverse methods for the Q3 loess and the S1 palaeosol have not been studied.

In this paper, three experiments with different scales were carried out to investigate the hydrological processes of the loess–palaeosol sequence during infiltration and the corresponding water-blocking mechanism was analysed. The main contents of this paper are:

• Exploring the variation in the water content profile in the layered loess–palaeosol sequence through the field undisturbed soil column infiltration test;
• Studying the variation in cumulative infiltration with time in the layered soils and infiltration models through the laboratory soil column experiment;
• Deriving the hydraulic properties of different soils through direct and indirect methods and further discussing the impact of the layered structure on hydrological processes under infiltration conditions.

2. Materials and Methods

2.1. Experimental Materials and Methods

2.1.1. Geological Environmental Conditions at the Study Area

The South Jingyang Platform is located between the city of Xianyang and Jingyang country, with an area of about 70 km². The elevation difference between the top and bottom of the platform varies from about 50 to 100 m and its slope ranges from 40° to 90° (Figure 1). The angle of the upper part of the slope is approximately 60° with a height difference of 10–15 m and the angle of the lower slope is approximately 50° with a height difference of 40–80 m. Malan loess (Q3 loess) and Lishi loess (Q2 loess) are the main constituents of the carcass. In the study area, the typical loess irrigated area of Shutangwang Village was selected (the geographical coordinates are 108.84° E, 34.49° N) for the field irrigation infiltration test and the acquisition of the indoor experiment soil samples.

2.1.2. The Field Irrigation Experiment

The purpose of this field experiment was to observe the water infiltration process under the irrigation condition. As shown in Figure 2a, we manually excavated the soil column with 200 cm in both length and width and 255 cm in height, sealed the surrounding of the soil column with bentonite and wrapped plastic film around the soil. One vertical side of the soil column was left uncovered to embed the sensors. The soil column was combined with two different soils, the upper 105 cm Q3 loess and the lower 150 cm S1 palaeosol. As shown in Figure 3, 8 Time-Domain Reflectometry (TDR)-3 moisture sensors were embedded along the depth of the soil column and the corresponding depth was 15 cm, 45 cm, 75 cm, 105 cm, 135 cm, 165 cm, 195 cm, 225 cm from the top surface of the soil column. We connected the installed sensors to the CR3000 data logger, backfilled the sensor holes with bentonite, and finally sealed the boundary with bentonite. Figure 2b shows a water retaining plate arranged above the soil column to provide a 6 cm constant head boundary condition during the irrigation process. The data logger acquisition interval was 3 min in order to achieve continuous, real-time observations of the entire process of irrigation. The irrigation test stopped when the moisture sensor at the 2.25 m depth responded, which happened after 58.65 h.

2.1.3. Large Laboratory Soil Column Infiltration Experiment

The undisturbed loess column was collected by carefully introducing a rigid polyvinyl chloride (PVC) cylinder into the soil. The PVC cylinder was 100 cm in height and 30 cm in inner diameter, with a wall thickness of 1 cm. The specific sampling process was as follows: first, the soil around the target location was excavated and slowly the upper part of the clod reserved in the middle was cut into a small size that was the same as the inner diameter of the cylinder; second, the cylinder was placed on the top of the clod and the cylinder was slightly pressed to place the undisturbed soil within it; third, a spatula was used to cut the next soil section and the cylinder pressed again until the soil reached the top of it; finally, the soil column was shovelled from the root to obtain the undisturbed soil column for the laboratory infiltration experiment. Note that the inner wall of the PVC cylinder was coated with Vaseline to ensure close contact with the soil to prevent preferential flow during the experimental procedure. During the soil column excavation procedure, the undisturbed small soil samples were obtained to satisfy the requirement of the hydraulic experiments in the laboratory.

We removed 3 cm of soil from the top of the column, leaving 47 cm of loess and 50 cm of palaeosol, to eliminate the influence of the disturbed soil. In addition, a sheet of filter paper and a permeable stone were placed at the bottom of the soil column. Afterwards, the soil column was placed on the base, which had a hole in the centre for water drainage or air expulsion. We sealed the gap between the side wall of the PVC cylinder and the base with glass glue to ensure that no water flowed out from here during the test. As shown in Figure 4, the soil column was equipped with five tensiometers (T1 to T5), installed at the corresponding reserved holes 12, 27, 47, 67, and 87 cm from the top of the soil column, respectively. A sheet of filter paper was placed on the top of the soil specimen to avoid disturbing the soil sample during the test. The water supply container was connected and an overflow hole was provided 2 cm from the top of the PVC cylinder to maintain a 1 cm constant water head. The tensiometers were connected to the CR3000 data logger and we observed the corresponding readings using a computer to judge whether the test device was normal. The data acquisition interval was 1 min. The experiment ended when the response of the fifth tensiometer was terminated, which was after 575 min. All sensors responded well.

2.1.4. Test of the Basic Physical and Hydraulic Properties of the Study Soil

A Bettersize 2000 laser particle size analyser (Dandong Better Instrument Corporation, Dandong, China) was applied to determine the grain size distribution of the soil, which is shown in

Table 1. Grain size distribution.

Granular Group (Particle Size Range)		Content Ratio (%)
		Q3 Loess Palaeosol
Sand (0.05–2 mm)		3.91	8.95
Silt (0.005–0.05 mm)	Coarse silt (0.01–0.05 mm)	67.88	56.57
	Fine silt (0.005–0.01 mm)	1.76	2.78
Clay $(<$0.005 mm)		$26.44$	31.7

. We used the HYPROP device and the Ksat device produced by the Universiti Malaysia Sabah (UMS) Company in Germany to obtain the SWCC and the saturated hydraulic conductivity, respectively (Figure 5). The saturated hydraulic conductivity of the loess was 1.33 × ${10}^{-6}$ m/s, which was larger than that of the palaeosol (1.19 × ${10}^{-6}$ m/s).

2.2. Theory of the Soil Hydraulic Properties

2.2.1. Infiltration Model of the Layered Soil

The relationship between the cumulative infiltration and time for layered soils can be expressed by piecewise formulas, meaning that the nonlinear section is simulated by existing models for homogeneous soils and the linear section is simulated by a linear equation. It is important to determine the stable infiltration rate $f_p$ and the transition time $t_p$ in this method. A large number of experimental and theoretical analyses have proven that the infiltration transition of layered soil from a nonlinear process to a linear process exists [18]. The infiltration process of the upper soil layer can be seen as a superposition of two parts. One part is a stable infiltration. The other part can be regarded as a flushing process, which can be described by the exponential relationship using Equation (1). The cumulative infiltration is the superposition of the flushing process and the stable infiltration [33].The schematic diagram of the proposed model is shown below:

As can be seen from Figure 6, Line 1 is a linear model and Line 2 is a nonlinear model. Line 3 is a superposition model of Line 1 and Line 2, and Line 4 is an asymptote of Line 2. The model form of Line 2 is:

$$Y=Y_0\left(1-e^{-\frac{x}{x_0}}\right),$$

(1)

where $Y_0$ is the value of the asymptote of Line 2, which is the value of Line 4; $x_0$ is the time constant, which is generally 1/3 of the transition time. For the infiltration problem of layered soils, the water flushing process of the upper soil layer adopts the relationship of the Equation (1), and the stable infiltration part is simulated by the relationship of the Line 1. It can be seen from Figure 6 that the value of the Line 2 no longer increases after the transition time and the infiltration process becomes linear. The infiltration model can be obtained:

$$I=f_{p}t+f_b(1-e^{-\frac{t}{\tau }});$$

(2)

$$i=f_p+\left(\frac{f_b}{\tau }\right)e^{-\frac{t}{\tau }},$$

(3)

where I is the cumulative infiltration (cm); $f_p$ is the steady infiltration rate (cm ${\min }^{-1}$), determined by the test data fitting; t is the time (min); $f_b$ is the intercept of the linear part and the vertical axis (cm), and the value can also be determined by fitting; τ is the turning time constant (min); $i$ is the infiltration rate (cm ${\min }^{-1}$).

2.2.2. Description of the Basic Soil Hydraulic Properties

SWCC defines the functional relationship between pressure head and soil water content. Many empirical formulas have been proposed to describe SWCC, such as the Brooks–Corey, Gardner, van Genuchten (VG), and Gardner–Russo models. In this paper, we used the VG model to fit the experimental data and used Origin 9.1 software to fit the model parameters. In this study, HYDRUS-1D software was used to calculate the forward model. We adopted the Bayesian approach and quantified parameter uncertainty using the DREAM(zs) algorithm. Time series of the pressure head measured at five different depths of the laboratory soil column were used for soil hydraulic parameter inversion. In the Bayesian framework, the parameter posterior distribution p (m|$\tilde{d}$) can be derived from the experimental data $\tilde{d}$ using Bayes’ law, as follows:

$$\mathrm{p}(m|\tilde{d})=\frac{\mathrm{p}\left(m\right)\mathrm{L}(\tilde{d}|m)}{\mathrm{p}\left(\tilde{d}\right)},$$

(4)

where p (m) is the prior distribution of the parameter vector m, L (m|$\tilde{d}$) ≡ p ($\tilde{d}$|m) is the likelihood function that summarizes the level of agreement between the observed and simulated data, and $\mathrm{p}(\tilde{d})$ is the normalized constant of the probability density function. We resorted to the efficient DREAM(zs) algorithm to numerically estimate the posterior distribution p (m|$\tilde{d}$). The DREAM(zs) algorithm includes various statistical tests to determine when the convergence of the sampled chains to a limiting distribution has been achieved. The most powerful convergence test is the multi-chain $\hat{\mathrm{R}}$-statistic [34]. The convergence exactly occurs as soon as the $\hat{\mathrm{R}}$-statistic decreases to less than the critical threshold 1.2 for all the parameters of the target distribution; otherwise, one should increase the number of iterations and the chains run for a longer period.

3. Results and Discussion

3.1. The Field Infiltration Experiment

Figure 7a shows that soil water content at the position of moisture sensors 1–8 changed with the irrigation time. The variation in soil water content experienced four stages: stable state, sharply rising state, slowly rising, and stable state again. Figure 7b shows that when the saturated zone of the loess reached a certain depth, it would not move down. The water migrated downward in the unsaturated state in the soil, so the actual migration form of water in the loess was dominated by unsaturated forms. At the same time, Figure 7 and

Table 2. The measured water contents at the start and end of the measurement (%).

Sensor Number	1	2	3	4	5	6	7	8
T = 0 h	15	18	20	20	25	30	25	23
T = 58.65 h	38	31	30	29	35	39	34	29

demonstrate that the water content of the palaeosol was significantly greater than that of the loess, which could be explained by the high content of clay in palaeosol soil. This is consistent with the conclusion in [8]. Figure 7b shows that the water content above and below the position of sensor 4 was discontinuous. From the perspective of energy, when the water moves from one layer to another, the water potential at the interface should be continuous, inevitably causing the soil moisture in the upper and lower layers of the interface to be discontinuous.

3.2. Laboratory Soil Column Infiltration Experiment Results

3.2.1. Variance of the Matric Suction during the Infiltration Process

As we can see from Figure 8, the changes in the matric suction were divided into three stages, as follows: stable before the wetting front arrived, sharply decreasing when the wetting front arrived, and basically stable after a period of time. The matric suction values were relatively high for all sensors during the initial stage of the infiltration. When the wetting front arrived, the infiltration capacity of the soil was high and the matric suction rapidly decreased. As time passed, the decreasing rate of the matric suction slowed. After a certain period of time, the soil matric suction was stable. Comparatively speaking, the closer the tensiometer was to the top surface, the faster the suction changed. It can be seen from the figure that the migration velocity of the wetting front below 47 cm was significantly slowed down, which indicates that the infiltration rate of water in the lower layer of palaeosol was reduced.

3.2.2. Infiltration Models

The studied layered soil was the combination of the upper Q3 loess and the lower S1 palaeosol, and the infiltration test lasted for 575 min. Before the wetting front reached the interface between loess and palaeosol (t $\le$ 100 min), the change process of cumulative infiltration with time was the same as that of the homogeneous soil. The relationship between the cumulative infiltration and the time was nonlinear, which can be simulated by Kostiakov’s infiltration model, as follows:

$$I=0.1966\times t^{0.7336}(0<t\le 100)R^2=0.9968,$$

(5)

where I is the cumulative infiltration amount (cm), t is the time (min), and $R^2$ is the relationship coefficient.

When the wetting front reached the interface between the two soils (t $>$ 100 min), the infiltration process was obviously slow and the cumulative infiltration showed a significant transition. The relationship between cumulative infiltration and time changed from a nonlinear process to a linear process, indicating that the infiltration process became stable. The test data after the transition was fitted by a linear relationship, yielding the following formula:

$$I=0.022\times t+3.427(100<t\le 380)R^2=0.9997;$$

(6)

$$I=0.0145\times t+6.295\left(380<t\le 575\right)R^2=0.9988,$$

(7)

where I is the cumulative infiltration amount (cm) and t is the infiltration time (min).

The relationship between the infiltration rate and time can be obtained from the derivation of the relationship between the cumulative infiltration amount and time over time. The relationship between soil infiltration rate and time in the upper part of the layered soil can be expressed by the following formula:

$$i=0.144\times t^{-0.2664},$$

(8)

where $i({\mathrm{cm}\min }^{-1})$ is the infiltration rate.

The variation in soil infiltration rate with time is shown in Figure 9b. The infiltration rate of the layered soil was the same as that of the homogeneous soil before the wetting front reached the lower layer. After the wetting front entered the lower layer, the infiltration rate became constant, that is, the infiltration turned into a stable process. There were two different linear parts, 0.0221 cm/min and 0.0145 cm/min, meaning that the infiltration rate reduced again at t = 380 min when the wetting front reached the fifth sensor, for the reason that the palaeosol under the fifth sensor had many calcareous concretions and stronger water resistance. According to the variation characteristics of the infiltration rate, the whole infiltration process of the layered soils can be divided into nonlinear and linear stages. Moreover, the infiltration rates of stable states were significantly smaller than the instantaneous infiltration rate of homogeneous soil at the transition time $t_p$. This phenomenon fully illustrates the anti-seepage effect of layered soil, which is consistent with [18,22].

According to Equation (2), we can describe the cumulative infiltration as follows:

$$I=0.022\times t+3.427\times \left(1-e^{-\frac{t}{33}}\right)\left(0<t\le 380\right);$$

(9)

$$I=0.0145\times t+6.295\times \left(1-e^{-\frac{t}{33}}\right)(380<t\le 575).$$

(10)

It can be seen from Figure 10 that both the calculated values and the measured values are in good agreement. The layered soil infiltration model established by superposition principle can simulate the layered soil infiltration well. That is to say, the cumulative infiltration could be divided into the flushing process and the stable infiltration process before the wetting front arrived at the interface between the two soils. The cumulative infiltration became stable, with an infiltration rate smaller than that of the homogeneous soil under the same situation.

3.3. The Basic Physical and Hydraulic Properties

3.3.1. Laboratory SWCC Experiments

As shown in Figure 11a,b, the values of the correlation coefficient ($R^2$) of the fitting were 0.999 for the loess and 0.998 for the palaeosol, indicating that the laboratory results agree well with the VG model. Figure 11c indicates that the SWCCs of the loess and palaeosol intersect with each other at a suction value of approximately 155 cm, meaning that the water holding capacity of the loess layer was stronger than that of the palaeosol layer in the lower matric suction range and was weaker than that of the palaeosol layer when the matric suction was greater than 155 cm. We can see from Figure 11d that the unsaturated hydraulic conductivity of the loess and palaeosol was very low when volumetric water content was less than 0.3. The unsaturated hydraulic conductivity rapidly increased to the saturated hydraulic conductivity with the increase in volumetric water content when the volumetric water content was greater than 0.3. Furthermore, the unsaturated hydraulic conductivity of the loess was always higher than that of the palaeosol under the same volumetric water content.

Table 1 shows that the content of silt was higher than that of clay for both soils, which is consistent with the results of previous studies [22]. However, the loess had a higher quantity of coarse silt and a lower quantity of fine silt and clay than the palaeosol. The main reason is the difference in the pedogenesis due to different climates. That is, palaeosol developed in warmer/wetter climatic conditions was subjected to stronger pedogenesis than the loess, which was deposited during colder/drier times. As a result, the loess layer had a stronger water holding capacity than the palaeosol layer had in the lower suction section, and the water holding capacity of the loess was less than the palaeosol layer had in the higher suction section, as shown in Figure 11c. Moreover, water adsorption by clay also explains why the water holding capacity of the palaeosol was stronger than that of the loess in the higher suction section.

Combining the results of Figure 11c,d, we can conclude that when the water enters the palaeosol from the loess, the suction of the palaeosol is greater than that of the loess and the wetting front directly enters the palaeosol. That is to say, the infiltration of the lower palaeosol layer is less affected by the infiltration of the upper loess layer. The water-blocking effect of palaeosol is mainly due to the fact that the hydraulic conductivity of palaeosol is lower than that of loess. After the infiltrated water reaches the lower soil, the water supply from the upper layer soil is greater than the soil hydraulic conductivity of the lower layer soil and the infiltration is under the condition of sufficient water for the lower layer soil. Therefore, excess water may accumulate above the interface between the upper and lower soil layers and gradually generate temporary groundwater levels, which is consistent with [22].

3.3.2. Inverse Simulation Results of the DREAM Algorithm

The posterior probability density distributions of the soil hydraulic parameters were estimated using the Bayesian simulation with the DREAM(zs) algorithm. The parameters’ convergence was tested by the multi-chain $\hat{\mathrm{R}}$-statistic method mentioned above. As shown in Figure 12, the posterior probability distributions are skewed and the peak values are obvious, indicating that the parameter identification is good and the uncertainty is relatively small. The posterior probability distributions are irregular and some of the highest are the maximum or minimum values at the boundary of the parameter ranges. The parameter values derived from the fitting results are near the ranges of the posterior distribution. However, differences remain between them. The main reasons are the hysteretic effect and the size effect. Numerical models neglecting these differences may compromise the reliability of the predictions and misrepresent the responses. However, the fitting parameters can provide a reference for a decision of initial prior ranges in the MCMC optimization algorithm. In addition, the distinct multiple frontal distribution of the posterior distribution of some parameters means that there may be more than one optimal solution. Therefore, the MCMC algorithm can effectively avoid the problem of obtaining a local optimal solution.

4. Conclusions

Large soil column infiltration experiments under a constant water head were conducted on both laboratory and field scales to simulate the water infiltration process under agricultural irrigation. The hydraulic property differences between the Q3 loess and palaeosol were analysed by comparing the basic hydraulic parameters. The main conclusions of this paper are:

• According to the comparison of SWCC between loess and palaeosol, the water holding capacity of palaeosol is stronger than that of loess under the same matrix suction. Therefore, the water content of palaeosol is always higher than that of loess at the beginning and ending of the field test. The loess does not weaken the infiltration of the palaeosol. The reason for the water-resistance of the palaeosol is that its unsaturated hydraulic conductivity is lower than that of loess, which also leads to the appearance of a temporary water table at the interface between loess and palaeosol under certain infiltration conditions.
• In the laboratory infiltration test, before the wetting front reached the interface between loess and palaeosol, the relationship between the cumulative infiltration and time was nonlinear, which is consistent with that of homogeneous soil. However, after the wetting front reaches the interface, the relationship between the cumulative infiltration and time became linear. The superposition model based on the capacitance law can continuously simulate the nonlinear and linear relationships between the cumulative infiltration and time well.
• The water infiltration simulation experiments with two different scales illustrated that the process of water content change in the layered loess–palaeosol sequence can be divided into two stages: the first stage is the free infiltration stage of the loess under the constant water head; the second stage is the free infiltration stage of the palaeosol under the variable water head.