Yushu Zhou , Bin Zhang ,*, Wei Qin , Qingchun Deng , Jun Luo , Hui Liu ,Dan Yang , Houling Wang , Ying Zhao

a School of Land and Resources, China West Normal University, Nanchong, 637009, PR China

b State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing,100048, PR China

c Research Center on Soil and Water Conservation of the Ministry of Water Resources, Beijing,100048, PR China

d Sichuan Provincial Engineering Laboratory of Monitoring and Control for Soil Erosion on Dry Valleys, China West Normal University, Nanchong, 637009,PR China

Keywords:Gully distribution Environmental factors Kernel density estimation Correlation analysis Large scale Northeastern China

ABSTRACT Gully erosion is the main cause of global land degradation. The factors controlling gully erosion at watershed scale have been extensively studied, but the spatial pattern of gullies and their primary environmental factors are unclear. In order to explore the primary factors at the regional scale of over thousands of square kilometers, 34 sample areas were selected in northeastern China and 10 environmental factors were proposed. Based on gully interpretations from high-resolution Google Earth images and spatial analysis, kernel density estimation (with polyline features), and correlation analysis, the results show that the maximum kernel density(KD)of the 34 sample areas ranged from 0.11 to 87.16,and there were significant positive correlations between the average KD and the maximum value and standard deviation.At the regional scale, topography(relief amplitude, slope,surface roughness, surface incision,and elevation with correlation coefficients(ρ)of 0.60,0.59,0.57,0.48,and 0.41 respectively)was the primary environmental factor, with a significant positive correlation with KD. The correlation coefficients of the normalized difference vegetation index(NDVI)(ρ=0.32),annual precipitation(ρ=0.19),and soil (silt, sand, and clay content with ρ of 0.139, -0.117, and 0.085, respectively) decreased successively with the KD. At sample area scale, there were six areas with weak correlations between environmental factors and KD,20 areas with topographic factors,six areas with annual precipitation,and two areas with NDVI as the primary factor.The correlation between soil texture and KD was low.There were different positive or negative correlations between the environmental factors and KD. The spatial distribution of gullies in northeastern China is heterogeneous,and the primary environmental factors were diverse. This study will help to understand the spatial pattern and formation mechanism of gullies at large scale.

1. Introduction

Soil erosion is one of the most important global environmental issues and it has been proposed as the main cause of land degradation (Valentin et al., 2005). Gully erosion is a common type of hydraulic erosion, and it is “a channel or miniature valley cut by concentrated runoff through which water commonly flows only during and immediately after heavy rains or during the melting of snow” (Ankeny,1982). Gullies can be either ephemeral or permanent channels (Amare et al., 2019; Poesen et al., 2003), which can provide morphological evidence of past intense soil erosion and reflect local environmental changes (Ionita et al., 2015; Luo et al.,2019).

Gully initiation is a critical geomorphological process and is controlled by multiple factors (Poesen et al., 2003), including terrain morphology, soil, rainfall, vegetation, runoff, lithology, and human activities(Valentin et al.,2005).Geographical differences in environmental conditions can lead to heterogeneity in the distribution of gullies(Sonneveld et al.,2005).Under the influence of an irregular boundary,runoff will rapidly cause gully erosion(Zhang&Zhong,1998).The development of a gully requires sufficient runoff,while the horizontal concentration of runoff is limited,and only the increase of slope length can achieve the required conditions(Wu&Liu, 2000). Slope is also important for gully formation, in mountainous or hilly areas, high-speed runoff occurs on steep slopes,making them prone to gully erosion (Valentin et al., 2005), with extensive erosion patterns developing on the larger slopes (Prelat et al., 2015).

The influence of slope on gully density varies with region, and gully density is positively correlated with slope in some areas(Jiang et al., 2017; Yan et al., 2005a), but in other regions it may initially increase and then decrease with increasing slope length (Xu et al.,2014). Slope direction is an important factor affecting gully distribution, due to a series of complex factors such as freezing and thawing, snow melt differences and rain direction; the gully density on sunny slopes is slightly larger than on shady slopes (Yan et al., 2005b). Soil mechanical composition, organic matter content, and land vegetation conditions are also important factors affecting gully erosion(Poesen et al.,2003;Xia et al.,2019).The clay content of soil is important in determining its susceptibility to gully erosion(Amiri et al.,2019).The depletion of soil organic matter will reduce the stability of the soil structure, and is conducive to crusting, runoff production, and erosion (Valentin et al., 2004).Vegetation can intercept rainfall and so reduce damage to the soil structure and the erosion capacity of runoff (Pathak et al., 2005),thereby delaying gully development(Li et al.,2015).In addition,the plant root system enhances the ability of soil to resist runoff erosion by improving soil structural stability and infiltration (Li, 1995).Precipitation and runoff generated by precipitation have an impact on gully erosion; under the same terrain and rainfall duration, an increase in the intensity of rainfall will generate more runoff,often leading to intensified gully erosion(Yao,2014).Human activity also affects the distribution of gullies. In recent decades, the black soil area of northeast China has experienced high levels of human activity, including land reclamation, and inappropriate farming systems and land uses, which have weakened the anti-erodibility of soil and rapidly increased gully density (Fan et al., 2005; Wang et al., 2014).

Geographical differences play a leading role in gully formation in some locations, but may not be significant in other areas(Sonneveld et al., 2005). At the regional scale, the spatial differentiation of gully density in the Loess Plateau is mainly affected by geological structure, rainfall, vegetation, and ground materials(Tian et al., 2013), however, the gully density in Heilongjiang Province generally increases with precipitation(Fan et al.,2015).At the watershed scale,soil properties,vegetation,and terrain are the main factors that contribute to gully erosion in small watersheds in the dry-hot valleys of China (Deng et al., 2015a, 2015b; Ding et al.,2017; Xu et al., 2018), while extreme rainstorms in the Majiagou watershed on the Loess Plateau have a large effect on gully formation (Jia et al., 2014).

The conditions that influence gully development and distribution are complex,and the relationship between gully development and environmental conditions at different watershed scales in northeast China has been investigated, but the spatial distribution pattern of gullies and their primary factors at the regional scale are still unclear. The objectives of this study were to determine the spatial distribution law of gullies and to identify the primary factors affecting gully distribution on the regional scale in northeast China.The results will not only help to reveal the regional gully distribution law and formation mechanism,but also provide a scientific basis for controlling gullies according to local conditions.

2. Material and methods

2.1. Description of study area

The study area(40°36′～53°13′N,118°12′～132°15′E)is located in northeastern China,and includes parts of Heilongjiang Province,Jilin Province, Liaoning Province, and eastern Inner Mongolia Autonomous Region (Fig. 1). The study area lies in the midtemperate and cold temperate zones and has a temperate monsoon climate with four distinct seasons, warm and rainy summers,and cold and dry winters. From southeast to northwest,the annual precipitation drops from 1000 mm to less than 300 mm.Surrounded by the mountains and hills of the Great Khingan Range,Lesser Khingan Range,and the Changbai mountains,the Northeast Plains (including Song-nen Plain, Liao River Plain, and Sanjiang Plain) lie in the middle of the study area. Forests (Larix gmelinii,Pinus koraiensis,andBetula costata)are the dominant land cover in the mountains,with temperate forest steppe,meadow steppe,and steppe distributed in the plain areas. Several drought-tolerant plants, such asAneurolepidium chinenseandStipa baicalensis, are typical of the desertified land that extends to the southwest(Wang et al.,2009).The average vegetation coverage is above 70%,with the highest coverage in the east(92%)and lowest coverage in the west(6%).The soil types in the study area are complex,with dark brown soil mainly distributed in mountainous areas, brown soil in the Liaodong Hilly area, and black soil and chernozem in the plain areas. As the main grain production area in China, the area has experienced serious gully erosion due to the environmental conditions and human activities.

2.2. Selection of sample areas

It would be difficult to interpret all of the gullies in the whole area due to the vast size of the territory. The selection of appropriate sample areas enabled an effective alternative method to be developed. Belt transects were established running either parallel or perpendicular to the mountain ranges,and were coincident with the prevailing wind direction in the area. Five belt transects were established in the northwest-southeast direction and seven in the northeast-southwest direction. The intersection point of the two transects was the central location of the sample area. The intersection point was used to make a square, with a side length of 30 km and trending in the north-south or east-west directions.A total of 35 samples were obtained.The numbering of sample areas was based on theSijformat, whereiis the northwest-southeast transect number (1-5), andjis the northeast-southwest transect number(1-7).Because sample area S51 lay outside of China,it was excluded from the study and only the other 34 sample areas were considered.

2.3. Data sources

10 factors of 4 types, including topography (slope, relief amplitude, surface roughness, and surface incision and elevation),vegetation(average normalized difference vegetation index,NDVI),annual precipitation,and soil texture(silt, sand,and clay content),were selected to characterize natural conditions affecting gullies.NDVI data(July 2018),vegetation zoning data of China(July 2018),annual precipitation data (2015), and soil texture of the topsoil layer at a depth of 0-20 cm in this study were obtained from the Resources and Environmental Data Cloud Platform of the Chinese Academy of Sciences (http://www.resdc.cn/). The NDVI data uses China’s monthly vegetation index spatial distribution data set products,the data on China’s vegetation zoning are from vegetation maps at a scale of 1:1000000 and the precipitation data was derived from China’s annual spatial precipitation interpolation data set since 1980. The terrain data was obtained from the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM),with a resolution of 90 m,which was originally produced by NASA(http://srtm.csi.cgiar.org/). All the above data are raster files with cell size of 1000 m. In this study, the GCS_WGS_1984 geographic coordinate system was used to describe the spatial data.

Fig.1. The sample areas distributed in the study area.

2.4. Gully interpretation

Google Earth images attained in the last five years with spatial resolution of over 2.0 m were used in this study. http://www.bigemap.com/. At the scale of 1:4000, gullies were visually interpreted in ArcGIS software (version 10.3) according to the shape,texture, tone and size. The gully shapefile was created with the polyline feature type and the geographic coordinate system of GCS_WGS_1984 in ArcCatalog,and then added to the ArcMap.After the image was also added, the gully thalweg can be visually vectorized following the principle from the main gully to branches,and from gully head to mouth (Fig. 2). The “length” field was added to the attribute table of the gully shapefile, and then calculated with the Calculate Geometry tool in ArcGIS.Gullies with a length of less than 200 m were removed,35,463 gullies were identified and used.This hard work cost our team of 32 people more than a month. It must be pointed out that a large gully can be interpreted whether it is a permanent gully or ephemeral one. Based on the images, it is difficult to distinguish an ephemeral gully from a permanent one.In order to improve the accuracy of the gully interpretation, a crosscheck was made after interpretation.Based on high-resolution Google Earth images, an accurate gully data set was obtained by visual interpretation in ArcGIS (Wang et al., 2016; Yibeltal et al.,2019).

2.5. Topographic factor extraction

Topographic factors were very important to gully formation because they can control the velocity and volume of the concentrated surface flow (Garosi et al., 2018). The DEM reflected the distribution of elevation across the study area.Based on the spatial analysis toolbox of the ArcGIS software,the parameters of slope(S),relief amplitude (Ra), surface roughness (Sr), and surface incision(Si) for each sample area could be obtained.

(1) Relief amplitude (Ra)

Rais an important parameter for quantitatively describing geomorphological features and their evolution (Yin et al., 2010),and can reflect the height difference in a certain area(Cheng et al.,2016; Zhang & Dong, 2012).Rais the difference between the maximum (Hmax) and the minimum (Hmin) elevation within a given area(statistical window):

The mean change-point method was then used to calculate the optimal statistical window(Han et al., 2012).

(2)Surface incision (Si)

Sican be used to quantitatively characterize the intensity of surface erosion and incision. It is the difference between the average (Hmean) and minimum (Hmin) elevation in a given area(Hu & Nie, 2015; Liu & Tang, 2006):

Fig. 2. Gullies in the Google Earth image. (left: gullies in hillslope; right: gullies in mountains).

(3)Surface roughness (Sr)Sr

is a parameter used to reflect the degree of surface fragmentation and erosion, and terrain features at different scales(Grohmann et al., 2011; Hobson,1972). It refers to the ratio of the surface area (As) to the plane area (Ap) in a given area:

2.6. Methods

2.6.1. Kernel density evaluation (KDE)

KDEis recognized as an accurate and effective method for estimating gully distribution(Ran,2018).It is an effective technique for the measurement of the probability density function of a random variable (Lopez-Novoa et al., 2016), and considered to be a fundamental data smoothing process in statistics (Lopez-Novoa et al.,2016).KDEcan measure local density changes and explore the continuous distribution of events, which can avoid the issue of prior knowledge or human subjectivity.KDEmainly uses a moving cell (or a window) to calculate a magnitude-per-unit area from point or polyline features using a kernel function to fit a smoothly tapered surface to each point or polyline. The Rosenblatt-Parzen formula is used to estimate theKD(Parzen, 1962; Rosenblatt,1956; Ruppert & Cline,1994):

wherek(x) is the kernel function andhis the bandwidth. (x-xi) is the distance from the estimated pointxto the samplexi,andnis the number of points whose path distance from the locationxiis less than or equal toh. In ArcGIS, a default bandwidth is automatically calculated based on the spatial configuration and number of input points.The algorithm is as follows:1)to calculate the mean center of the points, the distance from the mean center of all points, the median of these distances(Dm),and the standard distance(SD);2)to calculate the bandwidth:

The kernel density with polyline features calculated in ArcGIS is output as a raster file with cell size of 1000 m. Fig. 3 shows the kernel density map of a sample area.

We can make the minimum and maximum value ofKDof a raster, which characterizes the relative sparse,dense gullies in the area. The averageKDcan show the overall level of gully distribution.The standard deviation(SD)and coefficient of variation(CV)ofKDindicate the heterogeneity of gully distribution in an area.

2.6.2. Correlation analysis

As the Pearson correlation coefficient is applicable to paired,continuous, and normally distributed data, the data used in this study do not meet the above conditions and the Spearman correlation coefficient is preferred.The Spearman correlation coefficient was used to measure the correlation between environmental factors andKD:

wherenis the sample size and d is the ranking difference set ρ∈[-1,1], in which the larger the |ρ| is, the stronger the correlation is.According to the correlation coefficient, the strength of the correlation could be:|ρ|≥0.6,strong correlation;0.4≤|ρ|＜0.6,moderate correlation; 0.2≤|ρ| ＜0.4, low correlation; and |ρ| ＜0.2, weak correlation.A total of 30600 vector points(representing a pixel)of 34 sample areas extracted by function of fishnet points in ArcGIS are used as the data of correlation analysis. The attributes of vector points include the value ofKDand every environmental factor.The correlation of regional scale is calculated in SPSS software(version 17.0) with 30600 points (all sample areas), and the correlation of sample scale is calculated with 900 points (each of sample areas).

3. Results

3.1. Non-uniformity of gully distribution

The spatial distribution of gullies in northeastern China was found to be extremely unbalanced (Table 1). The maximumKDvaried substantially among the sample areas. There were three sample areas with maximum values less than 1.0, accounting for 8.8%of the total number of sample areas,which indicated that there were very few gullies in this area. The maximumKDof 14 sample areas was in the range of 2.0-10.0,accounting for 41.2%of the total.17 sample areas had a maximumKDgreater than 10.0, accounting for 50.0%of the total area.Three sample areas had a maximumKDgreater than 50.0, with sample area S47 having a maximumKDof 83.96.As to averageKD(Fig.4),there were four sample areas with a value of less than 0.1, accounting for 11.8% of the total, indicating few gullies in these areas;there were 18 areas with an averageKDin the range of 0.1-1.0,accounting for 52.9%of the total;there were six, five sample areas with an averageKDin the range of 1.0-2.0 and five with 2.0-3.0;there was only one area with an averageKDof over 3.0,accounting for 2.94%.The skewness of the averageKDis 1.13, which shows a positive skew distribution that there were more samples with smaller averageKD.The kurtosis of the average kernel density is 0.40, showing a fat tails distribution. The 95%confidence interval around the averageKDwas 0.63-1.30, the values of averageKDof six sample areas (S11, S12, S16, S27, S43,S56) fell within the confidence interval; but those of the other 28 sample areas went beyond the interval, which shows an obvious difference in gully density and distribution. For the standard deviation ofKD,there were nine sample areas with a value of less than 0.50,accounting for 26.5%of the total area.This indicated that the gully distribution in these areas was relatively uniform. The standard deviations of nine sample areas were in the range of 0.5-1.0,accounting for 26.5%. The standard deviation of 16 sample areas was greater than 1.0,accounting for 47.1%,which indicated that the gullies in these sample areas were very unevenly distributed.

Fig. 3. Kernel density distribution of gullies of a sample area (S16).

While thep-value between the standard deviation and coefficient of variation was 0.32,the tests between the other parameters were very significant(p＜0.01).The correlation coefficients among the maximum, average, and standard deviation ofKDwere 0.81,0.90, and 0.97; their relationship can be expressed by linear function with R2 of 0.65, 0.81, and 0.94, respectively (Fig. 5a-c). This highly positive correlation indicates that with an increase in the maximumKDof the sample area,spatial differences and clustering among gullies also increased significantly. The correlation coefficient between the coefficient of variation and averageKDwas-0.23,showing a weak trend that the spatial difference is less as there is more clustering of gully distribution (Fig. 5- d). The correlation coefficients between the maximum,standard deviation,and coefficient of variation were less than 0.1, showing no correlation between them (Fig. 5- e, f).

3.2. Correlations at the regional scale

The Spearman correlation coefficient was used to measure the correlation between environmental factors and theKDof all sample areas (Fig. 6). Among the natural environmental factors that affected the spatial distribution of gullies in northeastern China,topographical conditions were ranked first, followed in order by vegetation conditions, precipitation, and soil properties. The correlation coefficient between the five topographic factors andKDwas higher than that of the other environmental factors. For the topographical factors, the ranking of the correlation coefficients was listed as:Ra(ρ = 0.60) ＞S(ρ = 0.59) ＞Sr(ρ = 0.57) ＞Si(ρ = 0.48) ＞elevation (ρ = 0.41), which indicated that the relief amplitude, steepness, and surface roughness had the most significant impacts on gully distribution.At the regional scale,there was a positive correlation between topographic factors andKD, with strong positive correlations betweenRa,S, andKD, and moderate positive correlations betweenSr,Si, and elevation andKD. The correlation coefficient between the NDVI andKDwas 0.32, i.e., a low positive correlation. The correlation coefficient between precipitation andKDwas 0.19, also indicating a weak positive correlation.Soil texture had little influence on the distribution of gullies,with the absolute values of the correlation coefficients of the clay,sand,and silt content andKDall being less than 0.2.The correlation between the silt content andKD(ρ = 0.14) was relatively high,whereas the clay content(ρ=0.09)had almost no relationship withKD. There was a weak negative correlation between the sand content andKD(ρ=-0.12).At the regional scale,large-size gullies have undergone a long evolutionary stage, in which process and topographic conditions were the most fundamental and important factors, while the contributions of vegetation, precipitation, and soil conditions were relatively weak.

Table 1 Statistics of kernel density (KD)of the sample areas.

Fig. 4. Histogram of average KD.

3.3. Correlations at the local scale

At the local (or sample area) scale, the absolute values of the correlation coefficients between environmental factors and theKDof the sample areas(S25,S52,S41,S34,S24,and S14)did not exceed 0.20(Table 2).Sample areas S25,S34,S24,and S14 were located in the plain. The gullies were extremely underdeveloped, with averageKDvalues of less than 0.07; however, S52 and S41 had higher averageKDof over 1.5, but the correlation with environmental factors was very low. There were very few gullies in the sample areas in the plain, and they were not closely related to natural conditions, but were rather greatly affected by human activities. The correlations between the environmental factors and theKDin the remaining 28 sample areas indicated that there was a large difference in the gully distribution among the different areas.

3.3.1. Terrain

There were 20 sample areas in which the correlation coefficients between topographic factors andKDwere ranked highest,indicating that topographic factors were still the primary factors(Table 3).There were nine sample areas(S17,S57,S14,S16,S22,S46,S13, S32, and S56) with high correlation coefficients betweenKDand elevation,accounting for 32.14%of the 28 sample areas.For the relationship betweenKDand elevation in the 28 sample areas,the absolute correlation coefficient of S56, S13, S15, and S16 was greater than 0.60, for seven sample areas (S57, S44, S43, S37, S45,S12,and S42)it was in the range of 0.4-0.6,for eight sample areas(S46, S22, S17, S32, S27, S26, S35, and S36) it was in the range of 0.2-0.4, for nine sample areas (S23, S31, S32, S33, S36, S44, S46,S53, S55) it was less than 0.2, and for five sample areas (S53, S47,S21,S33,S23)it was less than 0.10.Except for five sample areas(S37,S35, S55, S11, and S33), there were positive correlations between KD and elevation.

Fig. 5. Relationships between statistical parameters of kernel density.

Fig. 6. Spearman correlation coefficient between environmental factors and kernel density (KD) (p ＜0.01).

Table 2 Correlation coefficients between environmental factors and kernel density (KD).

Table 3 The number of sample areas in different correlation coefficient categories.

The sample areas of S45, S12, S55, S23, S21, S53, and S11 had significant correlation coefficients between averageKDandRa.The correlation coefficient for S56 was 0.64;S45,S15,S13,S44,S57,S12,S16,S42,and S55 had correlation coefficients between 0.4 and 0.6;and S57, S33, and S47 had correlation coefficients of less than 0.2.The correlation coefficients of the other 15 sample areas were between 0.2 and 0.4. The averageKDandRaof S37 and S35 were negatively correlated,while the others were positively correlated.

There were four sample areas (S54, S31, S36, S44) with the largest correlation coefficient between averageKDandSi.For nine sample areas, the correlation coefficients were in the range of 0.4-0.6,while 16 sample areas had values of 0.2-0.4.The values in the other three sample areas were all less than 0.1. Except for S37 and S35, the correlation coefficients were positive, indicating that there was a positive relationship betweenKDandSi.

AlthoughSandSrwere not the primary factors, they were still important environmental factors influencing gully distribution.ForS, the correlation coefficient of 10 sample areas was greater than 0.40, and the maximum was 0.64. ForSr, there were nine sample areas with correlation coefficients greater than 0.40. The sample areas with correlation coefficients less than 0.2 for these two factors accounted for less than 15% of the total area.

Of the 28 sample areas, there were nine in which topographic factors ranked in the top five factors controlling gully distribution(S56, S45, S57, S13, S44, S15, S12, S46, and S36). There were 12 sample areas (S42, S55, S23, S43, S21, S22, S54, S53, S11, S17, S31,and S26) with four topographic factors among the top five factors,only 1 sample area(S33) without topographic factors in the first 5 factors, 3 sample areas (S16, S35, S47) with 2 topographic factors,and 3 sample areas (S27, S32, S37) with 3 topographic factors.Therefore,topographic factors were considered to play a key role in gully distribution.

3.3.2. Precipitation

With regard to the relationship between precipitation and theKDof gullies, the correlation coefficients of six sample areas (S42,S43, S37, S26, S27, and S47) were ranked first. The correlation coefficients of S42 (ρ = 0.67) and S43 (ρ = 0.61) exceeded 0.6. Four sample areas (S37, S56, S16, and S57) had values between 0.4 and 0.6, all four sample areas are located in the southeast of the study area and have more precipitation than other areas.Thirteen sample areas (S26, S27, S45, S13, S12, S22, S11, S15, S35, S17, S47, S21, and S54) had absolute correlation coefficients between 0.2 and 0.4,while the remaining nine areas (S23, S31, S32, S33, S36, S44, S46,S53, S55) had values less than 0.2, of which the correlation coefficients of S23, S36, and S44 were less than 0.1. The correlation coefficients of eight sample areas(S23,S33,S55,S21,S35,S11,S26,and S37) were negative, with the value of S37 being -0.58, i.e., a strong negative correlation. TheKDof the other sample areas tended to increase with an increase in precipitation.

3.3.3. Vegetation

Considering the correlation coefficients between NDVI andKD,the gully distribution of two sample areas (S35 and S33 with low values of 0.28 and 0.21 respectively) was dominated by the NDVI,which belongs to the grassland meadow subzone in the northern temperate zone.The maximum correlation coefficient betweenKDand NDVI was 0.59 at S16, which is located between mixed coniferous and deciduous broadleaf in the south of the temperate zone and deciduous oak forest in the north of the temperate zone.There were eight sample areas(S56,S37,S27,S35,S55,S42,S11,and S33)where the absolute correlation coefficient was in the range of 0.2-0.4. There were 19 sample areas (S17, S46, S15, S44, S21, S36,S57,S54,S13,S22,S45,S23,S53,S32,S12,S26,S31,S43,S47)where the correlation coefficients were less than 0.2, with 11 (S13, S15,S21, S22,S23, S36, S44,S45, S46, S53, S57) being less than 0.1.TheKDof six sample areas (S15, S46, S12, S31, S11, and S37) was negatively correlated with the NDVI,while in the other areas there was a positive correlation. Overall, the correlation between NDVI and gullyKDwas not strong.

3.3.4. Soil conditions

The correlations between soil texture andKDwere not significant,and the maximum absolute value of the correlation coefficient was only 0.38 (S56). The number of sample areas with correlation coefficients for the relationship between the sand, silt, and clay content andKDin the range of 0.2-0.4 was seven (S17, S15, S56,S16, S11S, S55, S32), two (S13, S56), and seven (S55, S21, S32, S16,S11,S42,S56),respectively,while at other sites the values were all less than 0.2. The correlation coefficient for the relationship between the silt content andKDwas distributed in the range of -0.26-0.18, with the values of 16 samples being negative. The correlation coefficient for the relationship between clay andKDwas in the range of-0.38 and 0.26,with a negative correlation in 50%of the sample areas. The correlation coefficient for the relationship between the sand content andKDwas in the range of-0.26 to 0.36,of which the correlation coefficients of 13 sample areas were positive. Soil texture made a lesser contribution to gully distribution than the other environmental factors investigated in this study.

4. Discussion

4.1. Variability of gully distribution

At a large scale, the spatial pattern of gully distribution tended to be clustered. In northeastern China, theKDin a sample area typically differed in spatial distribution, and also the averageKDand standard deviation ofKDvaried significantly among all sample areas, showing that the spatial distribution of gullies is nonuniform or clustered (Ran, 2018). In addition, the gully density in the Middle Volga Region of Russia(with a total area of 188 000 km2)and in the Loess Plateau of China(with an area of 640 000 km2)had significant variations(Golosov et al.,2018;Tian et al.,2013).In the Middle Volga Region(Russia),gully density is mainly correlated to topographical factors and land use,with most gullies closely linked to agricultural activity (Golosov et al., 2018). However, the distribution trend of gully was controlled by geological structure in the Loess Plateau(Tian et al.,2013). In northeastern China, there were characteristics of a plain in the middle(Northeast China Plain)and three mountains surrounding the plain.By comparison,averageKDand relief have the same trend in spatial distribution; the correlation coefficient between topographic factor and kernel density was greater than 0.4, which is higher than other factors, so the terrain may control the gully distribution in the study area. According to the power relation for channel initiation between slope (S) at the gully head and upstream drainage area (A) (Patton & Schumm,1975):

wherea,bare coefficients, it showed that there is a significant negative trend betweenSandA(Vanwalleghem et al., 2005): a gully can be triggered where the slope is bigger with relative low drainage area,therefore, there is a higher gully density in the hills and mountains than in the plains in northeastern China.

4.2. Scale difference of controlling factors

There were different environmental factors controlling gully development at different scales.At the runoff plot or hillslope scale,differences in site conditions(such as slope and aspect)may lead to different rates of gully development (Deng et al., 2020; Ran et al.,2018). At the small watershed scale, in a typical karst small watershed in the Lijiang River(SW China),elevation (300-500 m)and slope [15-25°] were reported to have the greatest impact on soil erosion (Qin et al., 2018). At the medium and large watershed scale, in the Wuyur River watershed (NE China), the dominant factors affecting gully development were found to be precipitation,topography, and soil (Hu & Wang, 2018); however, in the Wuding River Basin in northern Shaanxi(China),the main factors affecting gully density were geological factors, dynamic conditions, and ground surface properties (Yang & Liu, 2014).

In northeastern China, elevation,S, slope direction, and slope length were reported to have the largest impact on gully development in Liaoning Province (Huang & Fan, 2017), while gully density in Jilin Province increases with the increase of slope,and a strong linear relationship was identified between them (Xu et al.,2017). In a hilly area, slope has been shown to be an important factor controlling gully erosion on hillslope farmland and has a strong linear correlation with gully density. This contrasts with gully formation in the Loess Plateau, which is strongly correlated with vegetation, rainfall intensity, and surface material composition (Tian et al., 2013). In this study, gully distribution in northeastern China was primarily related to topographic factors at the regional scale, while the correlation between precipitation, soil factor,and kernel density was less than 0.2,meaning that there was less effect. At the sample scale, topographic factors were not the major causes of gully development in some sample areas. For example,in S42 and S43,the correlation between precipitation and kernel density is greater than 0.6, which is greater than that of topographic factors.In the sample area with high precipitation,the correlation between precipitation and kernel density is also relatively high. This shows that the primary factors controlling gully development in different regions and at different scales may be quite different.

Fig. 7. The active gully destroying the farmland in Heilongjiang Province (NE China).

4.3. Combined effect of human activities

Human activities mainly affect gully development from the aspects of vegetation destruction, agricultural land uses, and construction work. Plant roots reduce gully erosion in improving soil physical properties such as structural stability and infiltrability(Valentin et al., 2005). The original forests such as Great Khingan Range,Lesser Khingan Range,and the Changbai mountains,formed obvious high vegetation coverage zones.NDVI below medium level is mainly distributed in the arid and semi-arid regions in the west and southwest of northeast China. The vegetation degradation in these regions is also relatively obvious, which is closely related to people’s weak awareness of vegetation protection, overgrazing,deforestation, and other acts of vegetation destruction (Hu et al.,2017). Due to deforestation and grassland reclamation, natural erosion can be accelerated to a rate hundreds or even thousands of times the natural erosion rate under the same conditions(Fan et al.,2013). In recent decades, the intensity of human activities has increased and the vegetation destruction has been accelerated in northeastern China, which has weakened soil anti-erodibility and promoted gully development (Deng & Wang, 2019; Dong et al.,2019; Yu & Wu, 2003). With the rapid urbanization, the increased urban population has brought about a change of land use pattern,the destruction of vegetation ecological areas, and the significant decline of vegetation coverage (Xiao et al., 2017). The buffer area close to human settlements is subject to serious damage, due to destruction of the soil surface and vegetation.Therefore,the closer areas are to settlements,the higher the gully density will be(Wang et al., 2014). Road construction in mountainous areas will also damage the stability of the soil, and the erosion rate at a construction site is 2-40 000 times higher than before the land was disturbed(Harbor,1999).The gullies in this study area were formed both over previous geologic periods and due to modern human activities; with the former being relatively stable, while the latter were mostly active (Fig. 7). The distribution of active gullies and their formation mechanism should be the focus of future research,which will be especially important for gully control.

5. Conclusions

Based on kernel density evaluation of 34 sample areas, the relationship between the spatial distribution of gullies and natural environmental factors in northeastern China was investigated.The results showed that:(1)the gully distribution was very different in the region with averageKDof 0.01-3.51,maximumKDof 0.11-83.9,SD of 0-17.38, and CV of 0-11.0; (2) at the regional scale, the topographic factor had a higher correlation coefficient withKDthan that of other environmental factors, the relief amplitude was the primary factor controlling gully distribution (ρ = 0.60); (3) at the scale of the sample area, the gully distribution in the plain has a weak correlation with environmental factors.There were 58.82%of sample areas with topographic variables as primary factors (9 sample areas with elevation,7 with relief amplitude,4 with surface incision). Precipitation was ranked first in six sample areas, accounting for 18.75% of the total sample areas. Vegetation is the most relevant to gully distribution in two sample areas. However,soil contents cannot be the primary factor of gully distribution;(4)as to different areas, there were either positive or negative correlations between environmental factors andKD.In conclusion, the primary factors controlling the heterogeneity of gully distribution at the scale of the sample areas were more complex than those at the regional scale, with topographic factors still playing a key role in gully distribution.

Declaration of competing interest

The authors declared that they have no conflicts of interest to this work.

Acknowledgements

This research was funded by the IWHR Research&Development Support Program (Grant No. SE0145B132017), the National Key Research and Development Program of China (Grant No.2018YFC0507002),and the National Natural Science Foundation of China (Grant No. 41971015, 41807075), the Project of Science &Technology Department of Sichuan Province (Grant No.2018SZ0337, 2017JY0189), the Project of Sichuan Provincial Department of Education(Grant No.18TD0025,17AZ0385),and the Meritocracy Research Funds of China West Normal University(Grant No.17YC134,17YC105),the Fundamental Research Funds of China West Normal University (Grant No.17C032,16A001).