Net greenhouse gas balance in U.S. croplands: How can soils be part of the climate solution?

Agricultural soils play a dual role in regulating the Earth's climate by releasing or sequestering carbon dioxide (CO2) in soil organic carbon (SOC) and emitting non‐CO2 greenhouse gases (GHGs) such as nitrous oxide (N2O) and methane (CH4). To understand how agricultural soils can play a role in climate solutions requires a comprehensive assessment of net soil GHG balance (i.e., sum of SOC‐sequestered CO2 and non‐CO2 GHG emissions) and the underlying controls. Herein, we used a model‐data integration approach to understand and quantify how natural and anthropogenic factors have affected the magnitude and spatiotemporal variations of the net soil GHG balance in U.S. croplands during 1960–2018. Specifically, we used the dynamic land ecosystem model for regional simulations and used field observations of SOC sequestration rates and N2O and CH4 emissions to calibrate, validate, and corroborate model simulations. Results show that U.S. agricultural soils sequestered 13.2±1.16$$ 13.2\pm 1.16 $$ Tg CO2‐C year−1 in SOC (at a depth of 3.5 m) during 1960–2018 and emitted 0.39±0.02$$ 0.39\pm 0.02 $$ Tg N2O‐N year−1 and 0.21±0.01$$ 0.21\pm 0.01 $$ Tg CH4‐C year−1, respectively. Based on the GWP100 metric (global warming potential on a 100‐year time horizon), the estimated national net GHG emission rate from agricultural soils was 122.3±11.46$$ 122.3\pm 11.46 $$ Tg CO2‐eq year−1, with the largest contribution from N2O emissions. The sequestered SOC offset ~28% of the climate‐warming effects resulting from non‐CO2 GHG emissions, and this offsetting effect increased over time. Increased nitrogen fertilizer use was the dominant factor contributing to the increase in net GHG emissions during 1960–2018, explaining ~47% of total changes. In contrast, reduced cropland area, the adoption of agricultural conservation practices (e.g., reduced tillage), and rising atmospheric CO2 levels attenuated net GHG emissions from U.S. croplands. Improving management practices to mitigate N2O emissions represents the biggest opportunity for achieving net‐zero emissions in U.S. croplands. Our study highlights the importance of concurrently quantifying SOC‐sequestered CO2 and non‐CO2 GHG emissions for developing effective agricultural climate change mitigation measures.


| INTRODUC TI ON
Contemporary agriculture is facing multiple challenges such as feeding a growing world population and mitigating climate change (Chen et al., 2014;Foley et al., 2011;Pittelkow et al., 2015).By 2050, global food production may need to increase by 25%-70% (Hunter et al., 2017), or even double current production levels (Tilman et al., 2011), to meet projected food demands.This would inevitably lead to substantial increases in greenhouse gas (GHG) emissions due to nitrogen (N) fertilizer use and cropland expansion (Cavigelli et al., 2012;Molotoks et al., 2018;Thompson et al., 2019;Zabel et al., 2019).To date, agriculture has been a major force in anthropogenic global warming, contributing about 25%-30% and 35%-50% of global land biogenic emissions of nitrous oxide (N 2 O) and methane (CH 4 ), respectively (Tian et al., 2016).This constitutes a great challenge to achieve the Paris climate goal of limiting global warming to well below 2°C by the end of this century (Tian, Xu, Canadell, et al., 2020).Therefore, reducing GHG emissions from the agricultural sector is an imminent need for mitigating climate change.Global croplands account for about 10% of the terrestrial soil organic carbon (SOC) stock (IPCC, 2019;Watson et al., 2000) and could potentially sequester 0.90-1.85Pg C year −1 in the top 0.3 m of soils, which is equivalent to 26%-53% of the soil carbon sequestration target of 3.5 Pg C year −1 established by the 4p1000 Initiative for climate mitigation (Zomer et al., 2017).Increasing SOC stock is considered to be the most important countermeasure for GHG mitigation in agriculture (Mosier et al., 2006;Smith et al., 2010).Besides sequestering atmospheric CO 2 , enhancing SOC stocks can also provide multiple co-benefits, such as reducing soil erosion, strengthening climate resilience, and improving soil fertility and health (Lal, 2018;Sohi, 2012).Thus, advancing our understanding of the magnitude and spatiotemporal variations of net GHG balance (i.e., sum of SOC sequestration of CO 2 and emissions of N 2 O and CH 4 ) in agricultural soils, as well as drivers of their change, is critical for measuring the cumulative radiative forcing of non-CO 2 GHG emissions and CO 2 uptake (Robertson & Grace, 2004) and developing effective agricultural climate change mitigation strategies.Meanwhile, improved understanding could contribute to the achievement of Sustainable Development Goals including "Climate Action" and "Zero Hunger".
As one of the most important agricultural producers in the world, U.S. agriculture contributed a significant portion of global agricultural GHG emissions.Numerous studies have measured and quantified N 2 O and CH 4 emissions and SOC sequestration in U.S. croplands (Del Grosso et al., 2010;EPA, 2021;Linquist et al., 2018;Lokupitiya et al., 2012;Lu et al., 2021;Ogle et al., 2010;Penman et al., 2000;Yu et al., 2018).For example, it is reported that agriculture emitted ~10% of the national total GHG emissions in 2019 in the U.S. and was the largest source of N 2 O emissions (~75%; EPA, 2021).Nonetheless, these estimates are highly uncertain and wide-ranging, largely due to differences in quantification methods and data sources (Ogle et al., 2010;Tian et al., 2018Tian et al., , 2019;;Xu et al., 2012).Specifically, uncertainty in SOC stock changes in U.S. croplands could range from − 4.6 Tg C year −1 to + 4.9 Tg C year −1 (Ogle et al., 2006(Ogle et al., , 2010)), and uncertainty in U.S. agricultural soil N 2 O emissions could range from− 0.07 Tg N year −1 to + 0.1 Tg N year −1 (Cavigelli et al., 2012;Del Grosso et al., 2010).In addition, most of these studies have focused on estimating either individual GHG fluxes or SOC sequestration rate, while much less work has been done in quantifying magnitude and spatiotemporal variations in net soil GHG balance in U.S. croplands (EPA, 2021).Due to possible trade-offs between SOC sequestration and GHG emissions under different agricultural management practices (Guenet et al., 2021;Tian et al., 2011Tian et al., , 2015)), simultaneous quantification of SOC-sequestered CO 2 and non-CO 2 GHG emissions is crucial to accurately assess the overall climate abatement potential of mitigation measures.Furthermore, whether SOC sequestration of CO 2 in U.S. croplands can offset non-CO 2 net GHG emission rate from agricultural soils was 122.3 ± 11.46 Tg CO 2 -eq year −1 , with the largest contribution from N 2 O emissions.The sequestered SOC offset ~28% of the climate-warming effects resulting from non-CO 2 GHG emissions, and this offsetting effect increased over time.Increased nitrogen fertilizer use was the dominant factor contributing to the increase in net GHG emissions during 1960-2018, explaining ~47% of total changes.In contrast, reduced cropland area, the adoption of agricultural conservation practices (e.g., reduced tillage), and rising atmospheric CO 2 levels attenuated net GHG emissions from U.S. croplands.Improving management practices to mitigate N 2 O emissions represents the biggest opportunity for achieving net-zero emissions in U.S. croplands.Our study highlights the importance of concurrently quantifying SOC-sequestered CO 2 and non-CO 2 GHG emissions for developing effective agricultural climate change mitigation measures.GHG emissions and how far we are from achieving carbon-neutral agriculture remains unclear.
Climate-smart agriculture (CSA) management practices (e.g., reduced tillage, optimized N fertilizer use, and alternate wetting and drying irrigation) have been advocated to reduce GHG emissions without compromising crop yield (FAO, 2013;Miralles-Wilhelm, 2021).Various field investigations and meta-analyses have explored the effects and efficacy of these practices (Bai et al., 2019;Gerber et al., 2016;Shang et al., 2021;Sun et al., 2020).However, most existing work assessing the impact of CSA measures on GHG emissions has focused on a single management practice and one or two GHG fluxes (e.g., CO 2 or N 2 O) at a time (Huang et al., 2022;Lu et al., 2022;Yu et al., 2020).Relatively few studies have simultaneously quantified the integrated effects of multiple management practices on net soil GHG balance, especially at large spatial scales (e.g., national and continental scales).Notably, EPA (2013) and Tian et al. (2016) comprehensively quantified all three major GHG emissions and the abatement potential of non-CO 2 GHGs on both national and global scales.Considering that some CSA management practices may have antagonistic effects on SOC sequestration and non-CO 2 GHG emissions (Guenet et al., 2021), and the resulting effects of different practices typically have large variations and may be non-additive (Yue et al., 2019), studies that fail to combine SOC sequestration and non-CO 2 GHG emissions (as well as multiple practices together) may lead to inconsistencies when making comparisons that would not provide effective assessments (Shang et al., 2021).
Global environmental changes such as climate change, elevated atmospheric CO 2 concentration, and N deposition have also substantially affected agricultural GHG emissions (Ren et al., 2011(Ren et al., , 2020)).These factors vary over space and time in a highly heterogeneous geographical environment (e.g., diverse soil types and cropping systems) that can affect the effectiveness of CSA practices (Abdalla et al., 2013;Sun et al., 2020).This means a mitigation practice that is effective in one location or under certain conditions may not be effective elsewhere or under other conditions (Shang et al., 2021).In an example illustrating the importance of considering interactions between environmental factors and agricultural management practices, Huang et al. (2018) found that conversion from conventional tillage to no-tillage reduced GHG emissions in dry but not in humid climates.However, relatively few studies have quantitatively attributed changes in the net soil GHG balance of U.S. croplands to different drivers (including multiple management practices and environmental factors) over long-term periods (Moore et al., 2023), although such factorial contribution analyses are essential for accurately assessing impacts of these CSA practices and developing effective climate mitigation measures.
Field experiments provide feasible and reliable means of elucidating complex relationships of agricultural management practices and net soil GHG balance under multiple environmental changes (Plaza-Bonilla et al., 2018).However, directly extrapolating sitespecific findings to large spatial areas is difficult due to unique environmental and management conditions of each site (Huang et al., 2022).Process-based terrestrial biosphere models (TBMs), with well-represented crop growth processes and agricultural management practices (e.g., N fertilization, tillage, irrigation, and rotation), as well as detailed hydrological, biophysical, and biogeochemical processes, can account for effects of spatial and temporal variations in environmental and management conditions on net soil GHG balance at large scales (Bondeau et al., 2007;McDermid et al., 2017;You et al., 2022).However, the simulation performance of TBMs is largely limited by the availability of high-quality model forcing datasets (i.e., introducing uncertainties in input data), the lack of sufficient data for model calibration and validation (i.e., introducing uncertainties in model parameters), and the inadequate representation of relevant processes in the model (i.e., introducing uncertainties in model structures; Fisher et al., 2014;Gurung et al., 2020;Ogle et al., 2010).In view of respective strengths and weaknesses of field observations and TBMs simulations, the integration of model and data would provide promising means to overcome these bottlenecks (Fer et al., 2021;Peng et al., 2011).
In this study, we quantified the combined effects of multiple management practices and environmental changes on the magnitude and spatiotemporal variations of net soil GHG balance in U.S. croplands using a model-data integration approach.The model used here is the Dynamic Land Ecosystem Model v4.0 (DLEM v4.0), which is a highly integrated process-based TBM.DLEM v4.0 is capable of simultaneously depicting biosphere-atmosphere exchanges of CO 2 , N 2 O, and CH 4 , driven by multiple environmental forcings and management factors across site, regional, and global scales (Pan et al., 2021;Tian, Xu, Canadell, et al., 2010;Tian, Xu, Pan, et al., 2020;Yao et al., 2020;You et al., 2022).High-resolution model forcing datasets were developed to drive the DLEM.Field observations of SOC sequestration rates and non-CO 2 GHG emissions under various management practices and environmental conditions on U.S. croplands were compiled to calibrate, validate, and corroborate model simulations.The objectives of this work were (1) to estimate the net soil GHG balance of U.S. croplands as driven by changes in multiple management practices (e.g., N fertilization, tillage, and irrigation), climate conditions, historical land use, atmospheric CO 2 concentration, and N deposition spanning from 1960 to 2018, (2) to examine the relative contributions of SOC sequestration of CO 2 and non-CO 2 GHG emissions to the net soil GHG balance of U.S. croplands, and (3) to quantify factorial contributions of different drivers to the spatial and temporal variations in net soil GHG balance across the country.

| Metadata collection
A comprehensive literature search was conducted to identify peerreviewed publications reporting in-situ soil GHG emissions from U.S. croplands using several databases including Google Scholar, Web of Science, and Scopus.Search keywords included "cropland or crop or corn or maize or soybean or wheat or rice", "the United States or America or U.S. or USA", "soil organic carbon or SOC", "nitrous oxide or N 2 O", "methane or CH 4 ", and/or "greenhouse gases or GHG".To ensure the quality of compiled datasets, papers identified were further refined by the following criteria: (1) measurements were made in the field rather than in the laboratory; (2) ancillary information such as cropping systems, experimental year and duration, and applied management practices (e.g., N fertilizer use rate, tillage type, and irrigation) were provided; and (3) replicated field experiments were performed.
This search identified a total of 576 site-years of data representing 79 locations from 91 peer-reviewed publications (Figure 1 and

| Model descriptions
DLEM v4.0 is a highly integrated TBM that couples major biophysical, biogeochemical, and hydrological processes to quantify daily, spatially explicit carbon, water, and nutrient stocks and fluxes in terrestrial ecosystems and inland water systems at site, regional, and global scales (Pan et al., 2021;Tian, Xu, Canadell, et al., 2010;Tian, Xu, Pan, et al., 2020;You et al., 2022).The simulation of terrestrial carbon, water, and nutrient dynamics is driven by multiple environmental forcings (e.g., climate change, atmospheric CO 2 concentration, and N deposition) and various management factors (e.g., N fertilizer use rate, irrigation, and manure application).Nutrient loading and export from land to rivers and oceans are coupled with processes of leaching, runoff, and sedimentation (Bian et al., 2023;Liu et al., 2013;Tao et al., 2014).These nutrients are then used as inputs for the aquatic module of the DLEM to simulate aquatic biogeochemical processes (Tian, Xu, Pan, et al., 2020).For instance, simulated N loads (including nitrate, ammonium, dissolved organic N, dissolved inorganic N, and particulate organic N) serve as substrates for nitrification and denitrification processes within the aquatic module (Yao et al., 2020).Subsequently, these N loads are either released into the atmosphere as N 2 O or transported to the estuary.
To meet cross-scale agricultural application needs (e.g., management guidance, agricultural climate change mitigation and adaptation), DLEM v4.0 also includes mechanistic representations of dynamic crop growth and development processes, such as cropspecific phenological development, carbon allocation, yield formation, and biological N fixation (You et al., 2022).Additionally, agricultural management practices such as N fertilizer use, irrigation, tillage, manure application, dynamic crop rotation, cover cropping, and genetic improvements are also included.Regarding tillage practices, three aspects of tillage impacts on agroecosystems are represented in DLEM v4.0, including (1) changes in surface residue coverage and subsequent redistribution of soil organic matter and nutrients within tilled soil layers due to tillage mixing, (2) changes in litter interception, bulk density, soil moisture, and other waterrelated effects on processes such as nitrification, denitrification, and leaching, and (3) changes in the soil decomposition rate (You et al., 2022).The vertical soil profile in DLEM is discretized into 10 layers down to a depth of 3.5 m, in which the layer thickness increases  Dataset (Brown & Pervez, 2014;Pervez & Brown, 2010) as a base map.The annual manure N application dataset from 1860 to 2018 was acquired from Bian et al. (2021).The state-level earliest and latest crop planting dates were obtained from the NASS survey report (NASS, 2010), which provides planting and harvesting windows in most historical years.

Land use and land cover change (LULC):
We developed a spatially explicit annual LULC dataset at a spatial resolution of 1 × 1 km 2 over the contiguous U.S. during 1630-2020 using machine learning and geospatial modeling approaches (Li et al., 2023).
Multi-source datasets such as satellite-derived land cover maps, national inventories, topographical data, and model-based land use change data were used for the reconstruction.Furthermore, vegetation type datasets during 1860-2020 were derived from Chen et al. (2006).
3. Natural environmental changes (climate conditions, atmospheric CO 2 , and N deposition): The historical daily climate dataset (including precipitation, solar radiation, maximum, minimum and mean temperatures) from 1860 to 2018 was reconstructed from the North American Land Data Assimilation System product (Mitchell et al., 2004;Xia et al., 2012), the Climate Research Unit-National Centers for Environmental Prediction dataset (Mitchell & Jones, 2005), and the IPSL Climate Model dataset (Boucher et al., 2020) using the delta downscaling method (Liu et al., 2013).4. Soil properties and other auxiliary information: Soil physical and chemical properties were obtained from the ISRIC-WISE Harmonized Global Soil Profile dataset (Batjes, 2008).Other auxiliary information such as topography and river network was obtained from our previous studies (Tian, Xu, Canadell, et al., 2010;Tian, Xu, Pan, et al., 2020).

| Model calibration, validation, and sensitivity analysis
DLEM has been widely validated and applied to estimate N 2 O and CH 4 emissions and SOC stocks at multiple sites and in large-scale regions including the U.S. (Huang et al., 2020;Lu et al., 2021;Tian et al., 2012;Yu et al., 2018), North America (Tian, Xu, et al., 2010;Xu et al., 2012), China (Ren et al., 2011;Zhang et al., 2020), and across the globe (Friedlingstein et al., 2020;Ren et al., 2020;Saunois et al., 2020;Tian, Xu, Canadell, et al., 2020).In this study, we rigorously calibrated and validated DLEM, as driven by the forcing datasets developed in Section 2.3, to better simulate SOC stock, and N 2 O and CH 4 emissions in U.S. croplands using field observations compiled by the metadata collection described in Section 2.1.We calculated SOC sequestration rates as the differences in SOC stocks between two adjacent years.The changes in national SOC-sequestered CO 2 observed in this study can be attributed to changes in cropland area, natural environmental factors, and agricultural management activities.Additionally, since the collected measurements of SOC sequestration rates and SOC stocks were reported at various soil depths (e.g., 0-20 cm or 0-100 cm), we outputted simulation results at the corresponding soil depths for validation.Several metrics were used to quantify model performance, including coefficient of determination (R 2 ), root mean square error (RMSE), and normalized root mean square error (NRMSE).
In total, 576 site-year measurements representing 79 U.S. cropland sites covering major cropping systems were used to calibrate, validate, and corroborate model simulations (Figure 1).
Values of major parameters related to N 2 O, CH 4 , and SOC processes were determined through model calibration within a reasonable range of values reported in the literature (Table S2).
Specifically, we first used default parameter values to run the model, and then we manually tuned the parameters within the reported ranges to obtain a close match between observed and simulated values.During the calibration process, we utilized 10, 6, 4, and 6 site-year datasets for N 2 O, CH 4 , SOC sequestration rate, and SOC stock, respectively.We adopted the parameter set that obtained the minimal bias between the simulated and measured values across all calibration datasets as the optimal parameters and used it for the regional simulation.Additionally, apart from calibrating parameters related to N 2 O, CH 4 , and SOC processes, we calibrated the model with a focus on crop yields using data  .91, .46, and .64,respectively (Figure 2).Furthermore, given the significant differences in CH 4 emissions between rice and non-rice crops, we also calculated their individual R 2 values, which were .75 for rice and .54for non-rice crops.
We also quantified the sensitivity of simulated regional SOC sequestration rate and N 2 O and CH 4 fluxes in U.S. croplands to variations in model parameters (Tian et al., 2011).Specifically, we first conducted a variance-based global sensitivity analysis for each major crop type to quantify the relative importance of model parameters in simulating SOC sequestration rate and N 2 O and CH 4 emissions using the Sobol' method (Sobol, 1993).We then identified parameters having a significant impact on simulated SOC sequestration rate and N 2 O and CH 4 fluxes (Figures S1-S5).Next, we used the Monte Carlo sampling scheme to generate an ensemble of 100 sets of these sensitive parameters by randomly varying their values within a 20% range of their calibrated values based on their respective probability distribution functions (Tian et al., 2011;You et al., 2022).Finally, we used the generated parameter sets as inputs for DLEM to simulate regional SOC sequestration rate and N 2 O and CH 4 emissions from U.S. croplands.The sensitivity of simulated SOC sequestration rate and N 2 O and CH 4 fluxes to model parameter variations was expressed as ± 1 standard deviations (SDs) derived from these simulations.

| Model implementation and experimental design
Implementation of DLEM v4.0 included three major steps: an equilibrium run, a spin-up run, and a transient run.The equilibrium run was driven by average annual climate data during the 1860s and other factors (e.g., vegetation type, atmospheric CO 2 concentration, and N deposition) in 1860.The equilibrium state was assumed to be reached when changes in carbon, N, and water pools between two consecutive 20-year periods were less than 0.5 g C m −2 year −1 , 0.5 g N m −2 year −1 , and 0.5 mm year −1 , respectively.LULC turned on and off (Lu et al., 2021).Thus, the factorial contribution of LULC was calculated as the difference between S9 and S10 (Table 1).In addition, because our analysis focused on the period 1960-2018, we calculated the factorial contribution of each factor relative to the average state in the 1950s.That means, the factorial contribution of each factor from 1960 to 2018 was calculated by subtracting the average state of the 1950s from its original factorial contribution.

| Global warming potential calculation
The global warming potential (GWP) is an index to measure the integrated radiative forcing from the emission of 1 kg of a trace gas relative to that of CO 2 (Myhre et al., 2013).In GWP conversions, CO 2 is typically considered the reference gas with a GWP constant of 1.Note that the SOC sequestration rates and SOC stocks here were reported by studies at various soil depths (e.g., 0-20 cm or 0-30 cm), thus we outputted simulation results at corresponding soil depths for validation.SOC sequestration rates represent SOC stock changes and are calculated as the differences in SOC stocks between two adjacent years, in which a negative SOC sequestration rate represents a decrease in SOC stock (indicating decomposition), and a positive SOC sequestration rate represents an increase in SOC stock (indicating sequestration)."Others" represent all crop types not listed in the legend (e.g., barley, cotton, and sorghum).Additionally, due to the small amount of N 2 O and CH 4 emissions (close to 0) at many sites, these datapoints are stacked near the origin, resulting in the number of datapoints on the graph appearing to be fewer than the actual number of datapoints involved in the statistical analysis.
obtain a comprehensive assessment of the climatic impact of net soil GHG balance, we adopted the following equation to calculate the combined GWPs for SOC sequestration of CO 2 and N 2 O and CH 4 emissions: where in terms of their CO 2 equivalents, and this study used the GWP values integrated over a time horizon of 100 years for CO 2 , N 2 O, and CH 4 , which were 1, 265, and 28, respectively (Myhre et al., 2013).

| National budget and dynamics of net GHG balance in U.S. croplands
Our simulations showed that U.S. croplands acted as a net carbon sink during 1960-2018 with an average SOC sequestration rate of 13.2 ± 1.16 Tg C year −1 (at a depth of 3.5 m), and acted as a net source of N 2 O and CH 4 with average emission rates of 0.39 ± 0.02 Tg N year −1 and 0.21 ± 0.01 Tg C year −1 , respectively (Figure 3).Both SOC sequestration and N 2 O and CH 4 fluxes in U.S. croplands exhibited large interannual variations during 1960-2018, but showed overall significant increasing trends (with respective rates of 0.429 Tg C year −2 , 0.003 Tg N year −2 , and 0.001 Tg C year −2 ).
Using the GWP100 metric, sequestered SOC in U.S. agricultural soils reduced national net GHG balance at an average rate of 48.4 ± 4.25 Tg CO 2 -eq year −1 during 1960-2018, whereas N 2 O and CH 4 emissions contributed to the net GHG balance at average rates of 162.76 ± 8.33 Tg CO 2 -eq year −1 and 8.01 ± 0.37 Tg CO 2 -eq year −1 , respectively (Figure 4).Thus, non-CO 2 GHG emissions (i.e., sum of N 2 O and CH 4 emissions) from U.S. croplands surpassed SOC sequestered, indicating that U.S. croplands acted as a net source of GHGs.Statistically, sequestered SOC offset ~28% of climate-warming effects resulting from non-CO 2 GHG emissions during 1960-2018, and the proportion of the offset increased over time.When considering both SOC sequestration and non-CO 2 GHG emissions, the average net GHG balance during 1960-2018 was estimated to be a GHG source of 122.3 ± 11.46 Tg CO 2 -eq year −1 and exhibited a substantial decadal variability, increasing from 98.32 ± 7.13 Tg CO 2 -eq year −1 in the 1960s to 155.28 ± 11.48 Tg CO 2 -eq year −1 in the 1980s and gradually decreasing to 115.87 ± 16.18 Tg CO 2 -eq year −1 in the 2010s. (1)

| Spatial patterns of net GHG balance in U.S. croplands
With an average SOC sequestration rate of ~10.5 g C m −2 year −1 , simulation results over the study period indicated that most U.S. croplands acted as carbon sinks, with the Midwest, Southeast, and Northwest regions having relatively high SOC sequestration rates (Figure 5).The spatial pattern of N 2 O emissions varied substantially across the country, with hotpots in the Midwest region having peak N 2 O emission rates as high as 0.8 g N m −2 year −1 .In contrast, the distribution of CH 4 flux was polarized, with a high CH 4 emission rate of ~13 g C m −2 year −1 in the Mississippi Delta and the Sacramento Valley regions due to rice cultivation and rates approaching zero in other areas.
When simultaneously taking SOC sequestration and non-CO 2 GHG emissions into account, we found that the distribution of net soil GHG balance showed large spatial heterogeneity, with hotspots in the Midwest and Mississippi Delta regions where peak net soil GHG emissions were estimated to be higher than 2.5 Mg CO 2 -eq ha −1 year −1 (Figure 6).In contrast, some U.S.

| Factorial contributions of multi-driver changes to net GHG balance in U.S. croplands
We further quantified the factorial contributions of key drivers, experiments (Table 1).Our results revealed that increased use of N fertilizer was the dominant factor driving changes in the net GHG balance of U.S. croplands in comparison to the average state in the 1950s.This increase contributed to a net GHG balance increase of 79.9 Tg CO 2 -eq year −1 and roughly explained 47% of the total changes (Figure 8).Increased atmospheric N deposition and manure application also contributed to the increase in net GHG balance, with average rates of 8.3 Tg CO 2 -eq year −1 and 2.7 Tg CO 2 -eq year −1 , respectively.Conversely, LULC resulted in a substantial reduction in the net GHG balance in U.S. croplands, with an average mitigation rate of 42.2 Tg CO 2 -eq year −1 , which explained ~23% of the total changes.Rising atmospheric CO 2 concentration was the secondlargest mitigator, reducing the net GHG balance at an average rate of 20.1 Tg CO 2 -eq year −1 and accounting for ~9% of the total changes.
Reduced tillage and increased irrigated area were also effective in mitigating the net GHG balance in U.S. croplands, with average mitigation rates of 3.1 Tg CO 2 -eq year −1 and 1.9 Tg CO 2 -eq year −1 , respectively.Additionally, compared to the 1950s, climate change initially reduced net GHG balance in U.S. croplands by an average rate of 12.5 Tg CO 2 -eq year −1 from the 1960s to the 1990s, but later increased net GHG balance at an average rate of 36.5 Tg CO 2 -eq year −1 .

| Comparison with previous studies
We compared our estimates of SOC sequestration rate and N 2 O and CH 4 emissions in U.S. croplands with other regional estimates published (Table 2).Over the past six decades, our estimated SOC sequestration rate in U.S. croplands ranged from 4.9 ± 0.66 Tg C year −1 in the 1960s to 22.2 ± 1.91 Tg C year −1 in the 2010s and fell within the rate change reported by others (Ogle et al., 2006(Ogle et al., , 2010;;Pacala et al., 2001;Sleeter et al., 2018;West et al., 2008;Zhang et al., 2015).crop residues would have also significantly increased during the same period.Although the DLEM have accounted for the removal of a certain percentage of crop residues from the agroecosystem through tillage practices (50%, 25%, and 0% for conventional tillage, reduced tillage, and no tillage, respectively), we acknowledge that large uncertainties still exist concerning these residue removal percentages, potentially biasing the simulated residue inputs into soils and SOC sequestration rate.In terms of N 2 O emissions, our estimates ranged from 0.27 ± 0.02 Tg N year −1 in the 1960s to 0.45 ± 0.03 Tg N year −1 in the 2010s, which were also comparable to previous studies (Chen et al., 2016;Del Grosso et al., 2006;EPA, 2018;Griffis et al., 2013;Lu et al., 2021;Mummey et al., 1998;Tian et al., 2019).Tg C year −1 (EPA, 2015; Sass et al., 1999;Tian et al., 2015).
Overall, our individual estimates of SOC sequestration rate and N 2 O and CH 4 emissions in U.S. croplands were in similar magnitudes and variation ranges as other regional estimates, but differences still exist, possibly due to uncertainties in forcing data and differences in estimation methods.However, unlike previous studies that focused solely on individual GHGs or SOC, our work quantified both

| Impacts of agricultural management factors and environmental changes on net GHG balance
Leveraging agricultural management practices to curb net GHG emissions from croplands has recently come under sharp focus due to their large mitigation potential, low cost and accompanying co-benefits such as improved soil and water quality and biodiversity maintenance (Fargione et al., 2018).As seen in other studies (Christopher & Lal, 2007;Gerber et al., 2016;Lu et al., 2021), our factorial analysis indicated that N fertilization (including both synthetic N fertilizer use and manure application) contributed significantly to net GHG emissions from U.S. croplands (Figure 8).studies (Huang et al., 2022;Lu et al., 2022;Yu et al., 2020).For instance, Yu et al. (2020) found that reduced tillage intensity in U.S.  1960-2018 (Figures S7e andS8e).Irrigation possibly increased aboveground and belowground biomass and led to higher soil carbon inputs and SOC content (Bai et al., 2019;Blanco-Canqui et al., 2011).Additionally, we found that LULC reduced net GHG emissions from U.S. croplands (Figure 8) at an average mitigation rate of 42.2 Tg CO 2 -eq year −1 , primarily due to the increased national total amount of SOC-sequestered CO 2 in recent decades compared to the average state in the 1950s (Figure S6).Specifically, when compared with the 1950s, there was a decrease in cropland area during recent decades (Figure S8b), coinciding with a notable increase in the SOC sequestration rate.However, it is worth noting that the national total amount of SOC-sequestered CO 2 is determined by the multiplication of cropland area and the average SOC sequestration rate.Consequently, the recent reduction in cropland area was counterbalanced by the concurrent increase in SOC sequestration rate, to the extent that the total amount of CO 2 sequestrated in U.S.
croplands during recent decades still surpassed that of the 1950s.
This dynamic interplay has ultimately led to a reduction in net GHG Changes in environmental factors, including elevated atmospheric N deposition and climatic conditions over the past two decades, also contributed significantly to net GHG emissions in U.S. croplands, while rising atmospheric CO 2 concentration served to mitigate GHG emissions (Figure 8).(Carey et al., 2016;Pärn et al., 2018;Weier et al., 1993;Yvon-Durocher et al., 2014).( Xu et al., 2020;Yang et al., 2021), where increased N availability can promote nitrification and denitrification processes and thereby N 2 O emissions.Similar to findings of others (Ren et al., 2011;van Groenigen et al., 2011;Xu et al., 2020), our study showed a negative response to rising atmospheric CO 2 levels, which can be ascribed to increased soil carbon inputs from CO 2 enhanced crop biomass production (van Groenigen et al., 2011).

| Limitation and future work
Although our estimates of the net GHG balance are generally con- Accurately quantifying the uncertainty associated with model structures is a challenging task.Some studies have proposed using Bayesian methods to address this issue through conditioning model behavior on measurements (Engeland et al., 2005;Gurung et al., 2020).A notable example is the work of Gurung ventions (Fontaine et al., 2007), there is a pressing need for future research to intensify measurements on deep soil carbon and its dynamics.Such efforts will provide a more robust foundation for modeling deep soil carbon dynamics accurately.
Finally, our current estimates only include direct GHG emissions from agricultural soils, while neglecting indirect emissions due to N leaching in agroecosystems (e.g., stimulated N 2 O and CH 4 emissions in aquatic systems).Considering the high anthropogenic nutrient loading in stream and river networks caused by N fertilizer use (Tian, , Xu, Pan, et al. 2020;Yao et al., 2020), future work could quantify these indirect GHG emissions due to carbon and nutrient leaching to complement our current net GHG budget.We will address these limitations in future studies to provide an overall estimate of both direct and indirect emissions from U.S. croplands.

| CON CLUS ION
This study quantified the magnitude and spatiotemporal variations of the net soil GHG balance in U.S. croplands from 1960 to 2018 using a model-data integration approach.We found that , DLEM, greenhouse gas balance, methane, nitrous oxide, SOC sequestration, U.S. croplands | 3 of 20 YOU et al.
geometrically from top to bottom with values of 0.05, 0.05, 0.1, 0.2, 0.2, 0.3, 0.3, 0.5, 0.8, and 1 m, respectively.The SOC is distributed F I G U R E 1 Spatial distribution of field sites included in this study.Different colors represent different crop types, in which "others" represent crop types not listed in the legend (e.g., barley, cotton, and sorghum).Point size represents the number of replications/ observations at each site.| 5 of 20 YOU et al. belowground following an exponential function that characterizes the specific root density profile for each plant functional type.The well-represented crop growth processes and management practices enable DLEM v4.0 to simulate crop state variables (e.g., leaf area index and tissue biomass) across growth stages and biogeochemical fluxes and pools of carbon, N, and water related to agroecosystems (e.g., crop yield, GHG emissions, SOC, and nutrient leaching) across various spatial and temporal scales.More details about the representation of crop growth processes and different agricultural management practices in DLEM v4.0 are provided in You et al. (2022).
Four long-term datasets distributed over the U.S. at 5 × 5 arc-min grid spacing were developed to provide forcing conditions for DLEM v4.0: 1. Agricultural management practices (N fertilizer use rate, crop rotation, tillage, irrigation, and manure application): The annual crop-specific N fertilizer use rate dataset from 1910 to 2018 was reconstructed using the state-level N fertilizer use rates from USDA-NASS and the national-level commercial N fertilizer consumption data from Mehring et al. (1957) and USDA-ERS (2019) following Cao et al. (2018).The annual crop rotation dataset from 1910 to 2018 was developed by combining the USDA Cropland Data Layer (CDL) product and National Agricultural Statistics Service (NASS) survey of countylevel crop planting areas using the spatialization method of Yu et al. (2018).The annual tillage intensity map from 1960 to 2018 was reconstructed from the county-level tillage practices survey data obtained from the National Crop Residue Management Survey (CRM) of the Conservation Technology Information Center (https:// www.ctic.org/ CRM) following You et al. (2022).Tillage maps for missing years were kept the same as the nearest years with available data.The original five tillage practices in the CRM dataset were reorganized into four types by combining ridge and mulch tillage types to conservation tillage.The county-level CRM dataset was combined with the CDL-derived crop rotation map and the USDA-NASS crop planting area to estimate historical spatial distributions of tillage practices.The annual crop-specific irrigation dataset from 1950 to 2018 was downscaled from the county-level irrigation reanalysis (McManamay et al., 2021) and NASS irrigated cropland area survey, using the MODIS Irrigated Agriculture 2 concentration variations during 1860-2018 were from the NOAA GLOBALVIEW-CO 2 dataset derived from atmospheric and ice core measurements (www.esrl.noaa.gov).Monthly atmospheric N deposition variations during 1860-2018 were acquired from the International Global Atmospheric Chemistry (IGAC)/Stratospheric Processes and Their Role in Climate (SPARC) Chemistry-Climate Model Initiative (CCMI) (Eyring et al., 2013).
collected from the AmeriFlux Network, the Resilient Economic Agricultural Practices Project, and the United States Department of Agriculture-National Agricultural Statistics Service (further details can be found in You et al., 2022).After model calibration, field observed N 2 O, CH 4 , and SOC data (excluding data used for model calibration) were utilized to evaluate model performance.Overall, DLEM can well-simulate emissions of N 2 O and CH 4 , SOC sequestration rates, and SOC stocks compared with field observations from the metadata collection, in which RMSE (NRMSE) values were 0.16 g N m −2 (9.6%), 1.5 g C m −2 (4.4%), 156.9 g C m −2 (19.3%), and 1929.1 g C m −2 (17.6%), respectively, and the R 2 values were .6, The spin-up run was driven by detrended climate data during the 1860s to eliminate fluctuations due to the transition from equilibrium run to transient run.Finally, the transient run was driven by historical data from 1860 to 2018.For model forcing data not spanning the entire 1860-2018 period, we kept the data for the missing years the same as those of the nearest available years.For example, during the transient run, tillage data before 1989 and after 2011 were assumed to be constant at the levels of 1989 and 2011, respectively.The simulation results from the transient run are considered to represent actual changes in the real world due to the comprehensive inclusion of changes in all driving factors (e.g., LULC, climate conditions, and N fertilizer use).In this study, we will analyze these results from 1960 to 2018, and the model run between 1860 and 1959 serves as a spin-up for the slow soil biogeochemical cycles.We designed 11 simulation experiments to distinguish factorial contributions of different drivers to spatial and temporal variations in the net soil GHG balance of U.S. croplands (Table1).Attribution factors includedN fertilization, tillage, irrigation, manure application, climate change, atmospheric CO 2 concentration and N deposition, and LULC.A reference run (S0) was performed by keeping all factors at the 1860 level to examine model fluctuations resulting from internal system dynamics.The S0 run yielded background emissions with little human perturbation, utilizing the vegetation type map derived from the LULC dataset.An allcombined run (S1) was implemented by driving the model using all historically varying input forcings during 1860-2018 to represent the "best estimates" of SOC sequestration rate and N 2 O and CH 4 emissions from U.S. croplands.Net changes in SOC sequestration rate and N 2 O and CH 4 emissions driven by all factors were calculated as the difference between S1 and S0 simulations.Meanwhile, we performed seven additional simulations (S2-S8) to investigate individual contributions of changes in each factor to annual variations in SOC sequestration rate and N 2 O and CH 4 fluxes.Specifically, in each simulation one particular factor was kept constant at the 1860 level, while all other factors were set to vary over time, and the factorial contribution of this fixed factor was obtained by subtracting the simulation from the "allcombined" simulation (S1).Since LULC is usually accompanied by changes in the total input of management practices (e.g., manure and mineral fertilizer application), we calculated the factorial contribution of LULC by keeping all management factors constant at the 1860 level while varying other environmental factors with CH 4 and N 2 O emissions can be converted to 'CO 2 -equivalents' based on their respective GWP constants over a specified time horizon.To F I G U R E 2 Site-scale comparisons of model estimates and field observations of N 2 O (a), CH 4 (b), SOC sequestration rate (c), and SOC stock (d) across different crop types.Dashed line is the regression of observed data and modeled results, and the solid line is the 1:1 line.

F
CO 2 −C = − SOCSR, TA B L E 1 Factorial experiments to quantify the relative contributions of different drivers to changes in soil organic carbon sequestration rate and nitrous oxide and methane emissions from U.S. croplands.

3. 3 |
Relative contributions of SOC sequestration and N 2 O and CH 4 emissions to net GHG balance Given that non-CO 2 GHG emissions have surpassed the amount of SOC sequestered in U.S. croplands (Figure 4), we conducted a more comprehensive analysis of the spatial distribution of the relative contribution of SOC sequestration and N 2 O and CH 4 emissions to the net GHG balance of U.S. croplands (Figure 7).Over the study period, soil N 2 O emissions played a dominant role in controlling the net GHG balance of most croplands (e.g., the Midwest, Northern and Southern Plains hubs), followed by SOC sequestration, while CH 4 emissions only controlled the net GHG balance in the Mississippi Delta and Sacramento Valley regions, which are major rice growing areas.Meanwhile, the proportion of U.S. croplands dominated by SOC sequestration increased over time, indicating an increasing role of SOC sequestration in controlling the net GHG balance across the country.For example, most Midwest croplands were dominated by N 2 O emissions (red color) in the 1960s, but was controlled jointly by N 2 O emission and SOC sequestration (yellow and green colors) in the 2010s.
including multiple agricultural management practices and environmental forcings, to changes in the net soil GHG balance of U.S. croplands from 1960 to 2018 by setting up a series of simulation F I G U R E 3 Temporal variations in national SOC sequestration rate (a) and fluxes of N 2 O (b) and CH 4 (c) in U.S. agricultural soils from 1960 to 2018.Shaded area denotes the sensitivity of SOC sequestration rate and N 2 O and CH 4 fluxes to model parameter variations (represented by ±1SD).
For instance, Pacala et al. (2001) estimated the range of SOC sequestration rate in U.S. croplands to be 0-40 Tg C year −1 , Ogle et al. (2010) estimated a net increase in SOC of 17.5 ± 2.9 Tg C year −1 during 1995-2000, and Sleeter et al. (2018) estimated a net accumulation of SOC in U.S. croplands of 32 Tg C year −1 during 1973-2010.Our estimated SOC sequestration rate was generally consistent with these studies.In addition, our estimated SOC sequestration rate showed an increasing trend during 1960-2018, partly due to the increasing atmospheric CO 2 concentration (Walker et al., 2021), improved crop management practices (e.g., N fertilizer use and irrigation; Christopher & Lal, 2007), and advancements in crop technologies (e.g., genetics and breeding).These factors have resulted in significant increases in biomass production and grain yield, leading to greater crop residue inputs into soils.For instance, productions of major crops such as corn and soybean have nearly tripled to quadrupled in the U.S. over the past several decades (USDA, 2018).Despite the carbon allocation process in DLEM is driven by photosynthetic carbon supply and subject to multiple environmental stresses (You et al., 2022), the simulated ratio of yield to whole plant biomass has remained relatively stable over time, suggesting that the simulated F I G U R E 4 Temporal variations in national net greenhouse gas balance of U.S. croplands from the 1960s to the 2010s.Note that error bars represent ±1SD of net greenhouse gas balance in each decade.F I G U R E 5 Spatial patterns of average annual SOC sequestration rate (a) and fluxes of N 2 O (b) and CH 4 (c) in U.S. agricultural soils from 1960 to 2018.Note that negative values in soil fluxes represent uptake and positive values represent release.Therefore, negative soil CO 2 flux indicates SOC sequestration.| 11 of 20 YOU et al.
For example, the direct N 2 O emissions from U.S. agricultural soils estimated by the Global N 2 O Model Inter-comparison Project (NMIP) ranged from 0.3 Tg N year −1 to 0.62 Tg N year −1 during 2007-2016 (Tian et al., 2019; Xu et al., 2021), and our estimate of 0.44 Tg N year −1 fell well within this range.Additionally, our estimated CH 4 emissions (averaged to 0.23 Tg C year −1 over the last two decades) was also close to previous estimates centered on an annual emission rate of ~0.25

F
Spatial pattern of average annual net greenhouse gas balance of U.S. croplands from the 1960s to the 2010s.Note that error bars in insets represent ±1SD of the net greenhouse gas balance in each decade.NW, NP, MW, NE, SW, SP, and SE represent the Northwest, Northern Plains, Midwest, Northeast, Southwest, Southern Plains, and Southeast climate hubs, respectively.Regional Climate Hubs were formed and classified by the Agriculture Research Service, Forest Service, and Natural Resources Conservation Service of the USDA to help land managers adapt to a changing climate(Steele & Hatfield, 2018), with data sourced from https:// www.clima tehubs.usda.gov/ .SOC-sequestered CO2 and non-CO 2 GHG emissions in U.S. croplands simultaneously.In this way, we could gain a more complete picture of the magnitude and spatiotemporal variations of the net soil GHG balance in U.S. croplands and uncover new insights.For instance, by disentangling changes in the relative contribution of SOC sequestration and N 2 O and CH 4 emissions in the net GHG balance of U.S. croplands over the past few decades, we found that N 2 O emission and SOC sequestration jointly controlled the net GHG balance over most croplands in recent decades (Figure 7), highlighting the importance of considering both non-CO 2 GHG emissions and SOC sequestration when developing climate mitigation strategies.Meanwhile, our results revealed that the proportion of regions dominated by SOC sequestration increased over time.Possible explanations for this trend may include a gradual flattening of soil N 2 O emissions due to stabilized N fertilization amounts, as well as increased crop biomass due to improved management practices and breeding and genetic technologies that in turn enhanced residue inputs into soils and promoted SOC sequestration.However, given that meeting growing food demand will inevitably lead to increased non-CO 2 GHG emissions due to fertilization(Cavigelli et al., 2012;Molotoks et al., 2018), and significant increases in crop biomass may be limited by advances in crop breeding and genetic technologies, our findings therefore re-emphasize the necessity of concurrently reducing non-CO 2 GHG emissions and enhancing SOC sequestration to achieve the goal of carbon-neutral agriculture, which is particularly important since non-CO 2 GHG emissions have already surpassed the amount of SOC-sequestered CO 2 in U.S. croplands.Additionally, our results showed that the adoption of reduced tillage practices led to a decrease in net GHG balance, which underscores the potential of leveraging CSA management practices (e.g., notillage and cover cropping) to mitigate net GHG balance in U.S. croplands(Huggins & Reganold, 2008;Prokopy et al., 2019).Considering that some CSA practices could have quite different impacts (or even antagonistic) on SOC sequestration and N 2 O and CH 4 emissions F I G U R E 7 Spatial distributions of the relative contribution of SOC sequestration and emissions of N 2 O and CH 4 to the net greenhouse gas balance of U.S. croplands from the 1960s to the 2010s.| 13 of 20 YOU et al. (Guenet et 2021), our estimation of net GHG balance also allows for more comprehensive and effective mitigation efforts to combat climate change.

For
instance, Lu et al. (2021) reported that N fertilization was the dominant driver contributing to N 2 O emissions from U.S. agricultural soils, which increased N 2 O emissions by 0.33 Tg N year −1 since 1900.Synthetic N fertilizer use in U.S. croplands increased substantially (Figures S7a and S8a), from 2.48 Tg N year −1 in 1960 to 11.8 Tg N year −1 in 2015, which greatly promoted non-CO 2 GHG emissions (especially N 2 O) and exacerbated global climate warming.Although N addition could simultaneously stimulate SOC ac-cumulation in croplands, we found that SOC climate benefits were largely offset by non-CO 2 GHG emissions (FigureS6a,b).Therefore, optimizing N fertilizer use rates is an imminent need for achieving overall maximum benefits among enhancing SOC sequestration, improving crop yields, and curbing non-CO 2 GHG emissions(Gerber et al., 2016;Xia et al., 2017).According to CRM's survey data, the proportion of U.S. croplands adopting no-tillage practices increased significantly over the past three decades (FiguresS7c and S8c).Our factorial analysis suggested that the tillage intensity reduction suppressed net GHG emissions, which was consistent with other U.S.
corn-soybean cropping systems contributed to a net SOC accumulation of 1.0 Tg C year −1 during 1998-2008.Lu et al. (2022) reported a reduction rate of −5.5 Tg CO 2 -eq year −1 in GHG emissions from U.S. corn-soybean cropping systems during 1998-2008 as a result of tillage intensity reduction.Huang et al. (2022) suggested attenuated soil CO 2 and N 2 O emissions from Kentucky croplands under no-tillage compared to conventional tillage.Reduced net GHG emissions were also associated with an increase in irrigated U.S. cropland acreage over

F
Factorial contributions of multiple agricultural management practices and environmental forcings to changes in the net greenhouse gas balance of U.S. croplands from the 1960s to the 2010s, in comparison to the average state in the 1950s.Nfer represents nitrogen fertilizer use; Ndep represents atmospheric nitrogen deposition; LULC represents land use and land cover change (reflecting both cropland abandonment and expansion, as well as interannual crop rotation changes); and CO 2 represents atmospheric carbon dioxide concentration.Note that the sum of factorial contributions of individual drivers (i.e., stacked bars) does not equal to net changes in the net greenhouse gas balance (i.e., black line) due to interaction effects.emissions due to LULC relative to the average state of the 1950s.Overall, our factorial analysis of multiple agricultural management practices indicated large potential for CSA practices (e.g., optimized N fertilizer use and reduced tillage) to mitigate non-CO 2 GHG emissions and enhance SOC sequestration in agricultural soils, thereby leading to an overall reduction in net GHG emissions.
Our factorial analysis of each individual gas also indicated a positive response of climate warming on SOC decomposition, and N 2 O and CH 4 emissions (FigureS6).Additionally, the contrasting impacts of climate change on net GHG balance before and after the 1990s can be illustrated by the different response of N 2 O emissions to climate change (FigureS6b).Positive effects of increased N deposition on net GHG emissions were also reported in other studies sistent with observations and previous studies, some limitations remain in this study.These limitations include uncertainties arising from model forcing datasets, parameters, and structure, the lack of sufficient SOC-related measurements to constrain model simulations, and the neglect of indirect net GHG emissions due to N leaching in agroecosystems.First, model forcing datasets could introduce some uncertainties.For example, the crop-specific N fertilization data were reconstructed based on state-level surveys, which could not reflect actual spatial variations of fertilizer use in both magnitude and timing.Tillage intensity data were only available at the countylevel for recent decades, which lacked detailed spatial information and would inevitably introduce extrapolation errors in earlier years.Additionally, the assumption of crop residue removal percentage associated with tillage practices could also result in large uncertainties in the amount of crop residue inputs into soils, as it may diverge substantially from actual patterns.Therefore, joint community efforts to further improve model forcing datasets are needed.Second, the simplifications or omissions of real-world biophysical, biogeochemical, and hydrological processes in the DLEM, along with the under-representation of individual gas responses to various environmental factors, may also cause simulation biases.For example, the current DLEM representation of groundwater and irrigation practice (i.e., without considering irrigation amount and frequency) is relatively simple, which could lead to biased simulated soil moisture that, in turn, could affect GHG emission predictions.Furthermore, several studies have shown that soil freezing and thawing events induced nonnegligible amounts of N 2 O emissions(Del Grosso et al., 2022;Wagner-Riddle et al., 2017); nevertheless, the effects of soil freeze-thaw cycles on GHG emissions were not included in our current simulations and therefore constituted a possible source of deviation in our results.It is important to acknowledge that our model structures are inherently incomplete and uncertain, contributing to the overall uncertainty in model simulations.
et al.(2020), who applied a Bayesian model analysis framework incorporating the sampling importance resampling scheme to assess uncertainties in SOC estimates caused by both model parameters and structures.However, we concur withMarshall et al. (2007) that a robust estimation of model structural uncertainty requires the use of multiple models.Initiatives such as the Coupled Model Intercomparison Project(Eyring et al., 2016) and the Global N 2 O Model Intercomparison Project(Tian et al., 2018) provide valuable templates for this purpose.Therefore, we call for the initiation of multiple model inter-comparison projects, with a specific focus on net GHG balance, to comprehensively quantify uncertainty from model structures.Third, model parameters may also introduce biases into our results.Although we have performed a simplified estimation of parameter-induced uncertainty by assessing the sensitivity of the simulated net GHG balance to model parameter variations, the random parameter samples were generated based on their prior probability distribution functions, which may not accurately reflect their real distributions and could potentially cause biases(Wang & Chen, 2012).As such, a more thorough parameter uncertainty analysis is still needed, using random parameter samples generated based on their real or posterior probability distributions conditioned on measurements.Next, the lack of available spatialized and temporal datasets on SOC sequestration rate to constrain our model simulations over space and time could also result in significant uncertainty in our results, as SOC sequestration rates are typically much smaller than the magnitude of SOC pools.Moreover, most measurements of SOC sequestration rate are concentrated in shallow soils (e.g., 0-20 cm), while measurements in deeper soil layers remain sparse.Given that deep soil carbon represents over half of the global SOC stocks and is becoming more susceptible to climate change and human inter-

Table S1
and F CH 4 −C were annual fluxes of CO 2 , N 2 O, and CH 4 , respectively; SOCSR was SOC sequestration rate; molecular weight conversion fractions 44/12, 44/28, and 16/12 were used to convert the mass of CO 2 -C, N 2 O-N, and CH 4 -C into CO 2 , N 2 O, and CH 4 , respectively; GWP CO 2 , GWP N 2 O , and GWP CH 4 were GWP constants indicating radiative forcing of CO 2 , N 2 O, and CH 4 Comparisons of SOC sequestration rate and fluxes of N 2 O and CH 4 from other studies.
TA B L E 2