2. METHOD
2.1 Study area
The black soil region (BSR) is located in northeast China (42°57′~50°14′N, 122°5′~128°12′E) and covers an area of 20.6*104 km2(Figure 4A). The BSR produces about 25% of China’s corn. Related provinces include Heilongjiang, Jilin, Inner Mongolia and Liaoning. The climate is temperate and semi-humid monsoon climate, with a cold-dry winter and mild-wet summer. Mean annual precipitation ranges from 364 mm to 656 mm (Figure 4B). Mean annual temperature ranges from -0.5 ℃ to 7.1 ℃. The topography is a quasi-basin, with elevation of 79~776 m (Figure 4A). Geomorphology types include plain (56%), table-land (36%) and hill (8%). The table-land is not flat, but rather undulating. According to genetic soil classification system of Chinese Soil Taxonomy, main soil great groups of the BSR include meadow soil (36%), black soil (26%), chernozem (22%), dark-brown soil (5%) and aeolian sand soil (5%, Figure 4C). Meadow soils locate in the plain and valley bottom while black soils and chernozems mainly locate on the table-land. The original vegetation for black soils and chernozems were mainly steppe, with ground biomass of 1.1 kg/m2 and underground biomass of 1.7 kg/m2 (wet weight, CAS-FSI, 1980). Decomposition of the organic matter was historically rather slow due to the cold weather. As a result, most black soils and chernozems have fertile A horizons thicker than 30 cm and can be classified as mollisols in U.S. soil taxonomy (Figure 1). Parent materials are mainly loess and fluvial-lacustrine sediments. The difference between black soil and chernozem is that chernozem has calcium carbonates in the upper 100 cm due to its relatively drier climate. Dark-brown soils are mainly located in hill lands and classified as alfisols in U.S. soil taxonomy with parent materials of either slope deposits or residuum from granite and metamorphic rock. The BSR also includes Aeolians and soils which are deposits of wind eroded materials. Their clay content are less than 10% and SOM less than 0.5%. They are mainly located in the western BSR, where the climate is drier and topography is mainly plains and gently sloping sand dunes. Some aeolian sands have an A horizon but are generally thinner than 15 cm.
Nearing et al. (2017) provided a detailed description of the history of agricultural development in the BSR. Present landuse of the BSR is dominated by croplands (69%, Figure 4D) that are only recently cultivated (150 a). Sloping croplands cover an area of 7.2*104 km2, with black soil, chernozem and dark-brown soil as the main soil types. Crops on sloping croplands are dominated by corn. The conventional tillage system is highly mechanized and includes moldboard ploughing, disking, ridging, harrowing and seeding in April, followed by fertilizing and spraying throughout the growing season, and harvesting in October. Ridge-furrow systems are common with 65~70 cm wide and 10 cm high ridges. China statistical yearbooks suggest mean corn yield per unit area is about 6500 kg/ha in the BSR. The value was larger in south and smaller in north.
Water erosion (sheet-rill, ephemeral gully, gully) and wind erosion are all prevalent in the BSR. Runoff plot data indicates that sheet-rill erosion rate is 2~4 mm/a on bare soils and 0.7 mm/a on a 100 m long cropland (Liu et al., 2008). 137Cs data indicates that mean sheet-rill erosion rates range from 1.2 mm/a to 1.9 mm/a on three small watersheds (Fang et al., 2005; Liu et al., 2008; Fang et al., 2012). A sheet-rill erosion rate as large as 4~7 mm/a was reported for a 3°~7°slope (CAS-FSI, 1980). Density of ephemeral gullies and classic gullies was 1.1 km/km2 and 0.7 km/km2 in typical small watersheds, respectively (Zhang et al., 2016). Soil erosion in the BSR can be significant because contour farming is rare, whereas, longitudinal ridge system (LRS), in which long straight rows oriented according to the parcel’s borders, is the conventional tillage system in the BSR (Xu et al., 2018). Given that rows in the BSR tend to be extremely long with slope lengths that range from 200 m to 1000 m (Liu et al., 2011; Li et al., 2016), LRS involves undulating rows with sections inevitably aligned with the contour and other sections aligned more up and down slope. Thus, ridge directions are a combination of contour direction and up-down slope direction. The combination of the LRS on the undulating terrain inevitably results in significant slope gradients along furrows, which combined with long rows provides conditions for excessive erosion. Detachment of soils by water can be rather strong in these sloping furrows, especially in the long slope cases (Xu et al, 2018). Secondly, plant residue is rarely returned to soils. Straw is used as fuel or feed and remaining residue is generally burned in the field. Soils are essentially bare through winter until spring plant emergence during which time strong winds and storms often happened. Thirdly, high intensity mechanized tillage without subsoiling results in compacted layers for many croplands.
2.2 Extraction of sloping cropland
The GlobeLand30 dataset (30 m resolution, Chen et al., 2015) and the 1:1000000 geomorphologic atlas of the People’s Republic of China (Cheng et al., 2011) was used to extract sloping croplands. Total area of sloping croplands was found to be 71641 km2.
2.3 A horizon thickness
A horizon thickness (\(\text{Th}_{A}\)) was derived from the soil species record (SSR) of each province (Soil Survey office of Heilongjiang province; 1990; Master Station of soil fertilizer of Liaoning province, 1991; Soil Survey office of Inner-Mongolia autonomous region, 1994; Master Station of soil fertilizer of Jilin province; 1997). SSR is a dataset describing soil species, including thickness, area and main soil properties (color, bulk density, SOM and particle size etc.) of each horizon. The SSRs were published in 1990s based on data collected in a national soil survey in 1980s. The BSR includes 15 great soil-groups, 58 sub soil-groups and 720 soil species.\(\text{Th}_{A}\) and area of each soil species were extracted from the SSR. The national 1:1000000 scale soil map was used for mapping. As the minimum mapping unit of the available soil map is sub soil-group (a level higher than soil species), \(\text{Th}_{A}\) of each sub soil-group was calculated by weighing \(\text{Th}_{A}\) of each soil species by area and then assigned to the corresponding polygon in the soil map. Extrapolation across scales can cause unintentional bias, as soil thickness and properties may change considerably over short distances (Hartemink et al., 2007) but is the best approach given the data limitations. As SSR data were mainly measured in1980s, the\(\text{Th}_{A}\) in SSR was assumed to represent conditions in 1985. The \(\text{Th}_{A}\) in later years were calculated by subtracting eroded thickness from \(\text{Th}_{A}\) in 1985.
2.4 Water erosion rate
2.4.1 Soil erosion model
Only sheet-rill erosion was considered in this study and was based upon the China Soil Loss Equation (CSLE) which is similar to the USLE (Liu et al., 2002). The CSLE is:
\(E_{r}=RKLSBET\) (2)
where Ɛr is soil erosion rate (t/ha⋅a) and converted to a per area basis (mm/a) assuming a soil bulk density of 1.35 g/cm3; R is mean annual rainfall erosivity (MJ·mm/ha·h); K is soil erodibility (t·hm2·h / ha·MJ·mm); L, S, B, E, and T are dimensionless factors of slope length, slope gradient, biological-control, engineering-control, and tillage practices, respectively.
2.4.2 Rainfall erosivity
R was calculated from daily precipitation of 84 meteorological stations (1981 to 2010, China Meteorological Data Service Center) using an equation supplied by Xie et al. (2016). Daily precipitation lager than 12 mm was considered as erosive rainfall. R was calculated with following equations.
\({R_{\text{day}}=\alpha\bullet P_{\text{day}}}^{1.7265}\) (3)
\(R_{\text{year}}=\sum R_{\text{day}}\) (4)\(R=R_{\text{year}}\) (5)
where \(R_{\text{day}}\), \(R_{\text{year}}\) and R was rainfall erosivity of day, year, and mean annual rainfall erosivity, respectively, (MJ·mm/ha·h); \(\alpha\) is a parameter that equals 0.3937 in the warm season (May-September) and 0.3101 in the cool season (October-April). The R values of all 84 stations were interpolated by Kriging.
2.4.3 Soil erodibility
SOM and particle size distributions of the A horizon of each soil species were extracted from SSR. Soil texture was determined by particle size distribution according to USDA standard. Soil-structure code and permeability class were determined by soil texture (Wischmeier and Smith, 1978; Wall et al., 2002). From these data, KNvalues were calculated with nomograph supplied by Wischmeier and Smith (1978). KN values were multiplied by a coefficient to derive K values, as runoff plot data showed that there is a deviation between KN and K for some soils in the BSR (Zhang et al., 2008). The coefficient was 0.14, 0.52, 0.82 and 0.80 for brown soil, dark-brown soil, albic soil and chernozem, respectively. The coefficient is 1 for other soils.
2.4.4 Slope gradient factor and slope length factor
The S was calculated with the following equation:
\(\mathrm{S=}\left\{\par \begin{matrix}\mathrm{10.8\ }\sin{\operatorname{\theta}\mathrm{+0.03\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \theta<5}}\\ \mathrm{16.8}\sin\mathrm{\theta}\mathrm{-0.5\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 5\leq\theta<10}\\ \mathrm{21.9}\sin\mathrm{\theta}\mathrm{-0.96\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \theta\geq 10}\\ \end{matrix}\right.\ \) (6)
where \(\theta\) is slope gradient (°). The shuttle radar topography mission (SRTM) data were used to calculate \(\theta\). It is more difficult to calculate L than S. It is impractical to measure the slope lengths for all croplands in BSR. The slope length of terrain (\(\lambda_{T}\)) is rather long in the BSR, averaging 500~1000 m. But \(\lambda_{T}\) is not equivalent to the slope length of soil erosion (\(\lambda\)). The \(\lambda_{T}\) is often cut into smaller segments by features at cropland parcel boundaries in the BSR, such as windbreaks, roads and ditches, etc. Therefore, the\(\lambda\) is very close to the topographic length of each cropland parcel (\(\lambda_{P}\)). The \(\lambda_{P}\) of 3218 cropland parcels in the BSR and Heilongjiang province were measured by Lv et al. (2013). The mean \(\lambda_{P}\) was 260 m in table-land area and 152 m in hill area. These values were considered as \(\lambda\) and used to calculate the slope length factor (L) with the following equation:
\(L=\left(\frac{\lambda}{22.13}\right)^{m}\) (7)
where m is the slope-length exponent. The value of m is related to slope gradient (\(\theta\)). The m equals 0.2, 0.3, 0.4 and 0.5, when\(\theta\) range <1°, 1~3°, 3~5° and ≥3°, respectively.
2.4.5 Biological-control, engineering-control, and tillage practices
Values of B, E and T were determined from runoff plot data and field surveys. B was set to 0.5 for croplands (Fan et al., 2011), and E set to 1, as engineering-control practices are very rare in the BSR. Because directions of ridge-furrow are mostly between contour and up-down slope in the BSR, T was assumed to be the mean value of T for contour farming (TC) and T for up-down slope farming (1). TC was calculated according to a regression equation based on data supplied by Wischmeier and Smith (1978).
\(\mathrm{T}_{\mathrm{C}}\mathrm{=0.003}\mathrm{\theta}^{\mathrm{2}}\mathrm{-0.011\theta+0.553}\)R2=0.92 (8)
2.5 Soil life expectancy, erosion hazard degree and relative crop yield
The SLE of A horizon (SLEA) was calculated with Equation (1), with ThCR as 20 cm, based upon the initial thickness and the calculated Ɛr from Equation (2), according to the flow chart in Figure 5. The erosion hazard degree was determined according to SLEA and the industry standard of China (Table 1, Ministry of Water Resources of China, 2008). The relative crop yield, ratio of future crop yield to current value, was calculated with \(\text{Th}_{A}\) and a regression function (Figure 2). The flow chart of calculation is shown in Figure 5. All data layers were in a 900 m2 resolution for calculation and analysis, including \(\text{Th}_{A}\), R, K, L, S, B, E, T, Ɛr, SLEA and relative crop yield etc. However, a 900 m2 resolution is too small for visualization in this article. Therefore, a resolution of 100 km2 is used for displaying results with the value of each 10*10 km grid being the mean value of all 30*30 m grids within it. The 100 km2resolution is also the approximate size of townships in China, which is an administrative unit smaller than county.
2.6 Soil conservation scenario
To identify optimum soil conservation practices, the SLEA, \(\text{Th}_{A}\) and relative crop yield were predicted under six soil conservation scenarios (Present, Contour, Straw, Combo 1, No-till and Combo 2) which cover major soil conservation practices in the BSR. Detailed description of each scenario was shown in Table 2. The tillage practice factor (T) value was derived from runoff plot data in the BSR (Yang, 2019) or in USLE handbook (Wischmeier and Smith, 1978). The relative crop yield was predicted with future\(\text{Th}_{A}\) under each scenario and Equation (1). Residues were found to decrease soil temperature, emergence rate and crop yield in cool regions (Defelice et al., 2006). In the BSR, the intensity of this adverse effect varies spatially. In the southern (warmer) BSR, the effect of residue cover is negligible (Zhang et al., 2005; Gao et al., 2011), but in the northern (cooler) BSR, the adverse effect is significant. Compared to conventional tillage, corn yield can be reduced by 16% under No-till (Chen et al., 2011) and 11% under Straw (Zou et al., 2016). These ratios were used to estimate future crop yield under Straw, No-till and Combo 2 in the cooler BSR areas, including Heilongjiang province and Inner Mongolia (Figure 4A).The relationships incorporated in the model for NT and Straw may not fully represent their effectiveness as the database may be skewed towards short term data. However, given the data available for the CSLE and crop yield response, the model predictions in this study reflect the best estimates currently available.