Modeling of driving factors and headcut rates of ephemeral gullies in the loess plateau of China using high-resolution remote sensing images

ABSTRACT

Ephemeral gully headcut erosion contributes significantly to global land degradation and increased sediment yields, but the underlying driving factors and prediction models remain poorly understood. We conduct a comprehensive quantitative analysis of ephemeral gully headcut erosion in the Loess Plateau using an optimal parameters-based geographical detector (OPGD) model, leveraging high-resolution remote sensing images. Our findings reveal a varied ephemeral gully head advance rate spanning 0.04–5.54 m yr⁻¹ between 2009 and 2021 (average 1.37 m yr⁻¹), with over 58% of the erosion rates falling between 0.50 and 2.00 m yr⁻¹. Catchment area emerges as the primary driving factor, with an explanatory power of 61%. Moreover, the interactions between catchment area and slope degree, rainfall erosivity, and fractional vegetation coverage (FVC) had explanatory powers exceeding 80%. Furthermore, we developed a robust prediction model for ephemeral gully head advance rates based on the results from the OPGD model, incorporating the FVC factor. The validation of our model yielded a high coefficient of determination (R² = 0.92 m yr⁻¹) and low root mean square error (RMSE = 0.31 m yr⁻¹). Our study offers new insights into ephemeral gully headcut erosion control in the Loess Plateau and serves as a valuable reference for loess regions worldwide.

KEYWORDS:

• High-resolution remote sensing images
• headcut erosion
• optimal parameters-based geographical detector model
• prediction model

1. Introduction

Soil is a vital resource underpinning global food production, human habitation, and crucial ecosystem functions. However, soil erosion directly threatens land sustainability, exacerbating flood risks and increasing the frequency of natural disasters, jeopardizing global food security and the achievement of Sustainable Development Goals (Amundson et al. 2015; Bhardwaj, Wuest, and Wuest 2017; Borrelli et al. 2017; FAO 2019; Montanarella et al. 2016; Pal et al. 2021). Ephemeral gully erosion, a distinct form soil erosion prevalent on sloping farmlands in the Loess Plateau of China, poses a particular challenge. It renders land uncultivable and readily transitions into permanent gullies, severely impacting tillage, depleting fertile topsoil, and impeding the region’s sustainable socioeconomic development.

With the advancements in remote sensing technology, high-resolution satellite images have emerged as a valuable tool in studying ephemeral gully erosion. Research in this domain has primarily focused on identifying and delineating ephemeral gullies (Dai et al. 2020; Liu et al. 2022; Yang et al. 2017), mapping their susceptibility (Fiorucci et al. 2015; Mohebzadeh et al. 2024), assessing of soil losses (Karydas and Panagos 2020; Nachtergaele and Poesen 1999; Reece et al. 2023; Teasdale and Barber 2008), investigating influencing factors (Tang et al. 2022; Vandaele et al. 1997), and understanding their morphological characteristics and topographical thresholds (Gudino-Elizondo et al. 2018; Maugnard et al. 2014). The ephemeral gully head is particularly important, as it is at the dynamic forefront of gully formation, where understanding the primary driving factors and predicting advancement rates are critical for implementing effective erosion control measures.

Extensive studies have identified natural and anthropogenic factors driving ephemeral gully headcut erosion, including soil, climate, terrain, vegetation, and land use. In sub-humid tropical areas, the dominant cause of ephemeral gully head advancement and expansion is high rainfall, increasing subsurface water levels above the ephemeral gully bottom and reducing soil strength and erosion durability (Pal et al. 2022; Rengers and Tucker 2015; Roy et al. 2022; Tebebu et al. 2010; Zegeye et al. 2016). In arid and semi-arid areas, the main causes of ephemeral gully head advancement are higher temperatures and infiltration-excess runoff under extreme rainfall, resulting in summer cracks in clayey soils next to the ephemeral gully heads (Roshan and Negahban 2015; Xu et al. 2019). Vegetation can reduce the power of runoff erosion by increasing the soil infiltration capacity and runoff resistance and reducing the runoff flow velocity. Moreover, vegetative root growth can significantly enhance the anti-scourability of soil (Feng et al. 2022; Guo et al. 2019; Mukai 2017; Sun et al. 2016). Topography (e.g. profile curvature, catchment area, slope length, slope aspect, slope degree, and altitude) also controls the initiation and development of ephemeral gully heads (Arabameri et al. 2021; Chaplot 2013; Chowdhuri et al. 2021; Conoscenti et al. 2013; Jiang et al. 2021; Kariminejad et al. 2019; Saha et al. 2022; Torri et al. 2018; Vanmaercke et al. 2016; Wang et al. 2022).

The complexity and heterogeneity of environmental factors influencing ephemeral gully headcut erosion pose challenges for traditional statistical models, which often struggle to accurately capture non-linear and interactive effects. Geographical detector models offer a promising avenue to address this challenge by detecting spatial heterogeneity and uncovering underlying driving forces (Wang et al. 2010), an approach that has been applied widely in the natural sciences (Xie et al. 2021; Zhang et al. 2022; Zhou et al. 2021; Zhu, Meng, and Zhu 2020). However, traditional geographical detector models often rely on professional experience and previous research when discretizing spatial data, ignoring the possible impacts of data discretization methods and the break number of spatial strata on detection results. Song et al. (2020) developed an optimal parameter-based geographical detector (OPGD) model to enhance the reliability of geographic detection results, allowing for more effective extraction of spatial attributes and information. Despite these advancements, existing prediction models for ephemeral gully head advance rates mainly rely on empirical approaches, using the Thompson (Thompson 1964), SCS Ⅰ, and SCS Ⅱ (United States Department of Agriculture 1966) models. The Thompson model considers four factors – catchment area, slope degree, rainfall, and clay content – while the SCS I and SCS II models are based on catchment area and rainfall, with the SCS II model adding the $R_p$ factor (past average annual rate of gully head advancement). However, these models lack specificity for unique topography and climatic conditions, necessitating the development of an ephemeral gully head advance rate prediction model incorporating the relationships between ephemeral gully headcut erosion and terrain, soil, climate, and vegetation variables.

The Loess Plateau of China is the world’s largest and deepest loess deposit, covering an area of 640,000 km². Its unique characteristics, including high porosity, loose texture, strong water permeability, and concentrated rainfall patterns, have rendered it one of the most severely eroded regions globally (Fu et al. 2017; Maher 2016; Tang et al. 2015; Wu et al. 2018). Within this expansive plateau, the hilly and gully region accounts for 22.22% of its total area and represents the most extensive geomorphic area (Zhao et al. 2013). Persistent challenges of ephemeral gully headcut erosion have plagued this region for a long time, removing valuable nutrients from surface soils and forming numerous gullies and fragmented terrain that have significantly impeded economic and ecological development of the Loess Plateau. In response to these challenges, the Chinese government has implemented various measures, including comprehensive watershed management initiatives and the ‘Grain for Green’ project, which have effectively mitigated soil erosion and yielded significant economic and ecological benefits across the Loess Plateau. Despite these efforts, critical aspects such as the spatial distribution of ephemeral gully headcut erosion and the interactions among its influencing factors remain poorly understood.

This study leverages high-resolution remote sensing imagery from 2009 and 2021 to evaluate the impact of various independent variables (slope degree, catchment area, profile curvature, slope aspect, vegetation, and rainfall erosivity). By identifying key driving factors, we seek to develop a robust prediction model for ephemeral gully head advance rates in the hilly and gully regions of the Loess Plateau. The study aimed to (1) map the spatial distribution characteristics of ephemeral gully headcut erosion, (2) explore the contributions and relationships of driving factors to the ephemeral gully head advance rate, and (3) develop a prediction model for the ephemeral gully head advance rate in the hilly and gully regions of the Loess Plateau. Ultimately, this research aims to enhance the understanding and management of ephemeral gully headcut erosion, offering insights into effective erosion control strategies in the region.

2. Material and methods

2.1. Study area

The Zhoutungou watershed is located in Ansai County, Shaanxi Province, China (36°42′10″ – 36°47′10″ N, 109°09′00″ – 109°13′45″ E, altitude 1,026–1,439 m, area 33.6 km²) (Figure 1(a)). The study area is a typical small watershed in the hilly and gully region of the Loess Plateau with severe ephemeral gully erosion. There are over 2,300 ephemeral gullies in this area, with an average length of 45.84 m and width of 1.01 m, contributing to a total volume of ephemeral gully erosion reaching nearly 2.00 × 10⁴ m³ in 2021 (Liu et al. 2024). The predominant soil type in the study area is loessial soil, characterized by its loose structure and vertical joints, facilitating rapid water infiltration. Consequently, loessial soil is highly susceptible to erosion, with its weak erosion resistance making it prone to gully formation. The watershed experiences an average annual temperature of 10°C and an average annual precipitation of 550 mm. Notably, 75% of the precipitation occurs between July and September, with intense rainfall exacerbating soil erosion (Figure 2).

Figure 1. Location of the study area. (a) Geographical location, (b) remote sensing image, and (c) ephemeral gully samples in the Zhoutungou watershed on the Loess Plateau.

Figure 2. Monthly cumulative precipitation (orange bars) and average temperature (red curves) from 2009 to 2021 in the Zhoutungou watershed.

2.2. Remote sensing data collection

The high-resolution remote sensing data for 2009 was acquired from the QuickBird satellite, operating at an orbital altitude of 450 km. This data was processed into an RGB three-band Digital Orthophoto Map (DOM) with a spatial resolution of 0.5 m after data fusion. The high-resolution remote sensing data for 2021 was obtained from unmanned aerial vehicle (UAV), flying at an altitude of 220 m, with a spatial resolution of 0.5 m.

Before drone flight, image control points were marked in the basin using a STONEX S3+ RTK. These control points correct the imagery during unmanned aerial vehicle (UAV) image processing. Specifically, 18 image control points were designated within the Zhoutungou watershed. The selection criteria for these points encompassed a circular perimeter around the region, with additional points distributed as evenly as possible within the basin. The UAV launch site was strategically chosen in an open area with elevated altitude to maximize coverage, and minimize the number of UAV takeoff sorties. A total of 11,720 UAV images were captured and subsequently processed using Agisoft Metashape Professional software (Windows version 1.7.3) to generate a digital elevation model (DEM) with 0.15 m resolution and a DOM with 0.5 m resolution.

2.3. Research methodology

The research methodology, depicted in Figure 3, commenced with the spatial distribution analysis of ephemeral gully head advance rates by overlaying the distribution of ephemeral gullies from 2009 to 2021 using a DEM. Subsequently, multicollinearity analysis was conducted to assess the interrelationships among the driving factors. Following this, the OPGD model was used to identify the key driving factors influencing ephemeral gully head advance rates. Lastly, a prediction model for ephemeral gully head advance rates was developed and validated.

Figure 3. Flowchart of the methodology used in this study.

2.3.1. Extracting ephemeral gully head location

Spatial locations of ephemeral gullies in 2009 and 2021, as extracted by Liu et al. (2022), were used in this study. Ephemeral gullies exhibiting no spatial changes from 2009 to 2021 were identified based on high-resolution remote sensing images from QuickBird satellite and UAV, both with 0.5 m resolution. The process involved combining the layer of ephemeral gullies without spatial changes with the DEM layer in ArcGIS software to pinpoint the positions of ephemeral gully heads. Field investigations were conducted to verify the accuracy of these positions, ultimately determining the final positions of ephemeral gully heads. Subsequently, 92 ephemeral gullies uniformly distributed across the watershed and existing from 2009 to 2012 were selected to calculate ephemeral gully head advance rates.

2.3.2. Acquisition of driving factors

The driving factors considered in this study include terrain attributes, fractional vegetation coverage (FVC), and rainfall erosivity. Four terrain attributes, namely slope degree (Sl, °), catchment area (Ca, m²), profile curvature (Curv, m m⁻²), and slope aspect (Sa), were extracted from the DEM data obtained from the Shaanxi Geomatics Center of Ministry of Natural Resources of China (http://sxgis.cn/), with a 5 m resolution. These terrain factors were estimated using the terrain analysis extension in ArcGIS 10.2 (ESRI, USA). Sl refers to the slope or steepness above the ephemeral gully head, calculated using spatial analysis tools in ArcGIS software and defined by a pixel and eight adjacent pixels as the steepest slope in a plane. Ca refers to direct or indirect cumulative flow into the grid at the ephemeral gully head. Curv refers to the curvature characteristics of terrain in the direction perpendicular to the horizontal plane, showing the slope gradient with negative and positive values. Sa values in this region range from 0 to 360 degrees, including nine categories: flat, north, north east, east, southeast, south, southwest, west and northwest.

FVC – a quantitative indicator of surface conditions covered by plant communities – was calculated through dimidiate pixel model inversion on Google Earth Engine, as follows (Gutman and Ignatov 1998): (1) $\mathrm{FVC}={\displaystyle \frac{\mathrm{NDVI}-\mathrm{NDV}{I}_{\mathrm{soil}}}{\mathrm{NDV}{I}_{\mathrm{veg}}-\mathrm{NDV}{I}_{\mathrm{soil}}}}$(1) where $\mathrm{NDV}{I}_{\mathrm{soil}}$ is the $\mathrm{NDV}{I}_{}$ value of regions completely covered by bare soil, and $\mathrm{NDV}{I}_{\mathrm{veg}}$ is the $\mathrm{NDVI}$ value of regions completely covered by vegetation. After repeated comparisons, $\mathrm{NDVI}$ values with a cumulative frequency of 5% and 95% were taken as $\mathrm{NDV}{I}_{\mathrm{soil}}$ and $\mathrm{NDV}{I}_{\mathrm{veg}}$, respectively.

Rainfall erosivity was calculated and interpolated from precipitation data (National Meteorological Center of China, http://data.cma.cn/) in the study area and surrounding meteorological stations using the Kriging interpolation method in ArcGIS 10.2 (ESRI, USA). Rainfall erosion was calculated as follows (Zhang and Fu 2003): (2) $R_j=0.0534P_j^{1.6548}$(2) where R_j is the rainfall erosivity in the j_th year, and P_j is the rainfall amount in the j_th year. All data were resampled to a 5 × 5 m cell size in ArcGIS 10.2 (ESRI, USA) to ensure uniform spatial resolution.

2.3.3. Multi-collinearity analysis

Before using the OPGD model, multicollinearity among driving factors was assessed to identify and remove highly correlated variables to prevent misinterpretation of the model results (Garosi et al. 2019; Pourghasemi et al. 2017). b Tolerance and variance inflation factor (VIF) were calculated in SPSS V. 25 (IBM, USA), with a tolerance value <0.20 and a VIF value ≥5, indicating multicollinearity (Garosi et al. 2019; O’Brien 2007).

2.3.4. Optimal parameters-based geographical detector (OPGD) model

The OPGD model was constructed using Excel 2007 and comprises five components: parameter optimization, factor detector, interaction detector, risk detector, and ecological detector (Song et al. 2020). Before using the OPGD model, continuous variables were discretized using optimal parameters (Figure 4). The optimal parameter combinations for discretization methods and break numbers varied for different explanatory variables, with a quantile break with ten intervals for catchment area and FVC, an SD break with six intervals for profile curvature, an equal break with ten intervals for rainfall erosivity, and a quantile break with nine intervals for slope degree.

Figure 4. Spatial discretization parameter optimization for each continuous variable.

Factor detection involves assessing the spatial differentiation of the dependent variable (Y) and the explanatory power of the interpretation variable (X) using the q-statistic: (3) $\mathrm{Q}=1-{\displaystyle \frac{{\sum }_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}}$(3) where h is the stratification (classification or partition) of the dependent variable Y or the independent variable X, N_h and N are the units for layer h and region, respectively, ${\sigma }_h^2$ and ${\sigma }^2$ are the variances of layer h and region, respectively. Q values can be transformed to a non-central F distribution and tested for significance using Eqs (4) and (5): (4) $F={\displaystyle \frac{N-L}{L-1}{\displaystyle \frac{Q}{1-Q}\sim F(L-1,N-L\lambda )}}$(4) (5) $\lambda ={\displaystyle \frac{1}{{\sigma }^2}\left[\sum _{h=1}^L\overline{Y_h^2}-{\displaystyle \frac{1}{N}{\left(\sum _{h=1}^L\sqrt{N_h}{Y}_h\right)}^2}\right]}$(5) where λ is the non-central parameter, and Y_h is the mean value of layer h.

The interaction detector compares the Q values (Q (X₁) and Q (X₂)) to the interaction q value (Q (X₁ ∩ X₂)) to identify their interactions. The relationship between the variables can be categorized as: non-linear weakness (Q (X₁ ∩ X₂) < Min (Q (X₁), Q (X₂))), single-factor non-linear weakness (Min (Q (X₁), Q (X₂)) < Q (X₁ ∩ X₂) < Max (Q (X₁), Q (X₂))), two-factor enhancement (Q (X₁ ∩ X₂) > Max (Q (X₁), Q (X2))), independent (Q (X₁ ∩ X₂) = Q (X₁) + Q (X₂)), or non-linear enhancement (Q (X₁ ∩ X₂) > Q (X₁) + Q (X₂)).

The risk detector evaluates whether the spatial patterns represented by mean values differ significantly among sub-regions classified by a categorical or stratified variable using the t-test, whereas the ecological detector assesses whether one explanatory variable has a greater impact than another using the F-statistic. We used the Kolmogorov–Smirnov test in SPSS software to test the normality of each variable before performing the t-test and F-statistic, with the results showing that each variable followed a normal distribution.

2.3.5. Construction and verification of the prediction model for ephemeral gully head advance rate

For constructing the prediction model, 70% of the data were selected randomly, while the remaining 30% were used for verification. Using multiple regression and correlation analysis methods, the OPGD model was used to identify the main driving factors affecting the ephemeral gully head advance rate by constructing an empirical formula of the ephemeral gully head advance rate with the main driving factors. The validity, goodness of fit, and efficiency of the model were validated using the coefficient of determination (R²) and root mean square error (RMSE), calculated as follows: (6) $\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{(P_i-M_i)}^2}}$(6) (7) $R^2={\left\{{\displaystyle \frac{\left[{\sum }_{i=1}^n(P_i-\overline{P_i})(M_i-\overline{M_i})\right]}{\sqrt{{\sum }_{i=1}^n(P_i-\overline{P_i})}\sqrt{{\sum }_{i=1}^n(M_i-\overline{M_i})}}}\right\}}^2$(7) where $P_i$, $M_i$, $\overline{P_i}$, and $\overline{M_i}$ indicate the predicted, measured, mean predicted, and mean estimated ephemeral gully head advancement rates (m yr⁻¹), respectively.

3. Results

3.1. Spatial distribution of ephemeral gully head advance rate

Ephemeral gully headcut erosion within the study area was concentrated primarily in the southern and western regions of the watershed (Figure 5(a)). These areas are prone to soil erosion due to steep slopes and significant human disturbance. Advance rates were calculated for 92 ephemeral gully heads (Table S1). Notably, two typical ephemeral gullies, ephemeral gully 68 and ephemeral gully 85, experienced advancements of 30.93 m (advance rate 2.58 m yr⁻¹) and 40.50 m (advance rate 3.38 m yr⁻¹), respectively, from 2009 to 2021 (Figure 6). Overall, ephemeral gully head advance rates ranged from 0.04–5.54 m yr⁻¹ from 2009 to 2021 (average 1.37 m yr⁻¹). Figure 5(b) illustrates that rates between 0.50–1.00 m yr⁻¹ and 1.00–2.00 m yr⁻¹ accounted for 30.43% and 28.26% of the total, respectively, while rates >3.00 m yr⁻¹ accounted for less than 10%.

Figure 5. (a) Spatial distribution and (b) the relative frequency of ephemeral gully head advance rates from 2009 to 2021.

Figure 6. Ephemeral gully head advancement from 2009–2021 at two ephemeral gully sites: (a) ephemeral gully 68 and (b) ephemeral gully 85.

3.2. Multi-collinearity analysis of driving factors

A multicollinearity analysis used the six driving factors corresponding to the 92 ephemeral gully head advance rates (Table S1). Table 1 presents the VIF and tolerance values for each factor. Notably, the highest VIF (1.25) and lowest tolerance (0.80) values were associated with FVC. Generally, a tolerance value <0.20 and a VIF value ≥5 indicate a multicollinearity issue (Garosi et al. 2019; O’Brien 2007). However, none of the factors exhibited values suggesting multicollinearity. Consequently, all driving factors were included in the modeling process.

Table 1. Multi-collinearity testing indices for the independent variables in the study area.

Driving factors	Collinearity statistics
	Tolerance	VIF
Slope degree	0.82	1.22
Slope aspect	0.87	1.15
Profile curvature	0.85	1.17
Catchment area	0.84	1.19
FVC	0.80	1.25
Rainfall erosivity	0.89	1.13

Table 2. Ephemeral gully head advance rates in other global sites.

Continent	Region	Ephemeral gully head advance rate (m yr^–1)	References
Africa	Morocco	0–0.31	Marzolff and Ries (2007)
	Ethiopia	0.53–18.40	Yibeltal et al. (2019), Frankl et al. (2012), Nyssen et al. (2002)
	South Africa	0.33–1.67	Grellier et al. (2012)
	Kenya	0.8–1.5	Oostwoud Wijdenes and Bryan (2001)
	Tunisia	0.29–0.46	Maaoui et al. (2012)
Europe	Spain	0.07–0.8	Marzolff and Ries (2007), Martı́nez-Casasnovas (2003)
	Romania	0.92	Ionita (2006)
	Poland	0.21–0.52	Malik (2008)
Asia	Iran	0.1–2.57	Nazari Samani et al. (2010), Nazari Samani et al. (2018)
	Loess Plateau in China	0.23–2.39	Jiang et al. (2021), Xu et al. (2019)

3.3. Contributions of driving factors to ephemeral gully head advance rate

Using the OPGD model, we analyzed the contributions of driving factors to the ephemeral gully head advance rate. According to the factor detector results (Figure 7(a)), catchment area exhibited the greatest explanatory power (61.00%), followed by FVC (37.43%), rainfall erosivity (24.50%), slope aspect (21.26%), profile curvature (18.12%), and slope degree (10.06%). Interactions between different driving factors significantly enhanced the explanatory power of ephemeral gully head advance rates (Figure 7(b)). Bi-variable mutual enhancement was observed between catchment area and slope aspect, catchment area and FVC, and catchment area and rainfall erosivity, with all other interactions non-linear. Notably, the interaction between catchment area and slope degree emerged as the dominant factor, with an explanatory power of 81.64%, followed by catchment area and FVC (81.12%), and catchment area and rainfall erosivity (80.57%). Steeper slopes are more easily to generate surface runoff, and larger catchment areas accumulate more runoff. Therefore, the interaction between slope degree and catchment area has a greater impact on ephemeral gully erosion than catchment area with FVC and rainfall erosivity.

Figure 7. OPGD-based explanatory variable exploration of ephemeral gully head advance rate: (a) factor detector, (b) interaction detector, (c) ecological detector, and (d) risk detector.

The ecological detector revealed significant differences in the spatial distribution of the ephemeral gully head advance rate between the catchment area and other driving factors and between FVC and other driving factors, except for slope aspect, using the F-statistic at a 0.05 significance level (Figure 7(c)). Furthermore, the risk detector results identified high-risk areas of ephemeral gully head advance rate based on t-test results at a 0.05 significance level (Figure 7(d)). These high-risk areas were characterized by catchment areas >375 m², FVC ranging from 0 to 34%, profile curvature from –17 to –9 m m^–2, cumulative rainfall erosivity from 1,905 to 1,908 MJ mm hm^–2 h^–1 yr^–1, slope aspect from 112.5 to 127.5°, and slope degree from 26 °to 27.5°.

3.4. Relationship between main driving factors and ephemeral gully head advance rate

The interactions between catchment area and slope degree, rainfall erosivity, and FVC explained more than 80% of the ephemeral gully head advance rate in the OPGD model. Thus, we conducted regression analysis to explore these relationships (Figure 8). The results revealed a strong power function relationship between ephemeral gully head advance rate and catchment area, slope degree, FVC, and rainfall erosivity. Specifically, the ephemeral gully head advance rate positively correlated with catchment area, slope degree, and rainfall erosivity but negatively correlated with FVC.

Figure 8. Regression function fitting curves for ephemeral gully head advance rate and (a) catchment area, (b) FVC, (c) rainfall erovisity, and (d) slope degree.

3.5. Modeling and verification of ephemeral gully head advance rate

We developed a prediction model for the ephemeral gully head advance rate using catchment area, slope degree, rainfall erosivity, and FVC as independent variables. The model was constructed based on multiple regression and correlation analyses, using 70% of the ephemeral gully heads for model building. The developed model is expressed as: (8) $R_{\mathrm{EGH}}=9.26\times {10}^{-4}{S}^{0.187}{C}_A^{0.597}(R-1855.87{)}^{0.498}(1+\mathrm{FVC}{)}^{-0.565}(R^2=0.94)$(8) where $R_{\mathrm{EGH}}$ is the ephemeral gully head advance rate (m yr^–1), S is the slope of the approach channel above the ephemeral gully head (°), $C_A$ is the catchment area above the ephemeral gully head (m²), R is the annual average rainfall erosivity (MJ mm hm^–2 h^–1 yr^–1), and FVC is the fractional vegetation cover (%). The model exhibited a high coefficient of determination ($R^2$ = 0.94), indicating a strong correlation between the predicted and measured values.

Based on the prediction model, we randomly selected 30% of the ephemeral gully heads to calculate the error between predicted and measured values (Figure 9), resulting in a good linear fit ($R^2$ = 0.92) and an acceptable error (RMSE = 0.31 m yr⁻¹).

Figure 9. Ephemeral gully head advance rates from the field survey (measured values) and the prediction model (predicted values).

4. Discussion

4.1. Ephemeral gully head advance rate

The study revealed a wide range of ephemeral gully head advance rates (0.04–5.54 m yr⁻¹), with an average rate of 1.37 m yr⁻¹ from 2009 to 2021, highlighting the significant erosion dynamics occurring in the study area. Comparatively, ephemeral gully head advance rates in other semi-arid regions worldwide typically range from 0–3 m yr⁻¹, with rates less than 1 m yr⁻¹ recorded in Morocco, Tunisia, Romania, and Poland (Table 2). However, Ethiopia stands out with the highest recorded rate of 18.40 m yr⁻¹. The notably higher ephemeral gully head advance rates observed in the hilly and gully areas of the Loess Plateau in China compared to other semi-arid regions can be attributed to unique combination of soil properties and climatic conditions in the region. The loose and porous nature of the loessial soil, coupled with intense and concentrated rainfall events, contributes to the heightened susceptibility to erosion. Additionally, the large catchment area above the ephemeral gully heads exacerbates erosive processes, resulting in severe headcut erosion (Hu et al. 2019; Li et al. 2022; Wang et al. 2020; Zhao et al. 2021).

4.2. Driving factors of ephemeral gully head advance rate

The OPGD model identified catchment area as the most influential single factor impacting ephemeral gully head advance rate, with the interaction between catchment area and slope degree emerging as the primary driving force. High-risk areas, characterized by elevated ephemeral gully head advance rates, were mostly distributed in regions with larger catchment areas (>375 m²), low vegetation coverage (0–34%), and specific ranges of cumulative rainfall erosivity (1,905–1,908 MJ mm hm^–2 h^–1 yr^–1) or slope degree (26–27.5°). The catchment area plays a fundamental role in generating headcut erosion, with its size directly influencing the volume of runoff and flow kinetic energy, thus affecting the severity of headcut erosion. As a direct erosion factor, precipitation exacerbates soil erosion, particularly during high-intensity and high-frequency rainfall events (Liang et al. 2020; Zhao et al. 2020). Positive correlations between ephemeral gully head advance rate and catchment area and rainfall erosivity align with previous studies (Capra and Spada 2015; Han, Zheng, and Xu 2017; Li et al. 2016; Tang et al. 2022; Xu et al. 2017). Despite the uncontrollable nature of the monsoon climate on the Loess Plateau, larger catchment areas can be managed effectively through measures such as hedging or other plant engineering measures to disrupt catchment area connectivity, reduce runoff, and mitigate the impact of high-flow kinetic energy on gully head extension.

The longstanding ‘Grain for Green’ project on the Loess Plateau has yielded remarkable results (Fu et al. 2011; Ran et al. 2020; Sun et al. 2014), as evidenced by the negative correlation between ephemeral gully head advance rate and vegetation coverage in the study. Denser vegetation cover intercepts raindrops, reduces their velocity, and shields the soil surface from direct impact. Moreover, dense vegetation slows overland water flow, and root systems contribute to runoff reduction by enhancing soil characteristics (Chen, Guo, and Wang 2020; Dou et al. 2020; Hayas, Poesen, and Vanwalleghem 2017; Ran, Lu, and Xu 2013; Wang et al. 2014; Wang et al. 2021). Additionally, slope gradient significantly affects soil erosion, with steeper slopes dramatically increasing soil losses, particularly in areas with sparse vegetation cover (Liang and Fang 2021; Sun et al. 2014; Xu and Zhang 2020).

4.3. Prediction model for ephemeral gully head advance rate

Ephemeral gully headcut erosion is an indicator of severe land degradation, disrupting surrounding agricultural areas, altering watershed hydrology, and increasing sediment production (Guan et al. 2021). Thus, it is imperative to identify the driving factors behind ephemeral gully headcut erosion and develop a prediction model for ephemeral gully head advance rate to offer practical guidance for water and soil conservation management. While previous models such as SCS Ⅰ and SCS Ⅱ exclusively incorporated catchment area and rainfall variables, and the Thompson model included catchment area, slope degree, rainfall, and soil clay percentage to calculate ephemeral gully head advance rate (Table 3), our study revealed that four key variables – catchment area, slope degree, rainfall erovisity, and FVC – influence ephemeral gully head advance rate in the hilly and gully regions of the Loess Plateau. We integrated FVC into our prediction model, recognizing its significance as vegetation restoration efforts have substantially progressed in the hilly and gully areas of the Loess Plateau since 1999, rendering FVC an essential variable that cannot be ignored.

Table 3. Summary of global prediction models for ephemeral gully head advance rate

Study region	Prediction model for ephemeral gully head advance rate	References
USA	$R=0.15A^{0.49}{S}^{0.14}{P}^{0.74}E$ where $R$ is ephemeral gully head advance rate (m yr^–1), $A$ is drainage area above gully head (m^2), $S$ is slope of the approach channel above the gully head (%), $P$ is the sum of rainfall from 24-hour rains ≥12.7 mm, and $E$ is the soil clay percentage	Thompson (Thompson 1964)
USA	$R=1.5A^{0.46}{P}^{0.2}$ where $R$ is ephemeral gully head advance rate (m yr^–1), $A$ is drainage area above gully head (ha), and $P$ is the sum of rainfall from 24-hour rains ≥12.7 mm during gully development (mm)	SCS Ⅰ (USDA 1966)
USA	$R_A=R_1{\left(\frac{A_1}{A_2}\right)}^{0.46}{\left(\frac{P_2}{P_1}\right)}^{0.2}$ where $R_A$ is ephemeral gully head advance rate (m yr^–1), $R_1$ is past average annual ephemeral gully head advance rate (m yr^–1), $\frac{A_1}{A_2}$ is the ratio of the mean drainage area of the reaches to the average drainage area of the reach through which the gully has moved, and $\frac{P_2}{P_1}$ is the ratio of the amount of rainfall from 24-hour rains ≥12.7 mm during the study period	SCS Ⅱ (USDA 1966)
Hilly and gully region of Loess Plateau in China	$R=9.26\times {10}^{-4}{S}^{0.187}{C}_A^{0.597}(R-1855.87{)}^{0.498}(1+F{)}^{-0.565}$ where $R$ is ephemeral gully head advance rate (m yr^–1), S is slope of the approach channel above the ephemeral gully head (°), $C_A$ is catchment area above the ephemeral gully head (m^2), R is annual average rainfall erosivity (MJ mm hm^–2 h^–1 yr^–1), and F is fractional vegetation cover (%)	This research

Between 2000 and 2020, concerted efforts to implement water and soil conservation measures such as the ‘Grain for Green’ project, sloping land conversion program, and small watershed management led to significant vegetation recovery on the Loess Plateau (Fu et al. 2017; Jia et al. 2019; Ouyang et al. 2016; Zhao et al. 2013), with NDVI values increasing by 35.66% (Chen et al. 2023). This vegetation recovery has reduced water and soil losses in the region, as evidenced by the significant decrease (44%) in sediment transport at the main hydrological control stations in the Yellow River in 2020, as reported in the Chinese River Sediment Official Gazette 2020 (http://xxzx.mwr.gov.cn/xxgk/gbjb/zghlnsgb/), compared to the average value recorded from 1954–2020. Our findings demonstrate a negative correlation between vegetation coverage and ephemeral gully head advance rate, affirming the pivotal role of vegetation in mitigating erosion. By incorporating FVC into the prediction model, we observed improved accuracy, highlighting the significance of accounting for vegetation dynamics in ephemeral gully headcut erosion control strategies. Furthermore, these research insights offer valuable guidance for ephemeral gully headcut erosion control on the Loess Plateau and for addressing similar challenges encountered in other regions worldwide.

4.4. Limitations and future perspectives

While our study identified key natural factors influencing ephemeral gully headcut erosion, it is important to acknowledge certain limitations and consider future research directions. Firstly, factors such as land use land cover, soil type, and geology are known to play significant roles in affecting ephemeral gully headcut erosion (Sun et al. 2014; Xu and Zhang 2020; Zhang et al. 2017). However, these factors were not included in our analysis due to a lack of significant differences in their distribution across ephemeral gully headcut erosion location. Future research could explore integrating these factors into the modeling framework to provide a more comprehensive understanding of ephemeral gully headcut erosion dynamics.

Secondly, our study did not explicitly consider the impact of frequent and high-intensity rainstorm events on ephemeral gully headcut erosion. Given the increasing occurrence of extreme weather events in recent years, investigating the influence of rainstorm events on ephemeral gully headcut erosion dynamics represents an important avenue for future research.

Lastly, our prediction model for ephemeral gully head advance rate is based on data from a single watershed in the hilly and gully region of the Loess Plateau. Validating and refining the model using data from multiple regional watersheds would improve its robustness and generalizability.

5. Conclusion

Ephemeral gully erosion, particularly at the gully head, poses significant challenges to soil and water conservation efforts, especially in regions like the Loess Plateau. Our study focused on ephemeral gully headcut erosion in the hilly and gully regions of the Loess Plateau, shedding light on its spatial distribution and key influencing factors. Our analysis highlighted catchment area as the most significant single factor influencing ephemeral gully head advance rates. Additionally, our OPGD model underscored the importance of interactions between catchment area and slope degree, rainfall erosivity, and vegetation coverage in predicting ephemeral gully head advance rates. Incorporating vegetation coverage (FVC) into our prediction model further enhanced its accuracy. By developing a prediction model that incorporates these driving factors, including FVC, we aim to provide land managers with a valuable tool for implementing effective soil and water conservation measures. Strategies like hedges or other bioengineering measures can increase vegetation coverage, divide catchment areas, reduce runoff accumulation, and mitigate ephemeral gully head advancement. Our research contributes to a better understanding of ephemeral gully headcut erosion and offers practical guidance for land management planning in semi-arid and arid regions like the Loess Plateau. Moreover, the insights gained from our study can inform conservation efforts in other loess regions worldwide, aiding in the sustainability of these ecologically sensitive areas.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (41977064) and National Key R&D Program of China (2021YFD1900700). The authors express their gratitude to their colleagues in their research group for their help with the experiments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by National Key Research and Development Program of China [grant number: 2021YFD1900700]; National Natural Science Foundation of China [grant number: 41977064].