Evaluating the Water Holding Capacity of Multilayer Soil Proﬁles Using Hydrus-1D and Multi-Criteria Decision Analysis

: In semi-arid climate regions of China, vegetation restoration on open pit mining lands is limited by soil moisture. However, multi-layered soil proﬁles can impede water inﬁltration into deeper underground, leaving more water stored in the root zone. Here, three types of soils with contrasting texture, sandy loam (SL), sand (S), and silt loam (SiL), were used to construct four multilayer proﬁles: SL-SiL, SL-S, SL-S-SiL, and SL-SiL-S. Silt loam was taken from the humus layer, which is more conducive to plant growth than other layers, and it was allocated to the ﬁrst layer in the four proﬁles, while sand and silt loam underlay the silt loam layer. Column experiments and Hydrus-1D simulation of the vertical inﬁltration and drainage process were performed: (1) The simulated results showed that when the sand layer underlay the sandy loam layer (SL-SiL and SL-S-SiL), the sandy loam layer could hold more water than the silt loam layer underlaying the sandy loam layer (SL-SiL and SL-SiL-S). The water content of the sandy loam layer in SL-SiL (95 cm) and SL-S-SiL (95 cm) was 28.3% higher than SL-SiL (74 cm) and 10.5% higher than SL-SiL-S (86 cm). (2) Both the measured and simulated cumulative inﬁltration and wetting front penetration time were positively related to the thickness of the silt loam layer and negatively related to the thickness of the sand layer. (3) The simulated inﬁltration rate, accumulation inﬁltration, and wetting front of the ﬁrst layer were una ﬀ ected by the texture of the underlying layer. According to multi-criteria decision analysis, SL-S-SiL had the best water holding capacity and was suggested for land reclamation in the open pit mine in our research.


Introduction
Inner Mongolia, which is characterized by a semi-arid climate, is currently an important open pit mining area in China. Open pit mining causes 2-11 times the land damage as underground mining [1]. Consequently, large tracts of land have been destroyed that must be reclaimed every year. One crucial task of soil reclamation is to reconstruct a new profile that is conducive to plant growth. In mining operations, the excavation and transport activities degrade the soil structure and hasten nutrient loss, so the soils typically used for reclamation purposes are generally characterized by high bulk density, poor water holding capacity, and low available nutrients [2][3][4]. A further challenge is that water deficit is unfortunately a major factor limiting the survival of planted vegetation in arid and semi-arid climates. It follows that the main purpose of soil reconstruction should be to make the reconstructed profile better able to retain more water in the root layer. 2

of 21
To trap more water for plant growth, layered soil profiles have drawn much attention and investigation by many researchers [5]. Studies have shown that layered soil profiles can retain more water than do homogeneous soil profiles [6][7][8][9]. For example, Zettl et al. [10] compared the field capacities within a 1 m depth of reclaimed soil profiles and demonstrated that soils with textural layering had higher field capacities than those with homogenous layering. According to work by Huang and colleagues [11][12][13], layered soil can improve soil water storage capacity and reduce nutrient loss due to the strong textural contrast created. Soils with textural layering may hinder vertical water movement and allow the soil to store more water than those lacking any textural layering, during both the infiltration and drainage processes. Increased water storage in layered soils was attributed to two phenomena: the capillary barrier and the hydraulic barrier [14].
With advances in computing, more and more research on soil moisture has begun to apply numerical calculation to simulate soil moisture movement [15][16][17]. Hydrus-1D has been proven to be a promising mechanistic model for simulating water movement in variably saturated or unsaturated media [18,19]. Previous research highlighted that layered soils could hold more water and impede water into deeper layers. In a layered soil, the water dynamics are affected by the inner layer properties and the thickness of the layers, as well as their spatial configuration [20,21]; hence, the movement of layered soil is a complex process. Many indices can be used to characterize the water holding capacity of the soil profile, such as infiltration rate, accumulation rate, and field capacity, but evaluation indices are usually incommensurable (i.e., criteria with different units) and contradictory. For example, Wang found that the fine texture of silt clay loam or clay loam layers beneath the loam layer were more conducive to increase the amount of water infiltration, but a coarse texture of loamy sand below the loam layer allowed more water to be retained in the topsoil; however, its cumulative infiltration was lower than the former. The cumulative infiltration and water content of topsoil are generally considered as positive indicators for choosing a soil profile (i.e., the higher cumulative infiltration or water content represents the better water holding capacity of the soil profile). It is hard to evaluate which is better when these criteria are conflicting. In this case, multi-criteria decision analysis (MCDA) has been suggested to tackle this problem.
To our knowledge, few studies have focused on evaluating the water holding capacity of the soil profile through the MCDA method. It is important to design and select an appropriate soil profile suitable for reclamation before commencing land reclamation. In the study, the laboratory experiments and Hydrus-1D simulation of the vertical infiltration and drainage process were performed on four soil profiles characterized by different soil properties, layer thicknesses, and spatial configurations. The main objectives were two-fold: (1) to evaluate the effect of soil texture, thickness, and layer ordering on water holding capacity by column experiments and numerical modelling; (2) to use MCDA to evaluate the soil profiles and provide decision making concerning profile optimization.

Research Area
The study area (Figure 1) was in ShengLi Open pit Mine, located in Xilinhot, Inner Mongolia province of northeastern China. The area had a semi-arid continental climate, with mean annual precipitation of less 300 mm and mean annual evaporation of about 1746 mm. The mean annual temperature was 18.5 • C with a short plant growing season from June to October. A typical undisturbed profile above the bedrock is shown in Figure 1, for which the corresponding soil texture can be found in Figure 2, and three samples were taken from each layer. The soil profile (5.7 m) was divided into five layers (Layers A-E) according to the properties of the soil, and the ratio of each layer is also shown in Table 1. Three types of soil taken from Layer A, Layer D, and Layer E were used for filling columns to reconstruct the soil profiles. Layer B was not considered a single layer because it was adjacent to the humus layer (Layer A), which was thin (about 30 cm), the single minimum stripping thickness being more than 50 cm in ShengLi open pit mine, so their separation by the mining machine in the actual striping process was quite difficult. Layer C was too thin (<30 cm) to be worth selecting and using. It should also be pointed out that the soil profile had spatial heterogeneity, in that the thickness of each soil layer (Layers A-E) may vary across different source sites. Field investigation did show that Layers D and E were ripe for land reclamation, and since Layer A was a humus layer, which was more conducive to plant growth than other layers, it was allocated to the first layer in the subsequent profile design schemes.

Soil Column Experiments
The column experiments were carried out in the laboratory, and the experimental apparatus is displayed in Figure 3. The selected soils (Layer A, Layer D, and Layer E) were air dried to a mass constant, then sifted through a 2 mm sieve. The column experiments were conducted to simulate the infiltration and drainage process through different soil profiles. Considering the soil properties, the thickness of the layers, as well as their spatial configuration, four profiles suitable for the mining area in Inner Mongolia were designed ( Figure 4). Soil with a sandy loam texture was taken from Layer A, sand from Layer D, and silt loam from Layer E. These soils were compacted to maintain the original bulk density level. To perform this experiment, four columns with a diameter of 15 cm and a height of 120 cm were used. The soil was piled to a height of 100 cm, with an empty layer of 20 cm on top of the column to let water stand without allowing surface runoff. The EC-5 (METER Group, Inc. USA) was used for water content measurements, which was based on time-domain reflectometry (TDR) with a ±2% error in any porous medium. To enable the probes to measure the soil moisture content at each soil layer of the soil columns accurately, the probes were inserted in each layer at 10 cm intervals as shown in Figure 3. The data of the wetting front, infiltration rate, and accumulation infiltration were recorded manually at 5 min intervals. These data were collected before the wetting front reached the bottom of each column. After the infiltration ended, the constructed columns were soaked for 48 h to allow them to drain freely. Meanwhile, the soil was covered with a polyethylene sheet to prevent soil evaporation. During the drainage process, the change of water content of the whole soil column was measured by EC-5 at 1 hour (1h), 2 hours (2h), 4 hours (4h), 8 hours (8h), 12 hours (12h), 24 hours (1d), 48 hours (2d), 96 hours (4d), and 168 hours (7d). The above experiment for each profile was replicated three times. The average moisture data were used for parameter inversion and the following analysis.

Soil Column Experiments
The column experiments were carried out in the laboratory, and the experimental apparatus is displayed in Figure 3. The selected soils (Layer A, Layer D, and Layer E) were air dried to a mass constant, then sifted through a 2 mm sieve. The column experiments were conducted to simulate the infiltration and drainage process through different soil profiles. Considering the soil properties, the thickness of the layers, as well as their spatial configuration, four profiles suitable for the mining area in Inner Mongolia were designed ( Figure 4). Soil with a sandy loam texture was taken from Layer A, sand from Layer D, and silt loam from Layer E. These soils were compacted to maintain the original bulk density level. To perform this experiment, four columns with a diameter of 15 cm and a height of 120 cm were used. The soil was piled to a height of 100 cm, with an empty layer of 20cm on top of the column to let water stand without allowing surface runoff. The EC-5 (METER Group, Inc. USA) was used for water content measurements, which was based on timedomain reflectometry (TDR) with a ±2% error in any porous medium. To enable the probes to measure the soil moisture content at each soil layer of the soil columns accurately, the probes were inserted in each layer at 10 cm intervals as shown in Figure 3. The data of the wetting front, infiltration rate, and accumulation infiltration were recorded manually at 5 min intervals. These data were collected before the wetting front reached the bottom of each column. After the infiltration ended, the constructed columns were soaked for 48 h to allow them to drain freely. Meanwhile, the soil was covered with a polyethylene sheet to prevent soil evaporation. During the drainage process, the change of water content of the whole soil column was measured by EC-5 at 1 hour (1h), 2 hours (2h), 4 hours (4h), 8 hours (8h), 12 hours (12h), 24 hours (1d), 48 hours (2d), 96 hours (4d), and 168 hours (7d). The above experiment for each profile was replicated three times. The average moisture data were used for parameter inversion and the following analysis.

Hydrus-1D Simulation
Hydrus-1D was developed to solve the Richards equation based on finite element methods, which is widely used to simulate water and solute transport in soil [22].

Hydrus-1D Simulation
Hydrus-1D was developed to solve the Richards equation based on finite element methods, which is widely used to simulate water and solute transport in soil [22].

Governing Equation
The vertical soil column infiltration and drainage experiment could be simulated by Richards' equation describing one-dimensional saturated and unsaturated soil water movement, which is expressed as follows: where t is time (T); θ is volumetric water content (L 3 L −3 ); h is the soil water pressure head (L); z is the vertical spatial coordinate (L) and taken positive upward; K is the unsaturated hydraulic conductivity (L T −1 ). Richards' equation was solved numerically by Hydrus-1D [22].

Soil Hydraulic Parameters
The relationship between soil matrix potential, hydraulic conductivity, and moisture content was fitted by the following equations: where θ is volumetric water content (L 3 L −3 ); θ r and θ s respectively represent the residual moisture content and saturated water content (L 3 L −3 ); α (L −1 ), n, and m are van Genuchten's equation parameters, and m = 1 − 1/n; K s is the saturated hydraulic conductivity (L T −1 ); S e is the effective saturation.

Particle Swarm Optimization
The inversion method is usually used to obtain the parameters, but Hydrus-1D is incapable of solving multiple parameters using its built-in least squares algorithm; therefore, a multi-parameter optimization algorithm, particle swarm optimization (PSO), was used here [23]. PSO is a global search algorithm, and the average root mean squared error (RMSE) of the probes in each layer was taken to be an objective function to calculate the inverse parameter values. In this study, the inversed parameters derived from one column were validated at the corresponding layer with the same texture of the remaining soil columns, and those parameters were selected for subsequent simulations.

Initial and Boundary Conditions
Initial and boundary conditions for the infiltration process of the experiment consisted of the following: where h i is the initial pressure head (cm) in the soil profile; h 0 is the pressure head at the soil surface, equal to 3 cm in this study; L is the depth coordinate of the soil surface and is equal to 100 cm based on the column depth. Initial and boundary conditions for the drainage process from a saturated condition with zero flux across the surface were as follows:

Model Performance Evaluation
The performance of Hydrus-1D was evaluated using the (1) Nash-Sutcliffe efficiencies (NSE) index, (2) the root mean square of the prediction error (RMSE), and (3) the mean absolute error (MAE). The mathematical definition of these statistical indices is as follows: where E i is the i th estimated value, M i is the i th measured value, M is the mean of the observed values, and N is the number of observations.

Multi-Criteria Decision Analysis Methods
Multiple criteria decision analysis (MCDA) is a framework for evaluating limited alternatives (e.g., various soil profiles in our research) against multiple criteria. There is a number of approaches available for addressing an MCDA problem. The MCDA methods tested in our research included: the analytic hierarchy process (AHP); the technique for order preference by similarity to ideal solution (TOPSIS), and Grey relation analysis (GRA).

AHP
The different procedural steps in the method were as follows: (1) Build a multi-level hierarchical structure with a goal at the top level, criteria (and sub-criteria) at the intermediate levels, and alternatives at the lowest level.
(2) Construct a judgment matrix using pairwise comparisons for all alternatives. The relative importance of each indicator is demonstrated in Equation (14): Water 2020, 12, 773 where A is the judgment matrix, n is the size of the matrix, and a ij is the relative importance of indicator i to indicator j, which ranges from 1 to 9. a ij = 1/a ij and a ij = 1 when i = j. The scale of relative importance is shown in Table 2. Very strong importance 9 Extreme importance 2, 4, 6,8 Intermediate values between the two adjacent judgments (3) Normalize the judgment matrix by dividing a ij of the judgment matrix with the sum of its column; the sum of each column is 1. Then, the weight vector can be obtained by averaging across the rows. Since it is normalized, the sum of all elements in the weight vector is 1. The weight vector shows relative weights among the things that we compare.
(4) The consistency ratio (CR) is used to evaluate the consistency of the judgment matrix. For computing CR, the following Equation (18) is applied: where CI represents the consistency index computed according to Equation (6) and λ max is the largest eigenvalue of the judgment matrix, which can be calculated from Equation (20).
According to Vaidya [24], if the CR is >0.1, then the judgment matrix is unreasonable and must be reset.

TOPSIS
TOPSIS is a method based on how close a limited number of evaluation criteria is to the idealized target [25]. The procedures for applying the TOPSIS method are described as follows: (1) Create a performance matrix with m alternatives and n criteria. The performance matrix is represented as: Water 2020, 12, 773 9 of 21 (2) The values in the performance matrix are normalized into the range [0, 1]. The normalized values of each element in the performance matrix can be calculated as follows: Positive indicator (the larger the better), Negative indicator (the smaller the better), The normalized matrix is constructed as: The weighted value of the normalized indicator (r ij ) is calculated by the following equations: where w i is the weight of the indicators and X mn is the normalized value of indicators.
(3) The ideal point (A + ) is a composite of the best performance values, while the negative ideal point (A − ) is a composite of the worst performance values. They are determined by the following equations (separately): (4) The Euclidean distances from r ij to the ideal point (A + ) and negative ideal point (A − ) are calculated by the following equations (5) The value of the closeness coefficient (C j ) is calculated by Equation (32). A larger value of closeness indicates the better performance of a profile:

GRA
The GRA method was proposed by Deng [26]. The process consists of the following steps (1) The process of normalization is addressed in the above Section 2.3.1.
(2) Generate the reference sequence (R j ) from the normalized matrix by taking the largest normalized value of each criterion as: (3) Construct the difference matrix by calculating the difference between a normalized term and its reference value as: (4) The Grey correlation coefficient for each term is determined as: where χ generally takes the value of 0.5. (5) A Grey correlation degree is calculated as: where ϕ i is the Grey correlation degree that indicates the magnitude of correlation measured between the reference sequence and the i th data sequence.

Criteria Weights
The weight value of each MCDA method can affect the evaluation results. In addition to AHP, TOPSIS and GRA do not specify how the weights are calculated. As mentioned above, the calculation process of AHP includes the weight vector, which is also a widely used subjective weight calculation method [27,28], and one objective weight calculation method, the entropy method, was also adopted. Two kinds of weight vectors were applied to the MCDA above. In total then, five hybrid methods (AHP, AHP-TOPSIS, AHP-GRA, entropy-TOPSIS, entropy-GRA) were applied to evaluate the four soil profiles. The evaluation and inversion steps involved in this study are presented in Figure 5.
TOPSIS and GRA do not specify how the weights are calculated. As mentioned above, the calculation process of AHP includes the weight vector, which is also a widely used subjective weight calculation method [27,28], and one objective weight calculation method, the entropy method, was also adopted. Two kinds of weight vectors were applied to the MCDA above. In total then, five hybrid methods (AHP, AHP-TOPSIS, AHP-GRA, entropy-TOPSIS, entropy-GRA) were applied to evaluate the four soil profiles. The evaluation and inversion steps involved in this study are presented in Figure 5.

Calibration and Validation of Hydrus-1D
The Hydrus-1D model was calibrated using the measured water content of SL-S-SiL during the infiltration and drainage process and validated using the corresponding SL-S, SL-SiL, and SL-SiL-S column data. The results (Table 3) showed that Hydrus-1D successfully captured the moisture movement of the infiltration and drainage process in the four soil profiles, which was also reported in [11]. In order to reduce the computational cost, a set of hydraulic parameter values was roughly estimated by the trial-and-error method first [29][30][31], and these hydraulic parameters would be further optimized as the initial input values of the particle swarm optimization (PSO). The initial hydraulic parameters are given in Table 4, and the searching space was within a 20% fluctuation of the initial value. The searching space was the effective range of inversion parameters when executing the PSO. In the calibration process, the root mean squared error (RMSE), Nash-Sutcliffe efficiencies (NSE), and mean absolute error (MAE) values of simulated water content of SL-S-SiL were 0.0170, 0.9840, and 0.0118, respectively. The results showed that the validation results were slightly worse than the calibration results, but the results remained at an accurate level, which ranged from 0.0107 to 0.0265 (cm 3 cm −3 ) for RMSE, from 0.9662 to 0.9900 (dimensionless) for NSE, and from 0.0066 to 0.0192 (cm 3 cm −3 ) for MAE. The low MAE, RMSE, and high NSE values indicated that Hydrus-1D performed well. This was consistent with previous studies [17,19], which reported an excellent performance of Hydrus-1D in simulating the water movement, and similar accuracy was also obtained. Consequently, we concluded that the error of the simulation was low enough so that the performance of the model was deemed acceptable for further research. Further, the PSO algorithm accurately converged at the optimal value (Table 4) without falling into a local optimization space; PSO reportedly had the ability to invert multiple parameters at the same time [22,23].

Water Infiltration Rate and Cumulative Infiltration of Different Soil Profiles
The simulated maximum cumulative infiltration occurred in the SL-SiL profile, at 44.0 cm, followed by SL-S-SiL (39.2 cm), SL-SiL-S (36.6 cm), and SL-S (28.4 cm), which was close to the measured cumulative infiltration ( Table 5). The measured cumulative infiltration of SL-SiL, SL-S-SiL, SL-SiL-S, and SL-S was 38.5 cm, 34.7 cm, 32.9 cm, and 27.0 cm, respectively. The cumulative infiltration was positively related to the thickness of the silt loam layer (i.e., second layer of SL-SiL and SL-SiL-S, third layer of SL-S-SiL) and negatively related to the thickness of sand layer (i.e., second layer of SL-S and SL-S-SiL, third layer of SL-SiL-S). The infiltration rate was high and unstable in the initial stage of the infiltration process (0-10 min), but it soon leveled off. In the initial stage, compared with the estimated infiltration rate with the same values, the measured infiltration rate ranged from 0.5 and 0.9 cm min −1 .
The measured infiltration rate at the initial stage was significantly higher than its simulated value ( Figure 6), which was attributed to the low bulk density layer at the surface and the instability of the ponded water depth at the beginning of the infiltration process [10]. Before the water flow entered the second layer, the infiltration rate trend appeared stable. Although it had not reached a completely stable state, the subsequent changes were difficult to capture in the experimental data records because the differences were simply too small. For example, before the water flow entered the second layer (t = 96 min), the infiltration rate of SL-SiL was 0.07 cm min −1 , and after 300 min, the infiltration rate was 0.05 cm min −1 , which only decreased by 0.02 cm min −1 .   The measured steady infiltration rate of the different profiles ranged from 0.04 to 0.08cm min −1 .
The simulation results showed that SL-S reached a stable infiltration rate (0.07 cm/min) at 110.4 min, being the fastest profile to do so, whereas it took at least three times longer (496.8 min) for SL-SiL to reach a stable infiltration rate (0.05 cm/min). The infiltration rate of SL-S-SiL was also maintained at 0.07 cm/min after 110.4 min, but after the wetting front moved across the whole soil column, subsequent simulation results showed its infiltration rate reduced at 501.6 min, where it reached a steady infiltration rate (0.05 min/cm). The reason was that the low Ks of the silt loam layer hindered the infiltration of sandy and sandy loam layers, resulting in a decreased infiltration rate. This phenomenon is called the hydraulic barrier when a coarse-textured soil The measured steady infiltration rate of the different profiles ranged from 0.04 to 0.08 cm min −1 . The simulation results showed that SL-S reached a stable infiltration rate (0.07 cm/min) at 110.4 min, being the fastest profile to do so, whereas it took at least three times longer (496.8 min) for SL-SiL to reach a stable infiltration rate (0.05 cm/min). The infiltration rate of SL-S-SiL was also maintained at 0.07 cm/min after 110.4 min, but after the wetting front moved across the whole soil column, subsequent simulation results showed its infiltration rate reduced at 501.6 min, where it reached a steady infiltration rate (0.05 min/cm). The reason was that the low K s of the silt loam layer hindered the infiltration of sandy and sandy loam layers, resulting in a decreased infiltration rate. This phenomenon is called the hydraulic barrier when a coarse-textured soil overlies a fine-textured one [14,20]. The infiltration rate of SL-SiL-S reached a steady state when t = 230.4 min, but the infiltration rate did not decrease further like SL-S-SiL, because the infiltration rates through the sandy loam and silt loam layer were both very slow, and the water content of both layers could reach a saturation state before water flow entered the third layer, thereby enabling SL-SiL-S to achieve a steady infiltration rate quickly. In addition, the first and second layer had similar hydraulic properties, resulting in a weak flow barrier (Figure 7). When the water flow of SL-SiL-S entered the third layer, the water hydraulic barrier did not occur due to the higher K s of the sandy layer than the silt loam layer, resulting in the water flux being allowed to enter the third layer less than its loss. The low water supply of the silt loam layer and the fast infiltration rate of the sandy layer prevented the latter from reaching saturation (Figure 8).
Water 2020, 12, x FOR PEER REVIEW 16 of 25 textured one [14,20]. The infiltration rate of SL-SiL-S reached a steady state when t = 230.4 min, but the infiltration rate did not decrease further like SL-S-SiL, because the infiltration rates through the sandy loam and silt loam layer were both very slow, and the water content of both layers could reach a saturation state before water flow entered the third layer, thereby enabling SL-SiL-S to achieve a steady infiltration rate quickly. In addition, the first and second layer had similar hydraulic properties, resulting in a weak flow barrier (Figure 7). When the water flow of SL-SiL-S entered the third layer, the water hydraulic barrier did not occur due to the higher Ks of the sandy layer than the silt loam layer, resulting in the water flux being allowed to enter the third layer less than its loss. The low water supply of the silt loam layer and the fast infiltration rate of the sandy layer prevented the latter from reaching saturation (Figure 8).

Wetting Front of Different Soil Profiles
There were no significant differences among profiles in the penetration time of their first layer ( Figure 9 and Table 6). It may be concluded that the same texture had the same infiltration characteristics [32]. These results showed that the texture of underlying layers did not affect the

Wetting Front of Different Soil Profiles
There were no significant differences among profiles in the penetration time of their first layer ( Figure 9 and Table 6). It may be concluded that the same texture had the same infiltration characteristics [32]. These results showed that the texture of underlying layers did not affect the advancing speed of the wetting front at the first layer. The same patterns of results were obtained in terms of cumulative infiltration and infiltration rate. The four profiles had the same initial infiltration rate due to the same soil properties of the first layer. However, the conclusion was difficult to draw from the measured results, because according to the measured data of Table 5, the initial infiltration rates of different profiles differed greatly (ranging from 0.04 to 0.08 cm min −1 ). The simulation results showed that after passing through the first layer, the advancing speed of the wetting front in each profile began to differ. Because of the sandy texture of the second layer in the SL-S, the advance speed of the wetting front in this profile remained relatively fast, and its penetration time was thus the shortest of all four profiles: it took just 303.4 min for the water to reach the bottom of SL-S (Table 6). SL-SiL had the longest penetration time, at 546.1 min. The penetration times of SL-S-SiL and SL-SiL-S were between those of SL-S and SL-SiL. Compared with SL-S-SiL (444.8 min), the wetting front penetration time of SL-SiL-S was slightly shorter (430.4 min). In SL-SiL, after passing through the first layer, the wetting front advance speed in the silt loam layer was always higher than that of SL-SiL-S in the silt loam layer, but the simulation results showed no such difference. The simulation performance of Hydrus-1D with respect to the wetting front was not as robust as that for cumulative infiltration and infiltration rate, which underestimated the wetting front of all soil profiles. This may have arisen from preferential flow. The numerical model presumed a uniformly advancing wetting front in a soil profile, but in the natural soil profile, it was non-uniform. This could therefore explain why the recorded (i.e., observed) wetting front data values were larger than the estimated (i.e., predicted) ones from Hydrus-1D.   Note: "-" means no data, SL-SiL and SL-S are a two-layer structure.

Water Distribution with Different Soil Profiles
The EC-5 probe was used to record the changes in volumetric water content at different positions of each profile, and this analysis was carried out to distinguish profiles achieving better water retention performance. The optimized parameters could successfully simulate soil water dynamics in the drainage process. Figure 8 shows the simulated volumetric water content profiles of each profile at 7d of the drainage process. The simulation results accurately captured the real-time water content of the multi-layered soil columns. For each profile, the simulated water profiles basically matched those obtained from the experimental measurements. The moisture content profile of soil columns showed obvious discontinuity at the interface of different soil layers, as the vertical variation in soil texture led to an abrupt change in the soil moisture content profile ( Figure 10).

Water Distribution with Different Soil Profiles
The EC-5 probe was used to record the changes in volumetric water content at different positions of each profile, and this analysis was carried out to distinguish profiles achieving better water retention performance. The optimized parameters could successfully simulate soil water dynamics in the drainage process. Figure 8 shows the simulated volumetric water content profiles of each profile at 7d of the drainage process. The simulation results accurately captured the realtime water content of the multi-layered soil columns. For each profile, the simulated water profiles basically matched those obtained from the experimental measurements. The moisture content profile of soil columns showed obvious discontinuity at the interface of different soil layers, as the vertical variation in soil texture led to an abrupt change in the soil moisture content profile ( Figure  10). At a drainage time of 7d, the total water holding capacity (Figure 11) was unlikely to further change, and the moisture content could be considered to have reached field capacity [10]. The main water loss came from the sandy layer, because the water content ratio of the sandy layer to the whole profile significantly declined. The sandy layer was located in the second layer of SL-S and SL-S-SiL and the third layer of SL-SiL-S. In SL-SiL, there was no significant change in the At a drainage time of 7d, the total water holding capacity ( Figure 11) was unlikely to further change, and the moisture content could be considered to have reached field capacity [10]. The main water loss came from the sandy layer, because the water content ratio of the sandy layer to the whole profile significantly declined. The sandy layer was located in the second layer of SL-S and SL-S-SiL and the third layer of SL-SiL-S. In SL-SiL, there was no significant change in the proportion of water content in the sandy loam layer (first layer) and silt loam layer (second layer), but the total water holding capacity was constantly diminished, showing that the sandy loam and silt loam layers had similar drainage rates, which can be seen from Figure 8. The total water holding capacity of SL-SiL was largest (285 mm), due to the high field capacity of silt loam. However, the water holding capacity of the first layer in SL-SiL (74 mm) was the poorest of the four soil profiles. The first layer of SL-S-SiL and SL-S (95 mm) had the highest water holding capacity. A plausible explanation was that in both SL-S-SiL and SL-S, the sandy loam layer (first layer) lied above the sandy layer (second layer), so upon reaching equilibrium (when the matrix potentials at the interface are equal in Figure 8), the water content of the sandy loam layer was now much higher than the underlying layers, thus retaining more water than its field capacity due to the capillary barrier [20]. Since the first layer of each profile was taken from the humus layer, it was important to vegetation restoration, being the main active layer of plant roots in the initial stage of land reclamation. Compared with the total water holding capacity of a profile, the water holding capacity of its first layer was more important for early stages of vegetation restoration. The water holding capacities of the first layer of SL-S-SiL and SL-S were same, being 28.4% higher than that of SL-SiL (74 mm) and 10.5% higher than that of SL-SiL-S (86 mm). The results indicated that the water holding capacity of SL-SiL was the worst.
Water 2020, 12, x FOR PEER REVIEW 21 of 25 but the total water holding capacity was constantly diminished, showing that the sandy loam and silt loam layers had similar drainage rates, which can be seen from Figure 8. The total water holding capacity of SL-SiL was largest (285 mm), due to the high field capacity of silt loam. However, the water holding capacity of the first layer in SL-SiL (74 mm) was the poorest of the four soil profiles. The first layer of SL-S-SiL and SL-S (95mm) had the highest water holding capacity. A plausible explanation was that in both SL-S-SiL and SL-S, the sandy loam layer (first layer) lied above the sandy layer (second layer), so upon reaching equilibrium (when the matrix potentials at the interface are equal in Figure 8), the water content of the sandy loam layer was now much higher than the underlying layers, thus retaining more water than its field capacity due to the capillary barrier [20]. Since the first layer of each profile was taken from the humus layer, it was important to vegetation restoration, being the main active layer of plant roots in the initial stage of land reclamation. Compared with the total water holding capacity of a profile, the water holding capacity of its first layer was more important for early stages of vegetation restoration. The water holding capacities of the first layer of SL-S-SiL and SL-S were same, being 28.4% higher than that of SL-SiL (74 mm) and 10.5% higher than that of SL-SiL-S (86 mm). The results indicated that the water holding capacity of SL-SiL was the worst.

Evaluation of Soil Profiles
The selected indices are shown in Tables 7 and 8, and their calculated weight values are shown in Figure 12. The weights of each index calculated by the AHP and entropy method were assigned different values. The entropy method assigned more weight to Index 1, Index 2, and Index 4 and

Evaluation of Soil Profiles
The selected indices are shown in Tables 7 and 8, and their calculated weight values are shown in Figure 12. The weights of each index calculated by the AHP and entropy method were assigned different values. The entropy method assigned more weight to Index 1, Index 2, and Index 4 and had similar values to AHP in Index 3 and Index 5, while the weight value of Index 7 calculated by the AHP was obviously higher than that calculated by the entropy method. Generally, the entropy weight method, as an objective weight calculation method, was more convincing because it excluded anthropogenic interference, but expert knowledge also had an important role in determining such weight values. In this study, since based on prior knowledge, Index 7 was judged to be more reliable than the others for characterizing the water holding capacity of soil profiles, more weight was assigned to it, and the reason is explained in Section 3.3. Table 7. Selected indices for soil profiles.

Indices Description Unit
Index 1 Steady infiltration rate cm min −1 Index 2 Accumulation infiltration amount cm Index 3 Breakthrough time of wetting front min Index 4 Total profile moisture mm Index 5 Total profile moisture at 7d mm Index 6 Water moisture of first layer at 7d mm had similar values to AHP in Index 3 and Index 5, while the weight value of Index 7 calculated by the AHP was obviously higher than that calculated by the entropy method. Generally, the entropy weight method, as an objective weight calculation method, was more convincing because it excluded anthropogenic interference, but expert knowledge also had an important role in determining such weight values. In this study, since based on prior knowledge, Index 7 was judged to be more reliable than the others for characterizing the water holding capacity of soil profiles, more weight was assigned to it, and the reason is explained in Section 3.3. Total profile moisture mm Index 5 Total profile moisture at 7d mm Index 6 Water moisture of first layer at 7d mm  Entropy-TOPSIS indicated that SL-SiL-S had the best performance, and AHP-TOPSIS indicated SL-S-SiL was the best, whereas L30S70 ranked the worst in Entropy-TOPSIS and AHP-TOPSIS. This result arose from different weight assignments, but the GRA methods (i.e., entropy-GRA and AHP-GRA) maintained the same result under the two weight calculation methods used. In this study, the results appeared less affected by the weight values. According to the five methods, SL-S-SiL was considered to have a better water holding effect, followed by SL-SiL-S, SL-SiL, and  Entropy-TOPSIS indicated that SL-SiL-S had the best performance, and AHP-TOPSIS indicated SL-S-SiL was the best, whereas L 30 S 70 ranked the worst in Entropy-TOPSIS and AHP-TOPSIS. This result arose from different weight assignments, but the GRA methods (i.e., entropy-GRA and AHP-GRA) maintained the same result under the two weight calculation methods used. In this study, the results appeared less affected by the weight values. According to the five methods, SL-S-SiL was considered to have a better water holding effect, followed by SL-SiL-S, SL-SiL, and lastly, SL-SiL ( Figure 13).
Except for entropy-TOPSIS, the remaining four MCDA methods had the same ordering, and the overall performance of SL-S-SiL was the best profile. Therefore, the profile of SL-S-Sil was suggested for land reclamation in ShengLi open pit mine.
Water 2020, 12, x FOR PEER REVIEW 23 of 25 lastly, SL-SiL ( Figure 13). Except for entropy-TOPSIS, the remaining four MCDA methods had the same ordering, and the overall performance of SL-S-SiL was the best profile. Therefore, the profile of SL-S-Sil was suggested for land reclamation in ShengLi open pit mine.

Conclusion
This study assessed the performance of four constructed soil profiles in terms of two key water processes: infiltration and drainage, and further elaborated the water holding mechanism of the tested multi-layered profiles. Based on those results and according to the five MCDA methods (AHP, AHP-TOPSIS, AHP-GRA, entropy-TOPSIS, entropy-GRA) used to evaluate the water holding effect of the four profiles, we were able to draw three conclusions: 1. The wetting front, cumulative infiltration, and infiltration rate of the first layer in each profile was unaffected by the texture of the underlying layer. Flow barriers caused by contrasting textured soil would hamper water from entering the underlying layer. When a fine textured layer overlay a coarse-textured layer, water stagnation ensued because of the capillary barrier, yet when a fine textured layer underlay a coarse-textured layer, it was instead caused by the hydraulic barrier.
2. In the drainage process, the results showed that the first layer had a water holding capacity ranking among profiles as follows: SL-S = SL-S-SiL > SL-SiL-S > SL-SiL. However, the water holding capacity of the whole soil profile took the order of SL-SiL > SL-S-SiL > SL-SiL-S > SL-S. The SL-SiL

Conclusions
This study assessed the performance of four constructed soil profiles in terms of two key water processes: infiltration and drainage, and further elaborated the water holding mechanism of the tested multi-layered profiles. Based on those results and according to the five MCDA methods (AHP, AHP-TOPSIS, AHP-GRA, entropy-TOPSIS, entropy-GRA) used to evaluate the water holding effect of the four profiles, we were able to draw three conclusions: 1. The wetting front, cumulative infiltration, and infiltration rate of the first layer in each profile was unaffected by the texture of the underlying layer. Flow barriers caused by contrasting textured soil would hamper water from entering the underlying layer. When a fine textured layer overlay a coarse-textured layer, water stagnation ensued because of the capillary barrier, yet when a fine textured layer underlay a coarse-textured layer, it was instead caused by the hydraulic barrier.
2. In the drainage process, the results showed that the first layer had a water holding capacity ranking among profiles as follows: SL-S = SL-S-SiL > SL-SiL-S > SL-SiL. However, the water holding capacity of the whole soil profile took the order of SL-SiL > SL-S-SiL > SL-SiL-S > SL-S. The SL-SiL profile had the largest total water holding capacity, but the poorest water holding capacity in the first layer.
3. The four profiles were comprehensively evaluated using five MCDA methods. In addition to AHP, the results showed a ranking among profiles of SL-S-SiL > SL-SiL-S > SL-SiL > SL-S. The profile of SL-S-SiL had the best soil water holding capacity and was suggested for land reclamation in ShengLi open pit mine.
Author Contributions: Formal analysis, X.W. and H.L.; investigation, X.W., S.C. and H.L.; methodology, X.W.; project administration, Y.Z. and W.X.; supervision, Y.Z. and W.X.; writing, original draft, X.W.; writing, review and editing, Y.Z. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare that they have no conflict of interest.