The Gutun catchment is located in Ganguyi Town, northeast of Yan'an city in Shaanxi province. It is a typical hilly and gully loess region spanning an area of 24.30 km² and stretching for about 12 km from the northwest to the southeast. The main landforms are gully slopes, in addition to ridges and the beds of gullies, with an altitude of 900-1265 m. The area is characterized by a semi-arid continental monsoon climate, with an average precipitation of approximately 530.0 mm that is concentrated in the period from July to September of each year. According to data on daily meteorological observations recorded at the Yan'an station (a national primary standard weather station), and obtained from the China Meteorological Data Service Center (http://data.cma.cn/), the area received 568.0 mm of rainfall in July 2013, far exceeding those in the previous and subsequent years (114.1 mm each year on average in 1960-2012, and 141.0 mm on average in 2014-2022) (Figure 1). There were 20 rainy days in July 2013, including six days of torrential rain and two days of heavy rain.¹(1 The China National Standard on the Grade of Precipitation (GB/T 28592-2012) has conditional standards of precipitation within 12 and 24 hours. In this study, we classified the grade of precipitation only according to the 24-hour standard based on the collected data on the daily precipitation.) In addition, the daily precipitation on 12 days in July 2013 exceeded the average daily precipitation values in July 1960-2012 and 2014-2022. Moreover, there were nine days in July 2013 when the volume of daily precipitation exceeded the average daily precipitation plus the standard deviation in July in 1960-2012 and 2014-2022. There were 31 days in 2013 on which the daily precipitation exceeded the monthly average daily precipitation on rainy days in 1960-2012 and 2014-2022, and 15 days when it exceeded the monthly average daily precipitation plus its standard deviation on rainy days in 1960-2012 and 2014-2022 (Table 1). Therefore, the prolonged heavy rainfall in Yan'an in July 2013 provides a valuable period in which to examine topographical changes triggered by the extreme precipitation event.

With the development of high-resolution satellite remote sensing, stereoscopic photogrammetry based on the parallax analysis of stereo-pair images has emerged as an effective and efficient method for the generation of data for DEM (Li et al., 2018). Given the availability of data covering the entire catchment, we used DEM data on the Gutun catchment at a spatial resolution of 5 m that were processed from stereo-pair images, captured on May 29 and November 27, 2013, that were acquired from a three-line camera (TLC) sensor on the ZY-03 satellite platform. ZY-03 is China's first civilian high-resolution stereoscopic Earth-mapping satellite. Launched in January 2012, it has been widely used for surveying and mapping land resources, analyzing agricultural and forest irrigation and drainage, and for emergency mapping (Tang and Wang, 2016). Data on the Normalized Difference Vegetation Index (NDVI) were computed by using Landsat 7 imagery from June 2013, and were resampled to a spatial resolution of 5 m.

The relative elevation, reflecting the height of the slope with respect to the valley, is used to determine the gravitational potential energy that is relative to the stability of the slope and soil erosion. The Height Above the Nearest Drainage (HAND) index, which captures the vertical path of flow to the nearest drainage (Cao et al., 2019), was used to quantify the relative elevation. The position of the slope is the profile of the landscape on the slope and is related to its stability and soil erosion on it. The Relative Position Index (RPI) is an estimate of how far a given location is from a ridge or a valley. Its value is the ratio of the distance to the nearest ridge to the sum of the distances to the nearest streamline and ridge (Qin et al., 2009). The position of a given slope is classified as a valley, lower mid-slope, mid-slope, upper mid-slope, or ridge according to the value of the RPI (Skidmore, 1990). The HAND and RPI were extracted from the DEM data by using the AutoFuzSlpPos algorithm (https://github.com/lreis2415/AutoFuzSlpPos) (Qin et al., 2018; Zhu et al., 2018).

We used two classifications of the slope aspect: eight directions (i.e., north, northeast, east, southeast, south, southwest, west, and northwest), and four orientations (i.e., shady slope, half-shady slope, half-sunny slope, and sunny slope) (Pan et al., 2017). The shape of the slope was classified as convex if it had a positive curvature and concave if it had a negative curvature. The aspect, shape, and inclination of slopes in the Gutun catchment were calculated based on DEM data by using the Aspect, Curvature, and Slope tools in ArcGIS software. The slope length of the catchment was calculated from DEM data by using the Slope Length tool of SAGA software. The changes in the elevation, inclination, and aspect of slopes were calculated based on the corresponding raster data from two periods through the Raster Calculator tool of ArcGIS software. The elevation, relative elevation, RPI, length, inclination, aspect, and shape of slopes, and their changes in the subareas were calculated from the corresponding raster data by using MATLAB.

The correlation coefficients between the changes in elevation and topographical factors can help clarify the relative importance of the latter and reveal the mechanism of the topographical changes. They were calculated as the Pearson coefficient by using the corr(X, Y) function in MATLAB. To obtain the correlation coefficients, we quantified the categorical factors (i.e., slope shape and aspect, and gully land consolidation). The value of the curvature was used to represent the shape of the slope. The slope aspect along the eight directions was assigned values ranging from one to eight in order of the azimuth, while the slope aspect with respect to the four orientations was assigned values ranging from one to four based on the intensity of solar radiation. Gully Land Consolidation (GLC) was assigned the value one if it was within the areas of the land consolidation project and was otherwise set to zero.

To better highlight the relationship between the changes in elevation and factors of the terrain before the rainy season, we designed a Comparative Advantage Coefficient (CAC), which is the ratio of the percentage of the class of a given terrain-related factor (e.g., elevation or slope) in a given class of changes in elevation to the percentage of the class occupied by the given factor in the catchment. The CAC was calculated by using the following equation:

where CAC_i,j is the comparative advantage coefficient of class i of the terrain-related factor before the rainy season within class j of the changes in elevation, A_ij is the area of class i of the terrain-related factor before the rainy season within class j of the changes in elevation, m is the number of the classes of terrain-related factors, and n is the number of the classes of changes in elevation. If the CAC was greater than one, this means that the particular class of changes in elevation was more concentrated in the particular class of terrain-related factors. The CAC was designed to answer the question of where the changes in elevation had mainly occurred.

The sunny or shady properties of the slope aspect represented the key factor influencing topographical changes, followed by the height (i.e., elevation and relative elevation), inclination, and shape of the slope. The slope length, RPI, NDVI, and GLC played minor roles in topographical changes. A strong correlation (r = 0.75, p < 0.01) was observed between the elevation and the relative elevation owing to the small drop in elevation in valleys over a small spatial scale. The relative elevation was more strongly correlated with the other parameters (except for slope inclination) than the elevation, especially with the RPI, NDVI, and GLC. The length, shape, and aspect along eight directions of the slope had weaker correlation with one another than with the other parameters, except for the correlation between slope length and shape (r = −0.19, p < 0.01) as well as the correlation between the two classifications of the slope aspect (r = 0.20, p < 0.01). The RPI and slope inclination were moderately correlated with the NDVI and GLC (r was in the range of ±0.22-0.28, p < 0.01), while the NDVI was strongly correlated with the GLC (r = −0.54, p < 0.01) (Table 2).

We statistically analyzed the characteristics of topographical changes in the Gutun catchment before and after the rainy season in terms of the elevation, slope shape, and slope inclination. The main height, inclination, and morphology of the slope at points where the elevation changed were roughly consistent with the results of statistical analyses based on field surveys (Zhang et al., 2007; Wang et al., 2015; Zhu et al., 2017; Zhuang et al., 2018).

First, we used the relative elevation to approximately reflect the height of the slope. Our results showed that the rainy season reduced the elevation of the catchment above a relative elevation of 70 m. Zhang et al. (2007) used remote sensing data and a field survey to determine that landslides mainly occurred in the range of relative elevation of 50 m and 120 m in Baota in Yan'an.

Second, we found that the peak slope inclination was mainly distributed in the range of 24°-36° when the elevation decreased over different ranges. Wang et al. (2015) conducted a field investigation in Yan'an in July 2013, and found that most loess landslides occurred on slopes with angles in the range of 35°-50°. Zhuang et al. (2017) also analyzed the characteristics of loess landslides through a field investigation in Baota after a rainstorm and found that the highest density of landslides occurred over a range of slope angles of 25°-45°.

Third, we found that concave slopes dominated the morphology of areas with decreased elevation after rainfall while convex slopes occurred in areas with increased elevation. Different views have been proposed regarding the relationship between the slope shape and topographical stability. Wang et al. (2015) investigated geological disasters after the July 2013 rainfall event in Yanchuan county of Yan'an and found that linear and convex slopes (about 81% of the investigated loess slopes) were more prone to loess landslides than bench-shaped and concave slopes. Zhuang et al. (2018) found that concave slopes were prone to collapse and filling because their surfaces had steep upper-slopes with a higher potential energy. Our statistical analysis showed that the collapse of concave slopes reduced elevation while deposition at the bottom of convex slopes increased it in the Gutun catchment.

We used four elevation-related factors (i.e., elevation, relative elevation, slope length, and RPI) to represent slopes. The elevation is the height above sea level, the relative elevation is the height above the nearest drainage, the slope length is the distance to the ridge, and the RPI is the position between the nearest ridge and drainage. The elevation is an absolute indicator, while the relative elevation, slope length, and RPI are relative indicators of the slope.

The elevation can vaguely reflect the position of a given location on the slope and in the catchment. The relative elevation represents the height of a location on the slope but not in the catchment and is related to the gravitational potential energy that determines the stability of the slope. The relative elevation is more suitable than the elevation for analyzing topographical changes in the slope. The values of the correlation coefficient showed that the relative elevation was more strongly correlated with the other parameters considered here than the elevation except for the inclination of the slope (Table 2). This is possibly because slopes in the upper catchment were steeper than those in the lower catchment.

The RPI represents the position of a given location on the entire profile of the slope and is related to the positions of the slope (e.g., ridge, shoulder slope, backslope, foot slope, and valley). This indicator can be used to quantify the position of the slope while ignoring the impact of its height. The RPI did not exhibit a strong correlation with topographical changes but was suitable for analyzing the NDVI and GLC (Table 2). The slope length is among the main and most variable components of empirical models of soil erosion (Liu et al., 2000). It exhibited a weak correlation with topographical changes (Table 2), probably due to deposition on the foot slope.

The soil on shady slopes on China's Loess Plateau is moister than that on sunny slopes due to the effects of solar radiation and evaporation, because of which the vegetation on shady slopes is denser (Yang et al., 2012). Shady and half-shady slopes are thus more likely to be the habitat for arbor vegetation (Singh, 2018). The areas occupied by shady and half-shady slopes with values of the NDVI above 0.7 were larger than those occupied by sunny and half-sunny slopes (Figure 14). The correlation between the NDVI and the slope aspect along four orientations was −0.12 (p < 0.01), thus confirming that vegetation on sunny slopes was sparser than that on shady slopes (Table 2).

The vegetation root, which creates preferential paths of flow, plays a leading role in erosion that leads to landslides, and determines its depth and intensity (Guo et al., 2020; Zhuang et al., 2022). The depth of infiltration of moisture in the soil is positively correlated with the depth of roots (He et al., 2020). The topsoil of shady and half-shady slopes, with denser vegetation and moister content, in the rainy season in the study area in 2013 was more easily saturated such that it slid, and this led to landslide mass being deposited at the bottom of the slopes. The areas of increased elevation were located in the south, near the valley lines (Figure 3), forming the changed slope profile with the decreased elevation of ridges/upper mid-slopes and the increased elevation of mid-slopes/lower mid-slopes (Figure 7). There are two possible reasons for the smaller deposition on sunny and half-sunny slopes: (1) Such slopes might have undergone surface erosion owing to their relatively sparse vegetation and dryer soil. (2) Sunny and half-sunny slopes might have been subjected to shallower landslides owing to the shallow-rooted vegetation on them, such as herbs and shrubs.

Gully land consolidation can reduce the risk of landslides. The Gully Land Consolidation Project was initiated in some counties of Yan'an city in 2010 (Liu et al., 2015). Land consolidation was implemented in the main gully of the Gutun catchment to increase the area of arable land, and to balance economic development with eco-environmental protection by cutting off hilltops and filling gullies in 2011 and 2012 (Cao et al., 2023). The project included three types of areas: the main gully-filling area, the lateral gully-filling area, and the area of vegetation restoration (Figure 15). The main gully was wider and more important than the lateral gully, and this area was thus equipped with wider channels for flood drainage, and complete cutting-slope protection engineering was applied to it. The area, and its ratio of levels of changes in elevation showed that the area percentage of the main gully-filling area was the smallest, while the area percentage outside the project area was the largest in terms of a reduction in elevation greater than 0.5 m. The area percentage of decreasing elevation below 0.5 m in the main/lateral gully-filling area was larger than the area outside. This might have been related to the excavation of unstable slopes and the filling of gullies in the project. The area of vegetation restoration had a greater ratio of sedimentation than the area outside, where this might have been related to the downslope segment in which the area was located. Values of the correlation coefficient showed that the project had little effect on topographical changes in the region (r = 0.02, p < 0.01) (Table 2). Notably, the vegetation was negatively influenced by the project (r = −0.54, p < 0.01), probably because it was designed to consolidate the valleys and barren land.

We have explored the characteristics of topographical changes based on DEM data from two periods, before and after the rainy season in the Loess Plateau. We plan to improve our method and model in a number of ways in our future research.

We plan to use a greater volume of relevant data in future work. We focused here on the topographical changes triggered by heavy rainfall in July 2013, and thus the DEM data should have been considered from the periods before and after rainfall. However, we used DEM data processed from stereo-pair images obtained on May 29 and November 27, 2013 as this was the closest available dataset before and after the rainy season. In addition, owing to the lack of fine-grained data on rainfall, we were able to classify the grade of precipitation only according to the 24-hour standard. We can obtain more accurate results in future research by using more fine-grained datasets.

We also need to provide a more sophisticated analysis. The vegetation cover is an important factor influencing the occurrence and movement of rain-triggered landslides (Glade, 2003). We conducted only a coupling analysis of the NDVI and the shape of the slope based on a remote sensing dataset. On the one hand, more datasets can be coupled with topographical changes, including land cover/land use data and data from field surveys. On the other hand, a more detailed examination than that based on a two-factor coupled analysis could be carried out in future work to obtain more accurate and useful findings. Furthermore, we need to better distinguish between the (interactive) contributions of different drivers and factors on topographical changes.