Long-Term Gully Erosion and Its Response to Human Intervention in the Tableland Region of the Chinese Loess Plateau

The gully erosion process is influenced by both natural conditions and human activities on the tableland region, the Chinese Loess Plateau, which is a densely populated agricultural area with unique topography. For the purpose of assessing long-term gully growth rates, the influencing factors and potential of gully growth, KH-4B satellite images, Quickbird-2 images, and unmanned aerial vehicle (UAV) images were used to assess gully erosion from 1969 to 2019. The effects of runoff, topography and human activities were analyzed with information derived from historical and present images. Ninety-five investigated gullies were classified into four types: 45 growing, 25 stable, 21 infilled and four excavated gullies. The rates (RA) of 45 growing gullies ranged from 0.50 to 20.94 m2·yr−1, with an average of 5.66 m2·yr−1 from 1969 to 2010. The present drainage area, local slope, average drainage slope, annual runoff, and ratio of the terraced area were all significantly different between the stable and growing gullies. The long-term gully growth rate could be estimated using a nonlinear regression model with annual runoff (Qa) and the slope of the drainage area (Sd) as predictors (RA = 0.301Qa0.562Sd, R2 = 0.530). Based on the Sg-A and Sg-Qa relationship that was used to reveal the threshold conditions for gully growth, all growing gullies still have the potential to keep growing, but soil and water conservation measures, including terraces, could change the threshold condition by reducing the effective drainage area. The results of this study could be helpful for preventing further gully erosion by dealing with gullies far above the threshold line.


Introduction
Gully erosion is the process whereby water accumulates and often recurs in narrow channels and, over short periods, removes the soil from this narrow area to considerable depths, as defined by the Soil Science Society of America (https://www.soils.org/ publications/soils-glossary/#, accessed on 1 October 2021). Gully erosion represents a minimum of 10% to a maximum of 94% of the total sediment yield caused by water erosion [1], which is a major driver of land degradation on the global scale [2][3][4][5]. To control gully erosion, understanding gully head initiation and growth is necessary [1][2][3][4]6,7]. Gullies were developed only when the intensity of concentrated (overland) water flows during a rain event exceeded the threshold [1]. However, flow hydraulics for gully initiation are constrained in field measurements [1,3]. The critical environmental conditions of gullies were assessed in terms of precipitation [8], topography [7,9], soils and land use [10,11]. These factors control either the runoff hydraulics (e.g., precipitation, topography) or the resistance of the soil surface to incision (e.g., soils) or both (e.g., land use) [1].
Frankl et al. [12] found that the average headcut growth rate during the 2010 rainy season was 0.34 m·yr −1 , compared with 3.8 m·yr −1 during medium to long time periods

Study Area
The Zhongduo tableland (110.67-110.75°E, 35.97-36.04°N; 32 km 2 ) and the Shiquan tableland (110.60-110.66°E, 35.97-36.03°N; 20 km 2 ) on the southeastern Loess Plateau were selected as study areas; they are located in Jixian County, Shanxi Province, China ( Figure 1). The loess cover is approximately 150 m thick, and the main soil type is loessial carbonate cinnamon soil with weak alkalinity in the study area. There is a typical temperate continental monsoon climate with an average annual precipitation of 507.47 mm and an average annual erosive precipitation (precipitation > 12 mm·d −1 ) of 314.71 mm (1957-2019) (based on meteorological records from the China Meteorological Data Service Centre [CMDC], 9 August 2020, http://data.cma.cn/). Approximately 70% of precipitation events occur from June to September. People mainly live on the flat surface of the tableland, which has agricultural advantages. Most of the tableland area was reclaimed as farmland in the study area. Before the implementation of the Grain for Green Project (GFGP) in 2000, the main land use types were croplands for growing wheat (Triticum aestivum), corn (Zea mays), soybeans (Glycine max), and flue-cured tobacco (Nicotiana tabacum). Since 2000, most of the croplands have been converted into apple orchards for both

Study Area
The Zhongduo tableland (110. .75 • E, 35.97-36.04 • N; 32 km 2 ) and the Shiquan tableland (110.60-110.66 • E, 35.97-36.03 • N; 20 km 2 ) on the southeastern Loess Plateau were selected as study areas; they are located in Jixian County, Shanxi Province, China ( Figure 1). The loess cover is approximately 150 m thick, and the main soil type is loessial carbonate cinnamon soil with weak alkalinity in the study area. There is a typical temperate continental monsoon climate with an average annual precipitation of 507.47 mm and an average annual erosive precipitation (precipitation > 12 mm·d −1 ) of 314.71 mm (1957-2019) (based on meteorological records from the China Meteorological Data Service Centre [CMDC], http://data.cma.cn/, accessed on 9 August 2020). Approximately 70% of precipitation events occur from June to September. People mainly live on the flat surface of the tableland, which has agricultural advantages. Most of the tableland area was reclaimed as farmland in the study area. Before the implementation of the Grain for Green Project (GFGP) in 2000, the main land use types were croplands for growing wheat (Triticum aestivum), corn (Zea mays), soybeans (Glycine max), and flue-cured tobacco (Nicotiana tabacum). Since 2000, most of the croplands have been converted into apple orchards for both economic benefits and mitigating soil erosion. The valley area is dominated by forestland and shrubland due to natural vegetation restoration. The Gully Stabilization and Tableland Protection Project (GSTP) was launched in the 2010s to stabilize the channel and control gully erosion around the tableland, which includes water barriers, terraces, check dams, etc.

Base Data Collection and Gully Sampling
KH-4B satellite images, Quickbird-2 satellite images and UAV images were used to acquire information about topography, land use and land cover in 1969, 2010 and 2019, respectively (Table 1). KH-4B (where KH is Key Hole), a space reconnaissance satellite launched during the cold war era, offered the best ground resolution of approximately 1.8 m [38]. During the operational phase, panchromatic images of large areas of the world were recorded on film. These satellite images were declassified in 1995 and have been available to the public since then and are used for environmental and resource research [38,39] (KH-4B images can be downloaded from the U.S. Geological Survey website, https://earthexplorer.usgs.gov/, accessed on 19 May 2020). The Quickbird-2 images contained four multispectral bands (MULs) with a 2.44 m spatial resolution and one panchromatic band (PAN) with a 0.61 m spatial resolution, which were fused to create a 0.61 m pansharpened image after the images were preprocessed (Table 1). UAV images were acquired in 2019 with a 0.14 m ground sampling distance (GSD) with a flying height of 500 m. Pix4Dmapper 4.4.12 software was used to automatically mosaic and process UAV images, to generate digital orthograph models (DOMs) and digital elevation models (DEMs) with a spatial resolution of 0.14 m. KH-4B and Quickbird-2 images were georeferenced to UAV images with a set of ground control points (GCPs) visible in all three image phases. Gullies along the tableland boundary were investigated by even sampling, but the selected gullies needed to be easy to identify in the 1969, 2010, and 2019 images. The gully edges were interpreted visually from the KH-4B, Quickbird-2, and UAV images, and linear and areal gully growth rates (R L , m·yr −1 and R A , m 2 ·yr −1 ) were calculated. The present drainage area (A, m 2 ) of each investigated gully head, slope gradient at the gully head (S g , m·m −1 ), and average slope of the drainage area (S d , m·m −1 ) were extracted based on the DEMs generated from the UAV images. The land use types, including cropland, apple orchards, forestland, shrubland, grassland, bare land, fallow fields, paved roads, unpaved roads and construction sites as well as terraces in the drainage area in 1969, 2010, and 2019 were interpreted visually based on those images. All map layers were georeferenced to the UTM Zone 49N coordinate system. Precipitation data from 1969 to 2019 were obtained from the CMDC (Figure 2). The precipitation records were divided into three periods (1969-1999, 2000-2014, and 2015-2019), representing records before the GFGP, after the GFGP, and after the GSTP. To explore the relationship between gully growth rates and precipitation characteristics, precipitation parameters were calculated, i.e., average annual precipitation, average annual erosive precipitation (>12 mm·d −1 ) [40], rainfall amount in storm events (>50 mm·d −1 ), and rainfall amount in heavy rainstorm events (>100 mm·d −1 ) [41].

Topographic and Hydraulic Parameter Calculation and Data Analysis
The stream power index (SPI) is a parameter related to surface and subsurface runoff transport processes [42] and has been widely used in gully erosion susceptibility research [43][44][45][46]. The SPI at the gully head was calculated as follows: where A is the present drainage area (m 2 ) and Sd is the average slope in the drainage area (m·m −1 ). The runoff onto one gully head (Qgh, m 3 ) from the drainage area with certain land use scenarios was calculated as follows: where Qi is the runoff depth (mm) of land use type i, Ai is the area (m 2 ) of land use type i in the drainage area of each gully head, and n is the number of land use types in the drainage area of each gully head. The NRC runoff curve number method (CN method, USDA-National Engineering Handbook, Part 630, Chapter 10, 13 January 2021, www.nrcs.usda.gov/wps/portal/nrcs/detailfull/national/water/manage/hydrology) was used to estimate the event runoff of land use type i onto the gully head from the drainage area (Qi), which takes the impact of land use and precipitation into account [5,47,48].
where P is rainfall (mm), Si is the maximum potential loss to runoff for land use type i (mm), λ is the fraction of S representing the initial abstraction coefficient (λ is 0.2 in this paper), and CNi is the curve number of land use type i. Key impact factors of the curve number (CN) include soil properties, land use type, slope, antecedent moisture conditions, vegetation coverage, and land management practices [49,50]. Different land uses

Topographic and Hydraulic Parameter Calculation and Data Analysis
The stream power index (SPI) is a parameter related to surface and subsurface runoff transport processes [42] and has been widely used in gully erosion susceptibility research [43][44][45][46]. The SPI at the gully head was calculated as follows: where A is the present drainage area (m 2 ) and S d is the average slope in the drainage area (m·m −1 ). The runoff onto one gully head (Q gh , m 3 ) from the drainage area with certain land use scenarios was calculated as follows: where Q i is the runoff depth (mm) of land use type i, A i is the area (m 2 ) of land use type i in the drainage area of each gully head, and n is the number of land use types in the drainage area of each gully head. The NRC runoff curve number method (CN method, USDA-National Engineering Handbook, Part 630, Chapter 10, www.nrcs.usda.gov/wps/portal/nrcs/detailfull/ national/water/manage/hydrology, accessed on 13 January 2021) was used to estimate the event runoff of land use type i onto the gully head from the drainage area (Q i ), which takes the impact of land use and precipitation into account [5,47,48].
where P is rainfall (mm), S i is the maximum potential loss to runoff for land use type i (mm), λ is the fraction of S representing the initial abstraction coefficient (λ is 0.2 in this paper), and CN i is the curve number of land use type i. Key impact factors of the curve number (CN) include soil properties, land use type, slope, antecedent moisture conditions, vegetation coverage, and land management practices [49,50]. Different land uses were assigned to different CN values (Table 2) based on previous research on similar land uses on the Loess Plateau [51]. Runoff in 3 periods (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), and 2015-2019) generated by annual average total precipitation (P a69 , P a10 , and P a19 ), erosive precipitation (P ae69 , P ae10 , and P ae19 ), accumulated rainfall in rainstorm events (P r69 , P r10 , and P r19 ), and accumulated rainfall in heavy rainstorm events (P hr69 , P hr10 , and P hr19 ) were calculated via Equation (3) and corresponding land use scenarios (1969, 2010, and 2019) as follows: where Q a (m 3 ·yr −1 ) and Q ae (m 3 ·yr −1 ) are the annual average runoff in 51 years from 1969 to 2019 caused by total precipitation and erosive precipitation, respectively; Q r (m 3 ) and Q hr (m 3 ) are the accumulated runoff generated by rainstorms (39 events) and heavy rainstorms (2 events), respectively, from 1969 to 2019. The same method can be used to calculate the annual average runoff (Q a6910 and Q ae6910 , m 3 ·yr −1 ) and accumulated runoff (Q r6910 , m 3 ) generated by rainstorms (38 events) from 1969 to 2010.

Statistical Analysis and Assessment Methods
Classical statistics, such as measures of center, variation and shape, were used to assess the change of gully dimensions. The Mann-Whitney test and Pearson chi-square test (χ 2 ) were used to compare the differences in A, S d , S g , and Q gh between the 2 types of gullies (growing and stable). Pearson correlation coefficients (r) were used to evaluate the relationships of A, S d , S g , and Q gh to the gully head growth rates. Nonlinear regression was used to fit models for the gully growth rates. The coefficient of determination (R 2 ) and significant level of F test were used to evaluate the performance of the regression models. SPSS 24 software was used for statistical calculation and model parameter calibration.

Gully Growth from 1969 to 2019
According to changes in the gully edges from 1969 to 2019, 95 investigated gullies were classified into four types: stable, growing, infilled, and excavated ( Figure 3). Twentyfive gullies were classified as stable and had no obvious changes between the 1969 and Remote Sens. 2021, 13, 5053 7 of 17 2019 images. Infilled and excavated gullies showed obvious gully area changes caused by human activities. Twenty-one gullies were infilled with a decrease in gully area ranging from 30.22 to 708.17 m 2 , with an average of 257.15 m 2 . These infilled gullies enlarged the tableland area and prevented further gully headcuts locally. Four gullies were excavated with increases in gully areas of 71.74, 185.41, 244.16, and 355.16 m 2 . Two of the four excavated gullies were excavated for cave dwellings, and another two gullies were cut on slopes to prevent landslides. The average present drainage area, local slope, and average annual linear/areal gully growth rate of the four types of gullies are listed in Table 3. Among the 95 investigated gullies, 51 gullies were located on the Zhongduo tableland, and 44 gullies were located on the Shiquan tableland ( Figure 4).   A-Average present drainage area; S g -average slope gradient near the gully head; R L -linear gully growth rate; R A -areal gully growth rate. 1 (Figure 5b).

Influencing Factors of and Modeling of Gully Growth Rates
Infilled and excavated gullies were controlled by human activity; therefore, they were not taken into consideration for further analysis of the factors that influence gully growth rates.
growing from 2010 to 2019. The RL values from 1969 to 2010 of the 45 growing gullies ranged from 0.04 to 2.78 m·yr −1 , with an average of 0.30 m·yr −1 , and the RL values of 34 gullies were less than 0.30 m·yr −1 . The highest RL value was 2.78 m·yr −1 , which was probably due to the unpaved road on the gully head as a concentrated channel for runoff leading to gully growth (Figure 5a). The RA values from 1969 to 2010 of the 45 growing gullies ranged from 0.50 to 20.94 m 2 ·yr −1 , with an average of 5.66 m 2 ·yr −1 . Fifty percent of the RA values ranged from 2.64 to 6.72 m 2 ·yr −1 (Figure 5b).  Gullies without obvious human intervention were analyzed to determine the influencing factors from the aspects of topography and hydraulics. The Mann-Whitney test indicated that the differences in four influencing factors (A, S g , S d , and Q gh ) between growing gullies and stable gullies were statistically significant (p < 0.05). The average A values of growing and stable gullies were 2.96 ha and 2.08 ha, respectively. The average S g values of the growing and stable gully heads were 0.07 m·m −1 and 0.05 m·m −1 , respectively. The average S d values of the growing and stable gully heads were 0.11 m·m −1 and 0.16 m·m −1 , respectively. The Q a values of the growing and stable gully heads were 12942.49 m 3 ·yr −1 and 8351.19 m 3 ·yr −1 , respectively. The Pearson chi-square test indicated that gullies were growing or stable and depended significantly on the presence/absence of terraces in 1969 (χ 2 = 14.731, p < 0.001). Among the 25 stable gullies from 1969 to 2019, 15 were drained areas with terraces since 1969, while among the 45 growing gullies, 7 were drained areas with terraces since 1969. This result showed that terraces, as an effective soil and water conservation measure, had strong effects on gully growth in the past 50 years in the tableland region.
Pearson's correlation tests were used to examine the relationship between 45 growing gullies and their influencing factors. No significant correlations were found between S d or S g and gully growth rates from 1969 to 2010, while A, Q a6910 , Q ae6910 , Q r6910 , and Q hr correlated significantly with the gully head growth rates from 1969 to 2010 (Table 4). Compared to the relationship between R L and its influencing factors, larger correlation coefficients were found between R A and its factors. The correlation coefficients between R A and Q a6910 , Q ae6910 , and Q hr were larger than those between R A and A. These results indicate that gully growth rates were significantly influenced by runoff.
Nonlinear regression models for the gully growth rate from 1969 to 2010 were fitted, as summarized in Table 5. Drainage areas were widely used as surrogates for surface or subsurface runoff in previous studies of gully erosion [52], but the low R 2 indicated that they did not work well in this study. Compounding the average annual runoff (Q a6910 ) and average slope of the drainage area (S d ) presented a better performance in modeling the gully growth rate, as shown in Equation (12) in Table 5. The measured R A of each gully and the corresponding R A simulated by Equation (12) is plotted in Figure 6, suggesting the good performance of Equation (12). A-present drainage area in m 2 ; S g -local slope gradients near the gully head in m·m −1 ; S d -slope gradients of the drainage area in m·m −1 ; Q a6910 and Q ae6910 -annual average runoff in 41 years from 1969 to 2010 calculated with total precipitation and erosive precipitation, respectively, in m 3 ·yr −1 ; Q r6910 and Q hr -accumulated runoff generated by rainstorms (38 events) and heavy rainstorms (2 events), respectively, in m 3 . SPI-stream power index. ** Correlation is significant at the 0.01 level. * Correlation is significant at the 0.05 level.   (1)

Threshold Conditions of Gully Head Growth
The topographic slope-area (Sg-A) relationship has been widely used to understand the threshold conditions for gully initiation and development [7,52] and is commonly expressed in critical Sg-A relationships, i.e., Sg = kA −b . The Sg-A relationship of stable gullies was fitted with the lowermost data points with Sg = 0.006 −0.289 , and the threshold condition for growing gullies was Sg = 0.024A −0.354 (Figure 7a

Threshold Conditions of Gully Head Growth
The topographic slope-area (S g -A) relationship has been widely used to understand the threshold conditions for gully initiation and development [7,52] and is commonly expressed in critical S g -A relationships, i.e., S g = kA −b . The S g -A relationship of stable gullies was fitted with the lowermost data points with S g = 0.006 −0.289 , and the threshold condition for growing gullies was S g = 0.024A −0.354 (Figure 7a On the other hand, Figure 7a shows that 17 stable gullies were above the line Sg = 0.024 A −0.354 , among which 10 gullies developed terraces in their drainage area since 1969, indicating that the terraces could change the threshold condition by reducing the effective drainage area. Then, the critical relationship Sg-Qa can reveal the threshold conditions for gully growth better than the Sg-Qa relationship because gully erosion is driven mainly by accumulated runoff, which is usually substituted with drainage area (Figure 7b

Long-Term Gully Growth Rates in the Tableland of the Loess Plateau
Gully growth rates in different studies are listed in Table 6. Compared with the hilly gully region in the medium-or long-term, the average RL (0.30 m·yr −1 ) in the tableland region in this study was lower. The first reason may be the much smaller slope gradients of the drainage area on the tableland than in the hilly gully region. When the hydraulic radius is constant, the shear stress of the shallow flow is proportional to the slope gradient [53]; therefore, gentler slopes may produce a smaller shear stress of the shallow flow. At the same time, the slope gradient is related to hydrological connectivity [35][36][37]. Hydrological connectivity can be quantitatively described by topographic parameters, including the slope gradient and upslope area [54], so that gentler slopes may result in lower hydrological connectivity. In addition, soil and water conservation measures are very common in the study area and have an effect on hydrological connectivity. The tableland boundary was ridged with a shrub belt, and small reservoirs in the drainage area could reduce direct runoff into gullies. Terraces can also lead to a further decrease in hydrological connectivity [36,55]. The hydrological connectivity of the drainage area on the tableland should be lower than that in the hilly gully region owing to the gentler slope, soil and water conservation measures, and piping erosion in the drainage area. Lower hydrological connectivity may be the main reason for the lower gully growth rates in the tableland region. The On the other hand, Figure 7a shows that 17 stable gullies were above the line S g = 0.024 A −0.354 , among which 10 gullies developed terraces in their drainage area since 1969, indicating that the terraces could change the threshold condition by reducing the effective drainage area. Then, the critical relationship S g -Q a can reveal the threshold conditions for gully growth better than the S g -Q a relationship because gully erosion is driven mainly by accumulated runoff, which is usually substituted with drainage area (Figure 7b). Critical conditions for stable and growing gullies were S g = 0.483Q a −0.359 and S g = 0.070Q a −0.286 , respectively. The minimum annual average runoff onto the investigated growing gullies was approximately 1000 m 3 ·yr −1 , compared with 80 m 3 ·yr −1 for the investigated stable gullies. Among 45 growing gullies, 17 gullies that still grew from 2010 to 2019 are marked as gray dots in Figure 7.

Long-Term Gully Growth Rates in the Tableland of the Loess Plateau
Gully growth rates in different studies are listed in Table 6. Compared with the hilly gully region in the medium-or long-term, the average R L (0.30 m·yr −1 ) in the tableland region in this study was lower. The first reason may be the much smaller slope gradients of the drainage area on the tableland than in the hilly gully region. When the hydraulic radius is constant, the shear stress of the shallow flow is proportional to the slope gradient [53]; therefore, gentler slopes may produce a smaller shear stress of the shallow flow. At the same time, the slope gradient is related to hydrological connectivity [35][36][37]. Hydrological connectivity can be quantitatively described by topographic parameters, including the slope gradient and upslope area [54], so that gentler slopes may result in lower hydrological connectivity. In addition, soil and water conservation measures are very common in the study area and have an effect on hydrological connectivity. The tableland boundary was ridged with a shrub belt, and small reservoirs in the drainage area could reduce direct runoff into gullies. Terraces can also lead to a further decrease in hydrological connectivity [36,55]. The hydrological connectivity of the drainage area on the tableland should be lower than that in the hilly gully region owing to the gentler slope, soil and water conservation measures, and piping erosion in the drainage area. Lower hydrological connectivity may be the main reason for the lower gully growth rates in the tableland region. The average R L (0.30 m·yr −1 ) in this study was also lower (Table 6) than the gully growth rate in the Moldavian Plateau of Romania, where the topographical characters are similar to the hilly gully region of the Loess Plateau and the higher R L could be owed to the massive deforestation in 1990. Although the gully growth rates in the tableland regions were relatively low, unpaved roads with high erodibility, as the key components of hydrological connectivity [36], could lead to gully growth [56,57]. An unpaved road to the gully head may be the main reason for the largest erosion rate. Thus, soil and water conservation measures should focus on unpaved roads in tableland regions in the future. R L -Linear gully growth rate; R A -areal gully growth rate; P a -annual precipitation; MP-monitoring period.
Runoff is the main driver of gully headcuts because gully heads can only develop if concentrated (overland) flow intensity during a rain event exceeds a threshold value [1]. However, measuring critical flow hydraulics (e.g., runoff discharge and velocity) for gully erosion is constrained under field conditions [1,3,5]. The CN method provides an approach to estimate direct runoff in a certain land use type [47][48][49]. Likewise, combining a pixelbased CN approach with flow-routing algorithms makes it possible to account for the effect of the spatial patterns of topography, soil conditions and land cover [5]. Resent research has also found that runoff estimated by the CN method could be a good predictor for gully erosion at the rainstorm event scale [57]. In this study, the annual average runoff estimated by the CN method (Q a ) was used to fit a nonlinear regression model for the long-term areal gully growth rate (R A =0.301Q a 0.562 S d , R 2 = 0.530), suggesting that the long-term gully growth rate could be estimated with runoff and the slope of the drainage area as predictors. Soil and water conservation measures affect gully erosion by reducing runoff or slope gradients, which could be taken into account with the model presented in this study. However, it is difficult to accurately estimate runoff from the drainage area in the tableland of the Loess Plateau because soil and water conservation measures, including small reservoirs, could reduce runoff into the gullies, and there is sometimes piping erosion in the drainage areas. This can explain why the gully growth rate model built in this study with runoff as a predictor did not work as well as the models in previous studies [56,57]. This study mainly focused on the effects of topography, land uses and conservation practices with a reasonable assumption that soil properties varies slightly on the loess covered area [60]. The main soil type in the study area is the loessial soil that has a composition of more than 50% silt (0.002-0.05 mm) and less than 20% clay (<0.002 mm), and a porosity of approximately 50% [61], leading to very high infiltration capacity [61]. Nevertheless, further researches should be carried to investigate the variation of the soil properties and the infiltration rates caused by land use changes.

Topographic Thresholds of Gully Growth
The topographic threshold is the most common approach to understanding the critical condition for gully initiation and development [5,7,52]. Based on 45 growing gully heads and 25 stable gully heads, S-A relationships were fitted in this study, following the methods in previous research [7,52]. Parameters k and b showed very large variabilities in different environments [7], and variations in k seem mainly attributable to differences in land cover under the assumption that b is relatively constant [5,7]. The low k value indicated that the land cover condition in the study area was more prone to erosion [7]. For exponent b, the results of this study (b = 0.289) are comparable to those of several previous studies (Table 7). The S g -Q a relationship can express the threshold conditions better than the S g -A relationship [57] because the drainage area is only a substitute for runoff [52]. However, runoff from the drainage area was not estimated accurately in the tableland of the Loess Plateau (Section 4.1) so that the S g -Q a relationship is not much better than the S g -A relationship in this study.

Potentiality of Gully Growth
Gully erosion occurs when threshold conditions are exceeded and is controlled by the integrated effects of rainfall, flow hydraulics, topography, soil type, and land use [62,[65][66][67]. In other words, a gully can be stable if the threshold is not exceeded. Topographic thresholds can partly explain the development of gullies. Figure 7a shows that growing gullies are far above the critical line, S g = 0.006A −0.289 , indicating that, in terms of topographic conditions, these growing gullies have the potential to continue growing and gullies farther away from the threshold line have greater potential to grow. However, 25 gullies remained stable during 1969-2019, among which 15 gullies had terraces in their drainage area according to the images of 1969 and 2019, and 24 gullies had terraces in their drainage area according to the images of 2019. Terraces can change the effective drainage area by reducing runoff from the drainage area and then influencing gully erosion.
Since the GFGP and GSTP were put into practice in the tableland region of the Loess Plateau, soil and water conservation measures, including small reservoirs, terraces, shrub belts, and water barriers, have been widely adopted (Figure 8). Among 45 growing gullies from 1969 to 2010, only 17 continued to grow from 2010 to 2019, suggesting that the GFGP and GSTP reduced gully erosion to some extent. Gully erosion in the study area has been intervened with greatly by human activity. Some gullies were infilled and excavated by human activities directly, which are easily recognized. However, in most cases, gully development is influenced indirectly by soil and water conservation measures, land use changes, roads constructions, water conservancy facilities, etc. Quantifying of these effects need to be integrated with reliable runoff simulation model, which should be taken into consideration for future studies. On the other hand, studies have reported that several single rainfall events with high rainfall intensities promote the formation and development of gullies [56,65,[68][69][70]. With the gully growth potential in terms of topographic conditions, whether soil and water conservation measures are effective for gully erosion under extremely high rainfall intensity needs further investigation. However, the results of this study are still helpful for preventing further gully erosion by dealing with gullies far above the threshold line.
consideration for future studies. On the other hand, studies have reported that several single rainfall events with high rainfall intensities promote the formation and development of gullies [56,65,[68][69][70]. With the gully growth potential in terms of topographic conditions, whether soil and water conservation measures are effective for gully erosion under extremely high rainfall intensity needs further investigation. However, the results of this study are still helpful for preventing further gully erosion by dealing with gullies far above the threshold line.

Conclusions
In this study, high-resolution images were used to investigate gully growth rates from 1969 to 2019. Forty-five of 95 investigated gullies grew from 1969 to 2010 at rates ranging from 0.50 to 20.94 m 2 ·yr −1 (0.04 to 2.78 m·yr −1 ), with an average of 5.66 m 2 ·yr −1 (0.30 m·yr −1 ). Twenty-five gullies were stable, 21 gullies were infilled, and four were excavated by human intervention. The present drainage area, local slope, average drainage slope, annual runoff, and ratio of terraced area were all significantly different between the stable

Conclusions
In this study, high-resolution images were used to investigate gully growth rates from 1969 to 2019. Forty-five of 95 investigated gullies grew from 1969 to 2010 at rates ranging from 0.50 to 20.94 m 2 ·yr −1 (0.04 to 2.78 m·yr −1 ), with an average of 5.66 m 2 ·yr −1 (0.30 m·yr −1 ). Twenty-five gullies were stable, 21 gullies were infilled, and four were excavated by human intervention. The present drainage area, local slope, average drainage slope, annual runoff, and ratio of terraced area were all significantly different between the stable and growing gullies. For growing gullies, runoff is an important factor that influences the gully growth rates. The long-term areal gully growth rate (R A ) could be estimated using a nonlinear regression model with the annual average runoff (Q a ) and average slope of the drainage area (S d ) as predictors (R A = 0.301Q a 0.562 S d , R 2 = 0.530). Based on 45 growing gullies and 25 stable gullies, the S g -Q a relationship was able to reveal the threshold conditions for gully growth. Based on the threshold conditions, all growing gullies still have the potential to keep growing, while soil and water conservation measures could change the threshold condition by reducing the effective drainage area. An increasing number of gullies would possibly remain stable in the future since GFGP and GSTP were launched. Many soil and water conservation measures (including terraces, small reservoirs, shrub belts, etc.) and change of cropland to apple orchards have strong effects on gully erosion mitigation. The results of this study could be helpful for preventing further gully erosion by dealing with gullies far above the threshold line.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.