Integrated Use of Gcm, Rs, and Gis for the Assessment of Hillslope and Gully Erosion in the Mushi River Sub-catchment, Northeast China

The black soil region of Northeast China has suffered from severe soil erosion by water. Hillslope and gully erosion are the main erosion types. The objective of this research was to integrate the assessment of hillslope and gully erosion and explore spatial coupling relations between them in the Mushi River sub-catchment using geographical conditions monitoring (GCM) including remote sensing (RS) and geographic information system (GIS) techniques. The revised universal soil loss equation (RUSLE) model and visual satellite image interpretation were used to evaluate hillslope and gully erosion, respectively. The results showed that (1) the study area as a whole had slight erosion due to rill and sheet erosion, but suffered more serious gully erosion, which mainly occurs in cultivated land; (2) GCM contributed to the overall improvement of soil erosion assessment, but the RUSLE model likely overestimates the erosion rate in dry land; (3) the hillslope and gully erosion were stronger on sunny slopes than on shady slopes, and mainly occurred at middle elevations. When the slope was greater than 15 degrees, the slope was not the main factor restricting the erosion, while at steeper slopes, the dominant forest land significantly reduced the soil loss; (4) trends of gully erosion intensity and density were not consistent with the change in soil erosion intensity. To our knowledge, this study was one of the first that attempted to integrate gully erosion and hillslope erosion on a watershed scale. The findings of this study promote a better understanding of the spatial coupling relationships between hillslope and gully erosion and similarly indicate that GCM, RS, and GIS can be used efficiently in the hilly black soil region of Northeast China to assess hillslope and gully erosion.


Introduction
Soil erosion is regarded as a major environmental problem for sustainable development of human society because it seriously threatens natural resources, agriculture, and the environment [1].It is also one of the most significant forms of land degradation, causing 80% of the world's land degradation [2].Soil erosion has resulted in a series of environmental concerns, such as undermining land resource quality, leading to significant environmental disasters, and influencing sustainable regional development [3].Relevant and major global research programs and international organizations regard soil erosion as an important research topic [4].For example, the Future Earth Sustainability 2016, 8, 317; doi:10.3390/su8040317www.mdpi.com/journal/sustainability(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022)(2023) program focuses on studying the shared challenges of particular places and regions such as soil degradation [5].
The black soil region in Northeast China is considered one of the typical areas of the world that have been characterized by a high-tension man-land relationship within short-term time scales [6] and suffered severe soil erosion due to the mass reclamation in the last 100 years [7].Soil loss caused by water erosion including hillslope erosion and gully erosion in this region is mainly derived from farmlands in sloping areas [8].Previous studies on soil erosion in the black soil region were carried out from different perspectives [9][10][11][12][13][14][15].The conclusions confirmed that hillslope and gully erosion closely interact with each other both in space and time.However, few studies focus on exploring the relationship between hillslope and gully erosion on a large scale [16].Meanwhile, limitations of timeliness, accuracy, and precision of the information about the natural environment and human interference have impeded soil erosion assessment in black soil regions.
To obtain an objective and comprehensive geographical view of the spatial distribution and variation of natural resources, the National Administration of Surveying, Mapping, and Geo-information (NASG) started the Geographical Conditions Monitoring (GCM) project in 2012.The geographical conditions information census is the first national census for China, making full use of high-resolution remote sensing images including space-borne sensors (Pleiades-1A, Worldview 2, etc.), airborne sensors, and ground data.Important results obtained from the census, such as topography, land cover, and other information, will provide a basis for a more sustainable environment [17] and also provide a possibility for soil erosion assessment in black soil regions.
Splash erosion, sheet erosion (the previous two types are also called interrill erosion), rill erosion, and gully erosion are generally recognized as the main types of soil erosion by water [18].Splash erosion marks the beginning of other types of soil erosion involving raindrop impact, splash of soil particles, and formation of craters [19].It is called sheet erosion when a thin layer of soil is removed by raindrop impact and shallow surface flow from the whole slope.When rainfall exceeds the rate of infiltration, water accumulates on the surface of gently sloped land and moving water concentrates along tiny channels called rills that are less than 30 cm deep [18].Gully erosion begins when runoff concentrates in channels and rills develop, which may grow into gullies deeper than 30 cm over short time periods [20].
Spatial and quantitative information on soil erosion on a watershed scale contributes significantly to soil conservation management, erosion control, and catchment environment management [21].Previous studies showed that assessment of soil erosion methods mainly included empirical models [22,23], physically based models [15,[24][25][26], nuclear tracing [27], and spatially distributed multivariate models [28,29].Among these models, the revised universal soil loss equation (RUSLE) [23] is a revised version of the universal soil loss equation (USLE) model [22].RUSLE proved to be the most commonly applied model to estimate soil loss [30][31][32].Recently, many studies have applied remote sensing (RS) and Geographic Information System (GIS) methodologies to assess soil erosion [29,31,[33][34][35].In particular, satellite images with very high spatial resolution have offered possibilities for landscape feature mapping and quantitative evaluation of farm practices in soil erosion assessment [36,37].
RUSLE is designed to model sheet (interrill) and rill erosion and evaluate soil loss from agricultural areas caused by overland flow [38][39][40].However, it does not assess gully erosion [41].Gully erosion causes more significant soil loss and soil degradation compared with hillslope erosion [42].There are numerous techniques for monitoring and quantification of gully erosion, such as ruler, tape, microtopographic profilers [43], poles, total stations [44], pins [45], differential GPS [11,46], terrestrial laser scanners [47], stereoscopic photogrammetry on ground [48], aerial photographs [49,50], visual analysis of satellite images [51,52], object-oriented analysis [53], 3D photo-reconstruction [44,54], and unmanned aerial vehicles [55].Analysis of satellite images combined with field validation has proven to be an excellent and practical approach for mapping of gully features and assessing gully erosion over large areas [10,56].However, the spatial resolution of satellite images in previous studies was mostly higher than 5 m.According to a field survey, a considerable proportion of erosion gullies have widths less than 5 m, especially ephemeral gullies in the black soil region, and the analysis results are insufficient to estimate the potential hazard of gully erosion [57].Higher temporal and spatial resolution images provided new possibilities for studying gully erosion [53,58,59].
This study aims to integrate the assessment of hillslope and gully erosion.RUSLE was used to calculate soil loss caused by hillslope erosion, and visual satellite image interpretation with field verification was used to map geometric gully features and evaluate gully erosion.Another objective is to elucidate the spatial coupling relationship between hillslope and gully erosion using spatial overlay statistics based on the previous assessment.

Study Area
The study area shown in Figure 1e mainly included the middle and lower area of Mushi River basin with serious soil and water loss covering 253.13 km 2 .The target area (126 ˝1'39"E-126 ˝17'14"E, 44 ˝9'18"N-44 ˝22'20"N) is located in the Jiutai District (Figure 1d), Jilin Province, Northeast China (Figure 1a,b).This area is situated at the southeastern edge of the black soil region (Figure 1b).As the dominant soil type, black soil is classified as Udic Isohumisol in the Chinese Soil Taxonomy or Luvic Phaeozem in the FAO (Food and Agriculture Organization of the United Nations) system [60].The parent materials are Quaternary lacustrine and fluvial sand beds or loess sediments [61], and the main textural classes of the topsoil are silt-clay-loam to clay-loam (8%-27% sand, 29%-66% silt and 26%-40% clay) [62].Other soil types include meadow soil, albic soil, dark brown soil, and new alluvial soil.The catchment area has a continental climate with four distinctive seasons.The annual temperature is 5.3 ˝C, the highest average monthly temperature is 23.3 ˝C, and the lowest average monthly temperature is ´16.3 ˝C.The annual precipitation ranges from 550 to 600 mm.Sloped farmland is the main land use type because this area is located in the hilly area between Changbai Mountain and Songliao plain with large terrain undulation and little vegetation cover.Therefore, gully erosion is widely distributed in this area, which is one of the typical black soil erosion regions of Northeast China.

Data Sources
The sources and products of GCM, including high-resolution satellite images, digital elevation model (DEM), land cover data, and other important geographical condition data (road, water, residential sites, etc.), were provided by the Jilin Province Geomatics Center.
The high-resolution Pleiades image was used for the mapping of farming practices such as contour tillage and the monitoring of gully erosion.Pleiades-1A imagery shows a 0.7-m vertical viewing panchromatic band and four 2.8-m multi-spectral bands (blue, green, red, and near-infrared) (Figure 1e).The images were collected on 12 October 2013, and have been processed into DOM by the data providers.
The DEM of 5-m pixel size was produced using 1:10,000 scale topographical maps for terrain elevation and other topography-related analyses.The land use/land cover (LULC) map derived from Pleiades-1A, aerial images, and ground data were used as the basis for C-and P-factor determination.The classification accuracy of the LULC map from GCM was notably higher (Figure 2) compared with the land-use map with a scale of 1:100,000 (derived from Landsat8 imagery that was acquired on 16 October 2013, and by the Chinese Academy of Sciences).The soil spatial distribution map was derived from the 1:1,000,000 scale nationwide soil database, which was provided by the Chinese Academy of Sciences.Detailed soil attributes including mechanical composition, grain sizes, and organic content were taken from the second national soil survey in China.
where A is the mean soil loss per year (t¨hm ´2¨a ´1), R is the rainfall erosivity factor (MJ¨mm¨hm ´2¨h ´1¨a ´1), K is the soil erodibility factor (t¨hm ´2¨h ¨MJ ´1¨mm ´1¨hm ´2), LS is the topographical factor, C is the cover and management factor, and P is the conservation practices factor; L, S, C, and P are dimensionless.

Rainfall Erosivity (R)
Rainfall is the main driving force of soil erosion.Rainfall erosivity represents the potential of soil erosion caused by raindrop splash runoff erosion [63].The product of rainfall kinetic energy (E) and the maximum 30 min rainfall intensity (I30) is a classical calculation method for R.However, it is very difficult to reflect continuous records of rainfall data with this method.Therefore, many scholars have proposed more mature and simpler calculation methods of rainfall erosivity.In this study, we used a method based on monthly and annual rainfall proposed by Wischmeier [22] because of its simple equation, wide application range, and the ease of obtaining data.The equation model is shown as follows: R " where i is the month, P i is the monthly rainfall (mm), and P is the annual precipitation (mm).The R value of the formula is given in U.S. unit 100 ft¨tonf¨in/(acre¨h¨year); hence, it should be multiplied by a coefficient of 17.2 for conversion into the international unit MJ mm/(hm 2 ¨h¨a) [64].
The R factors of meteorological stations around the study area were calculated by using the monthly precipitation data from 1958 to 2014.Spatial distribution of the R factor in the study area was interpolated using the ordinary Kriging method.

Soil Erodibility Factor (K)
The soil erodibility factor (K) was used to measure the sensitivity of soil erosion to water and indicate the resistance of soil to erosion.In this study, the EPIC (erosion-productivity impact calculator) model was selected to estimate the soil erodibility factor [65].The formula is as follows: where S a is the sand (0.05-2 mm) content (%); S i is the silt (0.002-0.05 mm) content (%); C l is the clay content (%); C is the organic carbon content (%); and S n = (1 ´Sa )/100.The K value of the formula is given in U.S. unit t¨acre¨h/(100 acre¨ft¨tonf¨in).Thus, it should be multiplied by a coefficient of 0.1317 for conversion into the international unit t¨hm 2 ¨h/(MJ¨mm¨hm 2 ) [64].
Note that the statistical standard of soil particle size in the equation is based on the U.S. system, but the second national soil survey of China uses the international system.According to the graphic method [66], the sizes were converted by statistical software.The conversion equation correlation coefficients R 2 are all above 0.95 (p < 0.05).

Topographic Factor (LS)
The topographic factor (LS) indicates the impact of slope length and slope steepness on soil erosion.LS was also considered to be a crucial factor for surface runoff and soil erosion [23].Slope steepness has a significant effect on the rates of soil erosion.The steeper the slope is, the faster overland flow runs.Thus, the steeper slope will then increase the shear force of the soil particles.In addition, surface runoff and flow rate will have a steady growth with increasing slope length, leading to a greater erosion force on the soil surface [67].Usually the product of length and slope of terrain is used as the terrain factor.This study estimated the LS factor using Equation 4 [67].Usually the product of length and slope of terrain was used as terrain factors.This paper estimated the LS factor by using the formula [9]: where L is the slope length factor; S is the slope steepness factor; λ is the field slope length (m); and α is the slope angle in percentage.This paper selected the procedure derived from Hickey [68] and Van Remortel [69], and modified the corresponding parameter to generate the LS factor raster.

Cover Management Factor (C)
The C factor indicates the effect of vegetation cover and management on soil erosion.It is an important factor to evaluate the ability of vegetation canopy and ground covers to resist soil erosion, and its value is between 0 and 1.Using remote sensing image interpretation to acquire land cover type maps [70] or collecting land use data [67,71] directly and then assigning C values from the existing literature to the corresponding land use/land cover types is the main means of accessing C factors on the river basin and regional scale [72].This paper assigned C values to different LULC classes derived from GCM (Table 1) according to existing research in the black soil region [9,13,16].The practice factor (P) is defined as the ratio of soil loss after a specific preservation practice to the corresponding soil loss after upward and downward cultivation under otherwise equal conditions, ranging from 0 to 1, in which 0 represents the areas with no soil erosion due to good protective measures and 1 reflects the areas with no support practice.At the watershed level, farming practices such as terrace and contour ridge cannot be reflected in a land use map [73] but are clear in QuickBird images [37].
Ridge tillage is the main farming practice in the black soil region.This paper used the angle between ridge direction and slope aspect to calculate the P value, which is based on previous studies [37].The P equation is given as follows (Equation 5): where P is the value of the Water Conservation Measures factor and α is the angle between ridge direction and slope aspect in degrees (0 ˝-90 ˝).P h is 0.352 when the angle is 90 ˝ [9], i.e., Contour Ridge, which means that the ridge direction is parallel to the contour line and perpendicular to the slope aspect.P s is 1 when the angle is 0 ˝, i.e., Downslope Ridge, which means that the ridge direction is perpendicular to the contour line and parallel to the slope aspect.The specific steps are as follows: (1) according to the visual interpretation of the Pleiades image, this paper segmented the farmland into plots with the same ridge and assigned the ridge direction azimuth to each plot by outlining the ridge direction line (Figure 3); (2) DEM was used to generate the slope aspect layer; (3) the ArcGIS grid calculator was applied to calculate the angle between the ridge and slope direction.p-values of arable land on the flat ground and other LULC types were assigned as 1, usually without any soil and water conservation measures.

Calculation of Soil Erosion
According to the RUSLE equation, the R, K, LS, C, and P data layers (Figure 4) were considered to generate the hillslope erosion map in ArcGIS (Figure 5).The calculated soil loss values were divided into six erosion classes according to the soil erosion rate standard, Technological Standard of Soil Erosion Grade and classification SL190-96, issued by the Ministry of Water Resources of China [32].

Mapping and Assessment of Gully Erosion
In this research, visual 0.7-m resolution Pleiades-1A imagery interpretation was used to obtain gully distribution data of the study area.The gullies were classified as either ephemeral or permanent gully [42].Stream-channels and gullies in forests were not objects of this research.With the help of expert knowledge and field surveys that were conducted in May 2012, gully interpretation signs on Pleiades images were determined (Table 2).The gullies were mainly discriminated according to their spectral and spatial characteristics.Based on the topographical map, hydrographical net map, traffic map, and related information, the data of gully distribution were mapped including a polyline layer and a polygon layer representing gully length and area, respectively.ArcGIS 10.0 software was combined with expert experience.Gullies with a width less than 3 pixels were only outlined for length, not area.Field validation work was conducted in May 2015, based on which the gully data were modified and the final gully distribution data were obtained (Figure 6a).According to the field validation, the precision of gully interpretation was 88.9%.Some field paths are easily mistaken for ephemeral gullies.

Calculation of Soil Erosion
According to the RUSLE equation, the R, K, LS, C, and P data layers (Figure 4) were considered to generate the hillslope erosion map in ArcGIS (Figure 5).The calculated soil loss values were divided into six erosion classes according to the soil erosion rate standard, Technological Standard of Soil Erosion Grade and classification SL190-96, issued by the Ministry of Water Resources of China [32].

Mapping and Assessment of Gully Erosion
In this research, visual 0.7-m resolution Pleiades-1A imagery interpretation was used to obtain gully distribution data of the study area.The gullies were classified as either ephemeral or permanent gully [42].Stream-channels and gullies in forests were not objects of this research.With the help of expert knowledge and field surveys that were conducted in May 2012, gully interpretation signs on Pleiades images were determined (Table 2).The gullies were mainly discriminated according to their spectral and spatial characteristics.Based on the topographical map, hydrographical net map, traffic map, and related information, the data of gully distribution were mapped including a polyline layer and a polygon layer representing gully length and area, respectively.ArcGIS 10.0 software was combined with expert experience.Gullies with a width less than 3 pixels were only outlined for length, not area.Field validation work was conducted in May 2015, based on which the gully data were modified and the final gully distribution data were obtained (Figure 6a).According to the field validation, the precision of gully interpretation was 88.9%.Some field paths are easily mistaken for ephemeral gullies.As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area.The distribution of gully density was estimated by using the line density tool in ArcGIS (Figure 6b).The density was calculated as the ratio of total length of the gully within the circular kernel (50-m search radius) and the total area of the circular kernel.Similarly, gully intensity was calculated by using the measurement of total gully area per area of total subcatchment.The gully intensity was used as the main indicator in the analysis of topographical differentiation of hillslope and gully erosion.As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area.The distribution of gully density was estimated by using the line density tool in ArcGIS (Figure 6b).The density was calculated as the ratio of total length of the gully within the circular kernel (50-m search radius) and the total area of the circular kernel.Similarly, gully intensity was calculated by using the measurement of total gully area per area of total subcatchment.The gully intensity was used as the main indicator in the analysis of topographical differentiation of hillslope and gully erosion.As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area.The distribution of gully density was estimated by using the line density tool in ArcGIS (Figure 6b).The density was calculated as the ratio of total length of the gully within the circular kernel (50-m search radius) and the total area of the circular kernel.Similarly, gully intensity was calculated by using the measurement of total gully area per area of total subcatchment.The gully intensity was used as the main indicator in the analysis of topographical differentiation of hillslope and gully erosion.As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area.The distribution of gully density was estimated by using the line density tool in ArcGIS (Figure 6b).The density was calculated as the ratio of total length of the gully within the circular kernel (50-m search radius) and the total area of the circular kernel.Similarly, gully intensity was calculated by using the measurement of total gully area per area of total subcatchment.The gully intensity was used as the main indicator in the analysis of topographical differentiation of hillslope and gully erosion.As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area.The distribution of gully density was estimated by using the line density tool in ArcGIS (Figure 6b).The density was calculated as the ratio of total length of the gully within the circular kernel (50-m search radius) and the total area of the circular kernel.Similarly, gully intensity was calculated by using the measurement of total gully area per area of total subcatchment.The gully intensity was used as the main indicator in the analysis of topographical differentiation of hillslope and gully erosion.As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area.The distribution of gully density was estimated by using the line density tool in ArcGIS (Figure 6b).The density was calculated as the ratio of total length of the gully within the circular kernel (50-m search radius) and the total area of the circular kernel.Similarly, gully intensity was calculated by using the measurement of total gully area per area of total subcatchment.The gully intensity was used as the main indicator in the analysis of topographical differentiation of hillslope and gully erosion.
As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area.The distribution of gully density was estimated by using the line density tool in ArcGIS (Figure 6b).The density was calculated as the ratio of total length of the gully within the circular kernel (50-m search radius) and the total area of the circular kernel.Similarly, gully intensity was calculated by using the measurement of total gully area per area of total subcatchment.The gully intensity was used as the main indicator in the analysis of topographical differentiation of hillslope and gully erosion.

Soil Loss Due to Hillslope Erosion Using RUSLE
The rainfall erosivity (R) factor values range from 3568.45 to 3717.97 MJ¨mm¨hm -2 ¨h-1 ¨a-1 and increase from northwest to southeast (Figure 4a).This result indicates that the R factor was affected by the distribution of meteorological stations outside the area when missing sites in the study area.The soil erodibility (K) factor values were obtained using Equation 3 depending on different soil types and range from 0.029 to 0.038 t¨hm 2 ¨h¨MJ -1 ¨mm -1 ¨hm -2 (Figure 4b).Dark brown soil and black soil have higher erodibility values (K was 0.038 t¨hm 2 ¨h¨MJ -1 ¨mm -1 ¨hm -2 ), and meadow soil is less erodible (0.029 t¨hm 2 ¨h¨MJ -1 ¨mm -1 ¨hm -2 ) than other soil types.The topographic (LS) factor values range from 0 (flat areas, mainly in the valley plain between the hills) to 58.25 (steep hills, in the northern and southern part of the study area) (Figure 4c).As the main types of LULC in the study area, dry land occupies 65% of the total study area, while 26% of the area is forest.The C values range from 0.005 to 0.07.Desert and bare ground have the highest C values (1), while the C factor was 0 for water areas, housing construction areas, structures, and roads (Figure 4d).
The slope of the 59% slope farmland is higher than 5 degrees, where soil erosion is more serious due to unfavorable terrain factors.Concerning the angle between the ridge and slope direction, the P factor values range from 0.352 to 1 (Figure 4e).With the increase in the angle, the proportion of slope farmland is gradually reduced, and more than 75% of the whole slope farmland is cultivated land with askew ridges (angle between 10 ˝and 80 ˝).Remarkably, askew ridges are extremely common.
The hillslope soil erosion map of the study area is shown in Figure 5.The mean annual soil erosion due to rill and sheet erosion is approximately 18.07 t ha -1 a -1 , and the average soil loss amounts to 45.60 ˆ10 4 t in 2013 in the study area.The statistical results of soil loss (Table 3) show that 4.52% of the watershed area belongs to the severe or violent erosion class.Intense erosion occurs in 5.57% of the study area, while 12.42% of land is affected by moderate erosion.Overall, tiny and slight erosion dominated the region by covering 33.77% and 43.71% of the study area, respectively.However, with respect to the total amount of soil loss, only 1.42% of the total soil erosion occurred in tiny erosion areas, while slight erosion contributed to 23.56% of total soil loss.The moderate erosion area shows the largest amount of soil loss.In the study area, 50.07%and 10.09% of the total soil loss was observed under intense to violent erosion, respectively.

Estimation of Gully Erosion
The status of gully erosion in the study area is shown in Table 4.There are 1674 erosion gullies in the study area, among which 859 (more than half of the total) were erosion gullies with a width less than 3 pixels.This result indicated that the gully still relatively active in the study area.We assumed that these narrow gullies would further develop and become intensified erosion without any protection.In 2013, the gully density of the whole area was 1.07 km/km 2 (gully length per unit area), and the total area and length were 1,409,219 m 2 and 271.08 km, respectively.Figure 6b shows that the gully density ranged from 0 to 48.70 km/km 2 , and denser gully distribution was found in the northeastern part, where gully erosion was most serious.In addition, the erosion gullies, located in the cultivated land, occupied 79% of the total length, and the gully density was 1.65 km/km 2 for the whole farmland.Places within a certain distance from the erosion gully were not suitable for farming, which would accelerate the erosion [14].Therefore, the cultivated land areas engulfed, destroyed, and affected by gully erosion were much larger than 60.78 hm 2 , which is the gully area in farmland.These results indicate that the study area suffered more serious gully erosion, which mainly occurs in cultivated land.

Implementation of the RUSLE Model and Visual Image Interpretation
RUSLE has enhanced soil loss predication capabilities and can be applied to natural environmental conditions [23], compared with the USLE model, which was developed initially for field crops in agricultural areas [22].RUSLE estimates long-term average annual soil erosion amounts in a variety of environments, such as agriculture, forest, rangeland, mining sites, and construction sites [74].The USLE/RUSLE model is considered adaptable for estimation of soil erosion by water in the black soil region of Northeast China, where long slope length and gentle slope are the main topographical features [9,13,75,76].However, RUSLE may provide extremely high soil erosion predictions [31].
To validate the RUSLE model, this study selected the existing observations from previous research to examine the results.Zhang et al. (1992) monitored soil loss of different plots for many years in the experimental areas of Keshan and Binxian, which are also located in the black soil region [9].This study obtained grids consistent with the conditions of four plots in the experimental areas using raster calculator and zonal statistics tools in ArcGIS.The average calculated values of soil loss of these grids were compared with the observed values in the plots (Table 5).The calculated values were significantly higher than the observed ones, with an average deviation of 72.85% in the last three plots, while the calculation was close to the observation in the first plot.The results suggest that the RUSLE model likely overestimates the erosion rate in dry land even though the rainfall erosivity (R) factor values in the Mushi River subcatchment were significantly higher than those in Keshan and Binxian [77].The seven meteorological stations were distributed around the study area evenly (Figure 1c).With the existence of spatial dependence, the rainfall erosivity (R) factor can be estimated for the non-sampled locations through spatial interpolation.As one of the most effective interpolation methods, the ordinary Kriging method used the point data of these stations to produce a continuous surface covering the study area.However, it was difficult to obtain a very high resolution and accurate spatial distribution map of the rainfall erosivity (R) factor with raster format by spatial interpolation because the meteorological stations were scarce within the basin [78].Integration of multi-source rainfall data to improve the accuracy of spatial interpolation needs further consideration in the future.In addition, the very poor resolution of the soil map constrained the results significantly.A more precise soil map would improve the accuracy of the model calculation results in the future.Nevertheless, the high accuracy of LULS maps and high-resolution images contribute to the overall improvement of the erosion model [29].Compared with other methods for calculating the P factor in the black soil region of Northeast China [9,13,79], the calculation accuracy of the P factor for ridge tillage in this study was improved.
Pleiades imagery showed a good feasibility for the recognition of ephemeral gully in this study.The results indicate that sub-meter-resolution remote sensing images have great potential for gully mapping [57].However, this study only focuses on the length and area of the gully, ignoring gully volume, which is also an important factor for the determination of sediment yield from gully erosion.Application of remote sensing to extract the gully volume should be given further consideration.With respect to high spatial resolution satellite images (Pleiades-1A), although the visual interpretation method with field verification obtained satisfactory results in mapping gullies, (semi-)automatic mapping of gullies at a regional scale will be a future research direction [53,80,81].

Topographic Differentiation of Hillslope and Gully Erosion
The differentiation of hillslope and gully erosion on various elevations is demonstrated in Figure 7a.In the vertical direction, the hillslope and gully erosion showed a layered distribution, increasing first and then decreasing as the altitude increases.This result indicates that the geomorphic types of lower elevations (e.g., valley plain) were dominated by modern deposition.No significant erosion occurs at higher elevations (e.g., hill) due to the small catchment area, while slope land between the abovementioned regions was the main level of the occurrence and distribution of soil erosion [66].Moreover, the magnitude of gully erosion intensity (at 225-250 m) was lower than that of hillslope erosion intensity (at 250-275 m).We assume that, compared to slope erosion, gully erosion is not only affected by the slope but also by the catchment area and other topographic factors.The lower slope land is more conducive to the formation and development of gully erosion due to a longer slope and larger catchment area [10,56].
Sustainability 2016, 8, 317 14 of 20 soil loss compared with farmland.Although with increasing slope gravity-caused soil erosion increased and very steep slopes did not hold much soil [31], the dominant position of forest land gradually increased such that the proportion even exceeded 80% at 20-25 degrees (Figure 8), and forest land reduced the soil loss.When the slope was greater than 15 degrees, the gully intensity decreased significantly as the slope increased.Although the steeper slope will increase the shear force of the soil particles [67], the slope is not the main factor restricting erosion; it is more affected by slope length, catchment area, slope shape, and other terrain factors [16] when the slope is greater than 15 degrees.The relation of hillslope and gully erosion depending on different slope aspects was analyzed (Figure 7c).The results indicate that slope aspect has a consistent influence on gully erosion intensity and hillslope erosion.The hillslope and gully erosion on a sunny slope (S, SE, SW) were stronger than those on a shady slope (N, NE, NW) considering the monsoons, total solar radiation, diurnal temperature difference, and freezing and thawing [10].

The Relationship between Gully and Hillslope Erosion
The relationship between hillslope and gully erosion is illustrated in Figure 9. Gully erosion density and intensity represent gully length and area per unit area, respectively.The variation of gully erosion intensity and density is not consistent with soil erosion intensity.Gully intensity increases with slope erosion intensity, in particular rapidly increasing in the severe-violent erosion area.Gully erosion density reaches a peak at moderate erosion and then begins to decrease.This finding indicates that the control measures of hillslope erosion can effectively prevent the formation and development of erosion gully in the slight-moderate erosion area [16].However, the control measures need to be combined with other prevention and cure measures for gully erosion in moderate-violent erosion regions, where the lateral erosion of the gully erosion is significant, which results in notable widening of ditch sides and an increase in the area.The relation of hillslope and gully erosion on different slope is shown in Figure 7b.When the slope was less than 10-15 degrees, the hillslope erosion and gully erosion intensity both increased significantly with increasing slope.However, the hillslope erosion slightly decreased in steeper slopes at 20-25 degrees with an inverse trend.Forest land plays a predominant role in the control of soil loss compared with farmland.Although with increasing slope gravity-caused soil erosion increased and very steep slopes did not hold much soil [31], the dominant position of forest land gradually increased such that the proportion even exceeded 80% at 20-25 degrees (Figure 8), and forest land reduced the soil loss.When the slope was greater than 15 degrees, the gully intensity decreased significantly as the slope increased.Although the steeper slope will increase the shear force of the soil particles [67], the slope is not the main factor restricting erosion; it is more affected by slope length, catchment area, slope shape, and other terrain factors [16] when the slope is greater than 15 degrees.

The Relationship between Gully and Hillslope Erosion
The relationship between hillslope and gully erosion is illustrated in Figure 9. Gully erosion density and intensity represent gully length and area per unit area, respectively.The variation of gully erosion intensity and density is not consistent with soil erosion intensity.Gully intensity increases with slope erosion intensity, in particular rapidly increasing in the severe-violent erosion area.Gully erosion density reaches a peak at moderate erosion and then begins to decrease.This finding indicates that the control measures of hillslope erosion can effectively prevent the formation and development of erosion gully in the slight-moderate erosion area [16].However, the control measures need to be combined with other prevention and cure measures for gully erosion in The relation of hillslope and gully erosion depending on different slope aspects was analyzed (Figure 7c).The results indicate that slope aspect has a consistent influence on gully erosion intensity and hillslope erosion.The hillslope and gully erosion on a sunny slope (S, SE, SW) were stronger than those on a shady slope (N, NE, NW) considering the monsoons, total solar radiation, diurnal temperature difference, and freezing and thawing [10].

The Relationship between Gully and Hillslope Erosion
The relationship between hillslope and gully erosion is illustrated in Figure 9. Gully erosion density and intensity represent gully length and area per unit area, respectively.The variation of gully erosion intensity and density is not consistent with soil erosion intensity.Gully intensity increases with slope erosion intensity, in particular rapidly increasing in the severe-violent erosion area.Gully erosion density reaches a peak at moderate erosion and then begins to decrease.This finding indicates that the control measures of hillslope erosion can effectively prevent the formation and development of erosion gully in the slight-moderate erosion area [16].However, the control measures need to be combined with other prevention and cure measures for gully erosion in moderate-violent erosion regions, where the lateral erosion of the gully erosion is significant, which results in notable widening of ditch sides and an increase in the area.

Conclusions
To our knowledge, this study was one of the first studies to attempt to integrate gully erosion and hillslope erosion.The following main conclusions can be drawn: (1) although the study area as a whole had slight erosion due to rill and sheet erosion, the catchment suffered more serious gully erosion, which mainly occurred in the cultivated land; (2) GCM contributed to the overall improvement of the RUSLE model soil erosion assessment, particularly the C and P factor mapping and the recognition of ephemeral gully, in hilly black soil regions of Northeast China in which slope farmland is dominant.However, the RUSLE model likely overestimates the erosion rate in dry land; (3) topographic differentiation characteristics of hillslope and gully erosion are not completely similar.The hillslope and gully erosion were stronger on a sunny slope than on a shady slope, and occurred and distributed mainly at slope land of middle elevations.When the slope was greater than 15 degrees, the slope was not the main factor restricting the erosion, while at steeper slopes, the dominant forest land significantly reduced the soil loss; (4) trends of gully erosion intensity and density are not consistent with the change in soil erosion intensity.The control measures of hillslope erosion can effectively prevent the formation and development of erosion gully in the slightmoderate erosion area, but in moderate-violent erosion regions the measures need to be combined with other prevention measures for gully erosion.
The findings of this study contribute decisively towards improving our ability to better understand the spatial coupling relationships between hillslope erosion and gully erosion.The present study also indicates that using GCM, RS, and GIS technologies simultaneously results in an effective assessment of hillslope and gully erosion in a considerably short time and at low cost for large watersheds.The results might be helpful for decision makers developing sustainable land use plans and comprehensive soil and water conservation management in the study area.

Conclusions
To our knowledge, this study was one of the first studies to attempt to integrate gully erosion and hillslope erosion.The following main conclusions can be drawn: (1) although the study area as a whole had slight erosion due to rill and sheet erosion, the catchment suffered more serious gully erosion, which mainly occurred in the cultivated land; (2) GCM contributed to the overall improvement of the RUSLE model soil erosion assessment, particularly the C and P factor mapping and the recognition of ephemeral gully, in hilly black soil regions of Northeast China in which slope farmland is dominant.However, the RUSLE model likely overestimates the erosion rate in dry land; (3) topographic differentiation characteristics of hillslope and gully erosion are not completely similar.The hillslope and gully erosion were stronger on a sunny slope than on a shady slope, and occurred and distributed mainly at slope land of middle elevations.When the slope was greater than 15 degrees, the slope was not the main factor restricting the erosion, while at steeper slopes, the dominant forest land significantly reduced the soil loss; (4) trends of gully erosion intensity and density are not consistent with the change in soil erosion intensity.The control measures of hillslope erosion can effectively prevent the formation and development of erosion gully in the slight-moderate erosion area, but in moderate-violent erosion regions the measures need to be combined with other prevention measures for gully erosion.
The findings of this study contribute decisively towards improving our ability to better understand the spatial coupling relationships between hillslope erosion and gully erosion.The present study also indicates that using GCM, RS, and GIS technologies simultaneously results in an effective assessment of hillslope and gully erosion in a considerably short time and at low cost for large watersheds.The results might be helpful for decision makers developing sustainable land use plans and comprehensive soil and water conservation management in the study area.

Figure 1 .
Figure 1.The location map of the study area.(a) Location of Northeast China; (b) Jiutai District located at the southeastern edge of the black soil region; (c) spatial distribution of national meteorological stations; (d) location map of Mushi River sub-catchment with DEM data of the Jiutai District; (e) the detailed Pleiades-1A RGB image with 3, 2, 1 band combination covering the study area.

Figure 3 .
Figure 3. Schematic diagram for the ridge direction.(The black line represents the boundary between the plots and the red line indicates the direction of the ridge.)

Figure 3 .
Figure 3. Schematic diagram for the ridge direction.(The black line represents the boundary between the plots and the red line indicates the direction of the ridge).

Figure 4 .
Figure 4. Spatial distribution of factors of RUSLE.(a) The R factor map; (b) the K factor map; (c) the LS factor map; (d) the C factor map; (e) the P factor map.

Figure 5 .
Figure 5.The hillslope soil erosion map of the study area.

Figure 6 .
Figure 6.The distribution of (a) erosion gullies and (b) erosion gully density in the Mushi River subcatchment.

Figure 5 .
Figure 5.The hillslope soil erosion map of the study area.

Figure 5 .
Figure 5.The hillslope soil erosion map of the study area.

Figure 5 .
Figure 5.The hillslope soil erosion map of the study area.

Figure 5 .
Figure 5.The hillslope soil erosion map of the study area.

Figure 5 .
Figure 5.The hillslope soil erosion map of the study area.

Figure 5 .
Figure 5.The hillslope soil erosion map of the study area.

Figure 7 .
Figure 7. Relation of hillslope and gully erosion depending on different (a) elevations, (b) slopes, and (c) slope aspects.

Figure 8 .
Figure 8. Land use types in each slope class.

Figure 7 .
Figure 7. Relation of hillslope and gully erosion depending on different (a) elevations, (b) slopes, and (c) slope aspects.

Figure 7 .
Figure 7. Relation of hillslope and gully erosion depending on different (a) elevations, (b) slopes, and (c) slope aspects.

Figure 8 .
Figure 8. Land use types in each slope class.

Figure 8 .
Figure 8. Land use types in each slope class.

Table 1 .
C factor for different LULC classes.

Table 2 .
Interpretation signs on Pleiades.
˝E 44.356 ˝N Ephemeral gully

Table 2 .
Interpretation signs on Pleiades.

Table 2 .
Interpretation signs on Pleiades.

Table 2 .
Interpretation signs on Pleiades.

Table 2 .
Interpretation signs on Pleiades.

Table 2 .
Interpretation signs on Pleiades.

Table 2 .
Interpretation signs on Pleiades.

Table 3 .
Area and soil loss under different hillslope erosion classes in the study area.

Table 4 .
Status of gully erosion in the study area.

Table 5 .
[9]parison between observed data from Zhang et al.[9]and the calculation of RUSLE in this study.