An uncertainty-aware framework for solid-phase measurement and verification of enhanced weathering

Reliable verification of enhanced weathering as a carbon dioxide removal strategy requires accurate quantification of feedstock dissolution in amended soils. However, spatial heterogeneity introduces significant uncertainty, particularly in sampling designs that rely on sparse or repeated measurements at fixed locations. Here, we develop a probabilistic framework to evaluate how spatial uncertainty in solid-phase geochemical measurements influences the precision of feedstock dissolution estimates derived from an element-element mixing model. We first quantify how variance in soil compositions affects errors in modeled feedstock dissolution and apply distance-based sensitivity analysis to identify the measurement variance thresholds required to achieve desired uncertainty levels. Next, we simulate spatially heterogeneous soil conditions and various composite sampling configurations to identify the optimal sampling strategy likely to meet specified uncertainty criteria. Our findings underscore the necessity of accurately estimating field-scale variance in baseline soil concentrations prior to developing sampling plans. Analysis of data from existing high-density soil sampling campaigns indicates that geochemical variance is likely too high for element-element mixing models to serve as effective near-term constraints on feedstock dissolution. The framework presented here can be further extended to other solid- and multi-phase measurement models for enhanced weathering verification.

1 Introduction

Open-system carbon dioxide removal (CDR) technologies, such as enhanced weathering (EW), are being evaluated as potential keystone climate mitigation strategies. For EW, however, accurate and cost-effective verification remains a major obstacle to assessing large-scale effectiveness (Calabrese et al., 2022; Levy et al., 2024; Lebling et al., 2024; Mercer and Burke, 2023). EW involves amending soils, usually agricultural, with a reactive “feedstock,” such as basalt or Mg-silicate, to shift the alkalinity of the system and effectively dissolve additional CO₂ (Hartmann et al., 2013). Global projections (IPCC, 2018; NASEM, 2019; Beerling et al., 2020) anticipate removals of 0.5–4 Gt CO₂/year through EW, but only a fraction of these removals may be resolvable due to the geochemical “noise” inherent in soil systems (Reershemius and Suhrhoff, 2024; Suhrhoff et al., 2024; Derry et al., 2025).

Quantifying CDR in soil systems is complicated by spatial heterogeneity across scales—from mineral surfaces to landscapes (Lin, 2010; Vereecken et al., 2016)—that hinders attribution of the geochemical changes arising from EW activities. Verifying CDR through EW involves constraining both the rate of feedstock dissolution and the attendant increase in alkalinity and dissolved inorganic carbon (DIC) fluxes to the groundwater and stream network (Almaraz et al., 2022; Clarkson et al., 2024). While aqueous measurements are critical for constraining dissolved fluxes (e.g., Knapp et al., 2023; Hasemer et al., 2024; Rieder et al., 2024) and will thus constitute a major portion of measurement and verification, complementary approaches based on solid-phase mass balance (Reershemius et al., 2023) estimate cumulative feedstock dissolution using assumptions about stoichiometry and baseline compositions to compute potential CDR (Renforth, 2012; Sutherland et al., 2025; puro.earth, 2024). However, given the low initial enrichment of feedstock mass relative to native soil, confidently detecting a dissolution signal beyond a spatially heterogeneous baseline becomes a challenging problem (Suhrhoff et al., 2024), and accurately quantifying the magnitude of such signals can be even more challenging (Wasserstein et al., 2019; Bradford et al., 2023).

One widely considered approach to constraining feedstock dissolution relies on ratios of base cations to immobile elements in the solid phase (Reershemius et al., 2023). The resulting depletion and mixing equations (Equations 1, 2; Supplementary material 1) require analysis of multiple soil samples in sequence, from soil (baseline) to the initial mixture (soil + feedstock) to weathered compositions over multiple time points (Brimhall and Dietrich, 1986; Reershemius et al., 2023). These measurements are used to estimate the true fraction of feedstock dissolved (f_d), calculated as

$\begin{array}{ll}{f}_{\text{d}}=1-\frac{{\left[\text{M}\right]}_{\text{mix}}-{\left[\text{M}\right]}_{\text{bl}}}{{\left[\text{M}\right]}_{\text{mix}}^0-{\left[\text{M}\right]}_{\text{bl}}}& (1)\end{array}$

where [M] is base cation concentration of the baseline (_bl), initial mixture (${(}_{\text{mix}}^0)$), and weathered mixture (_mix). f_d gives the potential CDR before any correction for losses deeper in the soil profile or from down-gradient atmospheric equilibration (Renforth, 2012; Bertagni and Porporato, 2022; Dietzen and Rosing, 2023), hence an accurate estimate ($\widehat{f_{\text{d}}}$) is the focus here. In this approach, ${\left[\text{M}\right]}_{\text{mix}}^0$ is calculated by measuring an immobile tracer (T) and using the following element-element mixing equation with baseline and feedstock (_fs) endmembers,

$\begin{array}{ll}{\left[\text{M}\right]}_{\text{mix}}^0={\left[\text{M}\right]}_{\text{bl}}+\frac{({\left[\text{T}\right]}_{\text{mix}}-{\left[\text{T}\right]}_{\text{bl}})({\left[\text{M}\right]}_{\text{fs}}-{\left[\text{M}\right]}_{\text{bl}})}{{\left[\text{T}\right]}_{\text{fs}}-{\left[\text{T}\right]}_{\text{bl}}}& (2)\end{array}$

Here, heterogeneity can violate the basic assumptions of the mixing model if samples are not representative of the same geochemical “system” (Figure 1). Uncertainty in these assumptions was initially considered in the context of analytical variance (Reershemius et al., 2023). Building on this, Suhrhoff et al. (2024) evaluated how overall measurement variance (aggregate sampling and analytical variance), soil and feedstock compositions, and application rate influence the detectability of T upon amendment, as well as how the magnitude of f_d affects the ability to detect changes in M. They concluded that measurement variance in T and the contrast between soil and feedstock compositions determine whether EW signals can be detected. However, given the inherent spatial heterogeneity in soils, it remains unclear whether scalable sampling plans can yield the low measurement variances required. Moreover, the detection of changes in elemental abundances does not guarantee accurate mixing calculations: even small measurement variances can yield apparent mixtures that fall outside of the theoretical mixing space (Verma, 1998).

Figure 1

Figure 1. Illustration of how spatial uncertainty introduces error into feedstock dissolution calculations. The underlying mixing model assumes that baseline and mixture samples represent the same system or control volume. In reality, it is difficult to re-sample at the exact same location. For a presumed location (blue cross or blue circle) the actual sample location (red square or red circle) may be slightly different due to, e.g., positioning error. Furthermore, even if the exact same location is re-sampled, tillage or erosion may create a transient baseline that induces error, as might inconsistent sampling, preparation, or analytical techniques. Composite sampling can help reduce error from measurement variance by making it more likely to sample and re-sample similar control volumes.

These concerns highlight the need for a systematic approach to account for spatial heterogeneity in the design of EW verification protocols. To address this, we present a framework for incorporating spatial uncertainty into solid-phase measurements, starting with definition of overall uncertainty requirements (e.g., a maximum error and minimum confidence in $\widehat{f_{\text{d}}}$), followed by variance propagation and sensitivity analysis to help define corresponding measurement variance requirements (e.g., a maximum measurement variance and minimum confidence). We then use stochastic sampling simulations to infer a sampling approach that minimally meets these requirements and conclude with estimate and uncertainty reporting. This is detailed through the following steps:

1. Defining the uncertainty requirements and measurement model; the latter includes explicit relationships between input and output uncertainty using hierarchical Bayesian principles.

2. Determining the maximum measurement variances that fulfill the overall uncertainty requirements using variance propagation and sensitivity analysis.

3. Defining the measurement variance requirements and sampling model; the latter involves stochastic simulation of spatial fields and composite sampling plans.

4. Designing a sampling plan that minimally meets the measurement variance requirements using the sampling model and sensitivity analysis; if infeasible, reconsider the overall uncertainty requirements or measurement model.

5. Reporting the final estimate and overall uncertainty, with traceable and reproducible uncertainty quantification.

2 Integrated methods and results

Agricultural fields used for EW are typically selected for logistical reasons such as accessibility, rather than a detailed understanding of soil properties and heterogeneity. Still, the quantification agent (e.g., project developer, research scientist, practitioner) needs to estimate dissolved feedstock and resulting CDR with sufficient accuracy, ideally using as little sampling as possible. Soil sampling might include individual point samples with a defined footprint, or composite samples comprised of multiple point samples homogenized to represent a specified area. EW verification requires these points or areas to be re-sampled multiple times, but imprecise positioning devices and shifting baselines (e.g., due to tillage or erosion) make it difficult to precisely re-sample the same location or control volume, introducing error into the CDR calculation. Early characterization of soil variability can therefore play a central role in site selection and monitoring design, increasing the likelihood of obtaining precise and defensible CDR estimates.

In this example, we illustrate how such early characterization informs the full workflow of uncertainty assessment and sampling design. Although we assume a basaltic feedstock, the approach is generalizable to any amendment. For mixing-based solid-phase verification of EW, the measurement model is grounded in Equations 1, 2, which are solved for f_d based on measured baseline concentrations and either (i) the soil-feedstock mixture immediately after application to constrain application rates or (ii) the mixture after weathering has occurred (Reershemius et al., 2023); throughout this framework, we use spatial averages as inputs to the f_d calculation, rather than using individual measurements as inputs and then spatially averaging f_d. Field trials (Kantola et al., 2023; Beerling et al., 2024; Guo et al., 2023; Dietzen and Rosing, 2023; Larkin et al., 2022; Wang et al., 2024; Skov et al., 2024; Holzer et al., 2023; McDermott et al., 2024; Dupla et al., 2024; Sokol et al., 2024) report application rates of 5–100 tons of feedstock per hectare (t/ha), producing only 0.1%–3% mass enrichment when mixed in the top 0.15 m of the soil (Supplementary material 2), potentially introducing challenges in quantifying small dissolution signals amid heterogeneous baselines. Another key factor is the chemical differentiation between the feedstock and the baseline, which we express using feedstock-to-baseline ratios of mean cation (${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}$) and tracer (${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}$) concentrations.

The sampling model specifies how soil will be collected for analysis, through discrete point samples or homogenized composite samples, and provides the basis for simulation of the measurement variance associated with different sampling approaches. Because the measurement model depends on the baseline, and the sampling plan is typically fixed after feedstock application (e.g., through paired sampling), it is essential to establish a robust baseline sampling strategy and/or modify approaches as more data are collected. Finally, the uncertainty requirements reflect operational constraints, such as the need to support a compensatory claim or regulatory requirement, and define the probability that the estimated f_d will lie within an acceptable error margin of its true value.

Together, these elements motivate the five-step framework that follows, which links uncertainty requirements, measurement modeling, spatial heterogeneity, and sampling design into a coherent process for evaluating, optimizing, and reporting EW monitoring strategies.

2.1 Defining the uncertainty requirements and measurement model (Step 1)

Quantifying carbon removal from EW begins with establishing how much uncertainty is acceptable and expressing the measurement process in a form that makes those uncertainties explicit. This foundation is critical because the dominant source of error, spatial variability in soils, is unknown before sampling and can easily exceed what a measurement approach is designed to handle. As a result, a deployment cannot assume it will meet prescribed accuracy or confidence thresholds; it must demonstrate that these thresholds are achievable under the site's actual variability.

Step 1 sets up this demonstration. It defines the uncertainty bounds that a credible requirement must satisfy and formulates the measurement approach in a way that reveals how spatial variability and operational choices drive error. This structure supports Step 2, where variance is propagated through the measurement function and we evaluate when a given configuration can realistically meet the uncertainty requirements.

2.1.1 Uncertainty requirements

The uncertainty requirements, which might be defined by a standards development organization, here encapsulate:

• ϵ_max, the maximum relative error in $\widehat{f_{\text{d}}}$ and thus potential CDR.

• p_min, the minimum probability (i.e., confidence) that the relative error in $\widehat{f_{\text{d}}}$ is less than ϵ_max.

For instance, a quantification standard might require 90% confidence (p_min = 0.9) that the reported CDR is within 10% of the true value (ϵ_max = 0.1) (e.g., Verra, 2012). There is no a priori guarantee, however, that any particular field deployment can meet these requirements for a specific site, due to the inherent variability in measurement conditions and system parameters.

We use ϕ to represent the outcome where the uncertainty requirements are fulfilled. ϕ is achieved when p(ϵ ≤ ϵ_max), the probability that the error is sufficiently small, meets the specified confidence threshold p_min. Formally:

$\begin{array}{ll}p(ϵ\le {ϵ}_{\text{max}})\ge p_{\text{min}}\Rightarrow \varphi ,p(ϵ\le {ϵ}_{\text{max}})<p_{\text{min}}\Rightarrow \overline{\varphi },& (3)\end{array}$

where $\overline{\varphi }$ indicates the requirements are not fulfilled.

2.1.2 Measurement model and parameter set

Quantifying whether a measurement approach can meet the uncertainty requirements requires more than an error bound; it requires a model that makes each source of uncertainty explicit. In this context, the feedstock dissolution calculation (Equations 1, 2) is influenced not only by measured concentrations, but also by how spatial heterogeneity, measurement variance, and operational choices combine to perturb those measurements. We therefore define a measurement model that captures these effects and can be systematically explored in Step 2.

To calculate p(ϵ ≤ ϵ_max), the probability of sufficiently small error, for different measurement approaches, we define a measurement model with parameter space Θ. This model includes:

• a measurement function (e.g., feedstock dissolution calculation, Equations 1, 2),

• input parameters (e.g., M and T concentrations in the baseline, feedstock, and mixture),

• the function response (e.g., fraction of feedstock dissolved, f_d),

• measurement variance of each input parameter (e.g., spread of possible M and T measurement values given the point or composite sampling scheme), and

• operational parameters (e.g., feedstock-baseline differentiation, application rate, true fraction dissolved).

For this measurement function (Equations 1, 2), the input parameters are [M]_bl, [M]_fs, [M]_mix, [T]_bl, [T]_fs, and [T]_mix. Due to spatial heterogeneity, baseline and mixture measurements may be highly variable, interrupting their assumed comparability (Figure 2). To account for the impact of spatial heterogeneity, we consider each input parameter to have a distribution of possible measurement values, characterized by a mean (μ) and relative variance (υ), which we will propagate through the measurement function in Step 2. For brevity, mean μ values are written as natural logarithms (ln); a reference table for ln values is provided in Supplementary material 2. While measurement variance represents aggregate spatial and analytical uncertainty, we only consider spatial uncertainty in this study, as analytical uncertainty can be made negligible if necessary (Reershemius et al., 2023). Accordingly, we define μ and υ for the measurement distributions as follows:

Figure 2

Figure 2. Illustration of how measurement variance (υ) in cation (M) or tracer (T) concentrations may interrupt the assumed comparability of baseline and mixture samples to differing degrees depending on sampling approach. Composite approaches tend to result in re-sampling of more similar control volumes than point sampling does, reducing the measurement variance υ. The goal of Step 2 is to calculate how small υ_M and υ_T must be, or how narrow each probability density function of measurement values must be, so that the control volumes that comprise baseline samples are not too different from the control volumes that comprise mixture samples (i.e., the baseline sample is similar enough to the baseline that actually comprises the mixture sample).

Measurement means (μ): we set ${\mu }_{\text{M}}^{\text{bl}}$ and ${\mu }_{\text{T}}^{\text{bl}}$ (baseline means) as constants at the simulation scale. We then specify feedstock-baseline concentration ratios ${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}$ and ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}$ to obtain mean feedstock concentrations. A uniform application rate r_app and uniform fraction dissolved f_d together determine the mean mixture concentrations after amendment and weathering. This f_d also serves as the “true” fraction against which estimation errors are calculated.

Measurement variances (υ): we specify υ_M and υ_T, the (relative) measurement variances for baseline M and T. Due to feedstock mass enrichments of < 3%, we assume the measurement variances of the soil-feedstock mixture are equal to the baseline υ_M and υ_T. The feedstock itself is assumed homogeneous (negligible variance).

Collectively, these parameters comprise the parameter space Θ. Since we want to test the impact of different measurement distributions and operational parameters on fulfilling the uncertainty requirements, we initially consider a wide range of possible values for each parameter (Table 1). Mean baseline cation concentrations (${\mu }_{\text{M}}^{\text{bl}}$) range approximately 3–160,000 ppm ($ln({\mu }_{\text{M}}^{\text{bl}})$ range 1–12), and mean baseline tracer concentrations (${\mu }_{\text{T}}^{\text{bl}}$) range approximately 3–8,100 ppm ($ln({\mu }_{\text{T}}^{\text{bl}})$ range 1–9), considering the 10–90th percentiles of A-horizon soil elemental concentrations across the United States (Smith et al., 2019). Ranges for feedstock-baseline concentration ratios (${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}$, ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}$) are calculated based on six different basalt compositions (Lewis et al., 2021) relative to United States soils (Smith et al., 2019), with a minimum ratio of 1 representing equal elemental concentrations in the feedstock and soil. These ranges are used as bounds for the parameter space Θ. The next step involves random sampling of Θ to rigorously evaluate the individual and joint impacts of each parameter on whether the uncertainty requirements are fulfilled (ϕ) or not ($\overline{\varphi }$).

Table 1

Table 1. Exploratory Θ, uniform distributions used for variance propagation and sensitivity analysis to determine measurement variance requirements.

2.2 Determining maximum measurement variances (Step 2)

With the uncertainty requirements and measurement model defined in Step 1, Step 2 evaluates how measurement variance propagates through the model and determines the maximum variances under which a deployment can still meet the uncertainty requirements. This involves three components: (i) propagating measurement variances through the measurement model using nested Monte Carlo simulation (Step 2.1), (ii) quantifying how individual parameters influence the probability of meeting the uncertainty requirements (Step 2.2), and (iii) identifying deployment-specific limits on measurement variance, υ_max, given operational constraints (Step 2.3).

2.2.1 Propagating measurement variances through the measurement model (Step 2.1)

To evaluate whether a deployment can meet prescribed uncertainty requirements, we must first understand how spatially driven measurement uncertainty affects the calculation of $\widehat{f_{\text{d}}}$. Even small increases in baseline or mixture variability can compound through the measurement model, shifting the relative error in ways that are not obvious from the parameter definitions alone. Step 2.1 quantifies this compounding effect directly by propagating measurement variances through the model under thousands of plausible measurement scenarios.

We use nested Monte Carlo simulation to sample broadly across the parameter space Θ defined in Step 1 and to generate distributions of measurement outcomes for each scenario. This approach captures the full range of interactions between baseline means, measurement variances, feedstock-baseline differentiation, application rate, and true dissolution. The process begins with generating 10⁴ parameter realizations, or samples of Θ. Θ is a uniform multivariate distribution, meaning each parameter range in Table 1 is sampled from uniformly, and each realization represents a possible combination of baseline means (${\mu }_{\text{M}}^{\text{bl}}$, ${\mu }_{\text{T}}^{\text{bl}}$), measurement variances (υ_M, υ_T), and operational parameters (${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}$, ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}$, r_app, f_d). For a given parameter realization, we use these values to construct Gaussian measurement distributions for [M]_bl, [M]_mix, [T]_bl, and [T]_mix, and we sample from these distributions to generate 10⁴ measurement realizations. For each measurement realization, we compute $\widehat{f_{\text{d}}}$ and its relative error (ϵ), such that p(ϵ ≤ ϵ_max) for each parameter realization is the fraction of its measurement realizations where ϵ ≤ ϵ_max.

Results of these simulations indicate that ϕ is highly dependent on keeping the measurement variances below critical thresholds (Figure 3), while other parameters, such as ${\mu }_{\text{M}}^{\text{bl}}$ and ${\mu }_{\text{T}}^{\text{bl}}$, have less impact. Specifically, the distribution of p(ϵ ≤ ϵ_max) shows a clear divide (Figure 3A), indicating that while many realizations achieve ϕ, a significant number fail. The p₁₀-p₉₀ gray shaded regions in Figure 3B illustrate the spread of simulation outcomes across each parameter range; shaded regions that extend above p_min (black dashed line) indicate parameter values for the realizations that achieved ϕ in Figure 3A. Conversely, the unshaded regions above p_min indicate that high υ_M and υ_T and low r_app values will likely not result in ϕ.

Figure 3

Figure 3. Exploratory sensitivity analysis of the influence of means (μ), measurement variances (υ), and operational parameters on p(ϵ ≤ 0.1), the probability that the relative error in $\widehat{f_{\text{d}}}$ is no greater than 10%, using the parameter ranges in Table 1. (A) Shows the response across all 10⁴ realizations, with black dashed lines representing various confidence levels (p_min of 90%, 80%, 70%) of achieving ϕ; the black dashed lines separate the realizations that do (ϕ) and do not ($\overline{\varphi }$) fulfill uncertainty requirements of p(ϵ ≤ 0.1) ≥ 90%, 80%, 70%. (B) Provides the conditional response distribution for each parameter; black solid p₅₀ lines that extend above the black dashed p_min lines indicate where half of the realizations meet the uncertainty requirements and achieve ϕ; areas where the entire gray shaded p₁₀-p₉₀ region is below the black dashed p_min lines indicate where 90% of the realizations did not achieve ϕ, and the opposite of that would be true where the entire gray shaded p₁₀–p₉₀ is above the black dashed lines, but that is not achieved in this exploratory analysis.

The correlations among the expected p(ϵ ≤ ϵ_max) (black p₅₀ lines in Figure 3B) and the individual parameters provide additional insight into sensitivities. Specifically, the conditional distributions of p(ϵ ≤ ϵ_max) show a strong negative correlation with υ_M and υ_T, meaning high measurement variances make it unlikely to achieve ≤ 10% error in $\widehat{f_{\text{d}}}$ (Figure 3B). Intuitively, ${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}$, ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}$, r_app, and true f_d show moderate positive correlations with p(ϵ ≤ ϵ_max), indicating that greater values tend to increase the expected accuracy in $\widehat{f_{\text{d}}}$. Overall, the wide range of outcomes here emphasizes the importance of considering all possible outcomes early in site selection and monitoring design.

2.2.2 Quantifying parameter influence on fulfilling the uncertainty requirements (Step 2.2)

While Step 2.1 shows which combinations of parameters succeed or fail, it does not reveal why certain deployments fulfill the uncertainty criteria and others do not. In practice, quantification agents need to know which parameters most strongly affect the probability of satisfying the requirements so they can prioritize sampling effort, improve operational choices, or redesign deployment strategies. Step 2.2 isolates the influence of each parameter by comparing the conditions under which the requirements are fulfilled to those under which they are not.

To rigorously compare the sensitivity of ϕ to different parameters, we separate the 10⁴ realizations of Θ into one group that does fulfill the uncertainty requirements and one group that does not. This can be represented by partitioning Θ into two conditional distributions, Θ|ϕ and $\Theta |\overline{\varphi }$, and we can analyze the differences between these distributions to determine which parameters most significantly influence the outcome. A common way to quantify such sensitivities (Spear, 1980; Saltelli et al., 2004; Scheidt and Caers, 2009; Fenwick et al., 2014) is to compute the “distance” between Θ|ϕ and $\Theta |\overline{\varphi }$ for each parameter (Figure 4A), normalizing the parameter ranges to [−1, 1] so they do not influence comparison of the distances. The resulting sensitivity rankings (Figure 4B) highlight that ϕ is most influenced by υ_M, υ_T, r_app, and true f_d, and less so by ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}$ and ${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}$. Collectively, this emphasizes the dominant role of measurement variances in determining success.

Figure 4

Figure 4. Distance-based sensitivity calculations for the exploratory sensitivity analysis in Figure 3. (A) Shows threshold-conditional CDFs, i.e., partitions of the entire set of realizations (gray dashed line) into realizations that did (green line) and did not (red line) fulfill uncertainty requirements, with shaded areas to visualize distances between CDFs. (B) Provides a ranking of the parameters according to their influence on ϕ using this distance-based sensitivity metric.

2.2.3 Identifying measurement variance limits using deployment-specific constraints (Step 2.3)

Real deployments do not operate across the full exploratory parameter space. Project developers select feedstock, choose application rates, and have expectations about the likely dissolution behavior at a site. Once these operational choices are fixed, the question becomes much sharper: given the deployment's actual conditions, what is the maximum level of measurement variance that still ensures the uncertainty requirements can be met? Step 2.3 addresses this question by applying the sensitivity framework under deployment-specific constraints.

As an example, for our theoretical deployment we constrain ${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}$ to 38 and ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}$ to 25 (midpoints from Table 1), r_app to 40 tons/ha (median from Table 1), and f_d to 0.7-0.9 (high dissolution scenario). Performing the sensitivity analysis with these constraints (Figure 5) reveals that the expected, or median, p(ϵ ≤ ϵ_max) exceeds p_min for υ_M and υ_T less than approximately 2% (ln(υ) < −4) (Figure 5B). In contrast, the conditional response distributions for other parameters in Figure 5B do not show an expected p(ϵ ≤ ϵ_max) greater than p_min, as each distribution assumes values for all other parameters are randomly chosen from their respective ranges, thus incorporating effects from the entire ranges of υ_M and υ_T. While measurement variances are the primary control here, the application rate r_app and true f_d will become the two most significant parameters after constraining υ_M and υ_T (Figure 5C). For a commercial setting, this suggests that delaying intensive sampling until greater feedstock dissolution has occurred, though also delaying return on investment to the project developer, could be a key feature of profitable operations.

Figure 5

Figure 5. Deployment-specific sensitivity analysis where, relative to the exploratory sensitivity analysis in Figures 3, 4 and parameter ranges in Table 1, we apply constraints to soil-feedstock differentiation (${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}=38$, ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}=25$), application rate (r_app = 40 tons/ha), and expected dissolution (f_d=0.7–0.9) for confidence level p_min = 90%. (A) Shows the updated response across all 10⁴ realizations, (B) the updated conditional response distributions, and (C) the updated ranking of the parameters according to their influence on ϕ.

To determine specific measurement variance limits, we need to account for potential interactions between υ_M and υ_T by analyzing their joint conditional distribution (Figure 6). It is also important to consider which Θ is being used when determining such limits. Using an exploratory Θ, this analysis shows almost no combinations of υ_M and υ_T that achieve ϕ (Figure 6A). Using the constrained Θ, however, indicates the expected outcome is ϕ when both υ_M and υ_T are less than approximately 1.5% (ln(υ) < −4.2) (Figure 6B).

Figure 6

Figure 6. Combinations of base cation measurement variance (υ_M) and immobile tracer measurement variance (υ_T) that result in fulfilling uncertainty requirements (ϕ) of at least 90% likelihood of ≤ 10% error in $\widehat{f_{\text{d}}}$ for (A) loosely constrained, exploratory parameter ranges (Table 1) and (B) constrained parameter ranges for a theoretical deployment where ${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}=38$, ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}=25$, r_app = 40 tons/ha, and f_d=0.7–0.9.

Since, aside from measurement variances, the application rate and true dissolution exert the most influence on error, we infer υ_max for high dissolution (f_d 0.7–0.9) and low dissolution (f_d 0.1–0.3) scenarios with low (10 t/ha) and high (40 t/ha) application rates. For each scenario, we infer υ_max given different error margins (10%, 20%, 30%) and confidence levels (90%, 80%, 70%) (Table 2). Given these four different scenarios of application rate and dissolution, simulated maximum measurement variance υ_max values range from 0.1% to 5%.

Table 2

Table 2. Maximum measurement variances υ_max (%) for a deployment-specific parameter space with constrained feedstock-baseline concentration ratios (${\mu }_{\text{M}}^{\text{fs}}/{\mu }_{\text{M}}^{\text{bl}}=38$, ${\mu }_{\text{T}}^{\text{fs}}/{\mu }_{\text{T}}^{\text{bl}}=25$), high and low application rates (r_app = 10, 40 tons/ha), and high and low dissolution scenarios for different error margins (ϵ_max = 10%, 20%, 30%) and confidence levels (p_min = 90%, 80%, 70%).

The remaining analysis provides insight on combinations of inherent site characteristics and sampling designs likely to adhere to υ_max of 1% and 5% using stochastic simulations of spatial variability and composite sampling.

2.3 Defining the measurement variance requirements and sampling model (Step 3)

Quantifying whether a deployment can meet the uncertainty requirements ultimately depends on whether its sampling design can keep measurement variance below the maximum limits identified in Step 2. However, the true fine-scale spatial distribution of elemental abundances at a field site is never known in advance, and the structure of that heterogeneity, not just its magnitude, strongly shapes the variance obtained from any sampling plan. Step 3 therefore shifts from the measurement model to the sampling problem: given realistic spatial variability, what sampling designs are capable of achieving a measurement variance less than υ_max with at least probability p_min?

If the concentrations everywhere in a field were known at high resolution, one would expect to observe structured lenses or patches, rather than a completely random distribution. Synthetic fields can approximate this structure by drawing on geostatistical models of spatial variability, an approach widely used in contaminant remediation (Christakos and Olea, 1992; Ma and Chang, 2008; Demougeot-Renard et al., 2004; Verstraete and Van Meirvenne, 2008) and in estimating soil organic carbon stocks (Lark, 2009; Bradford et al., 2023; Potash et al., 2025). By sampling these synthetic landscapes, we can evaluate and compare composite sampling schemes across a broad range of plausible spatial conditions.

We simulate fields using a 1,000-by-1,000-cell grid, corresponding to, for example, 0.01 ha (100 m²) at 1-cm resolution, 1 ha (10,000 m²) at 10-cm resolution, or 100 ha (10⁶ m²) at 1-m resolution.

2.3.1 Measurement variance requirements

To illustrate this workflow, we begin by considering a maximum allowable measurement variance υ_max = 1% and report outcomes for both 1% and 5% at the end. We can generally denote [M] or [T] as an arbitrary spatial variable Z. Here, the measurement variance requirements for Z are defined by:

• υ_max, the maximum measurement variance in Z (e.g., from Table 2),

• p_min, the minimum probability (i.e., confidence) that the measurement variance in Z is below υ_max.

For a given spatial field and sampling plan, p(υ ≤ υ_max) is the likelihood that the resulting measurement variance (υ) will be less than the specified maximum measurement variance (υ_max). The measurement variance requirements are fulfilled when p(υ ≤ υ_max) exceeds the probability threshold p_min, and the corresponding outcome is denoted ϕ_Z, analogous to that for ϵ and ϕ. Formally:

$\begin{array}{ll}p(\upsilon \le {\upsilon }_{\text{max}})\ge p_{\text{min}}\Rightarrow {\varphi }_Z,p(\upsilon \le {\upsilon }_{\text{max}})<p_{\text{min}}\Rightarrow {\overline{\varphi }}_Z.& (4)\end{array}$

2.3.2 Sampling model and parameter set

To compute p(υ ≤ υ_max), the probability of sufficiently small measurement variance, for different combinations of spatial field and sampling plan, we first define a sampling model with parameter space Θ_Z that encompasses stochastic simulation of heterogeneous spatial fields and composite sampling plans.

A field's spatial heterogeneity can be characterized by its spatial covariance, or strength of correlation between values at different locations depending on the physical distances separating them, often analytically represented by a covariance or semivariogram function (Webster and Oliver, 2007). These functions involve distribution parameters, here μ and CV expressed as natural logarithms, and a correlation length, λ, which describes how distant two locations can be and still have correlated values, or the “size” of the heterogeneities (Figure 7). Different analytical forms (e.g., exponential, circular, Gaussian) are distinguished by the “smoothness” of the heterogeneities (Figure 7). We consider correlation lengths λ that ranges 0%–50% of the width of the grid (e.g., a 1-ha square grid has dimensions 100 m by 100 m, so λ ranges 0–50 m).

Figure 7

Figure 7. Examples of simulated spatial fields with increasing correlation lengths and different analytical covariance models; example scales shown for different relative variances CV_Z given mean concentration 1,100 mg/kg (ln(μ_Z) = 7).

The parameter space Θ_Z encompasses these spatial field parameters, as well as parametrization of a composite sampling plan, including the number of composite samples (n) and sub-samples (n_sub), radius of each composite sample (r_c) expressed as % of the width of the grid, and margin of error intrinsic to the positioning device (e_pos) also expressed as % of grid width (Table 3). Given this setup, point sampling would be simulated by specifying 1 sub-sample (n_sub = 1) and a composite radius of 0 m (r_c = 0%).

Table 3

Table 3. Exploratory Θ_Z, uniform distributions used in stochastic spatial sampling simulations and sensitivity analysis to determine combinations of spatial field and sampling plan that result in measurement variances less than a target maximum.

2.4 Designing a sufficient sampling plan (Step 4)

Even if a deployment has favorable chemistry and operational choices, it can only meet the uncertainty requirements if its sampling design can reliably keep measurement variance below υ_max. Step 4 therefore moves from understanding how sampling behaves (Step 3) to identifying which sampling plans are actually sufficient. This is where the statistical requirements from Steps 1–3 inform concrete operational decisions: How many samples are needed? How large should each composite area be? And under what spatial conditions can any design succeed?

Because spatial heterogeneity governs whether a sampling design can achieve υ_max, Step 4 evaluates how heterogeneity interacts with sampling parameters and determines which combinations yield a high probability of satisfying the variance requirement. It then incorporates deployment-specific constraints to narrow these possibilities to sampling plans that are both sufficient and operationally realistic. In short, Step 4 answers the practical question: given the field's spatial structure and the project's constraints, what is the minimum sampling effort needed for the measurement approach to succeed?

2.4.1 Quantifying influence of spatial heterogeneity and sampling parameters on fulfilling measurement variance requirements (Step 4.1)

The ability of any sampling plan to meet the variance requirement depends not only on how many samples are collected, but on the underlying spatial heterogeneity of the field. A design that works under weak heterogeneity may fail completely when variability is structured or patchy. Step 4.1 quantifies this dependence by evaluating how combinations of spatial properties and sampling parameters jointly affect the probability of achieving υ ≤ υ_max.

Using the spatial fields and sampling model defined in Step 3, we compute p(υ ≤ υ_max) for thousands of realizations of heterogeneity and sampling configurations. This allows us to partition the parameter space into conditions where the variance requirement is fulfilled (ϕ_Z) and where it is not (${\overline{\varphi }}_Z$), and then identify which aspects of heterogeneity and design most strongly determine success. This step reveals whether sufficient sampling is even possible under the field conditions a project might encounter.

To partition the parameter space Θ_Z into Θ_Z|ϕ_Z and ${\Theta }_Z|{\overline{\varphi }}_Z$ for sensitivity analysis, we use nested Monte Carlo simulations to compute p(υ ≤ υ_max) for 10⁴ realizations of Θ_Z. After first generating the spatial field, we choose random locations for the n composite samples. For a single configuration of random locations, we simulate 30 rounds of paired composite sampling, computing the mean Z each round and the relative variance υ of the 30 means. Considering a quantification agent would only sample a handful of times throughout the course of a deployment, this represents 30 paired sampling rounds for a given configuration and the theoretical variability introduced by random positioning error and inconsistent orientation of sub-samples over a heterogeneous field. In the context of solid-phase EW verification, this assumes calculation of f_d uses the mean of all n samples rather than each sample individually. Altogether, for a single parameter realization of Θ_Z, we simulate 100 different configurations of random locations, and p(υ ≤ υ_max) is the portion of configurations where the inferred measurement variance υ is less than υ_max.

Results of these nested simulations indicate that ϕ_Z is determined by CV_Z, the relative field-scale variance of Z (Figure 8). Most realizations show an extremely low or extremely high likelihood of achieving a sufficiently small measurement variance υ (Figure 8A), and there is a clear CV_Z threshold around 3% (ln(CV_Z)≈−3.5) that dictates this behavior (Figure 8B). In practice, quantification agents will already have or can readily maintain an estimate of CV_Z, which will leave n as the most influential parameter and indicate the minimum number of samples needed to achieve ϕ_Z. Overall, this highlights that spatial heterogeneity not only needs to be accurately constrained before designing a sampling plan, but may also determine whether any monitoring strategy can succeed.

Figure 8

Figure 8. Exploratory sensitivity analysis of the influence of spatial heterogeneity and sampling plan on p(υ ≤ 0.01), the probability (p) that the inferred measurement variance (υ) is less than a maximum measurement variance (υ_max) of 1%, using the parameter ranges in Table 3. (A) Shows the response across all 10⁴ realizations, with a black dashed line separating the realizations that do and do not give at least 90% (p_min) likelihood of measurement variance less than 1%. (B) Provides the conditional response distribution for each parameter; parameter values where gray shaded p₁₀-p₉₀ regions are entirely above the black dashed p_min line indicate where at least 90% of realizations fulfill uncertainty requirements; black p₅₀ lines that extend above the black dashed p_min line indicate parameter values where at least half of the realizations fulfill the uncertainty requirements. (C) Provides the ranking of the parameters according to their influence on ϕ_Z, the event where uncertainty requirements are fulfilled, using the distance-based metric illustrated in Figure 4.

2.4.2 Refining to a specific sampling plan using deployment-specific constraints (Step 4.2)

After identifying which parameters matter most for meeting the variance requirement, the next step is to determine what sampling plan is actually sufficient and feasible for a given deployment. Real projects operate under constraints: expected ranges of field-scale variance, limits on sampling effort, or prior knowledge about heterogeneity derived from preliminary measurements. Step 4.2 incorporates these constraints to refine the broad sensitivity results of Step 4.1 into specific, actionable sampling plans.

To narrow down to a specific sampling plan using preliminary empirical data, we first constrain relative field-scale variance CV_Z, the major control on ϕ_Z determined through the sensitivity analysis in Step 4.1 (Figure 8C). Spatially distributed point sampling is typically necessary to capture the true CV, and such data are sparse for soil elemental composition at ha-scales. Field studies that have employed very high-density point sampling (10²-10⁴ samples/ha) of soil elemental concentrations show ha-scale variances of 5% ≤ CV ≤ 44% for base cations (Ca, Mg, Na, K) and select tracer elements (Ti, Ni) in semi-arid grassland, scrubland, forested, and agricultural settings (Table 4). For a quantification agent interested in constraining site-specific variance, further stochastic point-sampling simulations indicate that, given observed ranges (Table 4), at least 20 point samples are needed to estimate CV to the nearest ln with 90% confidence. This suggests it would be feasible to collect the preliminary measurements needed to infer operational scalability for a robust array of potential empirical constraints.

Table 4

Table 4. Plot-scale spatial variances of cation and tracer concentrations for four land-use types in semi-arid climates; sampling density ranges 220-14,500 samples/ha (Bahri et al., 1993; Gallardo, 2003; Gallardo and Paramá, 2007).

While at the lower end of field-scale variances we could find reported in literature, we perform a constrained sensitivity analysis assuming CV_Z = 10% to demonstrate next steps in sampling design for a theoretical deployment (Figure 9). Given this constraint, ϕ_Z becomes most sensitive to n (Figure 9C), indicating at least 130 composite samples will result in >90% likelihood, 90 samples >50% likelihood, that υ will be sufficiently small (Figure 9B). It may also be important to constrain λ, which shows secondary interactions with sample design parameters (Figure 10), potentially requiring alternative technologies like remote sensing (Chadwick and Asner, 2018) to control costs.

Figure 9

Figure 9. Deployment-specific sensitivity analysis where, relative to the exploratory sensitivity analysis in Figure 8 and parameter ranges in Table 3, we constrain relative field-scale variance to 10% (ln(CV_Z) = −2.3) and consider up to 200 samples. (A) Shows the updated response across all 10⁴ parameter realizations, (B) the updated conditional response distributions, and (C) the updated ranking of parameters according to their influence on ϕ_Z.

Figure 10

Figure 10. Secondary parameter interactions relative to their influence on p(υ ≤ υ_max), the likelihood of sufficiently small measurement variance, using second-order Sobol indices (Saltelli et al., 2004; Scheidt and Caers, 2009; Fenwick et al., 2014).

This example considers a maximum measurement variance υ_max of 1% and relative field-scale variance CV_Z of 10%. Simulating minimum n for three other combinations of υ_max and CV_Z shows that n quickly increases with smaller υ_max (Table 5), generally indicating infeasible n for υ_max < 1%.

Table 5

Table 5. Number of samples n required to satisfy uncertainty requirements for theoretical deployment given field-scale variances CV_Z of 5% and 10% and maximum measurement variances υ_max of 1% and 5%; these values assume point to composite sampling with up to five sub-samples and a constrained parameter space as described in Step 4.2.

Altogether, the analysis here indicates that spatial heterogeneity in soil concentrations should be the foremost consideration when designing sampling plans for solid-phase verification of EW. Determining sufficient sampling plans requires preliminary constraints on relevant field-scale variances, and even minimally sufficient plans determined a priori may be operationally infeasible, pointing toward reconsideration of the overall uncertainty requirements or measurement model.

2.5 Reporting the final estimate and uncertainty (Step 5)

In principle, reporting $\widehat{f_{\text{d}}}$ and its uncertainty is relatively straightforward following sample collection and analysis, as we have predetermined the margin of error and confidence in the calculation. However, transparency, reproducibility, and traceability are critical for both scientific and compensatory applications. In particular, procedural compliance should be separated from the inherent scientific and environmental uncertainty in evaluating whether a project achieved a result within quantifiable confidence bounds. Thus, in addition to providing the underlying measurement data and appropriate metadata, ideally under an IGSN (International GeoSample Number) framework (Damerow et al., 2021), the reporting framework should systematically capture key uncertainty targets and measurement distributions as constrained with each sampling event. Supplementary material 3 provides an example reporting format, where all sources of uncertainty are documented, including deviations from initial estimated field variance, and the final estimate is presented with clearly defined error and confidence.

3 Discussion

The methodology we evaluate here shows how constructing a measurement model that is directly informed by uncertainty requirements can help ensure monitoring approaches will deliver defensible CDR. This is in contrast to the way EW field studies are traditionally designed, where the statistics are largely handled ex post. The probabilistic framework is sequential in nature and propagates input variability through a measurement model, highlighting the critical parameters that influence the overall uncertainty (Steps 1–2; Figure 4B). For EW, once application rate and expected dissolution range are constrained, the measurement variances of base cation and tracer concentrations set the uncertainty in $\widehat{f_{\text{d}}}$, indicating variance thresholds above which the uncertainty requirements will not be met. Stochastic spatial simulations (Steps 3–4) identify composite sampling approaches that adhere to the measurement variance thresholds, and empirically constrained sensitivity analysis (Figure 9) points to a minimally sufficient plan. The suggested reporting format (Step 5; Supplementary material 3) documents key constraints and outcomes of this approach and allows for iterative update as more data become available.

3.1 Model assumptions and limitations

The framework as presented relies on several assumptions that could potentially influence simulated outcomes (e.g, Tables 2, 5):

• Spatial variance is the main contributor to measurement variance, though analytical uncertainty is an additional factor for some low-abundance chemical tracers. To account for this, future implementations can subtract analytical variance from the maximum measurement variances identified in Step 2, resulting in a lower target measurement variance in the sampling simulations.

• Heterogeneity with depth can be implicitly captured in a two-dimensional representation of a sampling volume, which may not be true if stratigraphic variations with depth are widely inconsistent over a field.

• No covariance between cation and tracer concentrations, no anisotropy (i.e., spatial correlation length is equal in all directions), and a stationary mean (i.e., Z has spatially uniform mean at the scale of simulation).

• Uniform feedstock dissolution (i.e., same f_d in all 10⁶ grid cells) and constant bulk density (e.g., negligible error due to mass loss upon dissolution, negligible density change with tillage).

3.2 Importance of spatial heterogeneity and correlation length

Our results demonstrate that without thorough assessment of spatial variance, even well-composited samples may not adequately constrain measurement uncertainty, potentially resulting in intensive data collection that fails to achieve the confidence needed to assign CDR. To evaluate the scale of variability expected for agricultural and grassland soils, we reviewed studies involving very high-density point sampling (10²–10⁴ samples/ha) of various land-use types in semi-arid climates (Table 4). The relative variances in soil elemental abundances all exceed 5%, with some up to 44%, suggesting that baseline soil heterogeneity may be too large for a mixing-model approach to feasibly produce reliable dissolution estimates (Table 5). While alternative solid-phase models, such as isotope mixing or bulk cation stocks, may be less sensitive to spatial heterogeneity, the analysis here poses an important consideration for EW, namely whether soil property distributions, and our ability to capture them with measurements, will fall within the requirements needed for verification frameworks.

Correlation length λ shows not only primary influence on achieving ϕ_Z (Figure 9), but also secondary interactions with operational parameters n, n_sub, composite sampling radius r_c, and positioning error e_pos (Figure 10). n_sub is essentially a multiplier of n, but not always an “efficient” multiplier since sub-samples are correlated when r_c < λ. Sensitivity analysis shows that λ near zero is best (Figure 9B), presumably because the sub-samples are independent of each other for r_c>λ, maximizing the amount of information gleamed from the sub-samples. Positioning error (e_pos up to 10% of grid width) also shows secondary interactions with λ, potentially related to whether e_pos is large relative to λ and would thus include paired sampling schemes where e_pos might shift a re-sampling location such that it is spatially uncorrelated with the preceding sample at that location, increasing the measurement variance; r_c small relative to λ means re-sampling will occur over a spatially correlated space for a given location (purple shaded areas in Figure 1). Altogether, this shows how characterization of correlation length can provide statistical justification for composite sampling design.

3.3 Future work

Spatially coordinated field-scale solid-phase geochemical data is sparse in existing literature. Samples are usually collected in vertical configuration at single disperse sites, rather than laterally across the surface in a spatially explicit configuration, making it difficult to use the wealth of existing data to assess spatial variability. While the studies evaluated here show consistent variances over four land types within a geographical region (Table 4), high-density soil data from other climates and geomorphic settings would be valuable in determining where field-scale geochemical variances are low enough for accurate solid-phase verification of EW. Spatial distributions of certain elements may also be correlated to textural parameters, such as clay content, or reflected in aboveground biomass (Chadwick and Asner, 2018). Rigorous confirmation and ground-truthing of potential proxies in early deployments could enable more efficient heterogeneity estimates in later deployments and thus more precise baseline sampling strategies at scale.

A related challenge is the identification of intensively measured plots, simulated here at 0.01–100 ha, but can be smaller or larger, that are adequate to test research hypotheses or represent project areas spanning tens-of-thousands to millions of hectares. Quantifying the degree to which small-scale trials capture relevant spatial variability in larger landscapes remains a critical open question.

4 Conclusion

Collectively, simulating spatial variability and realistic sampling strategies can reduce logistical inefficiencies for EW and minimize the risk of failing to meet uncertainty requirements for CDR quantification. Effective verification of open-system CDR requires balancing standardized sampling guidelines with the flexibility to accommodate diverse field conditions, as well as managing tradeoffs between measurement costs and the certainty of outcomes. The proposed framework addresses these challenges by clearly defining uncertainty targets, allowing sampling plans to be rigorously evaluated using preliminary empirical constraints prior to intensive field deployment. This transparent, multi-stage approach helps differentiate procedural compliance from the inherent scientific and environmental uncertainties affecting removal performance.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://github.com/rogersdb/ew-spatial-uq.

Author contributions

BR: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. KM: Conceptualization, Funding acquisition, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. The authors gratefully acknowledge support from the Stanford Sustainability Accelerator (GHG-0012). BR acknowledges support from the Department of Energy Computational Science Graduate Fellowship (DE-SC0021110) and Stanford Data Science.

Acknowledgments

The authors appreciate the time taken by two reviewers to help make this manuscript more relevant and clear.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fclim.2025.1688361/full#supplementary-material

References

Almaraz, M., Bingham, N. L., Goertzen, H., Holzer, I. O., Geoghegan, E. K., Sohng, J., et al. (2022). Methods for determining the CO₂ removal capacity of enhanced weathering in agronomic settings. Front. Clim. 4:970429. doi: 10.3389/fclim.2022.970429

Crossref Full Text | Google Scholar

Bahri, A., Berndtsson, R., and Jinno, K. (1993). Spatial dependence of geochemical elements in a semiarid agricultural field: I. Scale properties. Soil Sci. Soc. Am. J. 57, 1316–1322. doi: 10.2136/sssaj1993.03615995005700050026x

Crossref Full Text | Google Scholar

Beerling, D. J., Epihov, D. Z., Kantola, I. B., Masters, M. D., Reershemius, T., Planavsky, N. J., et al. (2024). Enhanced weathering in the US Corn Belt delivers carbon removal with agronomic benefits. Proc. Natl. Acad. Sci. USA. 121:e2319436121. doi: 10.1073/pnas.2319436121

PubMed Abstract | Crossref Full Text | Google Scholar

Beerling, D. J., Kantzas, E. P., Lomas, M. R., Wade, P., Eufrasio, R. M., Renforth, P., et al. (2020). Potential for large-scale CO₂ removal via enhanced rock weathering with croplands. Nature 583, 242–248. doi: 10.1038/s41586-020-2448-9

PubMed Abstract | Crossref Full Text | Google Scholar

Bertagni, M. B., and Porporato, A. (2022). The Carbon-capture efficiency of natural water alkalinization: implications for enhanced weathering. Sci. Total Environ. 838(Pt 4):156524. doi: 10.1016/j.scitotenv.2022.156524

PubMed Abstract | Crossref Full Text | Google Scholar

Bradford, M. A., Eash, L., Kuebbing, S. E., Polussa, A., Jevon, F. V., Hammac, W. A., et al. (2023). Testing the feasibility of quantifying change in agricultural soil carbon stocks through empirical sampling. Geoderma 440:116719. doi: 10.1016/j.geoderma.2023.116719

Crossref Full Text | Google Scholar

Brimhall, G. H., and Dietrich, W. E. (1987). Constitutive mass balance relations between chemical composition, volume, density, porosity, and strain in metasomatic hydrochemical systems: results on weathering and pedogenesis. Geochim. Cosmochim. Acta 51, 567–587. doi: 10.1016/0016-7037(87)90070-6

Crossref Full Text | Google Scholar

Calabrese, S., Wild, B., Bertagni, M. B., Bourg, I. C., White, C., Aburto, F., et al. (2022). Nano- to global-scale uncertainties in terrestrial enhanced weathering. Environ. Sci. Technol. 56, 15261–15272. doi: 10.1021/acs.est.2c03163

PubMed Abstract | Crossref Full Text | Google Scholar

Chadwick, K. D., and Asner, G. P. (2018). Landscape evolution and nutrient rejuvenation reflected in Amazon forest canopy chemistry. Ecol. Lett. 21, 978–988. doi: 10.1111/ele.12963

PubMed Abstract | Crossref Full Text | Google Scholar

Christakos, G., and Olea, R. A. (1992). Sampling design for spatially distributed hydrogeologic and environmental processes. Adv. Water Resour. 15, 219–237. doi: 10.1016/0309-1708(92)90008-P

Crossref Full Text | Google Scholar

Clarkson, M. O., Larkin, C. S., Swoboda, P., Reershemius, T., Suhrhoff, T. J., Maesano, C. N., et al. (2024). A review of measurement for quantification of carbon dioxide removal by enhanced weathering in soil. Front. Clim. 6:1345224. doi: 10.3389/fclim.2024.1345224

Crossref Full Text | Google Scholar

Damerow, J. E., Varadharajan, C., Brodie, E. L., Burrus, M., Chadwick, K. D., Crystal-Ornelas, R., et al. (2021). Sample identifiers and metadata to support data management and reuse in multidisciplinary ecosystem sciences. Data Sci. J. 20, 1–19. doi: 10.5334/dsj-2021-011

Crossref Full Text | Google Scholar

Demougeot-Renard, H., De Fouquet, C., and Renard, P. (2004). Forecasting the number of soil samples required to reduce remediation cost uncertainty. J. Environ. Qual. 33, 1694–1702. doi: 10.2134/jeq2004.1694

PubMed Abstract | Crossref Full Text | Google Scholar

Derry, L. A., Chadwick, O. A., and Porder, S. (2025). estimation of carbon dioxide removal via enhanced weathering. Glob. Change Biol. 31:e70067. doi: 10.1111/gcb.70067

PubMed Abstract | Crossref Full Text | Google Scholar

Dietzen, C., and Rosing, M. T. (2023). Quantification of CO₂ uptake by enhanced weathering of silicate minerals applied to acidic soils. Int. J. Greenhouse Gas Control 125:103872. doi: 10.1016/j.ijggc.2023.103872

Crossref Full Text | Google Scholar

Dupla, X., Claustre, R., Bonvin, E., Graf, I., Le Bayon, R. C., and Grand, S. (2024). Let the dust settle: impact of enhanced rock weathering on soil biological, physical, and geochemical fertility. Sci. Total Environ. 954:176297. doi: 10.1016/j.scitotenv.2024.176297

PubMed Abstract | Crossref Full Text | Google Scholar

Fenwick, D., Scheidt, C., and Caers, J. (2014). Quantifying asymmetric parameter interactions in sensitivity analysis: application to reservoir modeling. Math. Geosci. 46, 493–511. doi: 10.1007/s11004-014-9530-5

Crossref Full Text | Google Scholar

Gallardo, A. (2003). Spatial variability of soil properties in a floodplain forest in Northwest Spain. Ecosystems 6, 564–576. doi: 10.1007/s10021-003-0198-9

Crossref Full Text | Google Scholar

Gallardo, A., and Paramá, R. (2007). Spatial variability of soil elements in two plant communities of NW Spain. Geoderma 139, 199–208. doi: 10.1016/j.geoderma.2007.01.022

Crossref Full Text | Google Scholar

Guo, F., Wang, Y., Zhu, H., Zhang, C., Sun, H., Fang, Z., et al. (2023). Crop productivity and soil inorganic carbon change mediated by enhanced rock weathering in farmland: a comparative field analysis of multi-agroclimatic regions in central China. Agric. Syst. 210:103691. doi: 10.1016/j.agsy.2023.103691

Crossref Full Text | Google Scholar

Hartmann, J., West, A. J., Renforth, P., Köhler, P., Wolf-Gladrow, D. A., Rocha, C. L. D. L., et al. (2013). Enhanced chemical weathering as a geoengineering strategy to reduce atmospheric carbon dioxide, supply nutrients, and mitigate ocean acidification. Rev. Geophys. 51, 113–149. doi: 10.1002/rog.20004

Crossref Full Text | Google Scholar

Hasemer, H., Borevitz, J., and Buss, W. (2024). Measuring enhanced weathering: inorganic carbon-based approaches may be required to complement cation-based approaches. Front. Clim. 6:1352825. doi: 10.3389/fclim.2024.1352825

Crossref Full Text | Google Scholar

Holzer, I. O., Nocco, M. A., and Houlton, B. Z. (2023). Direct evidence for atmospheric carbon dioxide removal via enhanced weathering in cropland soil. Environ. Res. Commun. 5:101004. doi: 10.1088/2515-7620/acfd89

Crossref Full Text | Google Scholar

IPCC (2018). Global Warming of 1.5°C: IPCC Special Report on Impacts of Global Warming of 1.5°C above Pre-industrial Levels in Context of Strengthening Response to Climate Change, Sustainable Development, and Efforts to Eradicate Poverty, 1 Edn. Cambridge: Cambridge University Press.

Google Scholar

Kantola, I. B., Blanc-Betes, E., Masters, M. D., Chang, E., Marklein, A., Moore, C. E., et al. (2023). Improved net carbon budgets in the US Midwest through direct measured impacts of enhanced weathering. Glob. Change Biol. 29, 7012–7028. doi: 10.1111/gcb.16903

PubMed Abstract | Crossref Full Text | Google Scholar

Knapp, W. J., Stevenson, E. I., Renforth, P., Ascough, P. L., Knight, A. C. G., Bridgestock, L., et al. (2023). Quantifying CO₂ removal at enhanced weathering sites: a multiproxy approach. Environ. Sci. Technol. 57, 9854–9864. doi: 10.1021/acs.est.3c03757

PubMed Abstract | Crossref Full Text | Google Scholar

Lark, R. M. (2009). Estimating the regional mean status and change of soil properties: two distinct objectives for soil survey. Eur. J. Soil Sci. 60, 748–756. doi: 10.1111/j.1365-2389.2009.01156.x

Crossref Full Text | Google Scholar

Larkin, C. S., Andrews, M. G., Bellamy, J., Goring-Harford, H., Sadekar, S., James, R. H., et al. (2022). Quantification of CO₂ removal in a large-scale enhanced weathering field trial on an oil palm plantation in Sabah, Malaysia. Front. Clim. 4:959229. doi: 10.3389/fclim.2022.959229

Crossref Full Text | Google Scholar

Lebling, K., Riedl, D., and Leslie-Bole, H. (2024). Measurement, Reporting, and Verification for Novel Carbon Dioxide Removal in US Federal Policy. Available online at: https://www.wri.org/research/measurementreporting-and-verification-novel-carbon-dioxide-removal?auHash=rLtKBvDkX0auKYcYp_RiaIDTIxpQAIy2kLH21DTMNZs (Accessed April 18, 2025).

Google Scholar

Levy, C. R., Almaraz, M., Beerling, D. J., Raymond, P., Reinhard, C. T., Suhrhoff, T. J., et al. (2024). Enhanced rock weathering for carbon removal-monitoring and mitigating potential environmental impacts on agricultural land. Environ. Sci. Technol. 58, 17215–17226. doi: 10.1021/acs.est.4c02368

PubMed Abstract | Crossref Full Text | Google Scholar

Lewis, A. L., Sarkar, B., Wade, P., Kemp, S. J., Hodson, M. E., Taylor, L. L., et al. (2021). Effects of mineralogy, chemistry and physical properties of basalts on carbon capture potential and plant-nutrient element release via enhanced weathering. Appl. Geochem. 132:105023. doi: 10.1016/j.apgeochem.2021.105023

Crossref Full Text | Google Scholar

Lin, H. (2010). Earth's critical zone and hydropedology: concepts, characteristics, and advances. Hydrol. Earth Syst. Sci. 14, 25–45. doi: 10.5194/hessd-6-3417-2009

Crossref Full Text | Google Scholar

Ma, H. W., and Chang, C. C. (2008). Assessment of the value of reducing uncertainty by sampling in a groundwater remediation system. Sci. Total Environ. 402, 9–17. doi: 10.1016/j.scitotenv.2008.04.034

PubMed Abstract | Crossref Full Text | Google Scholar

McDermott, F., Bryson, M., Magee, R., and Van Acken, D. (2024). Enhanced weathering for CO₂ removal using carbonate-rich crushed returned concrete; a pilot study from SE Ireland. Appl. Geochem. 169:106056. doi: 10.1016/j.apgeochem.2024.106056

Crossref Full Text | Google Scholar

Mercer, L., and Burke, J. (2023). Strengthening MRV Standards for Greenhouse Gas Removals to Improve Climate Change Governance. London: Grantham Research Institute on Climate Change and the Environment and Centre for Climate Change Economics and Policy; London School of Economics and Political Science.

Google Scholar

NASEM (2019). Negative Emissions Technologies and Reliable Sequestration: A Research Agenda. Washington, DC: National Academies Press, 25259.

Google Scholar

Potash, E., Bradford, M. A., Oldfield, E. E., and Guan, K. (2025). Measure-and-remeasure as an economically feasible approach to crediting soil organic carbon at scale. Environ. Res. Lett. 20:024025. doi: 10.1088/1748-9326/ada16c

Crossref Full Text | Google Scholar

puro.earth (2024). Enhanced Rock Weathering Methodology for CO₂ Removal. Helsinki.

Google Scholar

Reershemius, T., Kelland, M. E., Jordan, J. S., Davis, I. R., D'Ascanio, R., Kalderon-Asael, B., et al. (2023). Initial validation of a soil-based mass-balance approach for empirical monitoring of enhanced rock weathering rates. Environ. Sci. Technol. 57, 19497–19507. doi: 10.1021/acs.est.3c03609

PubMed Abstract | Crossref Full Text | Google Scholar

Reershemius, T., and Suhrhoff, T. J. (2024). On error, uncertainty, and assumptions in calculating carbon dioxide removal rates by enhanced rock weathering in Kantola et al., 2023. Glob. Change Biol. 30:e17025. doi: 10.1111/gcb.17025

PubMed Abstract | Crossref Full Text | Google Scholar

Renforth, P. (2012). The potential of enhanced weathering in the UK. Int. J. Greenhouse Gas Control 10, 229–243. doi: 10.1016/j.ijggc.2012.06.011

Crossref Full Text | Google Scholar

Rieder, L., Amann, T., and Hartmann, J. (2023). Soil electrical conductivity as a proxy for enhanced weathering in soils. Front. Clim. 5:1283107 doi: 10.3389/fclim.2023.1283107

Crossref Full Text | Google Scholar

Saltelli, A., Tarantola, S., Campolongo, F., and Ratto, M. (2004). Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models. Hoboken, NJ: John Wiley and Sons, Ltd.

Google Scholar

Scheidt, C., and Caers, J. (2009). Representing spatial uncertainty using distances and kernels. Math. Geosci. 41, 397–419. doi: 10.1007/s11004-008-9186-0

Crossref Full Text | Google Scholar

Skov, K., Wardman, J., Healey, M., McBride, A., Bierowiec, T., Cooper, J., et al. (2024). Initial agronomic benefits of enhanced weathering using basalt: a study of spring oat in a temperate climate. PloS ONE 1:e0295031. doi: 10.1371/journal.pone.0295031

PubMed Abstract | Crossref Full Text | Google Scholar

Smith, D., Solano, F., Woodruff, L., Cannon, W., and Ellefsen, K. (2019). “Geochemical and mineralogical maps, with interpretation, for soils of the conterminous United States,” in USGS Scientific Investigations Report (Reston, VA: U.S. Geological Survey), 2017–5118, doi: 10.3133/sir20175118

Crossref Full Text | Google Scholar

Sokol, N. W., Moreland, K., Slessarev, E., Pett-Ridge, J., Sohng, J., Goertzen, H., et al. (2024). Reduced accrual of mineral-associated organic matter after two years of enhanced rock weathering in cropland soils, though no net losses of soil organic carbon. Biogeochemistry 167, 989–1005. doi: 10.1007/s10533-024-01160-0

Crossref Full Text | Google Scholar

Spear, R. (1980). Eutrophication in peel inlet—II. Identification of critical uncertainties via generalized sensitivity analysis. Water Res. 14, 43–49. doi: 10.1016/0043-1354(80)90040-8

Crossref Full Text | Google Scholar

Suhrhoff, T. J., Planavsky, N. J., Reershemius, T., Wang, J., Jordan, J. S., and Reinhard, C. T. (2024). A tool for assessing the sensitivity of soil-based approaches for quantifying enhanced weathering: a US case study. Front. Clim. 6:1346117. doi: 10.3389/fclim.2024.1346117

Crossref Full Text | Google Scholar

Sutherland, K., Holme, E., Savage, R., Gill, S., Matlin-Wainer, M., He, J., et al. (2025). Enhanced Weathering in Agriculture. Available online at: https://registry.isometric.com/protocol/enhanced- weathering-agriculture#determination-2 (Accessed August 19, 2025).

Google Scholar

Vereecken, H., Schnepf, A., Hopmans, J. W., Javaux, M., Or, D., Roose, T., et al. (2016). Modeling soil processes: review, key challenges, and new perspectives. Vadose Zone J. 15, 1–57. doi: 10.2136/vzj2015.09.0131

Crossref Full Text | Google Scholar

Verma, P. S. (1998). Error propagation in geochemical modeling of trace elements in two-component mixing. Geofis. Int. 37, 327–338. doi: 10.22201/igeof.00167169p.1998.37.4.517

Crossref Full Text | Google Scholar

Verra (2012). VCS Module VMD0021: Estimation of Stocks in the Soil Carbon Pool. Available online at: https://verra.org/methodologies/vm0021-soil-carbon-quantification-methodology-v1-0/ (Accessed August 19, 2025).

Google Scholar

Verstraete, S., and Van Meirvenne, M. (2008). A multi-stage sampling strategy for the delineation of soil pollution in a contaminated brownfield. Environ. Pollut. 154, 184–191. doi: 10.1016/j.envpol.2007.10.014

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, F., Zhu, F., Liu, D., Qu, Y., Liu, D., Xie, J., et al. (2024). Wollastonite powder application increases rice yield and CO₂ sequestration in a paddy field in Northeast China. Plant Soil. 502, 589–603. doi: 10.1007/s11104-024-06570-5

Crossref Full Text | Google Scholar

Wasserstein, R. L., Schirm, A. L., and Lazar, N. A. (2019). Moving to a World Beyond “p < 0.05”. Am. Statist. 73, 1–19. doi: 10.1080/00031305.2019.1583913

Crossref Full Text | Google Scholar

Webster, R., and Oliver, M. A. (2008). Geostatistics for Environmental Scientists, 2nd Edn. Hoboken, NJ: John Wiley and Sons, Ltd, 1–315. doi: 10.1002/9780470517277

Crossref Full Text | Google Scholar