Interacting Effects of Land Use Type, Soil Attributes, and Environmental Factors on Aggregate Stability

Abstract

Soil erosion and surface pollution near reservoirs can adversely affect water quality and safety. Soil aggregate stability is an important predictor of soil water loss and erosion resistance that is strongly influenced by land use. This study therefore aimed to identify factors affecting soil aggregate stability near reservoirs to provide empirical and theoretical insights that could guide the development of management measures to increase land quality, optimize land use, and maximize sustainability. This study focuses on the land around the Shitoukoumen Reservoir in China and examines the effects of six land use types, eleven soil physicochemical properties, and five environmental factors. Ninety-four sets of soil samples were collected in 2021 for analysis of soil aggregates and properties. Particle size classification of soil aggregates was carried out using the wet sieve method and four indicators were calculated to evaluate the effects of land use, soil physicochemical properties, and environmental factors on soil aggregate stability: water stable aggregates (WSA), mean weight diameter (MWD), geometric mean diameter (GMD), and fractal dimension (D). Descriptive statistics and geostatistics were used to explore the spatial distributions of soil aggregate stability around the reservoir and the influence of soil properties was studied using correlation analysis and path analysis. The conclusion indicates that land use type significantly affects aggregate stability. The most stable aggregates were found in paddy fields (WSA = 0.77, MWD = 0.76, GMD = 0.57) and forests (WSA = 0.75, MWD = 0.76, GMD = 0.55), followed by an orchard, irrigated land, and grassland. Aggregate stability was worst in upland sites (WSA = 0.61, D = 2.28), where soil aggregates were highly fragmented. There were clear spatial correlations between all four stability indicators. The environmental factors and soil physicochemical characteristics with the strongest influence on aggregate stability were soil organic matter, pH, soil clay content, total nitrogen, and temperature changes. Path analysis revealed that some soil properties affect aggregate stability indirectly, with particularly complex relationships between clay, soil organic matter, and pH. In conclusion, land use type, soil organic matter, pH, soil clay content, total nitrogen, these soil physicochemical properties, and environmental factors, especially temperature, significantly affect soil aggregate stability around reservoirs. In the future, it is necessary to appropriately change upland into paddy land, increase forest land, and appropriately add organic fertilizer to improve soil quality.

1. Introduction

Soil aggregates are structural units <10 mm in diameter that are formed from soil particles under natural conditions [1]. They strongly affect the structure and fertility of the soil, so their formation is essential for maintaining the quality of arable land. The properties of soil aggregates depend on the properties of the soil-forming parent material as well as the soil fauna, microorganisms, plant roots, and inorganic binding agents that are present. Their basic structure develops gradually as fresh organic matter enters the soil from plant roots and is decomposed by animals and microorganisms, becoming an effective binder of mineral soil components [2]. Additionally, they provide habitats for animals and microorganisms, and their presence increases the suitability of arable land for cultivation [3,4,5]. Soil erosion is a major environmental problem worldwide and an important constraint on food security, ecological optimization, and social sustainability [6]. The formation of soil aggregates enhances soil structural stability and thus increases resistance to soil erosion. Different degrees of erosion in each land use type lead to differences in the stability of soil aggregates: land types with no artificial disturbance and more complex vegetation are more stable [7]. Therefore, favorable soil agglomeration characteristics facilitate the preservation of high-quality arable land. Previous studies have shown that overuse of agricultural land significantly affects soil structure and composition. This in turn reduces the supply of nutrients for plants and soil microorganisms, which accelerates soil erosion and reduces crop yields [8]. Maintaining a favorable soil structure with high aggregate stability is therefore essential for improving land fertility and achieving sustainable agriculture.

Analyses of soil aggregate stability commonly divide aggregates into four size classes that are separated by wet sieving: large macro-aggregates (>2 mm), small macro-aggregates (0.25–2 mm), micro-aggregates (0.053–0.25 mm), and the silt and clay fraction (<0.053 mm) [9,10]. Indicators of stability used in such analyses include soil water-stable aggregates (WSA) and the mean weight dimension (MWD), geometric mean diameter (GMD), and fractal dimension (D) of the aggregates [11,12]. Studies on factors affecting soil aggregate stability have investigated the influence of soil carbon pools, nitrogen pools, and land management practices. Xu et al. [13] suggested that organic carbon retention has a stronger effect on soil aggregate stability than the abundance of recalcitrant carbon and argued that reducing soil disturbance will be vital for future forest management. Additionally, Zhong et al. [14] suggested that vegetation restoration promotes the physical protection of soil organic carbon by increasing the abundance of soil aggregates, and that the formation of macro-aggregates is important for maintaining steady soil carbon levels. Finally, Ameer et al. [15] reported a strong correlation between soil organic carbon and soil aggregates, leading them to conclude that fallow periods can enhance soil fertility. The stability of soil aggregates depends on both environmental factors and the soil’s physicochemical properties, and the results of previous studies suggest that this fertility enhancement may be primarily due to changes in the organic carbon content that indirectly affect aggregate stability [16]. It is also well established that plant roots act as cementing agents that promote the aggregation of fine soil particles into stable macro-aggregates by adding organic matter to the soil [17]. In addition, temperature and moisture condition changes can influence soil aggregate stability [18]. Previous studies have also found that elevated rainfall contributes to the stability of soil aggregates [19]. Dramatic temperature changes (freeze–thaw cycles) also have a significant impact on soil aggregate stability [20]. Understanding the interplay between these factors will be crucial for maintaining soil fertility and enhancing agricultural productivity.

Reservoirs are areas exposed to severe soil erosion. Soil conditions in their surrounding areas are vital for sustainable management of water resources and maintenance of local ecological stability, and also play central roles in regional development. China has more reservoir dams than any other country in the world [21]. Reservoirs, formed by reservoir construction, are characterized by complex topographic conditions, concentrated rainfall, significant regional temperature, distinct soil properties, and severe soil erosion [22]. They also play a key role in the functioning of local ecosystems. Consequently, the impact of human activities on soil structure in reservoir areas has been studied extensively [23]. A deep understanding of the effects of soil properties, land use types, environmental factors, and their interactions on soil aggregate stability is vital for understanding soil quality and erosion resistance, and also for predicting water loss from land in ecologically sensitive sites such as reservoir areas. However, most studies on factors affecting soil aggregate stability in such areas have investigated the influence of individual factors using traditional statistical methods or spatial correlation analysis [18]. There have been few analyses of reservoir-scale spatial correlations, the effects of interactions between different factors, and the degree to which different factors act directly or indirectly on aggregate stability.

In this study, we hypothesize that soil aggregate stability changes periodically as a dynamic property of soil under the control of various mechanisms. Soil aggregate stability is affected by a combination of soil properties, natural factors, and human factors. Natural and human factors can indirectly affect soil aggregate stability by changing soil environmental properties, and soil itself can also indirectly affect soil aggregates through changes in other soil properties. The study’s objectives were to (a) clarify the spatial distribution of soil aggregate stability and land use in the Shitoukoumen reservoir area; (b) investigate the effect of land use on soil aggregate stability; and (c) identify the driving factors affecting soil aggregate stability. The results obtained will provide valuable support for future efforts to improve the quality of arable land in the reservoir area, optimize patterns of land development and protection, promote rational use of land resources in the reservoir area, and achieve sustainable development. In addition, they will provide a useful theoretical and practical reference for similar efforts in other reservoir areas.

2. Materials and Methods

2.1. Study Area and Soil Samplings

The Shitoukoumen Reservoir watershed is located near Changchun City in Jilin Province, northeastern China (between 43°06′~44°52′30″ N and 125°44′40″~126°18′20″ E) and covers an area of 8255 km², excluding the Yitong River. The soil types are mainly dark brown loam and albic soil, which are suitable for cultivation. Moreover, this area is in the midtemperate continental semihumid monsoon climate zone, with a large temperature difference between cold and summer seasons. It is a major source of drinking water for Changchun City, which has over nine million residents. The reservoir itself is located in the middle reaches of the Yinma River, which is a tributary of the Songhua River. Drastic hydrothermal changes in the area have caused frequent droughts and floods that have severely affected the local soil’s physicochemical properties (Soil pH, organic matter, total nitrogen, total phosphorus, and total potassium range from 4.1 to 7.3, 8.7 to 100.4 g/kg, 0.2 to 5.2 g/kg, 0.00004 to 2.6 g/kg and 10.0 to 102.5 g/kg, respectively), and the reservoir’s water quality has fallen in recent years. Reliable information on the stability of the watershed’s soil aggregates is therefore needed to guide efforts to control and prevent further soil erosion.

To investigate the stability of soil aggregates around the reservoir, 94 sampling points with a minimum spacing of 1000 m were chosen randomly in July 2021. Their locations were then overlaid on a land use map for the studied area based on data from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences to determine the land use type at each sampling point. Six land use types were found in the studied area: upland (UP), paddy field (PF), irrigated land (IL), grassland (GL), forest land (FL), and orchard (OR). The land use data come from the third national land resource survey conducted by the Ministry of Natural Resources, China. Soil samples were collected within 1 m of each sampling point at a depth of 0–20 cm. In each point, five soil samples were taken randomly and mixed into one composite sample. These samples, which weighed 1–1.5 kg, were thoroughly mixed before being analyzed. Figure 1 shows the land use types present in the study area at the time of sampling.

2.2. Analysis of Soil Properties

2.2.1. Analysis of Soil Aggregate Properties

Soil aggregates in each sample were separated into >2 mm, 1~2 mm, 0.25~1 mm, 0.05~0.25 mm, and <0.05 mm size fractions and the abundance of each fraction was recorded. Aggregate particle size classification was performed using the wet sieving method (refer to [24]). In this method, an 80 g soil sample is weighed on a 2 mm sieve that is then placed in an aluminum box containing deionized water and submerged for 5 min, keeping the soil below the water level. The sieve is then shaken vertically for about 3~5 min and any material remaining on the sieve is rinsed with deionized water and retained. The filtered water–soil mixture is then subjected to the same procedure using 1 mm, 0.25 mm, and 0.05 mm sieves in succession.

The indices used to evaluate soil aggregate stability in this work were WSA, MWD, GMD, and D. WSA, MWD, and GMD reflect the structural quality of the soil; higher values of these indices indicate greater structural stability, while lower values suggest unstable soil structures and an elevated likelihood of soil nutrient loss. D relates more to the characteristics of the aggregates themselves; lower values of D indicate more stable soil aggregate structures. The equations used to calculate these indices are:

$$\begin{array}{c}\mathrm{WAS}=\frac{M_{i>0.25}}{M_T}\end{array}$$

(1)

$$\begin{array}{c}\mathrm{MWD}=\frac{{{\displaystyle \sum }}_{i=1}^n\left(x_i{w}_i\right)}{{{\displaystyle \sum }}_{i=1}^n{w}_i}\end{array}$$

(2)

$$\begin{array}{c}\mathrm{GMD}=\exp \left[\frac{{{\displaystyle \sum }}_{i=1}^n{w}_i\ln x_i}{{{\displaystyle \sum }}_{i=1}^n{w}_i}\right]\end{array}$$

(3)

$$\begin{array}{c}\mathrm{D}=3-\frac{{\log }_{10}\left(\frac{M_{r<R_i}}{M_T}\right)}{{\log }_{10}\left(\frac{\overline{R_i}}{{\overline{R}}_{max}}\right)}\end{array}$$

(4)

Here, $M_{i>0.25}$ is the weight of aggregates larger than 0.25 mm, $M_T$ is the total mass of aggregates, $w_i$ is the weight percentage of soil aggregates in each aggregate fraction, and $x_i$ is the average diameter of each aggregate fraction [16,25]. $\overline{R_i}$ is the average diameter of the aggregates in the ith aggregate fraction, ${\overline{R}}_{max}$ is the average diameter of the largest aggregate fraction, and $M_{r<R_i}$ is the cumulative mass of aggregates with a particle size below $R_i$ [26]. Soil fractal dimensions are obtained by calculating linear fits for different particle size classes. For more detailed formula information refer to [4].

2.2.2. Soil Properties and Supplemental Data

Soil properties measured in this work include the soil sand content (Sand), soil silt content (Silt), soil clay content (Clay), soil pH, soil organic matter (SOM), total phosphorus (TP), available phosphorus (AP), total nitrogen (TN), available nitrogen (AN), total potassium (TK), and available potassium (AK). The procedures used to measure these properties are detailed in Soil Sampling and Analysis Methods, edited by Li et al. [27].

A DEM was generated using 30 m Shuttle radar mapping data from SRTM project images; elevation (ELE), slope, and aspect data were generated from this the Shuttle data using the ArcGIS10.7 software package. Temperature data were obtained from Landsat8OLI satellite images downloaded from the geospatial data cloud (https://www.gscloud.cn (accessed on 11 February 2023)); the annual average land surface temperature (LST) and seasonal land surface temperature difference (ΔLST) for the area of interest were determined by selecting 11 remote sensing images at 30 m resolution with low cloudiness. The LST data were obtained by inversion of the atmospheric correction method, which is performed in Figure 2.

2.3. Statistical Analysis

A one-way ANOVA test was used to calculate the effect of land use types on the stability of aggregates. Multiple comparisons use the Tukey’ HSD test. Spatial correlation analysis was used to characterize spatial correlations between the observed distributions of soil aggregate stability and land use patterns. Spatial autocorrelation models for WSA, GMD, MWD, and D were established using semivariance functions in GS+ 9.0, and ordinary Kriging interpolation was performed using ArcGIS 10.7. To enhance the spatial correlation of the data, the initial MWD and GMD values were logarithmically transformed.

The effects of soil properties on soil aggregate stability were evaluated first using Pearson correlation analysis and then by using path analysis to clarify the direct and indirect effects of soil properties on aggregate stability. It is worth noting that the accuracy of the interpolation process used in this study requires further discussion, limited by the large ratio of the study area to the number of samples and the uneven distribution of sample points.

Path models are based on multiple linear regression equations that decompose the correlations between independent and dependent variables into direct and indirect effects to clarify the interactions between independent, intermediate, and dependent variables [28]. The calculation procedure used in this work was described in detail by Song et al. [29].

3. Results

3.1. Soil Aggregate Size and Stability

The soil samples were separated into five size fractions: >2 mm, 2~1 mm, 1~0.25 mm, 0.25~0.05 mm, and <0.05 mm (

Table 1. Soil aggregate properties. Multiple comparisons use the Tukey’ HSD test. Different letters indicate significant difference between land use types at p = 0.05 level. PF: paddy field; IL: irrigated land; FL: forest land; UP: upland; OR: orchard; GL: grassland.

Land Use Type	>2 mm	2–1 mm	1–0.25 mm	0.25–0.05 mm	<0.05 mm
PF	13.67 ± 17.28 a	7.39 ± 11.11 ab	56.1 ± 19.76 a	15.81 ± 8.21 b	7.02 ± 6.4 a
IL	6.76 ± 9.83 ab	2.38 ± 1.46 b	62.44 ± 13.87 a	21.73 ± 12.41 ab	6.69 ± 4.2 a
FL	11.19 ± 12.87 a	11.27 ± 12.86 a	53 ± 16.88 a	19.96 ± 8.73 b	4.58 ± 2.81 a
UP	3.33 ± 3.45 b	3.06 ± 2.44 b	55.12 ± 16.5 a	30.69 ± 15.24 a	7.81 ± 7.14 a
OR	7.49 ± 5.68 ab	5.94 ± 2.83 ab	60.95 ± 9.73 a	22.34 ± 9.69 ab	3.73 ± 1.63 a
GL	3.83 ± 2.48 ab	5.85 ± 5.28 ab	64.34 ± 4.62 a	19.75 ± 5.04 b	6.22 ± 3.34 a

).

The 1~0.25 mm size fraction dominated, accounting for 58.66% of the total aggregate mass on average, followed by the 0.25~0.05 mm fraction (21.71% of the total aggregate mass). The <0.5 mm size fraction was least abundant overall, accounting for just 6.01% of the total aggregate mass. However, the 2–1 mm size fraction was least abundant in irrigated land and uplands, accounting for 2.38% and 3.06% of the total aggregate mass, respectively. Orchards had the lowest content of >2 mm aggregates (3.83% of the total aggregate mass; Figure 3). Macro-aggregates were more abundant in paddy and forest soils, while silt and clay were more abundant in uplands.

Inspection of normal Q-Q plots revealed that the aggregate stability index data satisfied the criteria for normality (

Table 2. Soil aggregate stability. WSA: water stable aggregates; MWD: mean weight diameter; GMD: geometric mean diameter; D: and fractal dimension (D).

	Minimum	Maximum	Average	Standard Deviation	Skewness	Kurtosis
WSA	0.09	0.95	0.70	0.16	−1.52	2.99
MWD	0.22	1.59	0.63	0.22	1.56	3.75
GMD	0.12	1.33	0.46	0.17	1.78	6.74
D	1.15	2.78	2.22	0.20	−1.18	8.70

).

As shown in Figure 4, paddy fields had the highest WSA, MWD, and GMD values, followed by forest land, indicating that these land use types were associated with the greatest soil aggregate stability.

The proportion of macro-aggregates was also higher in irrigated land and forests, suggesting that the presence of water-stable macro-aggregates can help maintain soil structure and enhance overall aggregate stability. There were some notable differences between the values of the soil aggregate stability indices for different land use types, suggesting that land use type has an important effect on aggregate stability. However, the mean values of D for different land use types other than orchards were relatively similar (ranging from 2.06 to 2.28), indicating a lack of clear fractal characteristics.

3.2. Spatial Distribution of Soil Aggregate Stability

The parameters of the semivariance function for each index were calculated using GS + 9.0 (

Table 3. Geostatistical analysis results.

Variables	Models	Nugget (C0)	Partial Sill (C + C0)	Nugget/Partial Sill (%)	Range (km)	R2	RRS
WSA	Gauss	0.016	0.032	0.5	14,809.03	0.862	0.00005
MWD	Gauss	0.066	0.133	0.496	12,124.36	0.825	0.00105
GMD	Gauss	0.082	0.166	0.494	14,514.59	0.735	0.00234
D	Gauss	0.009	0.046	0.196	1014.98	0.347	0.00383

).

Gaussian models were generated for each soil aggregate stability indicator, giving regression coefficients (R²) ranging from 0.347 to 0.862. The models for all indices other than D had high R² values. The residual sum of squares (RSS) for each index was small, indicating acceptably reliable data fits. The nugget (C0) parameter represents the undetectable experimental error and microscale variation. The C0 values for the aggregate stability indices ranged from 0.009 to 0.082, indicating that the error between sampling points was small. The partial sill (C + C0) values reflect the spatial variability due to natural and socioeconomic factors and were similarly low for all indices, indicating that environmental factors were the dominant source of variation in soil aggregate stability. The nugget:partial sill ratio indicates the strength of the system’s spatial correlation, which is generally considered to be strong if the ratio is <25%, intermediate if it is between 25% and 75%, and weak otherwise. The ratios for WSA, MWD, and GMD were close to 50%, indicating intermediate spatial correlation, while that of D was smaller, indicating a strong spatial correlation. The range variable indicates the distance over which spatial autocorrelation exists for the index in question. All indices other than D had relatively high range values. These results indicated that Kriging interpolation analysis could be performed for all four indices. However, the large ratio of study area to the number of samples and the uneven distribution of samples reduces the accuracy of the calculation of the semivariance function, which has an impact on the interpolation results.

Ordinary Kriging interpolation analysis was therefore performed in ArcGIS 10.7 by combining the semivariance functions discussed above to characterize the spatial distributions of WSA, MWD, GMD, and D in the study area (Figure 5). In general, higher values of WSA, MWD, and GMD were associated with smaller values of D. The soil WSA in the reservoir area ranged from 0.087 to 0.950, with higher values in the south (0.697 to 0.829), followed by the north, and lower values in the west, northwest, and central regions.

The spatial distributions of WSA, MWD, and GMD were all similar and indicate that the southern and western parts of the reservoir area have relatively high proportions of large soil aggregates and good soil stability whereas the central and northwestern parts have fewer macro-aggregates and worse soil stability. The values of D vary between 1.14 and 2.78 due to the low variance of the corresponding semivariance function. Together with the low spatial correlation distance of D, this meant that the distribution obtained for this index by Kriging interpolation was relatively fragmented. Overall, D was low in the reservoir area’s south-central and north-western regions but high in the central and eastern parts, exhibiting a negative spatial correlation with WSA, MWD, and GMD.

In conclusion, the soil aggregates are most stable in the southwest, west, and northeastern regions of the studied area, followed by the northern and southeastern regions, while the central region’s soil aggregates are less stable. By comparing these findings to the land use map shown in Figure 1, it can be seen that the regions of the study area with high soil aggregate stability are dominated by paddy fields and forests, while those with poor aggregate stability are dominated by construction land and uplands. This confirms the general conclusion that soil aggregates in paddy fields and forests are more stable than those in upland sites.

3.3. Relationships between Soil Properties and Soil Aggregate Stability

The semivariance functions and Kriging interpolation results discussed in the preceding section can also be used to analyze the relationships between the spatial distribution of soil aggregate fractions and soil properties. The spatial heterogeneity of aggregate stability indices is influenced by natural environmental factors and human activities, the latter of which can affect both the internal properties of the soil and socio-economic factors such as topographic conditions and land use patterns. Topographic factors influence soil properties by affecting local regional climatic conditions, vegetation cover, water, and thermal conditions as well as land use. The environmental factors whose influence on soil aggregate stability was evaluated in this work were ELE, aspect, slope, LST, ΔLST, sand, silt, clay, pH, SOM, TP, AP, TN, AN, TK, and AK (

Table 4. Statistical data on soil properties in the studied area. ELE: elevation; ΔLST: seasonal land surface temperature difference; LST: temperature; Sand: soil sand content; Silt: soil silt content; Clay: soil clay content; SOM: soil organic matter; TP: total phosphorus; AP: available phosphorus; TN: total nitrogen; AN: available nitrogen; TK: total potassium; AK: available potassium.

	Minimum	Maximum	Average	Standard Deviation	Variance	Skewness	Kurtosis
ELE	167.00	363.00	213.11	33.77	1140.46	1.74	4.30
aspect	0.00	352.88	164.01	104.64	10,948.84	0.17	−1.09
slope	0.00	17.21	4.91	3.54	12.51	1.11	1.09
ΔLST	3.85	9.46	7.18	1.02	1.05	−0.88	1.27
LST	5.71	14.36	10.68	1.70	2.89	−0.54	0.26
Sand	19.53	71.95	40.09	11.30	127.69	0.58	0.32
Silt	4.01	46.15	27.51	7.28	52.98	−0.39	1.03
Clay	18.04	53.25	32.40	7.07	49.93	0.45	0.42
pH	4.18	7.36	6.17	0.76	0.58	−0.62	−0.17
SOM	9.14	74.35	27.75	12.70	161.35	1.66	3.20
TP	0.11	1.55	0.44	0.30	0.09	1.95	4.17
AP	4.36	93.83	24.25	17.82	317.44	1.82	3.73
TN	0.03	4.30	1.07	0.90	0.80	1.56	2.49
AN	54.60	735.00	199.64	86.40	7464.23	2.68	14.86
TK	2.00	24.30	14.22	6.09	37.08	−0.51	−0.78
AK	13.20	632.50	142.18	89.78	8061.02	2.65	10.02

). The direct and indirect effects of these factors on soil aggregate stability were investigated using Pearson correlation analysis and path analysis.

Pearson correlation analysis was used to obtain the correlation matrix for the soil aggregate stability indices and soil physicochemical properties (Figure 6). There were strong and significant positive correlations between WSA, MWD, and GMD, and all of them had significant negative correlations with D, in accordance with the Kriging interpolation results. Sand had significant negative correlations with MWD and GMD, while clay, pH, SOM, and TN had significant positive correlations with WSA, MWD, and GMD. Clay, pH, SOM, and TN also exhibited significant positive correlations with WSA, MWD, and GMD, while D was significantly and negatively correlated with pH, SOM, slope, and aspect. Finally, the temperature variables LST and ΔLST had significant negative correlations with WSA, MWD, and GMD. These results show that multiple soil physicochemical properties and environmental factors significantly affect soil aggregate stability.

Path analysis was used to determine the direct and indirect effects of soil physicochemical properties on soil aggregate stability by using stepwise linear regression method in SPSS (Figure 7). The route analysis included WSA, MWD, and GMD, but no significant influential elements were found for D, so it was impossible to perform a path analysis for this index. Indicators with significant effects on WSA were SOM, AN, pH, and clay. The indirect effect of pH was substantial: 37.7% of its total impact on WSA resulted from indirect effects on SOM and AN. AN had a negative direct effect on WSA but also had a positive indirect effect via SOM and Clay. Indicators with significant effects on MWD and GMD included SOM, Clay, ELE, and pH. The latter two of these indicators had relatively strong indirect effects; ELE had a negative indirect effect on SOM via clay and pH as well as a strong positive direct effect on MWD and GMD. Conversely, pH had a positive indirect effect on MWD and GMD via SOM. In this context, it should be noted that the residual factors of WSA, MWD, and GMD in the path analysis were 0.84, 0.72, and 0.75, respectively; the magnitude of these values indicates that the path analysis neglected some important influential factors [29].

To verify the influence of land use patterns on soil aggregate stability, we investigated the relationships between land use types with the environmental factors and land use types with soil physicochemical properties found to correlate significantly with soil aggregate stability pH, Clay, SOM, TN, ELE, and Silt (Figure 8). The observed distribution of ELE is explained by the fact that higher parts of the reservoir area were dominated by forests and orchards, while lower parts were dominated by paddy fields and irrigated land. Paddy fields and forests were the land use types with the highest soil aggregate stability, followed by orchards and irrigated lands; grasslands and uplands had the lowest aggregate stability. In accordance with the previously discussed correlations, pH, Clay, SOM, and TN were highest in paddy fields and forests also had high SOM values. The lowest values of TN, SOM, and pH occurred in uplands, and grasslands also had low SOM and TN values. This shows that land use type significantly affects the soil elements and leads to differences in the stability of soil aggregates. It is worth noting that the large error in the data may be due to the sample size setting and the measurement accuracy, which needs to be regarded in the future study (Figure 8).

Soil texture distribution is shown in Figure 9, revealing that clay is predominant in paddy fields where soil aggregates were more stable. Sandy clay loam is also one of the main soil types, and the soil type with worse aggregate stability than clay is dominated by dry land and orchards, which is consistent with our conclusion.

4. Discussion

4.1. The Effects of the Environment on Soil Aggregate Stability

The topographic slope of land is known to correlate positively with the degree of surface soil erosion and the resulting loss of organic carbon and macro-aggregates. Conversely, elevation, and aspect are positively associated with soil aggregate stability. An earlier study on the hilly areas of the Lodgen watershed in Iran confirmed the conclusion that soil structure stability is higher in areas with shallower slopes [30], while another study found that topographic factors such as slope and aspect affect aggregate stability indirectly via their effects on vegetation activity (flat terrain and sunny areas have more lush vegetation), soil microorganisms (microbial carbon content is higher in flat land), hydrothermal conditions, and land use [31]. This is broadly consistent with our finding that elevation affected MWD and GWD indirectly via pH and Clay.

Land use changes under specific climatic conditions can affect soil aggregate stability. We also found that increases in LST and ΔLST had a significant negative impact on soil aggregate stability. This may be because seasonal temperature changes affect soil bulk density, soil moisture, and soil permeability, leading to higher erodibility [32]. It has also been suggested that intense temperature change process reduce soil aggregate stability [33]. However, another study suggested that aggregates are compressed when ice crystals grow in the pore spaces upon freezing, increasing their stability. However, upon thawing these crystals become liquid pore water, which disfavors macro-aggregate formation [34]. These results indicate a need to clarify the effects of freeze–thaw cycles on soil aggregate stability.

4.2. The Effects of Soil Properties on Soil Aggregate Stability

Previous studies have shown that soil texture, organic matter, and soil oxides dominated by iron and aluminum have important effects on soil aggregate stability. Our results indicate that soil organic matter has a significant positive effect on soil aggregate stability, which may be mainly due to the enhancement of soil structural stability resulting from the colloidal effect of organic matter [35,36]. The addition of organic matter promotes the formation of more stable multistage aggregates by increasing matrix molecular polarity and strengthening attractive electrostatic forces [37]. Additionally, interactions between organic material and minerals can strengthen the Van der Waals forces between adjacent particles in soil aggregates and thus mitigate electrostatic repulsion, leading to improved soil aggregate stability [38]. This suggests that soil structure stability and the quality of farmland around the Shitoukoumen Reservoir could be improved by treatment with organic fertilizer.

Soil texture is another important factor affecting soil aggregate stability; previous studies have found that stickier soil textures are associated with an increased abundance of macro-aggregates [39]. Aluminum and clay were identified as the most important binders for soil macro-aggregates, but elemental Fe also has a notable binding effect [39]. Additionally, Angers [40] suggested that mucilage enhances soil sorption by providing a direct surface for soil contact, leading to the formation of larger aggregates. The influence of clay particles is related to the soil-forming parent material, which is dominated by dark brown loam, black soil, white pulpy soil, and meadow soil in the research area. All of these materials are rich in expansive clay minerals such as montmorillonite and kaolinite. Taken together with the findings of Brubaker et al. [41], these results clarify the effects of clay on soil aggregate stability and highlight the beneficial effects of expansive minerals.

We also found that aggregate stability correlated positively with soil pH, TN, TK, TP, and AK, with the correlations for TN and pH being statistically significant. These results suggest that soil AP and AK promote the growth of small water-stable aggregates into larger ones and thus increase aggregate stability [42,43]. The amount of nitrogen in the soil is a key measure of its fertility, and soil macro-aggregates are one of the most significant nitrogen reservoirs. Additionally, nitrogen is one of the main aggregate binders [44]. Therefore, the more macro-aggregates there are, the more nitrogen is present in the soil. The positive effect of pH on macro-aggregates may be due to its beneficial effect on the stability of plant roots and their surface residues, which promote microbial activity and thus enhance the soil’s content of organic carbon and thus favor macro-aggregate formation [45]. This hypothesis is consistent with our observation that pH had positive indirect effects on WSA, MWD, and GMD via SOM [11,45]. However, other studies suggest that high pH suppresses macro-aggregate formation, possibly by reducing the thickness of the organic matter layer on the mineral surface or by causing clay particles to acquire negative surface charges and to thus repel one-another, causing clay dispersion [30,46]. Efforts should therefore be made to identify optimal pH ranges that maximize soil aggregate stability for different soil textures.

The path analysis revealed some of the indirect relationships between soil properties and aggregate stability. Clay, pH, and SOM were found to have indirect effects on all three aggregate stability indicators included in the path analysis, with pH being the most influential. The indirect effect of clay on SOM was attributed to the retention of SOM on clay, in accordance with previous reports indicating that clay is a more important reservoir of soil carbon than sand and silt [47]. AN also has an indirect effect on WSA via SOM because soil carbon cycling and N uptake are mutually coupled processes [48], and maintaining a reasonable C:N ratio has a positive effect on soil fertility. The soil N content in the research area is low, so adding organic N fertilizer could enhance soil fertility and promote macro-aggregate formation. However, SOM had a negative indirect effect on WSA via alkaline N, which is inconsistent with the previously discussed effect; the reason for this discrepancy is unclear. In summary, the influence of soil properties on soil aggregates is a complex and integrated mechanism. It is necessary to improve soil texture, add SOM, N, and P to improve the quality of soil aggregates.

4.3. The Effects of Land Use Type on Soil Aggregate Stability

The distribution of soil aggregate size fractions across different land use types observed in this work is consistent with previous reports [49]. Land use types significantly affected soil aggregate stability, with paddy fields and forests having the highest stability, followed by orchards, irrigated land, and grasslands; uplands had the lowest aggregate stability and the most fragmented aggregates, in accordance with previous findings [24,50,51]. In general, natural undisturbed soils have higher proportions of macro-aggregates and more stable aggregates than those subjected to artificial disturbance [52]. The agricultural soils in research area are mainly located in uplands and their soil aggregates are less stable due to suboptimal usage of chemical fertilizers and pesticides under artificial tillage conditions that cause mechanical and structural damage [52,53]. In contrast, forest land is rich in natural organic matter and plant root residues, suffers less soil disturbance, and has a more well-developed microbial community. Consequently, its soils have a high content of organic matter that favors macro-aggregate formation and increases aggregate stability [54,55]. Our results confirmed the beneficial effects of the high organic matter content of forest soils. Similar benefits appear to exist in paddy fields despite their extensive artificial disturbance because they have unique water circulation systems that increase their soil organic carbon levels and thus enhance aggregate stability [56]. The fertility of arable land near reservoirs in semi-arid or arid zones could thus be maintained by establishing paddies. It has been suggested that natural grasslands have higher organic matter input than uplands and that this organic matter decomposes slowly due to limited disturbance, leading to high soil aggregate stability [7]. However, our results do not support this position.

In conclusion, soil aggregate stability, as one of the key characteristics of land quality, has been a hot spot in the field of land quality evolution and land quality sustainability. For a long time, scholars at home and abroad have made many excellent systematic studies on the effects of natural environment, anthropogenic activities, and soil essential characteristics on soil aggregates, and detailed mechanisms of individual key influencing factors have been elaborated. The presented findings in this study highlight the interacting effects of land use type, soil attributes, and environmental factors on aggregate stability to improve soil quality.

5. Conclusions

The >0.25 mm fraction was the dominant size fraction among the water-stable soil aggregates in the studied area. Additionally, land use type significantly affected soil aggregate stability, which was highest in forests and paddy fields followed by orchards, irrigated land, and grassland; upland soil had the least stable soil aggregates. Geostatistical and interpolation analyses revealed high spatial autocorrelation of soil aggregate stability indices, confirming the effects of land use. Correlation analyses revealed that soil aggregate stability had significant positive correlations with SOM, pH, Clay, and TN but correlated negatively with environmental factors such as LST and ΔLST. Several factors indirectly affected soil aggregate stability via their effects on SOM. Overall, our results show that land use type has significant effects on soil properties and thus on soil aggregate stability. However, further research is needed to clarify the mechanisms determining soil aggregate stability and to identify additional variables with significant effects on aggregate structure. The findings of this study can provide certain theoretical reference and scientific support for soil and water conservation, resistance to soil erosion and soil erosion prevention and control, and optimal allocation of land resources in Shitoukoumen reservoir area.