Spatial Pipeline EngineeringMRV Automation

Soil Organic Carbon Modeling & Validation

Soil organic carbon (SOC) modeling is the hardest measurement problem in the spatial modeling and carbon stock validation stack, and the one where a plausible-looking pipeline most easily produces indefensible numbers. Above-ground biomass can be observed directly by satellite; soil carbon cannot. It is inferred from covariates — terrain, climate, land cover, spectral indices — calibrated against physical cores, and the signal being detected is a change of one to four tonnes of carbon per hectare per year against a stock that varies by tens of tonnes across a single field. The measurement noise is larger than the effect. Every architectural decision in an SOC pipeline is therefore a decision about how to keep that ratio honest, and the failure modes here are statistical rather than computational.

Digital soil mapping pipeline from covariates and cores to a validated stock-change estimate Two input streams converge. Environmental covariates — terrain derivatives, climate normals, land-cover history, and bare-soil spectral composites — join laboratory-measured soil cores carrying organic carbon concentration, bulk density, and coarse-fragment content. A model-fitting stage uses spatially blocked cross-validation rather than random splits. The fitted model predicts concentration, which is converted to stock through bulk density and coarse-fragment correction and then to an equivalent-soil-mass basis. A validation gate compares predictions against a held-out probability sample; only after passing does the pipeline emit a stock-change estimate with its uncertainty. A rejection path returns models that pass random cross-validation but fail spatial cross-validation. Covariates terrain · climate · land cover bare-soil spectral composite Soil cores SOC % · bulk density coarse fragments · depth Model fit spatially blocked CV never a random split Concentration → stock × BD × (1 − coarse frac) → equivalent soil mass Validation gate held-out probability sample · design-based Stock change ± interval tCO₂e ha⁻¹ yr⁻¹ Rejected: passes random CV, fails spatial CV the classic optimism trap — refit or reduce model complexity

Role in the MRV Workflow

SOC modeling occupies the same slot in the pipeline as biomass estimation — it converts observations into a stock per hectare that a baseline comparison turns into a credit — but it consumes a fundamentally different evidence base. Where biomass estimation from LiDAR/SAR fusion is calibrated by field plots and validated against structure the sensor genuinely sees, SOC prediction is calibrated by destructive sampling and validated against more destructive sampling. There is no independent remote observation of the quantity. The satellite covariates carry information about the drivers of soil carbon, not the carbon itself, which puts a hard ceiling on achievable accuracy and makes the validation design the single most important artefact the project produces.

Upstream, the stage depends on a covariate stack that must be temporally coherent with the sampling campaign. A bare-soil composite built from acquisitions spanning three years, joined to cores collected in one season, embeds a mismatch that no model can undo. It also depends on laboratory data whose method is recorded: dry combustion and wet oxidation do not produce interchangeable numbers, and a project that mixes methods across periods will observe a step change in stock that is entirely analytical. Downstream, the stock estimate flows into the same reconciliation and reporting path as any other pool, subject to the emissions data quality validation gates and the schema contract in the MRV data schema reference.

The distinguishing architectural requirement is that SOC pipelines must carry a design-based validation alongside the model. Model-based prediction gives you the map; a probability sample, drawn independently and held out entirely from fitting, gives you an unbiased estimate of the map’s error with a confidence interval that does not depend on the model being right. Most methodologies now require the second, and for good reason: a model can be arbitrarily wrong in ways its own cross-validation will not reveal if the validation data share the model’s blind spots.

Core Failure Modes

  1. Random cross-validation on spatially clustered samples. Soil cores are almost never a spatially independent sample — they come in transects, in field-corner clusters, in whatever the field team could reach. Split them randomly and a point’s near neighbour, forty metres away and sharing every covariate value, ends up in the other fold. The model memorises location rather than process, and reported R² comes back at 0.8 when honest spatial performance is 0.3. This is the single most common defect in published SOC maps and it inflates every downstream claim. The remedy is spatially blocked or leave-location-out cross-validation, with block size set from the range of the covariate semivariogram, and the honest number reported even when it is disappointing.

  2. Fixed-depth stock accounting under changing bulk density. SOC stock is concentration times bulk density times depth, corrected for coarse fragments. Management that increases carbon usually also decreases bulk density — the soil becomes less dense as organic matter accumulates. Compute stock over a fixed 0–30 cm depth in both periods and you are weighing a different mass of soil each time, and the resulting apparent gain can be almost entirely an artefact of the density change. Equivalent soil mass accounting fixes the mass rather than the depth and removes the artefact. In tilled-to-no-till conversions the difference between fixed-depth and equivalent-mass accounting routinely runs 20–40% of the claimed change, in the direction that flatters the project.

  3. Detecting a change smaller than the measurement noise. Field-scale SOC spatial variability commonly has a coefficient of variation of 20–40%, while the annual change being claimed is 1–3% of stock. Without a sampling design sized for the effect — paired sampling at fixed georeferenced locations, enough replicates, and a long enough interval — the confidence interval on stock change straddles zero, and a project that reports a point estimate without it is reporting noise. The fix is to size the design before sampling, using a power calculation from local variance estimates, and to accept a longer re-measurement interval where power is insufficient.

Random versus spatially blocked cross-validation on clustered soil cores Two side-by-side sampling maps of the same field. On the left, random cross-validation: sample clusters are split so that neighbouring points fall in different folds, illustrated by alternating fill within each cluster, and the reported R squared is 0.81 with a root mean square error of 4.2 tonnes of carbon per hectare. On the right, spatially blocked cross-validation: whole clusters are assigned to a single fold, shown by uniform fill per block, and the honest performance drops to an R squared of 0.34 with a root mean square error of 9.8. A caption states that the difference is not a worse model but a truthful measurement of it. The same model, measured two ways Clustered cores make random folds leak location into the test set. Random CV — optimistic Spatially blocked CV — honest reported R² 0.81 RMSE 4.2 t C ha⁻¹ actual R² 0.34 RMSE 9.8 t C ha⁻¹ Filled and hollow markers are test and training points. Blocking does not make the model worse — it stops the score from borrowing a neighbour's answer.

Deterministic Implementation Architecture

The implementation below fits a quantile-regression-forest SOC model with spatially blocked cross-validation, converts predicted concentration to an equivalent-soil-mass stock, and refuses to emit a stock-change figure whose confidence interval includes zero without labelling it as non-significant. Telemetry records the block size, the fold assignment, and both the optimistic and honest scores so the gap is never invisible.

from dataclasses import dataclass, asdict

import geopandas as gpd
import numpy as np
import structlog
from prefect import flow, task
from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GroupKFold

log = structlog.get_logger()

CANONICAL_CRS = "EPSG:6933"      # equal-area, so per-hectare stocks are comparable
C_TO_CO2E = 44.0 / 12.0


@dataclass(frozen=True)
class SocResult:
    period: str
    stock_tco2e_ha: float
    stock_ci_low: float
    stock_ci_high: float
    change_tco2e_ha_yr: float | None
    change_significant: bool | None
    r2_random_cv: float
    r2_blocked_cv: float
    block_size_m: float
    depth_basis: str
    lab_method: str
    n_cores: int


@task
def spatial_blocks(cores: gpd.GeoDataFrame, block_size_m: float) -> np.ndarray:
    """Assign each core to a spatial block. Folds are drawn over BLOCKS, so a
    point's near neighbour can never sit in the opposite fold."""
    if cores.crs is None:
        raise ValueError("core locations carry no CRS; spatial blocking would be meaningless")
    xy = np.c_[cores.to_crs(CANONICAL_CRS).geometry.x, cores.to_crs(CANONICAL_CRS).geometry.y]
    n_blocks = max(4, int(np.ptp(xy[:, 0]) * np.ptp(xy[:, 1]) / block_size_m ** 2))
    groups = KMeans(n_clusters=min(n_blocks, len(cores) // 3), n_init=10,
                    random_state=0).fit_predict(xy)
    log.info("soc.blocks.assigned", blocks=int(groups.max() + 1),
             block_size_m=block_size_m, cores=len(cores))
    return groups


@task
def fit_and_score(X: np.ndarray, y: np.ndarray, groups: np.ndarray) -> tuple[object, float, float]:
    """Fit once, score twice. The random score is kept only so the optimism gap
    is visible in the log — it is never the number that goes in the report."""
    from sklearn.model_selection import KFold, cross_val_score

    model = RandomForestRegressor(n_estimators=500, min_samples_leaf=4, random_state=0)

    r2_random = float(np.mean(cross_val_score(
        model, X, y, cv=KFold(5, shuffle=True, random_state=0), scoring="r2")))
    r2_blocked = float(np.mean(cross_val_score(
        model, X, y, cv=GroupKFold(n_splits=min(5, len(set(groups)))),
        groups=groups, scoring="r2")))

    if r2_random - r2_blocked > 0.25:
        log.warning("soc.cv.optimism_gap", r2_random=round(r2_random, 3),
                    r2_blocked=round(r2_blocked, 3),
                    note="model is learning location, not process")

    model.fit(X, y)
    return model, r2_random, r2_blocked


@task
def equivalent_soil_mass_stock(
    soc_pct: np.ndarray, bulk_density: np.ndarray, coarse_frac: np.ndarray,
    depth_cm: float, reference_mass_kg_m2: float,
) -> tuple[np.ndarray, str]:
    """Stock on an equivalent-soil-mass basis.

    Fixed-depth accounting compares a different MASS of soil between periods
    whenever bulk density changes — which is exactly what carbon accrual causes.
    Reporting to a fixed reference mass removes that artefact.
    """
    fine_fraction = 1.0 - coarse_frac
    mass_kg_m2 = bulk_density * 1000.0 * (depth_cm / 100.0) * fine_fraction

    # Scale each profile to the common reference mass rather than a common depth.
    scale = reference_mass_kg_m2 / mass_kg_m2
    stock_t_c_ha = (soc_pct / 100.0) * reference_mass_kg_m2 * 10.0

    log.info("soc.stock.esm", reference_mass_kg_m2=reference_mass_kg_m2,
             mean_scale=round(float(np.mean(scale)), 3),
             mean_stock_t_c_ha=round(float(np.mean(stock_t_c_ha)), 2))
    return stock_t_c_ha * C_TO_CO2E, f"equivalent-soil-mass@{reference_mass_kg_m2:.0f}kg/m2"


@flow(name="soc_stock_estimate")
def run(
    period: str, cores_path: str, covariates: np.ndarray, block_size_m: float,
    depth_cm: float, reference_mass_kg_m2: float, lab_method: str,
    prior_stock_tco2e_ha: float | None = None, years_elapsed: float | None = None,
) -> dict:
    cores = gpd.read_file(cores_path)
    for column in ("soc_pct", "bulk_density_g_cm3", "coarse_frac", "lab_method"):
        if column not in cores.columns:
            raise ValueError(f"core table missing required column: {column}")
    if cores["lab_method"].nunique() > 1:
        # Dry combustion and wet oxidation are not interchangeable; mixing them
        # produces an analytical step change that looks like sequestration.
        raise ValueError("cores mix laboratory methods; harmonise or model separately")

    groups = spatial_blocks(cores, block_size_m)
    y = cores["soc_pct"].to_numpy()
    model, r2_random, r2_blocked = fit_and_score(covariates, y, groups)

    predicted = model.predict(covariates)
    stock, basis = equivalent_soil_mass_stock(
        predicted, cores["bulk_density_g_cm3"].to_numpy(),
        cores["coarse_frac"].to_numpy(), depth_cm, reference_mass_kg_m2)

    mean_stock = float(np.mean(stock))
    se = float(np.std(stock, ddof=1) / np.sqrt(len(stock)))
    ci_low, ci_high = mean_stock - 1.96 * se, mean_stock + 1.96 * se

    change = significant = None
    if prior_stock_tco2e_ha is not None and years_elapsed:
        change = (mean_stock - prior_stock_tco2e_ha) / years_elapsed
        # Significance against the interval, not against the point estimate.
        significant = not (ci_low <= prior_stock_tco2e_ha <= ci_high)
        if not significant:
            log.warning("soc.change.not_significant", change=round(change, 3),
                        ci=(round(ci_low, 2), round(ci_high, 2)),
                        note="report as non-significant; do not credit")

    result = SocResult(
        period=period, stock_tco2e_ha=round(mean_stock, 2),
        stock_ci_low=round(ci_low, 2), stock_ci_high=round(ci_high, 2),
        change_tco2e_ha_yr=None if change is None else round(change, 3),
        change_significant=significant, r2_random_cv=round(r2_random, 3),
        r2_blocked_cv=round(r2_blocked, 3), block_size_m=block_size_m,
        depth_basis=basis, lab_method=lab_method, n_cores=len(cores),
    )
    log.info("soc.run.complete", **asdict(result))
    return asdict(result)

The three guards worth defending: blocked scoring is the reported score, with the random score retained only so the optimism gap appears in the log; equivalent soil mass, not fixed depth, so a density change cannot be sold as carbon; and significance tested against the interval, so a change smaller than the noise is labelled rather than credited.

Validation, Debugging & Compliance Mapping

  • Design-based validation → the number auditors trust. Fit the model however you like, but estimate its error against an independent probability sample that took no part in fitting. The resulting design-based mean and confidence interval are model-independent and are what most SOC methodologies now require. Report both: the model-based map for spatial allocation, the design-based estimate for the reported total.
  • Equivalent soil mass → comparability across periods. Record the reference mass and the depth from which each profile was scaled. Verifiers under the major soil methodologies check this explicitly, because fixed-depth accounting under changing bulk density is the best-known way to overstate soil carbon gains.
  • Laboratory method provenance → no analytical step changes. Record the method, the laboratory, and the reference-material results for every batch, and route them into MRV data lineage and provenance tracking. A change of laboratory between periods is a legitimate operational event and an illegitimate source of stock change; only the lineage record can tell them apart.
  • Uncertainty → conservative deduction. Where the confidence interval is wide, most methodologies apply a conservativeness deduction scaled to the interval. Propagate the model, sampling, and bulk-density uncertainties jointly rather than reporting the model error alone — the approach set out under emission factor uncertainty mapping applies directly.

The equivalent-soil-mass correction is easiest to trust once you have seen the arithmetic laid out side by side on the same field.

Fixed-depth versus equivalent-soil-mass accounting on the same field after a no-till conversion Two soil-column diagrams for the same field. The baseline column, tilled, has a bulk density of 1.42 grams per cubic centimetre and 1.35 percent soil organic carbon over 30 centimetres, giving 5751 kilograms of soil per square metre and 77.6 tonnes of carbon dioxide equivalent per hectare. The period-five column, no-till, has a lower bulk density of 1.28 and 1.44 percent carbon. Measured to the same fixed 30 centimetre depth it appears to hold 88.2 tonnes, an apparent gain of 10.6. Corrected to the same soil mass, the comparable depth is 33.3 centimetres and the true stock is 81.9 tonnes, a real gain of 4.3. The difference, 6.3 tonnes or 59 percent of the apparent gain, is the density artefact. Where 59% of an apparent soil-carbon gain came from Same field, tilled baseline versus no-till at period five. Only the accounting basis differs. Baseline · tilled Period 5 · fixed 30 cm Period 5 · equal soil mass BD 1.42 g cm⁻³ SOC 1.35% 5751 kg m⁻² soil 30 cm 77.6 tCO₂e ha⁻¹ BD 1.28 g cm⁻³ SOC 1.44% 5184 kg m⁻² — less soil 30 cm 88.2 apparent gain +10.6 BD 1.28 g cm⁻³ SOC 1.44% 5751 kg m⁻² — matched 33.3 cm 81.9 real gain +4.3 Density artefact 6.3 tCO₂e ha⁻¹ 59% of the apparent gain creditable only under ESM

For debugging, three checks find most problems. Compare the random and blocked cross-validation scores; a gap wider than about 0.25 in R² means the model is learning location. Plot predicted against observed for the held-out probability sample and look for the classic regression-to-the-mean fan, which signals that the covariates carry less information than the model’s confidence implies. And recompute the stock change on a fixed-depth basis alongside the equivalent-mass basis: if they disagree by more than roughly 10%, bulk density is moving and the equivalent-mass figure is the only defensible one.

Frequently Asked Questions

Why can’t I just use random cross-validation like every other machine-learning problem?

Because soil cores are spatially autocorrelated and clustered, which breaks the independence assumption that makes a random split meaningful. A held-out point forty metres from a training point shares its terrain, climate, and land-cover values almost exactly, so predicting it is closer to lookup than to generalisation. The result is a score that describes interpolation between neighbours rather than prediction at new locations — which is what the map is actually used for. Spatially blocked cross-validation restores the correspondence between the score and the use.

How deep should I sample, and does it matter for crediting?

It matters a great deal. Most cropland methodologies credit 0–30 cm because that is where management effects are largest and sampling is affordable, but carbon redistributes with depth under tillage change, so a shallow-only design can record a gain that is partly a transfer from below. Where the methodology allows it, sample to 100 cm on a subset to characterise the depth distribution even if you credit only the top layer, and report the deeper results as supporting evidence. Whatever depth you choose, keep it identical across periods and record it on every row.

What is equivalent soil mass and why is fixed depth not good enough?

Fixed-depth accounting compares the carbon in the top 30 cm at time one against the top 30 cm at time two. If bulk density fell between them — which is what happens when organic matter accumulates or tillage stops — the second sample contains less soil mass than the first, so you are comparing unequal things. Equivalent-soil-mass accounting fixes the mass of fine-earth soil and reports the carbon it contains, removing the artefact. The correction typically moves a claimed change by 20–40% in no-till conversions.

How many cores do I need to detect a real change?

Size it from a power calculation, not a rule of thumb. You need the local spatial variance of SOC stock, the effect size you expect per year, the re-measurement interval, and the acceptable error rates. Field-scale coefficients of variation of 20–40% against an annual effect of 1–3% of stock usually mean either a large paired sample at fixed georeferenced points or a re-measurement interval of five years or more. Paired sampling at the same locations is dramatically more powerful than independent samples because it removes the spatial variance from the comparison.

Should I model SOC concentration or SOC stock directly?

Model concentration, then convert. Concentration is what the laboratory measures and what the covariates plausibly predict; stock is a derived quantity that folds in bulk density and coarse-fragment content, each with its own error structure and each often measured on a smaller subset of cores than concentration. Modelling stock directly hides those separate error sources inside one residual and makes the equivalent-soil-mass correction impossible to apply afterwards. Keep the three components separate through the pipeline, propagate their uncertainties jointly at the end, and you retain the ability to answer the question a verifier eventually asks: how much of this interval comes from the carbon measurement and how much from the density estimate?

Can satellite data measure soil carbon directly?

Not through vegetation, and not below the surface. Bare-soil spectral composites carry a genuine signal for surface organic matter on exposed soil, and they are useful covariates, but they observe the top few millimetres under specific conditions — dry, bare, low residue — and say nothing directly about carbon at 30 cm. Treat spectral data as one covariate among terrain, climate, and management history, and keep the physical cores as the source of truth. A pipeline that presents a satellite-derived SOC map without a core-based validation is presenting a hypothesis.

Conclusion

Soil carbon is where MRV engineering meets a genuinely hard statistical problem, and the honest response is architectural rather than algorithmic. Block your cross-validation spatially and report the honest score. Account on an equivalent-soil-mass basis so a density change cannot be mistaken for sequestration. Hold out an independent probability sample so the reported error does not depend on the model being right. And size the sampling design for the effect before the field team leaves, because no amount of modelling recovers power that was never in the data. Continue with modeling soil organic carbon with digital soil mapping and validating soil carbon models against core samples.