Inﬂuence of Gully Land Consolidation on Phreatic Water Transformation in the Loess Hilly and Gully Region

: Gully Land Consolidation (GLC) is a proven method to create farmlands and increase crop yields in the Loess Hilly and Gully Region, China. However, GLC inﬂuences phreatic water transformation and might cause the farmlands water disasters, such as salinization and swamping. For exploring the inﬂuence of GLC on phreatic water transformation and mitigating disasters, a series of indoor experiments were conducted in the artiﬁcial rainfall hall. Then, we simulated the phreatic water transformation patterns under more conditions with HYDRUS-3D. Finally, an engineering demonstration in the ﬁeld was performed to validate our research. The indoor experiments indicated that GLC could increase phreatic water outﬂow rate 4.39 times and phreatic water coefﬁcient (PWC) 2.86 times with a considerable delay. After calibration and validation with experimental data, the HYDRUS-3D was used to simulate phreatic water transformation under more soil thickness and rainfall intensities. Accordingly, we summarized the relationship among PWC, rainfall intensities, and soil thickness, and therefore suggested a blind ditch system to alleviate farmlands disasters. Field application showed that a blind ditch system could avoid disasters with 3.2 times the phreatic water transformation rate compared to loess. Our research provides implications for sustainable land uses and management in the region with thick soil covers.


Introduction
The growing population and developing urbanization have increased the demands for quality arable land around the world, which even causes food security, especially in developing countries [1,2]. Taking China as an example, the proportion of urban areas' population has increased from 30.4% to 52.6% from 1998 to 2012, and total farmland reduced by 7.93 million ha simultaneously [3]. Using 7.63% of the world's arable land to feed China's 19.78% of the world's population is one of China's top priorities [4]. The modern land consolidation project, which originated in Germany [5], has gained global attention as a pivotal role in expanding arable land and improving agricultural productivity [6,7].

Research Area
The simulated gully in our indoor experiment is located in Yangjuangou catchment of Baota District in Yan 'an City (109°31'17.91"E, 36°41 '48.31"N), which is in the central part of the Loess Plateau and belongs to the northern Shaanxi Loess Hilly and Gully region (as shown in Figure 1). The Yangjuangou catchment covers an area of 2.02 km 2 , with elevation ranging from 1018 to 1268 m. This catchment has the typical climate and geomorphological characteristics of the Loess Hilly and Gully region. The climate is semiarid, the geomorphic type is dominated by loess beam, and loess gully with the valley's degree is 2.74 km/km 2 . The annual mean temperature is 8.8 °C in the region, with the mean lowest temperature of −6.9 °C and the mean highest temperature of 22.6 °C. The mean annual precipitation is about 500 mm, with 70% falling from July to September [39,40]. The undisturbed gully (UG) in this region comprises a loess layer, a weathered fissure layer, and a bedrock layer from top to bottom. The loess layer's thickness is 3-15 m with the bulk density ranging from 1.3 to 1.5 g/cm 3 , and the permeability coefficient is 0.2-0.5 m/d. The thickness of the weathered fissure layer is about 20 m, and the permeability coefficient is 20-80 m/d. The bedrock layer's permeability coefficient is about 0.03-0.22 m/d, which can be considered impervious compared with the upper two layers. In the GLC area, the thickness of surface cushion loess is generally 2.2 m higher than the gully phreatic water [21].
Simultaneously, the land consolidation project in the Yangjuangou catchment is one of the subprojects of GLC from 2013, with a construction scale of 27.84 hectares. The GLC aims to build a platform by stages as the gully farmlands. At the end of 2020, a series of engineering measures, including land leveling engineering, irrigation, and drainage engineering, significantly improved the local farmland production capacity, farmland quality, and land-use efficiency in the region. Selecting the gully in such a typical catchment as a research area can show the GLC's effect on the phreatic water transformation and related disasters prevention in the Loess Hilly and Gully Region.

Experimental Conditions and Equipment
The indoor simulated rainfall experiments were carried out in the Artificial Rainfall Simulation Hall of the State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau from June to October of 2018 (as shown in Figure 2b). With an automatic The undisturbed gully (UG) in this region comprises a loess layer, a weathered fissure layer, and a bedrock layer from top to bottom. The loess layer's thickness is 3-15 m with the bulk density ranging from 1.3 to 1.5 g/cm 3 , and the permeability coefficient is 0.2-0.5 m/d. The thickness of the weathered fissure layer is about 20 m, and the permeability coefficient is 20-80 m/d. The bedrock layer's permeability coefficient is about 0.03-0.22 m/d, which can be considered impervious compared with the upper two layers. In the GLC area, the thickness of surface cushion loess is generally 2.2 m higher than the gully phreatic water [21].
Simultaneously, the land consolidation project in the Yangjuangou catchment is one of the subprojects of GLC from 2013, with a construction scale of 27.84 hectares. The GLC aims to build a platform by stages as the gully farmlands. At the end of 2020, a series of engineering measures, including land leveling engineering, irrigation, and drainage engineering, significantly improved the local farmland production capacity, farmland quality, and land-use efficiency in the region. Selecting the gully in such a typical catchment as a research area can show the GLC's effect on the phreatic water transformation and related disasters prevention in the Loess Hilly and Gully Region.

Experimental Conditions and Equipment
The indoor simulated rainfall experiments were carried out in the Artificial Rainfall Simulation Hall of the State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau from June to October of 2018 (as shown in Figure 2b). With an automatic simulation device of under sprinklers, the simulated rainfall system could ensure the kinetic energy of simulated precipitation close to the natural rainfall for the mean fall-height of 18 m [41]. m, at the Rainfall Simulation Hall. The slope of the model could be adjusted manually from 0° to 35°. There were two water tanks, 0.15 m × 1 m × 1 m, in front and back of the model flume, to regulate the underflow level. Above the front side of the water tank, there was one surface water groove, and there was one drainage pipe of groundwater of 0.4 m high at the front side of the water tank (as shown in Supplementary Figure S2b). Moreover, two neutron probes access tubes down to a depth of 0.9 m in experimental flume for soil moisture control (as shown in Supplementary Figure S2a and Figure 2).

Gully Generalization and Artificial Rainfall Design
According to field survey results and sand-box model size, we generalized the gully into a straight shape, and the experimental flume was fixed at an angle of 3° in this research. The impervious steel bottom of the flume indicated the impervious bedrock layer of the gully. The air-dried and sieved coarse sand in the lower layer was 0.9 m thick (0.25 mm ≤ grain size ≤ 0.5 mm) with a bulk density of 1.5 g/cm 3 , which simulated the fissured layer for phreatic water transportation in the gully [44]. The top layer of the filling medium in the flume was 0.1 m local loess with a bulk density of 1.28 g/cm 3 and a water content of 13%, which simulated the gully's surface soil (grading curve is shown in Figure  3). The initial phreatic water level was parallel to the flume bottom, and the distance between the two lines was 0.4 m. The experiments were conducted in the sand-box model [42,43], 5.3 m × 1 m × 1 m, at the Rainfall Simulation Hall. The slope of the model could be adjusted manually from 0 • to 35 • . There were two water tanks, 0.15 m × 1 m × 1 m, in front and back of the model flume, to regulate the underflow level. Above the front side of the water tank, there was one surface water groove, and there was one drainage pipe of groundwater of 0.4 m high at the front side of the water tank (as shown in Supplementary Figure S2b). Moreover, two neutron probes access tubes down to a depth of 0.9 m in experimental flume for soil moisture control (as shown in Supplementary Figure S2a and Figure 2).

Gully Generalization and Artificial Rainfall Design
According to field survey results and sand-box model size, we generalized the gully into a straight shape, and the experimental flume was fixed at an angle of 3 • in this research. The impervious steel bottom of the flume indicated the impervious bedrock layer of the gully. The air-dried and sieved coarse sand in the lower layer was 0.9 m thick (0.25 mm ≤ grain size ≤ 0.5 mm) with a bulk density of 1.5 g/cm 3 , which simulated the fissured layer for phreatic water transportation in the gully [44]. The top layer of the filling medium in the flume was 0.1 m local loess with a bulk density of 1.28 g/cm 3 and a water content of 13%, which simulated the gully's surface soil (grading curve is shown in Figure 3). The initial phreatic water level was parallel to the flume bottom, and the distance between the two lines was 0.4 m.
The GLC dam could ensure that gully farmlands would remain stable in rainstorms [20]. So, the dam system in GLC was generalized as an aluminum-plastic plate with 10 mm thickness and 0.05 m height. The upstream area of the dam in the flume was simulated gully farmlands with a length of 1.5 m, soil height of 0.1 m, and slope gradient of 0 • . The entire simulated gully had two generalized dams, and the control ratio was 60% (as shown in Figure 2d). Before we performed the simulated gully surface with loess, three performance tests with pure coarse sand had been conducted to ensure the accuracy of the amount and flow field of phreatic water in the flume (as shown in Figure S2 in Supplementary Materials). As for the subsequent experiments with loess covered, they could be regarded as controlled experiments. The GLC dam could ensure that gully farmlands would remain stable in rainstorms [20]. So, the dam system in GLC was generalized as an aluminum-plastic plate with 10 mm thickness and 0.05 m height. The upstream area of the dam in the flume was simulated gully farmlands with a length of 1.5 m, soil height of 0.1 m, and slope gradient of 0°. The entire simulated gully had two generalized dams, and the control ratio was 60% (as shown in Figure 2d). Before we performed the simulated gully surface with loess, three performance tests with pure coarse sand had been conducted to ensure the accuracy of the amount and flow field of phreatic water in the flume (as shown in Figure S2 in Supplementary Materials). As for the subsequent experiments with loess covered, they could be regarded as controlled experiments.
Considering the influences of rainfall and duration on the GLC project, we designed 120 mm as the total artificial rainfall in this research, and three constant rainfall intensities were set as 60, 90, and 120 mm/h, with rainfall occurrence periods of 1.5, 2, and 5 years, respectively [45]. The size distribution of raindrops is an important factor in affecting rainfall infiltration, and the stain method was used to measure raindrop size in this research [46]. We spread the dye powder on the filter paper thinly and evenly, and the filter paper was pale pink under dry conditions. We put the filter paper in the rain, and then there was an approximately circular red stain formed where the raindrops hit (as shown in Figure S3 in Supplementary Materials). After measuring the average diameter of stains, we used the relationship between stain diameter and raindrop diameter to calculate the raindrop diameter [47]. Furthermore, after our measurement and statistics, the rainfall information of this research is shown in Table 1, and the calibration tests showed more than 85% uniformity of the rainfall.  Considering the influences of rainfall and duration on the GLC project, we designed 120 mm as the total artificial rainfall in this research, and three constant rainfall intensities were set as 60, 90, and 120 mm/h, with rainfall occurrence periods of 1.5, 2, and 5 years, respectively [45]. The size distribution of raindrops is an important factor in affecting rainfall infiltration, and the stain method was used to measure raindrop size in this research [46]. We spread the dye powder on the filter paper thinly and evenly, and the filter paper was pale pink under dry conditions. We put the filter paper in the rain, and then there was an approximately circular red stain formed where the raindrops hit (as shown in Figure S3 in Supplementary Materials). After measuring the average diameter of stains, we used the relationship between stain diameter and raindrop diameter to calculate the raindrop diameter [47]. Furthermore, after our measurement and statistics, the rainfall information of this research is shown in Table 1, and the calibration tests showed more than 85% uniformity of the rainfall.

Monitoring Methods and Data Collection
The experiments started after adjusting the initial phreatic water level up to 0.4 m. During this research, the main monitoring items included surface runoff amount, surface runoff in the process, phreatic water amount, and outflow process of phreatic water. Surface runoff was collected by collecting buckets at the beginning of runoff yield. Collecting buckets were used to collect surface runoff at the beginning of runoff yield for 30 s at 5 min intervals during experiments. Besides, phreatic water was collected by measuring cylinder per 2000 mL uninterrupted for about 24 h. According to the phreatic water condition, the monitoring interval was lengthened gradually until the phreatic water level dropped to 0.4 m, the initial control phreatic water level.

Water Transformation Model
The US Salinity Laboratory has developed HYDRUS to simulate water and solute transportation in variable saturated porous media [48,49]. HYDRUS-3D is a finite element computer model that can simulate the three-dimensional water flow movement in variable saturation porous media. The soil moisture module used by HYDRUS-3D follows Darcy's law and the continuity principle of mass conservation. It could reasonably simulate the transportation regular pattern of water and solute in soil [50][51][52]. Therefore, we chose the HYDRUS-3D to explore the influence of rainfall transformation to phreatic water under more conditions, such as more loess thickness and rainfall intensities with or without GLC.
According to the indoor experiment and calculation principle of HYDRUS-3D, the loess landfill layer with same size and slope had a tetrahedral grid (as shown in Figure 4). The calculation time unit of the model was minutes (min), and the calculation period was set according to the observation time of the indoor rainfall experiments. The impervious boundary was set around the soil layer. The rainfall was generalized as atmospheric boundary. The bottom soil layer was set as seepage boundary, and the dam measure was set as impervious boundary. Thus, the cumulative flow of seepage boundary in the simulation was accounted as a total amount of phreatic water. Initial water content and other parameters were set according to experiment and soil parameters.
The experiments started after adjusting the initial phreatic water level up to 0.4 m. During this research, the main monitoring items included surface runoff amount, surface runoff in the process, phreatic water amount, and outflow process of phreatic water. Surface runoff was collected by collecting buckets at the beginning of runoff yield. Collecting buckets were used to collect surface runoff at the beginning of runoff yield for 30 seconds at 5 min intervals during experiments. Besides, phreatic water was collected by measuring cylinder per 2000 ml uninterrupted for about 24 hours. According to the phreatic water condition, the monitoring interval was lengthened gradually until the phreatic water level dropped to 0.4 m, the initial control phreatic water level.

Water Transformation Model
The US Salinity Laboratory has developed HYDRUS to simulate water and solute transportation in variable saturated porous media [48,49]. HYDRUS-3D is a finite element computer model that can simulate the three-dimensional water flow movement in variable saturation porous media. The soil moisture module used by HYDRUS-3D follows Darcy's law and the continuity principle of mass conservation. It could reasonably simulate the transportation regular pattern of water and solute in soil [50][51][52]. Therefore, we chose the HYDRUS-3D to explore the influence of rainfall transformation to phreatic water under more conditions, such as more loess thickness and rainfall intensities with or without GLC.
According to the indoor experiment and calculation principle of HYDRUS-3D, the loess landfill layer with same size and slope had a tetrahedral grid (as shown in Figure 4). The calculation time unit of the model was minutes (min), and the calculation period was set according to the observation time of the indoor rainfall experiments. The impervious boundary was set around the soil layer. The rainfall was generalized as atmospheric boundary. The bottom soil layer was set as seepage boundary, and the dam measure was set as impervious boundary. Thus, the cumulative flow of seepage boundary in the simulation was accounted as a total amount of phreatic water. Initial water content and other parameters were set according to experiment and soil parameters.

Influence of Gully Land Consolidation on the Phreatic Water Based on Indoor Experiments
The implementation of GLC has controlled soil erosion and increased farmlands area over the Loess Plateau. However, various disasters have been reported in the gully farmlands created by the GLC, such as swamping and salinization [31,33]. These disasters have negatively affected local agriculture. To explore the main reasons causing these disasters,

Influence of Gully Land Consolidation on the Phreatic Water Based on Indoor Experiments
The implementation of GLC has controlled soil erosion and increased farmlands area over the Loess Plateau. However, various disasters have been reported in the gully farmlands created by the GLC, such as swamping and salinization [31,33]. These disasters have negatively affected local agriculture. To explore the main reasons causing these disasters, we first conducted indoor experiments in the artificial rainfall simulation hall. We attempted to reveal the influence of GLC on the phreatic water transformation under different gully conditions. The relevant results are shown in the following.
As shown in Figure 5, the phreatic water outflow processes with or without land consolidation treatment are entirely different. The hydrographs are similar under different rainfall intensities with the same gully treatments. Without land consolidation, the phreatic water hydrographs have weak responses to rainfall intensities, and the outflow rate continues to decrease with small fluctuations. Moreover, the magnitude of phreatic water is relatively small. The average phreatic water runoff under three rainfall intensities is 1.85 × 10 −5 m 3 /s, 1.82 × 10 −5 m 3 /s, and 1.67 × 10 −5 m 3 /s. In contrast, the phreatic water hydrographs with GLC have large peaks 6.2 h after the rainfall, and the outflow rate of phreatic water in GLC is 5.39 times that in UG. The average phreatic water runoff under three rainfall intensities is 9.91 × 10 −5 m 3 /s, 9.67 × 10 −5 m 3 /s, and 9.21 × 10 −5 m 3 /s. Meanwhile, we analyzed ratios among three water balance components, including surface runoff, soil water, and phreatic water (see Figure S4 in Supplementary Materials). Obviously, the GLC can adjust the average ratios from 79:9:12 to 40:14:46. The GLC effectively increases the proportion of precipitation transformed to soil water and phreatic water from 21% to 60%. In a word, GLC influenced the interaction mechanism between precipitation and phreatic water, and there was more phreatic water transformed with a long-prolonged time.
consolidation treatment are entirely different. The hydrographs are similar under different rainfall intensities with the same gully treatments. Without land consolidation, the phreatic water hydrographs have weak responses to rainfall intensities, and the outflow rate continues to decrease with small fluctuations. Moreover, the magnitude of phreatic water is relatively small. The average phreatic water runoff under three rainfall intensities is 1.85 × 10 −5 m 3 /s, 1.82 × 10 −5 m 3 /s, and 1.67 × 10 −5 m 3 /s. In contrast, the phreatic water hydrographs with GLC have large peaks 6.2 hours after the rainfall, and the outflow rate of phreatic water in GLC is 5.39 times that in UG. The average phreatic water runoff under three rainfall intensities is 9.91 × 10 −5 m 3 /s, 9.67 × 10 −5 m 3 /s, and 9.21 × 10 −5 m 3 /s. Meanwhile, we analyzed ratios among three water balance components, including surface runoff, soil water, and phreatic water (see Figure S4 in Supplementary Materials). Obviously, the GLC can adjust the average ratios from 79:9:12 to 40:14:46. The GLC effectively increases the proportion of precipitation transformed to soil water and phreatic water from 21% to 60%. In a word, GLC influenced the interaction mechanism between precipitation and phreatic water, and there was more phreatic water transformed with a long-prolonged time. Although the impact of GLC on the process of phreatic water transformation is demonstrated in Figure 5, it cannot directly express the influence of different rainfall intensities on the transformation quantity from precipitation to phreatic water. The runoff coefficient is widely used to describe the precipitation transformed to runoff [53][54][55]. Similarly, we introduced the phreatic water coefficient (PWC) to better evaluate the ability to transform rainfall into phreatic water under different rainfall intensities. The PWC can be expressed as follows: Where PWC is the phreatic water coefficient (dimensionless), is the phreatic water volume (m³), is the rainfall (mm), and is the catchment area (m 2 ). As shown in Table  2, with increasing rainfall intensity, PWC is decreasing regardless of the treatments. However, the average PWC of GLC is 0.46, which is 3.86 times that of UG. What's more, the range of PWC in GLC is smaller than UG. It indicates that the GLC has dramatically improved the phreatic water transformation and is less affected by the rainfall intensities. Although the impact of GLC on the process of phreatic water transformation is demonstrated in Figure 5, it cannot directly express the influence of different rainfall intensities on the transformation quantity from precipitation to phreatic water. The runoff coefficient is widely used to describe the precipitation transformed to runoff [53][54][55]. Similarly, we introduced the phreatic water coefficient (PWC) to better evaluate the ability to transform rainfall into phreatic water under different rainfall intensities. The PWC can be expressed as follows: where PWC is the phreatic water coefficient (dimensionless), Q re is the phreatic water volume (m 3 ), P is the rainfall (mm), and A is the catchment area (m 2 ). As shown in Table 2, with increasing rainfall intensity, PWC is decreasing regardless of the treatments. However, the average PWC of GLC is 0.46, which is 3.86 times that of UG. What's more, the range of PWC in GLC is smaller than UG. It indicates that the GLC has dramatically improved the phreatic water transformation and is less affected by the rainfall intensities. Through the above analysis, we can conclude that GLC impacts the process and quantity of phreatic water transformation. In the UG treatments, most precipitation directly flows out of the gully in surface runoff. As for GLC treatments, the dam of GLC impounds surface runoff firstly. Then, the impounded water would infiltrate into soil water and phreatic water after a long period (as shown in Figures 2d and 5). The more surface runoff impounded, the more phreatic water would be transformed under ideal circumstances. However, the GLC thickens the loess soil in the gully, and the frequently heavy rain accelerates the physical crust formation process of loess. The time of impounded runoff transformed to phreatic water would be longer. The long-term accumulation of impound surface runoff may cause farmland disasters. Therefore, it is necessary to explore the influence of different thickness loess soil on the phreatic water transformation for farmlands' disaster prevention.

Simulation Based on HYDRUS-3D
Based on the indoor experiments and analysis above, we know that the GLC significantly impacts gully phreatic water transformation, i.e., increase the transformation volume and extend the transformation time. The GLC creates farmlands by cutting the slope to fill the gully. It brings thick filling loess in the gully farmlands, which is why transformation time extends [33,56]. Due to the indoor experiments' limitations of time and workforce, it is impracticable to explore the influence of GLC on phreatic water transformation under various loess thicknesses. The HYDRUS-3D provides a way to reasonably simulate water and solute transportation regular patterns in soil with high efficiency. Therefore, we used HYDRUS-3D to explore the effects of different loess thickness on phreatic water transformation under more rainfall intensities.

Calibration and Validation
In HYDRUS-3D, the van Genuchten model and the Brooks-Corey model are commonly used to model soil moisture characteristics. The Brooks-Corey model is better for homogeneous and isotropic coarse-textured soil with narrow pore size distribution but has worse accuracy for fine-textured soils. The van Genuchten model, on the other hand, is continuous and universal for most soils in a wide range of moisture content, and its line shape fits well with the measured data curve. Simultaneously, it can make the suction force of saturated soil zero, which is in line with the characteristics of soil suction force change during the process of moisture absorption. Considering the soil's nature for this research, which is high clay content [39], therefore, the van Genuchten model is more suitable than another model [57,58]. The expression is: where S e is the relative saturation (dimensionless), θ s is the saturated water content of the soil (cm 3 /cm 3 ), θ r is the retained water content of the soil (cm 3 /cm 3 ), h is the soil negative pressure head (cm), K h is the unsaturated hydraulic conductivity (cm/min), K s is the saturated hydraulic conductivity (cm/min), a, m, n are fitting parameters, which are related to the physical properties of the soil, where m = 1 − 1 n , n > 1, and l is the pore connectivity parameter, which generally is 0.5.
Model calibration is the process whereby selected model input parameters are adjusted within reasonable limits to produce simulation results that best match the known or measured values. It is the most critical process in building a model because the calibration and validation quality inevitably determines the reliability of any conclusions and recommendations made using the simulation results. Thus, we chose the accumulation process of phreatic water in indoor UG treatments under the rainfall intensity of 60 mm/h as the calibrated data to obtain the van Genuchten model's final parameters. As for default parameters, we can determine them using the Rosetta neural network prediction module of HYDRUS-3D. The final parameters are shown in Table 3. After determining the final parameters by calibration, we took the indoor experiment's data with GLC treatments under the same rainfall intensity (60 mm/h) to validate the model. As shown in Figure 6, there are the results of calibration and verification. Figure 6a draws the accumulation process by observing and calculating data. Analyzing the accumulation process from experiments beginning at the phreatic water outflow stop under UG and GLC treatments. We can conclude that the model calculation curves fit well with the indoor experiment's phreatic water accumulative curve. The phreatic water accumulation of UG increases at a slow rate, which is 1.82 × 10 −5 m 3 /s on average and coincides with the above. The calculated accumulation of GLC is like the indoor experimental process. The phreatic water transformation rate was slow before rainfall stops (0-120 min). Then, the rate increased to a maximum of 1.81 × 10 −4 m 3 /s at 376 min. Later, the rate decreased until most surface impounds water transformed to phreatic water (376-647 min). Finally, the remaining surface water transformed to phreatic water at a similar rate as UG, which is 1.98 × 10 −5 m 3 /s. Additionally, we use the mean absolute percentage error (MAPE) to evaluate model accuracy. The MAPE of the calibration period is 8.03%, and the validation period is 12.08%. Figure 6b compares the calculated and observed quantities of phreatic water under all experimental conditions. The MAPE of the observed value and calculation is 7.21%. Therefore, the simulation of precipitation transformed to phreatic water performed well, and we consider that HYDRUS-3D can meet the simulation accuracy of phreatic water transformation in this research.

Model Simulation of Phreatic Water Transformation under More Conditions
After calibration and verification, the accuracy of the phreatic water transformation model constructed by HYDRUS-3D shows that we can use this model to simulate more

Model Simulation of Phreatic Water Transformation under More Conditions
After calibration and verification, the accuracy of the phreatic water transformation model constructed by HYDRUS-3D shows that we can use this model to simulate more conditions, such as more different soil thicknesses and rainfall intensities. Therefore, we constructed GLC and UG gully models with treatments of filling 0.1-1 m thickness engineering loess and adding two more rainfall intensities that have not been tested indoors (30 and 150 mm/h). Then, we count the phreatic water transformation in various conditions, as follows. Figure 7 shows the relationship between PWC under more loess thickness and rainfall intensities, including thickness range from 0.1 to 1 m and rainfall intensities from 30 to 150 mm/h. In this graph with a double logarithmic coordinate system, we can summarize the following laws. The average PWC of GLC is 4.12 times that of UG. Furthermore, the phreatic water transformation curve's changing pattern is similar under the same gully condition regardless of the rainfall intensities. The rainfall intensities determine how much the phreatic water quantities are negatively correlated. Thus, the quantity of rainfall transformed to phreatic water in the UG is more than GLC. For example, the phreatic water transformation of UG with 30 mm/h is a smaller portion more than that in GLC under 150 mm/h rainfall intensity, with the loess thickness increasing from 0.2 to 0.4 m.
Meanwhile, all curves have a characteristic: the PWC decreases quickly when the soil thickness exceeds 0.4 m, and the average decrease rate per meter of 0.4-0.5 m is 139.3 times that of 0.1-0.4 m. Moreover, the rate of the decrease per meter of GLC with loess thickness of 0.1-0.4 m is 3.87 times that of UG. It also shows that the accumulated runoff caused by GLC further inhibits the phreatic water transformation [40]. The area with thick soil cover and heavy rain is prone to disasters caused by more surface runoff accumulation and the tardy phreatic water transformation, especially in the Loess Plateau, where the loess has low water permeability and is easy to form the physical crust to impede the infiltration process.  According to Figure 7, the changing pattern of curves is similar, and to better summarize the phreatic water transformation laws, we can describe this pattern with such equation: where is the phreatic water coefficient (dimensionless), ℎ is the thickness of the loess soil (m), and and are parameters related to rainfall intensities. Then, summarizing the parameters in the Equation (4), we can establish the correlation among parameters , , and (rain intensity) as a linear relationship, as shown in Figure 8, and the equations describing the linear relationship between the parameters , , and are According to Figure 7, the changing pattern of curves is similar, and to better summarize the phreatic water transformation laws, we can describe this pattern with such equation: where PWC is the phreatic water coefficient (dimensionless), h is the thickness of the loess soil (m), and a and b are parameters related to rainfall intensities. Then, summarizing the parameters in the Equation (4), we can establish the correlation among parameters a, b, and i (rain intensity) as a linear relationship, as shown in Figure 8, and the equations describing the linear relationship between the parameters a, b, and i are shown as Equations (5)- (8): Next, we incorporate Equations (5)-(8) into (4), respectively. Equations (9) and (10) can be obtained as show in the following: PWC GLC = 4.24 i 1.18 h 0.0035i+0.96 (10) By comparing and analyzing the Equations (9) and (10), we can conclude that the PWC of GLC is more sensitive to loess thickness. The average standard deviation of GLC is 0.17 under the same rain intensity and different loess thickness, which is 3.74 times that of UG. Although, when the loess thickness is the same, the PWC of UG are stable within a small range, while rainfall intensities increase. Too little phreatic water transformation means more surface runoff and soil erosion. Therefore, there is a need to reduce the impact of loess thickness on PWC of GLC to increase gully phreatic water and reduce erosion. Meanwhile, considering that the cultivated soil thickness in the Loess Hilly and Gully Region is mostly 0-0.2 m [59,60], the average rainfall intensity is within the range of 30-140 mm/h [61]. We conclude that GLC transforms rainfall to phreatic water more efficiently in gully agricultural land by calculating Equations (9) and (10). The average PWC of the GLC is 4.15 times that of the UG under the condition that impounds surface rainfall transforms to phreatic water smoothly. On the one hand, the GLC dams protect farmland from floods and increase the phreatic water transformation volume; on the other hand, more surface runoff accumulation caused by GLC brings some hidden dangers of disasters, especially in the area with thick loess and frequent heavy rain. Region is mostly 0-0.2 m [59,60], the average rainfall intensity is within the range of 30-140 mm/h [61]. We conclude that GLC transforms rainfall to phreatic water more efficiently in gully agricultural land by calculating Equations (9) and (10). The average PWC of the GLC is 4.15 times that of the UG under the condition that impounds surface rainfall transforms to phreatic water smoothly. On the one hand, the GLC dams protect farmland from floods and increase the phreatic water transformation volume; on the other hand, more surface runoff accumulation caused by GLC brings some hidden dangers of disasters, especially in the area with thick loess and frequent heavy rain. .

Field Application and Gully Farmlands' Disaster Prevention
It is concluded that the farmlands near the dam with thick soil are more likely to accumulate surface runoff by field investigation, indoor experiments, and model simula-

Field Application and Gully Farmlands' Disaster Prevention
It is concluded that the farmlands near the dam with thick soil are more likely to accumulate surface runoff by field investigation, indoor experiments, and model simulation. In response to this problem and better utilization of GLC for phreatic water transformation, this research simulated phreatic water transformation and drainage measures in a sub-basin in the Yangjuangou GLC area, Yan'an, Shaanxi Province.
There is a GLC dam that is 30 m in length and 5 m high at the exit of this sub-basin. The farmland created by GLC is 30 m wide and 65 m long, with a longitudinal gradient of 3 • , and the soil properties are the same as the loess in the experiments. The simulated phreatic water drainage and transformation measure is the underground blind ditch filled with gravel, which adopts the main and branch blind ditch. Three blind ditches are connected by the horizontal blind ditch in the upper consolidation dam. In HYDRUS-3D, it was set as a "combined cave" in the simulated soil layer, and the buried depth designed according to Equation (10) is 0.4-2 m and varies with the slope. The surrounding interface was generalized as free drainage, and its internal permeability coefficient was 9.7 cm/s, which simulates filling gravel (as shown in Figure S5a in Supplementary Materials). Based on the simulation, we demonstrated the blind ditch transformation and drainage engineering project in this sub-basin, with the same blind ditch layout as the simulation (as shown in Figure S5b in Supplementary Materials). A cellar was arranged at the blind ditch terminal to collect, monitor, and efficiently utilize the phreatic water.
Three days after rainfall of 54.2 mm in August 2019, we conducted a field investigation on this sub-basin. We measured the cellar's water collection, and about 1.87 × 10 −1 m 3 of phreatic water was collected. By inputting the rainfall into the forecasting model above, we calculated the accumulated flow of free infiltration, which is 1.63 × 10 −1 m 3 . Compared with the actual measurement, the MAPE is 13.1%. Figure 9 depicts the comparison of the phreatic water process between the blind ditch and cellar collecting and compares the observed and simulated water depth of the cellar. By comparing the data at observation time, the MAPE is 9.43%. Later, we counted the average infiltration rate at the junction between the blind ditch and loess ( Figure S6 in Supplementary Materials shows the water velocity vector in different simulation periods). It can conclude that the maximum infiltration rate in the loess is 2.88 × 10 −4 m/d, which remains unchanged on the first day after the rainfall stopped. In contrast, the maximum rate in the blind ditch comes on day 1.5, which is 2.1 times that of the loess layer. The maximum drainage transformation efficiency of the blind ditch is 3.2 times that of the loess layer (as shown in Figure S7 in Supplementary Materials). Therefore, the blind ditch could effectively avoid swamping and salinization by increasing the infiltration rate to promote rainfall transformation to phreatic water. Compared with the actual measurement, the MAPE is 13.1%. Figure 9 depicts the comparison of the phreatic water process between the blind ditch and cellar collecting and compares the observed and simulated water depth of the cellar. By comparing the data at observation time, the MAPE is 9.43%. Later, we counted the average infiltration rate at the junction between the blind ditch and loess ( Figure S6 in Supplementary Materials shows the water velocity vector in different simulation periods). It can conclude that the maximum infiltration rate in the loess is 2.88 × 10 −4 m/d, which remains unchanged on the first day after the rainfall stopped. In contrast, the maximum rate in the blind ditch comes on day 1.5, which is 2.1 times that of the loess layer. The maximum drainage transformation efficiency of the blind ditch is 3.2 times that of the loess layer (as shown in Figure  S7 in Supplementary Materials). Therefore, the blind ditch could effectively avoid swamping and salinization by increasing the infiltration rate to promote rainfall transformation to phreatic water.

Impact of GLC on Gully Surface Runoff and Phreatic Water
This research strongly suggested that the GLC played an essential role in transforming rainfall to gully phreatic water. We confirm that as a soil and water conservation measure, the dam of GLC impounds the rainfall and forms reservoirs of gully phreatic water, which alleviate water shortage and increase the area of quality farmland [56,62]. On the other hand, as a unique and essential invention designed for both food production and environment in the Loess Plateau, dam land has thick loess and higher water retention

Impact of GLC on Gully Surface Runoff and Phreatic Water
This research strongly suggested that the GLC played an essential role in transforming rainfall to gully phreatic water. We confirm that as a soil and water conservation mea-sure, the dam of GLC impounds the rainfall and forms reservoirs of gully phreatic water, which alleviate water shortage and increase the area of quality farmland [56,62]. On the other hand, as a unique and essential invention designed for both food production and environment in the Loess Plateau, dam land has thick loess and higher water retention capacity [63]. The gully surface runoff could be formed in a short time after rainfall [64] and the characteristics of loess (clayey and easily form physical soil crusts) generate surface water ponding [65][66][67]. It takes more time for surface ponding to be transformed into phreatic water. As shown in Figure 6a, there are three stages in transforming rainfall to phreatic water in the GLC area.
In the first stage, the rainfall runoff is intercepted by the dam to form surface ponding, and the wetting peak in the loess does not reach the weathered fissure layer. The phreatic water infiltration rate is mainly affected by the nature of the loess. Besides, the pressure of the surface ponding would further increase the transformation rate of phreatic water, so the initial phreatic water transformation rate of GLC is higher than that of UG. Moreover, due to the loess's higher water retention capacity, the wetting peak tends to take a long time to reach the weathered fissure layer, and this is why there are more delays in phreatic water transformation in the GLC area [68,69]. Once the wetting peak reaches the weathered fissure layer, the pathway for transforming surface ponded water into phreatic water is formed [70]. Then, the first stage of the transformation process ends, and the second phase begins. During the second stage, the transformation rate of phreatic water in GLC increased significantly. However, when most surface ponding has transformed to phreatic water, the hydrostatic pressure of surface ponding is not enough to increase the transformation rate, resulting in the transformation rate that is only influenced by the loess's nature. Thus, it is the third stage of the phreatic water transformation process. In this stage, a small amount of surface ponding and mostly saturated soil moisture reaches the weathered fissure layer in the form of infiltration. The transformation rate of phreatic water in GLC is just the same as that of UG. When the soil moisture becomes unsaturated, the phreatic water transformation comes to an end.
On the positive side, GLC reduces gully soil erosion by accumulating rainfall-runoff and ensures the water use of gully agriculture, especially under the special rainfall conditions of the Loess Plateau (heavy rainfall in a short period is common) [20,71,72]. However, according to the above analysis, the first stage significantly affects the phreatic water transformation rate. Due to excessive phreatic water transformation delays and surface ponding from rainfall accumulation quickly, the farmlands would be prone to suffering disasters such as swamping and dam breaks [73,74]. What's more, there is a high evaporation rate and surface ponding caused by transformation slowly after a short-termed rainfall. It is easy to occur salinization in the near-dam farmlands. Even if the transformation is timely, the water table's rise led by excessive phreatic water may cause the same problem [33,36]. Therefore, measures must be taken to prevent gully farmlands' water disasters by reducing the phreatic water transformation delay.

Prevention of Gully Water Disasters
In response to the water disaster in the GLC, it seems there are two ideas to solve these problems. Firstly, we can drain surface runoff by building drainage trenches in the farmlands and the drainage canals on both sides of the reshaped gully [33]. This method applies to farmlands downstream of the gully, because if drainage trenches are placed on farmlands upstream of the gully, there is more surface runoff flow to downstream lands, increasing the risk of water disasters in the downstream gully. Moreover, the upper stream's landform is more twisted than that of the downstream and does not favor water drainage.
Secondly, another concept is to accelerate the transformation of surface runoff to phreatic water, accelerating the rate of wetting peak reaching the weathered fissure layer. Our idea shortens the distance of wetting peak to weathered fracture layer by constructing an artificially weathered fissure zone close to the surface-a blind ditch system filled with gravel. Our research also confirms that it is an effective measure to solve excessive surface ponding in the GLC farmland. We analyzed the soil samples collected in the field. Compared with farmland without a blind ditch, the increase of soil moisture content was reduced by 46.81%. The maximum soil moisture content's absolute value was reduced by 1.53% and maximum soil conductivity decreased by 15.41 µs/cm. Overall, we should lay out the comprehensive measures reasonably to avoid gully water disasters and ensure the sustainable use of GLC farmlands, especially in the upstream gully [75].
What's more, the combination of the empirical Formulas (1), (9), and (10) can not only provide a theoretical foundation for gully water disaster prevention in GLC areas but also provide a reference for groundwater recharge management in thick soil areas around the world. Taking this field demonstration as an example, we calculated the initial design buried depth of the blind ditch by using the above empirical formula. Simultaneously, the collection cellar laying at the end of the blind ditch also serves as a water source for farmland, especially in the Loess Plateau [76,77]. Besides, the blind ditch system is an effective measure to lower the groundwater level [78], reducing the risk of dam breaks and swamping [79].

Methodological Limits and Future Research
The method of combining the experiment and the mathematical model had an absolute advantage in this research. It used fewer experimental data to calibrate and verify the selected model, HYDRUS-3D, and we used the calibrated model to simulate and predict the other treatments. Then, the model established by experiments was used to design and simulate the field demonstration project, which improves the efficiency of blind ditch design.
However, we only simulated water transformation in a gully with or without GLC, which is limited in exploring the more specific process of GLC on soil moisture, soil salinity, and phreatic water in the Loess Plateau. It is necessary to evaluate the effectiveness of the GLC from the scale of the area [80]. Moreover, we focused on disaster prevention in this research and did not consider other ecosystem aspects such as crop yield. At the same time, a research to compare the accumulation of biochar in farmland with or without disaster protection measures (blind ditch) is underway. Although the gully with GLC is buried with consolidation loess, the buried loess is not mixed with the gully's original loess or the bedrock at the gully's bottom. The phreatic water will flow along the cracks between different soil layers, and how to solve this problem is the direction of the next simulation.

Conclusions
In this research, we tried to explore the reasons for gully water disasters by simulating the impact of GLC on the transformation of rainfall to phreatic water under different conditions and compared it with the field demonstration project. Limited by soil texture and climate of the Loess Hilly and Gully Region, the GLC with dam impounds too much surface runoff. On the one hand, it reduces gully erosion and replenishes phreatic water; on the other hand, it increases farmland swamping and salinization probability. Our research proposed that the blind ditch system could reduce the probability of disasters and maintain gully farmlands' agricultural productivity. It helps sustainable land uses and management in regions around the world with thick soils, especially in the Loess Hilly and Gully Region with GLC.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 441/13/4/538/s1, Figure S1: Primary process of GLC. Figure S2: Pure coarse sand experiments. (a) Soil flume filled with pure coarse sand, (b) performance testing. Figure S3: Raindrop diameter measurement based on stain method under rainfall intensities of 120 mm/h. Figure S4: Water balance components ratio under different rainfall intensities. Figure S5: Field simulation and demonstration. (a) Simulation of GLC with the blind ditch in Yangjuangou sub-basin by HYDRUS-3D, (b) demonstration of the blind ditch in Yangjuangou sub-basin. Figure S6: Compensation of water velocity vector in different simulation periods. Figure S7: Flow rate at the junction of blind ditch and the loess.
Author Contributions: Z.G. (Zihao Guo) performed the conceptualization, data processing, and manuscript preparation; J.G. provided constructive suggestions to the revised manuscript; P.S. collected indoor experiment data and S.D. collected field demonstration data; P.S., J.L., X.L., and H.W. helped to design and conduct the indoor experiment; J.G. and S.D. were responsible for the design of the field demonstration project; R.A. contributed to the enhancement of the manuscript language; Z.G. (Zhe Gao) provided suggestions for the manuscript preparation. All authors have read and agreed to the published version of the manuscript.