Understanding soil loss in Mollisol permanent gully head cuts through hydrological and hydromechanical responses

During permanent gully development, soil losses on steep slopes and in channel beds are primarily driven by the hydromechanical response and water storage within the soil mass. However, this aspect has been largely overlooked in previous studies on gully erosion in the Mollisol region of northeast China. In this study, erosion intensities during the 111 d of the rainy season and the 97 d of the snow-melting season were analyzed in relation to soil water storage, drainage capacity, and soil suction stress. This analysis was supported by monitoring soil moisture, temperature, and precipitation, as well as experimental investigations of soil hydromechanical properties. Under the same confining stress, Mollisols at the interrupted head cut of Gully no. II exhibited a more rapid increase and more effective dissipation of pore water pressure compared to those at the uninterrupted head cut of Gully no. I. The combination of the soil water characteristic curve and the hydraulic conductivity function revealed that the Mollisols in Gully no. II had a lower air-entry pressure and higher saturated hydraulic conductivity during wetting and drying cycles than those in Gully no. I. The head cut area of Gully no. II demonstrated a rapid water infiltration and drainage response coupled with high soil water storage capacity. The absolute suction stresses within the Mollisols of Gully no. II were lower than those in Gully no. I, potentially leading to greater erosion per unit of steep slope area. Notably, gravitational mass wasting on steep slopes was closely associated with soil suction stress, and a correlation was observed between erosion per unit in the gully bed area and soil water storage. Therefore, predicting soil loss in permanent gullies requires more emphasis on soil water storage and the hydromechanical response of the soil mass rather than solely on rainfall amounts. Specifically, considering the required water storage capacity to generate runoff intensity and reduce suction stress may enable more accurate predictions of soil loss at the permanent gully head cut.

1 Introduction

Gravitational mass wasting refers to the downward movement of rock, regolith, and/or soil caused by gravity along the sloping top layers of the Earth's surface (Evans, 2004; Allen et al., 2018). This process can be classified into four types based on the speed of material movement and moisture levels: falls and avalanches, landslides, flow, and creep (Bierman and Montgomery, 2014). Mass-wasting events occur in various sizes with undetermined failure planes and are influenced by both hydrological and hydromechanical responses (Stein and Latray, 2002; Rengers and Tucker, 2014). On the steep slopes of permanent gullies, gravitational mass wasting typically involves debris-free soil falling due to bed undercutting caused by intensive channelized flow or persistently high soil moisture (Harmon and Doe, 2001). Soil loss during the rainy season results from steep slopes losing support provided by debris deposits, while soil loss during the melting season may occur due to persistent low soil suction stress. In unsaturated soil mechanics, a high potential for or intensity of soil loss from gravitational mass wasting is associated with low soil suction stress (Lu and Godt, 2013). However, it remains unclear whether soil loss from gravitational mass wasting is consistently correlated with soil suction stress during these two stages.

Permanent gullies are initiated in areas where concentrated flows erode and transport bed sediments (Kirkby and Bracken, 2009; Sidle et al., 2017) and expand when gravitational mass wasting follows instantaneous or prolonged water infiltration (Poesen et al., 2010; Tebebu et al., 2010). The development of permanent gullies can be characterized by factors such as the topographical threshold and volumetric retreat rate of gully head cuts (Svoray et al., 2012; Torri and Poesen, 2014; Guan et al., 2021; Zare et al., 2022), the gully length–area–volume relationship (Li et al., 2015, 2017), and their role in upstream drainage areas during rainy days (Hayas et al., 2019). Soil loss from permanent gullies is largely governed by hydrological factors (Gómez-Gutiérrez et al., 2012), including the flow rate, total water volume, rainfall intensity and amount, and hydromechanical properties of the soil mass. These soil properties are influenced by land use, plant roots, texture, and structure. The hydrological processes near the head cut, the hydromechanical response of the soil mass to water infiltration, and their relationships with soil loss due to gravitational mass wasting remain poorly understood. Under natural conditions, water infiltration occurs following rainfall or snow-/ice-melting events. The infiltration rate is strongly influenced by the amount and intensity of precipitation, which determines soil water storage. However, the amount of stored water varies depending on the amount of rainfall, melting rate, and temperature. During the snow-/ice-melting season, prolonged soil saturation and extended periods of low soil suction stress result in longer water infiltration durations compared to rainfall events. This extended saturation may lead to increased soil loss due to gravitational mass wasting. In contrast, rain events typically generate intensive channelized flows that erode steep slopes and trigger gravitational mass wasting. Therefore, comparing soil loss between these two seasons is challenging. This issue can be addressed by considering the associated hydrological processes of head cuts and the hydromechanical responses within the soil mass.

In the Mollisol region of northeast China (MEC), over 296 000 permanent gullies have developed since 1960 (Yang et al., 2017; Dong et al., 2019). Gravitational mass-wasting processes have led to rapid gully widening due to overfarming and lack of maintenance (Wang et al., 2009). Various studies have examined hydrological processes affecting ephemeral gully development and volume disparities caused by rainfall and snowmelt (Tang et al., 2022; Jiao et al., 2023), tillage practices (Xu et al., 2018; Li et al., 2021), and morphology (Zhang et al., 2016). However, permanent gullies pose a greater threat to croplands than ephemeral gullies as soil loss from permanent gully erosion can account for 50 %–65 % of the total soil loss (Zhang et al., 2022). The relatively high area expansion ratio is influenced by the combination of permanent gullies with cropland use, large ridge orientation angles, and sunny slope orientations (Li et al., 2016; Liu et al., 2023). Tang et al. (2023) identified the rainfall threshold for permanent gully development, showing that the maximum 3 d cumulative rainfall best explained permanent gully bed erosion, while cumulative erosive rainfall was most strongly correlated with gravitational mass wasting. Gravitational mass wasting on the steep slopes of permanent gullies can occur during both the rainy season and the snow-melting season (Zhang et al., 2020; Zhou et al., 2023). Some studies have demonstrated that soil loss during the snow-melting season remarkably accounts for a large percentage of annual sediment production (Hu et al., 2007, 2009), with gully heads retreating faster during this season than in summer (Wu et al., 2008). Despite this, the hydrological processes near the gully head cut and the hydromechanical response of Mollisols to water infiltration during the two seasons have not been thoroughly documented. Additionally, the relationship between gravitational mass wasting and soil loss remains poorly understood. In the MEC, while the snow-/ice-melting season is shorter in duration than the cumulative rainy days (Wang et al., 2021a; Fan et al., 2023; Wen et al., 2024), meltwater infiltration persists for a significantly longer time than rainwater infiltration. Therefore, soil water storage may surpass drainage owing to continuous meltwater infiltration and limited water drainage pathways. In contrast, during the summer, rain infiltration temporarily increases but quickly diminishes once rainfall ceases and water drains. Stored water is primarily influenced by rainfall events and initial soil water content (Farkas et al., 2005; Xu et al., 2018). The duration of low soil suction stress, characterized by high soil moisture, differs substantially between the two seasons. Intensive rainstorms during the rainy season also generate channelized flow (Wen et al., 2021), which may erode the bed and result in gravitational mass wasting. Therefore, soil loss from gravitational mass wasting may coincide with low soil suction stress during the snow-/ice-melting season but not necessarily during the rainy season.

Soil loss from gravitational mass wasting on the steep slopes of permanent gullies remains poorly understood in the MEC. However, few studies have explored the hydrological and hydromechanical responses of the soil mass. This study investigated the effects of monitored soil water changes and suction stress on soil loss during the rainy and snow-melting seasons at the head cuts of two permanent gullies – one with no human activity and the other one experiencing human activity. Soil loss in the head cut areas during the two seasons was observed. Differences in the physical properties of Mollisols, such as pore water pressure dissipation under a given confining stress, the soil water characteristic curve (SWCC), and the hydraulic conductivity function (HCF), were compared. Soil loss per unit area on steep slopes and gully beds was analyzed in relation to soil water storage, drainage, and suction stress. The objective of this study was to characterize the relationship between soil loss intensity on steep slopes and the hydromechanical response of the soil mass, as well as the relationship between soil loss intensity in channel beds and water storage.

2 Study area

Northeast China is one of the three main Mollisol regions worldwide, covering a total area of 1 030 000 km². This region contributes 20 % of China's grain production and more than 40 % of its corn. Since the late 19th century, much of the Mollisol region has been gradually converted from native vegetation to cropland, which now constitutes 80 % of the total land area. The primary crops grown are soybean and corn. The study area lies in a typical severe gully erosion zone within the Mollisol region of northeast China, where native grasslands and forests were completely converted to croplands by 1968. This area is situated in a transitional rolling hilly region encompassing the area from the Songnen Plain to the Greater Khingan mountains in the west, of the Lesser Khingan mountains in the north, and near the Nen River (Fig. 1a). The farmland is characterized by a gently rolling landscape with a thick black organic soil layer overlying sandstone, mudstone, and sandy conglomerate.

The two permanent gullies examined in this study are located 1.4 km apart on south-facing and north-facing rolling slopes (Fig. 1b and c). The catchment area above Gully no. I is 0.22 km², with a relative relief of 25.85 m and a channel gradient of 3.3 %. In comparison, the catchment above the head cut of Gully no. II is 0.35 km², with a relative relief of 26.1 m and a channel gradient of 3.2 %. Gully no. I has a broader and deeper profile than Gully no. II (Fig. 2a and b). The mean depth of Gully no. I is 3.5 m, while that of Gully no. II is 1.23 m. The mean length and width of Gully no. I are 25.3 and 8.72 m, respectively, while those of Gully no. II are 28.2 and 5.61 m. The gully area and volume for Gully no. I are 199.3 m² and 863.6 m³, respectively. In contrast, Gully no. II has an area of 143.3 m² and a volume of 123.6 m³.

Both gullies are still expanding as they are connected to the river network that drains into the Nen River. Although grass covers the area near the sidewalls and ridges of the gullies, mass-wasting events occur frequently during the melting and rainy seasons. Differences in gully planform and depth suggest that mass-wasting processes at the sidewalls and head cuts occur at different rates and scales. The mass movement observed at the sidewalls of the two gullies differs in scale, as shown in Fig. 2c and d. Gully no. II has lower sidewall height and width than Gully no. I (Fig. 3). Notably, the head cut area of Gully no. II has been subjected to tillage activities, whereas the head cut of Gully no. I has not been subjected to these activities. Consequently, Gully no. II represents an early stage in the development of a large permanent gully.

The study area experiences a continental monsoon climate, with annual precipitation ranging from 347 to 775 mm and an average of 546 mm between 1971 and 2018 (Tang et al., 2023). Most rainfall occurs between June and August, contributing 70 %–90 % of the annual precipitation, with an average of 461 mm. Snowfall primarily occurs from November to April, accounting for 10 %–30 % of the total annual precipitation. The average temperatures in the coldest and warmest months are −22.5 and 20.8 °C, respectively, with an annual average temperature of 0 °C.

Figure 1 Location of the two permanent gullies in the Mollisol region of northeast China. (a) The red star marks the observation site in the study area (from ESRI). (b) Monitoring sites and ground controlling points at permanent Gully no. I. (c) Monitoring sites and ground controlling points at permanent Gully no. II. (The background of panel a is from ESRI. The area between the blue lines marks the gully bed and that between the mint green and blue lines marks the steep slope.

Figure 2 Close view of the steep slope and head cut of the two permanent gullies, with the (a) cross-section and upstream view of the permanent Gully no. I, (b) cross-section and downstream view of the permanent Gully no. II, (c) ground control points (blue dot circles) and the soil moisture–temperature monitoring site (yellow star) at permanent Gully no. I, and (d) ground controlling points and the soil moisture–temperature monitoring sites at permanent Gully no. II. The location of the head cut of the two gullies is shown in Fig. 1. The area between the blue lines marks the gully bed. The area between the mint green and blue lines marks the slope.

Figure 3 Difference between the two permanent gullies' cross sections. The location of the cross-section lines is shown in Fig. 1b and c.

3 Material and methods

3.1 Monitoring work

Near the gully head cut, frequency-domain reflectometry sensors were installed to monitor soil moisture and air temperature at depths of 20, 40, 60, and 80 cm (Fig. 2c). The two monitoring sites share the same rainfall records as Gully no. II (Fig. 2d). A trench was excavated to collect soil samples from these two monitoring sites. The soil samples were analyzed for pore water pressure dissipation using triaxial consolidated undrained (CU) compression tests with a GDS triaxial apparatus (GDS, UK). Unsaturated permeability was measured using the transient release and imbibition method (TRIM; Lu and Godt, 2013).

To observe the gravitational mass-wasting process during the rainy and melting seasons, the study area was scanned using numerous control points (indicated by dots in Fig. 1a and b and dashed circles in Fig. 2c and d) installed in and around the gully. An uncrewed aerial vehicle (UAV) was employed to improve the accuracy of the UAV-derived map and digital elevation models (DEMs), enabling the acquisition of highly accurate topographic data. Three UAV flights were conducted on 28 June and 17 October 2022 and 20 June 2023 following the same flight routine and image overlap settings. The first two flights in 2022 spanned 111 d during the rainy season, while the latter two covered the winter of 2022 and spring of 2023. As low soil moisture persists from October each year and snow cover in winter does not cause gravitational mass movement, the effective melting season in this study began on 15 March 2023 and lasted for 97 d. Pix4D software was used for image synthesis and gully topography generation. This software reallocates the point cloud and filters out vegetation-layer points. Since the vegetation layer, primarily composed of grass blades, varies in height, while ground points remain fixed, the vegetation layer was removed using the filtering tool. The DEM products were spatially registered in ArcGIS 10.2 using a standard orthoimage layer, ground control points, and spline functions (Table 1). The erosion depth at the head cut was determined by calculating the differences between the two DEMs. Using this erosion depth and the grid size, the linearity and erosion per unit area were calculated. Differences between the DEMs generated positive and negative terrain values, which reflected soil loss from gravitational mass wasting. The eroded soil volume per unit of steep slope surface area, referred to as erosion per unit area, was used to address the erosion caused by gravitational mass wasting.

Table 1 Detailed information on three UAV flights and the digital elevation models.

3.2 Tests of pore water pressure rising and dissipation

The consolidation module of the GDS triaxial apparatus was used to record the pore water pressure within the soil mass under a given confining stress. The soil samples were initially saturated in a vacuum pump and then consolidated in the chamber of the GDS apparatus at effective confining pressures of 100, 200, and 300 kPa with a 10 kPa backpressure. The consolidation process was completed when the pore water pressure decreased to the backpressure values.

For the pore water increasing stage, the following applies:

$$\begin{array}{}\text{(1)}& P_{\uparrow }=P_{\mathrm{0}}\times t^{b_{\uparrow }},\end{array}$$

where P_↑ is the recorded pore water pressure during the increasing stage (kPa), P₀ is the initial pore water pressure since loading (kPa), t is the time (s), and b_↑ is the rising proxy reflecting the steepness of the power-law curves of pore water pressure increase.

For the pore water dissipation stage, the following applies:

$$\begin{array}{}\text{(2)}& P_{\downarrow }={\displaystyle \frac{P_{\max }}{\mathrm{1}+b_{\downarrow }\times t}},\end{array}$$

where P_↓ is the recorded pore water pressure during the dissipation stage (kPa), P_max is the maximal pore water pressure since loading (kPa) and is the rollover point of the pore water pressure curve, t is the time (s), and b_↓ is the dissipation proxy reflecting the water drainage ability of soil mass at given confining pressure. It reflects the concavity of the pore water pressure dissipation curve.

3.3 Hydromechanical properties

TRIM was used to test the unsaturated permeability of the soil mass (Lu and Godt, 2013). The SWCC and HCF were obtained using Hydrus-1D (Wayllace and Lu, 2012). Using the models proposed by Mualem (1976) and van Genuchten (1980), the constitutive relations between the suction head (h), water content (θ), and hydraulic conductivity (K) under drying and wetting states can be represented by the following equation:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}={\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\mathit{\alpha }\left|h\right|\right)}^n}}\right]}^{\mathrm{1}-\frac{\mathrm{1}}{n}},\end{array}$$

and

$$\begin{array}{}\text{(4)}& K=K_{\mathrm{s}}{\displaystyle \frac{{\left\{\mathrm{1}-{\left(\mathit{\alpha }\left|h\right|\right)}^{n-\mathrm{1}}{\left[\mathrm{1}+{\left(\mathit{\alpha }\left|h\right|\right)}^n\right]}^{\frac{\mathrm{1}}{n}-\mathrm{1}}\right\}}^{\mathrm{2}}}{{\left[\mathrm{1}+{\left(\mathit{\alpha }\left|h\right|\right)}^n\right]}^{\frac{\mathrm{1}}{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}n}}}},\end{array}$$

where θ_r is the residual moisture content (%), θ_s is the saturated moisture content (%), α and n are empirical fitting parameters, α is the inverse of the air-entry pressure head, n is the pore size distribution parameter, and K_s is the saturated hydraulic conductivity (cm s⁻¹).

Based on the observed volumetric water content and the SWCC, the suction stress (σ^s, kPa) throughout the observation stage can be expressed as follows:

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}^{\mathrm{s}}=-{\displaystyle \frac{S_e}{\mathit{\alpha }}}{\left(S_e^{n/\left(\mathrm{1}-n\right)}-\mathrm{1}\right)}^{\mathrm{1}/n}.\end{array}$$

3.4 Soil water storage and drainage

In this study, the hydrological process of the steep slope is of utmost importance for analyzing gravitational mass wasting because of the varied soil water storage and drainage in the rainy and snow-melting seasons. Soil water is temporarily stored during rainstorms but drains after they cease. The drainage process during melting is not addressed herein because melting water constantly contributes to high soil moisture. Therefore, soil water storage (S_s) during rainstorms and the snow-melting season and drainage (S^d) after a rainstorm can be evaluated using the soil depth and the difference between the maximum soil moisture and antecedent soil moisture:

$$\begin{array}{}\text{(6)}& {\displaystyle }{S}_e={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}},\text{(7)}& {\displaystyle }{S}_{\mathrm{s}}=S_e^{\mathrm{w}}\mathrm{\Delta }{h}_i,\text{(8)}& {\displaystyle }{S}^{\mathrm{d}}=P-S_e^{\mathrm{d}}\mathrm{\Delta }h,\end{array}$$

where S_e is the degree of saturation, θ is the in situ observed volumetric moisture content being measured (%), Δh_i is the soil layer i (200 mm in this work; i=1, 2, 3, 4), $S_e^{\mathrm{w}}$ and $S_e^{\mathrm{d}}$ are the residual soil moisture in the wetting and drying processes (%), and P is the accumulated rainfall (mm) and equals 0 mm in the snow-melting season. To show the soil water storage during the rainy and snowmelt seasons and the water drainage after rainfall, all the information including rainfall amount, air temperature, soil moisture, and temperature in various soil layers was considered. The recorded rain events were categorized into four groups: light rain, moderate rain, torrential rain, and rainstorms, with rain amounts of < 10, 10–25, 25–25, and 50–100 mm, respectively.

4 Results

4.1 Erosion per unit area of gully bed and slope

The erosion per unit area in both the bed and slope areas during the snowmelt season was greater in Gully no. I than in Gully no. II (Fig. 4). This could be attributed to lower meltwater storage and higher meltwater runoff at the head cut of Gully no. I. In contrast, during the rainy season, the erosion per unit area in the bed of Gully no. II exceeded that of Gully no. I, likely due to rapid soil water storage and drainage generating intensive runoff at the head cut of Gully no. II. The primary cause of steep slope erosion in both gullies was gravitational mass wasting. For Gully no. II, the erosion per unit area during the snowmelt season was significantly higher than that during the rainy season. Additionally, during the snowmelt season, erosion per unit area on the slopes of Gully no. II exceeded that of Gully no. I. Although erosion per unit area during the rainy season was slightly higher for Gully no. I than for Gully no. II, this difference was negligible compared to the substantial variation observed during the snowmelt season. The steep slopes of the permanent gullies were primarily stabilized by soil suction stress, which is a function of the soil moisture and hydromechanical properties of the soil mass.

As channel bed erosion was closely correlated with hydrological processes and slope erosion was influenced by soil suction stress, further examination of the soil water storage, drainage, and hydromechanical properties of the soil mass in the two permanent gullies was conducted. One key difference in the hydrological processes at the head cut was that soil water storage and drainage occur during the rainy season, whereas water drainage was absent during the snowmelt season. These results could be attributed to the continuous infiltration of meltwater from snow and ice into macropores and fissures. Once the melting process was completed, soil water storage ceased, and water drainage began during the transition period between the snowmelt and rainy seasons.

Figure 4 Differences in the erosion per unit area for the gully bed and slope.

Table 2 Physical properties and pore water pressure changes in the soil mass.

4.2 Physical properties of Mollisols

4.2.1 Pore water pressure rising and dissipation

Under the same confining pressure, pronounced differences were observed in the rising and dissipation ratios of the pore water pressure within the Mollisols of the two gullies. The pore water pressure results during the consolidation process at effective confining pressures of 100, 200, and 300 kPa were compared (Fig. 5). The physical properties and the rising and dissipation ratios and proxies are listed in Table 2. The peak value of the pore water pressure within the Mollisols of Gully no. I was higher than that in Gully no. II. The peak value of the pore water pressure within the Mollisols of Gully no. II increased to 57.6, 139.0, and 141.7 kPa under the confining stresses of 100, 200, and 300 kPa, respectively. In contrast, the peak value of the pore water pressure within the Mollisols of Gully no. I increased to 87.9, 176.1, and 237.3 kPa, respectively.

Figure 5 Variation in pore water pressure under effective confining pressure of 100, 200, and 300 kPa by GDS triaxial shear tests (GDS Instruments, UK). The proxy for the pore water pressure rising and dissipation are calculated using Eqs. (1) and (2). The rising and dissipation ratio is calculated using the pore water pressure difference during a given time interval. The values of proxy and ratio are shown in Table 2.

The high peak pore water pressure indicates that the Mollisols in Gully no. II exhibited strong hydraulic conductivity, as reflected by the increased ratio. The dissipation ratio and proxy further demonstrated the connectivity of the soil pores. During the rising stage, the ratio of the Mollisols in Gully no. II was 2–18.6 times greater, and the rising proxy was 1.08–1.92 times larger than those observed in Gully no. I. In the dissipation stage, the ratios were 6.20–58.6 times greater, and the proxies were 2.65–8.0 times larger compared to the Mollisols in Gully no. I. The largest difference between the two gullies was observed under a confining stress of 200 kPa. These findings suggest that the increased pore water pressure and enhanced dissipation properties in Gully no. II are indicative of active hydrological processes at its head cut.

4.2.2 Hydromechanical properties of Mollisols

Figure 6 shows the results of the TRIM tests, SWCC, HCF, and estimated suction stress at varying degrees of saturation. Water outflow mass was measured at 10 min intervals during both the drying and wetting processes. The SWCC and HCF differed between the drying and wetting processes because water flow during the drying process relates to the applied suction level, while water flow during the wetting process was measured at a positive pressure head (Lu and Godt, 2013). The water outflow masses measured for the Mollisols in Gully no. II were generally higher than those in Gully no. I. During the drying tests, the water outflow masses for Mollisols from Gully nos. II and I were 0.0713 and 0.060 g per 10 min, respectively. In the wetting tests, the water outflow masses were 0.031 and 0.0208 g per 10 min, respectively (Fig. 6a). Overall, the permeability of the Mollisol Gully no. II was higher than that of the Mollisol Gully no. I. Similar results were obtained for pore water pressure increase, dissipation ratio, and proxy, as shown in Table 2.

Figure 6 Differences in the hydromechanical properties of the two soil masses. (a) Water flow mass in the drying and wetting process. (b) SWCC for soil mass of permanent Gully no. I. (c) SWCC for soil mass of permanent Gully no. II. (d) Suction stress–volumetric water content curves for the two soil masses. The mass of water outflow was recorded at 10 min for each test.

Using the parameters listed in Table 3, the SWCC and HCF curves of the Mollisols are plotted in Fig. 6b and c. Air-entry pressure and residual water content are two key parameters that describe the hydrological and mechanical characteristics of Mollisols. Air-entry pressure represents the critical point at which air enters saturated soil and drainage begins. The values of α^d and α^w indicated that the air-entry pressure required for Mollisols in Gully no. I was greater than that in Gully no. II, with differences of 79.4 and 28.0 kPa under drying and wetting conditions, respectively (Table 3). Therefore, water infiltration in Gully no. II, during both the rainy and snowmelt seasons, was more active compared with that in Gully no. I. Residual moisture did not vary markedly due to the similarity in soil type.

The saturated hydraulic conductivity of the Mollisols in Gully no. I was lower than that in Gully no. II during both drying and wetting processes. As shown in Table 2 and Fig. 5, the pore water pressure rising ratio and proxy, along with the dissipation ratio and proxy, further demonstrate that the permeability of the Mollisols in Gully no. II was higher than that in Gully no. I. These results suggest that the pore water pressure varied with confining stress, air-entry pressure, and saturated hydraulic conductivity under drying and wetting conditions. Consequently, it is more challenging for the Mollisols in Gully no. I to absorb and drain water compared to those in Gully no. II.

Figure 6b and c illustrates the matric suction and hydraulic conductivity at various soil moisture levels. However, direct comparisons of suction stress with other hydrological and mechanical parameters listed in Table 3 were not feasible. Hence, the suction stress at various soil moisture levels was determined (Fig. 6d). The absolute suction stress at specified soil moisture levels was higher for Mollisols in Gully no. I than for those in Gully no. II. This indicates a higher likelihood of gravitational mass wasting for the Mollisols in Gully no. II.

Table 3 Parameters describing the SWCC and the HCF from Hydrus-1D.

Note that the superscripts d and w indicate drying and wetting states.

4.3 Hydrological response

4.3.1 Monitoring results

In total, 24 light rain events, 2 moderate rain events, 5 torrential rain events, and 1 rainstorm event were recorded (Fig. 7a). During the snowmelt season, the air temperature began to rise above 0° on 20 March, with an initial gradient of 0.15 °C d⁻¹, which increased to 2.3 °C d⁻¹ after 23 April (Fig. 7b). Regarding soil moisture changes, the volumetric water content at a depth of 20 cm in Gully no. II showed a significant increase starting on 23 April, whereas only a slight increase was observed in Gully no. I. This suggests that the head cut of Gully no. II experienced higher soil moisture levels. Soil moisture patterns during the rainy and snowmelt seasons differed between the two sites. In the rainy season, the volumetric water content at a depth of 20 cm consistently remained at a lower level compared with those at the other three soil depths (Fig. 7c). In contrast, during the snowmelt season, the volumetric water content in the 40 cm soil layer was the highest (Fig. 7d). Overall, Gully no. II exhibited greater soil moisture fluctuations than Gully no. I in both seasons. This indicates that water infiltration from rainfall and snowmelt into the head cut of Gully no. II was more active than that in Gully no. I. The observed differences demonstrate that the amount of stored and drained water at the head cut of Gully no. II was significantly greater than that in Gully no. I.

To further analyze water infiltration differences during the rainy and snowmelt seasons, an in-depth comparison of the rate of soil moisture increase at a depth of 20 cm was conducted (Fig. 8). Among the four types of rain events, the mean rates of increase for Gully no. II were 0.027, 0.053, 0.102, and 0.356, respectively, which were 1.12, 1.35, 1.34, and 1.78 times higher than those for Gully no. I (Figs. 8a and 9a). During the snowmelt season, the soil moisture increase ratios at the initial, medium, and final stages for Gully no. II were 3.48, 1.60, and 1.66 times higher, respectively, than those in Gully no. I (Fig. 8b). Therefore, the water infiltration rates at the head cut areas of Gully no. II were consistently greater during both the rainy and snowmelt seasons.

Figure 7 Field-monitored rainfall conditions, air and ground temperature, and volumetric water content. (a) Rain events during the rainy season. (b) Soil, air temperature, and volumetric water content during the snow-melting season. (c, d) Monitored volumetric water content during the rainy and snow-melting seasons.

Figure 8 Volumetric water content increasing ratio in snow-melting ratio and the rainy season. (a) Rate of increase in VWC (volumetric water content) at varied rain events. (b) Rate of increase in VWC at three stages of temperature increase.

Figure 9 Hydrologic behavior for gully head cut during light rain events. (a) Lower rate of increase in VWC for Gully no. I. (b) Higher soil water storage for Gully no. II. The three crossing lines of the boxes show the 75th quantile (Q₃), median (Q₂), and 25th quantile (Q₁) from top to bottom. The length of the box is referred to as the interquartile range (IQR $=Q_{\mathrm{3}}-Q_{\mathrm{1}}$). The crossed square inside the box is the average value. The upper and lower limits of whiskers are Q₃+1.5IQR and Q₃−1.5IQR, respectively. The solid squares are the outliers.

4.3.2 Soil water storage and drainage

Figure 10 shows the stored and drained water in the soil column at the head cuts of the two gullies. During the snowmelt season, the water stored in Gully no. II was higher than that in Gully no. I. The stored water ratio was calculated by dividing the amount of water stored in Gully no. II by that in Gully no. I and was typically greater than 1.0 throughout the snowmelt season (Fig. 10a). This ratio increased sharply from 26 April, indicating that the amount of water stored in the head cuts of Gully no. II was higher.

For the four types of rain events, the mean water stored in the head cuts of Gully no. II during the 24 light rain events was greater than that in Gully no. I (Figs. 9b and 10b). The differences in stored water between the two gullies were 4.0, 8.1, 15.2, and 46.3 mm, respectively. These results show that the stored water, be it during the snowmelt or rainy seasons, was generally higher in the head cuts of Gully no. II. However, the amount of water stored in Gully no. II was not always greater. Between 26 August and 3 September 2022, the amount of water stored at the head cut of Gully no. II was lower than that in Gully no. I, which could be attributed to high temperatures and light rain events (Fig. 10c). During a torrential rainfall event on 22 September, the amount of water stored in Gully no. II exceeded that in Gully no. I. The soil water storage capacity of Gully no. II exhibited stronger fluctuations compared with that of Gully no. I. Rapid water infiltration was often followed by rapid water drainage. Figure 10d shows the water drainage and drainage ratios of the two gullies during the rainy season, where the amount of water drained from Gully no. II was higher than that from Gully no. I. This suggests that the head cut area of Gully no. II had better soil water storage capacity during both the snowmelt and rainy seasons, along with more efficient water drainage during the rainy season than Gully no. I.

In summary, rapid soil water storage and drainage in the head cuts of Gully no. II during torrential rains or rainstorms coincided with observed pore water pressure rise, dissipation, and the hydromechanical properties of Mollisols. The high permeability of Mollisols at the head cut of Gully no. II was responsible for more rapid soil water storage, drainage processes, and water retention. This could considerably influence the erosion intensity of the steep slope and gully bed in permanent gullies.

Figure 10 Hydrological response during the rainy and snow-melting season. (a) Soil water storage and the storage ratio during the snow-melting season. (b) Soil water storage at varied rain events. (c) Soil water storage and the storage ratio for the two permanent gullies. (d) Soil water drainage and the drainage ratio during the rainy season. During the rainy season, soil water storage and drainage synchronously change with the onset and end of rainfall.

4.4 Hydromechanical response and soil loss

The Mollisols in the head cut areas of the two permanent gullies differed in hydromechanical properties, leading to considerable variations in monitored soil moisture in the field. Suction stress was estimated based on field-monitored soil moisture at each site and the relationship between soil moisture and matric suction (Figs. 6d and 7c–d). During the rainy season, the absolute value of suction stress in the Mollisols of Gully no. II was lower than that of Gully no. I (Fig. 11a). Similarly, smaller absolute suction stress values were observed in Gully no. II during the snowmelt season (Fig. 11b). The lower suction stress during the snowmelt season likely contributed to strong erosion on the slopes of Gully no. II, as illustrated in Fig. 4.

As the hydrological processes in the head cut area are closely related to channel bed erosion, the hydromechanical response directly influences slope stability. It is crucial to analyze the relationships among erosion per unit area on the channel bed, soil water storage, and slope erosion with suction stress. Generally, a high absolute value of suction stress is associated with strong cohesive forces between soil particles, which enhances soil stability. Conversely, a low absolute value of suction stress indicates a higher potential for slope failure. Therefore, the relationship between the absolute value of suction stress and erosion per unit area is expected to be negative. Figure 11c shows the reciprocal relationship between the suction stress and erosion per unit area of the slope, indicating that gravitational mass wasting occurred on the slope and that the permanent gully expanded when suction stress remained relatively low for a prolonged period – approximately 5.6 kPa in this study area.

Erosion of the channel bed is closely associated with runoff discharge during erosive rain events. During such events, the amount of stored soil water decreases runoff amount and intensity. The less rainwater stored during erosive rain events, the higher the runoff amount or the more intensive the channeled flow. Consequently, the relationship between soil water storage and erosion per unit area of the channel bed is expected to be negative. Figure 11d shows the reciprocal relationship between erosion per unit area of the channel bed and soil water storage. It indicates that excessive rainwater during erosive rain events could create intensified channelized flow, eroding the channel bed once the stored water in the Mollisols reaches a threshold, which is 139.3 mm in this study area.

Figure 11 Relationship between hydrology and the hydromechanical state with the erosion per unit area over approximately 3 months. (a) Suction stress during the rainy season. (b) Suction stress during the snow-melting season. (c) Erosion per unit area on the slope decreases with suction stress. (d) The erosion per unit area on the channel bed decreases with the amount of soil water storage. The time for the monitored rainy and melting seasons was 111 and 97 d.

5 Discussion

The physical processes of permanent gully development can be categorized into gravitational mass wasting on steep slopes and sediment delivery on channel beds (Montgomery and Dietrich, 1992; van Beek et al., 2008; Luffman et al., 2015). Traditionally, most studies on gully erosion have focused on soil loss caused by water erosion and piping. Soil loss estimation is typically determined by several primary factors, including the upslope contributing area, topographic conditions, erosive rainfall, and land use (Li et al., 2015; Xu et al., 2017; Wang et al., 2021b; Tang et al., 2022). However, the physical mechanics of bed erosion and slope erosion differ, making it challenging to accurately predict soil loss on steep slopes. The gravitational mass-wasting process on a slope differs from rainfall-induced shallow landslides, particularly for those without failure planes (Poesen et al., 1998; Guo et al., 2020). Despite these differences, the two processes share similarities, such as reduced soil strength due to water infiltration (Guo et al., 2019). Therefore, a detailed mechanical analysis is necessary to understand gravitational mass wasting on slopes and sediment delivery on channel beds.

This study thoroughly investigated the effects of hydrological factors and hydromechanical properties on soil loss from both slopes and channel beds. Mass failure on hillslopes was primarily governed by suction stress, while erosion on channel beds was influenced by soil water storage and runoff amount. Therefore, hydrological factors related to soil water storage and drainage were analyzed (Fig. 10), along with volumetric changes during various rain events and snowmelt stages (Fig. 8). We also examined the hydromechanical properties and pore water pressure under a given confining stress (Table 2 and Fig. 5), relationship between the degree of saturation and suction stress (Fig. 6), and variation in suction stress during the rainy and snowmelt seasons (Fig. 11a and b). Field observations revealed two permanent gullies with distinct erosion patterns on their slopes and channel beds. Gully no. II showed signs of head cut disruption, in contrast to Gully no. I, resulting in notable disparities in erosion per unit area for both seasons and sites. The hydromechanical properties of the Mollisols differed distinctly between the two gullies, directly influencing water movement. This was evident from the observed increases in pore water pressure, dissipation ratio, and proxy. In the head cut of Gully no. II, the Mollisols were significantly disturbed, with the soil mass exhibiting higher permeability and lower suction stress at a given saturation degree. These findings indicate more active water infiltration in Gully no. II than in Gully no. I, triggered by changes in the soil's water storage and release capacity, as well as a higher ratio of volumetric water content. Consequently, the head cut area of Gully no. II experienced more intense hydrological processes. Additionally, the observed rainfall amount of 139.3 mm in this study was smaller than the 177 mm proposed by Tang et al. (2023). This discrepancy could be explained by differences in plant interception capacity and depression detention during the rainy season.

The soil water storage and drainage capacity at the head cut considerably influenced soil loss. This study primarily focused on soil water storage and its impact, and runoff was not directly addressed. From a water balance perspective, soil water storage and runoff depth were approximately equal to the rainfall depth. Consequently, the erosion per unit area of the channel bed was inversely proportional to soil water storage, as shown in Fig. 11d. Some researchers have identified factors leading to mass failures on steep slopes, including long-duration storms (Xu et al., 2020), initial soil moisture in the pre-winter season (Wen et al., 2024), tensile crack morphology (Zhou et al., 2023), and heaving and thawing (Thomas et al., 2009). The head cut of Gully no. II exhibited a high level of disturbance, resulting in greater permeability, quicker water pressure response, and higher soil moisture levels during the rainy and snowmelt seasons. Further, soil suction stress in Gully no. II was lower, leading to more intense slope erosion compared to Gully no. I. As the two gullies were only 1.4 km apart and experienced similar climatic conditions, soil properties appear to be the dominant intrinsic factor governing soil loss on gully slopes.

Gully bed erosion rates generally depend on runoff intensity. While some studies reported that runoff hydraulics during the rainy season were significantly higher than those during the snowmelt season, others have demonstrated that gully heads may retreat faster during snowmelt than in summer (Wu et al., 2008; Hu et al., 2009). In this study, the accumulated snowfall depth was high – reaching 49.6 mm – compared to the average snow depth of 30 mm. The snowfall melted intensively between 3 and 10 May 2023 (Fig. 7a and b). The heavy snowfall during the winter of 2022 and the intensive melting in early spring 2023 likely led to high soil moisture levels and intensive runoff, ultimately causing substantial bed erosion. Long-term soil saturation during the snowmelt season facilitated prolonged water infiltration and reduced suction stress. Therefore, the highest erosion per unit area occurred during the snowmelt season rather than the rainy season.

Dong et al. (2011) identified a critical soil water content for gravitational mass wasting, ranging from 31.0 % to 33.8 %, corresponding to a volumetric water content of 39.0 % to 48.0 % and a suction stress of 11.0 kPa. These findings also demonstrated that the direct shear apparatus had limitations in differentiating the contributions of effective cohesion and suction stress to total cohesion. As shown in Fig. 10b and supported by the findings of Xu et al. (2020), the high soil water storage in Gully no. II during the snowmelt season (Fig. 9a) and prolonged water infiltration lowered suction stress and increased erosion per unit area. This suggests a reciprocal relationship between absolute suction stress and erosion per unit area. The results shown in Fig. 11c and d are key findings and main contributions in the study domain of gully erosion as they clarify the role of suction stress of stored water in soil loss from steep slopes and gully beds. Additionally, our results indicate that soil water storage does not necessarily equal the rainfall amount during an event but is partially influenced by antecedent soil moisture. Figure 11 illustrates that antecedent soil moisture or precipitation substantially affects surface runoff depth and soil loss during permanent gully expansion in MEC – an aspect largely neglected in previous studies. That is, antecedent precipitation should be considered when predicting soil loss as it is closely related to soil water storage and indirectly affects runoff generation and intensity (Sachs and Sarah, 2017; Wei et al., 2007; Schoener and Stone, 2019; Wang et al., 2009). Notably, the theoretical framework underpinning this study posits that soil loss on steep slopes occurs through bank slope instability, while soil loss in gully beds results from the balance between shear forces from runoff water and soil erodibility. Therefore, soil loss in permanent gullies can be more accurately predicted using soil water storage and the hydromechanical response of the soil mass rather than relying solely on rainfall amount.

6 Conclusions

Permanent gully development is a hydrogeomorphic phenomenon, and its physical mechanics can be attributed to the hydrological and hydromechanical responses of the head cut. In the Mollisol region of northeast China, numerous studies on gully development have focused on soil loss in response to rainfall or snow depth. However, relatively few studies have addressed the physical mechanics of gravitational mass wasting. This study provides a comprehensive analysis of soil loss on steep slopes and channel beds in two permanent gullies. Our analysis considered key hydrological processes, such as infiltration, soil water storage, and drainage, as well as hydromechanical responses, including changes in suction stress levels. The following conclusions were drawn:

1. Mollisols in the head cut areas of Gully no. II exhibited a higher permeability than those in Gully no. I. This can be attributed to the elevated ratio and proxy for pore water pressure rise and dissipation. The TRIM test results confirmed that the saturated Mollisols in the Gully no. II drain faster than those in Gully no. I, owing to their higher air-entry pressure and saturated hydraulic conductivity during the wetting and drying cycles.

2. The head cut area of Gully no. II exhibited more intense hydrological processes than that of Gully no. I. This could be explained by the higher ratio of soil moisture increase observed during the four rain event types and three snow-melting stages. Soil water storage in Gully no. II experienced greater fluctuations during torrential rains and rainstorms. Overall, the absolute suction in Gully no. II remained lower than that in Gully no. I, potentially triggering greater erosion on the steep slopes.

3. The relationships between erosion per unit area on the steep slope and channel bed were analyzed for the suction stress and soil water storage. Our findings indicate that low suction stress and high soil water storage can increase gravitational mass wasting while reducing erosion on the channel bed. The two empirical relationships and their efficiency can be enhanced by incorporating data from ongoing monitoring efforts to enhance the prediction of future soil loss.