Spatial Properties of Soil Physical Quality Index S in Black Soil Croplands under Permanent Gully Erosion

: Soil physical quality (SPQ) is a limiting factor affecting crop production. However, the impact of gully erosion on the SPQ index S, deﬁned by Dexter as the inﬂection point of the soil water retention curve (SWRC), remains unclear, especially when considering different latitudinal regions. This study aimed to apply Dexter’s S-theory to evaluate the distribution of index S in black soils adjacent to various gully positions and investigate its relationship with bulk density (Bd), soil organic matter (SOM), and particle percentage. Soil properties (SWRC, Bd, SOM, and particle percentage) from nine gullies in croplands in three latitudinal regions (Harbin, Hailun, and Nenjiang in Heilongjiang province) were determined at the gully edge (GE 0 ) and 50 m beyond the edge into the croplands (GE 50 ) at the following gully units: head, mid-upper, middle, mid-lower, tail, and conjunctions between main gully and gully branch. The S-index was calculated using parameters such as n , θ s , and θ r , with SWRC data ﬁtted into the van Genuchten model. The results showed spatial variations in the S-index across latitudinal regions, with slightly higher S-values in Harbin than in Hailun and Nenjiang. The S-index also showed noticeable differences at GE 0 and GE 50 and at the junctions between the main gully and its branches. Approximately 51% of the samples at GE 0 and 28.2% of the samples at GE 50 had S-values below 0.035, which Dexter proposed as the boundary between good and poor SPQ, indicating a degradation of SPQ at the gully-surrounding areas. A decreased S-index in the gully vicinity was signiﬁcantly ( p < 0.05) associated with increased bulk density (1.33 vs. 1.21 g cm − 3 for GE 0 and GE 50 ) and decreased SOM (36.80 vs. 39.36 g kg − 1 for GE 0 and GE 50 ). In summary, this study indicates that gully erosion affects the farmland S-index at the gully-surrounding areas through SOM and Bd. Accordingly, measures suited to the increase in the S-index of the gully-surrounding areas may be implemented to maximize the crop yield of farmlands.


Introduction
Soil physical quality (SPQ) plays a crucial role in maintaining the sustainability of agro-ecosystems [1].SPQ encompasses various attributes, including field capacity, available water content, organic matter content, and structural stability, all of which collectively influence the water movement [2].However, evaluating SPQ when dealing with multiple parameters can be complex.Dexter introduced a single index, denoted as S, which represents the slope of the soil water retention curve (SWRC) at its inflection point.Determining the S-index requires utilizing parameters obtained by fitting SWRC data into the van Genuchten model or employing pedo-transfer functions [3][4][5][6].The S-index provides a comprehensive assessment of SPQ and has found widespread application in quantifying soil friability, compaction, hydraulic conductivity, and crop yield in various studies [7][8][9].
Soil organic matter (SOM) and bulk density (Bd) have been demonstrated to be key factors influencing the SWRC and S-value documented in the literature [5,10,11].Gully erosion is a well-known contributor to global cropland soil degradation [12], leading to Land 2023, 12, 1641 2 of 12 the depletion of SOM and an increase in Bd due to the loss of fertile topsoil and the transportation of organic matter-rich fine particles from the upper convex locations on both sides of the gully into the gully channel [13][14][15].Gullies in Israel have revealed their influence on SOM variations between fields, at the north/south slopes of the gully, and in the gully channel [16].Similarly, in Iran, SOM at the gully edge decreased by 18% compared to locations 50-70 m away from the edge [17].A transect study conducted parallel to the gully bank indicated lower SOM values in gully head areas compared to gully tail areas [18].Moreover, accelerated water flow through gullies resulted in thinning of topsoil and the exposure of subsoil layers, subsequently leading to soil compaction and increased bulk density [19].Therefore, it is important to investigate the changes in SOM and Bd resulting from gully development to enhance the understanding of soil water retention and the SPQ index S in the gully-surrounding areas.
Both permanent and temporary gully erosion have the potential to influence soil hydraulic properties, specifically the SWRC, hydraulic conductivity, and bulk density.However, their immediate effects might differ due to the differing soil disturbances they introduce, while their long-term impacts on SPQ remain a topic of interest.Many studies have primarily examined the effects of ephemeral gully erosion on the SPQ index S in various locations, such as Sicily in Italy, Ardabil province in Iran, and Mississippi in the USA [17,[20][21][22].Ephemeral gully erosion, characterized by small channels that can be easily filled through normal tillage practices, has been found to contribute to approximately 30% of the total soil loss in actively eroding areas in Navarra, Spain [23].In contrast, a permanent gully (with channel depth >0.5 m and resistance to normal tillage obliteration) located in Cape Bon (northeast Tunisia) showed an annual soil loss rate of 4.8 mm [24].Nevertheless, comparing SPQ across different permanent gullies becomes challenging when multiple soil indexes are employed, particularly in cases of varying soil textures and SOM content.To facilitate such comparisons across 12 FAO/USDA textural classes, the single index S has been demonstrated [3].However, using the single index S to assess SPQ in different croplands affected by permanent gully erosion raises uncertainties.
The black soil region in China is a crucial grain-producing area in northeastern China.Unfortunately, gully erosion has increased greatly in this region, expanding from approximately 68.53 to 928.99 km 2 from 1950 to 2013 in Heilongjiang province, China [12].This erosion has impaired topsoil quality and leads to crop yield loss.The impact of gully erosion on SPQ may exhibit regional variations in Heilongjiang province due to considerable differences in historical farming practices: the northern (Nenjiang), middle (Hailun), and southern (Harbin) lands were reclaimed on large scales in the 1930s, 1890s, and 1860s, respectively (Nenjiang Farm History, Hailun County Records, and Bayan County Records) [25].Therefore, it is necessary to conduct research to address the effect of these regional variations of permanent gullies on SPQ index S in order to provide region-specific recommendations for effective management practices.We hypothesize that gullies lead to a reduction in SOM content and an increase in bulk density in the gully-surrounding areas, subsequently resulting in a decrease in the S-value.The objective of this research was to investigate the spatial distribution of the S-index in different positions of permanent gullies across various latitudinal regions of black soil and to study the relationship between the S-index, SOM, clay percentage, and Bd under the influence of gully erosion.

Research Site and Research Design
The soil and water erosion zone in the northeastern black soil region of China covers 218,700 km 2 , constituting 20.11% of the entire black soil area.The typically black soil region is undergoing serious water erosion challenges, with the black soil layer experiencing an annual loss at a rate ranging from 0.1 to 0.5 cm, as reported by the 2019 China Soil and Water Conservation Bulletin.In this study, nine gullies were selected from croplands averagely inclined at 1%, 2%, and 3% across three distinct latitudinal regions in Heilongjiang province: Harbin (HB) at 46.3 • N, Hailun (HL) at 47.5 • N, and Nenjiang (NJ) at 48.8 • N (Figure 1a).The selected gullies represent a variety of dimensions, which include gully depths from 0.37 to 1.22 m, gully widths from 2 to 13.34 m, and gully lengths from 151 to 2132 m.This selection encompasses both permanent gullies, identified by gully depths exceeding 0.5 m and extending up to 25-30 m, as well as ephemeral gullies, defined by depths greater than 0.5 m according to criteria outlined in [23].The landscapes in the study regions are characterized as rolling landscapes featuring gentle slope gradients of 1 • -7 • and slope lengths of 500-4000 m.The study areas experience a continental monsoon climate with annual precipitation of 347-775 mm.The rainy season occurs from May to October, accounting for 85.6% of the annual precipitation, while snowfall mainly occurs from November to April.The soil type is clay loam, classified as Udic Argiboroll in the US Soil Taxonomy.The cropland in the region is farmed with a soybean-maize rotation, and ridge-furrow cultivation is the most common tillage practice.
Land 2023, 11, x FOR PEER REVIEW 3 of 13 an annual loss at a rate ranging from 0.1 to 0.5 cm, as reported by the 2019 China Soil and Water Conservation Bulletin.In this study, nine gullies were selected from croplands averagely inclined at 1%, 2%, and 3% across three distinct latitudinal regions in Heilongjiang province: Harbin (HB) at 46.3° N, Hailun (HL) at 47.5° N, and Nenjiang (NJ) at 48.8° N (Figure 1a).The selected gullies represent a variety of dimensions, which include gully depths from 0.37 to 1.22 m, gully widths from 2 to 13.34 m, and gully lengths from 151 to 2132 m.This selection encompasses both permanent gullies, identified by gully depths exceeding 0.5 m and extending up to 25-30 m, as well as ephemeral gullies, defined by depths greater than 0.5 m according to criteria outlined in [23].The landscapes in the study regions are characterized as rolling landscapes featuring gentle slope gradients of 1°-7° and slope lengths of 500-4000 m.The study areas experience a continental monsoon climate with annual precipitation of 347-775 mm.The rainy season occurs from May to October, accounting for 85.6% of the annual precipitation, while snowfall mainly occurs from November to April.The soil type is clay loam, classified as Udic Argiboroll in the US Soil Taxonomy.The cropland in the region is farmed with a soybean-maize rotation, and ridge-furrow cultivation is the most common tillage practice.

Soil Sampling and Measurement
In 2022, georeferenced sampling points were flagged for soil sampling across different gully units (Figure 1b).These gully units included the points of head, mid-upper, middle, mid-lower, tail, and intersection positions of the main and branch gullies or of the permanent and temporary gullies.Soil samples were collected at the gully edge (GE0) and 50 m further into the croplands (GE50) at these gully units (Figure 1b).At each georeferenced point, both undisturbed core soil samples (r = 5 cm, h = 5.1 cm) and disturbed soil samples (a composite of three samples) were collected at a depth of 0-20 cm.Control soil samples were collected in fields without gullies.The sampling methodology in this study totaled 10-16 samples per gully.

Soil Sampling and Measurement
In 2022, georeferenced sampling points were flagged for soil sampling across different gully units (Figure 1b).These gully units included the points of head, mid-upper, middle, mid-lower, tail, and intersection positions of the main and branch gullies or of the permanent and temporary gullies.Soil samples were collected at the gully edge (GE 0 ) and 50 m further into the croplands (GE 50 ) at these gully units (Figure 1b).At each georeferenced point, both undisturbed core soil samples (r = 5 cm, h = 5.1 cm) and disturbed soil samples (a composite of three samples) were collected at a depth of 0-20 cm.Control soil samples were collected in fields without gullies.The sampling methodology in this study totaled 10-16 samples per gully.

Soil Water Retention Curve and S Calculation
Soil water retention characteristics were determined at seven different matric potentials (0, −5, −10, −30, −100, −500, and −1500 kPa) using undisturbed core soil samples.The soil water retention curve (SWRC) was determined using a combination of the sandbox method and centrifuge method (GR21G, Hitachi, Japan) [26].The plants' available water content (AWC = θ 30 − θ 1500 ) was calculated as the difference between field capacity (FC) at −30 kPa and permanent wilting point (PWP) at −1500 kPa.After completing the measurements of the last matric potential, the soil cores were oven-dried to determine the bulk density (Bd).
The process of S calculation is shown in a flowchart in Figure 2. To obtain the S-value, the SWRC must be plotted as the logarithm (base e) of the water potential (cm) against the gravimetric water content (g g −1 ).In this study, the S-value was obtained following the method described by Dexter [3] (Equation (A8)), using parameters obtained from fitting the SWRC data into the van Genuchten equation (Equation ( 1)).The S-values are always negative, but the S modulus is used for presentation and discussion purposes.By employing the S-index, SPQ can be classified into distinct categories: S < 0.020 (very poor), 0.020 ≤ S < 0.035 (poor), 0.035 ≤ S < 0.05 (good), and S ≥ 0.05 (very good), based on data from soils with clay content ranging from 4% to 73% in seven countries [8,27].The S equation is defined by Equations ( 1) and ( 2) below: where θ g is the gravimetric water content and θ s and θ r are the saturated and residual water content, respectively; h (cm) represents the water potential; n is a dimensionless variable that describes the shape of the curve; and A (cm −1 ) is the reciprocal suction characteristic of the soil.The parameter m is calculated as method and centrifuge method (GR21G, Hitachi, Japan) [26].The plants' available water content (AWC = θ30 − θ1500) was calculated as the difference between field capacity (FC) at −30 kPa and permanent wilting point (PWP) at −1500 kPa.After completing the measurements of the last matric potential, the soil cores were oven-dried to determine the bulk density (Bd).The process of S calculation is shown in a flowchart in Figure 2. To obtain the S-value, the SWRC must be plotted as the logarithm (base e) of the water potential (cm) against the gravimetric water content (g g −1 ).In this study, the S-value was obtained following the method described by Dexter [3] (Equation (A8)), using parameters obtained from fitting the SWRC data into the van Genuchten equation (Equation ( 1)).The S-values are always negative, but the S modulus is used for presentation and discussion purposes.By employing the S-index, SPQ can be classified into distinct categories: S < 0.020 (very poor), 0.020 ≤ S < 0.035 (poor), 0.035 ≤ S < 0.05 (good), and S ≥ 0.05 (very good), based on data from soils with clay content ranging from 4% to 73% in seven countries [8,27].The S equation is defined by Equations ( 1) and (2) below: where θg is the gravimetric water content and θs and θr are the saturated and residual water content, respectively; h (cm) represents the water potential; n is a dimensionless variable that describes the shape of the curve; and ɑ (cm −1 ) is the reciprocal suction characteristic of the soil.The parameter m is calculated as m = 1 − 1/n.

Soil Organic Matter and Particles
The SOM content was determined by oxidation with potassium dichromate and titration with the FeSO4 method.The distribution of soil sand, silt, and clay fractions was measured using the pipetting method, with pretreatment involving the removal of SOC

Soil Organic Matter and Particles
The SOM content was determined by oxidation with potassium dichromate and titration with the FeSO 4 method.The distribution of soil sand, silt, and clay fractions was measured using the pipetting method, with pretreatment involving the removal of SOC using H 2 O 2 and dispersion using sodium hexametaphosphate (5%), followed by 16 h of shaking.

Statistical Analysis
Significant differences (p = 0.05) in S-values among gully positions (GE 0 and GE 50 ) and latitudinal regions were determined using the SPSS22 (IBM Corp., Chicago, IL, USA).The spatial distribution of the S-index values on the map was conducted using ArcGIS.Regression analysis between S-values and other soil parameters (i.e., Bd and SOM) was performed using Origin2021.

Soil Water Retention Curve Fitting Parameters
Figure 3 shows the water retention curves observed at different positions adjacent to the gully.At the gully head region, the GE 0 showed lower water retention from saturation to small negative pressure compared to the GE 50 areas.At the GE 50 areas, water retention tended to retain water to larger negative pressure, yet it exhibited a faster water loss rate once the air-entry threshold was exceeded.At the gully tail, the GE 0 areas exhibited lower water retention near saturation and wilting point compared to the GE 50 areas.This trend diverged from the patterns observed at the gully head and mid-gully areas.The fitting parameters for the soil water retention curve in the van Genuchten (VG) model, specifically θs and n, generally exhibited slightly lower values at the GE 0 sites compared to both the GE 50 and control (no gully) sites (Table 1).
Significant differences (p = 0.05) in S-values among gully positions (GE0 and GE50) and latitudinal regions were determined using the SPSS22 (IBM Corp., Chicago, IL, USA).The spatial distribution of the S-index values on the map was conducted using ArcGIS.Regression analysis between S-values and other soil parameters (i.e., Bd and SOM) was performed using Origin2021.

Soil Water Retention Curve Fitting Parameters
Figure 3 shows the water retention curves observed at different positions adjacent to the gully.At the gully head region, the GE0 showed lower water retention from saturation to small negative pressure compared to the GE50 areas.At the GE50 areas, water retention tended to retain water to larger negative pressure, yet it exhibited a faster water loss rate once the air-entry threshold was exceeded.At the gully tail, the GE0 areas exhibited lower water retention near saturation and wilting point compared to the GE50 areas.This trend diverged from the patterns observed at the gully head and mid-gully areas.The fitting parameters for the soil water retention curve in the van Genuchten (VG) model, specifically θs and n, generally exhibited slightly lower values at the GE0 sites compared to both the GE50 and control (no gully) sites (Table 1).     1 The soils were collected from 0-20 cm; 2 the soils were collected from 0-10 cm; 3 the soils were collected from 0 cm; 4 the soils were collected from 0-10 cm (R: Row; IR: inter-row) No data was found for θ s , θ r , α, n, S; 5 The soils were collected from 0-20 cm.

Gully Erosion in Our Study and Other Types of Erosion in the Literature
Land 2023, 12, 1641 6 of 12

S-Index Value, SOM, and Bd Spatial Distribution
The S-value displayed spatial heterogeneity near the gully positions and across different latitudinal regions (Figures 4 and 5).The S-value exhibited a trend of GE 0 < GE 50 < no gully, although the differences among them were not statistically significant (Figure 5a-c).For example, in areas inclined at 1% toward the tail of the gully, the S-values were 0.041, 0.065, and 0.043 for GE 0 , GE 50 , and no gully areas, respectively.Among the GE 0 areas, 49% of the samples had S-values > 0.035, whereas in the GE 50 areas, this proportion was 71.8% (Figure 5a).Furthermore, the S-value notably reached its nadir at the junctures where the main gully intersected with subsidiary gullies or at transitional zones between ephemeral and permanent gullies.For example, at positions marked by red crosses in Figure 4, the S-value was 0.018 in Hailun and 0.009 in Jiusan.The coefficient of variation (CV) substantiated the spatial fluctuations of the S-value across latitudinal regions, with CV ranging from 13.91% to 32.10% in Harbin, 32.56% to 49.13% in Hailun, and 28.12% to 64.95% in Jiusan.Overall, the S-values in this study were lower than those reported for other types of erosion in the existing literature (Table 1), indicating a degradation of the soil physical quality due to gully erosion.The study also revealed that the SOM content decreased at the GE 0 (average SOM = 36.80g kg −1 ) compared to the GE 50 areas (average SOM = 39.36 g kg −1 ) (Figure 6a).Additionally, Bd increased at the GE 0 (average of 1.33 g cm −3 ) compared to the GE 50 areas (average of 1.21 g cm −3 ), with a slight increase observed with latitude (Figure 6b).Lastly, the clay percentage also decreased at the GE 0 compared to the GE 50 areas (Figures 5c and 6c).Generally, the decline in S-values in the gully vicinity and in low latitudinal areas in Figures 4 and 5 can be attributed to an increase in bulk density and a decrease in SOM, which was confirmed by a significant correlation between S, Bd, and SOM (r = 0.57 and 0.42 for S vs. Bd and S vs. SOM, respectively) (Figure 7).The study also revealed that the SOM content decreased at the GE0 (average SOM = 36.80g kg −1 ) compared to the GE50 areas (average SOM= 39.36 g kg −1 ) (Figure 6a).Additionally, Bd increased at the GE0 (average of 1.33 g cm −3 ) compared to the GE50 areas

Influence of Gully Erosion on Spatial Distribution of S
This study investigated the spatial distribution pattern of the S-index concerning gully positions and latitude in the black soil of Heilongjiang province, China.Firstly, the findings demonstrate the varying detrimental impacts of gully erosion positions on the soil physical quality index S.A higher percentage of samples (51%) at the GE0 displayed S-values below 0.035 compared to areas further away from the permanent gully edge (GE50).This observation suggests that soils at the gully vicinity are characterized by poor physical quality according to Dexter's criteria [27].Within this subset of GE0 samples, 16% exhibited S-values below 0.02, indicating an extreme degradation of soil structure according to Dexter's criteria [3].Specifically, the convergence points where the main gully intersects with subsidiary gullies exhibited the lowest S-values, indicating the most unfavorable SPQ (Figure 3).Similarly, in Sicily (Italy), the lowest soil quality index values were found inside the ephemeral gully channel and its immediate vicinity compared to fields further away from the gully channel [21].Studies conducted in Sarcham (northwestern Iran) also observed markedly lower S-index values at the gully wall (0.04) compared to lands 50-70 m outside the gully (0.06) [17].Studies conducted in the Loess Plateau of

Influence of Gully Erosion on Spatial Distribution of S
This study investigated the spatial distribution pattern of the S-index concerning gully positions and latitude in the black soil of Heilongjiang province, China.Firstly, the findings demonstrate the varying detrimental impacts of gully erosion positions on the soil physical quality index S.A higher percentage of samples (51%) at the GE 0 displayed S-values below 0.035 compared to areas further away from the permanent gully edge (GE 50 ).This observation suggests that soils at the gully vicinity are characterized by poor physical quality according to Dexter's criteria [27].Within this subset of GE 0 samples, 16% exhibited S-values below 0.02, indicating an extreme degradation of soil structure according to Dexter's criteria [3].Specifically, the convergence points where the main gully intersects with subsidiary gullies exhibited the lowest S-values, indicating the most unfavorable SPQ (Figure 3).Similarly, in Sicily (Italy), the lowest soil quality index values were found inside the ephemeral gully channel and its immediate vicinity compared to fields further away from the gully channel [21].Studies conducted in Sarcham (northwestern Iran) also observed markedly lower S-index values at the gully wall (0.04) compared to lands 50-70 m outside the gully (0.06) [17].Studies conducted in the Loess Plateau of China revealed a decline in soil quality index from 36.6% to 10.6% in the vicinity of gullies compared to gully-free lands [30].
Secondly, the S-index exhibited discernible regional spatial variations in black soils in our study, with slightly lower S-values at higher latitudes (Nenjiang) compared to middle and lower latitudes (Hailun and Harbin).Similarly, the S index estimated via pedo-transfer function successfully predicted the spatial variability of soil physical quality with an Oxisol under a livestock farming system in Mato Grosso do Sul, Brazil [31].Furthermore, our study highlighted the distinct impact of permanent gully erosion on the S-index compared to ephemeral gully erosion and sheet erosion.The S-values after permanent gully erosion in our study were slightly lower than those reported for sheet erosion, ephemeral gully erosion, and tillage erosion in other studies (Table 1) [1,20,22].These differences likely arise from variations in the spatial redistribution of soils affected by diverse types of erosion.
The S-index has proven to be a reliable indicator of soil physical quality when assessing the effects of tillage erosion (no till vs. conventional till) [29], compaction [11], and varying degrees of sheet erosion [28].By employing the S-index, SPQ can be classified into distinct categories: S < 0.020 (very poor), 0.020 ≤ S < 0.035 (poor), 0.035 ≤ S < 0.05 (good), and S ≥ 0.05 (very good), based on data from soils with clay content ranging from 4% to 73% in seven countries [8,27].An S-index > 0.035 indicates favorable soil physical conditions for plant root growth [3,27] and is beneficial for high millet yield in Sahelian Arenosol [32].In our study, the low S-values at GE 0 areas were likely responsible for the lower average soybean yield of approximately 3351 kg ha −1 compared to 4281 kg ha −1 in GE 50 areas for the same nine gullies in our published papers [33].However, Wilson [20] determined that the overall soil quality index, based on sensitive soil indicators such as soil field capacity, saturated hydraulic conductivity, and soil penetration resistance, rather than index S alone, could predict the corn yield loss due to ephemeral gully erosion and soil-scrapping to fill the gully in Mississippi, USA.They attributed the limited performance of the S-index to the relatively small difference in soil erosion intensity within a single field site.

Influence of SOM, Clay Percentage, and Bd on Index S under Gully Erosion
This study also examined the spatial distribution pattern of Bd and SOM in relation to gully positions and latitudinal regions in black soil croplands.The results found a significant (p < 0.05) negative effect of Bd and a positive effect of SOM on the index S (Figure 6).The spatial heterogeneity of S-values in the gully-surrounding areas can be attributed to the alteration of preferential entrainment/deposition of soil particles by the presence of gullies [16,34], leading to a gradual change in soil properties, particularly in terms of Bd and SOM.
Bulk density primarily exhibited higher values at GE 0 and at higher latitudes (Nenjiang).In the vicinity of the gully, black soils with high clay percentages (over 30%) are subjected to compaction from heavy machinery as normal tillage machinery often needs to make multiple U-turns when unable to cross the gullies with depths >30-60 cm [35].This high compaction resulting from past and ongoing machinery traffic impacted the index S by reducing soil capillary porosity [1].In the black soil croplands of northeasten China, agricultural machinery engages in field operations 6-8 times each year throughout the growth phases of soybeans and maize, which encompass a range of tasks, such as land preparation, sowing, tillage, harvesting, and so on [36].Thereby, these machinery tasks increased the risk of soil compaction near gullies.Furthermore, regional differences in Bd were probably related with the different degrees of mechanized operations in regions of China and Brazil [37][38][39].Although the reclamation history in higher latitudes (Nenjiang) started in the 1930s, which was shorter than that of Hailun (1890s) and Harbin (1860s), the use of higher tractor power, faster tillage speed, and a higher degree of mechanized operations in Nenjiang increased the potential risk of soil compaction [36].This explained why higher Bd was observed in Nenjiang than in Hailun and Harbin in this study.Similarly, long-term integrated agricultural production systems in Brazil have also resulted in increased compaction and Bd and decrease in SOM [39].The degree of soil compactness showed a negative correlation with the index S using 12 FAO/USDA textural classes [11].
In compacted soils resulting from agricultural machinery, the runoff flow rates experience a substantial increase of 220% compared to those in uncompacted soils [37].This augmented runoff exacerbates the reduction of soil organic matter (SOM) to a significant degree [40].In undulating landscapes devoid of gullies, SOM tends to deposit in depressed areas due to the movement of water and tillage, often resulting in heightened yields in these locations relative to the more eroded higher shoulder sites [41].However, the presence of gullies disrupts the established patterns of soil translocation.Firstly, organic-enriched soil particles from the upper convex locations on both sides of the gully were transported and ultimately lost in the gully channel rather than being deposited within the concave sites without gullies, resulting in slightly lower SOM values at the gully edge (GE 0 ) in this study.Similarly, slightly lower SOM and fine particle percentages existed in the north/south slopes of the gully and its channel compared to the north/south fields further away from the gully channel in Israel [16].The lower SOM and fine particle percentages at GE 0 of gullies contribute to differences in soil water retention properties and S-values [3,11].Additionally, the relatively higher SOM content in Nenjiang compared to Hailun and Harbin was attributed to its shorter history of reclamation.Black soil organic carbon content showed a significant decrease from 38.57 to 18.53 g kg −1 with a gradual increase in cultivation duration from 20 to 100 years [42,43].The clay percentage in this study did not have a significant correlation with the index S.
In this study, index S revealed spatial variations due to gully erosion within black soil cropland.The slightly lower S-values near gullies were due to the lower SOM and higher Bd around these areas.Management practices specific to areas with S-values below 0.035 near gullies will be essential for future gully management efforts.However, a limitation of this study is that it does not consider the S-index at a catchment scale or temporal variations.However, prior studies have demonstrated that emergence of gullies in the three latitude regions from the south to the north (Harbin, Hailun, and Nenjiang) took place during the 1920s, 1940s, and 1950s, respectively, which is consistent with the chronological order of black soil area reclamation for farming purposes [44].Therefore, visualizing temporal changes in soil characteristics across the nine gullies can be achieved through spatial depiction.Identifying the main parameters that influence the temporal and spatial distribution of the index S and implementing management strategies may have potential to guide effective land management practices at a catchment scale.

Conclusions
This study investigated the spatial distribution of the S index across various gully positions and latitudinal regions within black soil croplands.It sought to establish the relationship between the S-index, SOM, clay percentage, and Bd in the gully-surrounding areas.The S-values exhibited a noticeable reduction at GE 0 compared to GE 50 , particularly evident at the convergence points where the main gully met its branch.Additionally, a discernible pattern emerged with lower S-values observed in regions of both higher and lower latitudes within Heilongjiang province.The S-values were consistently below the threshold value of 0.035, as proposed by Dexter, indicating compromised soil physical quality attributed to gully erosion.Moreover, a significant correlation (p < 0.05) was observed between the S-index, Bd, and SOM.The increased Bd and decreased SOM at GE 0 and in the higher latitude region (Nenjiang) were responsible for the lower S-values in those gully-surrounding areas.These findings enhance our understanding of how gullies affect the S-index of black soils in terms of soil physical quality.To improve the S-values in cropland, it is recommended to implement agricultural practices geared toward reducing Bd and increasing SOM content in the vicinity of gullies.

Figure 1 .
Figure 1.Gully locations.(a) Locations of nine gullies at different latitudes in Heilongjiang Province, China; (b) soil sampling protocol carried out adjacent to a gully in Nenjiang County-Jianshan.The gully photo was taken by an unmanned aerial vehicle.

Figure 1 .
Figure 1.Gully locations.(a) Locations of nine gullies at different latitudes in Heilongjiang Province, China; (b) soil sampling protocol carried out adjacent to a gully in Nenjiang County-Jianshan.The gully photo was taken by an unmanned aerial vehicle.

Figure 2 .
Figure 2. A flowchart showing the process of soil collection and soil physical quality index S calculation.

Figure 2 .
Figure 2. A flowchart showing the process of soil collection and soil physical quality index S calculation.

Figure 3 .Table 1 .
Figure 3. Soil water retention curves at six positions (GE0 and GE50) of the head, middle, and tail of the gully in Hailun-Dacheng.Table 1.The parameters of the van Genuchten model for different positions of gully in this study (i.e., Hailun-Dacheng) and in different types of erosion in the literature.

Figure 3 .
Figure 3. Soil water retention curves at six positions (GE 0 and GE 50 ) of the head, middle, and tail of the gully in Hailun-Dacheng.

and 2023 ,
11, x FOR PEER REVIEW 7

Figure 4 .
Figure 4. Spatial distribution of the S-index values adjacent to gullies, specifically in Ha Dacheng, Nenjiang-Jianshan, and Nenjiang-Nangangtun.The S-index values are represen bars, contour lines are depicted in red, and the red crosses highlight the positions of conjun between main gully and gully branches or the transition from ephemeral gullies to permane lies in the field.

Figure 4 .
Figure 4. Spatial distribution of the S-index values adjacent to gullies, specifically in Hailun-Dacheng, Nenjiang-Jianshan, and Nenjiang-Nangangtun.The S-index values are represented by bars, contour lines are depicted in red, and the red crosses highlight the positions of conjunctions between main gully and gully branches or the transition from ephemeral gullies to permanent gullies in the field.

Figure 4 .
Figure 4. Spatial distribution of the S-index values adjacent to gullies, specifically in Hailun-Dacheng, Nenjiang-Jianshan, and Nenjiang-Nangangtun.The S-index values are represented by bars, contour lines are depicted in red, and the red crosses highlight the positions of conjunctions between main gully and gully branches or the transition from ephemeral gullies to permanent gullies in the field.

Figure 5 .
Figure 5.A comparison of the S-index values among gully positions of GE0, GE50, and the no-gully field and among latitudinal regions.(a) distribution of S-index values with sequence of sampling; (b) comparison of S-index values among the gully units at the same land slope.Different lowercase letters indicate significant differences of S-index values among the gully units (head, mid, and tail) at the same land slope; (c) comparison of S-index values among the three latitudinal regions.

Figure 5 .
Figure 5.A comparison of the S-index values among gully positions of GE 0 , GE 50 , and the no-gully field and among latitudinal regions.(a) distribution of S-index values with sequence of sampling; (b) comparison of S-index values among the gully units at the same land slope.Different lowercase letters indicate significant differences of S-index values among the gully units (head, mid, and tail) at the same land slope; (c) comparison of S-index values among the three latitudinal regions.

Figure 6 .
Figure 6.Distribution of basic soil properties in the latitudinal regions.(a) properties of soil organic matter content in three latitudinal regions; (b)properties of soil bulk density in three latitudinal regions; (c)properties of soil particle percentage at different gully positions and latitudinal regions (i.e., land inclined at 1%).

Figure 6 .
Figure 6.Distribution of basic soil properties in the latitudinal regions.(a) properties of soil organic matter content in three latitudinal regions; (b) properties of soil bulk density in three latitudinal regions; (c) properties of soil particle percentage at different gully positions and latitudinal regions (i.e., land inclined at 1%).

Figure 7 .
Figure 7. Correlation and regression relationship between index S and soil parameters.(a) Correlation between S, sand, silt, clay, SOM, and Bd; (b) regression relationship between index S and soil organic matter and bulk density.S: soil physical index S; sand: sand percentage; silt: silt percentage; clay: clay percentage; SOM: soil organic matter; Bd: bulk density.

Figure 7 .
Figure 7. Correlation and regression relationship between index S and soil parameters.(a) Correlation between S, sand, silt, clay, SOM, and Bd; (b) regression relationship between index S and soil organic matter and bulk density.S: soil physical index S; sand: sand percentage; silt: silt percentage; clay: clay percentage; SOM: soil organic matter; Bd: bulk density.

Table 1 .
The parameters of the van Genuchten model for different positions of gully in this study (i.e., Hailun-Dacheng) and in different types of erosion in the literature.Gully Erosion