Validation and Modification of the Van Genuchten Model for Eroded Black Soil in Northeastern China

: The soil water characteristic curve is highly related to soil physical characteristics, which may be affected by the soil erosion degree. To explore the applicability of the soil moisture characteristic curve model in northeastern China, two erosion degrees, (1) lightly and (2) severely eroded black soil sampled from 15 sites, were collected. The soil water contents at eight soil water suctions were measured and the parameters of the Van Genuchten (VG) model were estimated. Then, two input sets—SSCBD (sand, silt and clay percentages and bulk density) and SSCBDθ 33 θ 1500 (SSCBD, and water contents at 33 and 1500 kPa suction) based on the Rosetta model were compared for the VG model prediction. The results showed that the parameters in the VG model had significant difference under the two eroded soils of the saturated water content (   ), but the opposite was true for the residual water content (   ), the scale parameter (  ) and the shape parameter (  ). In addition, the   and   had no significant differences but the opposite was true for the  and  under the two input sets. The simulated soil water content values of the VG model parameters derived from the Rosetta model underestimated the measured ones, except the water contents at 0 kPa. Therefore the relationships between m and n were modified for accuracy. The validation results showed that the VG model performed well when the sand content was less than 80% for the input set of SSCBD. Using the input set of SSCBDθ 33 θ 1500 can lead to higher simulation accuracy and wider applicability compared with SSCBD under black soil.


Introduction
The soil water characteristic curve describes the relationship between soil water potential and soil water content under unsaturated conditions, which is an important parameter for the study of soil water movement and retention [1]. The soil water characteristic curve can be measured directly using a pressure membrane gauge, centrifuge, or the tension meter method. However, these direct measurements are expensive, difficult, and time-consuming [2][3][4]. Some concise and indirect deducing models were successively proposed in the past several decades, such as those of Brooks-Corey [5], Gardner [6], Van Genuchten [7], Kosugi [8], and Seki [9], et al. These models serve as viable inexpensive and rapid ways to estimate soil hydraulic properties [10][11][12][13][14][15]. Among these models, the BC model (Brooks-Corey model) and VG model are widely used in related studies since the range of applicability is wider and accuracy is higher compared with other models. Some studies have shown that the VG model is more suitable to estimate the soil water content profile of soils with lower sand and higher clay content compared to the BC model [16,17]. Moreover, the VG model can describe more accurate soil water characteristic curves for farmland [18].
To determine the parameters of the VG model, undisturbed soil should be sampled, then the soil water contents under different soil suction conditions should be measured [9]. Since soil sampling might be affected by natural conditions and is more suitable for a small-scale study, the application of the VG model is limited. The soil water characteristic curve is closely related to soil physical and chemical characteristics [18][19][20][21][22][23][24]. Comparisons of soil water characteristic curves among different soil textures, landforms, and land uses have been conducted in several studies [1,[25][26][27]. Therefore, many researchers have proposed methods to estimate VG parameters with soil physicochemical properties [11,[28][29][30][31][32]. Among these studies, the Rosetta model is the most widely used currently [11,31]. The Rosetta model was developed based on artificial neural network analysis, which can generate parameters of the VG model through iterative correction between different input variables and output variables [4]. This method reduces the measurement workload and improves the efficiency of obtaining model parameters, especially for a large-scale study [33]. Recent publications related to the SSCBDθ33θ1500 to the Rosetta model to estimate the VG model parameters in the northeastern black soil. This study substantiates the application of the Rosetta model to estimate VG model parameters at a large scale for northeastern black soil.

Study Area
The study was conducted at the Jiusan Soil and Water Conservation Experimental Station of Beijing Normal University, which is located in the Heshan Farm (northwestern part of Heilongjiang Province, 125°16′ E-125°21′ E, 48°59′ N-49°03′ N). The study area is in the transitional zone from the southern foothills of the Xiaoxing'anling Mountains to the Songnen Plain. Elevation ranges from 310 to 390 m, with long and gentle slopes. The slope length is 800-1500 m and the slope angle is generally between 1° and 4° [49][50][51][52]. The climate is the mid-temperate continental monsoon climate, with an annual average air temperature of 0.4 °C and large temperature differences between winter and summer. Annual precipitation is around 500 mm, with 67% precipitation between June and August [50,52].
The primary soil type in this area is black soil, accounting for 64.2% of the total area, which is classified as Udic Isohumisols in the Chinese Soil Taxonomy, and Udic Argiboroll in the USDA Soil Taxonomy, or a Luvic Phaeozem in the FAO (Food and Agriculture Organization of the United Nations)/UNESCO (United Nations Educational, Scientific, and Cultural Organization) system [53]. The typical soil profiles in the study area can be divided into three layers, the topsoil which is porous and fertile (A), the illuvial horizon (B), and the parent material layer (C). The soil color of each layer is different, from layer (C) to layer (A) the soil color gets darker, and the clay content decreases as the gravel content increases with depth [54]. The average soil bulk density is 1.27 g cm −3 , with soil organic matter content ranging from 3% to 5% [53,55]. The area is divided by artificially planted windbreak forests, and ridges along the forests lead to three types of ridge direction, which are parallel, oblique or vertical to the contours. The corresponding planting patterns were planted on contour, cross slope or up-down slope.

Data Collection
Data were collected from 15 field sites in a watershed with an area of 28 km 2 . Soil samples were taken from all of the field sites at depths of 0-20, 20-40 cm, and 40-60 cm. Soil erosion occurring in the study area results in a thinner black soil layer, even disappearing in some locations. The field sites were classified by their erosion phases. Erosion phases were classified based on the structure of soil profiles (Table 1, Figure 1). The distribution of sampling sites appears in Figure 2. Table 1. Classification of soil erosion.

Eroded Soil Degree Indication
Lightly A significant portion of horizon A remaining Severely Exposure of horizon C or evidence of incorporation of horizon C material into exposed horizon B through cultivation A, mineral horizon; B, illuvial horizon; C, parent material.  The soil water characteristic curve, mechanical composition, and bulk density of the soil samples collected from all three layers of each field site were measured. (1) The soil water characteristic curve for a soil is defined as the relationship between water content and suction for the soil. Soil water contents under 8 suctions (0, 33, 50, 100, 300, 500, 1000, and 1500 kPa) were measured and produced a data set of 360 (15 soil sample sites × 3 layers × 8 soil water suctions) pairs of measured suction (h) versus water content. Soil water contents for unsaturated soil were measured using the pressure membrane apparatus; soil saturated water contents (0 kPa) were measured with several 100 cm 3cutting rings.

Model Calibration and Validation
The Van Genuchten (VG) model was applied in this study to describe the soil water characteristic curve. The expression is as follows, where is the effective water content; is the volumetric soil water content (cm 3 cm −3 ) when the matric potential is ℎ (kPa); is the saturated water content (cm 3 cm −3 ); is the residual water content (cm 3 cm −3 ), which is the soil water content in a very dry condition; is a scale parameter that is related to the inverse of the air entry suction (kPa −1 ), and , are shape parameters, and the default relationship between the two parameters is, According to classification of soil erosion in Table 1, 7 field sites (4 lightly eroded and 3 severely eroded soil field sites) of soil water contents were selected for fit, and the remaining 8 field sites of soil water contents (5 lightly eroded and 3 severely eroded soil field sites) were used for verification (Table 2). Firstly, 168 pairs of measured suction (h) versus water content from the 7 fitted field sites were used to fit a set of parameters in the VG model using the program under different eroded soil. Then, the derived set was used as input in the VG model equation and received another data set of soil water content. Finally, 192 pairs of measured suction (h) versus water content from the 8 validated field sites were used to compare with the results derived from the model. To determine the parameters of the VG model, undisturbed soil should be sampled, then the soil water contents under different soil suction conditions should be measured. The Rosetta model can predict the parameters , , and for the VG model, and the inputs of the Rosetta model are easier to collect. The Rosetta model has five options for input data: (a) input soil texture class; (b) input sand, silt and clay percentages; (c) input sand, silt and clay percentages and bulk density; (d) input sand, silt, and clay percentages, bulk density and water content at 33 kPa suction; (e) input sand, silt, and clay percentages, bulk density and water contents at 33 and 1500 kPa suction. This study used two options, sequence (c) with four inputs and sequence (e) with six inputs, to compare the fitting accuracy for the VG model (Table 3). Table 3. Five hierarchical models of the Rosetta model.

Model
Hierarchical Sequence Data Soil Texture Class Sand, Silt and Clay Percentages Bulk Density θ33 θ1500 θ33 and θ1500 represent water contents at 33 kPa and 1500 kPa suction, respectively.
The measured data set from 7 field sites (4 lightly eroded and 3 severely eroded soil field sites) were classified according to the Table 2. Hierarchical sequence data were obtained for 21 groups as the (c) and (e) model forms of the Rosetta model according to the Table 3. The 21 groups' fitted parameters as input of the VG model, including , , n and α, can obtain the 168 pairs of fitted suction (h) versus water contents by calculation. Model accuracy was verified by comparing with the measured values. If the model does not perform well according to the verification results, the relation between m and n in the VG model needs modification. Finally, 192 pairs of measured water contents of the 8 field sites (5 lightly eroded and 3 severely eroded soil field sites) in Table 2 were used to validate the performance of the revised VG model.

Accuracy Assessment Methods
The measured soil water contents were compared with the predicted soil water contents to assess the model performance by the mean of residuals (MR), Root Mean Square Error (RMSE), Akaike's information criterion (AIC) [56], and the 1:1 line regression method (one of the T test methods, by estimating whether the confidence interval of the slope and intercept of the regressed equation include the number of 1 and 0, respectively. If included, it indicates that there is no difference between the regressed curve and the 1:1 line. It further indicates that there is no difference between the simulated values and observed values) [39,56,57], which are calculated as follows, where N is the total number of events, is the measured value of soil water content, is the estimated value of soil water content, and q is the number of model parameters.

Validation of the VG Model
A total of 168 pairs of measured suction (h) versus water contents from the seven fitted field sites in Table 2 were used to estimate a set of parameters in the Van Genuchten (VG) model using the program under different eroded soil. T-tests showed that there were no significant differences for parameters , and while there was significant difference for parameter ( Table 4). The mean and were 0.00 cm 3 cm −3 and 1.11, respectively. The and of the lightly eroded soil were 0.41 cm 3 cm −3 and 0.27 kPa −1 , which were decreased by 26.83% and 81.48% for the severely eroded soil, respectively. Different letters indicate significant differences at p = 0.05.
The calculation soil water contents match well with measured values, indicating the VG model is suitable for the northeastern black soil (Figure 3). Compared with the severely eroded soil, the lightly eroded soil leads to better simulation quality using the VG model, except the saturated water content at 0 kPa ( Figure 3b, Table 5). The mean of residuals (MR) and Root Mean Square Error (RMSE), of the lightly eroded soils were lower by 147.73% and 30.9% than those of severely eroded soils, respectively (Table 5). In addition, the R 2 of lightly eroded soils was higher than that of severely eroded soils (Figure 3b).  For the lightly eroded soil, the measured soil water contents were higher than the simulated except the water content at 1500 kPa, while the measured soil water contents were lower than the simulated except the water content at 0 kPa for severely eroded soil (Figure 3a). The effects of fitting results of the model in medium to high suctions were better for lightly eroded soil, while the model predicted saturated water content better for severely eroded soil (Figure 3a). Figure 4 illustrates that there were no significant differences between parameters and obtained using two input sets, SSCBD (sand, silt and clay percentages and bulk density) and SSCBDθ33θ1500 (SSCBD, and water contents at 33 and 1500 kPa suction), but the differences between the scale parameters α and the shape parameters obtained using two input sets were significant at the = 0.05 significance level for the same erosion degree. The mean and for the lightly eroded soil samples were 1.29 and 1.12 times those for the severely eroded soil samples. Parameter obtained by SSCBD was 63.2% and 42.5% lower than SSCBDθ33θ1500 for different eroded soils. Parameter was 1.46 and 1.38 for the lightly and the severely eroded soils, respectively, which was 13.4% and 1.4% higher than that of SSCBDθ33θ1500. Comparison of the VG model parameters obtained using the Rosetta model with inputs of SSCBD and SSCBDθ33θ1500, respectively, for the lightly (L) and severely (S) eroded soil degrees. (a), (b), (c), and (d) is the residual water content ( ), saturated water content ( ), scale parameter ( ), and shape parameter ( ), respectively. Different letters indicate significant differences at p = 0.05 between output parameters using two inputs (SSCBD and SSCBDθ33θ1500) at the = 0.05 significance level with the t-test. Figure 5 shows that the input of SSCBDθ33θ1500 predicted better results of the soil water contents compared with SSCBD. The measured soil water contents were higher than the simulated contents except the water contents at 0 kPa for different eroded soils. Those results were further confirmed by the relatively low MR, RMSE and AIC values obtained with input SSCBDθ33θ1500, as shown in Table  6. For lightly eroded soil, using SSCBDθ33θ1500 significantly improved the quality of predicted soil water contents, while there were no improvements for severely eroded soil. Figure 5. Comparison of the measured and simulated soil water contents before correction, using the input of SSCBD and SSCBDθ33θ1500, respectively, for the lightly and severely eroded soils.  The results of Figure 5 show that the simulated water contents of the VG model using the Rosetta model were lower than the measured water contents. This may relate to the relationship between and described in Equation (2). To improve the accuracy of predicted soil water contents via the Rosetta model, the relationships between and were adjusted in this study for using SSCBD and SSCBDθ33θ1500, respectively. Since the measured soil water content was known, and parameters , , and could be obtained by the Rosetta model, new values of could be calculated according to Equation (1)  where Equation (6) was the modified relationship for SSCBD, and Equation (7) was for SSCBDθ33θ1500. Statistical t-tests indicated that the difference between the new linear relations and Equation (2) were significant at the = 0.05 significance level. Figure 6 showed that the adjustment of the relationship between m and n in the VG model will improve the accuracy of soil water contents predicted using parameters generated by the Rosetta model, except saturated water contents. Adding indexes of water contents at 33 and 1500 kPa suction can effectively improve the simulated accuracy of soil water contents predicted using parameters generated by the Rosetta model. Those results were further confirmed by the relatively low MR, RMSE and AIC values obtained after correction, as shown in Table 7. Figure 6. Comparison of the measured and simulated soil water contents after correction, using the input of SSCBD and SSCBDθ33θ1500, respectively, for the lightly and severely eroded soils. For lightly eroded soil, the water contents of fitting results of the model in low suctions were higher than those measured after correction. Compared using SSCBDθ33θ1500, the effects of fitting results when using SSCBD in saturated water contents were better. For severely eroded soil, the effects of fitting results when using SSCBDθ33θ1500 were improved after correction, while there were no significant improvements when using input sets of SSCBD. The water contents of fitting results of the model in low suctions were higher than those measured after correction using SSCBDθ33θ1500, while the fitting results when using SSCBD were all higher than measured.

Model Validation
After modification, the model validation samples showed that adding the soil suction of 33 kPa and 1500 kPa corresponding to soil water contents could effectively reduce the MR, RMSE and AIC between the simulated water content values of the model and the measured values, except for severely eroded soil at 0 kPa ( Table 8). The effects of the fitting results of lightly eroded soil, in terms of soil water contents, were better than those of severely eroded soil. For lightly eroded soil, the water contents of the fitting results for the model in low to medium suctions were higher than the measured contents, and the model was not good at fitting the suction at 0 kPa. For severely eroded soil, the water contents of the fitting results were higher than measured, except at 0 kPa using SSCBD and from 50 kPa to 1000 kPa using SSCBDθ33θ1500. Using SSCBD, the model only fitted better at 0 kPa. Using SSCBDθ33θ1500, the model was not good at fitting the suction at 0 kPa. Figure 7a showed that the simulated water contents' accuracy, excluding the samples with sand content less than 80%, was higher than the simulated accuracy with all samples. In addition, the difference between the two lines and the 1:1 line were significant at the = 0.05 significance level. The difference between the line using SSCBDθ33θ1500 and the 1:1 line was not significant at the = 0.05 significance level (Figure 7b). When it comes to using input of SSCBD, the simulated values have relatively higher precision for soil sand contents of less than 80%. When it comes to using the input of SSCBDθ33θ1500, all sample points were applied. Figure 7. Comparison of the measured and simulated soil water contents using the input of SSCBD (a) and SSCBDθ33θ1500 (b), respectively. * means there were significant differences at p = 0.05 between the 1:1 line and the fitting line with the t-test.

The Difference in Model Parameters
The soil water characteristic curve was found to be highly related to soil physical characteristics [58,59], which were affected by soil erosion degrees [17]. Previous research has shown that there were differences between soil moisture characteristic curves obtained from farmland, sand dune and meadow soil [60]. Bai [61] found that there was a negative correlation between soil's available water and sand content, whereas there was a significantly positive correlation between soil's available water and soil bulk density, sand content, clay content and organic matter content. As a result, the ability of soil water holding decreases after soil has eroded [62,63]. Therefore, the value of of the lightly eroded soil is higher than the severely eroded soil. The parameter of α is related to the inverse of the air entry suction. The higher the value is, the better the aeration of soil is. However, severe erosion will destroy the soil structure [60]. Therefore, the of the lightly eroded soil is better than the severely eroded soil. In addition, the scale parameters obtained using two input sets were significant at the = 0.05 significance level for the same erosion degree. Results achieved here were similar to the previous research [11] in that the correlation for increased considerably when one or two retention points were added to the predictors. The simulated water contents of the VG model using the Rosetta model were lower than the measured contents. This may relate to the relationship between and described in Equation (2). The equation was established by Mualem model [7,64] and selected the value of = 0.5 as a constant. However, some researchers found that for a given soil, the VG model performed most accurately when was considered as an unknown parameter instead of a constant [65]. Kosugi [66] found that (0 < < 1) was also related to the width of the pore radius distribution. Results showed that the m values calculated according to Equation (2) were overestimated. An overestimation of m was noted by Schapp's and Ghanbarian's results in [11] and [56].

Performance of Model
The effect of fitting results of lightly eroded soil in soil water contents was better than in severely eroded soil. The results achieved here are consistent with previous research in that the VG model is more suitable for sticky soil rather than sandy soil [67]. The effect of fitting results of the model in medium to high suctions was better. The results achieved here were similar to previous research that the soil water characteristic curve is closely related to soil texture, and the distribution of macropores has an effect on water holding capability under lower suction for soil [17,58].
The input of SSCBDθ33θ1500 predicted better results of the soil water contents compared with SSCBD. The results achieved here were similar to previous research [11] in that the RMSE values decreased for water retention when more predictors were used, and the estimations by model e are much better than the estimations by model c in Table 3. The simulated water contents of the VG model using the Rosetta model were higher than the measured contents after correction. This means that the water content overestimates values by the Rosetta model, which were inconsistent with Schapp's [11] and Liu's [68] results. The reason may be that the relationship between and was modified in this study. Although the range of applications for model using SSCBD was studied in this article, more research is needed to test whether the content of sand, silt and clay is suitable for the VG model in northeastern China using SSCBDθ33θ1500.
We can use the result of the model for the SWCC in this study to predict the soil water in the cropland in the black soil region. Then, the results could be used to compare the soil conservation measures for increasing the efficiency of water use in the future.

Conclusions
The applicability of the VG model for different eroded soils was explored and compared. The following results were achieved, (1) Parameters and of the VG model were influenced by soil erosion degrees. The and had no significant differences but the opposite was true for the and under the two input sets at the = 0.05 significance level.
(2) Using model e in Table 3 could improve the model fitting accuracy. The lightly eroded soil leads to better simulation quality than severely eroded soil using the VG model, except for the saturated water content. After modifying the relationships between m and n, the accuracy of the VG model was improved; the model fitting was better at 0 kPa using SSCBD for severely eroded soil. (3) With input SSCBD, the VG model performed well when the sand content was smaller than 80%.