Next Article in Journal
Bed-Sediment Transport Conditions along the Sagavanirktok River in Northern Alaska, USA
Previous Article in Journal
The Color Formation Mechanism of the Blue Karst Lakes in Jiuzhaigou Nature Reserve, Sichuan, China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Institute of Land Reclamation and Ecological Restoration, China University of Mining & Technology (Beijing), Beijing 100083, China
2
China ENFI Engineering Co., Ltd., Beijing 100036, China
3
Department of Land Management, Zhejiang University, Hangzhou 310058, China
4
School of Mines, China University of Mining & Technology, Xuzhou 221116, China
*
Author to whom correspondence should be addressed.
Water 2020, 12(3), 773; https://doi.org/10.3390/w12030773
Submission received: 9 January 2020 / Revised: 14 February 2020 / Accepted: 7 March 2020 / Published: 11 March 2020
(This article belongs to the Section Hydraulics and Hydrodynamics)

Abstract

:
In semi-arid climate regions of China, vegetation restoration on open pit mining lands is limited by soil moisture. However, multi-layered soil profiles can impede water infiltration 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 profiles: 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 first layer in the four profiles, while sand and silt loam underlay the silt loam layer. Column experiments and Hydrus-1D simulation of the vertical infiltration 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 infiltration 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 infiltration rate, accumulation infiltration, and wetting front of the first layer were unaffected 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.

1. 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.
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.

2. Materials and Methods

2.1. 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.

2.2. 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.

2.3. 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].

2.3.1. 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:
θ t = t [ K ( h ) ( h z + 1 ) ]
where t is time (T); θ is volumetric water content (L3 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].

2.3.2. Soil Hydraulic Parameters

The relationship between soil matrix potential, hydraulic conductivity, and moisture content was fitted by the following equations:
θ ( h ) = { θ r + θ s θ r [ 1 + | α h | n ] m h < 0 θ s h 0
K ( h ) = { K s   S e 1 / 2 [ 1 ( 1 S e 1 / m ) m ] 2 h < 0 K s h 0
S e = θ θ r θ s θ r
where θ is volumetric water content (L3 L−3); θ r and θ s respectively represent the residual moisture content and saturated water content (L3 L−3); α (L−1), n , and m are van Genuchten’s equation parameters, and m = 1 1 / n ; Ks is the saturated hydraulic conductivity (L T−1); S e is the effective saturation.

2.3.3. 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.

2.3.4. Initial and Boundary Conditions

Initial and boundary conditions for the infiltration process of the experiment consisted of the following:
h ( z , t ) = h i ( z ) ,   t = 0
h ( z , t ) = h 0 , z = L
h ( z , t ) z = 0 ,   z = 0
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:
h ( z , t ) = 0 ,   t = 0
h ( z , t ) = 0 , z = 0
K ( h z + 1 ) = 0 ,   z = L

2.3.5. 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:
N S E = 1 i = 1 N ( M i E i ) 2 i = 1 N ( M i M ¯ ) 2
R M S E = 1 N i = 1 N ( E i M i ) 2
M E = 1 N i = 1 N ( E i M i )
where E i is the ith estimated value, M i is the ith measured value, M ¯ is the mean of the observed values, and N is the number of observations.

2.4. 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).

2.4.1. 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):
A = ( a i j ) n × n = [ 1 a 1 n a n 1 1 ]
where A is the judgment matrix, n is the size of the matrix, and a i j is the relative importance of indicator i to indicator j, which ranges from 1 to 9. a i j = 1 / a i j and a i j = 1 when i = j . The scale of relative importance is shown in Table 2.
(3) Normalize the judgment matrix by dividing a i j 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.
W ¯ i = A i j = 1 n A i
w i = 1 n i = 1 n W ¯ i
W = ( w 1 w 2 w i w n ) T
(4) The consistency ratio (CR) is used to evaluate the consistency of the judgment matrix. For computing CR, the following Equation (18) is applied:
CR = C I R I
C I = λ m a x n n 1
where CI represents the consistency index computed according to Equation (6) and λ m a x is the largest eigenvalue of the judgment matrix, which can be calculated from Equation (20).
λ m a x = 1 n i = 1 n ( A W ) i w i = 1 n i = 1 n ( j = 1 n a i j w i w i )
According to Vaidya [24], if the CR is >0.1, then the judgment matrix is unreasonable and must be reset.

2.4.2. 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:
X m n = [ X 11 X 1 n X m 1 X m n ]
(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),
X i j = X i j m i n { X i j } m i n { X i j } m i n { X i j }
Negative indicator (the smaller the better),
X i j = X i j m i n { X i j } m i n { X i j } m i n { X i j }
The normalized matrix is constructed as:
X m n = [ X 11 X 1 n X m 1 X m n ]
The weighted value of the normalized indicator ( r i j ) is calculated by the following equations:
r i j = w i X m n
where w i is the weight of the indicators and X m n 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):
A + = { r 1 + , r 2 + , , r n + }
A = { r 1 , r 2 , , r n }
r j + = { m a x j { r i j } i f   i   i s   a   p o s i t i v e   i n d i c a t o r m i n j { r i j } i f   i   i s   a   n e g a t i v e   i n d i c a t o r }
r j = { m a x j { r i j } i f   i   i s   a   n e g a t i v e   i n d i c a t o r m i n j { r i j } i f   i   i s   a   p o s i t i v e   i n d i c a t o r }
(4) The Euclidean distances from r i j to the ideal point ( A + ) and negative ideal point ( A ) are calculated by the following equations
D j + = j = 1 n ( r j + r i j ) 2 ,   j = 1 , 2 , , n
D j = j = 1 n ( r i j r j ) 2 ,   j = 1 , 2 , , n
(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:
C j = D j D j + + D j

2.4.3. 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:
R j = m a x i = 1 m { r i j }
(3) Construct the difference matrix by calculating the difference between a normalized term and its reference value as:
Δ i j = | R j X i j |
Δ = [ Δ Δ 1 n Δ m 1 Δ m n ]
(4) The Grey correlation coefficient for each term is determined as:
ψ i j = m i n i = 1 m m i n i = 1 n Δ i j + χ m a x i = 1 m m a x i = 1 n Δ i j Δ i j + χ m a x i = 1 M m a x i = 1 n Δ i j
where χ generally takes the value of 0.5.
(5) A Grey correlation degree is calculated as:
φ i = 1 n j = 1 m [ ψ i j ]
where φ i is the Grey correlation degree that indicates the magnitude of correlation measured between the reference sequence and the ith data sequence.

2.4.4. 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.

3. Results and Discussion

3.1. 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 (cm3 cm−3) for RMSE, from 0.9662 to 0.9900 (dimensionless) for NSE, and from 0.0066 to 0.0192 (cm3 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].

3.2. 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.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 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 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 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).

3.3. 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.

3.4. 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).
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.

3.5. Evaluation of Soil Profiles

The selected indices are shown in Table 7 and Table 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.
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 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.

4. 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.

Funding

This work was supported by the National Key Research and Development Program of China (2016YFC0501103).

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

  1. Li, M.S. Ecological restoration of mineland with particular reference to the metalliferous mine wasteland in China: A review of research and practice. Sci. Total Environ. 2006, 357, 38–53. [Google Scholar] [CrossRef] [PubMed]
  2. Xiao, W.; Lv, X.; Zhao, Y.; Sun, H.; Li, J. Ecological resilience assessment of an arid coal mining area using index of entropy and linear weighted analysis: A case study of Shendong Coalfield, China. Ecol. Indic. 2020, 109, 105843. [Google Scholar] [CrossRef]
  3. Lv, X.; Xiao, W.; Zhao, Y.; Zhang, W.; Li, S.; Sun, H. Drivers of spatio-temporal ecological vulnerability in an arid, coal mining region in Western China. Ecol. Indic. 2019, 106, 105475. [Google Scholar] [CrossRef]
  4. Shrestha, R.K.; Lal, R. Ecosystem carbon budgeting and soil carbon sequestration in reclaimed mine soil. Environ. Int. 2006, 32, 781–796. [Google Scholar] [CrossRef] [PubMed]
  5. Yan, W.M.; Chiu, C.F.; Yuen, K.V. Prediction and modeling of permeability function and its application to the evaluation of breakthrough suction of a two-layer capillary barrier. Can. Geotech. J. 2017, 54, 778–788. [Google Scholar] [CrossRef]
  6. Miller, D.E.; Gardner, W.H. Water infiltration into stratified soil 1. Soil Sci. Soc. Am. J. 1962, 26, 115–119. [Google Scholar] [CrossRef]
  7. Alfnes, E.; Kinzelbach, W.; Aagaard, P. Investigation of hydrogeologic processes in a dipping layer structure: 1. The flow barrier effect. J. Contam. Hydrol. 2004, 69, 157–172. [Google Scholar] [CrossRef]
  8. Hübner, R.; Günther, T.; Heller, K.; Noell, U.; Kleber, A. Impacts of a capillary barrier on infiltration and subsurface stormflow in layered slope deposits monitored with 3-D ERT and hydrometric measurements. Hydrol. Earth Syst. Sci. 2017, 21, 5181–5199. [Google Scholar] [CrossRef] [Green Version]
  9. Tan, S.H.; Wong, S.W.; Chin, D.J.; Lee, M.L.; Ong, Y.H.; Chong, S.Y.; Kassim, A. Soil column infiltration tests on biomediated capillary barrier systems for mitigating rainfall-induced landslides. Environ. Earth Sci. 2018, 77, 589. [Google Scholar] [CrossRef]
  10. Zettl, J.D.; Barbour, S.L.; Huang, M.; Si, B.C.; Leskiw, L.A. Influence of textural layering on field capacity of coarse soils. Can. J. Soil Sci. 2011, 91, 133–147. [Google Scholar] [CrossRef]
  11. Huang, M.; Barbour, L.S.; Elshorbagy, A.; Zettl, J.D.; Si, B.C. Infiltration and drainage processes in multi-layered coarse soils. Can. J. Soil Sci. 2011, 91, 169–183. [Google Scholar] [CrossRef]
  12. Huang, M.; Barbour, L.S.; Elshorbagy, A.; Zettl, J.D.; Si, B.C. Water availability and forest growth in coarse-textured soils. Can. J. Soil Sci. 2011, 91, 199–210. [Google Scholar] [CrossRef]
  13. Huang, M.; Spies, J.; Barbour, S.L.; Si, B.C.; Zettl, J. Impact of textural layering on water retention within drained sand profiles. Soil Sci. 2013, 178, 496–504. [Google Scholar] [CrossRef]
  14. Si, B.; Dyck, M.; Parkin, G.W. Flow and transport in layered soils. Can. J. Soil Sci. 2011, 91, 127–132. [Google Scholar] [CrossRef] [Green Version]
  15. Jha, R.K.; Sahoo, B.; Panda, R.K. Modeling the water and nitrogen transports in a soil–paddy–atmosphere system using HYDRUS-1D and lysimeter experiment. Paddy Water Environ. 2017, 15, 831–846. [Google Scholar] [CrossRef]
  16. Ma, Y.; Feng, S.; Song, X. Evaluation of optimal irrigation scheduling and groundwater recharge at representative sites in the North China Plain with SWAP model and field experiments. Comput. Electron. Agric. 2015, 116, 125–136. [Google Scholar] [CrossRef]
  17. Karandish, F.; Šimůnek, J. An application of the water footprint assessment to optimize production of crops irrigated with saline water: A scenario assessment with HYDRUS. Agric. Water Manag. 2018, 208, 67–82. [Google Scholar] [CrossRef] [Green Version]
  18. Chen, M.; Willgoose, G.R.; Saco, P.M. Spatial prediction of temporal soil moisture dynamics using HYDRUS-1D. Hydrol. Process. 2014, 28, 171–185. [Google Scholar] [CrossRef]
  19. Azad, N.; Behmanesh, J.; Rezaverdinejad, V.; Abbasi, F.; Navabian, M. Developing an optimization model in drip fertigation management to consider environmental issues and supply plant requirements. Agric. Water Manag. 2018, 208, 344–356. [Google Scholar] [CrossRef]
  20. Li, X.; Chang, S.X.; Salifu, K.F. Soil texture and layering effects on water and salt dynamics in the presence of a water table: A review. Environ. Rev. 2014, 22, 41–50. [Google Scholar] [CrossRef]
  21. Wu, Q.F.; Fan, J.; Yang, X.L.; Pan, Y.W.; Wang, Y.F.; Qiao, Y.Q.; Wang, S. Experiment and simulation of Infiltration from layered soils in open pit mine in Jin-Shaan-Meng adjacent region. Acta Pedol. Sin. 2015, 52, 86–96. [Google Scholar]
  22. Brunetti, G.; Šimůnek, J.; Piro, P. A comprehensive numerical analysis of the hydraulic behavior of a permeable pavement. J. Hydrol. 2016, 540, 1146–1161. [Google Scholar] [CrossRef] [Green Version]
  23. Zambrano-Bigiarini, M.; Rojas, R. A model-independent Particle Swarm Optimisation software for model calibration. Environ. Model. Softw. 2013, 43, 5–25. [Google Scholar] [CrossRef]
  24. Vaidya, O.S.; Kumar, S. Analytic hierarchy process: An overview of applications. Eur. J. Oper. Res. 2006, 169, 1–29. [Google Scholar] [CrossRef]
  25. Behzadian, M.; Khanmohammadi Otaghsara, S.; Yazdani, M.; Ignatius, J. A state-of the-art survey of TOPSIS applications. Expert Syst. Appl. 2012, 39, 13051–13069. [Google Scholar] [CrossRef]
  26. Tscheikner-Gratl, F.; Egger, P.; Rauch, W.; Kleidorfer, M. Comparison of multi-criteria decision support methods for integrated rehabilitation prioritization. Water 2017, 9, 68. [Google Scholar] [CrossRef] [Green Version]
  27. Ramya, S.; Devadas, V. Integration of GIS, AHP and TOPSIS in evaluating suitable locations for industrial development: A case of Tehri Garhwal district, Uttarakhand, India. J. Clean. Prod. 2019, 238, 117872. [Google Scholar] [CrossRef]
  28. Zhu, S.; Li, D.; Feng, H.; Gu, T.; Zhu, J. AHP-TOPSIS-based evaluation of the relative performance of multiple neighborhood renewal projects: A case study in Nanjing, China. Sustainability 2019, 11, 4545. [Google Scholar] [CrossRef] [Green Version]
  29. Wang, J.; Huang, Y.; Long, H.; Hou, S.; Xing, A.; Sun, Z. Simulations of water movement and solute transport through different soil texture configurations under negative-pressure irrigation. Hydrol. Process. 2017, 31, 2599–2612. [Google Scholar] [CrossRef]
  30. Bonfante, A.; Basile, A.; Acutis, M.; De Mascellis, R.; Manna, P.; Perego, A.; Terribile, F. SWAP, CropSyst and MACRO comparison in two contrasting soils cropped with maize in Northern Italy. Agric. Water Manag. 2010, 97, 1051–1062. [Google Scholar] [CrossRef]
  31. Jiang, J.; Feng, S.; Huo, Z.; Zhao, Z.; Jia, B. Application of the SWAP model to simulate water-salt transport under deficit irrigation with saline water. Math. Comput. Model. 2011, 54, 902–911. [Google Scholar] [CrossRef]
  32. Hillel, D.; Talpaz, H. Simulation of soil water dynamics in layered soils. Soil Sci. 1977, 123, 54–62. [Google Scholar] [CrossRef]
Figure 1. The location of the research area and a typical natural profile.
Figure 1. The location of the research area and a typical natural profile.
Water 12 00773 g001
Figure 2. Soil texture distribution (USDA) of the soil samples collected from the research area.
Figure 2. Soil texture distribution (USDA) of the soil samples collected from the research area.
Water 12 00773 g002
Figure 3. Experimental apparatus used in the infiltration and drainage process.
Figure 3. Experimental apparatus used in the infiltration and drainage process.
Water 12 00773 g003
Figure 4. Schematic diagram of the four designed soil profiles. SL, sandy loam; SiL, silt loam; S, sand.
Figure 4. Schematic diagram of the four designed soil profiles. SL, sandy loam; SiL, silt loam; S, sand.
Water 12 00773 g004
Figure 5. The flowchart of soil profile evaluation.
Figure 5. The flowchart of soil profile evaluation.
Water 12 00773 g005
Figure 6. Measured and simulated infiltration rate and accumulation.
Figure 6. Measured and simulated infiltration rate and accumulation.
Water 12 00773 g006
Figure 7. Optimized soil water retention curves and hydraulic conductivity for all soils.
Figure 7. Optimized soil water retention curves and hydraulic conductivity for all soils.
Water 12 00773 g007
Figure 8. Measured and simulated water content at different depths in the infiltration process, (a) SL-SiL; (b) SL-S; (c) SL-S-SiL and (d) SL-SiL-S.
Figure 8. Measured and simulated water content at different depths in the infiltration process, (a) SL-SiL; (b) SL-S; (c) SL-S-SiL and (d) SL-SiL-S.
Water 12 00773 g008
Figure 9. Measured and simulated wetting fronts in the soil columns.
Figure 9. Measured and simulated wetting fronts in the soil columns.
Water 12 00773 g009
Figure 10. Simulated water content in the drainage process at (a) SL-SiL; (b) SL-S; (c) SL-S-SiL and (d) SL-SiL-S.
Figure 10. Simulated water content in the drainage process at (a) SL-SiL; (b) SL-S; (c) SL-S-SiL and (d) SL-SiL-S.
Water 12 00773 g010
Figure 11. The actual water holding capacity of each layer at different times, (a) SL-SiL; (b) SL-S; (c) SL-S-SiL and (d) SL-SiL-S.
Figure 11. The actual water holding capacity of each layer at different times, (a) SL-SiL; (b) SL-S; (c) SL-S-SiL and (d) SL-SiL-S.
Water 12 00773 g011
Figure 12. Weight values calculated by AHP and the entropy method.
Figure 12. Weight values calculated by AHP and the entropy method.
Water 12 00773 g012
Figure 13. The ranking of soil profiles based on five MCDA methods. GRA, Grey relation analysis.
Figure 13. The ranking of soil profiles based on five MCDA methods. GRA, Grey relation analysis.
Water 12 00773 g013
Table 1. Physical properties of the experimental soils.
Table 1. Physical properties of the experimental soils.
SampleClay (%)Silt (%)Sand (%)TextureBulk Density (g cm−3)Porosity (%)The Proportion of Total Profile (%)
Layer A4.5537.7357.72Sandy Loam1.230.514.6%
Layer B4.4334.0461.54Sandy Loam1.590.3921.8%
Layer C3.3221.874.87Loamy Sand1.720.332.2%
Layer D0.152.6997.16Sand 1.540.3735.6%
Layer E7.4471.7220.83Silt Loam1.220.5135.9%
Table 2. The fundamental scale.
Table 2. The fundamental scale.
Intensity of ImportanceDefinition
1Equal importance
3Moderate importance
5Strong importance
7Very strong importance
9Extreme importance
2, 4, 6, 8Intermediate values between the two adjacent judgments
Table 3. The simulation performance using optimized parameters by Hydrus-1D.
Table 3. The simulation performance using optimized parameters by Hydrus-1D.
DatasetSoil ProfileRMSENSEMAE
Validation dataSL-S-SiL0.0170 ± 0.00360.9840 ± 0.00900.0118 ± 0.0041
Calibration dataSL-SiL0.0265 ± 0.00450.9662 ± 0.02610.0192 ± 0.0052
SL-S0.0107 ± 0.00310.9900 ± 0.00870.0066 ± 0.0019
SL-SiL-S0.0207 ± 0.00440.9679 ± 0.01110.0135 ± 0.0025
Table 4. Optimized parameters used in PSO and the search space. Residual moisture content ( θ r ), saturation moisture content ( θ s ), a and n the shape parameters of the soil water characteristic curves, and saturated hydraulic conductivity ( K s ).
Table 4. Optimized parameters used in PSO and the search space. Residual moisture content ( θ r ), saturation moisture content ( θ s ), a and n the shape parameters of the soil water characteristic curves, and saturated hydraulic conductivity ( K s ).
Soil TextureParametersInitial
Values
Lower
Boundary
Upper
Boundary
Optimized Values
Sandy loam θ r (cm3 cm−3)0.0000.000.0500.003
θ s (cm3 cm−3)0.4400.3740.5060.450
a (cm−1)0.0200.0170.0230.019
n 1.5301.3011.7611.514
K s (cm min−1)0.5000.0510.0690.06
Sand θ r (cm3 cm−3)0.0000.000.0500.007
θ s (cm3 cm−3)0.3100.2640.3570.312
a (cm−1)0.0550.0470.0630.059
n 2.5002.1252.8752.599
K s (cm min−1)0.5000.2130.2880.238
Silt loam θ r (cm3 cm−3)0.000.0000.0500.016
θ s (cm3 cm−3)0.4700.4000.5410.479
a (cm−1)0.0200.0170.0230.021
n 1.4001.1901.6101.462
K s (cm min−1)0.5000.0340.0460.045
Table 5. Statistics of the infiltration characteristics of various soil profiles.
Table 5. Statistics of the infiltration characteristics of various soil profiles.
ProfilesTypeInitial Infiltration
(cm min−1)
Steady Infiltration Rate (cm min−1)Accumulation Infiltration (cm)
SL-SiLMeasured0.80.0438.5
Simulated0.30.0544.0
SL-SMeasured0.50.0827.0
Simulated0.30.0728.4
SL-S-SiLMeasured0.90.0634.7
Simulated0.30.0539.2
SL-SiL-SMeasured0.90.0632.9
Simulated0.30.0636.6
Table 6. The penetration times of the wetting front through the soil profiles.
Table 6. The penetration times of the wetting front through the soil profiles.
ProfileData TypeInfiltration Time of Each Layer (min)
First LayerSecond LayerThird LayerTotal Time
SL-SiLMeasured77.1469.0546.1
Simulated96.0546.0642.0
SL-SMeasured79.2224.2303.4
Simulated96.0235.5331.5
SL-S-SiLMeasured79.8105.8259.2444.8
Simulated96.099.0288.0483.0
SL-SiL-SMeasured83.1193.1154.2430.4
Simulated96.0204.0165.0465.0
Note: “—“ means no data, SL-SiL and SL-S are a two-layer structure.
Table 7. Selected indices for soil profiles.
Table 7. Selected indices for soil profiles.
IndicesDescriptionUnit
Index 1Steady infiltration rate cm min−1
Index 2Accumulation infiltration amount cm
Index 3Breakthrough time of wetting front min
Index 4Total profile moisture mm
Index 5Total profile moisture at 7d mm
Index 6Water moisture of first layer at 7dmm
Table 8. Evaluation indices of different soil profiles.
Table 8. Evaluation indices of different soil profiles.
ProfileIndex 1Index 2Index 3Index 4Index 5Index 6
SL-S0.0728.4331.535313795.0
SL-SiL0.0544.0642.048128574.0
SL-S-SiL0.0539.2483.042622295.0
SL-SiL-S0.0636.6465.040822486.0

Share and Cite

MDPI and ACS Style

Wang, X.; Zhao, Y.; Liu, H.; Xiao, W.; Chen, S. Evaluating the Water Holding Capacity of Multilayer Soil Profiles Using Hydrus-1D and Multi-Criteria Decision Analysis. Water 2020, 12, 773. https://doi.org/10.3390/w12030773

AMA Style

Wang X, Zhao Y, Liu H, Xiao W, Chen S. Evaluating the Water Holding Capacity of Multilayer Soil Profiles Using Hydrus-1D and Multi-Criteria Decision Analysis. Water. 2020; 12(3):773. https://doi.org/10.3390/w12030773

Chicago/Turabian Style

Wang, Xin, Yanling Zhao, Huifang Liu, Wu Xiao, and Shuzhao Chen. 2020. "Evaluating the Water Holding Capacity of Multilayer Soil Profiles Using Hydrus-1D and Multi-Criteria Decision Analysis" Water 12, no. 3: 773. https://doi.org/10.3390/w12030773

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop