From Weather to Enhanced Weathering: Weather Foundation Models for Carbon Dioxide Removal at Operational Scale

# From Weather to Enhanced Weathering: Weather Foundation Models for Carbon Dioxide Removal at Operational Scale
Source: [https://arxiv.org/html/2610.28973](https://arxiv.org/html/2610.28973)
Sasankh Munukutla††thanks:Corresponding author\.Affiliation:Terradot Soil, Inc\.Affiliation:San Francisco, CA, USAEmail:[sasankh@terradot\.earth](mailto:)Danrong ZhangAffiliation:Microsoft ResearchAffiliation:Redmond, WA, USAKenji TakedaAffiliation:Microsoft ResearchAffiliation:Cambridge, UKJuan M\. Lavista FerresAffiliation:Microsoft AI for GoodAffiliation:Redmond, WA, USA

###### Abstract

Enhanced Rock Weathering \(ERW\) is a carbon removal solution with the potential to remove gigatons of CO2per year, but quantifying that removal today requires dense soil sampling and laboratory analysis\. We evaluate five weather Foundation Models \(FMs\) as frozen encoders feeding lightweight decoders, under a leave\-one\-field\-out protocol on 5,192 chronosequences from 264 fields across Brazil and the US\. To our knowledge, this is the first holistic evaluation of weather FMs for ERW and the largest evaluation on operational ERW data\. We find that across weather FMs, representations transfer to ERW: our best model, TerraNova\-ERW, built on Aurora 1\.5, positively orders measurements within fields it has never seen\. Zero\-shot prediction of a new field’s mean removal is not yet feasible, and we attribute this to data scale: the same architecture excels on∼\\sim10510^\{5\}natural\-weathering sequences, anchoring the scale we expect ERW data must reach\. Few\-shot prediction, blending the model with a few samples per field, is supported today: TerraNova\-ERW recovers over 60% of the highest\-removal fields with two samples per field, roughly halves the samples needed to estimate a field’s mean removal, and can cut sampling on our densest fields by 28–98%\.

## 1Introduction

To avoid the most adverse effects of climate change, nearly every mitigation pathway requires Carbon Dioxide Removal \(CDR\) to balance hard\-to\-abate residual emissions, with∼\\sim10 gigatons of CDR required annually by 2050\[[1](https://arxiv.org/html/2610.28973#bib.bibx1)\]\.

Enhanced Rock Weathering \(ERW\) is a promising CDR solution, which accelerates Earth’s natural weathering process that draws down carbon over geological timescales\. By applying finely crushed silicate rocks, such as basalt, to agricultural soils, ERW accelerates this process to occur in years rather than millennia\. Scaled globally, ERW has the potential to remove gigatons of CO2per year while providing agronomic co\-benefits such as improved crop yields and soil health\[[2](https://arxiv.org/html/2610.28973#bib.bibx2)\]\.

The key obstacle to scaling up ERW today is the Measurement, Reporting, and Verification \(MRV\) required to quantify the CDR\. Quantifying removal currently requires dense physical ground\-truth samples, including soil samples, along with rigorous laboratory analysis at every site\. Process\-based models such as Reactive Transport Models \(RTMs\) have been used to model ERW, but they require explicitly defining and approximating all relevant processes, such as mineral kinetics, and can be challenging to parametrize and computationally expensive to scale\. Meanwhile, deep learning is increasingly applied across climate and Earth system problems\[[3](https://arxiv.org/html/2610.28973#bib.bibx3)\], and weather Foundation Models \(FMs\) pretrained on decades of reanalysis data offer rich representations of the atmospheric drivers, such as precipitation, temperature, and soil moisture, that govern weathering rates\.

Previous work by\[[4](https://arxiv.org/html/2610.28973#bib.bibx4)\], to our knowledge, the first application of deep learning to ERW, trained lightweight decoders on frozen representations from Microsoft’s Aurora\[[5](https://arxiv.org/html/2610.28973#bib.bibx5)\]to show that atmospheric FM representations transfer to rock weathering\. We build upon and advance this work by evaluating five weather FM backbones on our own operational ERW data spanning Brazil and the US, to our knowledge, the first holistic evaluation of multiple weather FMs for ERW and the largest evaluation on operational ERW data\. We then show how our best model, TerraNova\-ERW, can be used operationally to reduce the sampling MRV requires\.

## 2Methods

### 2\.1Task and Dataset

We formulate ERW modeling as a sequence regression task using atmospheric FM representations: given monthly historic atmospheric representations and site\-specific features \(static soil, dynamic site, amendment, spatial, temporal; see Appendix[A2](https://arxiv.org/html/2610.28973#A1.SS2)for the full list of features\) attit\_\{i\}, predict atti\+1t\_\{i\+1\}the CDR rate in t CO2e/ha/yr\.

Similar to\[[4](https://arxiv.org/html/2610.28973#bib.bibx4)\], we assemble a dataset of natural weathering data using the same sources and process\. We augment this with ERW data from our own deployments in the US and Brazil\. The data is assembled as chronosequences, such that a given training example is a pair of measurements from two time points\. In total, we have192,511192\{,\}511sequences of natural weathering data with global coverage and5,1925\{,\}192sequences of our own operational ERW data spanning 264 agricultural fields\. The two datasets play distinct roles: all main ERW results \(Tables[1](https://arxiv.org/html/2610.28973#S3.T1)and[5](https://arxiv.org/html/2610.28973#A1.T5)\) are on the ERW sequences\. The natural\-weathering data is used only in two separate experiments: a joint\-training experiment \(Appendix[A6](https://arxiv.org/html/2610.28973#A1.SS6)\) and a natural\-only scale study \(Appendix Table[7](https://arxiv.org/html/2610.28973#A1.T7)\)\.

### 2\.2Model Architecture

The general architecture of TerraNova\-ERW is summarized in Appendix[A1](https://arxiv.org/html/2610.28973#A1.SS1); key elements follow\.

#### Weather FM Backbones:

We consider five leading weather FMs to generate latents: Aurora 1\.5\[[6](https://arxiv.org/html/2610.28973#bib.bibx6)\], Pangu\[[7](https://arxiv.org/html/2610.28973#bib.bibx7)\], SFNO\[[8](https://arxiv.org/html/2610.28973#bib.bibx8)\], FCN3\[[9](https://arxiv.org/html/2610.28973#bib.bibx9)\], and GraphCast\[[10](https://arxiv.org/html/2610.28973#bib.bibx10)\]\. For each backbone, we precompute global monthly latent representations by running the pretrained, frozen model on primarily ERA5\[[11](https://arxiv.org/html/2610.28973#bib.bibx11)\], extracting the activations from the final encoder layer\. At each site, we pool the latent field at three fixed spatial scales \(local, regional, global\) and concatenate them into a singledmd\_\{m\}\-dimensional embedding per month, obtained for every month spanningtit\_\{i\}toti\+1t\_\{i\+1\}of a chronosequence; this keeps training lightweight\. Full details on the pooling scales and per\-backbone extraction are in Appendix[A3](https://arxiv.org/html/2610.28973#A1.SS3)\.

#### Temporal Aggregation:

The interval spanning a chronosequence can vary, and theTTmonthly embeddings spanningtit\_\{i\}toti\+1t\_\{i\+1\}must be reduced to a fixed\-size representation\. We aggregate by trapezoidal integration of the latent trajectory, selected over supervised and self\-supervised LTAE\[[12](https://arxiv.org/html/2610.28973#bib.bibx12)\]alternatives \(Appendix[A4](https://arxiv.org/html/2610.28973#A1.SS4)\)\.

#### Lightweight Decoders:

We build upon\[[13](https://arxiv.org/html/2610.28973#bib.bibx13)\]and\[[4](https://arxiv.org/html/2610.28973#bib.bibx4)\]to train lightweight decoders that fuse the aggregated weather representation with site features to predict the target\. Considering the backbone latents are high\-dimensional while labeled ERW data is scarce, decoder capacity is an important choice: we sweep the fusion mechanism and the per\-backbone latent\-projection width, favoring shallow decoders with wide hidden layers \(“light wide”\); the backbone latent is separately down\-projected before fusion\. Throughout,*TerraNova\-ERW*denotes the Aurora 1\.5 backbone paired with a light\-wide decoder and ad=64d\{=\}64latent projection \(“64d”\), with two refinements to its inputs and training: the three location inputs enter without per\-fold standardization, and early stopping monitors validation RMSE rather than MAE\. This is the configuration selected in Section[3](https://arxiv.org/html/2610.28973#S3)\.

### 2\.3Training and Evaluation

#### Leave\-One\-Field\-Out \(LOFO\):

We train and evaluate all models by rotating the held\-out field\. Fields sharing a 2\-decimal location key are merged into*components*\(264 fields→\\to101\), since splitting fields that are spatially close would leak\. Each fold trains with field\-disjoint ERW validation, and predicts the held\-out component \(full LOFO protocol details in Appendix[A5](https://arxiv.org/html/2610.28973#A1.SS5)\)\. We favor this evaluation approach over a traditional held\-out validation split given the dataset size\.

#### Point vs\. Site protocol:

*Point Level*scores every held\-out measurement, testing whether the model captures within\-field structure: the capability*model\-and\-measure \(few\-shot\)*MRV would use, where model predictions can help complement physical samples\.*Site Level*compares the predicted field mean to the true field mean, testing whether the model can place an unseen field’s overall removal without any on\-field measurement: the capability*fully model\-based \(zero\-shot\)*MRV would require\.

#### Metrics and Uncertainty:

We report MAE, RMSE, and Spearman rank correlation throughout; the row\-shuffled ERA5 experiment serves as the pooled no\-information reference\. Uncertainty \(±\\pm\) is a bootstrap over held\-out*fields*: resampling the scored components \(64 for treated\-row metrics\) with replacement 1000 times\. At the point level, we also report “% best”, the share of the closable error a model removes:

%best=100×constant−modelconstant−best\\%\\,\\text\{best\}\\;=\\;100\\times\\frac\{\\text\{constant\}\-\\text\{model\}\}\{\\text\{constant\}\-\\text\{best\}\}where the constant is each fold’s treated training mean and the best is each field’s own measured median \(MAE\) or mean \(RMSE\), the best any single value per field can achieve\. It is undefined at site level, where the best case is the exact field mean\.

#### Experimental Configurations:

All experiments use trapezoidal aggregation and vary a single design choice, so each contrast isolates one effect: one configuration per weather FM backbone \(which backbone transfers best\), raw ERA5 inputs in place of the latent \(the contribution of the FM encoding\), and row\-shuffled ERA5 \(isolating the value of weather information itself\), plus natural\-weathering co\-training, input\-group ablation, and tabular AutoML configurations detailed in Appendix[A6](https://arxiv.org/html/2610.28973#A1.SS6)\.

## 3Results

Table 1:LOFO results on treated rows\. Point level scores every held\-out measurement; site level compares predicted vs\. true field means\.Boldmarks the best value in each block; “% best” is defined in Section[2\.3](https://arxiv.org/html/2610.28973#S2.SS3)and reported at point level only\. Per\-country \(US, Brazil\) point and site breakdowns, backbone comparisons, and other configurations are in Appendix Table[5](https://arxiv.org/html/2610.28973#A1.T5)\.#### FM representations transfer to CDR prediction across architectures; Aurora performs best\.

We see that the weather FM representations transfer to CDR across all backbones, while shuffling the weather away destroys ranking, showing weather*information*is necessary \(all\-ERW point MAE 3\.06–3\.26 and Spearman\+0\.12\+0\.12to\+0\.29\+0\.29across backbones and Aurora variants, vs\. 3\.32 and−0\.01\-0\.01shuffled; raw ERA5 3\.21 and\+0\.22\+0\.22; Appendix Table[5](https://arxiv.org/html/2610.28973#A1.T5)\)\. The FM encoding matches raw ERA5 on error size at the base configuration, while TerraNova\-ERW more than doubles the share of closable error removed \(41% of MAE and 54% of RMSE vs\. 19% and 18% for raw ERA5; Table[1](https://arxiv.org/html/2610.28973#S3.T1)\), and is more robust\. Raw ERA5 ranking collapses on one of three seeds \(per\-seed Spearman−0\.02\-0\.02to\+0\.22\+0\.22\), whereas the Aurora \(base, 64d, 512d\), SFNO, FCN3, and GraphCast encodings hold\+0\.16\+0\.16to\+0\.28\+0\.28on every seed; only Pangu and the natural\-weathering co\-training variant share ERA5’s instability \(Appendix Table[6](https://arxiv.org/html/2610.28973#A1.T6)\)\. Aurora 1\.5 with a light\-wide decoder, ad=64d\{=\}64latent projection, unstandardized location inputs, and RMSE early stopping achieved the best value on every all\-ERW metric \(point MAE 3\.06, RMSE 4\.20, Spearman\+0\.29\+0\.29; site MAE 1\.78, RMSE 2\.57, Spearman\+0\.26\+0\.26\) while being stable on every seed, so we adopt it as our best model, which we call TerraNova\-ERW\. It was selected among about a dozen decoder and training variants evaluated on the same 68 folds\.

#### Model\-and\-Measure MRV is supported; fully model\-based MRV is not yet\.

At the point level on treated rows, weather FMs positively order held\-out measurements on fields they have never seen \(Appendix Table[5](https://arxiv.org/html/2610.28973#A1.T5)\) with TerraNova\-ERW’s pooled Spearman\+0\.29\+0\.29\(Brazil\+0\.40\+0\.40\), which is the capability model\-and\-measure MRV builds on\. At the site level the model also orders unseen fields \(pooled Spearman\+0\.26\+0\.26; US\+0\.17\+0\.17, Brazil\+0\.37\+0\.37\), but its error on an unseen field’s mean is only marginally below the no\-information reference: pooled site MAE is 1\.78±\\pm0\.18 vs\. 1\.96 for Shuffled ERA5, and on the US subset 1\.50±\\pm0\.16 vs\. 1\.57 \(Appendix Table[5](https://arxiv.org/html/2610.28973#A1.T5)\)\. Single fields with few measurements are high\-variance targets, so ranking is recoverable where absolute measurement is not\. Fully model\-based quantification of unseen fields is therefore not yet supportable\.

#### Natural Weathering Scale:

Trained and evaluated exclusively on natural weathering chronosequences, the same TerraNova\-ERW architecture reaches validation MAE 0\.014 and RMSE 0\.017, roughly 90% lower error than the deployable constant predictor \(Appendix Table[7](https://arxiv.org/html/2610.28973#A1.T7)\)\. Given∼\\sim105diverse sequences, this model family excels at natural weathering prediction\. This anchors our expectations for the needed data scale for fully model\-based MRV for ERW\. From our own scale\-up and the growth of the ERW industry, we estimate this scale will be reached in a few years\.

## 4Operational Deployment of TerraNova\-ERW

Figure 1:Sampling\-budget simulation \(Appendix[A8](https://arxiv.org/html/2610.28973#A1.SS8)\)\. \(a\) Share of the ten truly highest\-removal fields recovered when ranking 52 fields under a budget ofkkmeasured points per field\. \(b\) Median error of the estimated field mean vs\. budget\. \(c\) Percent of samples the hybrid saves on the three densest fields at two precision targets\.#### Site Selection with Limited Features:

Site\-specific features require physical samples taken at the field or, as a fallback, gridded products\. We test whether a reduced\-input variant of the model can operate before either is available, using only the weather FM latents, obtainable globally from a field’s coordinates alone, together with the planned amendment specification\. This Aurora\+\+Amendment configuration matches the full\-input model at the same decoder settings on point\-level error and orders unseen fields better \(site Spearman\+0\.32\+0\.32vs\.\+0\.13\+0\.13for Aurora 1\.5 \(base\); Appendix[A7](https://arxiv.org/html/2610.28973#A1.SS7)\), and it remains competitive with TerraNova\-ERW on site ranking \(\+0\.26\+0\.26\) without any site characterization\. The advantage lies in ranking, not error: this supports prioritizing sites rather than quantifying them, and lets sampling and site visits be deferred\.

#### Reducing Samples for MRV:

Since TerraNova\-ERW supports model\-and\-measure MRV, we quantify what the model is worth in samples by simulating sparser campaigns on the LOFO out\-of\-fold predictions: each Monte\-Carlo draw revealskkrandom measured intervals per field, and a*hybrid*estimator blends the revealed points with the model prediction at weights fitted out\-of\-sample \(Protocol and results in Appendix[A8](https://arxiv.org/html/2610.28973#A1.SS8)\)\. Three results stand out \(Figure[1](https://arxiv.org/html/2610.28973#S4.F1)\)\.*\(1\)*When ranking candidate fields, the hybrid recovers on average 63% of the ten highest\-removal fields with two samples per field, vs\. 49% for sampling alone at the same budget; the model alone, with no samples, already finds 50%, as many as two samples per field \(Figure[1](https://arxiv.org/html/2610.28973#S4.F1)a\)\.*\(2\)*At every budget below∼10\{\\sim\}10points per field, the hybrid estimates field means with lower median error than sampling alone \(1\.13 vs\. 2\.18 t CO2e/ha/yr at one point\), and the hybrid at three points matches what sampling alone needs five to ten points to reach \(Figure[1](https://arxiv.org/html/2610.28973#S4.F1)b\)\.*\(3\)*On our three densest fields, which hold 32% of all treated measurements, the hybrid saves 28–98% of samples \(Figure[1](https://arxiv.org/html/2610.28973#S4.F1)c\): at the±1\.0\\pm 1\.0t target one sample replaces 6–11 samples on all three fields, and at±0\.5\\pm 0\.5t Fields B and C need a single sample \(96% and 98% saved\) because the model alone already places their means within 0\.10 and 0\.32 t, while Field A still needs 44 of 62 \(28% saved\)\.

## 5Discussion

Pathway to climate impact:MRV accounts for a third to over half of the cost of ERW\[[14](https://arxiv.org/html/2610.28973#bib.bibx14)\]\. Since sampling and laboratory analysis dominate MRV spend, cutting samples by 28–98% attacks the largest term in a cost that can reach $275/tCO2\[[15](https://arxiv.org/html/2610.28973#bib.bibx15)\]\. TerraNova\-ERW complements rather than replaces physical sampling, enabling sparser, targeted sampling designs and helping prioritize site selection\. As latents can be precomputed globally and the decoder is lightweight, the marginal cost of modeling a new site is near zero, supporting cost\-effective MRV as ERW scales to gigaton removal\.

Limitations and future work:Our deployment data covers only two countries and two feedstock types; temporal horizons are limited to monitoring to date; LOFO does not cover entirely new regions\. Further validation across lithologies, climates, regions, and longer horizons along with training on more ERW data will enable models like TerraNova\-ERW to be more robustly deployed\.

## Acknowledgments and Disclosure of Funding

We gratefully acknowledge Katya Larina, Jackson Harris, Naomi Asimov Alley, Laisa de Assis Batista, Jono Finger, Jon Golla, Ella Holme, William Larsen, Jennifer Mills, Elaine Pegoraro, Silvia Palotti Polizel, Soroush Rezvan, and Adam Wolf for their invaluable contributions in assembling and curating the data\. We thank Prof\. David Lobell and Prof\. Scott Fendorf for their expert review and domain feedback, which significantly strengthened this work\. We thank the anonymous reviewers for their constructive feedback\. This work was supported by Terradot, which provided the operational ERW data\. Finally, we thank the Microsoft AI for Good and Aurora teams for providing compute resources and for their generous support and responsiveness as we built on Aurora\.

## References

- \[1\]IPCC“Climate Change 2022: Mitigation of Climate Change”In*Contribution of Working Group III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change*Cambridge University Press, 2022
- \[2\]David Beerling et al\.“Potential for large\-scale CO2removal via enhanced rock weathering with croplands”In*Nature*583\.7815Nature Publishing Group, 2020, pp\. 242–248
- \[3\]Grey Nearing et al\.“Global prediction of extreme floods in ungauged watersheds”In*Nature*627\.8004Nature Publishing Group, 2024, pp\. 559–563DOI:[10\.1038/s41586\-024\-07145\-1](https://dx.doi.org/10.1038/s41586-024-07145-1)
- \[4\]Sasankh Munukutla, Felipe Oviedo and Kenji Takeda“From Weather to Weathering: Foundation Models for Carbon Removal via Rock Weathering”In*ICLR 2026 Workshop on Foundation Models for Science: Real\-World Impact and Science\-First Design*, 2026
- \[5\]Cristian Bodnar et al\.“A foundation model for the Earth system”In*Nature*641, 2025, pp\. 1180–1187DOI:[10\.1038/s41586\-025\-09005\-y](https://dx.doi.org/10.1038/s41586-025-09005-y)
- \[6\]Jonathan Weyn et al\.“Aurora 1\.5: Fine\-Tuning a Foundation Model for Medium\-Range Ensemble Weather Prediction”, 2026URL:[https://www\.microsoft\.com/en\-us/research/wp\-content/uploads/2026/07/Aurora\_1\_5\_Paper\.pdf](https://www.microsoft.com/en-us/research/wp-content/uploads/2026/07/Aurora_1_5_Paper.pdf)
- \[7\]Kaifeng Bi et al\.“Pangu\-Weather: A 3D high\-resolution model for fast and accurate global weather forecast”In*Nature*619\.7970Nature Publishing Group, 2023, pp\. 533–538
- \[8\]Boris Bonev et al\.“Spherical fourier neural operators: Learning stable dynamics on the sphere”In*International conference on machine learning*, 2023, pp\. 2806–2823PMLR
- \[9\]Boris Bonev et al\.“FourCastNet 3: A geometric approach to probabilistic machine\-learning weather forecasting at scale”In*arXiv preprint arXiv:2507\.12144*, 2025
- \[10\]Remi Lam et al\.“Learning skillful medium\-range global weather forecasting”In*Science*382\.6677American Association for the Advancement of Science, 2023, pp\. 1416–1421
- \[11\]Hans Hersbach et al\.“The ERA5 global reanalysis”In*Quarterly Journal of the Royal Meteorological Society*146\.730Wiley, 2020, pp\. 1999–2049
- \[12\]Vivien Sainte and Loic Landrieu“Lightweight Temporal Self\-attention for Classifying Satellite Images Time Series”In*Advanced Analytics and Learning on Temporal Data: 5th ECML PKDD Workshop, AALTD 2020, Ghent, Belgium, September 18, 2020, Revised Selected Papers*12588, Lecture Notes in Computer Science, 2020, pp\. 171–181Springer Nature
- \[13\]Fanny Lehmann et al\.“Finetuning a Weather Foundation Model with Lightweight Decoders for Unseen Physical Processes”In*arXiv preprint arXiv:2506\.19088*, 2025URL:[https://arxiv\.org/abs/2506\.19088](https://arxiv.org/abs/2506.19088)
- \[14\]Leo Mercer, Josh Burke and Sue Rodway\-Dyer“Towards Improved Cost Estimates for Monitoring, Reporting and Verification of Carbon Dioxide Removal”, 2024URL:[https://www\.lse\.ac\.uk/granthaminstitute/publication/towards\-improved\-cost\-estimates\-for\-monitoring\-reporting\-and\-verification\-of\-carbon\-dioxide\-removal/](https://www.lse.ac.uk/granthaminstitute/publication/towards-improved-cost-estimates-for-monitoring-reporting-and-verification-of-carbon-dioxide-removal/)
- \[15\]Ian\. Power et al\.“Enhanced Rock Weathering Has Greater Promise as a Sustainable Farming Practice than a CO2Removal Technology”In*Proceedings of the National Academy of Sciences*123\.26, 2026, pp\. e2602421123DOI:[10\.1073/pnas\.2602421123](https://dx.doi.org/10.1073/pnas.2602421123)
- \[16\]L\. Poggio et al\.“SoilGrids 2\.0: producing soil information for the globe with quantified spatial uncertainty”In*SOIL*7\.1, 2021, pp\. 217–240DOI:[10\.5194/soil\-7\-217\-2021](https://dx.doi.org/10.5194/soil-7-217-2021)
- \[17\]M\. Rodell et al\.“The Global Land Data Assimilation System”In*Bulletin of the American Meteorological Society*85\.3, 2004, pp\. 381–394DOI:[10\.1175/BAMS\-85\-3\-381](https://dx.doi.org/10.1175/BAMS-85-3-381)
- \[18\]Balázs\. Fekete, Charles\. Vörösmarty and Wolfgang Grabs“High\-resolution fields of global runoff combining observed river discharge and simulated water balances”In*Global Biogeochemical Cycles*16\.3, 2002, pp\. 15\-1–15\-10DOI:[10\.1029/1999GB001254](https://dx.doi.org/10.1029/1999GB001254)
- \[19\]Cort\. Willmott and Kenji Matsuura“Terrestrial air temperature and precipitation: Monthly and annual time series \(1950–1999\)” Center for Climatic Research, University of Delaware, 2001URL:[http://climate\.geog\.udel\.edu/~climate/](http://climate.geog.udel.edu/~climate/)
- \[20\]Pasquale Borrelli et al\.“An assessment of the global impact of 21st century land use change on soil erosion”In*Nature Communications*8, 2017, pp\. 2013DOI:[10\.1038/s41467\-017\-02142\-7](https://dx.doi.org/10.1038/s41467-017-02142-7)
- \[21\]FAO and IIASA“Harmonized World Soil Database version 2\.0: Documentation and User Guide”, 2023URL:[https://www\.fao\.org/3/cc3823en/cc3823en\.pdf](https://www.fao.org/3/cc3823en/cc3823en.pdf)
- \[22\]Maosheng Zhao, Faith Heinsch, Ramakrishna\. Nemani and Steven\. Running“Improvements of the MODIS terrestrial gross and net primary production global data set”In*Remote Sensing of Environment*95\.2, 2005, pp\. 164–176DOI:[10\.1016/j\.rse\.2004\.12\.011](https://dx.doi.org/10.1016/j.rse.2004.12.011)
- \[23\]Robert Carver and Alex Merose“ARCO\-ERA5: an analysis\-ready cloud\-optimized reanalysis dataset”In*American Meteorological Society Meeting Abstracts*103, 2023, pp\. 4A–1URL:[https://ams\.confex\.com/ams/103ANNUAL/meetingapp\.cgi/Paper/415842](https://ams.confex.com/ams/103ANNUAL/meetingapp.cgi/Paper/415842)
- \[24\]Chi Wang, Qingyun Wu, Markus Weimer and Erkang Zhu“FLAML: A fast and lightweight automl library”In*Proceedings of Machine Learning and Systems*3, 2021, pp\. 434–447

## Appendix AAppendix

### A1Architecture Diagram

![Refer to caption](https://arxiv.org/html/2610.28973v1/fig_architecture.png)Figure 2:TerraNova\-ERW architecture\. A frozen weather foundation\-model encoder \(Aurora 1\.5 by default; Pangu\-Weather, SFNO, FCN3, and GraphCast are drop\-in alternatives\) maps ERA5 inputs to a global monthly latent, from which local \(latent cell at the field location\), regional \(3×33\{\\times\}3mean\), and global \(map mean\) features are pooled and aggregated over the observation window by a time\-weighted trapezoidal mean\. The aggregated 3072\-d vector is projected to 64\-d and concatenated with five 64\-d MLP\-encoded site\-specific feature groups; a three\-layer MLP head predicts the carbon dioxide removal \(CDR\) rate\. Only the projection, site\-specific feature encoders, and head are trained\.
### A2Feature Variables

Table 2:Summary of input features and target variables used, organized by category\. ERA5 variables are not consumed directly: they are fed to the frozen Aurora 1\.5 encoder, and only its pooled backbone latent enters the predictor\. The ERA5 block lists Aurora 1\.5’s complete input set: 18 surface fields plus a computed insolation channel, 36 static fields, and 5 atmospheric variables on 13 pressure levels \(50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000 hPa\), for 120 fields in total\. Each field is supplied at two history timesteps \(00:00 and 06:00 UTC on the first of the month\)\. Parenthesized codes are the ERA5 / Aurora short names\. If a given variable isn’t available from the original dataset, it is obtained from an external source as shown below\. Spatial and temporal variables are listed in physical units but are normalized before use \(latitude by90∘90^\{\\circ\}, longitude to\[−1,1\]\[\-1,1\], area and interval bylog⁡\(1\+x\)\\log\(1\+x\), month and day of year by sine/cosine pairs\)\.CategoryVariableUnitsExternal SourceERA5Surface \(18 fields \+ insolation\):2m Temperature \(2t\)KERA5\[[11](https://arxiv.org/html/2610.28973#bib.bibx11)\]2m Dewpoint Temperature \(2d\)KERA5Skin Temperature \(skt\)KERA5Soil Temperature, level 1 \(stl1\)KERA5Volumetric Soil Water, level 1 \(swvl1\)m3/m3ERA510m U\-wind \(10u\)m/sERA510m V\-wind \(10v\)m/sERA5100m U\-wind \(100u\)m/sERA5100m V\-wind \(100v\)m/sERA5Mean Sea Level Pressure \(msl\)PaERA5Surface Pressure \(sp\)PaERA5Total Column Water Vapour \(tcwv\)kg/m2ERA5Total Cloud Cover \(tcc\)fractionERA5Low Cloud Cover \(lcc\)fractionERA5Medium Cloud Cover \(mcc\)fractionERA5High Cloud Cover \(hcc\)fractionERA5Sea\-Ice Cover \(ci\)fractionERA5Snow Depth \(sd\)m w\.e\.ERA5InsolationW/m2Computed\[[6](https://arxiv.org/html/2610.28973#bib.bibx6)\]Static \(36 fields\):Land\-Sea Mask \(lsm\)fractionERA5Surface Geopotential \(z\)m2/s2ERA5Std\. Dev\. of Orography \(sdor\)mERA5Std\. Dev\. of Filtered Subgrid Orography \(sdfor\)mERA5Angle of Sub\-gridscale Orography \(anor\)radERA5Anisotropy of Sub\-gridscale Orography \(isor\)–ERA5Slope of Sub\-gridscale Orography \(slor\)–ERA5Lake Cover \(cl\)fractionERA5Lake Depth \(dl\)mERA5High Vegetation Cover \(cvh\)fractionERA5Low Vegetation Cover \(cvl\)fractionERA5Soil Type, one\-hot: slt\_0–slt\_7 \(8\)–ERA5High Vegetation Type, one\-hot: tvh\_\{0,3,4,5,6,18,19\} \(7\)–ERA5Low Vegetation Type, one\-hot: tvl\_\{0,1,2,7,9,10,11,13,16,17\} \(10\)–ERA5Atmospheric \(5 variables×\\times13 levels\):Temperature \(t\)KERA5U\-wind \(u\)m/sERA5V\-wind \(v\)m/sERA5Specific Humidity \(q\)kg/kgERA5Geopotential \(z\)m2/s2ERA5Static SoilSand Content%SoilGrids\[[16](https://arxiv.org/html/2610.28973#bib.bibx16)\]Silt Content%SoilGridsClay Content%SoilGridsDynamic SiteSoil pHpHSoilGridsSoil Cation Exchange Capacitycmol/kgSoilGridsSoil MoisturemmGLDAS\[[17](https://arxiv.org/html/2610.28973#bib.bibx17)\]Soil Runoffmm/dayUNH/GRDC\[[18](https://arxiv.org/html/2610.28973#bib.bibx18)\]Evapotranspirationmm/day\[[19](https://arxiv.org/html/2610.28973#bib.bibx19)\]Soil Organic Carbong/kgSoilGridsSoil Erosionkg/ha/dayGLOSEM\[[20](https://arxiv.org/html/2610.28973#bib.bibx20)\]Soil Porosity%GLDASNitrification Ratekg/ha/day–Base Saturation%HWSD\[[21](https://arxiv.org/html/2610.28973#bib.bibx21)\]Net Primary ProductivitygC/m2/dayMODIS\[[22](https://arxiv.org/html/2610.28973#bib.bibx22)\]Precipitationmm/yearERA5AmendmentOxides:CaO%–MgO%–Na2O%–K2O%–MnO%–SO3%–P2O5%–Minerals:Quartz%–Calcite%–Actinolite%–Piemontite%–Chlorite%–Titanite%–Muscovite%–Albite%–Olivine%–Plagioclase%–Augite%–Orthopyroxene%–Magnetite%–Ilmenite%–Apatite%–Physical:Reactive Surface Aream2/gCalculated if unavailableApplication Ratet/haCalculated if unavailableSpatialLatitude∘–Longitude∘–Field Areaha–TemporalObservation Intervaldays–Month \(sine/cosine\)––Day of Year \(sine/cosine\)––Yearyr–Amendment Agedays–ERW Trial Indicatorbinary–TargetCDR Ratet CO2e/ha/yrERW: measured∗; Natural Weathering: from source or calculated if unavailable∗The CDR rate labels for the operational ERW data are empirically derived and consider the near\-field zone\. Each georeferenced point is soil\-sampled over a near\-surface credit horizon at a pre\-application baseline and at one or more later credit campaigns, and the samples undergo geochemical analysis\. For each consecutive pair of samplings, removal is quantified by a solid\-phase cation mass balance: an immobile trace element \(e\.g\., Zr, Ti\) constrains the rock mass present at the point, and the control\-corrected loss of rock\-derived cations is converted to CO2removal at 2 mol CO2per mol divalent cation, with corrections for background weathering, strong\-acid \(fertilizer\-derived\) weathering, and plant uptake\. We use near\-field CDR to denote the initial CO2uptake from cation release within the sampled soil horizon\. Project emissions from quarrying, grinding, transport, and spreading are not deducted from the label: they do not vary at the sampling\-point scale and are accounted for at issuance rather than in the spatial rate surface\. Downstream loss and reversal factors are likewise not applied; unlike project emissions these do vary by catchment, so the labels are an upper bound on net durable removal\. Labels are annualized by the interval length\.

Table 3:Frozen weather encoders implemented and evaluated, expressed relative to the Aurora 1\.5 input set of Table[2](https://arxiv.org/html/2610.28973#A1.T2)\. Aurora 1\.5 is the reference: Pangu\-Weather, SFNO and FCN3 require strict subsets of its 120 fields and are therefore served from the same cached ERA5 tensors by variable remap alone\. GraphCast is the sole exception: it additionally requires vertical velocity \(ww\) on all 13 pressure levels, which Aurora 1\.5 does not use, so it is fetched separately from ARCO\-ERA5\[[23](https://arxiv.org/html/2610.28973#bib.bibx23)\]\. It also declares 6\-hourly total precipitation \(tp06\), but that field is target\-only for the operational checkpoint and never a network input, so it is zero\-filled\.
### A3Backbone Latent Extraction Details

Each backbone is run pretrained and frozen on primarily ERA5, plus any additional data the given model requires \(Table[3](https://arxiv.org/html/2610.28973#A1.T3)\)\. At each site, we pool the extracted latent field at three fixed spatial scales and concatenate them:*local*, the single latent cell containing the site \(site\-scale weather\);*regional*, the mean of the surrounding3×33\{\\times\}3window \(∼\\sim±\\pm112 km; mesoscale context\); and*global*, the whole\-map planetary scale; yielding a3​db=dm3d\_\{b\}=d\_\{m\}dimensional embedding per month \(e\.g\.,dm=3072d\_\{m\}\{=\}3072for Aurora 1\.5’sdb=1024d\_\{b\}\{=\}1024\)\. Fixed scales are used rather than field\-area\-adaptive windows because fields are sub\-cell at the latent grid resolution\. Each backbone thus contributes only admd\_\{m\}\-dimensional per\-month embedding \(Aurora 1\.5 = 3072, Pangu = 1152, GraphCast = 1536, SFNO = 1152, FCN3 = 168\)\.

Table[4](https://arxiv.org/html/2610.28973#A1.T4)details, for each weather FM backbone, the layer from which latents are extracted, the native latent dimensiondbd\_\{b\}, and the resulting per\-month embedding dimensiondmd\_\{m\}\.

Table 4:Latent extraction details for each weather FM backbone\.
### A4Temporal Aggregation Methods

TheTTmonthly embeddings spanningtit\_\{i\}toti\+1t\_\{i\+1\}, i\.e\.\{z1,…,zT\}\\\{z\_\{1\},\\dots,z\_\{T\}\\\}, must be reduced to a fixed\-size representation\. We explore three methods:*\(1\) Trapezoidal Integration*, a fixed, parameter\-free aggregationz¯=1T−1​∑t=1T−1zt\+zt\+12\\bar\{z\}=\\frac\{1\}\{T\-1\}\\sum\_\{t=1\}^\{T\-1\}\\frac\{z\_\{t\}\+z\_\{t\+1\}\}\{2\}approximating the time\-integral of the latent trajectory, imposing the inductive bias that every month contributes;*\(2\) Supervised LTAE*, a Lightweight Temporal Attention Encoder\[[12](https://arxiv.org/html/2610.28973#bib.bibx12)\]trained end\-to\-end with the decoder to learn month\-level importance weights \(e\.g\., up\-weighting wet seasons that dominate dissolution\); and*\(3\) Self\-supervised LTAE*, the same architecture trained with a masked\-month reconstruction objective on unlabeled latent sequences, testing whether temporal structure can be learned without weathering\-rate labels\.

Trapezoidal integration was selected for all final experiments after experimentation\.

### A5LOFO Protocol Details

Of the 101 location components, the 68 components with≥\\geq10 measurements each receive a fold \(5,062 of 5,192 rows\)\. Each fold trains with field\-disjoint ERW validation, early\-stops on validation MAE \(validation RMSE for TerraNova\-ERW\), ensembles seeds 42/43/44, and predicts the held\-out component\. Scalers are fit on training data only; TerraNova\-ERW’s three location inputs bypass the per\-fold scaler\. Four of the 68 components contain only control rows, so treated\-row metrics cover 64 components \(51 US, 13 Brazil\)\.

### A6Full Results

Each experimental configuration varies one design choice relative to the shared protocol, so each contrast isolates one effect: one configuration for each backbone \(which weather FM backbone transfers best\); joint training with natural weathering \(the value of cross\-regime pretraining\); raw ERA5 inputs in place of the latent \(the contribution of the FM encoding\); row\-shuffled ERA5 \(real weather vectors, zero information, so the raw\-vs\-shuffled contrast isolates the value of weather information itself\); input\-group ablation, withholding the amendment and site\-specific feature branches in turn \(which inputs are actually required, and what is still recoverable when site characterization is unavailable; Appendix[A7](https://arxiv.org/html/2610.28973#A1.SS7)\); and AutoML \(whether deep learning is needed\)\. Folds, scalers, head, seeds, early stopping, and scoring are identical across configurations, except for TerraNova\-ERW’s two refinements \(unstandardized location inputs and RMSE early stopping; Appendix[A5](https://arxiv.org/html/2610.28973#A1.SS5)\)\.

Table 5:Complete LOFO results on treated rows with a 1000\-resample marginal bootstrap over held\-out components on MAE\. US pipeline: 55 components \(51 with treated rows\); Brazil pipeline: 13 components; site==field\-mean prediction vs\. field\-mean truth\.Boldmarks the best value in each block\. For Aurora 1\.5, three variants are explored: the latents are projected to a narrower 64d and a wider 512d \(note the latents are usually down\-projected to 192d before fusing with site\-specific features\), and a variant with natural weathering co\-training on 192,511 natural \+∼\\sim4,193 ERW rows per fold, evaluated on the identical ERW\-only folds\. The 64d row is TerraNova\-ERW, which additionally uses unstandardized location inputs and RMSE early stopping \(Appendix[A5](https://arxiv.org/html/2610.28973#A1.SS5)\)\. FLAML\[[24](https://arxiv.org/html/2610.28973#bib.bibx24)\]: AutoML run with all site\-specific features\. Search over GBM and tree\-ensemble families; ExtraTrees was selected most often \(32/68 folds\), followed by LightGBM \(21/68\)\.Table[6](https://arxiv.org/html/2610.28973#A1.T6)breaks the all\-ERW treated point\-level results out per seed\. The FM encodings differ from raw ERA5 less in average performance than in reliability: raw ERA5 ranking collapses on one of its three seeds \(Spearman−0\.02\-0\.02on seed 43\), and Pangu shows the same failure mode on two seeds, while the Aurora \(base, 64d, 512d\), SFNO, FCN3, and GraphCast encodings stay within a narrow band \(\+0\.16\+0\.16to\+0\.28\+0\.28\) on all three\. Shuffled ERA5 sits stably near zero, as expected for the no\-information reference\. Since an operational MRV model is retrained as data accrues, a configuration whose ranking skill depends on the training seed is not deployable even if its ensemble looks competitive\.

Table 6:Per\-seed pooled results, all\-ERW treated rows, point level \(n=4244n\{=\}4244per seed; seeds 42/43/44\)\. MAE in t CO2e/ha/yr;ρ\\rhois the Spearman rank correlation and sd\(ρ\\rho\) is its standard deviation across the three seeds\.Boldmarks the best value in each column;underlinemarks seed\-level ranking collapses \(ρ≤0\.05\\rho\\leq 0\.05\)\. FLAML is omitted\.Table 7:Natural\-weathering validation performance of the same TerraNova\-ERW architecture \(64d decoder head on the frozen Aurora 1\.5 encoder\), trained and evaluated exclusively on natural weathering chronosequences: 173,070 training sequences and a location\-disjoint validation set of 19,441 sequences, with no ERW data seen during training or validation\. Values are the mean over five seeds; MAE and RMSE are stable across seeds \(±0\.001\\pm 0\.001\)\. The target \(CDR, t CO2e/ha/yr\) spans 0\.185 to 0\.686 over the training sequences and 0\.211 to 0\.681 over validation, with validation mean 0\.254 and standard deviationσ=0\.081\\sigma=0\.081\. The deployable constant predictor, the training\-set mean of 0\.375, gives the reference row; against it the model scores∼\\sim90% lower MAE and RMSE\. This run uses a similar configuration as the other experiments \(standardized location inputs, MAE early stopping\)\.
### A7Input Ablation: Site Selection with Limited Features

We rerun the identical 68\-fold LOFO protocol over four input configurations, holding folds, decoder, early stopping, and scoring fixed\. Evaluation covers operational US and Brazil ERW rows only \(5,062 rows; 4,244 treated, 818 control\); the configurations are Aurora latents alone \(no amendment and site\-specific features\), Amendment only \(24 variables\), Aurora\+\+Amendment, and Aurora\+\+full Site Specific Features\. All four configurations are the Aurora 1\.5 \(base\) configuration, so the full\-input corresponds to the Aurora 1\.5 \(base\) row of Appendix Table[5](https://arxiv.org/html/2610.28973#A1.T5), not TerraNova\-ERW; TerraNova\-ERW is included in Table[8](https://arxiv.org/html/2610.28973#A1.T8)as a reference row for comparison\.

Table 8:Input ablation on treated rows, same 68 LOFO folds as Table[1](https://arxiv.org/html/2610.28973#S3.T1)\. MAE and RMSE in t CO2e/ha/yr\. Site metrics average predictions and targets within each held\-out component before scoring\.Aurora latents alone are not deployable: point Spearman is−0\.12\-0\.12, so the model fails to order measurements within an unseen field, even though its site\-level error is competitive\. Adding the amendment variables is the single change that restores rank signal at both levels, lifting point Spearman to\+0\.11\+0\.11and site Spearman to\+0\.32\+0\.32\. Relative to the full stack, Aurora\+\+Amendment gives up little at the point level \(MAE 3\.27 vs\. 3\.21, Spearman\+0\.11\+0\.11vs\.\+0\.24\+0\.24\) while ordering held\-out fields far better \(site Spearman\+0\.32\+0\.32vs\.\+0\.13\+0\.13, site MAE 1\.80 vs\. 2\.04\); a plausible reading is that soil and spatial covariates are informative about within\-field variation but are partly absorbed as field identity, which does not transfer to an unseen component\. One caveat: the site\-level advantage is concentrated in ranking rather than error, with site MAE staying near 1\.8–2\.0 t CO2e/ha/yr across configurations, so these results support site prioritization rather than model\-based quantification\.

### A8Sampling\-Budget Simulation

#### Protocol\.

A*field*is one location component and a*point*is one measured CDR interval on that field\. Each Monte\-Carlo draw revealskkrandomly chosen points per field, forms an estimate of the field’s mean, and scores it against the mean of the points that were*not*revealed, so no estimate is graded on its own inputs\. Model predictions are the LOFO out\-of\-fold predictions of TerraNova\-ERW \(Aurora 1\.5, 64d latent projection\) \(68 folds, seeds 42–44 averaged\); priors and variance components are likewise fitted leaving the scored field out\. Estimators are:*sampling only*\(y¯k\\bar\{y\}\_\{k\}\);*sampling shrunk*, thekk\-point mean shrunk toward the leave\-one\-field\-out prior with the error\-optimal weightw=τ2/\(τ2\+σ2/k\)w=\\tau^\{2\}/\(\\tau^\{2\}\+\\sigma^\{2\}/k\)from the between\- and within\-field variance components, but with no model, isolating the value of shrinkage itself;*model only*\(f^\\hat\{f\}on the unrevealed points\);*hybrid \(error weights\)*, the same shrinkage applied to the model’s residuals, i\.e\. the error\-optimal blend of measurements and model; and*hybrid \(rank weights\)*,α​z​\(y¯k\)\+\(1−α\)​z​\(f^\)\\alpha\\,z\(\\bar\{y\}\_\{k\}\)\+\(1\-\\alpha\)\\,z\(\\hat\{f\}\)with oneα\\alphaper budget tuned for the ranking objective on the even\-numbered Monte\-Carlo draws, with every estimator scored on the odd\-numbered draws only \(α\\alpharises from 0\.65 at one point to 0\.85 at five\)\. Ranking fixes the portfolio at the 52 fields with≥13\\geq 13treated points so every budget orders the same set; error tiers use every field with enough unrevealed points left to grade, sonnfalls withkk\(64 fields atk=1k\{=\}1, 63 atk=2k\{=\}2, 62 atk=3k\{=\}3, 52 atk=5k\{=\}5, 39 atk=10k\{=\}10, 25 atk=20k\{=\}20\) and the fields surviving at highkkare systematically the denser ones\.

Three caveats bound these numbers\. First, the truth is itself an estimate \(the mean of the unrevealed points\), and within\-field noise carries 64% of total variance, which puts a floor under every error and caps attainable rank agreement\. Second, rows at different budgets are not computed over identical field sets \(see the*fields*column of Table[10](https://arxiv.org/html/2610.28973#A1.T10)\)\. Third, the residual “poor” tail is where measurements matter most: with no measured points the model alone leaves 13% of fields in the “poor” tier \(error\>4\{\>\}4t\), and a single measured point per field brings that to 5%, because even one point of local calibration removes a stable per\-field offset in the model\.

Table 9:Ranking 52 fields for deployment potential under a measurement budget ofkkpoints per field \(mean over 100 held\-out Monte\-Carlo draws,±\\pmacross draws\)\. “Top\-10 found” is the share of the 10 genuinely highest\-removal fields appearing in the estimator’s top 10; picking at random gives 19%\.α\\alphais the tuned weight on the measured mean\. The rank\-weighted hybrid leads on Spearman and leads or ties on top\-10 at every budget; the error\-weighted hybrid is strongest on the top\-quartile screen\.Table 10:Share of fields by accuracy of their estimated mean CDR under a budget ofkkpoints per field: typical\|\|error\|\|within 1, 1–2, 2–4, and beyond 4 t CO2e/ha/yr, against field means spanning−5\-5to\+10\+10\. The no\-model, no\-point baseline is thek=0k\{=\}0row of “sampling shrunk”, at 41% great\. The hybrid beats sampling alone at every budget through ten points and matches or beats shrinkage without a model through three points \(shrinkage overtakes on the “great” tier atk=5k\{=\}5, 69% vs\. 63%, while the hybrid retains the lower median error through ten points\)\. What measurements add to the hybrid is mainly the tail, with “poor” falling from 13% atk=0k\{=\}0to 5% with a single point\.Estimatorpointsgreat≤1\\leq 1good 1–2fair 2–4poor\>4\>4fieldsSampling only19%36%41%14%64216%56%24%5%63329%52%19%0%62542%52%6%0%521072%23%5%0%392096%4%0%0%25Sampling shrunk \(no model\)041%22%27%11%64139%31%25%5%64243%25%27%5%63348%32%18%2%62569%25%6%0%521077%18%5%0%392092%8%0%0%25Model only045%22%20%13%64144%23%20%13%64246%22%19%13%63347%21%21%11%62548%23%21%8%521046%31%13%10%392052%36%4%8%25Hybrid \(error weights\)044%25%19%13%64145%30%20%5%64248%29%19%5%63352%27%18%3%62563%27%8%2%521077%18%5%0%392092%8%0%0%25Table 11:Samples required to estimate the field mean within a precision target on the three densest fields, which together hold 32% of all treated points\. Values are interpolated from per\-field Monte\-Carlo error curves \(300 draws per budget\): samples needed by sampling alone vs\. the hybrid, and the resulting share of samples saved\. On Fields B and C the±0\.5\\pm 0\.5t target is met with one sample because the model alone already places their means within 0\.10 and 0\.32 t; this rests on two fields and should be read per field rather than as a typical saving\. Field names are anonymized\.