Effect of Shelterbelt Construction on Soil Water Characteristic Curves in an Extreme Arid Shifting Desert

Abstract

To explore the impact of artificial shelterbelt construction with saline irrigation on the soil water characteristic curve (SWCC) of shifting sandy soil in extreme arid desert areas, three treatments including under the shelterbelt (US), bare land in the shelterbelt (BL) and shifting sandy land (CK) in the hinterland of the Taklimakan Desert were selected. The age of the shelterbelt is 16, and the vegetation cover is mainly Calligonum mongolicum. The soils from different depths of 0–30 cm were taken keeping in view the objective of the study. The SWCCs were determined by the centrifugal method and fitting was performed using various models such as the Gardner (G) model, Brooks–Corey (BC) model and Van Genuchten (VG) model. Then, the most suitable SWCC model was selected. The results showed that electrical conductivity (EC) and organic matter content of BL and US decreased with the increasing soil depth, while the EC and organic matter content of CK increased with the soil depth. The changes in soil bulk density, EC and organic matter of 0–5 cm soil were mostly significant (p < 0.05) for different treatments, and the differences in SWCCs were also significant among different treatments. Moreover, the construction of an artificial shelterbelt improved soil water-holding capacity and had the most significant impacts on the surface soil. The increase in soil water-holding capacity decreased with increasing soil depth, and the available soil water existed in the form of readily available water. The BC model and VG model were found to be better than the G model in fitting results, and the BC model had the best fitting result on CK, while the VG Model had the best fitting result on BL with higher organic matter and salt contents. Comparing the fitting results of the three models, we concluded that although the fitting accuracy of the VG model tended to decrease with increasing organic matter and salinity, the VG model had the highest fitting accuracy when comparing with BC and G models for the BL treatment with high organic matter and salinity. Therefore, the influence of organic matter and salinity should be considered when establishing soil water transfer function.

1. Introduction

Soil moisture is an important factor affecting plant growth and is a major driving force for the sustainability of many terrestrial ecosystems. Moisture changes have significant impacts on vegetation and soil properties [1,2,3]. Especially in arid and semi-arid regions, soil moisture availability is one of main factors limiting the type and quantity of vegetation, and the water deficiency can lead to severe degradation of vegetation and reduce vegetation cover [4,5]. Therefore, it is important to understand the soil moisture change pattern for the maintenance of vegetation in arid desert areas [6].

The Taklamakan Desert is located in the hinterland of the Eurasian continent in Chinese southern Xinjiang, and is the largest desert in China and the second largest shifting desert in the world. Southern Xinjiang is one of the poorest regions with an extreme drought climate in China. To accelerate the development of regional social economy, the Taklimakan Desert Highway (TDH) was completed in 1995. TDH is 522 km long running across the Taklimakan Desert from north to south and is the longest highway crossing a moving desert in the world. However, serious sand disasters are a great threat to TDH. Thus, in 2003, a shelter–forest belt (the Taklimakan Desert Highway Shelterbelt (TDHS)) having a length of 436 km was constructed on both sides along the highway, which was dominated by shrubs to protect the highway from shifting sand [7]. Low-quality saline groundwater is used for irrigation in order to ensure the survival of shelterbelt plants [8,9,10,11]. The regional ecological environment on both sides of the highway has been greatly improved by the fixation of the shifting sand dunes with the artificial shelterbelt [12]. However, saline groundwater irrigation may aggravate soil salinization in future and is harmful for shelterbelt plants [13,14]. Therefore, it is necessary to study the water and salt transport of shelterbelt soils which is helpful for the sustainable utilization and management of TDHS [15].

SWCC, as a component of the unsaturated soil mechanics framework, provides the information needed to characterize the properties of unsaturated soils [16], is an interpretation of the basic constitutive relationships of unsaturated soil phenomena [17], and is an important tool for studying the properties of unsaturated soils [18]. SWCC is one of basic hydraulic properties that simulate water and solute transport under unsaturated conditions [19]; is universally used in agriculture [20], soil physics [21], soil chemistry [22], mineralogy research [23], geotechnical engineering [24]; and is widely used in the soil–plant–atmosphere continuum (SPC) [18,25] and other fields. Due to its importance in soil hydrodynamics and solute transport modeling, many SWCC models, both numerical and theoretical, have been developed [21]. A good SWCC model should have simple and clear parameters and be easy to use. It can satisfy the three characteristics of accuracy, universality and simplicity as much as possible at the same time [26]. The VG model [27], BC model [28] and G model [29] have relatively few parameters, and can accurately describe the SWCC of various soil textures [30,31,32,33,34,35,36,37,38]. SWCC is influenced by soil texture, bulk density, organic matter, salinity, temperature, etc. [19,30,31,32,33,39,40,41,42]. Therefore, the fitting results of the models are often different in various study areas. For example, some scholars [43,44] pointed out that the G model can accurately fit SWCC; however, others [45] came to the opposite conclusion, pointing out that the G model cannot accurately fit SWCC. Matlan et al. [26] compared four models and found that the BC model has the most accurate description of the SWCC of sandy soils. Li et al. [46] also proved this, pointing out that the BC model is more suitable than the VG model on soils with high sand content and low clay content. For the SWCC of aeolian sand covered by biocrust, the fitting effect of the VG model is better than the BC model [47].

The literature on SWCC in arid and semi-arid regions is relatively limited as compared to farmland and forest ecosystems. The influence of artificial shelterbelts and long-term saline irrigation on SWCC of sandy soils and the applicability of SWCC models is still unclear. This hinders our understanding towards the soil water-holding capacity and water availability of shelter forests. In this study, we firstly assumed that artificial shelterbelts and long-term saline irrigation have impacts on the water-holding capacity and water availability of different soil layers. In addition, we assumed that the VG model, BC model, and G model have different accuracies in fitting SWCC. Therefore, our study collected soils from 0 to 30 cm layers under the shelterbelt (US), bare land in the shelterbelt (BL) and shifting sandy land (CK) in the hinterland of the Taklimakan Desert, and their SWCCs were determined by the centrifugal method. Combined with its bulk density, organic matter, salinity and other properties, SWCC, pore distribution and soil moisture were analyzed, and the VG model, BC model and G model were used and compared to fit the SWCCs. We aim to reveal the impact of artificial shelterbelt construction on SWCC on shifting sandy soil in extreme arid deserts under saline drip irrigation, so as to provide a basis for desert shelterbelt construction and sustainable management.

2. Materials and Methods

2.1. Study Area

The sampling area is located at the Taklamakan Desert Research Station of the Chinese Academy of Sciences in the hinterland of Taklamakan Desert (39°01′ N, 83°36′ E, 1100 m.a.s.l.), Xinjiang, China (Figure 1). The study area belongs to a warm temperate arid climate. The annual average temperature is 12.4 °C. December is the coldest month with an average monthly temperature of −8.1 °C and July is the hottest month with an average monthly temperature of 28.2 °C. Annual precipitation is 24.6 mm, average relative humidity is 29.4% and annual potential evaporation is up to 3638.6 mm. Annual average wind speed is 2.5 m·s⁻¹, and the maximum instantaneous wind speed is up to 20 m·s⁻¹. The soil type is Xeric Quartzipsamments [48] (Soil Survey Staff 2014), derived from shifting eolian sand, and the basic properties are shown in

Table 1. Basic physiochemical properties of the collected soil samples.

Soil Depth (cm)	Location	Bulk Density (g/cm3)	EC (µS/cm)	pH	Organic Carbon (g/kg)	Soil Mechanical Composition (%)			Soil Texture
						Sand (0.05–2 mm)	Silt (0.002–0.05 mm)	Clay (0–0.002 mm)
0–5	CK	1.55 a	451 c	9.39 a	0.208 c	89.05	2.83	8.12	Sandy
	BL	1.43 b	19,215 a	8.73 a	1.765 b	82.38	9.01	8.61	Loamy sand
	US	1.25 c	10,175 b	9.15 a	3.222 a	65.33	25.46	9.21	Sandy loam
5–10	CK	1.49 b	762 c	8.69 a	0.379 a	89.34	2.73	7.93	Sandy
	BL	1.53 a	4340 b	8.8 a	1.15 a	86.92	4.81	8.27	Loamy sand
	US	1.44 c	6680 a	8.95 a	1.715 a	87.94	4.36	7.7	Loamy sand
10–20	CK	1.47c	2455 a	8.62 b	0.788 a	90.17	1.84	7.99	Sandy
	BL	1.51a	854.5 c	8.71 b	1.058 a	86.09	4.63	9.28	Loamy sand
	US	1.5b	1157 b	9.6 a	1.182 a	88.48	3.34	8.18	Loamy sand
20–30	CK	1.47 b	3240 a	8.82 c	0.83 ab	87.97	4.87	7.16	Sandy
	BL	1.6 a	265 c	9.52 b	0.412 b	87.49	2.7	9.81	Loamy sand
	US	1.44 c	951.5 b	10.11 a	0.912 a	88.09	1.31	10.6	Loamy sand

Note: CK: shifting sandy land; BL: bare land without vegetation cover in the shelterbelt; US: under the shelterbelt. Different lowercase letters represent the significant differences between different treatments (p < 0.05). (Soil texture was classified according to USDA standards based on actual measured soil mechanical composition).

.

2.2. Sample Collection and Determination

In August 2020, undisturbed soil samples were collected with a cutting ring of 100 cm³ from the soil layers of 0–5 cm, 5–10 cm, 10–20 cm and 20–30 cm under US, BL and CK, and all samples had three replicates (the pictures of the sampling points are shown in Figure 2). The ring knife (with soil) samples were soaked in deionized water for 24 h until saturation, weighed, and then the SWCCs were determined using a high-speed centrifuge (CR 22G III model, Hitachi, Japan). Suction values were determined in lab using a centrifuge with different speeds and time settings. The selected suction values were 10.2 cm (310 r/min, 10 min), 30.6 cm (540 r/min, 12 min), 51 cm (690 r/min, 17 min), 71.4 cm (820 r/min, 21 min), 102 cm (980 r/min, 26 min), 204 cm (1200 r/min, 28 min), 612 cm (2190 r/min, 49 min), 1020 cm (3100 r/min, 57 min), 4080 cm (6200 r/min, 77 min) and 8160 cm (8770 r/min, 87 min), respectively. The weight of the ring knife sample was weighed after the completion of centrifugation at each suction value. After centrifugation, the ring knife was dried in an oven (105 °C) and then weighed, and the water content at the corresponding suction value was obtained from the difference in weights. The moisture under different suction values was plotted according to the SWCC. Meanwhile, soil samples of corresponding layers were also collected with a soil drill and brought to the laboratory and air dried. Part of the samples was passed through a 2 mm sieve to determine soil pH, EC (with an EC500 pH/Conductivity Meter (ExStik, Boston, MA, USA)), and soil mechanical composition was determined using the hydrometer method and soil texture was classified according to USDA standards. Soil particles were divided into sand particles of 2–0.05 mm, silt particles of 0.05–0.002 mm and clay particles of ≤0.002 mm, and the left soils were used to measure soil organic matter content with potassium dichromate external heating method after passing through a 0.25 mm sieve. The physical and chemical properties of the soil samples are shown in Table 1.

2.3. Calculations of Soil Equivalent Pore Sizes and Moisture Constants

The equivalent pores are divided into six levels based on their diameters, including narrow micropores (≤0.3 μm), micropores (0.3–5 μm), fine pores (5–30 μm), medium pores (30–75 μm), macropores (75–100 μm) and interstices (≥100 μm) [47,49].

Based on the measured and fitted SWCC parameters, saturated water θ_s, field capacity θ_f, wilting coefficient θ_r, available water content θ_a and readily available water content θ_ra were obtained. While θ_s is the soil water content when water suction is zero, θ_f is the soil water content of when pF = 1.8 (pF is expressed as the logarithm of the centimeter height of the water column of the soil water potential), and θ_r is the soil water content when pF = 4.2, pF = 3.8 represents temporary wilting coefficient. The water content in the range of pF 1.8–4.2 is θ_a, and the water content in the range of pF 1.8–3.8 is θ_ra [50].

2.4. SWCC Modeling

VG model, G model and BC model were used to describe the SWCC.

VG model:

$${\theta }_{(h)}={\theta }_r+\frac{{\theta }_s-{\theta }_r}{{(1+{\left|\alpha h\right|}^n)}^m}$$

(1)

where θ_(h) is the volume water content, cm³ cm⁻³; θ_s is the saturated water content, cm³ cm⁻³; θ_r is the residual water content, cm³ cm⁻³; h is the matric suction, cm; α is a scale parameter that is related to the inverse of the air entry suction, cm⁻¹; m, n are the fitting parameters; m = 1 − 1/n (n > 1).

G model:

$$\theta =\alpha h^{{}^{-b}}$$

(2)

where θ is the volume water content, cm³ cm⁻³; h is the matric suction, cm; a and b are fitting parameters.

BC model:

$$S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}=\left\{\begin{array}{ll}{\left(\alpha h\right)}^{-\lambda }& \alpha h>1\\ 1& \alpha h\le 1\end{array}\right.$$

(3)

where S_e is the effective saturation; λ is the curve shape parameter; the physical meanings of other parameters are the same as above.

2.5. SWCC Model Fitting Accuray Assesments

The coefficient of determination (R²), root mean square error (RSME) and relative error (RE) were used to quantitatively evaluate the fitting effect of the models. Grey correlation analysis allows ranking the importance between different influencing factors [51,52]. The soil physiochemical properties and model parameters were quantitatively analyzed by grey correlational method, and the correlational degree was sorted. The calculation formulas are as follows:

$$R^2=\frac{{\displaystyle \sum _{i=1}^N{({\theta }_i-\overline{{\theta }_i})}^2}}{{\displaystyle \sum _{i=1}^N{({\beta }_i-\overline{{\theta }_i})}^2}}$$

(4)

$$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^N{({\theta }_i-{\beta }_i)}^2}}{N}}$$

(5)

$$RE=\frac{\left|{\theta }_{\mathrm{i}}-{\beta }_i\right|}{{\theta }_{\mathrm{i}}}\times 100$$

(6)

$${\gamma }_{0i}=\frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{m+\rho M}{{\Delta }_{\mathrm{i}}(k)+\rho M}}$$

(7)

In Equations (4)–(6), N is the total number of samples of matric suction; ${\theta }_i$ represents the measured value of soil moisture corresponding to the ith pressure value; $\overline{{\theta }_i}$ represents the average value of measured soil moisture; ${\beta }_i$ is the fitted value of soil moisture corresponding to the ith pressure value. In Equation (7), ${\gamma }_{0i}$ is the correlation degree, Δ_i(k) is the difference sequence, M is the maximum difference sequence, m is the minimum difference sequence, and $\rho$ is the resolution coefficient, which is generally 0.5 in the models [51,52].

RETC software was used to solve and fit the parameters, Excel2010 and Origin2018 were used for data processing and mapping and SPSS18.0 was used for one-way ANOVA and multiple comparisons (a = 0.05, LSD).

3. Results

3.1. Effects of Artificial Shelterbelt Construction on Soil Physiochemical Properties

Soil physiochemical properties were improved after artificial shelterbelt construction. From Table 1, according to USDA system, it is clearly stated that all soil layers of CK are sandy soil, all soil layers of BL treatment are loamy sandy soil, 0–5 cm soil of US treatment was sandy loam soil and 5–30 cm soil was loamy sandy soil. There were significant differences in bulk density and EC among CK, BL and US in the same soil layer (p < 0.05). Except for CK, the soil bulk density under BL and US treatments increased with the increasing soil depth, and the soil bulk density of 0–5 cm under BL and US treatments was significantly lower than CK; the bulk density of the 5–30 cm soil layer was the highest in BL, except for 10–20 cm; the bulk density of US treatment was lower than that of CK and BL. The bulk density of the 0–5 cm soil layer in US treatment was the smallest (1.25 g/cm³), and the highest bulk density at all was 1.6 g/cm³. With the increase in soil depth, the EC of CK gradually increased, while EC of BL and US treatments decreased with the increasing depth. The difference in EC between BL, US and CK in the 0–5 cm soil layer was the largest, which was 18,764 µS/cm and 9724 µS/cm, respectively. The EC of CK at 10–30 cm was higher than that of BL and US.

3.2. Screening of Soil Water Characteristic Curve Models

As listed in

Table 2. R² and RMSE of different fitting models.

Location	Soil Depth (cm)	VG Model		BC Model		G Model
		R2	RMSE (%)	R2	RMSE (%)	R2	RMSE (%)
CK	0–5	0.994	0.0102	0.997	0.0079	0.919	0.0682
	5–10	0.996	0.0082	0.998	0.0071	0.906	0.0862
	10–20	0.998	0.0150	0.996	0.0090	0.917	0.0710
	20–30	0.998	0.0631	0.998	0.0059	0.908	0.0636
BL	0–5	0.987	0.0165	0.977	0.0259	0.944	0.0395
	5–10	0.997	0.0067	0.992	0.0117	0.884	0.0787
	10–20	0.997	0.0078	0.998	0.0054	0.935	0.0704
	20–30	0.995	0.0100	0.997	0.0082	0.948	0.0431
US	0–5	0.996	0.0086	0.998	0.0050	0.934	0.0765
	5–10	0.993	0.0089	0.984	0.0147	0.934	0.1133
	10–20	0.997	0.0073	0.995	0.0095	0.913	0.0822
	20–30	0.997	0.0060	0.998	0.0048	0.908	0.1576
	mean	0.9955	0.0140	0.9941	0.0096	0.9206	0.0792

, R² of the fitting values of VG, BC, and G models was ranged between 0.884 and 0.998, which showed the correlations of three models with the measured data were high. As shown in Figure 3, the fitted results of the G model were all higher than the measured points, and the fitted values of the VG model and BC model for BL in the 5–20 cm soil layer were smaller than the measured values.

The fitting effect of the VG model and BC model for different soil layers of various treatments were always better than that of the G model. The fitting effect of the VG model and BC model was similar for the 10–30 cm soil layers of each treatment. The fitting effect of the BC model for the 0–10 cm soil layer of CK was better than the VG model, while the results of BL were opposite. For US, the BC model had a good fitting effect for the 0–5 cm soil layer, but the fitting effect was contradictory for the 5–10 cm soil layer. For all treatments, the fitting errors of three models increased with the increasing water potential and tended to be stable (Figure 4). The relative errors of the VG model and BC model for different treatments and soil layers were always lower as compared with relative error of the G model.

In addition, R² and RMSE of the VG and BC models showed that the simulating results of CK were better than those of BL and US (Table 2). The results indicated that the increase in salinity and organic matter content may affect the fitness of each model.

Table 3. Parameters in the modeling of soil water characteristic curves.

Location	Soil Depth (cm)	VG Model				BC Model				G Model
		θr (%)	θs (%)	a (cm−1)	n	θr (%)	θs (%)	a (cm−1)	λ	a	b	a × b
CK	0–5	0.048	0.411	0.032	2.841	0.044	0.405	0.044	0.728	0.375	0.728	0.273
	5–10	0.044	0.409	0.036	2.805	0.041	0.402	0.047	0.663	0.373	0.663	0.247
	10–20	0.041	0.409	0.033	2.428	0.036	0.400	0.053	0.708	0.379	0.708	0.268
	20–30	0.042	0.409	0.033	2.690	0.039	0.403	0.046	0.733	0.383	0.733	0.281
BL	0–5	0.063	0.430	0.022	2.670	0.043	0.421	0.040	0.990	0.366	0.990	0.362
	5–10	0.047	0.355	0.022	3.291	0.040	0.355	0.037	0.710	0.363	0.710	0.258
	10–20	0.041	0.423	0.044	2.302	0.036	0.409	0.062	0.711	0.379	0.711	0.269
	20–30	0.047	0.467	0.044	2.257	0.040	0.452	0.062	0.826	0.381	0.826	0.315
US	0–5	0.080	0.461	0.031	2.160	0.073	0.451	0.045	0.774	0.290	0.774	0.224
	5–10	0.038	0.347	0.040	2.074	0.016	0.343	0.101	0.616	0.352	0.616	0.217
	10–20	0.051	0.411	0.034	2.295	0.046	0.403	0.049	0.711	0.339	0.711	0.241
	20–30	0.045	0.333	0.035	2.552	0.043	0.326	0.047	0.521	0.328	0.521	0.171

Note: θ_s: saturated water content; θ_r: wilting coefficient; α is a scale parameter that is related to the inverse of the air entry suction; n is the fitting parameter; λ is the curve shape parameter; a and b are fitting parameters; a × b is the specific water capacity when the soil water suction is 1 Bar.

lists the parameters of the three models in SWCC modeling. From grey correlation calculation (

Table 4. Correlation degree analysis between VG model parameters and soil basic physiochemical parameters.

Correlation Degree	pH	EC (µS/cm)	Organic Carbon (g/kg)	Bulk Density (g/cm3)	Sand (0.05–2 mm)	Slit (0.002–0.05 mm)	Clay (0–0.002 mm)
a	0.925	0.684	0.783	0.930	0.926	0.744	0.917
n	0.939	0.694	0.794	0.952	0.948	0.781	0.922

), it could be concluded that VG model parameters a and n have a higher degree of correlation with the soil physiochemical parameters (BC and G models were not included, because the VG model had the best fitting effect). The grey relational degrees between pH, EC, organic matter, bulk density, sand content, silt content, clay content and model parameters (a and n) were all larger than 0.6, and the order of correlation showed as: bulk density > sand content > pH > clay content > organic matter > silt content > EC.

3.3. Effects of Artificial Shelterbelt Construction on Soil Water Retention Performance

Soil porosity was significantly changed after artificial shelterbelt construction. As listed in Figure 5, few micropores were found under all treatments. Compared with CK, the contents of interstices and macropores in the 0–10 cm soil layer under BL and US were much lower, and the contents of fine pores and micropores were higher. The medium pores of the BL increased the most, with 0–5 cm increased by 3.74% and 5–10 cm increased by 4.94%. Fine pores of US increased the most: 0–5 cm increased by 4.54% and 5–10 cm increased by 2.32%.

Soil water parameters of each treatment are shown in

Table 5. Soil water parameters calculated from the VG models.

Depth	Location	θs (%)	θf (%)	θr (%)	θa (%)	θra (%)
0–5	CK	40.92	13.89	4.75	9.13	9.13
	BL	42.73	23.46	6.31	17.15	17.14
	US	46.11	23.67	8.05	15.61	15.56
5–10	CK	40.67	12.21	4.40	7.80	7.80
	BL	35.80	16.46	4.74	11.72	11.72
	US	34.32	14.39	3.82	10.56	10.51
10–20	CK	40.36	15.90	4.06	11.84	11.83
	BL	41.41	13.65	4.13	9.53	9.51
	US	41.20	17.51	5.14	12.38	12.35
20–30	CK	40.57	14.07	4.24	9.83	9.83
	BL	45.67	15.84	4.70	11.14	11.12
	US	39.59	12.38	4.54	8.59	8.52

Note: θ_s: saturated water content; θ_f: field capacity; θ_r: wilting coefficient; θ_a: available water content; θ_ra: readily available water content. (BC and G model were not included because the VG model had the best fitting effect).

. In the 0–5 cm soil layer, θ_s, θ_f, θ_r, θ_a and θ_ra under BL and US were higher than those under CK. Compared with CK, θ_s of BL and US increased by 4.42% and 12.67%, θ_f increased by 68.9% and 70.41%, θ_r increased by 32.84% and 69.47%, θ_a increased by 87.84% and 70.97% and θ_ra increased by 87.73% and 70.43%, respectively. In the 5–10 cm soil layer, θ_s of CK was the largest, but θ_f, θ_r, θ_a and θ_ra of US and BL were higher than those of CK. In the 10–20 cm soil layer, higher θ_f, θ_r, θ_a and θ_ra were observed under US than in BL and CK. In the 20–30 cm soil layer, higher θ_s, θ_f, θ_r, θ_a and θ_ra were observed under BL than in US and CK. Meanwhile, the contents of available water in each treatment were almost equal to that of readily available water content.

As shown in Figure 6, the shapes of SWCC of each soil layer were similar, but the variations in soil moisture under per unit suction were clearly different. Under 0–10 cm and 1000–10,000 cm suctions, soil water was lost slowly with the increase in suction. However, under the suction of 10–1000 cm, the curve trended to be steep, and the soil water decreased rapidly with the increase in suction. Figure 6 clearly showed the differences in the course of the water retention curves of the three sampled areas. Soil water-holding capacity of 0–5 cm soil layer was highest under US, followed by BL and CK (Figure 6a). The difference in water-holding capacity among treatments in the remaining soil layers gradually became smaller. At the suction value corresponding to the field water-holding capacity (p_F = 1.8, i.e., 63 cm water column), the moisture of the 5–10 cm soil layer was as follows: BL ≥ US≥ CK; 10–20 cm soil layer was as follows: US ≥ CK ≥ BL; and 20–30 cm soil layer was as follows: BL ≥ CK ≥ US. When reaching the suction value corresponding to the temporary wilting coefficient (p_F = 3.8, i.e., 6309 cm water column), the water content of the 5–10 cm soil layer behaved as follows: BL ≥ CK ≥ US; the water content of the 10–20 cm soil layer behaved as follows: US ≥ BL ≥ CK; and the water content of the 20–30 cm soil layer behaved as follows: BL ≥ CK ≥ US.

4. Discussion

4.1. Artificial Shelterbelt Construction Greatly Changed the Soil Physiochemical Properties

SWCC is affected by various soil properties such as texture, bulk density, porosity, organic matter and salinity [42,53,54,55]. In our study, compared with shifting sandy land, the soil properties of BL and US changed obviously (Table 1), which caused the transformation of SWCC (Figure 6). Soil physiochemical properties of the 0–5 cm soil layer were changed most significantly. Shelterbelt construction under saline irrigation significantly decreased the soil bulk density, increased EC and organic matter content, and the soil texture changed from sandy soil to loamy sand and sandy loam (Table 1). After shelterbelt construction, the soil bulk density of 0–5 cm decreased significantly as compared with CK, which is primarily due to the continuous input of plant litters in the surface soil that lead to loosening the soil particles [56,57]. The bulk density of BL was the highest in the 5–30 cm soil layer, indicating that long-term saline irrigation would increase the soil bulk density [58,59]. The change in soil texture was mainly reflected by the increase in silt and clay particles. The main reason is that the shelterbelts reduced the wind speed, which promotes the precipitation of sand and the accumulation of dust fall [60,61]. The second reason is the accumulation of litter and the role of microorganisms. The volume and average radius of soil macropores increased with the increase in volumetric rock fragment content [62,63]. This directly led to the reduction in large pores and increase in small pores under BL and US (Figure 5a), which improved soil water-holding capacity and changed the SWCC. In addition to the influence of rock fragments, soil macropores are also controlled by biological factors [64]. Therefore, when the content of rock fragments is similar, the distribution of soil macropores will also be different (Figure 5b–d).

Under long-term saline irrigation, the salts were added into the soil and resulted in surface accumulation due to strong evaporation rate [13,65]. Therefore, the soil ECs of BL and US in the 0–10 cm soil layer were significantly higher than shifting sandy land, and the EC of BL and US decreased with the increasing soil depth. Vegetation cover resulted in less evaporation under US than BL, and more intense surface salt accumulation under BL. Therefore, soil EC of BL was higher than US in the 0–5 cm soil layer, and EC of BL was lower than US in the 0–10 cm soil layer. Wind erosion has an important impact on the cycle of soil organic carbon. The fine particles in the sand and dust adsorb soil organic carbon. Under the action of wind erosion, soil organic carbon is redistributed along with the movement of sand and dust [65,66]. Therefore, the increase in soil organic matter was mainly due to the accumulation of litters, and atmospheric dustfall also played a certain role in promoting it [67].

4.2. Artificial Shelterbelt Construction Increased the Soil Water-Holding Capacity

Soil water moves in pores and its transfer rate is directly determined by the size and distribution of pores. Soil bulk density is negatively correlated with soil porosity. Changes in soil primary particles such as sand, silt and clay also affect the distribution of pores. Organic matter and salinity have a direct impact on soil structure and adsorption. These soil physical properties may directly or indirectly affect soil water conductivity. Therefore, SWCC is affected by bulk density, texture, organic matter, porosity, aggregate stability, salinity and other properties [19,30,31,32,33,40,41,42].

Compared with CK for the 0–5 cm soil layer, vegetation coverage resulted in the decline in bulk density and increase in the salinity, organic matter and silt content under BL and US. On the one hand, salt contents in soil will occupy the pore space, and it will cause some soil particles to flocculate together, increase soil pores and enhance soil water-holding capacity [68]. Plant litters will reduce soil bulk density and increase soil organic matter content, soil saturated water content and water conductivity [69,70]. Therefore, the water contents of BL and US were higher than those of CK under the same suction, and with the increase in suction, the water contents of BL and US decreased, which means the SWCC integrally moved upwards and the trend slowed down. The contents of soil interstices and macropores in the 5–10 cm of CK were more than that in BL and US (Figure 5b). SWCC reflects the dehumidification process of soils, during which water is stored in pores and water retained in macropores is preferentially expelled as suction increases. Therefore, the change in water content per unit suction of CK was higher than that of BL and US. Therefore, the curve was the maximum under 0–10 cm suction, and then decreased rapidly to the minimum. The soil physiochemical properties of each treatment in the 10–20 cm soil layer had little difference, and their water-holding capacities were similar. For the 20–30 cm soil layer, the bulk density of BL was significantly higher than that of CK and US, and the water-holding curve was at its peak. Our results are consistent with Lipiec et al. [54]. Comparing the soil water changes over the entire suction range, the water-holding capacity of the surface 0–5 cm soil increased the most significantly, and the water-holding capacity of BL and US increased over the entire suction range compared to CK. The increase in water retention capacity of shelterbelt soils under long-term saline water irrigation is mainly due to the following reasons. Soil texture controls the physical, hydrological and chemical properties of the soil and has a strong influence on water flow paths, residence times and the magnitude and location of salt accumulation. Soil texture governs the water and solute transport. Under irrigation conditions, finer soils limit water infiltration, and coarse-grained soils retain significantly less water than fine-grained soils [71,72]. The accumulation of soil salts is dominated by sodium salts originated from the irrigation water. Excessive concentration of sodium ions in the soil solution will disperse and swell the soil structure, leading to the reduction and blockage of connected pores, thus reducing the permeability and hydraulic conductivity of the soil [73,74]. The salinity of pore water also affects the development of the diffusion double layer (DDL) around soil particles, which controls the microstructural changes in soil particles during hydration. As the salt content in irrigation water increases, the interlayer space between DDLs expands and soil aggregates are disrupted by swelling and clay dispersion [41,75]. The more significant salinity damages the soil structure and the salt stress results in higher absorption of soil water [76].

4.3. Screening of SWCC Models for Artificial Shelterbelt

The fitting R² results of VG model and BC model were significantly higher than those of the G model, and the RE was generally lower than that of G model. Therefore, the fitting results of the VG and BC models were better for SWCC, while the BC model had better fitting effect for CK. However, the fitting effects of the VG model were better than the BC model for bare BL and US, especially for the surface soil. This is because the VG model considers more influencing factors when predicting. Therefore, it has a higher accuracy on soils with more complex physiochemical properties. The soil salt content and organic matter content of BL and US were maximum, so the fitting effect of the VG model is better than the BC model. Therefore, comparing the fitting results of the three models, the results concluded that the VG model is the best choice in this regard.

When using RETC software to predict the parameters of the VG model, only the influence of bulk density and texture can be considered. In our study, it can be found that the increase in organic matter and salinity reduce the fitting accuracy of the VG model, and through grey correlation analysis, we found that the soil bulk density, sand content, pH, clay content, organic matter, silt content, EC and other physiochemical indicators of the study area had a larger correlation degree with the parameters a and n of the VG model (Table 4). Meanwhile, the basic parameters such as organic matter, pH and EC were used as input variables to predict the parameters of the VG model, and our results found that the fitting effect of the VG model could be improved [77,78]. Therefore, pH, EC, organic matter and other indicators should also be used as input variables for calculation when predicting the parameters of the VG model.

5. Conclusions

The construction of an artificial shelterbelt with long-term saline water irrigation increased the water retention capacity and the content of fine pores in the soils, and reduced the content of macropores. Artificial shelterbelt construction had the greatest impact on the surface soil, and the impact gradually decreased with the increasing soil depth. Compared with the shifting sandy land (CK) in the 0–5 cm soil layer, the saturated water content θ_s of the soil under BL and US increased by 4.42% and 12.67%, the field capacity θ_f increased by 68.9% and 70.41%, and the available water content θ_a increased by 87.84% and 70.97%. In the 5–10 cm soil layer, the θ_f, θ_r, θ_a and θ_ra of US and BL were higher than those of CK. In the 10–20 cm soil layer, the soil of US had the best water-holding performance. In the 20–30 cm soil layer, the soil water-holding performance of BL was the best. The available water content of each treatment was in the form of readily available water content. Comparing the coefficient of determination (R²), root mean square error (RSME) and relative error (RE) of the three models, it was found that the G model always overestimated the soil water content and had a lower prediction accuracy, while the BC and VG models had higher prediction accuracies. Although both BC and VG models are suitable for fitting SWCC of the shelterbelt, the VG model is more effective. The parameters a and n of the VG model had a higher degree of correlation with the soil EC and organic matter. In summary, the increase in organic matter and salinity reduced the fitting accuracy of the model. When predicting the parameters of the VG model and establishing the soil transfer function, soil EC, organic matter and other indicators should be calculated as input variables.