Assessment of soil erosion at multiple spatial scales following land use changes in 1980-2017 in the black soil region, (NE) China

Impact of land use and land cover (LULC) change on soil erosion is still imperfectly understood, especially in northeastern China (NEC). Based on the Revised Universal Loss Equation (RUSLE), the variability of soil erosion at different spatial scales following land use changes in1980, 1990, 2000, 2010, and 2017 was analyzed. The regionally spatial patterns of soil loss coincided with the topography, rainfall erosivity, soil erodibility, and use patterns, and around 45% soil loss came from arable land. Regionally, soil erosion rates increased from 1980 to 2010 and decreased from 2010 to 2017, ranging from 3.91 to 4.45 t ha-1 yr-1 with an average of 4.22 t ha-1 yr-1 in 1980-2017. The rates of soil erosion less than 1.41 t ha-1 yr-1 decreased from 1980 to 2010, and increased from 2010 to 2017, and opposite changing patterns occurred in higher erosion classes (i.e., above 5 t ha-1 yr-1). At a provincial scale, Liaoning Province experienced the highest soil erosion rate of 9.43 t ha-1 yr-1, followed by Jilin Province, the east Inner Mongolia, and Heilongjing Province. Arable land continuously increased at the expense of forest in the high-elevation and steep-slope areas from 1980 to 2010, and decreased from 2010 to 2017, resulting in increased areas with erosion rates higher than 7.05 t ha-1 yr-1. At a county scale, around 75% of the countries had soil erosion rate higher than its tolerance level. The county numbers with higher erosion rate increased in 1980-2010 and decreased in 20102017, resulting from the sprawl and withdrawal of arable land. The results indicate that appropriate policies can control soil loss through limiting arable land sprawl in areas of unfavorable regions in the NEC.

mountainous areas, and broad leaved forest scatter in the plains. Arable land in the plains covers about one third of the total area. The main crops are maize, soybean, and rice, occupying around 87% of the total. A single tillage operation is typically used with a cultivation harrow with 0.25 m depth in late-September and in early-May. From October to April, arable land is left fallow with no vegetation cover during cold weather.
The main soil types are Luvisols, Phaeozems, Chernozems, and Gleysols, occupying 37.5%, 23.3%, 7.4%, 6.8% of the total area in the NEC (Figure 1d). The parent materials are Quaternary lacustrine and fluvial sand beds or loess sediments [20]. Fuvisols are widely distributed in mountainous areas, while Phaeozems and Chaeozems which are called black soils by local people mainly appear in the Songnen Plain [9]. The textural classes of the topsoil are silt-clay-loam to clay-loam (8-27% sand, 29-66% silt and 26-40% clay), with soil organic carbon usually above 4%. The hilly regions with long slope in the mountain-plain transition zone are prone to soil erosion. Rill erosion dominates the study area. However, gullies also widely occur [16]. Figure 1 is about here.

Data collection
The RUSLE (Revised Universal Soil Loss Equation) was employed to estimate soil erosion rate. Several data layers are required to run this model: a Digital Elevation Model (DEM), land use maps, a rainfall erosivity map, a soil erodibility map, the land use management and soil conservation practice maps.
The DEM was prepared from Shuttle Radar Topography Mission (SRTM) with 90-m resolution. It was provided by Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn).
The Harmonized World Soil Database (HWSD) Version 1.2 was obtained from Food and Agriculture Organization of the United Nations (FAO-UN; http://www.fao.org). Soil properties such as percentages of sand, clay and silt, depth of soil, organic carbon are available from the HWSD.
The daily precipitation data in 1970-2018 at the 120 meteorological stations were acquired from National Climate Centre of the China Meteorological Administration (http://data.cma.cn/).
The land use datasets with 30-m resolution were acquired from different sources. The datasets in 1980, 1990 and 2010 were obtained from the Resource and Environmental Sciences, Chinese Academy of Sciences (RES-CAS). The datasets in 2010 and 2017 were downloaded from the GlobeLand30 (http://www.globallandcover.com) and from the research team of Gong Peng in Tsinghua University (http://data.ess.tsinghua.edu.cn/fromglc2017v1.html).

Description of the RUSLE and data treatment
In the present study, soil erosion intensity was calculated using RUSLE. It combines five factors, including rainfall erosivity, soil erodibility, slope length and gradient, land use management, and soil conservation practices: where A is the soil erosion rate (t ha -1 yr -1 ), R is the rainfall-runoff erosivity factor (MJ mm ha -1 h -1 ), K is the soil erodibility factor (t ha h ha -1 MJ -1 mm -1 ), LS is a combination of slope length L and slope gradient S factor (-), C is the crop/cover management factor (-), and P is the soil conservation factor (-).
The classic method to calculate R factor is to multiply rainfall energy and the maximum 30-min rainfall intensity for rainfall events [21]. However, such datasets are not available for many places in the world. In such instances, alternative method has been proposed to calculate R factor using readily available precipitation data [22][23]. In the current study, the daily precipitation data in 1970-2018 was used to estimate mean annual R-factor values, using the method proposed by Zhang et al. [23]: Where Ri is the half-month R-factor (MJ mm ha -1 h -1 yr -1 ), and Dj is the erosive rainfall day j. Dj is equal to actual rainfall when it is greater than 12 mm. Otherwise it equal to zero. Here, a and b are empirical parameters. Detailed description was given by Fang and Sun (2017). A Co-Kriging interpolation was used to derive an R-factor value map in ArcGIS 10.5 software. In the NEC, mean annual R-factor values ranged from 475 in western NEC to 5320 MJ mm ha -1 h -1 yr -1 in southeastern NEC (Figure 2a). Figure 2 is about here. Soil texture and soil organic matter of the top 30 cm soils were available from the HWSD, and the K factor values were calculated by using EPIC method. This method has been successfully used in the NEC 10,17 : where SAN is the sand content (%), SIL is the silt content (%), CLA is the clay content (%), SC is the soil organic carbon content (%) and SN = 1-SAN/100. The K factor values of the study region varying from 0 to 0.045 t ha -1 h MJ -1 mm -1 ha -1 (Figure 2b).
The L-factor measures the slope length, and S-factor is proportional to the local slope. This study used the 2D approach to calculate LS factor as proposed by Van Rompaey et al. 24 , using the methods by McCool et al. 25 where Li,j is the slope length factor of the grid cell with coordinates i and j, Ai,j is the contribution area at the inlet of a pixel (m 2 ), D is the grid cell size (m); , is the slope aspect direction for a pixel; m is slope length exponent, β is empirical factor (-), and is the slope angle. The derived LS factor values ranged from 0.03 to 70.02. Higher values occurred in the mountainous regions, and lower values in the plains (Figure 2c).
The C and P factors closely link to land use types [27]. In order to obtain reasonable C-and P-factor values, a detailed survey of the literatures published in the NEC was first carried out [19,[28][29][30][31]. Then, C-factor values in forest, grassland, shrub land, cultivated land, and residential area were assigned to 0.004, 0.043, 0.07, 0.228, 0.03, and those in bare land, wetland and water bodies, and other unused land were 1, 0, and 0.06, respectively. For P-factor values, cultivated land and water body were set to 0.352 and 0, and other lands were unit in the study area.

LULC changes in the NEC
The land uses in 1980-2017 were shown in Figure 3 and numerically in Table 1. Spatially, forests were distributed in mountainous areas, and arable lands in the three plains. Among the land use classes, forest was predominant and covered 37% of the whole area, followed by arable-and grass-land covering 30% and 22%, respectively. These three kinds of lands occupied around 90% of the total. Shrub land, wetland, water body, and unused lands each covered 2% of the total area. Table 1 is about here. Figure 3 is about here. During the past decades, great changes in arable land and forest occurred. Before 2010, arable land percentages increased from 27% in 1980, 29% in 1990, 30% in 2000, to 35% in 2010, and then decreased to 32% in 2017, fluctuating from 41.34 million km 2 in 1980 to 45.63 million km 2 in 2010. In contrast, forest area percentages decreased from 38% in 1980, 36% in 1990, 37% in 2000, to 34% in 2010, and then increased to 37% in 2017. The area percentages of grass land decreased too before 2000, ranging from 21% in 1980, 20% in 1990, and 19% in 2000, and then increased sharply thereafter, occupying around 25% in 2010-2017. Shrub land percentages occupied around 4% before 2000, and less than 0.1% thereafter. The percentages of residential area in general increased from 1% in 1980 to 2% in 2017. On the contrary, the area percentages of wetland and water body decreased from 3% and 3% in 1980 to zero and 1% in 2017, respectively. Unused lands kept unchanged through time. Since 2010, bare land has almost disappeared. Noticeably, great LULC occurred in 2000-2010 (i.e., P3 in Table 1), except for those of the wetland and bare land in 2010-2017.

Changes in soil erosion
In 1980-2017, spatial patterns of annual soil erosion rates changed little ( Figure 4). Annual soil erosion rates in the NEC ranged from 0 in the plains to over 4000 t ha -1 yr -1 in the mountainous areas, with an average of 4.22 t ha -1 yr -1 (Figure 5a). It was three times the soil loss tolerance (T=1.41 t ha -1 yr -1 ; [6]). Annual soil loss reached around 5.24 million ton from the NEC. The soil erosion rates were classified into five erosion classes from low (0-T), moderate (1-5T), high (5-10T), very high (10-20T) to critical (>20T). Figure 5b indicated that 57% of the NEC belonged to low erosion class, followed by moderate erosion with area covering around 30%. The high, very high and extreme erosion areas covered 6%, 3.5%, and 2.7% of the total, respectively. Figure 4 is about here. Figure 5 is about here. The area percentages with different erosion classes differed with time ( Figure 5). The area percentages in the low erosion class decreased from 59% in 1980 to 56% in 2010, and increased to 57% in 2017. Opposite changing patterns occurred for the area percentages with high, very high, and extreme erosion severities. The threshold values occurred in 2010. In contrast, the area percentages with moderate erosion severity increased continuously with time.
The rates of annual soil erosion presented an obvious spatial heterogeneity among the provinces and the EIM ( Figure 5). The mean soil erosion rate was the highest in Liaoning Province (i.e., 9.43 t ha -1 yr -1 ) in 1980-2017, followed in Jilin Province, and the EIM. The lowest erosion rate occurred in Heilongjiang Province, with an average of 2.46 t ha -1 yr -1 (Figure 5c). Considering the area percentages, their annual soil losses accounted for 21%, 17%, 26%, and 35% of the total NEC, respectively.
Temporally, the soil erosion rates of the three provinces increased from 1980 to 2010 and then decreased from 2010 to 2017. The rates of soil erosion in Heilongjiang and Jilin provinces in 2017 were less than the counterparts in 1980. The annual erosion rate in the EIM fluctuated from 1980 to 2017 (Figure 5d).
Low and moderate erosion zones dominated the area, with similar temporal patterns to the NEC ( Figure 6). Specifically, the areas with low erosion class dominated Heilongjian Province, Jilin Province and the EMI, and the areas with moderate erosion class were dominant in Liaoning Province. In the three provinces, the area percentages with low erosion class decreased from 1980 to 2010, and increased from 2010 to 2017. In contrast, it continuously decreased with time in the EIM. The area percentages with moderate erosion class increased with time for Jilin Province and the EIM. In contrast, it increased in 1980-2010 and decreased from 2010 to 2017 for Heilongjiang Province, and opposite changing trend in Liaoning Province. Figure 6 is about here. The spatio-temporal characteristics of soil erosion severity were mapped for the 224 counties in the NEC (Figure 7). In 1980-2017, the number of the counties with soil erosion rates higher than one T accounted for around 75% of the total. The counties with soil erosion rates higher than 5T and 10T accounted for 24% and 8% of the total, respectively ( Figure 8). For the top 16 counties with erosion rate higher than 10T (i.e., 14.1 t ha -1 yr -1 ), 12 and 4 counties were distributed in the southeastern and southwestern NEC, respectively. Almost all the counties with soil erosion rate lower than one T were distributed in the plains (Figs. 1 and 9). Figure 7 is about here. Figure 8 is about here. Arable land continuously sprawled from south to the north in1980-2010, and a little withdraw occurred in 2017. In 1980-2010, among the 224 countries, there were 52 countries in low erosion zone (i.e., less than T), 119 countries with moderate erosion zone (1-5T), 37 countries with high erosion zone (5-10T), 14 and 2 countries with very high (10-20T) and extreme erosion (> 20T) zones, respectively (Figure 8a). The number percentages of the countries with low, moderate, high, very high, and extreme zones occupied 23%, 53%, 16%, 6%, and 1% respectively.
In 1980-2017, the number percentages of the counties with low erosion rate ranged from 20% to 25% of the total and kept unchanged after 2000. The number percentages of the counties in moderate soil erosion zone decreased from 54% in 1980 to 49% in 2017, among which the soil erosion rate with 1-2T decreased fastest, from 21% in 1980 to 17% in 2017. Noticeably, the numbers of counties with high and very high classes increased from 14% and 5% in 1980 to 18% and 8%, respectively ( Figure  8b).

Model performance
In the NEC, soil erosion rates have been identified using runoff plot [14], 137 Cs tracer technique [15,[32][33], and model estimation [10,29]. Comparisons of the RUSLE-derived rates of soil erosion in the present study to published results can evaluate the performance of the RUSLE.
Some published results are from runoff plots, natural slopes, and catchments. At plot and slope scales, the RUSLE-derived rates of soil erosion in the current study were less than the reported results (Table 3). This difference could result from two aspects. One aspect is that the 90-m resolution topographic data can smooth slope gradient [6,34], resulting in lower soil erosion rates. The other aspect is the difference in time spans between the studies in literature and the current study. The rates of soil erosion have also been reported at larger scales, including at a 28.5 ha catchment in Heshan Farmland [15], at a 916 km 2 Shuangyang river catchment 10 , and in the Sanjiang Plain [29]. A comparison at a catchment scale is more meaningful since this study was done at a large spatial scale. Table 3 demonstrated that the RUSLE-derived rates of soil erosion in the current study are basically consistent with the published results. Table 3 is about here.

Impact of LULC change on soil erosion
Land use greatly influences soil erosion. The rates of soil erosion on bare land, shrub land, and arable land were much higher than those on other land use types, and the total soil loss amounts from these lands changed from 75% in 1980, 96% in 2010, to 79% in 2017 of the total in respect years. Table 2 is about here. These changes can be attributed to the LULC with time, due to population pressure and national agricultural policies [5,11,35]. Since the mid-20 century, extensive settlement and land exploitation have begun in the NEC, and a large amount of woodland were converted into farmland 11 . From 1980 to 2000, around 4% forest and 3% wetland were lost, and around 3381 km 2 woodland was converted to arable land in 1990-2000 12 . As a result, soil loss from arable land accounted for 42% of the total, followed by grassland, forest land, and other land uses ( Table 2). The NEC is the Chinese most important grain production base, and a series of policies further stimulated arable land reclamation in 2000-2010. For example, the rescission of the tax from agriculture in 2005 and the household responsibility system had triggered farmland cultivation at the expense of woodland (i.e., forest and shrub lands), resulting in a 46% soil loss of the total from arable land in 2010. To decrease soil loss, the "Grain to Green" project was introduced, and reforestation was initiated in the NEC [12,36], including farmland shelterbelt, sand fixation forest, and conversion of farmland to forest on steep slopes [11]. From 2010 to 2017, forest area increased by 39,928 km 2 , with the largest change ratio of 9% during 2010-2017 (Table 2).

Impact of topography on soil erosion
The spatial distribution of soil erosion rate and their changes with time can be explained by topography in the NEC. First, steep mountainous areas usually have higher LS values 36 . This distribution characteristics has been found in literatures [1,[37][38]. Second, more rainfall and higher rainfall erosivity occurred in the mountainous areas (Figure 2a), and lower values in the plains ( Figure 3). The spatial pattern of R value in mountainous regions can also be found in Sudetes Mts., SW Poland 1 . In the NEC, soil erosion rates increased with increasing slope gradients (Figure 9ab). The rate of soil erosion above 300 m a.s.l. was four to five times that in areas below 200 m a.s.l. Due to steep slope and high R values, the rates of soil erosion were high in high-elevation and steep mountainous regions, and low in plains (Figure 9ab). Therefore, the spatial distributions of soil erosion rate at regional and county scales were quite related to topographic factors (Figs.1 and 5), which is consistent with the studies by Bakker et al. [39] and Latocha et al. [1]. This can explain the spatial distribution of soil erosion rate at regional and county scales in the NEC ( Figures.5 and 9).
In areas with elevation above 200 m a.s.l and slopes above 5 degree, the distribution of arable land showed similar change patterns in different years, i.e., the areas of arable land increased from 1980-2010 and decreased in 2017. In contrast, in the low-elevation and gentle-slope regions, arable land area increased continuously (Figure 9cd). Because arable land is the major sediment source in the NEC and the soil erosion rate is high in the high-elevation and steep-slope areas (Table 2)

Study limitations
The soil erosion rate was calculated through multiplying the five factors (i.e., R, K, LS, C, and P).The data accuracies, especially the values of the C and P factors, greatly influence the results. Without consideration of the varied C-and P-factor values with time could result in biased estimation [6]. Furthermore, the datasets with 90 m resolution could influence its estimation accuracy at plot or slope scales. Increase in computational capacity and high spatial resolution data could also improve estimation accuracy.
Except for data resolution, the preclusion of gully erosion from the RUSLE can also induce biased estimation. In the NEC, gullies are widely distributed, and over 295,000 (ephemeral) gullies were found in this region [40]. The soil loss from gullies is high. In a reservoir catchment in Baiquan County, the sediment from gullies occupied around 40% of the total [41]. Similar result was also obtained by Huang et al. [42] in Heshan Farmland, NEC. At a catchment scale, models like WaTEM/SEDEM [24] and TeTIS [16] have combined soil erosion module and sediment deposition module together to simulate soil sediment yield. The combination of such module with RUSLE could improve the estimation accuracy.
The uncertainty also comes from the applications of LS calculation equations. This kind of limitation has fully been explained by Alewell et al. [43], and thus detailed explanation was not given here. However, it should be stated that the application of the upslope area-based L factor calculation could improve the estimation accuracy.
However, although such limitations and uncertainty mentioned above could to some extent influence our estimation accuracy (Table 3), the comparisons with published results verified the feasibility of our estimated results. Furthermore, since the focus of this study was to estimate the impact of LULC on soil erosion with time, these uncertainties can be accepted.

Conclusions
Assessment of the impact of regional LULC change on soil loss is paramount importance to appropriately regulate land use. Soil loss characteristics and its changing patterns following LULC changes in 1980-2017 were explored using RUSLE method at regional, provincial, and county scales.
In the NEC, higher rates of soil erosion occurred in the mountainous regions, and lower values in the plains. The spatial patterns of soil loss coincided with those of RUSLE-LS, -C, -K, and -R factors in the NEC (Figs. 2 and 4). Impacted by the LULC in 1980-2017, threshold values of soil erosion rate occurred. The regional soil erosion rates increased from 3.91 t ha -1 yr -1 in 1980 to 4.45 t ha -1 yr -1 in 2010 and decreased to 4.22 t ha -1 yr -1 in 2017. The rates of soil erosion less than the tolerance value (i.e., 1.41 t ha -1 yr -1 ) decreased from 1980 to 2010, and increased in 2017. In contrast, opposite changing patterns occurred in higher erosion classes. At provincial scale, due to higher R values in the southeast and southwest of the NEC, the highest soil erosion rate of 9.43 t ha -1 yr -1 occurred in Liaoning Province, and the lowest value in Heilongjing Province. Due to the sprawl of arable lands at the expense of natural vegetation from 1980 to 2010, the area percentages with soil erosion rate higher than 7.05 t ha -1 yr -1 increased before 2010, and decreased after 2010. At a county scale, the counties in the southeast and southwest of the NEC experienced the highest soil loss. Around 75% of the counties suffered soil loss with soil erosion rate higher than its tolerance level. The areas with soil erosion rate less than T increased in 1980-2010 and decreased in 2017. Therefore, appropriate policies can control soil loss through limiting the sprawl of arable lands in unfavorable regions.     Table 3. The published rates of soil erosion in the black soil region (the study area locations of the references are labeled in Figure 1). Note: * The average bulk density 1.3 t m -3 was used when the unit mm yr -1 was converted to unit t ha -1 yr -1 according to Xie et al. (2019); CGWCR indicates the Complication Group of Water Conservancy Record in Baiquan County; UP recurrents up-down cultivated land. The first numbers in the column "Region" recurrent the counties labeled in Figure 1.