PaperlayerView exact revision on arXiv

Source-linked AI summary

Advancing global crop simulations with differentiable modelling

Yunan Lin, Sebastian Bathiany, Maximilian Gelbrecht, Jens Heinke, Christoph Müller, Niklas Boers

arXiv:2610.36827v1cs.CEphysics.geo-ph

TL;DR

Systematic calibration of global gridded crop models remains difficult, despite its importance for reliable crop projections. This paper uses automatic differentiation in a fully process-based model to calibrate parameters against site-scale carbon and water flux observations. For rice under SSP5–8.5, optimization reduces the projected end-of-century production increase from 55.1% to 22.1%.

  • Problem

    Systematic calibration of global gridded crop models remains difficult, and uncertain parameterizations can diverge under future climates.

    Source: Introduction · p. 2

  • Method

    The authors implement differentiable Agrocosm and use automatic differentiation to optimize selected model parameters against site-scale carbon and water flux observations.

    Source: Introduction · p. 3

  • Results

    For rice under SSP5–8.5, optimization reduces the projected end-of-century production increase from 55.1% to 22.1%, while similarly skilled historical parameterizations can produce markedly different future responses.

    Source: Discussion · p. 12

  • Takeaways & Limitations

    Food-security assessments based on a single parameterization may underestimate uncertainty in crop projections under future climate change.

    Source: Discussion · p. 13

    Flux-tower observations are limited and unevenly distributed, global simulations use a single parameter set per crop, and the 20-member bootstrap does not cover the full parameter range.

    Source: Discussion · p. 13

Abstract

from arXiv · show

Global gridded crop models (GGCMs) are crucial for assessing the performance of food producing systems, for example under climate change or their contribution to environmental challenges. However, systematic calibration of their parameters remains a long-standing challenge. Persisting uncertainties in these parameters propagate into future crop projections, thereby affecting and potentially biasing food security assessments. Here we present a fully process-based differentiable GGCM implementation that offers a fundamentally new approach to parameter calibration, uncertainty quantification, yield loss attribution, and management optimization. Using automatic differentiation, we optimize selected model parameters against site-scale carbon and water flux observations, reducing normalized root-mean-square errors at held-out sites by 24.8-36.2% compared to the default state-of-the-art model. Without calibration against yield data, the optimized model better captures detrended yield anomalies across 61.2-64.3% of global crop area and reduces low-yield errors by 3.4-43.5% for wheat, rice, maize and soybean. It also ranks among the best performers in comparison with GGCMI Phase 3 models. Future simulations show that calibration reduces the projected end-of-century rice production increase from 55.1% to 22.1% under SSP5-8.5. Parameter ensembles further reveal that parameterizations with comparable historical skill can produce markedly different future production responses. Beyond calibration, differentiability enables the attribution of crop losses under extreme events and supports the design of adaptive crop management strategies. For the 2018 European drought, we identify spring heat and water limitation as major contributors to wheat yield losses and develop more effective irrigation strategies for mitigating simulated regional yield losses.

Introduction

Uncertain process representations and parameters weaken climate-related crop-risk projections, while systematic calibration of global crop models remains difficult. Agrocosm applies automatic differentiation to calibrate a fully process-based crop model and examine crop responses, uncertainty, and extreme-weather impacts.

Introduction

  • Uncertainty in process representation and parameterization makes global crop projections diverge, while models often underestimate low-yield risks and losses during climatic extremes.These limitations reduce confidence in assessments of production risks and adaptation benefits under climate change.Source: Introduction · p. 2
  • Systematic calibration remains challenging because current approaches rely on expert judgement and repeated forward searches, often targeting end-of-season yield rather than sub-seasonal processes.Growing-season carbon and water flux observations are increasingly available but typically do not constrain those dynamics.Source: Introduction · p. 2
  • Parameterizations with similar historical performance may diverge under future climates, while forward-simulation ensembles and growing-season-mean emulators constrain extreme-event sensitivity analysis.Forward ensembles are computationally intensive and restricted to predefined perturbations; the emulators are less reliable outside their training ranges and do not resolve within-season timing.Source: Introduction · p. 2
  • Agrocosm introduces gradient-based parameter optimization in a fully process-based GGCM, using automatic differentiation to compute sensitivities of crop growth, soil dynamics, and yield to parameters and daily weather.The authors describe this as the first demonstration of gradient-based parameter optimization in a fully process-based GGCM.Source: Introduction · p. 3
  • The study calibrates Agrocosm against site-scale carbon and water fluxes, then evaluates historical yield variability, future responses, parameter uncertainty, and 2018 European wheat losses and irrigation strategies.For rice under SSP5–8.5, parameter uncertainty accounts for about twice as much end-of-century projection variance as climate-model uncertainty.Source: Introduction · p. 3

Results

Flux-based parameter calibration improves held-out carbon and water flux simulations and historical yield evaluations without using yield data for optimization. It also changes future crop projections, reveals parameter-driven uncertainty, and supports weather-loss attribution and irrigation design.

Calibration using automatic differentiation

  • Held-out-site nRMSE decreases by 34.9% for GPP, 24.8% for RECO, and 36.2% for ET after flux-based calibration.Across 27 sites, mean temporal correlation rises from 0.65 to 0.81 for GPP, 0.65 to 0.71 for RECO, and 0.54 to 0.74 for ET.Source: Figure 1 · p. 4, Calibration using automatic… · p. 5

Historical yield evaluation

  • 64.0%, 64.3%, 62.3%, and 61.2% of evaluated wheat, rice, maize, and soybean area, respectively, show improved yield-anomaly correlation without yield-based optimization.Grid-cell RMSE also decreases across 96.6%, 89.1%, 67.2%, and 85.1% of those crop areas, respectively.Source: Historical yield evaluation · p. 5, Figure 2 · p. 6
  • 43.5%, 21.0%, 3.4%, and 15.4% are the reductions in low-yield-year crop area-weighted nMAE for wheat, rice, maize, and soybean, respectively.With optimized parameters, Agrocosm ranks first by mean rank for temporal correlation and nRMSE among the compared models, improving from third with default parameters.Source: Historical yield evaluation · pp. 5, 7, Figure 2 · p. 6

Future projections across parameterizations

  • Under SSP5–8.5, end-of-century production changes shift from 12.8% to 17.9% for wheat, 55.1% to 22.1% for rice, and −15.5% to −12.5% for maize; soybean remains nearly unchanged at 37.4%.The same directions occur under SSP1–2.6, but with smaller changes.Source: Future projections across… · p. 7, Figure 3 · p. 8
  • Despite higher historical yield correlation than the default simulations, bootstrap parameterizations span end-of-century production changes of 21.83–49.35% for rice under SSP5–8.5.Parameter sampling accounts for 63% of crossed-ensemble variance in rice production responses.Source: Figure 4 · p. 9, Future projections across… · p. 10
  • A yield-change IQR of at least 10% occurs across 74.3% of rice area, compared with 3.2% for wheat, 27.4% for maize, and 13.2% for soybean.At least 17 of 20 parameter sets agree on response direction across 92.0–98.9% of crop area, while the magnitude of spread varies by crop and region.Source: Figure 4 · p. 9, Future projections across… · p. 10

Extreme event attribution

  • In the worst-affected quarter of European wheat area in 2018, May sensitivities indicate that increased rainfall raises yield while increased temperature and radiation lower it.Replacing weather in May with reference-year conditions produces the largest median yield gain, +0.16 t ha−1, followed by June at +0.11 t ha−1.Source: Extreme event attribution · pp. 10, 12, Figure 5 · p. 11
  • Across seasonal irrigation budgets of 10–100 mm, AD-guided irrigation recovers 24.4–113.4% of the simulated yield deficit, versus 22.6–92.0% for deficit-based and 12.7–85.2% for uniform irrigation.Its recovery advantage over deficit-based irrigation ranges from 3.4 to 21.4 percentage points.Source: Figure 5 · p. 11, Extreme event attribution · p. 12

Discussion

Systematic gradient-based calibration improves a process-based crop model beyond its flux-observation targets, while altering future projections and enabling weather attribution and irrigation analysis. The study also identifies limits from sparse observations and persistent difficulty simulating low-yield failures.

Model performance

  • Calibration moves Agrocosm from third to first among 13 AgMIP models for mean temporal-correlation and nRMSE ranks, though no model performs consistently well across crops.Optimized simulations also have lower low-yield nMAE than the AgMIP ensemble median for all four crops.Source: Discussion · p. 12

Future projections and uncertainty

  • Under SSP5–8.5, optimization reduces the projected end-of-century rice production increase from 55.1% to 22.1%, while similarly skilled parameterizations can yield markedly different future responses.For rice, parameter uncertainty accounts for 63.3% of projection variance under SSP5–8.5, exceeding the 31.6% attributed to driving climate models.Source: Discussion · pp. 12–13

Yield attribution and management

  • For the 2018 European wheat event, May temperature and June precipitation are leading contributors to simulated yield loss, and AD-guided irrigation alleviates simulated losses more effectively than deficit-based or evenly distributed irrigation.The diagnosis is consistent with independent evidence linking the 2018 European wheat losses to heat and drought.Source: Discussion · p. 13

Limitations

  • Unevenly distributed flux-tower observations limit robust spatially varying parameter estimation, so global simulations use one parameter set per crop despite possible regional variation in physiological thresholds.The 20-member bootstrap does not cover the full range of parameters, and attribution is limited to events the model can reproduce and represent.Source: Discussion · p. 13
  • Low-yield improvements remain modest for rice, maize and soybean, and accurately representing yield failure from heat, drought and excess water remains challenging.The authors identify improved stress representations, harmonized observations and independent validation as requirements for addressing this challenge.Source: Discussion · p. 13

Methods

The methods implement Agrocosm as a differentiable crop–soil model, initialize it with a target-constrained agricultural spin-up, and use flux observations to calibrate parameters. Historical and future simulations, parameter ensembles and event-level analyses evaluate model behavior and explore uncertainty and management strategies.

Model description

  • Agrocosm implements LPJmL crop and soil formulations in Julia with daily coupled dynamics and a CPU/GPU-compatible architecture supporting automatic differentiation.The global simulations use a 0.5° × 0.5° grid; the implementation supports gradient-based optimization without changing standard forward-simulation routines.Source: Model description · pp. 13–14

Model spin-up

  • Agricultural soil states are initialized from HWSD soil targets using a 600-year spin-up under fixed 2015 management and repeating 1901–1930 climate.SOC and total nitrogen are constrained during spin-up, while litter pools evolve freely; initial soil water is set to field capacity.Source: Model spin-up · p. 14

Gradient-based model optimization

  • Reverse-mode AD computes gradients of observation loss through selected continuous daily transitions, while initialization and discrete crop events remain outside differentiation or fixed during gradient evaluation.Calibration uses 100 growing-season records from 27 flux sites, optimizes crop parameters against GPP, RECO and ET, then soil parameters against RECO and ET with crop parameters fixed.Source: Gradient-based model optimization · pp. 14–15

Historical crop simulations

  • Historical simulations compare default and optimized parameters for four crops against detrended GDHY yields, evaluate low-yield years, and benchmark against 13 AgMIP models.Global simulations run at 0.5° resolution for 1901–2019; low-yield years are each grid cell’s lowest 20% of detrended GDHY anomalies.Source: Historical crop simulations · pp. 15–16

Future crop simulations

  • Future simulations compare default and optimized parameters under SSP1–2.6 and SSP5–8.5 using five climate models, with management held at 2015 levels.Production responses are measured relative to 1995–2014, and AD estimates yield sensitivities to CO2 and temperature across specified future periods.Source: Future crop simulations · pp. 16–17

Parameter uncertainty

  • A 20-replicate site bootstrap recalibrates parameters by resampling 27 physical sites, then evaluates historical skill and future production responses across five climate models.Variance decomposition separates contributions from parameter sampling, climate model and their interaction; grid-cell response robustness requires agreement from at least 17 of 20 parameter sets.Source: Parameter uncertainty · pp. 17–18

Yield loss attribution

  • For the 2018 European wheat event, AD measures yield sensitivity to daily weather, counterfactual simulations replace weather inputs, and irrigation schemes are compared under shared seasonal water budgets.The schemes are AD-guided, deficit-based and evenly distributed supplementary irrigation, evaluated using actual water use and yield recovery.Source: Yield loss attribution · pp. 18–19

Supplementary information

Agrocosm’s soil states are initialized from HWSD stocks and static soil properties, then refined through a target-constrained agricultural spin-up. The supplementary methods describe the stock remapping, initial pool allocations, and convergence checks used before production simulations.

A.1.1 HWSD soil inputs

  • HWSD stocks are remapped from seven source layers to five Agrocosm layers by conserving areal stocks under uniform within-layer stock density.Component shares are normalized, structural zeros are distinguished from unresolved values, and layers are retained only when at least 99% of component share is resolved.Source: §A.1.1 · p. 27
  • For the 2–3 m Agrocosm layer, stock density from HWSD’s deepest layer is extrapolated, so these carbon and nitrogen stocks are estimates rather than measurements.Source: §A.1.1 · p. 28

A.1.2 Construction of initial soil states

  • Initial SOC and organic nitrogen use a 40:60 fast-to-slow split; mineral nitrogen, litter, and vertical-routing allocations receive provisional starting values.Initial water is set to field capacity estimated from soil texture and SOC-derived organic matter; spin-up later replaces the provisional pool allocations.Source: §A.1.2 · p. 28

A.1.3 Spin-up and convergence criteria

  • A 600-year agricultural spin-up repeats 1901–1930 climate under fixed 2015 management, annually restoring soil pools to HWSD targets while litter pools evolve freely.Spin-up-derived pool fractions and vertical-routing fractions initialize production simulations, where annual target correction is not applied.Source: §A.1.3 · pp. 28–29

A.2 Spin-up of four crops

  • Convergence checks compare states at the same phase of successive 30-year cycles; Table S1 reports per-cell convergence for four crops under rainfed and irrigated conditions.Stability requires specified stock and correction changes no greater than 0.01 and fast-pool fraction changes no greater than 0.01 for three consecutive comparisons.Source: §A.1.3 · p. 29, Table S1 · p. 29

B FLUXNET observations

The FLUXNET calibration constraints are daily gross primary production, ecosystem respiration, and evapotranspiration, paired with meteorological forcing and crop-season information. The retained dataset covers 100 growing-season records at 27 physical sites across four crops.

B FLUXNET observations

  • Daily GPP, RECO, and ET are the observational constraints, with ET derived from gap-filled latent heat flux and air temperature.GPP and RECO use the night-time flux-partitioning products.Source: B FLUXNET observations · p. 30
  • Candidate seasons require complete daily coverage of five Agrocosm forcing variables: air temperature, precipitation, longwave and shortwave radiation, and atmospheric CO2.Short CO2 gaps are filled using the same no-leap calendar day from another available year.Source: B FLUXNET observations · p. 30
  • The final calibration set contains 100 growing-season records from 27 physical sites: 36 wheat, 22 maize, 26 rice, and 16 soybean records.Table S2 lists site identifiers, coordinates, and crop coverage; Fig. S1 maps the global and regional site distribution.Source: B FLUXNET observations · p. 30, Table S2 · p. 30, Figure S1 · p. 31

Station-flux metrics

Station-flux fit is summarized using Pearson correlation and normalized root-mean-square error. Both metrics are calculated at the growing-season and flux level.

Station-flux metrics

Global yield metrics

Global yield metrics compare simulated and observed detrended yield variability, including performance in low-yield years. Parameter selection uses sensitivity and optimization diagnostics to retain crop and soil parameters for calibration.

Global yield metrics

  • Combined grid-cell yield is the crop-area-weighted mean of rainfed and irrigated yields, and yield variability is evaluated with detrended-anomaly correlation and RMSE.The normalized RMSE uses the standard deviation of detrended GDHY yield anomalies as its scale.Source: Global yield metrics · p. 32
  • Low-yield years are the lowest 20% of detrended GDHY yield anomalies, and their error is summarized with normalized mean absolute error.Source: Global yield metrics · p. 33

D.1 Parameter selection

  • The candidate parameter set spans canopy, photosynthesis, respiration, nitrogen allocation, soil-carbon turnover, evaporation, and transpiration processes; discrete management decisions remain fixed during calibration.Source: §D.1 · p. 34
  • Parameter selection evaluates the normalized trajectory Jacobian, parameter correlations, multi-start optimization, and leave-one-parameter-out experiments, retaining parameters with reproducible gradients and independent information.Source: §D.1 · p. 34
  • The final calibration fits 12 crop parameters and four soil parameters per crop, with parameter vectors fitted independently for wheat, maize, rice, and soybean.Source: §D.1 · p. 34

D.2 Two-stage calibration

The calibration separates crop- and soil-parameter fitting because simultaneous optimization produced less stable, poorer fits. It balances observations across sites, seasons, and crops while keeping GPP as an independent soil-stage diagnostic.

D.2 Two-stage calibration

  • Two-stage optimization reduces the search dimension and limits compensation between crop and soil parameters after simultaneous fitting produced less stable, poorer results.Crop parameters are fitted independently for each crop, then held fixed while shared soil parameters are optimized across crops.Source: §D.2 · pp. 34–36
  • Calibration bounds define optimization search ranges, not universal limits, for crop and soil parameters.The crop- and soil-parameter bounds are listed in Tables S3 and S4.Source: Table S3 · p. 35, Table S4 · p. 36
  • The crop-stage objective gives each physical site equal total weight within a crop, preventing sites with more years from dominating.The weighting is defined over site–seasons for each crop.Source: §D.2 · p. 35
  • The soil-stage objective weights crops, physical sites within crops, and seasons within sites equally, while excluding GPP as an independent diagnostic.The soil stage optimizes four shared parameters with crop parameters fixed.Source: §D.2 · p. 36
  • Full-data calibration uses a learning rate of 0.01, 100 crop-stage epochs, 20 soil-stage epochs, and 30-day reverse-mode checkpointing blocks.Adam uses β1 = 0.9, β2 = 0.999, and ϵ = 10−8.Source: §D.2 · pp. 36–37

D.3 Grouped-site cross-validation

Grouped-site fivefold cross-validation evaluates parameter fits on physical sites withheld from calibration. The fits converge across folds, and held-out results show broad gains in GPP and ET, with less uniform RECO responses.

D.3 Grouped-site cross-validation

  • All seasons and crops from each physical site stay in one fold, preventing site-specific soil and climate information from crossing between training and validation data.Fold-specific fits are refitted from default values and evaluated on held-out sites without parameter updates.Source: §D.3 · p. 37
  • Across crops, mean training-loss reductions range from 38.4% to 62.2% for crop parameters and from 12.3% to 28.3% for soil parameters.Final reported losses remain within 0.05% of the best evaluated epoch on average.Source: §D.3 · p. 37, Figure S2 · p. 38
  • Strong GPP and ET gains occur across all four crops, while RECO responses are less uniform; wheat improves consistently across all fluxes.Maize RECO has the smallest nRMSE reduction (7.3%), and rice RECO correlation increases by ∆r = 0.067.Source: §D.3 · p. 38, Figure S3 · p. 39
  • Performance is summarized as an unweighted mean of 27 site-level values, so each physical site contributes equally regardless of its crop or season count.Table S6 reports the grouped-site cross-validation performance, and Table S7 separates crop- and soil-stage contributions.Source: §D.3 · p. 38, Table S6 · p. 38, Table S7 · p. 39

D.4 Full-data parameter optimization

After cross-validation, separate parameter sets are retrained on all 100 growing seasons for use in global simulations. The final crop and soil parameterizations are reported separately.

D.4 Full-data parameter optimization

  • Separate parameter sets for winter wheat, spring wheat, maize, rice, and soybean are retrained using all 100 growing seasons for global simulations.These full-data sets are retrained after cross-validation.Source: §D.4 · p. 40
  • The final full-data crop parameterization used in global simulations is reported in Table S8, with values shown to three significant digits.Simulations use full-precision values.Source: Table S8 · p. 40
  • The final full-data soil parameterization used in global simulations is reported in Table S9, with long values shown to three significant digits.Simulations use full-precision values.Source: Table S9 · p. 41

E Global historical yield evaluation

The section evaluates historical yield performance and future production projections across crop parameterizations. It also frames how parameter ensembles with similar historical skill can yield different future responses.

E Global historical yield evaluation

  • Observation-based calibration improves historical low-yield simulations, evaluated by nMAE for the lowest 20% of detrended yield anomalies during 1981–2016.The figure compares wheat, rice, maize, and soybean; lower nMAE indicates better performance.Source: Figure S4 · p. 43, Figure S5 · p. 44
  • Historical yield evaluation includes grid-cell temporal correlations, normalized RMSE, and changes in seasonal process–yield correlations for 1981–2016.The process variables include GPP, LAI-days, ET, water deficit, growing-season length, and aboveground biomass.Source: Figure S6 · p. 46, Figure S7 · p. 47, Figure S8 · p. 48

F Parameter-ensemble diagnostics

  • The ensemble tables frame future production differences across parameterizations and climate models, with a two-factor ANOVA partitioning variation into parameter, climate-model, and interaction components.Production responses cover 2080–2099; the supplied passages do not state the resulting numerical comparisons.Source: Table S11 · p. 48, Table S12 · p. 51, Table S13 · p. 52
  • For SSP1–2.6, parameterizations with similar historical performance are presented as producing different future projections.The diagnostics compare historical correlation, end-of-century global production changes, variance shares, trajectories, and grid-cell yield-change spread.Source: Figure S11 · p. 53

Related papers

Loading 2610.36827v1…