Sustainability of soil organic carbon in consolidated gully land in China’s Loess Plateau

Massive gully land consolidation projects, launched in China’s Loess Plateau, aim to restore 2667 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{km}^2$$\end{document}km2 agricultural lands in total by consolidating 2026 highly eroded gullies. This effort represents a social engineering project where the economic development and livelihood of the farming families are closely tied to the ability of these emergent landscapes to provide agricultural services. Whether these ‘time zero’ landscapes have the resilience to provide a sustainable soil condition such as soil organic carbon (SOC) content remains unknown. By studying two watersheds, one of which is a control site, we show that the consolidated gully serves as an enhanced carbon sink, where the magnitude of SOC increase rate (1.0 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{g\,C}/\mathrm{m}^2/\mathrm{year}$$\end{document}gC/m2/year) is about twice that of the SOC decrease rate (− 0.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{g\,C}/\mathrm{m}^2/\mathrm{year}$$\end{document}gC/m2/year) in the surrounding natural watershed. Over a 50-year co-evolution of landscape and SOC turnover, we find that the dominant mechanisms that determine the carbon cycling are different between the consolidated gully and natural watersheds. In natural watersheds, the flux of SOC transformation is mainly driven by the flux of SOC transport; but in the consolidated gully, the transport has little impact on the transformation. Furthermore, we find that extending the surface carbon residence time has the potential to efficiently enhance carbon sequestration from the atmosphere with a rate as high as 8 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{g\,C}/\mathrm{m}^2/\mathrm{year}$$\end{document}gC/m2/year compared to the current 0.4 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{g\,C}/\mathrm{m}^2/\mathrm{year}$$\end{document}gC/m2/year. The success for the completion of all gully consolidation would lead to as high as 26.67 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{Gg\,C}/\mathrm{year}$$\end{document}GgC/year sequestrated into soils. This work, therefore, not only provides an assessment and guidance of the long-term sustainability of the ‘time zero’ landscapes but also a solution for sequestration \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_2$$\end{document}CO2 into soils.

Soil and SOC transport within watersheds and the consolidated gullies Figure S1 shows soil thickness changes after 50 years evolution in the Reference watershed, GLC (gully land consolidation) watershed excluding the consolidated gully, and the consolidated gully, respectively. Here we assume a zero weathering rate for the soil column in the loess plateau region 1 , and, therefore, the thickness change is from the surface soil erosion or deposition only. Severe erosion happens at the sharp slope transition in the upland gullies for both watersheds. However, the soils detached from erosional sites are mostly deposited at a nearby location within the watersheds possibly due to micro-depressions. Hence, the soil thickness change has a relatively large spatial range compared to the spatial mean values. In general, the consolidated gullies is under a depositional environment. The strongest deposition zone is near the two edges of consolidated gullies along the upland. There is no clear strong erosion zone in the consolidated gullies. Moderate erosion and deposition happen throughout the consolidated gullies. Figure S1: The total soil depth change after 50 years of co-evolution in the (a1) Reference Watershed, (b1) GLC Watershed, and (c1) consolidated gully land. The corresponding probability distribution functions for the spatial variability of Reference Watershed, GLC Watershed, and consolidated gully are shown in (a2), (b2), and (c2), respectively.

Comparison of Soil Organic Carbon transport and transformation
The SOC stock change at each of the 2×2 m 2 grid box is the net of the surface SOC transport and SOC biogeochemical transformation. The surface SOC transport is due to lost through erosion or gained through deposition, hence referred as the lateral flux. The SOC biogeochemical transformation is due to the soil-atmosphere CO 2 exchange that includes indirect accumulation via plant residues or release to the atmosphere as CO 2 through microbial decomposition [2][3][4] , hence referred as the vertical flux. Figure S2a shows the stocks and flows of carbon through different reservoirs, which are air, plants, and soils. Figure S2b (expanded from Figure 4a in the main text) uses the data from Reference Watershed to show the lateral transport (x-axis) and vertical transformation (y-axis) relationship at each 2-D grid box. Within this domain, it can be divided into six groups, and each group has different legacies of SOC stock change: 1. x > 0, y > 0: ultimate sink (blue). SOC transport results in a net deposition and biogeochemical transformation results in a net carbon accumulation. 2.
4. x < 0, y < 0: ultimate source (red). SOC transport results in a net erosion and biogeochemical transformation results in a net carbon loss.
6. x > 0, y < 0, |y| > |x|: transformation dominated (yellow); SOC accumulation < decomposition; the site is a source to the atmospheric CO 2 . Figure S2: Illustration of the contribution of Soil Organic Carbon (SOC) transport and transformation at each grid point. a) Conceptual illustration of stocks and flows of carbon between atmosphere, soil, and plants. b) Illustration of six groups with different signs and relative magnitudes in the SOC transport and transformation dynamics.

Litter input estimation
The litter input is estimated from the Normalized Difference Vegetation Index (NDVI) ( Figure  S3a). The relationship between surface litter input and NDVI follows an exponential relationship. Figure S3: Normalized difference vegetation index (NDVI) and associated litter input. a) NDVI is processed from Landsat satellite bands for two years of record (2016-2017). The plot shows spatial average for the natural (or watersheds) area is spatially averaged for the natural area (a1) and consolidated gully land (a2), respectively. The black line is the mean for the 2 years record, which is then used through the simulation period repeatedly for each year. b) The estimated relationship between the above ground input of plant residues and NDVI.

Depositional zone and erosion zone
Within the PDF of the Reference watershed ( Figure S1(a2)), two zones (erosional and depositional zone) are chosen as shown in Figure S4. The spatial means of the fluxes of SOC transport and transformation are extracted to reveal the SOC intra-annual cycle under soil erosion and deposition conditions. The results are shown in Figure 5 Figure S4: Two zones from the probability density function of total soil depth changes. Each zone accounts for 5% of the total watershed area.

Land management practices in the consolidated gullies
Comparative study of set different scenarios which have different combinations of litter input and SOC mean residence time in the surface (5 cm) soil to test how SOC stock changes after a 50-yr evolution ( Figure S5).  Figure S6 shows the estimation of the initial SOC profile in the Gutun study sites. We collected twenty profiles which includes SOC content and bulk density. Among the twenty sites, eighteen sites are with natural plants and two sites are with corn field. Details of the empirical relationship between SOC concentration and soil depth can be found in the Methods Section. Figure S6: Illustration of the estimation of initial SOC profiles. Dots are from field sampling as shown in Figure S7a. Solid lines are fitted curve using the non-linear least square method for SOC(Z) = ae −bZ + c, where a, b, and c are parameters, and Z is the soil depth measured from the surface.

Surface SOC and forcing data
The spatial map of surface SOC and land cover are from the field survey in 2015. SOC is sampled at 0-10 cm and 10-20 cm at the locations as shown in Figure S7a. The surface values of SOC is based on the SOC profiles ( Figure S6) when the soil thickness is zero. The rainfall data with 10-yr record and 40-yr simulation are shown in the Figure S7c.  '30 N , 109 o 19'23 E), 44 km away on the northwest side of our study site (Gutun station), which is the closest available station. Soil organic matter parameters Carbon in fast (or litter) pool C l kg C/m 3 see footnote b Carbon in slow (or humus) pool C h kg C/m 3 see footnote b Carbon in biomass pool C b kg C/m 3 see footnote b a from surface to bottom of the initial 1 m soil depth for each layer. The distance among each other is the thickness for each layer. b The carbon concentration at each grid is a function of soil depth(Z) following the profile as shown in Figure S6.