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.