Monitoring and evaluation of gully erosion in China's largest loess tableland based on SBAS-InSAR

Gully erosion is widespread in central China's ecologically fragile loess plateau. However, research on the monitoring and evaluation of large-scale erosion is scarce. Here, we collected 16 pairs of Sentinel-1A images from the rainy season, spanning April 5-October 26, 2019, subsidence is then quantified and its results can indirectly reflect the intensity of erosion on the largest loess plateau in China, the Dongzhiyuan tableland, using small baseline subset synthetic aperture radar interferometry (SBAS-InSAR). The results showed that, the average erosion rate was 9 mm/year, throughout the study area, Surprisingly, the average erosion rate in the gully reached 160 mm/year, the sides of the gully correspond to the most extensive and intense erosion. Heavy rainfall, slope gradients below 15°, changes in groundwater levels due to freeze–thaw action and water infiltration during rainfall, and lack of vegetation due to continued gully side erosion, drive the intense erosion. These results demonstrate the reliable ability of the SBAS-InSAR method in accurately assessing large-scale gully erosion. The InSAR-calculated subsidence demonstrates a strong predictive capability for soil erosion and is generally consistent with previously observed levels of erosion, although much noise interference also occurs. Therefore, in the future InSAR-calculated sedimentation could be used as an indirect measure of erosion levels at multiple watershed scales.


Introduction
Gully erosion is the slow or abrupt subsidence of the surface of a gully caused by natural phenomena and human activities . Gully erosion is one of the most dangerous forms of soil erosion (Belayneh et al. 2020). The ecological resource crisis caused by gully erosion seriously hinders the social and economic development of the region (Shi et al. 2020). It has increased, for instance, the likelihood of other geological phenomena, such as mudslides and landslides, which can lead to human casualties and economic losses, especially among the most vulnerable populations . Gully erosion is considered a major challenge to land resource degradation and is becoming increasingly common (Guo et al. 2019). In recent decades, gully erosion has been recognized as an important environmental issue (Huang et al. 2023).
Especially located at the boundary of the Loess Plateau, with active geological tectonic activity, fragile regional stability, and porous geology, the Loess Plateau has long been considered as one of the most susceptible areas in the world to gully erosion (Zhuang et al. 2022). Determining methods for monitoring and analyzing gully erosion has been the focus of research, and traditional methods for monitoring gully erosion include field surveys, positional observations, and GPS methods. Traditional methods are labor-intensive, time-consuming, and only applicable to small areas (Addisie et al. 2017;Hu et al. 2022). Another option is to consider Light Detection and Ranging (LIDAR) (Castillo et al. 2012), which is more accurate and can be done quickly, but usually, it is done in small areas. It places the sensor less than one kilometer away from the target, resulting in a limited scale of evaluation. This is a key drawback as gully erosion can be governed by mesoscale processes (Haregeweyn et al. 2017), making the impact of localized instances of the study limited. Gully erosion has also been determined using DEMs, but this method has not been commonly used due to inherent accuracy limitations (Raj et al. 2022). 3D scanners have been successful in obtaining high-resolution DEMs that can be used to assess surface landforms and gully erosion. However, it is not suitable for assessing deformation and erosion on a fairly large regional scale . Instead, monitoring deformation near gullies over a large area (960 km 2 ) is crucial.
Compared to the gully erosion monitoring techniques described above, combining traditional aerial remote sensing monitoring methods (Unmanned aerial vehicles (UAVs)) with high-resolution satellite remote sensing imagery has distinct advantages, including the ability to rapidly monitor surface deformation and to produce large-scale and high-resolution products (Guan et al. 2021). In general, UAVs are more advantageous for short-term gully erosion monitoring, but they are also more sensitive to the environment. In addition, UAVs cannot be used to monitor the continuous and dynamic deformation of gully surfaces due to their limited flight range. Gully surface deformation contributes to gully erosion and thus is considered to be a measure of gully erosion. Monitoring the change of the gully surface with time is a direct way to obtain the deformation trend of the gully. Synthetic aperture radar interferometry (InSAR) is a high-resolution technique, 5 × 20 m spatial resolution, that is primarily used to monitor surface deformation over long time series and large areas ). In addition, multi-temporal InSAR is suitable for analyzing time series images (Wu et al. 2019). The small baseline subset InSAR (SBAS-InSAR) method (Berardino et al. 2002) and the permanent scatterer InSAR (PS-InSAR) method (Ferretti et al. 2000) have been widely used for surface deformation with millimeter-level measurement accuracy. Many studies have shown that SBAS-InSAR has higher data utilization aiming at deformation monitoring over large areas with distributed targets, and thus the SBAS-InSAR method is more suitable than other InSAR methods for solving low coherence problems in vegetated and gully regions (Pu et al. 2021). SBAS-InSAR has been used in a wide area for ground surface, especially urban surface subsidence (Pu et al. 2021), volcanic activity (Meng et al. 2022b), and glacial recession . However, it is rarely used to monitor gully erosion. the SBAS-InSAR method can fill the gap between historical hazard data and gully assessment on a regional scale (Su et al. 2022). Therefore, exploring the application of the SBAS-InSAR method in gully deformation can promote the spatial and temporal scales of gully deformation monitoring and expand the rapid assessment of gully deformation at the regional scale.

3
The Loess Plateau in China, covering an area of 640,000 square kilometers, is one of the most severely eroded areas in the world (Kou et al. 2020). Gully erosion on the Loess Plateau is a serious problem, mainly manifested in the widespread landslides on both sides of the gully (Wu et al. 2017). There is little research on gully erosion on a large-scale ). This study focuses on the Dongzhiyuan tableland, a borderline region on the easternmost edge of the Gansu Region in Central China. Gully erosion has fragmented the entire Loess Plateau, and the Dongzhiyuan tableland is the largest remaining loess tableland in China (Huo et al. 2020), and the continuous gully erosion and retreat here makes the tableland shrink year by year (Gao et al. 2016;Qin et al. 2022). Particularly, gully erosion has shown an increasing trend that has the sed attention of the local government (Guan et al. 2021) as it is becoming a notable factor for hindering economic development the environment conservation. Therefore, there is an expanding interest in monitoring the spatial and temporal distribution of gully erosion. This study aims at shedding some light on the time behavior of gully erosion in this region as a way to identify time-pervasive factors that enhance it, allowing for the structuring of remedial initiatives.
For this purpose, 16 SAR images taken by the Sentinel 1A satellite while it was in a descending orbit, were sourced. Then, vertical surface displacements were estimated using SBAS-InSAR over the whole study region. This outcome allowed for (1) computation of gully erosion rates through assessment of changes on the surface over time.
(2) Characterize locations where gully erosion is more predominant by considering slope and elevation as their most critical morphological conditioners. (3) Assessed the role of rainfall as a covariant in gully erosion rates, by comparing their time histories. This study will be beneficial to the study of gully erosion characteristics and deformation mechanism, and provide relevant reference for InSAR reliability verification.

Study area
The Dongzhiyuan plateau spans between 106°14′-108°42′E,34°50′-37°19′N, which corresponds to the southeastern edge of the Gansu Province, China (Fig. 1a), an arid and semi-arid zone. The topography of the study area is dominated by gully incised tablelands, on terrain that displays an elevation ranging between 882 and 1540 m.a.s.l with an average of 1300 m.a.s.l (Fig. 1b). The entire Dongzhiyuan covers an area of 2835 km 2 , while the tableland area is 960 km 2 . Its average annual temperature is 8.5 °C in 2019, and its annual accumulated precipitation is 672.3 mm (Fig. 2). Nevertheless, its average annual evaporation is 1600 mm (Xu et al. 2018). However, trends are not uniform either spatially or over time. Rainfall in the south of the plateau is larger than that in the north, and at least 70% of the yearly value occurs between July and September, the peak occurred in August with 159 mm of precipitation (Fig. 2). Furthermore, storms tend to be high-intensity and short-duration . As expected, vegetation cover mirrors geographical trends on rainfall, spices adapted to arid terrain are predominant in the north, while lush, more perennial vegetation is widespread in the south (Ling et al. 2022). The water table exhibits extensive variation due to the high time-clustering of rainfall. Particularly, it is usually deep, as evaporation exceeds ground percolation during most of the year, leading to low water content in the surface soil (Li et al. 2014). Top layers (150-200 m thick) in the plateau are comprised mostly of Quaternary Upper Pleistocene loess that displays an on-site density of 1.79-1.97 g cm −3 . Loess depth generally exceeds 100 m. Moreover, due to their loose nature, they are susceptible to collapse, leading to the formation of cavities and voids within them. The loess surficial strata are underlaid by a black kiln clay soil. This geomorphological arrangement is observed all over the plateau . It must be stressed that these geomorphological units showcase different erosion-resistant traits. The Heilu soil retains moisture and soil particles, while loess tends to decompose, especially in areas on both sides of the gully where loess crumbling and accumulation has occurred and is more likely to be washed away by rainwater, resulting in subsidence on the gully surface (Fig. 1d, e).

Generation of basic data with the aid of GIS
Diverse data sources were tapped for collecting the information required by this study. Firstly, Google Earth Pro was employed to collect images in the visual spectrum (0.4-0.7 µm) to get an overview of the study site and select features of interest. Furthermore, a Digital Elevation Model (DEM) with a 30 m resolution was sourced for U.S National Aeronautics and Space Administration (NASA) (http:// www. gsclo ud. cn/ search). It is useful for performing orbital interference de-duplication for InSAR analysis, providing ancillary data for geocoding and mapping sub-watershed in the study area. Both sources were integrated using ArcGIS 10.7, which was instrumental for mapping the study area, identifying gullies that were vectorized as lines, and allowing for the computation of geomorphological parameters using elevation data. Then, SAR images from the Sentinel satellite constellation were collected through the Copernicus hub (https:// scihub. coper nicus. eu/) to estimate displacement time histories using SBAS-InSAR. Further details about these datasets are presented in the following section. Obtained displacement data were interpolated to match the resolution of the DEM, allowing for exploration of dependency traits among surficial deformation, slope, and absolute elevation. The 2019 daily precipitation data from the Xifeng climate station located in the study area were provided by the National Climatic Data Center (NCDC) (https:// www. ncei. noaa. gov) to study the relationship between gully surface subsidence and precipitation.

Sentinel-1A data
As stated previously, SAR images were downloaded from the Copernicus hub, supported by ESA (https:// scihub. coper nicus. eu/). The analysis period ranged from April 2019 to October of the same year. Images selected for the analyses were taken each 12 days, after carefully verifying all shared the same incidence angle (37°). All were taken during the descending pass or the satellite and they correspond to the Vertical transmitting, Vertical receiving (VV) polarization mode, as shown in Table 1. The Interferometric Wide swath (IW) mode was used, which creates images with a 250 km swath at 5 × 20 m spatial resolution. Moreover, the accuracy was improved by considering orbital precision files ) (https:// scihub. coper nicus. eu/ gnss).

SBAS-InSAR Monitoring Principle
SBAS-InSAR is a procedure that uses SAR images to construct surface change time histories by grouping them into master and slave sets. This allows for effective data extraction from noisy sets (with low coherence) even when using a low number of images. Particularly, the following relationship holds: (Chen et al. 2020).
where N is the number of remote sensing images while M is generated considering N + 1 SAR images \ in such a way that M takes the value of 102 possible pairs in this study.
( Equation (2) represents the interferometric phase composition of image j (generated by two SAR images at times t A and t B , given that t A < t B ) within the pixel element at azimuth x and r distance from the satellite x, r. φ(t A , x, r) and φ(t B , x, r) are the phase values of the interferograms at time t A and t B at pixel point (x, r). Φ def , j (x, r) represents the phase difference between t A and t B while Φ topo,j (x, r) denotes the topographic residual phase after extracting the reference DEM elevation data. Φ atm,j (x, r) is the error caused by the atmospheric phase, while Φ orb,j (x, r) represents the orbital error, and Φ noise,j (x, r) represents random noise. Clearly, Φ atm,j (x, r) can be obtained from the interferometric phase by performing time domain high-pass filtering and spatial low-pass filtering. After removing the phase error caused by the atmosphere, the interferometric phase mainly consists of deformation and residual elevation.
After that, it is assumed that the velocity vector in a certain direction can be expressed by a linear model. To obtain the surface deformation information and the residual elevation DEM component, we can transform the Φ def,j (x, j) term in Eq. (2) as follows: where refers to the radar wavelength, B Vj refers to the vertical baseline of the interference relative to j, θ refers to the satellite incidence angle, and △z refers to the altitude difference. Also, it is possible to formulate the C vector as: Finally, the displacement can be computed in this manner: where D is an M × (N − 1) size matrix of known coefficients. Thus, the deformation rate can be obtained by performing singular value decomposition (SVD) on the result of expression 6. Figure 3 illustrates the flowchart of the data processing done in this study. The main steps are described as follows:

Data processing workflow
Step 1: We generated connection maps, determined the spatial and temporal baseline thresholds, divided the input 16 sentinel-1A SAR images into different subsets (spatial baseline of 93 m and temporal baseline of 220 days), and generated 102 high-quality interferometric pairs.
Step 2: The process of generating interferograms was the most critical step. The Goldstein filter was used in this process to improve the accuracy of phase unwrapping and measurement precision to maximize the signal-to-noise ratio of the interferogram. External DEM data were used to remove the topographic phase. The coherence threshold was 0.3 and the Minimum Cost Flow (MCF) is used for unwrapping.
Step 3: We used ground control points (GCPs) to eliminate the remaining terrain phase and satellite orbit faults. The GCPs should be located in the best-unwrapped phase region and away from the deformation region (GCPs should be chosen on the plane of the tableland in this study.).
Step 4: In this process, we estimated and removed the remaining constant phase and phase slopes that remained after the unwrapping. The deformation rate and residual topography were estimated using a linear model. Re-unwrapping was then performed using a minimum cost flow. A spatial-temporal filter mitigates the effects of atmospheric delay (spatial filter window: 1600 m; temporal filter window: 365 days) and obtained the average deformation rate and cumulative gully surface deformation along the line of sight (LOS). Considering the large area of subsidence and landslides on both sides of the ravine in the study area, we calculated the subsidence rate of the ravine area to indirectly reflect the intensity of erosion.

Gully erosion
The average elevation within the Dongzhiyuan tableland was 1247 ± 11 m, while the elevation of the bottom of the most severely eroded gully bottom within 880 m.a.s.l (Fig. 4a). Assessment of DEM data lead to an average slope of 29 ± 6° (Fig. 4c), which was slightly lower than expected. This outcome could be explained by the relatively low resolution of the DEM (30 m) and radar shadowing. Contrarily, the southern extent showed gentle slopes. The result indicated that the gully erosion rate in the Dongzhiyuan tableland reached a maximum of 160 mm/year while the maximum uplift rate was 79 mm/year (Fig. 4b). Also, there was a large variation in erosion rates between northern and southern regions. In the former it ranges between 160 and 36 mm/year and, as expected, it clustered around gullies (Fig. 4c). Particularly, erosion rates exceeded 100 mm/year within the Zhangjia mountain, Taoshu, and Liu gully (Fig. 4b). By comparing the different soil types, the results show that the Heilu soil type was mainly distributed in the tableland (Fig. 4e), where severe erosion rarely occurred, with an average deformation rate of 3 mm/yr (Fig. 4f). On the contrary, Fig. 3 Processing workflow of SBAS-InSAR method the average erosion rates of Alluvial soil, Lately deposited soil, and Huangmian soil were 16.4 mm/yr, 19.5 mm/yr and 12 mm/yr, respectively (Fig. 4e, f). Especially Alluvial soil and Lately deposited soil types were the most severely eroded, forming two huge gullies ( Fig. 4e; Fig. 4a). Besides, some uplifting was recorded, particularly in the center of the tableland, which may be related to intensive engineering activity. The average rate in the study area was 12 mm/year while the maximum rate on the tableland reached 79 mm/year. More details could be found by focusing on specific locations, such as Areas A, B, and C highlighted in (Fig. 4b). The most intense erosion in Area A were at the bottom and sides of the gully (Fig. 5a). The cumulative erosion from April to October at points P2 and P3 reached 52 and 41 mm, respectively, especially from August 15 to August 27, when P1, P2, P3, and P4 all observed the most severe erosion (Fig. 5d). Another Zone B with very severe erosion was located in the Qianyuan gully (Fig. 5b). Although the average cumulative erosion across the four observation sites from April to October reached 27 mm, its maximum deformation was observed between June 16 and June 28, when the average erosion was 16 mm (Fig. 5e). Area C is located west of Area B (Fig. 5b) and erosion occurred Topography and deformation of the study area. a DEM distribution in the Dongzhiyuan tableland; b average deformation rate map of the study area from April to October 2019; c distribution of gully landslides within the study area; d cumulative deformation rate of the study area from April to October 2019; e Soil types in the study area; f Distribution of average deformation rate density by soil types mainly in the gully area (Fig. 5c). From April to October, the average cumulative erosion at the four observed sites in Area C reached 32 mm (Fig. 5f). However, severe erosion of these points occurred mainly from April 5 to July 22. Thereafter, erosion gradually decreased from August to October (Fig. 5f). Notably, the trend reversed in August, with a slow increase observed in both. However, the gains were reversed due to the more pronounced erosion recorded on October 14. The results for the other trenches also showed these changes, indicating that erosion and uplift were not linear phenomena but followed each other, resulting in oscillatory changes in vertical displacement.
Contrarily, uplift was the major trend observed at the Xiapan gully head near Yima town (Fig. 6a), exceeding 30 mm from April to October (Fig. 6b). Particularly, cumulative uplift reached 32 and 25 mm at points A and B, between April 5 and May 23. This was the largest uplift rate recorded in this study. However, the uplift moderated afterward, and between June 16th and August 15th, it was almost zero. Then vertical movement resumed and eventually total uplifted exceeded 40 mm in both locations from June 16th (Fig. 6b). Anomalous uplift values were observed at the head of the trench close to the city and the deformation values were irregularly distributed. Considering that because the ditch head of  (Fig. 4b). a, b, c perspective views of the average deformation pixel values of Yinshan, Qianyuan, and Dongzhuang gully areas superimposed on the DEM, respectively; d, e, f time series deformation of the three gully severely eroded regions the Dongzhi tableland surface close to the city was intensively planted with vegetation by humans, the vegetation developed luxuriantly in summer (Fig. 6d,e), which led to the instability of the surface deformation results (Meng et al. 2022a). Finally, partial uplift appeared in the InSAR results.

Topography and erosion
The elevation within the Dongzhiyuan tableland ranges from 1200 to 1400 m.a.s.l (Fig. 7a), while the elevation of gullies within the study region is between 882 and 1540 m.a.s.l. Moreover, we are surprised to find that the higher the altitude, the smaller the erosion rate, which also indicates that the erosion is mainly driven by the downstream concentrated runoff. It was notable that erosion at elevation 900 m.a.s.l appeared to be more severe than that at this particular location, as the average cumulative deformation was 28 mm/yr and 17 mm/yr, respectively (Fig. 7c, e). In the slope interval of less than 5°, the larger the slope, the greater the erosion rate (Fig. 7d). This indicated that the slope of the gentle slope was an important factor to control the erosion rate of the gully. However, the local especially  Correlation between gully erosion rate and slope and elevation in the study area. a Elevation values in the study area; b slope distribution in the study area; c relationship between mean deformation and cumulative deformation and elevation, and the red dashed line marks the outliers of the relationship; d coupling relationship between slope and deformation rate; e The average deformation rate of the points with elevation at 900 m.a.s.l superimposed on the Google Earth satellite image steep terrain seemed to be unfavorable to erosion, which was related to the strong stability of this vertical loess .
There were also spatial effects on the erosion rates within gullies, as erosion tends to cluster in the northern and eastern regions of the Dongzhiyuan tableland. We speculated that this was related to the catchment area of the basin (Fig. 1b, c). A larger catchment area tended to collect more runoff during rainfall (Fig. 8a). Therefore, the focus was brought to the gully bottom of three large basins, as shown in Fig. 8a. The largest erosion rate (41 mm/year) was observed close to the bottom of the gully near the drainage outlet of the basin (Fig. 8b). Then there was a sudden reduction in erosion rates at a location 3.5 km away from the bottom, where it eventually turned to be less than 20 mm/year close to gully head. For the second case, erosion at the bottom was observed to be 37 mm/year at the bottom while it was observed to be 17 mm/year at the top (Fig. 8c). Finally, the erosion rate reached 47 mm/year at the bottom of Gully 3, decreasing to 26 mm/year. Besides, Gully 3 was similar to Gully 1 in that the erosion rate abruptly increased rapidly near the water outlet. This indicated that the increase in erosion rate at the ditch bottom was not linear, but pulse, which may be related to the flow collection of tributaries in different sub-basins (Fig. 8b, d, e).

Precipitation and gully erosion
Gully erosion typically responds to precipitation, which causes an increase in loose soil pore-water pressure and reduces the effective normal stress along the basal shear zone. To analyze the temporal correlation between gully erosion and precipitation, and taking into account the re-entry cycle of the Sentinel-1A satellite, we counted the time series of Fig. 8 Statistical analysis of the correlation between the average erosion rate and elevation changes in the gullies. a Altitude data in the study area and the spatial locations where gullies 1, 2, and 3 are located; b, c, d Altitude data obtained from gullies 1, 2, and 3 at intervals of 300 m, 450 m, and 380 m along the red line and the average erosion rate in the area; b, c, d The x-axis represents the collection distance, the y 1 -axis represents the altitude data, and the y 2 -axis represents the average gully erosion rate; e Study area subbasins gully surface deformation and rainfall in areas A, B, and C using a 12-day unit period. The cumulative rainfall in the study area between April and October was 672 mm. The 12-day average rainfall was 40 mm (Fig. 9). The observed rainfall showed a seasonal behavior. Most of the rainfall occurred during the summer months, from July to September, accounting for more than 70% of the total. At the same time, the deformation showed a cyclical oscillation similar to that of precipitation. The three regions exhibited two distinct erosion phases during the study period. The first phase was from April 5 to August 3, when intense erosion occurred in all three regions, with erosion amounts reaching 42, 38, and 34 mm, respectively. The rainfall and erosion rates were synchronized. The second stage was from August 15 to October 26, when the erosion rate was slower. In particular, the rainfall was 75 mm on September 8, reaching a recent peak, but the erosion rate of the gully did not change accordingly and still maintained a relatively slow rate and even a weak uplift was observed in areas A and B. Gully erosion, as a complex geological event, is not only related to the change of rainfall pore pressure but also to some other factors.

Investigation of gully erosion processes in four soil textures
Unsurprisingly, soil texture affects critical shear stress, and factors such as soil structure and chalk content affect the resistance and lift experienced by individual particles (Douglas-Mankin et al. 2020). In the present study area, the Heilu soil has an internal formation of agglomerates and agglomerates due to long-term cultivation resulting in soil mixing with large amounts of plant residues. These types of sediments are less susceptible to entrainment and have a more compact structure, reducing the exposure of larger particles to hydraulic forces. Huangmian, Alluvial soil and Lately deposited soil have high chalk Fig. 9 Average deformation and rainfall for a 12-day cycle, during the investigation content, loose soils with very weak impact and erosion resistance, and are highly susceptible to erosion (Xiong et al. 2023). In contrast, Alluvial and Lately deposited soil is mainly slope deposits and flood deposits originating from nearby terraces, and is distributed by hydraulic transport, gravitational accumulation mainly throughout the gully bottoms and hillslope, and has a very unstable soil structure, leading to easier erosion (Fig. 4e). Moreover, the increase in soil sand content makes the catchment more erosive, leading to an increase in gully erosion rate sedimentation and a decrease in gully head height . Soil cohesion decreases with increasing sand content; therefore, it seems reasonable that soils with higher sand content may be eroded at lower shear stresses, which also leads to differences in erosion rates across the gullies in the study area. Also, it is more urgent to control the erosion of Huangmian, Alluvial soil and Lately deposited soil, namely to enhance the management of erosion in gully areas.

Steep terrain drives intense gully erosion
The results indicated that the lowest elevation location between 884 and 900 m.a.s.l in this study had the highest erosion rate (Fig. 7c), which is consistent with previous studies (Zabihi et al. 2018). That meant potential gully erosion was more severe at lower elevations (Dickson et al. 2007). In addition,  concluded that the closer to the river, the more severe the erosion, which is also in general agreement with our findings (Fig. 8), and we speculate that this is related to water infiltration in the gully. However, the main erosion control measures currently remain in the form of vegetation greening and terracing (Bastola et al. 2018), while monitoring and management of key erosion areas on both sides of the gully are neglected.
The steeper the slope corresponds to more intense erosion (Wang et al. 2016), but the correlation between soil erosion rate and slope weakens when the slope exceeds a certain threshold, in particular, once the gradient of slope in our study exceeds 16°, which is a specific value depending on the resolution of the slope, steep slopes seem to hinder erosion (Fig. 7d), which may be related to the property that the steeper the loess, the more stable the properties . Soil structures on hillslopes are more fragile and unstable to withstand shear forces and have more dramatic deformation patterns compared to gully bottoms. On the other hand, soil erosion increases and then decreases with slope, with the most intense erosion occurring between 5 and 20° , and we obtained a similar result (Fig. 7d) In general, for our study area, the slope varies significantly, with the northern location within the study area nearly 35°, favoring further erosion, leading to the further development of the gully, which in turn accelerates erosion of the tableland margin.

Water infiltration drives gully erosion
Runoff from rainfall is considered to be the most intuitive factor influencing soil erosion (Sun et al. 2021). Before August 3, the cumulative precipitation increased and the stability of the slope decreased significantly, which means that the source material increased rapidly. In particular, the strong rainfall that occurred during the drought area had a greater impact on slope stability (Fig. 9). Then, the accumulated rainfall increased, but the influence on the stability of the slope gradually decreases. Our results showed that from August 27 to September 20, the cumulative amount of continuous heavy rainfall reached 140 mm, but the erosion rates did not increase with cumulative rainfall. We confirm that rainfall events in the study area have a more significant effect on erosion in the gully during the early rainfall than during the late rainfall, furthermore, the erosion rate is most accelerated by extreme rainfall (Fig. 9). Zhou also stated similar findings (Zhou et al. 2022). High-intensity rainfall, over-irrigation of agricultural fields, and freeze-thaw cycles induce gradual groundwater infiltration within the upper loess layer, and this water flow generally overflows from slopes (Zhou et al. 2014). The same mechanism leads to pothole erosion and increases the chances of hillslope failure of the walls of gullies within the study area. Freezing-thaw effects are more predominant at the end of winter in March, leading to increases in groundwater levels (Chou and Wang 2021;Pei et al. 2023), which accelerate gully erosion. Likewise, over-irrigation can also compromise slope stability within the study area, by allowing infiltration right into the loess layer through crevices , eventually leading to an overall loss of strength. Specifically, outcomes of this research indicate that erosion suddenly becomes severe at the hillslope (Figs. 5a, 10a), a location where seepage of infiltrated groundwater leaks out, as it is being pressed by most of the hillslope weight (Fig. 10a). This was further verified by analysis of high-resolution images that showed severe seepage precisely at that height (20 m) above the gully bottom (Fig. 10b, c, d). In particular, seepage was prevalent in July after extreme rainstorms (Fig. 9), as infiltration continues and the water table rises, this will further exacerbate the landslides on both sides of the gully and the new accumulation will be washed away in the coming rainy season, which is a major cause of severe erosion on both sides of the erosion gully. In order to alleviate this problem, we recommend that vegetation restoration, especially herbaceous planting with a very shallow root system that does not damage the structure of the loess, should be implemented with emphasis on the sides of the gullies, which, in addition to avoiding the direct erosion of the loess by scouring, herbaceous plants can also consume some of the water flow within the soil.

Conclusion
We demonstrated that the SBAS-InSAR method can quantify and assess erosion by monitoring subsidence. In particular, the sides of gullies are the most eroded locations, which are closely related to rainfall, topography, water infiltration, and lack of vegetation cover. The most worrying scenario is that erosion will be further enhanced by continued topographic evolution in the northern part of the study area, a vicious cycle in which erosion in the northern part of the study area could fall into a permanent state of high erosion rates, leading to further shrinkage of the largest loess plateau in China. Our analysis suggests that a focus on enhancing the planting of herbaceous vegetation on the sides of the gullies, rather than focusing on greening the plateau and terraces, prevents the soil from being moved away from its initial position by depleting the water table and direct scouring by water flow. Our recommendations are only a first step in the fight against widespread loess gully erosion and in the future should focus on nature-based defense strategies, taking into account the economic costs and the difficulty of implementing measures, issues that should be given high priority in the research agenda.
Author contributions Yuxiang Tao contributed to conceptualization; Xiaobo Luo contributed to data curation; Yongcheng Gou contributed to investigation; Haibo Tian and Pinglang Kou contributed to methodology; Haibo Tian and Pinglang Kou contributed to software; Haibo Tian contributed to visualization; Haibo Tian contributed to writing-original drafts; Yuxiang Tao, Pinglang Kou, Andres Alonso, Haibo Tian, Chenyu Gong, Yunpeng Fan, Changjian Lei contributed to writing-review and editing. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.