Knowledge-guided machine learning can improve carbon cycle quantification in agroecosystems

Domain: earth_science
Published: 2024-01-08
Editorial outcome: major_revision

Abstract
Accurate and cost-effective quantification of the carbon cycle for agroecosystems at decision-relevant scales is critical to mitigating climate change and ensuring sustainable food production. However, conventional process-based or data-driven modeling approaches alone have large prediction uncertainties due to the complex biogeochemical processes to model and the lack of observations to constrain many key state and flux variables. Here we propose a Knowledge-Guided Machine Learning (KGML) framework that addresses the above challenges by integrating knowledge embedded in a process-based model, high-resolution remote sensing observations, and machine learning (ML) techniques. Using the U.S. Corn Belt as a testbed, we demonstrate that KGML can outperform conventional process-based and black-box ML models in quantifying carbon cycle dynamics. Our high-resolution approach quantitatively reveals 86% more spatial detail of soil organic carbon changes than conventional coarse-resolution approaches. Moreover, we outline a protocol for improving KGML via various paths, which can be generalized to develop hybrid models to better predict complex earth system dynamics.

Introduction
Crop production systems and their interactions with the environment, known as agroecosystems, cover about one-third of the Earth’s land surface. As soil constitutes the largest single carbon reservoir on land, agroecosystems play a key role in the global terrestrial carbon cycle through crop interactions with soils and atmosphere1,2. Globally, agriculture is a significant source of greenhouse gasses (GHGs); yet, carbon uptake by crops also removes large amounts of carbon dioxide (CO2) from the atmosphere, some of which can be stabilized in soil3. Because most intensively cultivated soils are carbon-unsaturated, practices that increase soil organic carbon (SOC) represent a low-cost, large-scale strategy for reducing atmospheric GHG concentrations4–6. Thus, it is essential to accurately quantify carbon fluxes and changes in SOC in agroecosystems so that appropriate and effective conservation practices can be identified for any given location.

Increasing agricultural carbon sequestration is a key strategy for mitigating climate change. Significant efforts and investments have been made in the U.S. and across the globe to implement programs that incentivize SOC enrichment7,8. In light of these initiatives, it is important to develop robust and scalable methods for reliably quantifying field-level carbon sequestration, both to assess the climate mitigation effect and to ensure that mitigation actions by individual farmers are compensated fairly and accurately. Traditional carbon quantification methods that rely on soil sampling, emission factors, and process-based (PB) modeling entail inherent barriers to achieving the required levels of accuracy, scalability, and cost-effectiveness9–11. In particular, high spatial heterogeneity and seasonality due to variations in environmental conditions, crop types, and management practices present challenges for accurately quantifying carbon budgets12. While PB modeling approaches incorporate scientific knowledge, large uncertainties arise in PB models if local- and crop-specific parameters are not calibrated properly or if the underlying mechanisms are oversimplified or incompletely represented12,13. Additionally, PB models with detailed representations of existing scientific principles can be computationally prohibitive when applied to large regions at high spatial-temporal resolution (e.g., 250 m daily). On the other hand, data-driven machine learning (ML) approaches have the potential for high computational efficiency and accuracy14–17 but suffer from out-of-sample prediction failure in the absence of large training datasets, which are unavailable for most agricultural applications. Moreover, the results of ML models are often uninterpretable due to their black-box nature18. Therefore, new methods19 are needed to overcomethe limitations of PB and ML models, enabling cost-effective, accurate, and interpretable measurement and monetization of carbon outcomes at the individual field level. This will reduce errors in aggregated quantifications and promote more sustainable land management practices12,20.

The growing field of knowledge-guided machine learning (KGML)21,22 provides a promising methodology that combines the advantages of PB models, ML models, and multi-source datasets (e.g., in-situ and remote sensing data). Existing KGMLs can successfully model certain Earth systems for which dynamic processes are well-represented by established governing equations, such as hydrology and atmospheric sciences14,22–26. However, complex and crucial ecosystem processes such as biogeochemical cycling are mathematically non-linear and substantially more complicated. Furthermore, unlike surficial systems, soil interactions in agroecosystems largely cannot be directly observed by remote sensing, whereas in-situ direct measurements are often expensive and limited. Therefore new KGML approaches must be developed to incorporate sufficient biogeochemical knowledge and effectively assimilate indirect measurements (e.g., remote sensing and survey data) to capture terrestrial processes that are less directly observable27,28.

To address the existing gaps in carbon budget modeling capabilities, we developed a novel KGML framework that combines prior biogeochemical knowledge of carbon dynamics with a deep learning model to generate reliable predictions of agricultural carbon fluxes, crop yields, and changes in soil carbon stocks (KGML-ag-Carbon, Fig. 1). In-situ eddy covariance (EC) flux tower data, regional survey yield data, remotely-sensed gross primary production (GPP) data, and synthetic data generated by a PB model were assimilated into KGML-ag-Carbon. The model effectiveness is demonstrated here for corn and soybean production in the U.S. Midwest (Fig. S1), with highly accurate outputs for carbon fluxes, crop yields, and changes in soil carbon at high spatial (250 m) and temporal (daily) resolution, providing usable data for land managers. We also analyzed the improvement resulting from each KGML component and certain biogeochemical responses. The model design presented here exemplifies a solution to challenges in simulating dynamic heterogeneous systems, which will help advance broader applications of KGML for understanding Earth processes.

Results
Overview of the KGML-ag-Carbon framework
KGML-ag-Carbon is a novel framework combining process-based understanding and advanced AI approaches for simulating complex biogeochemical cycles under intensive management practices for agroecosystems. KGML-ag-Carbon distinguishes itself from previous KGML applications in other disciplines by its use of a well-validated PB model, ecosys29, and its ability to directly assimilate remote sensing data. The deep learning model based on the gated recurrent unit (GRU) mechanism30 was used to develop KGML-ag-Carbon’s architecture (Fig. 1). The ecosystem theory of carbon allocation is the basis for the ecosys model, which was used to design the hierarchical structure of KGML-ag-Carbon, including submodules related to plants, soil, atmospheric carbon exchange, and feedback of carbon from plants to soil (Fig. S2). Importantly, outputs from the ecosys model provided synthetic data on ecosystem carbon allocation, associated fluxes, and environmental responses that were used to pre-train the KGML-ag-Carbon model. This pre-training step confers a major advantage to the model by improving the predictive ability with a minimal amount of labeled samples and accelerating the convergence in model tuning using labeled samples. GPP data, which represents the dominant carbon input to agroecosystems, was retrieved from remote sensing observations and assimilated into KGML-ag-Carbon as spatial constraints.

KGML-ag-Carbon resolves the major carbon budget components, including autotrophic respiration (Ra), heterotrophic respiration (Rh), total ecosystem respiration (Reco, Ra+Rh), and net ecosystem carbon exchange (NEE) on a daily scale, and yield on an annual scale. As in natural ecosystems, changes in agroecosystem soil carbon storage are determined by the mass balance of input and output carbon fluxes31,32. Ecosystem carbon inputs originate from plant photosynthesis, i.e., gross primary production (GPP), while soil carbon inputs include both aboveground and belowground litter and root exudates. Carbon outputs occur through respiration, including Ra from plant shoots and roots and Rh from SOC decomposition by microbes and fungi. Disturbances such as harvesting also remove carbon from the ecosystem periodically. Based on the carbon fluxes and yield estimated from KGML-ag-Carbon, annual changes in SOC can be determined using the mass balance equation ΔSOC = -NEE - Yield12,13,33.

We systematically explored multiple paths for improving the prediction performance of KGML-ag-Carbon, including pre-training the model with synthetic data and incorporating knowledge-guided (KG) loss functions, which addresses broader issues about reducing uncertainty for hybrid modeling. Over 14 million synthetic data and various KG loss functions were used to pre-train KGML-ag-Carbon to learn the prior knowledge from the PB model (Fig. 1). Using synthetic data generated by a PB model is several orders of magnitude cheaper than the cost of collecting real-world observations. The KG loss functions include biogeochemical/physical constraints such as mass balance (Ra+Rh-GPP=NEE), prediction thresholds (e.g., 0<yield<GPP), and responses of outputs to inputs (e.g., Rh should monotonically increase with SOC content under other fixed conditions).

A subset of the observed yield data from USDA (320 out of 630 counties) along with Reco (Ra+Rh) and NEE data from 11 cropland EC flux tower sites were then used to fine-tune the pre-trained KGML-ag-Carbon model to improve its prediction ability for real-world carbon budgets (Fig. 1). KG loss functions with biogeochemical/physical constraints similar to those used in pre-training were included in the fine-tuning. However, to preserve sufficient pre-training knowledge, the synthetic data were merged with observational data during fine-tuning, while extra constraints were added to the KG loss functions to maintain the responses of outputs to inputs (e.g., changes in Ra responses to the environment remain within 10%). Details on the structural development, datasets utilized, and training strategies for the KGML-ag-Carbon model are provided in the Methods section.

Model performance in crop yield and carbon flux predictions
We evaluated the performance of KGML-ag-Carbon both before and after fine-tuning along with the sensitivity of model performance to the real-world training sample size (Figs. 2, S4, S6, S8). As an initial check on pre-training effectiveness, the pre-trained KGML-ag-Carbon model results for the test set of synthetic data (two years out of the 18-year period) were compared with ecosys simulations and found to be highly consistent, with R2 values of 0.99, 0.99, 0.97, and 0.97 for yield, Ra, Rh, and NEE, respectively (Fig. S4). R2 values of the ecosys model for corn and soybean yield predictions were 0.49 and 0.42, respectively, as benchmarked with observed county-scale crop yields, while values for daily Reco and NEE predictions were 0.67–0.89 and 0.59–0.88, respectively, compared with measurements from EC sites (green stars/boxes in Fig. 2, derived from Zhou et al.13). After fine-tuning, the R2 of KGML-ag-Carbon for corn and soybean yield predictions on a test set of 210 counties were 0.91 and 0.88, while values for daily Reco and NEE predictions tested on 2 years of out-of-sample data from 11 EC flux tower sites were 0.94 and 0.96, respectively (Fig. S8).

The robustness test (Fig. 2) reveals that compared with a GRU-based pure ML model using the same inputs, KGML-ag-Carbon both consistently outperforms the pure ML model and has much lower sensitivity to the number of real-world training samples. In summary, the major differences between KGML-ag-Carbon and pure ML lie in the additional pre-training process and the customized model structure and loss functions guided by known scientific knowledge. These advancements enhance the optimization process of the ML model and allow reliable predictions to be made with fewer labeled samples. The pure ML model performance approached that of KGML-ag-Carbon at large sample sizes but performed poorly with small training sets, particularly for crop yield (Fig. 2a, b). The reduced need for training samples is a central advantage of KGML-ag-Carbon for real-world crop yield and carbon flux estimates because available training data are usually limited and collecting data from physical sampling is costly. The improvements in carbon flux predictions with increasing sample size mostly arise from capturing the interannual carbon dynamics, which can be more easily learned from seasonal patterns of GPP and climate (Fig. 2c, d). Even without fine-tuning (a training sample size of 0), the KGML-ag-Carbon by assimilating the GPP data as input, can outperform both ecosys model and pure ML model trained with small training samples (Fig. 2).

We conducted several additional experiments to evaluate the performance of KGML-ag-Carbon under different training, validation, and testing dataset splits, including out-of-sample performance of yield predictions in the spatial and temporal domains (Fig. S9a-f). For example, we used data from Illinois for testing and data from other states for training and validation and used several continuous years of data for training and validation with other years for testing. We also examined the effect of using extreme years with exceptionally high or low yields for testing and other years for training and validation. We note that the KGML-ag-Carbon model outperforms pure ML and process-based models in predicting yield in extreme years (Fig. S9e, f) primarily because it is constrained by both observations and synthetic data generated from the PB model. Similarly, out-of-sample performance for Reco and NEE predictions was investigated on both annual and daily temporal scales and specifically examined for sites (EC flux towers US-NE1–3) with longer-term observations covering 2001–2019 (Fig. S9g–l). The results demonstrate consistently better performance of KGML-ag-Carbon compared to pure ML in all tested situations. The details of robustness tests on the KGML-ag-Carbon performance can be found in the Methods section.

Pathways to reduce KGML-ag-Carbon uncertainty
To understand the contribution of different strategies to improvements in the performance of KGML-ag-Carbon, we conducted full-factorial tests to include or exclude different model components and selected five representative models to use in interpreting the results (Fig. 3). The results reveal that using GPP data as an input and pre-training KGML-ag-Carbon with synthetic data contribute most to improving the performance of KGML-ag-Carbon relative to other strategies. When using larger real-world observations for model fine-tuning, the GPP data has the largest contribution to improving KGML-ag-Carbon performance; while pre-training with synthetic data is more important when using smaller real-world observation sets for model fine-tuning (Fig. 3a, b). This indicates that under data-limited situations, pre-training based on datasets generated by process models with sufficiently well-represented mechanisms can provide prior knowledge to significantly help improve the performance of ML. In contrast, when good-quality observational datasets are available, the ML model can learn complex relationships directly from the data so pre-training is less important.

The improvements to KGML-ag-Carbon provided by the hierarchical structure and KG loss functions are relatively small compared to those from GPP inputs and pre-training processes. One potential reason is that the model performance metrics (R2 here) were already very high after adding GPP inputs and pre-training (Fig. 3a, b). However, including the hierarchical structure and KG loss functions significantly increases the ability of KGML-ag-Carbon to capture complex carbon flux dynamics and the interpretability of the predictions.

The pre-training process significantly reduces the residual mass balance of carbon fluxes (i.e., GPP-Reco-NEE), while the hierarchical structure and KG loss functions further reduce the mass balance residual to near-zero (Fig. 3c), indicating that the inclusion of hierarchical structure and KG loss functions constrain the model to follow physical rules. Although KGML-ag-Carbon has an overall performance similar to the ML+GPP and ML+GPP+pretrain models, significant improvements in NEE predictions are achieved in winter and summer, especially over periods with complex dynamics (Fig. 3d). The advantage of KGML-ag-Carbon over other models for various time periods is mainly attributed to the incorporated knowledge of distinct temporal patterns in Ra and Rh, such as Rh equals Reco in winter when plant growth is absent, which can be utilized to separately improve the performance of the Ra and Rh submodels in KGML-ag-Carbon. Other results from full-factorial tests and mass balance tests can be found in Fig. S11 and Fig. S12, respectively. The details outlining how the contributions of KGML-ag-Carbon components were identified can be found in the Methods section.

High-resolution carbon flux estimates across the U.S. Midwest
Using the fine-tuned KGML-ag-Carbon, we predicted regional daily carbon fluxes and annual crop yields across the U.S. Midwest at a 250-m spatial resolution (smaller than a typical U.S. Midwest field, Fig. 4a–c). Inputs that drive KGML-ag-Carbon include daily climate data from NLDAS-2, topsoil properties from gSSURGO, remotely sensed daily GPP from the SLOPE product34, and crop rotation information. The high spatial resolution is facilitated by high-resolution soil information, crop rotation maps, and the GPP product. The high temporal resolution comes from the climate and GPP product data, which provide daily information on environmental conditions and ecosystem carbon inputs. The procedures for generating high-resolution predictions across the U.S. Midwest are given in the Methods section.

The multi-year-averaged SLOPE GPP data and the carbon fluxes generated by KGML-ag-Carbon are closer to the EC flux tower observations (same dataset used in Fig. 2) than estimates from Trendy35, a widely used carbon flux ensemble product generated by a suite of dynamic global vegetation models (Fig. 4). Although methods for evaluating KGML-ag-Carbon are somewhat limited at the regional scale, EC flux tower data and Trendy are suitable datasets for comparison in the absence of ideal benchmarks for this large region. The distributions of GPP values are similar across SLOPE, EC flux tower observations, and the Trendy ensemble (Fig. 4d). The distributions of Reco and NEE estimated by KGML-ag-Carbon are similar to those of EC flux tower observations, but the ensemble of Trendy models overestimated both Reco and NEE in the U.S. Midwest compared with flux tower observations (Fig. 4e, f). In addition, the estimated distributions of GPP, Reco, and NEE vary widely among individual Trendy models, which may arise from differences in structure and parameters among models. This reflects the large uncertainties remaining in current state-of-the-art PB models for carbon budgets, especially for regional-scale estimates. To summarize these comparisons, the carbon fluxes estimated by KGML-ag-Carbon demonstrate high spatial-temporal resolution and accuracy, providing a novel product for precise regional-scale carbon budget quantification down to a single field.

Methods
Synthetic pre-training data for the KGML model
We used the agroecosystem model ecosys to generate synthetic data for crop yield, ecosystem autotrophic respiration (Ra), ecosystem heterotrophic respiration (Rh), net ecosystem exchange (NEE), and gross primary production (GPP). This synthetic data was used to pre-train the Knowledge-Guided Machine Learning for the Agricultural Carbon budget model (KGML-ag-Carbon). Ecosys simultaneously simulates carbon, water, and nutrient cycles within the soil and plant system based on biophysical and biochemical principles29. Its ability to simulate carbon fluxes and crop yields has been extensively validated across midwestern U.S. cropping systems13. We conducted county-level simulations using ecosys for 293 counties in the states of Illinois, Iowa, and Indiana using climate data from the North American Land Data Assimilation System (NLDAS-2 and soil data from the Gridded Soil Survey Geographic Database (gSSURGO). The synthetic database contains 10,335 simulations whose inputs include soil information, planting and harvest dates, crop parameters, and crop rotation information randomly selected from among predefined ranges to ensure a representative synthetic database. Within each county, the soil information was randomly selected from among the top 10 dominant cropland soil types in each county. The predefined range of planting dates is from April 15 to June 10, and the harvest date is from October 31 to November 20, which represents the general crop calendar in this region. In the database, one-third of the total simulations have corn-soybean rotations, one-third are soybean-corn rotations, and the remaining had corn and soybean planted randomly from 2001 to 2018 to represent common rotation strategies within this region.

Datasets for fine-tuning, validation, and extrapolations
We fine-tuned and validated KGML-ag-Carbon for crop yield estimation over 637 counties and carbon fluxes estimation (i.e., Ra, Rh, NEE) at 11 cropland EC flux tower sites located within major U.S. corn and soybean production regions (Fig. S1). For fine-tuning and validation of the regional crop yield submodule, the 250-m daily GPP product derived from machine learning models based on Soil-Adjusted Near-Infrared Reflectance of vegetation (SANIRv)34, 0-30 cm gSSURGO soil properties, NLDAS-2 climate data, and crop type information (CDL and CSDL63 were used after and prior to 2008, respectively) were used as KGML-ag-Carbon inputs. County-scale corn and soybean yields from NASS, and USDA (https://quickstats.nass.usda.gov/) were used as benchmarks for fine-tuning and validation. For carbon flux estimation, we used the daily Ra, Rh, and NEE decomposed from observed NEE at EC flux tower sites using the ONEFlux tool61, observed climate data from EC flux towers, gSSURGO soil information, and CDL crop type, while observed NEE and Reco from EC flux tower data were used as benchmarks. Because the daily GPP data we used is an average of GPP decomposed from NEE using different daytime and nighttime partition methods, it may not preserve the mass balance among NEE, Reco, and GPP; thus, we corrected the EC flux tower daily GPP by replacing it with observed Reco plus NEE in the following analysis. For estimating carbon fluxes at the regional scale, we used SANIRv-based GPP, NLDAS-2 climate data, and gSSURGO soil information as the model inputs.

The structure of KGML-ag-Carbon
KGML-ag-Carbon uses a hierarchical structure64 to incorporate the causal relations between different variables and processes with ecosystem knowledge for guidance, as presented in Fig. S2. It contains five submodules, including (1) a GRU_Ra submodule for daily Ra estimation, (2) a GRU_Rh submodule for daily Rh estimation, (3) a GRU_NEE submodule for daily NEE estimation, (4) an attention module for crop yield estimation, and (5) a GRU_Basis submodule to connect and support the other four modules. We used a type of recurrent neural network called a Gated Recurrent Unit (GRU) as the basic machine learning module to develop our model structure. GRU has been proven to perform similarly to Long short-term memory (LSTM65) in using cell states as internal memories to preserve historical information; however, GRU uses a simpler structure with fewer hidden states compared to LSTM and thus often remains more stable with a small number of training samples28,30,66.

The recursive representations of GRU can be represented as:
ht = zt ⊙ nt + (1-zt) ⊙ ht-1 (1)
zt = g(Wxz xt + bxz + Whz ht-1 + bhz) (2)
nt = tanh(Wxn xt + bxn + rt ⊙ (Whn ht-1 + bhn)) (3)
rt = g(Wxr xt + bxr + Whr ht-1 + bhr) (4)
where ht is the hidden state at time t, xt is the input at time t, and ht-1 is the hidden state at time t-1 or the initial hidden state at time 0. z, n, and r are the update, reset, and new gates, respectively. g is the sigmoid function and ⊙ is the Hadamard product. Wxz, Whz, Wxn, Whn, Wxr, and Whr are learnable linear transformation matrices. bxz, bhz, bxn, bhn, bxr, and bhr are corresponding learnable bias vectors.

Each GRU cell in KGML-ag-Carbon represents a GRU with 64 hidden units (h vector dimension = 64), and each dense cell is a linear transformation layer, which can be represented as:
yflux,t = Why ht + bhy (5)
where yflux,t is the predicted flux target variables at time t, including Ra, Rh, and NEE. Why and bhy are the learnable weight and bias, respectively. The GRU_basic, GRU_Ra, and GRU_NEE submodules have one layer of GRU cells while GRU_Rh has two layers of GRU cells. 20% of the output hidden states from GRU cells are randomly dropped by replacing them with zero values (the so-called 20% dropout) to avoid overfitting.

The attention module in KGML-ag-Carbon is a modified version of the traditional LSTM attention model67, containing two layers: ATTN_Weight and ATTN_Densor. ATTN_weight can be represented as:
αt = exp(e(t)) / Σt=1 365 exp(e(t)) (6)
et = tanh(Whe,4 ReLU(Whe,3 ReLU(Whe,2 ReLU(Whe,1 ht + bhe,1) + bhe,2) + bhe,3) + bhe,4) (7)
where αt is the probability attention score calculated from a softmax function, representing the importance of time over the whole year. et is the weight score of ht at time t calculated from a 4-layer feedforward neural network (FNN) with a Rectified Linear Unit (ReLU) as the activation function for the first three layers and a hyperbolic tangent (tanh) function for the last layer. Whe,i and bhe,i are the learnable weight and bias for the ith layer in the FNN, respectively (i=1, 2, 3, and 4).

αt and ht are then multiplied in the ATTN_Densor layer to calculate the annual yield:
yyield = Wcy,4 ReLU(Wcy,3 ReLU(Wcy,2 ReLU(Wcy,1 c + bcy,1) + bcy,2) + bcy,3) + bcy,4 (8)
c = Σt=1 365 αt ht (9)
where yyield is the predicted yield for the input year, calculated from a 4-layer FNN with ReLU as the activation function for the first three layers. c is the self-weighted context vector, which has the same dimensions as the hidden state. Wcy,i and bcy,i are the learnable weight and bias for the ith layer in the FNN, respectively (i=1, 2, 3, and 4). The attention module for yield collects simulated information of each day from the GRU_basis submodule as input and weighs the contribution of each day’s information to the final yield prediction.

Crop annual residue yres can be expressed as:
yres = ReLU(Σt=1 Tx (GPP) - Σt=1 Tx (Ra,t) - yyield) (10)
where GPP and Ra,t are the GPP input and predicted Ra at time step t; yyield is the annual predicted yield; and Tx is the number of days in the input time series (in this study Tx = 365). The ReLU function is used to prevent a situation in which the sum of predicted annual yield and Ra is bigger than the annual GPP.

The KGML-ag-Carbon inputs (X) include seven daily climate variables: surface downward shortwave radiation (RADN, MJ m-2 day-1), maximum air temperature (TMAX_AIR, °C), the difference between the maximum and minimum air temperature (TDIF_AIR, °C), maximum humidity (HMAX_AIR, kPa), the difference between the maximum and minimum humidity (HDIF_AIR, kPa), wind speed (WIND, km day-1), and precipitation (PRECN, mm day-1). Additional inputs are daily GPP (gC m-2 day-1), year, crop type (corn/soybean), and nine soil properties averaged from 0 to 30 cm soil depth: bulk density (TBKDS, Mg m-3), sand content (TCSAND, g kg-1), silt content (TCSILT, g kg-1), water content at field capacity (TFC, m3 m-3), water content at wilting point (TWP, m3 m-3), saturated hydraulic conductivity (TKSat, mm h-1), soil organic carbon (TSOC, gC kg-1), pH (TPH), and cation exchange capacity (TCEC, cmol+ kg-1). To increase the efficiency of the training process, we used the Z-normalization method to normalize each variable separately on synthetic data. The Z-normalization method can be expressed as:
Z = (x - μ) / σ (11)
where Z is the normalized variable; x is the vector of a particular variable over all the data samples in the dataset; μ is the mean value of x; and σ is the standard deviation (STD) of x. The scaling factors derived from the ecosys synthetic data for each variable were used to normalize observed data into the same ranges as synthetic data. TDIF_AIR and HDIF_AIR were used instead of absolute minimums of temperature (TMIN_AIR) and humidity (HMIN_AIR) because TMIN_AIR and HMIN_AIR follow similar trends as TMAX_AIR and HMAX_AIR, causing Z-normalization to be poorly defined numerically. Using the difference between maximum and minimum values provides clearer information about daily air temperature and humidity variations.

Xt are the inputs to the submodules of GRU_Basis, GRU_Ra, GRU_Rh, and GRU_NEE. Additionally, the output hidden states from GRU_Basis are inputs to GRU_Ra, GRU_Rh, and the attention module. The predicted annual yield, daily GPP, and daily Ra are then used to calculate the carbon in crop annual residue after harvest in Residue_layer. The annual residue is fed back to the soil for Rh calculation by inputting it on the 300th day of the year to GRU_Rh to assess the relationship between soil and plant carbon pools. Finally, the GRU_NEE takes predicted daily Ra and Rh together with Xt as input to predict daily NEE to assess the contribution of different carbon fluxes to NEE.

Training strategies for KGML-ag-carbon
We used a five-step training method to train KGML-ag-Carbon with ecosys-generated synthetic data and observed data, including (1) pre-training yield and Ra submodules using synthetic data, (2) pre-training Ra, Rh, and NEE submodules using synthetic data, (3) fine-tuning the yield submodule using observed data, (4) retraining Ra, Rh, and NEE submodules using synthetic data, and (5) fine-tuning Ra, Rh, and NEE submodules using observed data (Table 1). We utilized an enhanced mini-batch learning strategy68 to effectively capture and maintain long-term temporal dependencies in the model. The best-performing submodules in the validation set at each step are saved for training in the next step.

Specifically, the KGML-ag-Carbon model was pre-trained using synthetic data to gain prior knowledge in steps 1-2, with two years of data randomly selected from the 18-year period of synthetic data for model validation, while the remaining 16 years of data were used for model training. In step 1, we trained the yield and Ra submodules together since they are crop-related variables and are used together for crop residue calculations (Eq. (10)), with the GRU_Rh and GRU_NEE submodules “frozen” by setting the learning gradient to zero. We used a mean-square-error (MSE)-based self-paced learning (SPL) method69,70 to build our training losses to train the model from “easier” samples to “harder” samples (Note S1).

In step 2, we further pre-trained the submodules for Ra, Rh, and NEE prediction together with the knowledge-based losses and responses by freezing the attention module and GRU_Basis module, considering the relationship of carbon fluxes. Besides the MSE loss, the loss function for step 2 also involves (1) the knowledge of mass balance (GPP-Ra-Rh = -NEE, considering the positive NEE direction to be from soil to atmosphere) to control the relationship between the input GPP and predicted Ra, Rh, and NEE, (2) the partial dependence plot (PDP, Fig. S3) to control the response of Rh to TSOC (Note S2). After two steps of pre-training, the KGML-ag-Carbon can successfully imitate ecosys for simulating yield, Ra, Rh, and NEE (Fig. S4).

In step 3, we fine-tuned the yield submodule with county-level crop yield data. The GRU_Ra, GRU_Rh, and GRU_NEE submodules were fully frozen and the GRU_Basis submodule was partially frozen by setting the learning rate to 20% of the original one. We included the knowledge-guided constraints in the loss function to control the range of yield (bigger than 0 and less than 0.5 times annual GPP) and maintain three key responses (i.e., yield responses to TSOC, GPP, and year) learned from the PB model (Note S3, Fig. S5). Detailed information on using coarse resolution (county-level) yield data to fine-tune our high-resolution model (250 m) is described in supplementary Note S3.

Step 4 is similar to step 2 in terms of using synthetic data to train the Ra, Rh, and NEE submodules to avoid too much prior knowledge loss after fine-tuning the yield submodule. An experiment comparing Ra, Rh, and NEE prediction performance after step 4 and models with and without step 2 demonstrated the effectiveness of step 2 (Fig. S6). We attempted to remove step 2 and trained the model only in step 4 for carbon fluxes with 80 maximum epochs. The results showed a performance drop for Ra, Rh, and NEE pre-training, especially at the annual scale (Fig. S6; with step 2: annual RMSE=13.9, 24.4, and 28.9 gC m-2 day-1 for Ra, Rh, and NEE, respectively; without step 2: annual RMSE=17.0, 29.3, and 34.4 gC m-2 day-1).

Finally, we fine-tuned KGML-ag-Carbon using the daily observed Reco (Ra+Rh), NEE, and GPP data from 11 EC flux tower sites throughout the U.S. Midwest, with the GRU_Basis and Attention modules frozen (Fig. S2). The learning rates of the GRU_Ra, GRU_Rh, and GRU_NEE submodules were set to 50%, 20%, and 50% of the original one at the fine-tuning stage, respectively, to avoid overfitting and losing too much prior knowledge. The loss function for step 5 involves a similar mass balance constraint as step 2 but contains a different response constraint to preserve Ra and Rh responses to environmental variables learned from the processes-based model in data-sparse regions (Note S4, Fig. S7). Additionally, we introduced a method to separate the Ra and Rh during winter by assuming that most Reco during winter is from Rh since the selected EC flux tower sites were fallow during winter (Note S4). At each site, two years of data were randomly selected from the whole observed period as validation data, and the remaining data were used as training data. The final fine-tuned out-of-sample testing results are presented in Fig. S8.

Robustness test for the performance of KGML-ag-Carbon
To investigate the robustness of KGML-ag-Carbon for yield, Reco, and NEE predictions, we conducted several experiments with different training sample sizes to compare the performance of KGML-ag-Carbon with a pure ML model under different conditions (Fig. 2, Fig. S9).

For yield predictions, we first randomly sampled 210 counties out of all 637 counties in the U.S. Midwest from NASS data for testing and kept 100 counties from the remainder for validation (Fig. 2a, b). Specifically, to conduct the yield robustness test, different sample sizes of 2, 5, 10, 20, 40, 80, 160, 240, and 320 counties were randomly selected from the remaining data, with ensemble times of 50, 20, 10, 10, 5, 4, 2, 1, and 1, respectively. The KGML-ag-Carbon model was trained following the 5-step training strategy described above, with a varying training sample size for fine-tuning in step 3. A sample size of zero for fine-tuning was also considered by skipping step 3 for the KGML-ag-Carbon training. The pure ML model for yield prediction is a 2-layer GRU model with attention, which is similar to GRU_Basis combined with the attention module presented in Fig. S2, with the same input features as KGML-ag-Carbon (including GPP). The pure ML models were trained with a similar method as KGML-ag-Carbon in step 3 with doubled maximum training epochs but without a knowledge-guided loss (only MSE loss) and pre-training (all other steps). In addition, we conducted further experiments with different training/testing split methods such as (1) using counties except for Illinois for training and Illinois for testing to detect spatial transferability, (2) training on the prior few years of data and testing on latter years to detect temporal transferability, and (3) training on the normal years and testing on extreme years (Fig. S9a–f). We used counties from Illinois (100) for testing and randomly sampled 100 counties from the remaining states for validation to test the KGML-ag-Carbon model in an independent out-of-sample testing data set (Fig. S9a, b). The training sample selection method was the same as the random sampling method. For detecting temporal transferability, we trained the model with all counties but split the 21-year data into training/validation/testing periods (Fig. S9c, d). We selected the front 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20 years for training and validation, with the last 1, 1, 1, 2, 2, 2, 3, 3, and 3 years of the selected periods as the validation sets, respectively. The remaining years of the 21-year period were used for testing the model. The prediction performance of each testing year was calculated separately and presented in Fig. S9c, d. To test the performance in extreme years, we trained the KGML-ag-Carbon model and pure ML model with data from all counties but excluded the selected extreme years of 2002, 2003, and 2012 (Fig. S9e, f). The extreme years were selected by detecting the outliers (outside the range of mean ± two times the STD) for each year based on a yield distribution calculated from the detrended yield for all counties and all years (Fig. S10). 2002, 2003, and 2012 have the top three numbers of outliers, with 98, 89, and 349 counties.

For Reco and NEE predictions, we divided the 11 EC flux tower sites into 6 testing groups based on the spatial distribution to detect the spatial transferability of the KGML-ag-Carbon at different temporal scales (Fig. 2c, d, Fig. S9g–l). We conducted the ensemble experiments, and each time, we selected one group on which to test KGML-ag-Carbon and the pure ML model, which were trained and validated by randomly selected sites from the remaining groups. Specifically, we first divided the 11 eddy-covariance sites into 6 testing groups based on the spatial distribution, with US-Bo1 and 2 as group 1, US-Br1 and 3 as group 2, US-IB1 as group 3, US-KL1 as group 4, US-NE1, 2 and 3 as group 5, and US-Ro1 and 5 as group 6 (Fig. S1). We selected one site as the validation data for each group and selected different sample sizes of 1, 2, 3, 4, 5, 6, and 7 sites as the training data from the remaining sites. The validation data traversed each of the remaining sites and training data of the same size would be forced to be different from each other. For example, group 1 has 2 sites for testing. If we would like to choose a training/validation sample for sample size 5, we would first select one site from the remaining 9 sites (excluding 2 test sites) and randomly sample 5 sites from the remaining 8 sites (excluding 2 test sites and 1 validation site). The 5 sampled sites would be compared with the previously selected 5-site training data and if they are the same, the sampling would be applied again. This process was conducted 9 times to cover all of the remaining sites so that the ensemble count for each sample size in group 1 was 9. Similarly, the ensemble times for each sample size in groups 2, 3, 4, 5, and 6 are 9, 10, 10, 8, and 9, respectively. The KGML-ag-Carbon model was trained following the 5-step training strategy described above with the training sample size varying for fine-tuning in step 5. A sample size of zero for fine-tuning was also examined by skipping step 5 in the KGML-ag-Carbon training. The pure ML model is a multitask 2-layer 64-unit GRU for Ra, Rh, and NEE simulation with the same input as the KGML-ag-Carbon model. The pure ML models were trained by a similar method as KGML-ag-Carbon in step 5 with doubled maximum training epochs but without a knowledge-guided loss (only MSE loss) and pre-training (all other steps). Finally, we investigated the overall performance of Reco and NEE prediction by combining results from all ensemble experiments at daily and annual scales (Fig. 2c, d, Fig. S9g, h), and investigated the performance at one representative location (the area containing US-NE1, 2, and 3 with 19-year data at each site) at daily and annual scales (Fig. S9i–l).

Detecting the contributions of KGML-ag-Carbon components
To investigate the contributions of different KGML-ag-Carbon components to the final ready-to-go KGML-ag-Carbon performance, we conducted full-factorial tests for each component in the model and tested the model performance on an out-of-sample dataset (Fig. S11). Specifically, we included or excluded four components: (1) using GPP data as an input (GPP for short), (2) pre-training the model with synthetic data, (3) incorporating the KGML-ag-Carbon structure, and (4) implementing KG loss functions and the 5-step training strategy (if structure is applicable). In total, 16 individual models were trained. The training and testing data are similar to the robustness experiment described above. Specifically, to determine the contributions to the yield (flux) predictions, we used training sets of 5 and 40 counties (1 and 7 sites) to train the models, referred to as small and large training sample sets, respectively. The optimized models were tested on out-of-sample datasets, which are NASS yields from 210 randomly selected counties and Reco and NEE from 6 groups of EC flux tower sites (the models tested on one group were trained and validated with data from sites chosen from other groups). We calculated the mean and STD of the prediction accuracy for all the models from ensemble experiments and detected the performance changes by comparing the models with and without each KGML-ag-Carbon component (Fig. S11). To illustrate the factors that contribute to the KGML-ag-Carbon model performance, we selected five representative models from the 16 trained models to showcase the direction of performance improvement. These models include (1) ML, (2) ML+GPP, (3) ML+GPP+pre-training, (4) ML+GPP+pre-training+KG structure, and (5) ML+GPP+pre-training+KG structure+KG loss (Fig. 3a, b). To further detect the influences of knowledge-guided components (i.e., pre-training, hierarchical structure, and KG loss functions) on improving the prediction performances, we compared three kinds of models, including an ML model with GPP data, an ML model with GPP and pre-training, and KGML-ag-Carbon, regarding the mass balance residues of predictions and the performance in capturing complex daily fluxes for a representative site-year (US-NE1-year 2016; Fig. 3c, d; Fig. S12).

High-resolution predictions across the U.S. Midwest
After fine-tuning KGML-ag-Carbon with county-scale corn and soybean yield as well as EC flux tower observations from agroecosystem sites (Table 1, Step 5), the model was used to simulate regional annual crop yields and daily carbon fluxes (i.e., NEE, Ra, Rh, Reco) with a spatial resolution of 250 m over the main corn- and soybean-producing region of the U.S. Midwest (Fig. S1) from 2000 to 2020. To evaluate the performance of regional-scale carbon flux estimates, we compared the model results with Trendy35, which was generated by a suite of dynamic global vegetation models at a monthly scale with spatial resolutions of 0.5°×0.5° or coarser. The carbon flux values from this study were regridded to 0.5° by averaging the value of pixels within a 0.5° grid for comparison. The distribution of annual accumulated GPP, NEE, and Reco from these two datasets and the observations from the selected EC flux tower sites were used for the comparison (Fig. 4). The wide range of variation observed in the Trendy models ensemble can be attributed to the inclusion of diverse processes and alternative parameterizations adopted by models from different research communities, as described by Sitch et al.35.

Investigating the benefits of high-resolution quantification
To generate 0.0025-degree-resolution ΔSOC estimates for the U.S. Midwest (Fig. 5a–c), we employed the mass balance equation ΔSOC = -NEE - crop yield12,13,33 over the period 2000–2020. Specifically, we regridded the 250-m-resolution NEE and crop yield estimations from KGML-ag-Carbon into 0.0025° estimations for use in the mass balance equation. To minimize the influence of undecomposed surface crop residues, which do not contribute to ΔSOC but are counted as part of our ΔSOC estimations through the mass balance approach, we selected the 21-year averaged value of ΔSOC. We then focused on regions where more than 50% of the area was planted with corn or soybean crops (Fig. S13a). The ΔSOC values were converted to percentage fractions (Fig. 5b) using ML-based SOC stocks derived from SoilGrids71 (Fig. S13b). Specifically, we used corn and soybean fractions from CDL and CSDL data (Fig. S13a) to exclusively identify corn and soybean agroecosystems (total fraction > 0.5). This alignment with our model’s current training scope helped reduce the mixed pixel effect resulting from inputting remotely sensed GPP data from other ecosystems. The SoilGrids SOC stock (Fig. S13b) was derived from organic carbon density (OCD) in each layer of the 200-cm soil depth at a 250-m resolution71.

To attribute the spatial patterns of estimated ΔSOC, we conducted Pearson correlation analyses between the input variables (including seven climate variables and nine soil variables) and the target variables (including GPP, NEE, Ra, Rh, Yield, Reco, Residue, and ΔSOC) (Fig. S14). In our approach, each variable was temporally aggregated to a 21-year scale and Z-normalized using Eq. (11). The Residue variable was computed as GPP-Ra-Yield, representing the net carbon return from plants to the soil. While GPP served as an input to the KGML-ag-Carbon model, we included it as a target variable in the correlation assessment due to its pivotal role in the carbon cycle. In addition, we conducted a multiple linear regression to assess the total influence of climate factors and soil factors on ΔSOC. For more comprehensive explanations, please refer to the Supplementary discussion.

To demonstrate the advantages of high-resolution carbon budget quantification, we produced a 0.5-degree-resolution ΔSOC estimation (Fig. 5d) and conducted a comparative analysis with the 0.0025-degree-resolution estimation (Fig. 5e, f). Specifically, we employed KGML-ag-Carbon at a 0.5-degree-resolution and applied the mass balance approach to derive ΔSOC using 0.5-degree estimates of NEE and crop yield (Fig. S15; Fig. 5d). To achieve this, we employed a mean aggregation approach for each input variable, converting from 250-meter resolution to 0.5-degree resolution. However, the aggregation of crop types from high resolution to coarse resolution was not straightforward. To better emulate realistic crop rotations in the coarse-resolution simulation, we conducted two simulations involving corn–soybean rotations (corn in even years and soybean in odd years) and soybean-corn rotations (soybean in even years and corn in odd years). Subsequently, we used the corn/soybean fractions from CDL and CSDL data to compute weighted averages of corn and soybean estimations for each year using those two simulations. The differences between the 0.5-degree-resolution and 0.0025-degree-resolution ΔSOC estimations were then calculated (Fig. 5e, f) by subtracting the 0.0025-degree-resolution estimation from the 0.5-degree-resolution estimation. To enhance clarity, regions outside the corn/soybean agroecosystem were excluded, and the estimated differences were converted to percentage fractions using methods similar to those employed in generating Fig. 5b–c. We have adopted the normalized root mean square error (NRMSE) to describe the overall differences:
NRMSE = RMSE / (Q3 - Q1) (12)
Where RMSE is the root mean square error between 0.5-degree- and 0.0025-degree-resolution estimations, and Q1 and Q3 represent the three quantiles and one quantile of 0.0025-degree-resolution estimation, respectively. This method was chosen to avoid the denominator becoming too small. We also applied a similar approach to calculate the differences in other variables, including GPP, Ra, Rh, NEE, and Yield (Fig. S15). We note that the relative differences in Fig. S15 were calculated as (0.5-degree-resolution estimation - 0.0025-degree-resolution estimation)/0.0025-degree-resolution estimation, representing the relative differences relative to each 0.0025-degree pixel.

We conducted an extensive literature review to gather available soil organic carbon (SOC) measurements in the U.S. Midwest. This effort involved data from 18 sites, each with multiple SOC measurements at the plot level (~10 m) after 2000, facilitating ΔSOC validation for the KGML-ag-Carbon model (Fig. S16, Table S1). Observed data such as bulk density, initial SOC stock in the top 30 cm, and rotation management were integrated into the input feature when applicable. Other inputs were directly derived from our 250-m-resolution regional database, such as weather forcings and GPP based on the sites’ geophysical locations. We have used an empirical equation72 to simulate the percentage fraction of SOC at different depths to total stock (assumed to be SOC in 0–100 cm), expressed as:
FSOC,Z = -0.011Z2 + 2.029*Z (13)
Where F is the estimated SOC percentage between 0 to Z cm depth. This conversion factor aided in translating observed SOC values to the entire profile or to the top 30 cm in cases where depth-specific data was unavailable. It is worth noting that all of the collected ΔSOC data pertain to plot-level (~10 m) experimental measurements that primarily focus on detecting the influences of management practices. Data from those plots often lack the requisite localized forcing data needed by our model. Consequently, we resort to utilizing field-level (250 m) forcings such as remotely sensed GPP and reanalysis of NLDAS weather forcing, which poses a scale mismatch when compared to the plot-level observations. To illustrate this scale mismatch, we selected two sites from different studies36,38 to compare the sizes of the experimental plots with the sizes of our predictions and neighboring real fields, as depicted in Fig. S16a, b.

Development environment description
We used Pytorch 1.6.0 (https://pytorch.org/get-started/previous-versions/, last access: 21 Oct 2023) and Python 3.7.11 (https://www.python.org/downloads/release/python-3711/, last access: 21 Oct 2023) as the programming environment for model development. Statistical analysis, such as linear regression, was conducted using Statsmodels 0.14.0 (https://github.com/statsmodels/statsmodels/, last access: 21 Oct 2023) In order to use a GPU to speed-up the training process, we installed the CUDA Toolkit 10.1.243 (https://developer.nvidia.com/cuda-toolkit, last access: 21 Oct 2023). A desktop with an NVIDIA 2080 super GPU was used for code development and testing. The training processes, which required extensive time and memory space, were conducted on the Mangi and Agate clusters (https://www.msi.umn.edu/mangi, last access: 21 Oct 2023) from the High-Performance Computing facility of the Minnesota Supercomputing Institute (HPC-MSI, https://www.msi.umn.edu/content/hpc, last access: 21 Oct 2023) with two-way NVIDIA Tesla V100 GPUs.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript \"Knowledge-based artificial intelligence significantly improved agroecosystem carbon cycle quantification\" uses multiple model structure and training strategies to achieve better carbon cyle quantification than purely-PB or purely-ML approaches. The authors demonstrate the relative contributions of various techniques to model performance and demonstrate the economic/policy value of high-resolution carbon flux information. In my opinion the novelty of this work is in the systematic evaluation of various forms of knowledge guidance and the connection of quantitative improvements in predictions to the carbon market.\nAs far as I can ascertain, the methodology and interpretation of results are sound. While the large number of training strategies can become somewhat convoluted, the authors justify each variation and make attempts to summarize clearly (e.g. Table 1).\nThe authors make several references to high-resolution carbon flux and crop yield predictions. It would be helpful to have the resolution contextualized earlier in the paper. I recognize this is covered to some extent in the discussion of basis risk.\nI believe this application of KGML will be interesting in many application fields beyond agroecosystem studies and therefore read by non-experts. As such, I think the paper would be improved by providing more explanation of the major carbon cycle components and building more intuition about how they interact with each other.\nMinor concerns:\nI suggest removing language like \"remarkable\" (L24) and \"outstanding\" (L174)\nCan you comment on why there is such significant spread in the ensemble of Trendy models?", "reviewer_2": "This work demonstrates the usefulness and applicability of combining process-based models with the data-adaptiveness of machine learning methods. The authors go an interesting route by simply pretraining on outputs from a process-based agricultural model, ecosys, to embed knowledge into a machine-learning framework. Additionally, the authors embed knowledge in constructing the machine learning model by considering the structure of their studied system by sensibly stacking stateful ML models. Furthermore, they designed so-called knowledge-guided loss functions to incorporate knowledge into the training process further.\nThe authors have developed, trained and tested this framework for modelling fine-scale carbon budgets and crop yields in the Midwest of the USA. They have done detailed tests and their line of argumentation is overall easy to follow. Overall, the authors make a sound case for their KGML-ag-Carbon framework.\nI see four areas where I would like to see substantial improvements and clarifications:\n1) Claim ‘help save over 500 million US$ per year for stakeholders’\nIt is unclear to me who these stakeholders are. Insurances? Would they overpay farmers because of coarse-scale carbon-budget accounting? Or are these stakeholders farmers who cannot sell their carbon sequestration services to the market because of a lack of reportability? My question is, who would lose 500 million US$, and who would gain 500 million US$ in your calculations? Would climate change mitigation be affected by this 500 million US$ example because funds are not going into something cost-effective?\n2) Realism of ΔSOC\nThe calculation of carbon credit basis risk relies on the simulation of NEE and yields since the authors define ΔSOC as NEE – yield. In my opinion, the readers should see histograms and maps of ΔSOC and how they compare with current SOC stocks, so we can assess the realism of these ΔSOC predictions. The authors should be able to compare their ΔSOC with repeated SOC measurements in this area, ideally at the flux towers. To me, this is paramount to avoid carbon double-accounting and ensure SOC is increasing. One point that the authors could discuss more thoroughly is the calculation of harvest residues that should remain on the field. Is it accurate to assume that harvest residues that remain on the field are GPP – Ra – yield? Is there uncertainty regarding where the residues go?\n3) Static SOC\nThe KGML-ag-carbon has TSOC as a static input feature (chapter 5.3). This is problematic since that way the KGML-ag-carbon model is not dynamic, i.e. the essential variable the model predicts, ΔSOC, does not feedback on SOC. Instead, the initial SOC content is treated as a static, not changing variable. First, this is a missed opportunity which the authors could easily remedy in KGML-ag-carbon by introducing a recurrent cell which updates as SOC[t+1] = SOC[t] + ΔSOC[t]. Second, by explicitly outputting SOC concentrations, it would be easy to compare the outputs of the model with the ground truth, i.e. if SOC contents are indeed increasing in these fields.\n4) The trade-off with climate costs of N fertilization and N2O\nWith the claim that stakeholders can save 500 million $ in the climate mitigation market, the authors neglect the nitrogen cost for permanently storing carbon as soil organic carbon. Harvest residues are generally depleted in nitrogen compared to the harvest product of interest, especially soil microbes, which decompose harvest residues to SOC with a carbon use efficiency that depends on the stoichiometric imbalance between harvest residues and soil microbes. Nitrogen inputs are needed for a high carbon use efficiency. The energy costs of these nitrogen fertilizers would need to be considered for these carbon credits. Similarly, KGML-ag-carbon would reward management practices that enhance carbon storage but would not consider a potential trade-off of high fertilization with high N2O emissions. The authors seem to have worked on N2O separately in 'KGML-ag: A Modeling Framework of Knowledge-Guided Machine Learning to Simulate Agroecosystems: A Case Study of Estimating N2O Emission using Data from Mesocosm Experiments' in GMD. In my opinion, it is very important to avoid incremental publishing, especially in high-impact journals such as Nature Communications and with claims like the 500 million $ costs saved for stakeholders, which may come from false accounting when neglecting N2O. The climate benefits of management practices have to be reported here, not financial implications for one group of stakeholders.\nMiscellaneous issues and questions\n• Why is an attention layer introduced for yield predictions? Can you explain why GRU modules were not used for this task?\n• Is mass conservation encoded into the model? Is GPP = Yield + Ra + Rh + NEE encoded or learned?\n• In Figure 2, the authors show how the number of counties and flux towers influences the performance of the different models – it would be good also to have the number of observations associated with each sample, i.e. the number of annual or daily measurements per training set.\n• For Figure 3, it would be better to have a full-factorial test of the contributions of GPP Data, Pretrain data from ecosys, and KG Loss. Maybe the authors did this already, and I misunderstood. For example, how well would a model perform if we combined ML with a knowledge-guided loss function? It seems the authors only added them in the last step.\n• The authors compare KGML-ag-carbon with potential runtimes of ecosys. However, they compare GPU and CPU runtimes. Nowadays, process-based models can be run in parallel on GPUs – the comparison may be unfair.\n• The GPP product SLOPE is not a pure remote-sensing product (line 241) but has machine learning models included – it would be good to note that down for the reader\n• While the methods are described in great detail in the supplement it would be better to provide the code associated with this study to ensure reproducibility. There is also no mention of which framework was used to code KGML-ag-carbon and what kind of GPU was used for training.\nSome typos\nLine 110: cheapercompared – space is missing\nLine 656: Fig. SX – number is missing", "reviewer_3": "In this manuscript, the authors design an AI framework, Knowledge Guided Machine Learning (KGML), to calculate carbon capture in agriculture sites using data from an ecosystem model, ecosys, to pretrain the KGML and observations of yield and eddy covariance to fine tune the model. They perform a variety of experiments to examine the performance sensitivity to the size of the training data and compared the results with a simple ML model. They also performed additional experiments to explore alternative training, out-of-sample and extreme year data, and different data-splitting to understand how to reduce uncertainty. The KGML method outperformed the ecosystem model and pure ML model for capturing yields and EC fluxes. Finally, the authors use this approach to estimate carbon capture in agriculture fields. I have a few suggestions and questions for the authors to address, but overall, I find this paper to be an interesting new application of machine learning that can have useful approaches for a variety of biogeochemical modeling.\nI’m not an expert in AI so I can’t comment much on that element. I did get a little lost in some of the details of the methods since I am not familiar with KGML, but I was able to understand overall what the authors did. I liked how the authors essentially did a sensitivity study of the approach, comparing with pure ML, and against the ecosystem model itself. One additional suggestion would be to look at how the method performs with pretraining from a different ecosystem model. I think it outside the scope of this manuscript, and the authors discuss the importance of ecosystem model choice in Section 4.2. I wonder how well the ecosys model works outside the U.S. or with other crop types, such as wheat or rice. This discussion might be expanded more, especially since the ecosys model doesn’t perform great compared with observations.\nI’m not convinced that the model would be accurate enough to estimate carbon capture in an agriculture field given the lack of validation data. Discussion section 4.1 seems out of place and unnecessary.\nThe field sites where data is collected (particularly the flux sites) don’t have much information about them in the manuscript. All those plots have different crops and management practices (e.g., rotation, tillage, fertilizer application, irrigation, etc.). The reader doesn’t find out until section 4.3 that these aren’t taken into account in the model. While I like that the authors suggested a future application could be to quantify management, it should be emphasized that management practices are strong drivers of carbon inputs. There should be some data at these sites that can be used as model input or at least include a comment on how that will influence the results.\nAnother consideration for this model is that ML can only learn from the past and extreme events are often not included but they can be important for future estimates of carbon fluxes. For example, carbon fertilization is generally not included in ecosystem models, nor is crop response to extreme heat or flash droughts. Including these in ML models could improve our prediction of crop models.\nThe methods section indicates that most of the data is used for pre-training and training purposes. Do the authors worry about overfitting the model? Also, what is the difference between “easier” and “harder” samples (Line 623)\nIn the abstract, line 23, do the authors mean remote data instead of satellite?\nLine 110: cheapercompared should be two words", "reviewer_4": "In this work, the authors proposed a Knowledge-Guided Machine Learning (KGML) framework and demonstrated that KGML can outperform existing process-based or ML models in quantifying carbon cycles. In general, none of the key parts in KGML framework (i.e., knowledge-guided loss, synthetic data pre-training, the use of indirect data) makes a novel contribution to existing machine learning techniques. Although the combined use of all those parts hasn’t been explored by existing works yet, this manuscript presents an overall performance of KGML without clearly studying the contribution of each part (e.g., whether the sequence of including different parts would impact the contributions in Figs. 3ab, why knowledge-guided loss always has a very small contribution in Fig. 3). In addition, the descriptions of some key parts (e.g., knowledge-guided losses, and the use of indirect data) are not clear. In terms of quantifying carbon cycles, this manuscript does not provide a clear description of the case studies and the general motivation for studying them. The quality of the figures in the manuscript should be improved as well. Therefore, the reviewer cannot support the publication of this manuscript in Nature Communications."}, "author_response": "Pleasebeawareoftheformattingofallresponses:\n1.Reviewercommentinblack,responseinblueandquotationfromthemaintextinred;\n2.Thereferredsectionsarebasedonthecleanversionoftherevisedmanuscript,nottheoldortrack\nchangeversion.\nThe manuscript \"Knowledge-based artificial intelligence significantly improved agroecosystem carbon\ncycle quantification\" uses multiple model structure and training strategies to achieve better carbon cycle\nquantification than purely-PB or purely-ML approaches. The authors demonstrate the relative\ncontributions of various techniques to model performance and demonstrate the economic/policyvalueof\nhigh-resolution carbon flux information. In my opinion the novelty of this work is in the systematic\nevaluation of various forms of knowledge guidance and the connection of quantitative improvements in\npredictionstothecarbonmarket.\nAs far as I can ascertain, themethodologyandinterpretationofresultsaresound.Whilethelargenumber\nof training strategies can become somewhat convoluted, the authors justify each variation and make\nattemptstosummarizeclearly(e.g.Table1).\nResponse: We greatly appreciate these positive comments onourmanuscript,particularlytherecognition\nthat our systematic evaluation of various forms of knowledge guidance is novel, as well as the\neconomic/policyvalueofhigh-resolutioncarbonfluxinformation.\nTheauthorsmakeseveralreferencestohigh-resolutioncarbonfluxandcropyieldpredictions.Itwouldbe\nhelpful to have the resolution contextualized earlier in the paper. I recognize this is covered to some\nextentinthediscussionofbasisrisk.\nResponse: Thank you for this suggestion, which we fully agree with. We have revised the second\nparagraph of the Introduction to emphasize the need for accurate, cost-effective, and high-resolution\ncarbon quantification methods, particularly at the field level. These revisions should help better\ncontextualizetheresolutionofourmethodinthepaper.Therevisedparagraphisbelow:\n“Increasing agricultural carbon sequestration is a key strategy for mitigating climate change, and\nsignificant effortsandinvestmentshavebeenmadeintheUSandacrosstheglobetoimplementprograms\nthat incentivize SOC enrichment7,8. It is thus more important than ever to develop robust and scalable\nmethods for reliably quantifying field-level carbon sequestration, both to assess the climate mitigation\neffect and to ensure that mitigation actions by individual farmers are compensated fairly and accurately.\nHowever, traditional carbon quantification methods, which rely on soil sampling, emission factors, and\nprocess-based (PB) modeling, encounter inherent challenges that hinder their ability to achieve the\nrequired levels of accuracy, scalability, and cost-effectiveness9–11. … Therefore, new methodsareneeded\nto overcome the limitations of PB and ML models, enabling cost-effective, accurate, and interpretable\nmeasurement and monetization of carbon outcomes at the individual field level, to reduce errors in\naggregatedquantificationsandpromotemoresustainablelandmanagementpractices12,19.”\nI believe this application of KGML will be interesting in many application fields beyondagroecosystem\nstudies and therefore read by non-experts. As such, I think the paper would be improved by providing\nmore explanation of the major carbon cycle components and building more intuition about how they\ninteractwitheachother.\nResponse: Thank you for this suggestion. We agree that it's important to provide more context and\nelaborate on carbon cycle components and their interactions. Thus, we have added more description to\nsection 2, paragraph 2, covering carbon inputs/outputs of agroecosystems, respiration types, and\n2\ndisturbance impacts to enhance the accessibility to non-experts and clarify the significance of\nKGML-ag-Carbon.Therevisedcontentisbelow:\n“KGML-ag-Carbon resolves themajorcarbonbudgetcomponents,includingautotrophicrespiration(Ra),\nheterotrophic respiration (Rh), total ecosystem respiration (Reco, Ra + Rh) and net ecosystem carbon\nexchange (NEE) on a daily scale, and yield on an annual scale. As in natural ecosystems, changes in\nagroecosystem soil carbon storage are determined by the mass balance of input and output carbon\nfluxes30,31. Carbon inputs originate from plantphotosynthesis,i.e.,grossprimaryproduction(GPP),while\ncarbon inputs to the soil include both aboveground and belowground litter and root exudates. Carbon\noutputs occur through respiration, including Ra from plant shoots and roots, and Rh from SOC\ndecomposition by microbes and fungi. Disturbances suchasharvestingmayalsoremovecarbonfromthe\necosystem periodically. Based on the carbon fluxes and yield estimated from KGML-ag-Carbon, annual\nchanges in SOC can be determined using the mass balance equation ΔSOC=-NEE-Yield12,13,32.Detailed\ninformationaboutthestructuraldevelopmentofKGML-ag-CarbonisgiveninMethods5.3.”\nMinorconcerns:\nIsuggestremovinglanguagelike\"remarkable\"(L24)and\"outstanding\"(L174)\nResponse:Followingthissuggestion,wehaveremovedsuchexpressionsfromtherevisedmanuscript.\nCanyoucommentonwhythereissuchsignificantspreadintheensembleofTrendymodels?\nResponse: The significant spread intheTrendymodelsensembleismainlyduetotheinclusionofdiverse\nprocesses and alternative parameterizations in the models from different research communities, as\ndescribed by Sitch et al. (2015). More details about the parameterizations and processes of the Trendy\nmodels are described in Table S1 of the supplementary material for that paper:\nhttps://bg.copernicus.org/articles/12/653/2015/bg-12-653-2015-supplement.pdf. We have also added a\nbriefexplanationinsection5.7:\n“... ThewiderangeofvariationobservedintheTrendymodelsensemblecanbeattributedtotheinclusion\nof diverse processes and alternative parameterizations adopted by models from different research\ncommunities,asdescribedbySitchetal.(2015)15.”\nReference:\nSitch,S.,Friedlingstein,P.,Gruber,N.,Jones,S.D.,Murray-Tortarolo,G.,Ahlström,A.,Doney,S.C.,\nGraven,H.,Heinze,C.,Huntingford,C.,Levis,S.,Levy,P.E.,Lomas,M.,Poulter,B.,Viovy,N.,\nZaehle,S.,Zeng,N.,Arneth,A.,Bonan,G.,…Myneni,R.(2015).Recenttrendsanddriversof\nregionalsourcesandsinksofcarbondioxide.Biogeosciences,12(3),653–679.\nThis work demonstrates the usefulness and applicability of combining process-based models with the\ndata-adaptiveness of machine learningmethods.Theauthorsgoaninterestingroutebysimplypretraining\non outputs from a process-basedagriculturalmodel,ecosys,toembedknowledgeintoamachine-learning\nframework. Additionally, the authors embed knowledge in constructing the machine learning model by\nconsidering the structure of their studied system by sensibly stacking stateful ML models. Furthermore,\nthey designed so-called knowledge-guided loss functions to incorporate knowledge into the training\nprocessfurther.\nThe authors have developed, trained and tested this framework for modelling fine-scale carbon budgets\nand crop yields in the Midwest oftheUSA.Theyhavedonedetailedtestsandtheirlineofargumentation\nisoveralleasytofollow.Overall,theauthorsmakeasoundcasefortheirKGML-ag-Carbonframework.\nIseefourareaswhereIwouldliketoseesubstantialimprovementsandclarifications:\n1)Claim‘helpsaveover500millionUS$peryearforstakeholders’\nIt is unclear to me who these stakeholders are. Insurances? Would they overpay farmers because of\ncoarse-scale carbon-budget accounting? Or are these stakeholders farmers who cannot sell their carbon\nsequestration services to the market because of a lack of reportability? My question is, who would lose\n500 million US$, and who would gain 500 million US$ in your calculations? Would climate change\nmitigation be affected by this 500 million US$ example because funds are not going into something\ncost-effective?\nResponse: The stakeholders we refer to are participantsinthecarbonmarket,whichincludefarmerswho\nadopt new practices, carbonprojectdevelopers(e.g.companiesorgovernments)thatpayfarmerstomake\nchanges and collect generated carbon credits,andtheentitiesthatpurchasecarboncredits.Overpayments\nto farmers represent a loss of money that could have funded more carbon removal, whereas\nunderpayments to farmers risk farmers not being able to recover the cost of implementing carbon\nsequestration practices and being discouraged from continuing them. We calculated the totallosstoboth\ncarboncreditbuyersandfarmers,whichrepresentsmisallocatedfunds.\nMoreover, our main point in mentioning the potentiallossistoemphasizetheimportanceofaccurate\ncarbon quantification in making cost-effective decisions for climate mitigation. The overall goal of\nclimate mitigation is to reduce greenhouse gas emissions and support carbon sinks, which requires a\nmulti-faceted approach and careful allocation of limited resources. Inaccurate carbon quantification can\nlead to an inefficient carbon market, which may result in wasted funds that could be better invested in\nmore cost-effective solutions. Therefore, we believe thatimprovingtheaccuracyandreliabilityofcarbon\nquantification is crucial for optimizing investments and achieving our collective climate goals. We have\naddednewtexttoemphasizethispointinthesecondparagraphoftheintroduction:\n“Increasing agricultural carbon sequestration is a key strategy for mitigating climate change, and\nsignificant effortsandinvestmentshavebeenmadeintheUSandacrosstheglobetoimplementprograms\nthat incentivize SOC enrichment7,8. It is thus more important than ever to develop robust and scalable\nmethods for reliably quantifying field-level carbon sequestration, both to assess the climate mitigation\neffect and to ensure that mitigation actions by individual farmers are compensated fairly and accurately.\nHowever, traditional carbon quantification methods, which rely on soil sampling, emission factors, and\nprocess-based (PB) modeling, encounter inherent challenges that hinder their ability to achieve the\nrequired levels of accuracy, scalability, and cost-effectiveness9–11. … Therefore, new methodsareneeded\nto overcome the limitations of PB and ML models, enabling cost-effective, accurate, and interpretable\nmeasurement and monetization of carbon outcomes at the individual field level, to reduce errors in\naggregatedquantificationsandpromotemoresustainablelandmanagementpractices12,19.”\nWe have removed the statement regarding the 500 million US$ savings to acknowledge the concern\nthat this finding is based on scenario simulations with many uncertainties. Instead, we have revised the\nabstract and added discussion in section 4.1 to emphasize the uncertainties in both observation and\nsimulationofsoilorganiccarbon:\n“The prediction improvement indicates that over 70% of the error in NEE and ΔSOCestimationscanbe\nreduced.”\n2)RealismofΔSOC\nThe calculation of carbon credit basis risk relies on the simulation of NEE and yields since the authors\ndefine ΔSOC as NEE – yield. In my opinion, the readers should see histograms and maps of ΔSOCand\nhow they compare with current SOCstocks,sowecanassesstherealismoftheseΔSOCpredictions.The\nauthors should be able to compare their ΔSOC with repeated SOC measurements in this area, ideally at\nthe flux towers. To me, this is paramount to avoid carbon double-accounting and ensure SOC is\nincreasing. One point thattheauthorscoulddiscussmorethoroughlyisthecalculationofharvestresidues\nthat should remain on the field. Is it accurate to assume that harvest residues that remainonthefieldare\nGPP–Ra–yield?Isthereuncertaintyregardingwheretheresiduesgo?\nResponse: We fully agree that validation of ΔSOC is crucial to ensure the use of this model for carbon\naccounting. The validation suggested by the reviewer would be an ideal approach, but one that is not\ncurrently feasible due to limitations in existing SOC data. An exhaustive literature search revealed three\nmainreasonsforthis:\n(i) SOC measurements have high uncertainty from various sources, including spatial heterogeneity,\nsampling design, variability in bulk density, and variation in soil processing methods and laboratory\nassays (Goidts et al., 2009; Stanley et al., 2023), which together could easily lead to uncertainty\nexceeding 50% in field-level SOC stockestimation,evenovermulti-yeartimescales.Theuncertainty\nof annual ΔSOC is even larger because the magnitude ofchangeisrelativelysmall.Thisislikelythe\nreasonwhythemostaccessibleSOCdataintheUSarestaticmeasurements(Batjesetal.,2020).\n(ii) Scale mismatch between point sampling and field-level estimation. Our literature search also found\nthat most ΔSOC measurements are from plot-level experiments (Al-Kaisi & Kwaw-Mensah, 2020;\nIbrahim et al., 2018; Jin et al., 2015; Khan et al., 2007; Olson et al.,2014;Poffenbargeretal.,2017;\nSchmer et al., 2014; Varvel, 2006; Venterea et al., 2006), which are not representative of the\nfield-level values needed for ΔSOC validation. There’re some dense sampling efforts for field-level\nSOC benchmarking by private companies, but unfortunately, data are not available for public\nresearch. A discussion of soil sampling methods that can be used for benchmarking field-level SOC\nchangeisgivenbyGuanetal.(2022).\n(iii) Other key variables needed to verify model performance are also missing. We have attempted to\ndownscale the KGML-ag-Carbon model from field scale to plot scale to validate the ΔSOC using\nplot-level data. However, such experiments often do not report precise latitude and longitude\n5\ncoordinates that are needed to extract remotely sensed GPP data to constrain the KGML-ag-Carbon\nmodel, or have too small plot size so that even the 30m resolution SLOPE-GPP dataistoocoarseto\ndistinguishbetweenplots.\nFor the KGML-ag-Carbon model, we therefore adopted a mass balance approach to indirectly\nestimate ΔSOC from NEE and yield, which have been well-validated against available individual\ncomponents. The ΔSOC estimates integrate the best available data, including static SOC, weather\nforcing, soil properties,croptypeandremotelysensedGPP.However,weacknowledgethatthisapproach\nhas limitations: Ra and Rh cannot be separately validated, as they are seldom partitioned in flux tower\nsites, and the influences of different management practices are not directly considered, which could\nfurther affect the accuracy of the calculated ΔSOC. This type of situation is a commonchallengeforthe\nmodeling community; thus, we emphasize the need for more comprehensive measurements of different\ncomponentsofthecarboncycle.\nSimilarly, we acknowledge that predicting residue fates using GPP-Ra-Yield alsohaslimitationsdue\nto the lack of direct observations of residues and Ra, and the absence of management data regarding\nresidueremovalfromthefield.ThecropresidueestimatedintheResidue_layerofKGML-ag-Carbonwas\nonly used as a prior information input into the GRU_Rh submodule. After pretraining using\necosys-generated synthetic data, KGML-ag-Carbon is able to determine the relationship between\nGPP-Ra-Yield and actual residue when there are disturbances like removal. The Rh was further\nconstrained together with Ra by Reco, which is derived from EC flux-tower measurements. While\npartitioned Rh data fromRecowouldbeamoreidealconstraint,suchpartitioneddataisoftenlimitedand\navailableonlyatcoarsetemporalresolution(seasonally~annually)(Jianetal.,2021).\nBecause we don't have direct ΔSOC validation, as pointed outbythereviewer,wehaveremovedthe\nbasis risk calculation using ΔSOC that involves policy-sensitive US$ values. Instead, we have added\ndiscussion on the need for accurate SOC measurements and focus on the benefits of high-resolution\nquantification of NEE (whichiswell-validatedinthisstudy)andΔSOC.Wehavealsoaddedadiscussion\nof the uncertainty of themodel’intermediatevariablesandtheneedformeasurementstoconstrainthem.\nWe have updated Section 4.1, Fig. 5 (shown here as Fig. R3), Fig. S13 (with the removalofpanelcthat\ndisplays basis risk), and Methods 5.8 to reflect these changes. While we recognize the importance of\naccurate field-level SOC measurements for validating carbon accounting and the need formeasurements\nof Ra, Rh, and crop residue to constrain the underlying processes, we believe that our study provides\nvaluable insights into the potential of the extensible KGML-ag-Carbon framework for reducing carbon\nbudgetquantificationerrors,whichcanaidfuturecarbonaccountingandclimatemitigationefforts.\nTherevisedFig.5(hereasR3)is:\nusingcoarseresolutionforcarbonbudgetestimations.(c-d)ΔSOCabsoluteerrorscausedbyusingcoarseresolution.\n0.0025-degreeresolutionistreatedas“field-level”groundtruth.\nReviseddiscussionsection4.1islistedbelow:\n“The field-level quantification of carbon budgets, crop yields, and ΔSOC produced using\nKGML-ag-Carbon (as demonstrated for the US Midwest) provides an accurate, cost-effective, and\nhigh-resolution product for potentially improving carbon sequestration assessments. We used the term\nabsolute error to quantify the advantages of high-resolution carbon budgets generated by\nKGML-ag-Carbon. Previous approaches often assess the carbon budget using factors that are averaged\nover a certain region, and absolute error describes the potential consequences resultingfrommismatches\nbetween local-scale factors and the regional baseline. As used here, absolute error represents the\ndifference of NEE or ΔSOC that can be attributed to individual fields using a high-resolution carbon\nbudget product compared to larger-scale estimates based on lower-resolution products (Methods 5.8).\nGenerally,theabsoluteerrordecreasesquicklywithincreasingresolutionofthecarbonbudgetproductfor\nthe US Midwest (Fig. 5 and Fig. S13). For example, when using a 1x1°-resolution product, the\nuncertainties in estimated NEE and ΔSOC of individual fields are up to 90 gC/m2/year and 101\ngC/m2/year (75% quantile), respectively(Fig.5).Basedonthiscalculation,thoseabsoluteerrorswilllead\nto a total 5.6 and 6.2 GgC/year (1 Gg = 109 g) errors for NEE and ΔSOC quantification,respectively,in\nthe US Midwest. However, if using a product with a relatively high resolution of 0.005°, maximum\n7\nestimate uncertaintiesdecreaseto27gC/m2/yearforbothNEEandΔSOC(75%quantile).Thetotalerrors\nwould be reduced by 70% and 73%purelybychangingtheresolutionofNEEandΔSOCquantifications,\nrespectively (Fig. 5). More detailed results about the benefits of high-resolution NEE, crop yield, and\nΔSOC qualifications are given in Fig. S13. To be noted, the KGML-ag-Carbon employs a mass balance\napproach to estimate ΔSOC from NEE and yield, which are estimated by integrating all available data,\nincluding weather forcing, soil properties (which include static SOC), crop type, and remotely sensed\nGPP, and are well-validated by observations. This approach allows us to make the best use of existing\ndata to estimate the regional ΔSOC at low cost and high resolution, even in the absence of field-level\nmeasurements. However, it is important to acknowledge that the lack of validation using field-level\nmeasurements of ΔSOC limits the certainty of our estimates, as such measurements are oftenscarceand\nassociated with large uncertainties. For instance, at 250-m scales, uncertainty in SOC arises from lab\nmeasurement error (up to 12%),spatialsamplingerror(upto50%),andresamplingerror(upto45%)35,36.\nThe uncertainty of ΔSOC over a long period of time will be even larger. Including measurements with\nsuch uncertaintiesintheKGML-ag-Carbontrainingdatamaydecreasethereliabilityofothercomponents\nsuch as NEE and yield, which have relatively lowermeasurementerrors.It'salsoworthnotingthatwhile\nthe NEE, Reco, and crop yield in the KGML-ag-Carbon are well-constrained, the intermediate variables\nlike Ra, Rh, and crop residue may stillcontainhighuncertaintyduetoalackofdirectobservationaldata.\nThese variables, however,arecrucialforabetterunderstandingoftheunderlyingmechanisms.Therefore,\nwhile theKGML-ag-Carbonoffersafeasiblemethodforreducingcarbonbudgetquantificationerrorsdue\nto its high resolution (250 m), high accuracy (Fig.2-4),andlowcomputationalcost,italsohighlightsthe\nneed for accurate field-level ΔSOC measurements to improve the reliabilityofΔSOCquantification,and\ntheneedforaccuratemeasurementsofRa,Rh,andcropresiduetoconstraintheunderlyingprocesses.”\nInaddition,wehavealsorevisedMethodssection5.8:\n“5.8Investigatingthebenefitsofhigh-resolutioncarbonbudgetquantification\nWe used the absolute error to quantify the uncertainties in NEE and crop yield quantification for\nindividual fields in coarse-resolution carbon budget estimations (Fig. 5, Fig. S13). Specifically, we first\nregridded the 250-m-resolution products fromKGML-ag-Carbonintoa0.0025-degreeproducttobeused\nas a baseline for field-level carbon credit quantification Then we aggregated the 0.0025-degree products\nto various coarse resolutions, including 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, and 1 degrees. The\ncoarse-resolution products were compared with the0.0025-degreeproducttocalculatetheNEEandyield\ndifferencesasabsoluteerrors.”\n3)StaticSOC\nThe KGML-ag-carbonhasTSOCasastaticinputfeature(chapter5.3).Thisisproblematicsincethatway\nthe KGML-ag-carbon model is not dynamic, i.e. the essential variable the model predicts, ΔSOC, does\nnot feedback on SOC. Instead, the initial SOC content is treated as a static, not changing variable.First,\nthis is a missed opportunity which the authors couldeasilyremedyinKGML-ag-carbonbyintroducinga\nrecurrent cell which updates as SOC[t+1] = SOC[t] + ΔSOC[t]. Second, by explicitly outputting SOC\nconcentrations, it would be easy to compare the outputs of the model with the ground truth, i.e. if SOC\ncontentsareindeedincreasinginthesefields.\nResponse: As suggested, we conducted extra simulation experiments using a recurrent cell to update\nSOC[t+1] as SOC[t] + ΔSOC[t] in the KGML-ag-Carbon model and tested this new approachwithboth\n8\nsynthetic and observed data. However, the results indicate that including a predicted variable with large\nuncertainty as an input in the model reduces the prediction performance of target variables in synthetic\ndata tests, especially those that are most relevant, such as Rh (Fig. R4). Additionally, adding dynamic\nSOC either has no significant influence on performance in real observations or slightly reduces the\nperformanceofRecoandNEE(Fig.R5).\nFig. R4. Performance comparison between KGML-ag-Carbon models using synthetic data with static SOC input\nfrom the gSSURGO database and dynamic SOCinputfromupdatingstaticSOCwithΔSOCpredictedfromthelast\ntimestep.\nFig. R5. Performance comparison between KGML-ag-Carbon models using observed data with static SOC input\nfrom the gSSURGO database and dynamic SOCinputfromupdatingstaticSOCwithΔSOCpredictedfromthelast\ntime step. Models were trained using either 5 or 40 counties (1 or 7 sites) for yield (flux) optimization and tested\nusingout-of-sampleyielddatafor 210countiesfromNASS(theout-of-samplegroupofsites).\nforcing (TF) approach in ML. However, this method has limitations in modeling time series, as it\nprimarily focuses on learning ΔSOC[t], hindering the model's ability to accurately accumulate SOC\nchanges over time. Additionally, due to temporal dependencies, ΔSOC[t] is often small, causing\nTF-trained models to overly rely on SOC[t] for predicting SOC[t+1]. In the testing phase, errors in\npredicting SOC[t] cannot be corrected, leading to error accumulation and poor SOC predictions. More\ndetailed explanation of the limitations of TF can be found in section 7.8 of Xu et al. (2023). Moreover,\nwhile using a recurrent cell to update SOC[t+1] is an ideal solution when sufficient accurate SOC\nmeasurements areavailableasconstraints,wecannotcurrentlydosoatthelandscapelevelduetothelack\nof reliable measurements of SOC change. Therefore, we choose to retain SOCasastaticinputfeaturein\nthisstudy.\n4)Thetrade-offwithclimatecostsofNfertilizationandN2O\nWith the claim that stakeholders can save 500 million $ in the climate mitigation market, the authors\nneglect the nitrogen cost for permanently storing carbon as soil organic carbon. Harvest residues are\ngenerally depleted innitrogencomparedtotheharvestproductofinterest,especiallysoilmicrobes,which\ndecompose harvest residues to SOC with a carbon use efficiency that depends on the stoichiometric\nimbalance between harvest residues and soil microbes. Nitrogen inputs are needed for a highcarbonuse\nefficiency. The energy costs of these nitrogen fertilizers would need to be considered for these carbon\ncredits. Similarly,KGML-ag-carbonwouldrewardmanagementpracticesthatenhancecarbonstoragebut\nwould not consider a potential trade-off of high fertilization with highN2Oemissions.Theauthorsseem\nto have worked on N2O separately in 'KGML-ag: A Modeling Framework of Knowledge-Guided\nMachine Learning to Simulate Agroecosystems: A Case Study of Estimating N2O Emission using Data\nfrom Mesocosm Experiments' in GMD. In my opinion, it is very important to avoid incremental\npublishing, especially in high-impact journals such as Nature Communications and with claims like the\n500 million $costssavedforstakeholders,whichmaycomefromfalseaccountingwhenneglectingN2O.\nThe climate benefits of management practices havetobereportedhere,notfinancialimplicationsforone\ngroupofstakeholders.\nResponse: Thank you for the valuable discussion. We have revised the discussion in section4.1tofocus\non the NEE and ΔSOC regional results, andremovedthelesswell-validatedbasisriskcalculationrelated\nto ΔSOC. Additionally, we have added new contentinsection4.3todiscussthelimitationsofthecurrent\nmodelandpotentialsolutionsforintegratingC-NinteractionsinfutureversionsofKGML-ag-Carbon.\nWe acknowledge that management practices that aim to enhance carbon storage may unintendedly\nincrease N fertilizer use, especially when new practices reduce yields,therebyincreasingN2Oemissions\nand partially offsetting climate mitigation effects. In practice, this is less of a concern because most\nagricultural carbon programs impose strict requirements on “additionality”, so they require farmers to\nreporttheirfertilizeruseandwilldiscountthecarboncreditiffarmersincreasefertilizerrate.\nTo attract more farmers who have concerns about potential yield reduction when adopting\nconservation practices without increasing fertilizer use, carbon programs indeed should have better\nquantification approaches that can evaluate the comprehensive GHG outcomes, especially the tradeoff\nbetween SOC sequestration and N2O emissions. However, most of current quantification approaches\nentail costly directmeasurementslikesoilsamplingtodetermineSOCchange(Wendt&Hauser,2013)or\neddy-covariance flux sensors for GHG emissions measurement (Baldocchi et al., 1988), or they use\n10\nsimple emission factor estimates to approximate 'carbon outcomes' based on various management\npractices(IPCC,2019).\nIn light of these considerations, the role of KGML-ag-Carbon will become pivotal in the existing\ncarbon program due to its high accuracy,efficiency,andcost-effectiveness.TheKGML-ag-carbonmodel\ndoes not incentivize any particular management practices. Instead, it provides a novel means of\nquantifying the carbon outcomes of varying practices. In the US, this might mean theimplementationof\ncovercropsandreducedtillagewithoutanincreaseinfertilizer.\nCurrently"}]

Decision letter