Numerically Simulated Water Movement in Reclaimed MultiLayered Soil Backfilled with Yellow River Sediments

Author contributions All the authors contributed to the study conception and design. Sampling, soil column preparation, data collection and numerically simulated were performed mainly by Xiaotong Wang and Yusheng Liang. The first draft of the manuscript was written by Zhenqi Hu, Xiaotong Wang and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. Abstract The Yellow River interlayer filling reclamation technology can effectively improve the quality of destroyed cultivated land in the condition of limited soil resources. However, it is a conundrum to choose the appropriate sandwich strategy according to the amount of soil that can be backfilled. This study using the Hydrus-1D model to simulate water movement in reclaimed multiple-layered soils were to understand the mechanism of interlayer, and predict the optimal profile for reclaimed with Yellow River sediments. Simulations were performed on 18 soil profiles that were divided into a control check (CK) and two general scenarios that the total thickness of soil were 50 and 60 cm. Treatments in both scenarios exhibited interaction of different positions and thicknesses of soil interlayer. Results showed that removing part of the subsoil overlying the sediment placed it between sediment layers will improve the infiltration character of the conventional reconstructed soil profiles (T50-0 and T60-0). Moreover, changing the thickness of the interlayers affected infiltration character and soil water-holding capacity more than changes in the position of the layers for same total thickness of native soil. The optimal reconstructed soil profiles for scenarios 1 and 2 were T50-6 (interlayer thickness of 20 cm and located at a position of 30 cm) and T60-9 (interlayer thickness of 30 cm and located at a position of 30 cm), which could have a better infiltration character that were more closed to the native farmland.


Introduction
Underground coal mining has caused severe ecological and environmental problems, including land subsidence and the degradation of farmland . China possesses the largest coal mining subsidence area in the world (Wang et al. 2016), currently 700,000 km 2 and is continually increasing at a rate of 130 km 2 every year ). The reclamation of subsided farmlands with locally available unconsolidated materials plays a crucial role in the existing reclamation strategies. Typical fill materials include coal gangue and fly ash, both limited in their availability and having heavy metals (Tang et al. 2018). In contrast, sediments from the river present a relatively low risk of polluting the soils with heavy metals (Wang et al. 2016), and are readily available, making them a potentially useful fill material for land reclamation. In the united states, for example, Illinois River sediment was successfully used to reclaim a brownfield along the south Chicago lakefront (Darmody et al. 2004). Thus, reclaiming damaged land by using river sediments has broad potential for many countries around the world.
The Yellow River is the fifth longest river in the world and there exists large number of mining subsided located along its river banks. The sediment load in the Yellow River is the world's largest (Wang et al. 2016). However, Yellow River sediment is coarse-textured, with limited soil water-holding capacity . When it is combined with conventional methods of soil reconstruction, in which soils are directly spread over a backfilled layer of Yellow River sediment, the reclaimed soils possess poor soil water holding-capacities in comparison to native soils (Wang et al. 2014). Therefore, the effective use of Yellow River sediment requires the development of a reclamation strategy to enhance its water holding capacity.
Recent studies have suggested that layered soils consisting of materials of varying texture help maximize plant-available water in reclaimed soils (Huang et al. 2011;Huang et al. 2013;Zettl et al. 2015). Hu et al. (2017), for example, designed multi-layered soil profiles for the reclamation of subsidence land in eastern China. In this case, subsoil interlayers were emplaced into a layer of Yellow River sediment, which resulted in good productivity of maize growth (Hu et al. 2017). The authors developed a conceptual model for introducing the function of interlayers in Yellow river sediment ).
However, the mechanism of the interlayered soil profile, particularly with regards to hydrologic processes, has not been fully explored.
Infiltration is an important hydrologic properties of soil reclamation, because it is instrumental in partitioning precipitation into surface runoff and soil water storage (Doerr et al. 1996;Li et al. 2012). In multi-layered soils, both infiltration and the water holding capacity of the constructed profiles varied considerably in response to differences in interlayer position, thickness and texture (Cheng et al. 2013;Miangoleh et al. 2014). The heterogeneity of soils, combined with complex boundary conditions associated with natural terrains, pose considerable problems to researchers attempting to characterize infiltration at field research and laboratory research. The Hydrus-1D numerical software package based on the Richards equation provides a fast and accurate way for investigating infiltration, and is widely used for simulating various hydrologic processes in the vadose zone (Wang et al. 2014). For instance, Tan et al. (2014) used Hydrus-1D to simulate multi-layered soil water flow under ponded conditions in paddy fields, whereas Jiang et al. (2010) simulated water flow in nonuniform soils under natural climatic conditions where spray irrigation was used on a dairy farm. Similarly, Bourgeois et al. (2016) used the model to estimate the hydraulic properties of a shallow soil, including its water content, which was generally close to the measured value.  simulated water movement in layered water-repellent soils and found that soil water repellency was more important than interlayer position in controlling water movement.
However, the Hydrus-1D model has not been used to simulating water movement within a multi-layered reconstructed soil comprised of Yellow River sediment.
The main objectives of this study were to: 1) assess the performance of the Hydrus-1D model in simulating infiltration through multi-layered reconstructed soil composed of Yellow River sediment under various interlayer configurations; 2) investigate the effect of interlayer thickness and position on water movement, and 3) determine the reconstructed soil profile for different thicknesses of native soil.

Infiltration Experiments
Experimental materials included native soils (topsoil, subsoil) and Yellow River sediment. The native soil was obtained from an area of subsided land in Qihe County, Dezhou City, Shandong Province (36°28´52.69´´N, 116°28´3.18´´E). Yellow River sediment was collected from the Jinzhuang diversion canal (36°29´40.05´´N, 116°28´54.3´´E), which is approximately 1.89 km from the soil sampled subsidence area. The experimental materials were air dried and then sieved through a 2 mm mesh. The particle size distributions of the sampled soils were then determined using a pipette method.
Soil bulk density was measured by the gravimetric method (Table 1).
The infiltration experiments were carried out using polymethyl methacrylate cylinders (120 cm in length, 19 cm inner diameter). The top 5 cm of the cylinder was used for water application. The following 110 cm of the soil cylinder was packed with soil and compacted at 5-cm increments until reaching the bulk density specified in the experimental design (Wang et al., 2014) (Table 1). The surface of each soil layer was corrugated to increase its roughness before adding and compacting the next layer (Ma et al.,2010). During the filling and compacting of the cylinder, ten soil moisture sensors (ECH2O EC-5) were installed at 10, 30,45,55,65,75,85,95, and 105 cm below the soil surface. These sensors were used to measure the soil water content every 5 min. A Mariotte bottle was used to maintain 3cm constant water head ( Fig. 1(a)).
During the experiment, the time required for the wetting front to advance about every 5 cm was recorded, along with changes in the water volume within the Mariotte bottle. The infiltration experiment ended when the wetting front reached The experimental soil profiles were all 110 cm long, and consisted of two control check (CK1, CK2) and three experimental treatments (T1-T3). CK1 consisted of 20 cm of topsoil overlying a continuous layer of subsoil respected native farmland. CK2 represented the conventional method for reconstructing a soil profile in the area; it consisted of 20 cm of topsoil, followed by 20 cm of subsoil that overlaid 70 cm of sediment. T1-T3 consisted of various combinations of subsoil interlayer thicknesses and positions within the sediment, all of which were overlain by 20 cm of topsoil and 20 cm of subsoil (1(b)). Interlayer thickness of subsoil included 10 cm (T1) and 20 cm (T2), whereas interlayer position varied from 20 cm (T1, T2) to 30 cm (T3) ( Fig. 1(b)).

Numerical modeling
The Hydrus-1D model, applied in this study, is based on the Richards equation, a commonly used equation for simulating the movement of water in unsaturated soils. The Richards equation is arguably one of the most reliable and accurate equations to solve questions in the hydroscience (Farthing et al. 2017). It is expressed as: where h is the soil water pressure head (cm); z is the vertical elevation above a given datum (coordinate), which is positive upward (cm); t is infiltration time (min); is the soil water content (cm 3 cm -3 ); and K(h) is the hydraulic conductivity (cm min -1 ).
The Van Genuchten (1980) model is commonly used to describe the relationships between K(h) and (ℎ), as follows: where is the residual water content (cm 3 cm -3 ); is the saturated water content (cm 3 cm -3 ); is saturated hydraulic conductivity (cm min -1 ); and , m, n and l (set at 0.5 for all the treatments) are empirical parameters in the model.
The initial conditions assume that the initial soil water content of each layer is uniformly distributed.
where L is the length of the simulated soil column (cm); and h0 is the initial soil water pressure head (cm) that was calculated according to the initial soil water content of each layer.
The upper boundary condition was set at constant pressure head of 3 cm, and the lower boundary condition was free drainage. The boundary condition for water flow was expressed as: The Hydrus-1D was used to simulate the infiltration dynamics of reclaimed soils composed of Yellow River sediments.
The control checks (CK1, CK2) were used to calibrate the hydraulic parameters α, and n. The experimental treatments (T1-T3) were used to validate the modeled flow through the multi-layered soil profiles (Table 2). Simulated treatments consisted of 18 theoretical soil profiles that were divided into CK and two scenarios (50 cm and 60 cm thick native-soil layers) that possessed a similar total thickness of reclaimed soil (Table 3). 18 soil profiles were set where the thickness of overlying soil was less than 30 cm and the number of interlayers was one layer, according to greenhouse experiment (Hu et al. 2017) reclamation technology ) and the quantity of native soil. (Fig.2)In the case of total thickness of native soil was 50 cm, 7 treatments included the thickness of interlayers that were 10 cm and 20 cm, and the position of interlayers (the thickness of the L3) that were 0 cm, 10 cm, 20 cm and 30 cm, respectively. In the case of total thickness of native soil was 60 cm, 10 treatments included the thickness of interlayers that were 10 cm, 20 cm and 30 cm, and the position of interlayers (the thickness of the L3) that were 0 cm, 10 cm, 20 cm and 30 cm, respectively.

Statistical Evaluation of model performance
The performance of the model performance was evaluated using the Nash-Sutcliffe efficiency coefficient (NSE), Relative Root Mean-Squared Error (RRMSE), and coefficient of determination (R 2 ). Their mathematical expressions are shown in equations (7)-(9) respectively: where n is the number of samples, Si is the simulated value of infiltration, and Mi is the measured value of infiltration or evaporation.

Calibration of parameters
The parameters and KS were obtained using the ring sampler method and permeameter method, respectively. Other model parameters, , and ,were obtained by water retention curves (Fig.3). Soil water retention curves were measured using a pressure-membrane instrument. They were then calibrated until they model's result fit the observed values of cumulative infiltration (I), depth to the wetting front (Zf), and water content volume (θv) from the infiltration experiments that represented the control check (CK1 and CK2) ( Table 4).
Comparison between the measured and simulated I, Zf, and θv from control check infiltration case CK1 and CK2 (Fig.   4) the RRMSE values were smaller than 0.06, and the R 2 values were larger than 0.99 (Table 5). These statistical comparisons indicate that the calibrated hydraulic parameters and n used in the simulations were representative of those associated with infiltration into CK1 and CK2.

Validation of parameters
The calibrated hydraulic parameters were used to simulated I, Zf, and θv within the multilayered soil colmns (T1, T2 and T3), in which subsoil interlayers were emplaced into Yellow River sediment. The rate of the Zf migration decreased when the wetting front entered a subsoil interlayer. Comparison between measured and simulated values of I, Zf, and from T1, T2 and T3, demonstrated that they were similar (Fig. 5). In the case of T1, the NSE values for I, Zf and were larger than 0.96, the RRMSE values were smaller than 0.07, and the R 2 values were larger than 0.96. In the case of T2, the NSE values were larger than 0.97, the RRMSE values were smaller than 0.10 and the R 2 values were larger than 0.98. In the case of T3, the NSE values were larger than 0.97, the RRMSE values were smaller than 0.06 and the R 2 values were larger than 0.98 (Table 5). Therefore, the simulation was considered to be excellent (0.96<NSE<1, RRMSE<0.1 and 0.96<R 2 <1) (Jamieson et al., 1991).
The above statistical analyses indicate that the Hydrus-1D could well simulate infiltration characteristics of multi-layered reconstructed soil composed of Yellow River sediment when using calibrated parameters.

Simulation result from Hydrus-1D
In order to analyze the influence of interlayer thickness and position on the infiltration characteristics of the reconstructed soils. Hydrus-1D was used for simulating I, Zf and .
During the early stages of infiltration when water was moving through the same overlying topsoil, changes in I with time was consistent in all treatments and increased rapidly (Fig.6). Each of the experimental treatments exhibited differences in the timing of when I decreased (i.e., the inflection point of I curve, Fig. 5). For treatments in scenario 1, the time inflection of the I curve occurred at 720 min for T50-0, 480 min for T50-1 to T50-3, and 300 min for T50-4 to T50-6.
In the scenario 2 treatments, the I curve was 1100 min for T60-0, 720 min for T60-1 to T60-3, 480 min for T60-4 to T50-6, and 300 min for T60-7 to T60-9. The observed trends indicate that the timing of the changes in I (the inflection point) was determined by the thickness of the subsoil layer that covered the sediment, which is in accord with previous studies (Zhang et al., 2007;Fan et al., 2016.). It is also worth noting that a second inflection point in the I-curve occurred when a second subsoil interlayer was present (T50-1 to T50-6 and T60-1 to T60-9). The I value of the treatments also differed because of these turning points. In the case of scenario 1 treatments, when I was 34 cm, the infiltration time was 3470, 3400, 3300, 3260, 3160, 2870 and 2720 min for T50-0 to T50-6, respectively. These data showed that the infiltration time was shortened to varying degrees by putting a subsoil interlayer into the sediment; the required time gradually decreased with an increase in subsoil interlayer thickness when depth under the total thickness of subsoil is 30cm. The transit time of T50-6 was most closely to CK (2700 min). For treatments in scenario 2, when I was 34 cm, the infiltration time was 3420, 3470, 3440, 3470, 3320, 3197, 3220, 3000, 2750, and 2685 min for T60-0 to T60-9. The infiltration time gradually decreased with an increase in the thickness of the subsoil interlayers under the total thickness of subsoil is 40 cm. The effect of position of interlayer got stronger as the thickness of subsoil interlayers increased. Using a subsoil interlayer thickness of 10 cm did not improve the infiltration over that of conventional reconstructed soil profile (T60-0), in addition, the transit time of the water in T60-9 was most closely to CK. The above data suggested that removing part of the subsoil overlying the sediment and placing it between sediment layers could reduce the surface runoff of conventional reclaimed soil. In addition, thickness of the subsoil interlayer improved water infiltration rate more than position under the same total thickness of native soil was used for reclamation.
The relationships between Zf and time for the various treatments (Fig.7) showed similar trends at the initial Zf through the profile where all the materials were the same. However, the changes in Zf over time differed in reconstructed soil profiles containing Yellow River sediment (all treatments except CK) when the wetting front had to move across a boundary between a sediment and subsoil layer. The Zf moved downward faster when it was in sediment, but moved slower when the Zf entered a subsoil layer. Therefore, reconstructed soil profiles containing Yellow River sediment exhibited irregular Zf curves characterized by varying rates of downward movement. In the case of scenario 1 treatments, the time required for Zf to reach the bottom of the soil columns in treatments T50-0 to T50-6 was 3470, 3560, 3575, 3560, 3420, 3345 and 3330 min, respectively. The rate of downward movement in T50-2 was the slowest which was 86.85 % of the rate of CK.
Different thicknesses and positions of subsoil interlayer had different effects on Zf, but the differences were not significant (P > 0.05) when the total soil thickness was consistently 50 cm. In scenario 2, the time required for the Zf to reach the bottom of the columns in treatments T60-0 through T60-9 was 3608, 3830, 3864, 3823, 3785, 3824, 3760, 3500,3420 and 3390 min, respectively. The rate at which Zf increased beyond T60-0 was most significant when the interlayer thickness was 10 cm (T60-1 to T60-3). Moreover, the rate of downward movement in T60-2 was the slowest which was 93.87 % of the rate of CK. This indicates that the thickness of the interlayer affected the Zf migration speed more than position when the total of quantity (thickness) of soil interlayers was consistently 60 cm.
The distribution of in the soil profiles has an important influence on the absorption of water by plant roots. In order to describe the dynamic characteristics of the within the soil column (profile) that was determine when Zf equaled depths of 10,25,35,45,55,65,75,85,95,110 cm. The simulate values with depth in the soil columns for the various treatments are shown in Fig. 8. At the beginning of infiltration, the simulated results fluctuated at the interface between the subsoil and sediment. These fluctuations were due to the fact that the designed initial could not satisfy the hydrostatic condition. The was saturated in the native soil layer, while it was unsatured in the underlying sediment layers until the Zf passed through. This phenomenon has also been observed by other researchers (Hammecker et al., 2003;Rao et al., 2006;Wang et al., 2014). In fact, Zf within the sediment became unstable and broke into narrow wetting columns or "fingers" during some treatments (Hillel, et al., 1988;. At the end of infiltration, the within the soil profile of CK was uniformly distributed. However, was discontinuously distributed in reclaimed soil with sediment. Large values in the soil layer, and lower values occurred in the sediment layer. This confirmed the presence of an unwetted zone in the sediment layer. For example, the averaged of the soil and sediment layers for T50-0 were 48.52 and 24.02 cm 3 /cm 3 , and for T60-5 were 48.57 and 20.19 cm 3 /cm 3 , respectively. In contrast to the conventional reconstructed soil profiles (T50-0 and T60-0), the values of increased within the first sediment layer (below subsoil interlayer) of the multi-layered reconstructed soil profiles. For instance, the for the first sediment layer in treatments T50-1 to T50-6 were 39. 06, 38.16, 36.15, 41.92, 41.87and 41.71 cm 3 /cm 3 , respectively. The values for the first sediment layer in treatments T60-1 to T60-9 were 36. 63, 34.45, 31.24, 41.12, 40.25, 36.08, 42.00, 41.98, 41.80 and 41.71 cm 3 /cm 3 , respectively. The results showed that the values of within the first sediment layer gradually increased with an increase in the total subsoil interlayer thickness, and gradually decreased as the depth of the subsoil interlayers increased. Moreover, the influence of interlayer position on water content decreased as the interlayer thickness increased.
The thickness of the subsoil interlayer affected water infiltration more than position when the total thickness of the soil profile was the same. A comparison of the change in I and Zf with time also showed that as the infiltration rate increased, the rate of migration of Zf did not correspondingly increase. Furthermore, showed that both the topsoil and subsoil layers were saturated while the sediment layers were unsaturated as the wetting front, Zf , passed through.
Given the system's water balance, variations in water content ∆θ (cm 3 cm −3 ) and infiltration rate (cm min −1 ) determined the depth (cm) to which Zf will advance in a given time, t (min). Based on a piston type assumption (Jia and Tamai, 1997), the relationship can be expressed as: Therefore, the amount of water required per unit distance of Zf can represent the water-holding capacity of soil column during infiltration. In scenario 1, I for treatments T50-0 through T50-6 were 34. 03, 34.89, 35.54, 35.77, 35.36, 35.62 and 37.77 cm when Zf advancement to the soil columns respectively. Compared with T50-0, the water holding capacity of the reconstructed soils in scenario 1 increased (improved) by 2.53%, 4.44%, 5.11%, 3.91%, 4.67% and 10.99% for T50-1 through T50-6, respectively. The soil water-holding capacity reached a maximum in T50-6, which was 91.3% of CK, 41.35).

Conclusion
Multi-layered soil profiles were created by sequentially backfilling subsided areas with sediment delivered by hydraulic pump, followed by the addition of native soil layers over the sediment. The Hydrus-1D model was able to effectively The thickness and position of the interlayers was shown to affect infiltration characteristics of the reclaimed soils. In general, the thicker the native soil layer that was used, the better the soil water-holding capacity of the reconstructed soil profile. However, the quantity of native soil that is available in most subsidence areas is insufficient to use. The data collected herein demonstrated that removing part of the subsoil overlying the sediment and placing it between sediment layers will effectively improve the infiltration character and water-holding capacity of the conventional reconstructed soil profiles (T50-0 and T60-0) with the limit quantity of native soil. This Multi-layered soil profiles could reduce the surface runoff of conventional reclaimed soil. The thickness of the soil interlayer affected I, the migration speed of Zf, and the vertical distribution of θv more than the position of the interlayer when used total native soil was the same. Treatments in scenario 1 demonstrated that the optimal reconstructed soil profile (T50-6) had a soil water-holding capacity of 91.3% of the CK and the infiltration rate was most close to CK. In scenario 2, the optimal reconstructed soil profile (T60-9) possessed a water-holding capacity of this soil was 94.7% of the CK and the infiltration rate was most close to CK.
Overall, Numerical simulation provides the groundwork for designing field-scale experiments given the specific quantity of native soil, represents an important way to move forward soil reclamation. The study was greatly meaningful to improve the water character of the reclaiming damaged land by using river sediments in the world.