Sustainability of soil organic carbon in consolidated gully land in China’s Loess Plateau
SCALE model
Combining the processes of biogeochemical transformation, soil erosion/deposition and resultant landscape evolution, and bioturbation by soil fauna, the governing equations, which conserves SOC mass, can be expressed as24:
$$begin{aligned} frac{ partial }{ partial t } int _{0}^{Z} mathbf{C } dz =int _{0}^{Z} mathbf{g } dz – nabla cdot mathbf{q }_C + int _{0}^{Z} nabla cdot big [D(z) nabla mathbf{C } big ] dz end{aligned}$$
(1)
where (mathbf{C }) is the SOC concentration ([M L^{-3}]), (mathbf{C } = [ C_l,C_h,C_b]^T) represents the fast (or litter), slow (or humus), and microbial biomass pool, respectively; (mathbf{g }) is the rate of the biogeochemical transformation process, which is a function of the fast (or litter), slow (or humus), and microbial biomass pool; (nabla cdot mathbf{q }_C) is the surface SOC flux associated with soil transport; and (D(z)) is the bioturbation diffusivity. When the vertical column is discretized into several layers, the equations can be written as:
$$begin{aligned} text {Top soil layer: }&frac{ partial big ( mathbf{C }_1 z_1 big ) }{ partial t } = mathbf{g }_1 z_1 – nabla cdot mathbf{q }_C end{aligned}$$
(2a)
$$begin{aligned} text {Sublayers: }&frac{ partial mathbf{C }_n }{ partial t } = mathbf{g }_n + nabla cdot big [D(z) nabla mathbf{C } big ] dz end{aligned}$$
(2b)
where the subscripts 1 and n denote the surface soil layer and the (n{mathrm{^{th}}}) layer below-surface, respectively. Other details of these equations can be found in24.
The mechanisms of soil transport and the resultant landscape evolution can be categorized into two groups: overland flow-driven transport and diffusion-driven transport from other disturbances (e.g., wind, animal, and raindrop splash). The 2-D mass conservation equation of soil transport and the resultant landscape evolution follows Exner equation:
$$begin{aligned} frac{partial eta }{partial t} = U -nabla cdot q_d – nabla cdot q_s end{aligned}$$
(3)
where (eta) is soil surface elevation [L]; U is the rate of tectonic uplift or glacial rebound ([L T^{-1}]); (q_d) is the volume flux of sediment per unit width by hillslope diffusion ([L^2 T^{-1}]); (q_s) is the volume flux of sediment per unit width by overland flow ([L^2 T^{-1}]).
The diffusion-driven transport (nabla cdot q_d) is a combination of wind erosion, animal disturbance, soil creep, raindrop splash, and biogenic transport. The 2-D equation of (q_d) is expressed as a linear relationship with slope47:
$$begin{aligned} q_d = – D_{x} frac{ partial eta }{partial x} – D_{y} frac{ partial eta }{partial y} end{aligned}$$
(4)
where (D_x) and (D_y) are the soil diffusion coefficient in x and y direction, respectively ([L^2 T^{-1}]). The overland flow-driven transport (nabla cdot q_s) is a combined form of the divergence of stream-power but limited by the detachment capacity48:
$$begin{aligned} nabla cdot q_s = min Bigg ( D_c, frac{q_{s,out} – sum q_{s, in}}{d_s} Bigg ) end{aligned}$$
(5)
where (D_c) is the detachment capacity, which is the upper limit of of local erosion rate ([L/S]); (q_{s,out}) is the sediment flux out of a cell and (sum q_{s, in}) is the total sediment flux into a cell assumed at sediment transport capacity.
The rate of change of SOC on the surface driven by soil transport is (nabla cdot q_s), which has a linear relationship with soil transport flux:
$$begin{aligned} nabla cdot mathbf{q }_{C} = nabla cdot ( k_{soc} mathbf{C }_{1} q_d) + nabla cdot ( k_{soc} mathbf{C }_{1} q_s) end{aligned}$$
(6)
where (mathbf{C } = [C_l , C_h , C_b ]^T), and the subscript 1 denotes the surface soil layer; (q_d) and (q_s) are soil transport flux of diffusion and overland flow; (k_{soc}) is an enrichment ratio, which represents a preferential transport (mobilization and deposition) of SOC. The SOC fluxes driven by diffusion and overland flow sediment transport are:
$$begin{aligned} nabla cdot ( k_{soc} mathbf{C }_{1} q_d)= & {} – frac{partial }{partial x} bigg (k_{soc} mathbf{C }_{1} D_{x} frac{ partial eta }{partial x} bigg ) – frac{partial }{partial y} bigg (k_{soc} mathbf{C }_{1} D_{y} frac{ partial eta }{partial y}bigg ) end{aligned}$$
(7)
$$begin{aligned} nabla cdot big (k_{soc} mathbf{C }_{1} q_s big )= & {} left{ begin{array}{ll} k_{soc} mathbf{C }_{1} D_c , &{}quad mathrm{if} D_c < frac{q_{s,out} - sum q_{s, in}}{d_s} \ frac{ k_{soc} mathbf{C }_{1,out} q_{s,out} - sum big ( k_{soc}mathbf{C }_{1,in} q_{s, in} big ) }{d_s} , &{}quad mathrm{otherwise} end{array} right. end{aligned}$$ (8) The initial conditions include land surface elevation, SOC profiles, soil moisture, and surface water depth for each grid box. Table S1 (in Supplementary Information) lists variables associated with initial conditions. We simulate the SOC surface transport and vertical transformation at a daily time step with 2 (times) 2 (mathrm{m}^2) spatial resolution on the surface and a range of vertical grid sizes varying from 5 to 60 (text{cm}). The initial elevation is from lidar (Light Detection and Ranging) DEM (digital elevation model) with a 2-m resolution, collected by a drone in 2015. The SOC profiles at each 2-D grid are estimated by combining the surface SOC contents from the field survey in 2016 and SOC profiles from soil cores sampled at a nearby site 2-km away (see Section below). We assume that the initial soil thickness is 1-m because the most active soil thickness49 regarding SOC and soil moisture are within the top 1 (text{m}), and the deeper SOC contents are quite uniform19,20, indicating a relatively stable condition. Also, we neglect bedrock weathering, and therefore, the soil thickness change is only from the surface soil erosion or deposition. The SOC stock and soil thickness change in the results are compared with the initial SOC stock in the top 1 m. Other initial values such as soil moisture profile and surface water depth (i.e., initial values assigned spatially uniform at each grid point) have a short-time memory in that the impacts only last for a few days to weeks and are predominantly determined by the external forcing and resulting dynamics (more information in Supplementary Information Table S1). Initial soil organic carbon profiles SOC in a natural setting exponentially decreases with soil depth. Here we assume the following relationship between SOC content and soil depth: $$begin{aligned} SOC(Z) = a e^{-bZ} + c end{aligned}$$ (9) where (text{SOC}(text{Z})) represents SOC content at depth (text{Z}); (text{Z}) is zero at the surface, and positive downward; (a), (b), and (c) are positive and constant parameters, where (a+c) represents surface SOC content, (b) represents the decay rate, and (c) represents the relative stable or immobile SOC at depth. Parameters (b) and (c) in this study are estimated from twenty deep cores—ten sites at hills of the Reference Watershed, eight sites at GLC Watershed in the hills, and two sites in the consolidated gully land (Figure S7a). The SOC contents are shown as dots in Figure S6 (in Supplementary Information). We assume the exponential decay rate ((b)) and the immobile SOC ((c)) are spatially uniform in the cornfield in the consolidated gully area and natural field, respectively. We use the least square non-linear method to fit the sample points, and the fitted curves are solid lines in Figure S6 (in Supplementary Information) with the corresponding relationships obtained as: $$begin{aligned} text {For trees/shrubs: }&SOC(Z) = 8.04 e^{-9.63Z}+2.55 end{aligned}$$ (10a) $$begin{aligned} text {For crops: }&SOC(Z) = 3.71 e^{-7.06Z}+2.27 end{aligned}$$ (10b) hence in the natural area, (b= 9.63) and (c = 2.55); and in the cornfield, (b = 7.06) and (c = 2.27). The parameter (a) varies spatially. Soil samples near the soil surface were collected across the whole area of both the GLC and Reference Watershed in 2016. To ensure that the sampling sites are uniformly distributed and represent all land cover types in each watershed, we superimposed an 80-m (times) 80-m grid on the DEM map. In each grid cell, we selected one representative site to collect soil samples; in consolidated gullies, adjacent sampling sites were spaced at an interval of 40-m because of the relatively narrow width. SOC was measured in a laboratory by using the dichromate oxidation method50. In the GLC Watershed, soil samples were collected at 89 locations with 178 samples (0–10 cm and 10–20 cm); in the Reference Watershed, soil samples were collected at 72 locations with 144 samples (Figure 2a, in Supplementary Information). We used Kriging51 to interpolate the SOC content in each 2-D grid box in the two watersheds. Then (a) is back-calculated with the given values of the two layers at 5 cm and 15 cm (the middle point of the two layers, respectively). The final surface SOC (which equals (a+c)) is shown in Figure S7a (in Supplementary Information). Forcing data To explore the co-evolution of the landscape and the vertical profiles of SOC over the decadal time scale, we target a 50-years simulation. The meteorological data is collected from the China National Field Observation Station in An’sai ((36^{circ } 51^{prime } 30^{prime prime } N), (109^{circ } 19^{prime } 23^{prime } E); data source: http://asa.cern.ac.cn), 44 km away to the northwest from the two watersheds, with a 10 years record from 2008 to 2017, which is the best available data for the simulation. The mean annual precipitation is 560 mm. These data are used to train a stochastic Weather Generator52, which is used to create an ensemble of another 40 years of data (Figure S7c in Supplementary Information). Landcover is also obtained from the field survey in 2015 (Figure S7b in Supplementary Information). It is essential for simulating surface water runoff and the input of SOC from plant residues. Different types of landcover provide different fractions that control the surface water runoff velocity, and such fraction is represented by Manning’s coefficient (Table S2 in Supplementary Information). The plant residues include dead leaves, roots, stems, and corn stover after harvest. Here, we estimate plant residues as a function of the Normalized Difference Vegetation Index (NDVI) (Figure S3a) which is obtained from Landsat satellite data (see more information in the Section below). Litter input estimation The NDVI is collected from Landsat satellite data for a full 2 years period, 2016 and 2017. It is spatially divided into three areas, the consolidated gully land, the GLC Watershed, and the Reference Watershed. The spatial distribution of NDVIs for the two watersheds are nearly identical, so we took the spatial means of NDVI for the two watersheds excluding the consolidated gully to represent the natural area (Figure S3a in Supplementary Information). During the growing season, the NDVIs in a natural area (Figure S3(a1) in Supplementary Information) are smaller than the one in the consolidated gully (Figure S3a2 in Supplementary Information). This is because the crop inside the consolidated gully has higher vegetation density. We assume the rate of litterfall has an exponentially increasing relationship with NDVI (Figure S3b in Supplementary Information). As the NDVI increases, the plant residues increase in general. When a plant’s growth slows down, the NDVI increase also slows down and is close to the maximum, but on the other hand, the litterfall rate increases much faster near the end of the growing season. This characterization allows us to fill in the gap due to the lack of litterfall data about the various land vegetation types, including the cornfield in the consolidated gully. Governing equations of SOC transformation The equation below shows the SOC transformation, which is directly affected by litter input and decomposition rate29: $$begin{aligned} frac{partial big ( C_l + C_h + C_b big ) }{partial t} = I_{litter} - big (r_r K_l C_l+r_r K_h C_h big ) end{aligned}$$ (11) where (C_l), (C_h), and (C_b) are defined in Eq. 1; (I_{litter}) is the litter input from the sum of above-ground litter fall and below-ground root-litter ([M L^{-2}T^{-1} ]); (r_r) defines the fraction of decomposed organic carbon to (hbox {CO}_2) ([-]) ((0 le r_r le 1-r_h)), which typically ranges between 0.6 and 0.8; (K_l) and (K_h) are rates of carbon decomposition in fast and slow pool, respectively ([T^{-1}]). They are regulated by soil moisture and carbon-nitrogen ((C/N)) ratio as shown below29: $$begin{aligned}&K_l = varphi (C/N) f_d(theta ) k_l C_b end{aligned}$$ (12a) $$begin{aligned}&K_h =varphi (C/N) f_d(theta ) k_h C_b end{aligned}$$ (12b) where (k_l) and (k_h) represent the rate of decomposition as a simplified term that encompasses different organic components in the litter and humus pool, respectively ([L^3 T^{-1} M^{-1}]); (varphi (C/N)) is a ratio that is from the reduction of the decomposition rate if the immobilization (controlled by nitrogen content) fails to meet the nitrogen demand by the microbes ([-]). (varphi approx 1) in agricultural fields where nitrogen supply is usually sufficient from fertilizers; (f_d(theta )) ([-]) represents the soil moisture effects on decomposition29. The optimal soil moisture condition is the field capacity which provides the highest (f_d)29, meaning that very dry or very wet conditions will result in a smaller (f_d), and hence reduce the decomposition rate. The decomposition rates for litter (or fast) pool and humus (or slow) pool from the equations are (r_r K_l) and (r_r K_h), respectively. In this study, we test the different mean residence times on the surface soil layer (5-cm) by assigning a new decay rate of the decomposition parameter (k_h) in humus (or slow) pool. A complete list of parameters can be found in Supplementary Information Table S2. More
