Forest Carbon Baseline & Additionality Modeling #

Forest Carbon Baseline & Additionality Modeling is the stage that decides how many tonnes a REDD+ or afforestation project may credibly claim, by constructing the counterfactual — what would have happened to the forest without the intervention — and subtracting it from the observed outcome. It is the most contested quantitative step in the Spatial Modeling & Carbon Stock Validation stack, because the counterfactual is unobservable by definition and every methodological choice moves the credited volume. This component consumes the calibrated stock surfaces and confidence envelopes produced by emission factor uncertainty mapping and hands its baseline surface to threshold tuning for carbon stock baselines, where the numeric cut-offs that separate forest from non-forest are calibrated against it.

The last five years have seen forest carbon baselining move decisively away from self-selected, static, project-level baselines towards dynamic performance benchmarks derived from matched control areas and jurisdictional deforestation rates — a shift crystallised by Verra’s consolidated REDD methodology VM0047. That transition is not a compliance footnote; it changes the GIS data a pipeline must ingest, the statistical machinery it must run, and the lineage it must retain. Where an older baseline was a single historical average frozen into a project design document, a modern baseline is a spatially explicit, periodically re-estimated surface that must survive the same scrutiny the emission factor uncertainty mapping stage applies to every other figure, and it must feed the same auditable threshold tuning logic downstream. This page covers how baselines and additionality are modelled defensibly: matched controls and synthetic controls, historical reference regions, deforestation-risk modelling, leakage belts, permanence buffers, and the statistical tests that make the credited difference something an auditor can reproduce from the pixels.

Role in the MRV Workflow #

Baseline and additionality modelling occupies the point in the carbon accounting pipeline where physical measurement becomes a financial claim. Everything upstream — cloud-masked reflectance, fused biomass, calibrated emission factors, propagated variance — produces the observed carbon stock inside the project boundary. This stage answers the only question a registry actually pays for: how much of that stock is genuinely additional, meaning it exists because of the intervention and would not have existed otherwise. Because the counterfactual can never be measured directly, it must be estimated, and the credibility of the estimate is entirely a function of how comparable the reference is to the project and how honestly the estimate is bounded and deducted.

The stage has hard upstream contracts. It requires a project boundary and a leakage belt as versioned, topology-valid vector layers; a stack of historical land-cover or forest-cover rasters spanning at least the reference period the methodology requires; and a set of covariate surfaces — elevation, slope, distance to roads and settlements, distance to prior forest edge, and baseline canopy cover — resampled to a common equal-area grid so that matching is done on like-for-like pixels. Downstream, it emits a spatially explicit baseline surface (expected forest loss or expected stock trajectory per pixel or per stratum), a matched-control diagnostic bundle, and a credited-difference figure already net of the deductions a verifier expects. That figure flows into the same append-only lineage layer every other MRV output uses, so a third party can trace an issued credit back to the specific control pixels and the specific counterfactual that justified it.

Three families of baseline construction dominate practice, and a production system should treat the choice as a data-and-defensibility trade-off rather than a preference. A historical-average baseline projects the project area’s own past deforestation rate forward; it is cheap and simple but self-referential and easy to game. A matched-control or synthetic-control baseline builds the counterfactual from areas that resemble the project on observable drivers of deforestation, which is far more defensible but demands rich covariate data and careful statistics. A jurisdictional or dynamic performance benchmark — the direction of travel under VM0047 — allocates an expected loss rate from a wider administrative region’s observed deforestation, updated periodically, which resists gaming but requires reliable jurisdiction-wide monitoring. The table below sets out how they compare on the axes that actually determine whether credits survive audit.

Modern pipelines rarely run one approach in isolation. A VM0047-aligned system typically combines a jurisdictional expected-loss allocation with a matched-control cross-check and an explicit uncertainty treatment, so that the dynamic benchmark constrains the aggregate while covariate matching confirms the project’s own risk profile is genuinely represented. Whichever combination is used, the baseline is no longer a number written once — it is a re-estimated surface with a defined validity window, and the pipeline must be built to regenerate it on schedule.

Core Failure Modes #

Three failure modes account for most of the credibility crisis that has engulfed forest carbon in recent years. Each has a specific root cause in the spatial and statistical construction of the baseline, and each has a measurable, well-documented impact on issued volume.

1. Non-comparable control area inflating additionality. The root cause is selecting a reference region — or synthetic control pixels — whose deforestation pressure is systematically higher than the project area truly faced, without demonstrating covariate balance. If the donor pool skews towards flatter, more accessible, closer-to-road land while the project sits on steeper or remoter terrain, the estimated counterfactual predicts forest loss the project would never have experienced, and the observed-minus-counterfactual difference is credited as avoided deforestation that was never at risk. This is not a marginal effect: independent reanalyses of REDD+ projects using rigorous synthetic controls found that a large share of projects had reference-region deforestation rates several times higher than matched counterfactuals justified, with headline baselines overstated by roughly two- to three-fold and many projects delivering additionality of only a small fraction of the credits issued. The failure is invisible in the project’s own reporting because the project never audits its control against itself — only an external matched-control reconstruction exposes it.

2. Ignoring leakage, causing over-crediting. The root cause is treating the project boundary as a closed system and crediting avoided emissions inside it while the deforestation activity simply relocates to unprotected forest outside it. Activity-shifting leakage (loggers and farmers move next door) and market leakage (reduced local supply raises prices and drives clearing elsewhere) both mean the atmosphere sees less benefit than the boundary suggests. If the pipeline does not monitor a leakage belt and deduct the loss observed there, the credited figure double-counts protection that did not occur at the landscape scale. Empirical leakage estimates for avoided-deforestation projects commonly fall in the range of roughly 20 to 80 percent of gross claimed reductions depending on the driver, so omitting the deduction can more than halve the true impact while leaving the issued volume untouched.

3. Static baseline drift as regional deforestation changes. The root cause is fixing a historical-average baseline at project design and holding it constant for a decade while the regional deforestation rate that the baseline was meant to represent moves underneath it. When regional clearing slows — because of enforcement, commodity-price shifts, or macroeconomic change — a frozen high baseline keeps crediting against a threat that has receded, and the gap between the assumed counterfactual and reality widens every year. Even a modest divergence of a few percentage points of forest-loss rate per year compounds: over a ten-year crediting period a static baseline can drift to double the genuine counterfactual, generating a steadily growing stream of non-additional credits. This is precisely the failure the shift to periodically re-estimated dynamic benchmarks under VM0047 is designed to eliminate, and it is the strongest engineering argument for building baselines as regenerable surfaces rather than constants in a config file.

Deterministic Implementation Architecture #

The implementation below constructs a matched-control counterfactual and computes additionality with an explicit uncertainty deduction. It selects donor pixels by covariate matching on elevation, slope, distance-to-road, and baseline forest cover; enforces an equal-area CRS so per-pixel areas are honest; gates match quality with a standardised-mean-difference balance test that raises before any credit is computed; estimates the counterfactual loss from the matched controls; and subtracts a conservative uncertainty deduction, a leakage-belt loss, and a permanence buffer contribution before emitting an issued figure. Structured telemetry records every gate so the run is reproducible.

from __future__ import annotations

from dataclasses import dataclass, asdict
from datetime import datetime, timezone

import numpy as np
import geopandas as gpd
import rioxarray  # noqa: F401  (registers the .rio accessor)
import xarray as xr
from scipy.spatial import cKDTree
import structlog

structlog.configure(
    processors=[
        structlog.processors.add_log_level,
        structlog.processors.TimeStamper(fmt="iso", utc=True),
        structlog.processors.JSONRenderer(),
    ]
)
log = structlog.get_logger()

EQUAL_AREA_CRS = "EPSG:6933"     # World Cylindrical Equal Area — area-true for pools
MAX_SMD = 0.10                   # max |standardised mean diff| per covariate after matching
CARBON_FRACTION = 0.47          # IPCC default C fraction of dry biomass
DEFAULT_UNCERTAINTY_DEDUCTION = 0.15   # conservative deduction on gross additionality
DEFAULT_BUFFER_RATE = 0.15             # non-permanence buffer-pool contribution


@dataclass
class AdditionalityResult:
    gross_avoided_tco2e: float
    uncertainty_deduction_tco2e: float
    leakage_deduction_tco2e: float
    buffer_contribution_tco2e: float
    issued_credits_tco2e: float
    n_matched_pairs: int
    max_covariate_smd: float


def _load_equal_area(path: str, name: str) -> xr.DataArray:
    """Open a raster, force the canonical equal-area CRS, refuse untagged input."""
    da = rioxarray.open_rasterio(path, masked=True, chunks={"x": 2048, "y": 2048})
    if da.rio.crs is None:
        raise ValueError(f"{name}: untagged CRS; refusing to guess a datum.")
    if da.rio.crs.to_epsg() != 6933:
        da = da.rio.reproject(EQUAL_AREA_CRS)
        log.info("reprojected_to_equal_area", layer=name, target=EQUAL_AREA_CRS)
    return da.squeeze(drop=True)


def _standardised_mean_diff(treat: np.ndarray, ctrl: np.ndarray) -> float:
    """Absolute standardised mean difference — the covariate-balance diagnostic."""
    pooled_sd = np.sqrt((treat.var(ddof=1) + ctrl.var(ddof=1)) / 2.0) or 1.0
    return float(abs(treat.mean() - ctrl.mean()) / pooled_sd)


def select_matched_controls(
    covariates: dict[str, xr.DataArray],
    project_mask: xr.DataArray,
    donor_mask: xr.DataArray,
    covariate_order: tuple[str, ...] = (
        "elevation", "slope", "distance_to_road", "baseline_cover"),
) -> tuple[np.ndarray, np.ndarray, float]:
    """Nearest-neighbour covariate matching of donor pixels to project pixels.

    Matching is done on z-scored covariates so no single driver (e.g. metres of
    distance-to-road) dominates the Euclidean distance. Returns project pixel
    indices, their matched donor indices, and the worst post-match SMD.
    """
    stack = np.stack([covariates[c].values.ravel() for c in covariate_order], axis=1)
    proj_sel = project_mask.values.ravel().astype(bool)
    donor_sel = donor_mask.values.ravel().astype(bool)
    finite = np.isfinite(stack).all(axis=1)

    mu, sd = np.nanmean(stack[finite], axis=0), np.nanstd(stack[finite], axis=0)
    sd = np.where(sd == 0, 1.0, sd)
    z = (stack - mu) / sd

    proj_idx = np.flatnonzero(proj_sel & finite)
    donor_idx = np.flatnonzero(donor_sel & finite)
    if proj_idx.size == 0 or donor_idx.size == 0:
        raise RuntimeError("empty project or donor pool after finiteness filtering.")

    tree = cKDTree(z[donor_idx])
    _, nn = tree.query(z[proj_idx], k=1)
    matched_donor = donor_idx[nn]

    worst_smd = max(
        _standardised_mean_diff(z[proj_idx, j], z[matched_donor, j])
        for j in range(z.shape[1])
    )
    log.info("covariate_matching_complete",
             n_project=int(proj_idx.size), n_matched=int(matched_donor.size),
             worst_smd=round(worst_smd, 4), covariates=list(covariate_order))
    return proj_idx, matched_donor, worst_smd


def estimate_additionality(
    forest_loss_start: xr.DataArray,
    forest_loss_end: xr.DataArray,
    covariates: dict[str, xr.DataArray],
    project_mask: xr.DataArray,
    donor_mask: xr.DataArray,
    leakage_belt_loss_tco2e: float,
    pixel_area_ha: float,
    biomass_density_t_ha: float,
    uncertainty_deduction: float = DEFAULT_UNCERTAINTY_DEDUCTION,
    buffer_rate: float = DEFAULT_BUFFER_RATE,
) -> AdditionalityResult:
    """Counterfactual avoided emissions from matched controls, net of deductions."""
    proj_idx, matched_donor, worst_smd = select_matched_controls(
        covariates, project_mask, donor_mask)

    # Hard balance gate: an unbalanced match manufactures additionality.
    if worst_smd > MAX_SMD:
        log.error("balance_gate_failed", worst_smd=round(worst_smd, 4), limit=MAX_SMD)
        raise RuntimeError(
            f"covariate balance SMD {worst_smd:.3f} exceeds {MAX_SMD}; "
            "control area is not comparable — additionality would be inflated.")

    # Counterfactual loss rate = observed loss in the matched controls.
    donor_loss = (forest_loss_end.values.ravel()[matched_donor]
                  - forest_loss_start.values.ravel()[matched_donor])
    proj_loss = (forest_loss_end.values.ravel()[proj_idx]
                 - forest_loss_start.values.ravel()[proj_idx])
    counterfactual_loss_px = np.clip(donor_loss, 0.0, None)
    observed_loss_px = np.clip(proj_loss, 0.0, None)

    # Avoided forest loss (ha) -> biomass -> CO2e.
    avoided_ha = float((counterfactual_loss_px - observed_loss_px).sum()) * pixel_area_ha
    avoided_ha = max(avoided_ha, 0.0)
    gross_tco2e = avoided_ha * biomass_density_t_ha * CARBON_FRACTION * (44.0 / 12.0)

    unc_ded = gross_tco2e * uncertainty_deduction
    net_after_unc = gross_tco2e - unc_ded - leakage_belt_loss_tco2e
    net_after_unc = max(net_after_unc, 0.0)
    buffer = net_after_unc * buffer_rate
    issued = net_after_unc - buffer

    log.info("additionality_estimated",
             gross_tco2e=round(gross_tco2e, 1),
             uncertainty_deduction_tco2e=round(unc_ded, 1),
             leakage_deduction_tco2e=round(leakage_belt_loss_tco2e, 1),
             buffer_tco2e=round(buffer, 1),
             issued_tco2e=round(issued, 1),
             worst_smd=round(worst_smd, 4))

    return AdditionalityResult(
        gross_avoided_tco2e=round(gross_tco2e, 1),
        uncertainty_deduction_tco2e=round(unc_ded, 1),
        leakage_deduction_tco2e=round(leakage_belt_loss_tco2e, 1),
        buffer_contribution_tco2e=round(buffer, 1),
        issued_credits_tco2e=round(issued, 1),
        n_matched_pairs=int(proj_idx.size),
        max_covariate_smd=round(worst_smd, 4),
    )


if __name__ == "__main__":
    cov = {name: _load_equal_area(f"data/{name}.tif", name)
           for name in ("elevation", "slope", "distance_to_road", "baseline_cover")}
    result = estimate_additionality(
        forest_loss_start=_load_equal_area("data/cover_2016.tif", "cover_2016"),
        forest_loss_end=_load_equal_area("data/cover_2026.tif", "cover_2026"),
        covariates=cov,
        project_mask=_load_equal_area("data/project_mask.tif", "project_mask"),
        donor_mask=_load_equal_area("data/donor_mask.tif", "donor_mask"),
        leakage_belt_loss_tco2e=42_000.0,
        pixel_area_ha=0.09,          # 30 m pixel on an equal-area grid
        biomass_density_t_ha=140.0,
    )
    log.info("baseline_run_sealed",
             result=asdict(result),
             generated_at=datetime.now(timezone.utc).isoformat(),
             method="matched_control_counterfactual_vm0047_aligned")

The design choices are deliberate. Matching is performed on z-scored covariates so that distance-to-road in metres does not swamp slope in degrees inside the Euclidean nearest-neighbour search; the standardised-mean-difference gate then refuses to proceed if any covariate remains unbalanced after matching, which is the single most effective guard against the non-comparable-control failure mode. The equal-area CRS makes the hectare arithmetic honest, and the deductions are applied in the order registries expect — uncertainty and leakage first, then the non-permanence buffer contribution on the net figure — so the issued number is already conservative when it leaves the function. In a full VM0047-aligned deployment the counterfactual loss rate would be reconciled against a jurisdictional allocation rather than taken from the donor pool alone, but the matched-control balance test remains the diagnostic that proves the project’s own risk profile is genuinely represented.

Validation, Debugging & Compliance Mapping #

Every output of this stage maps to a specific control in the standards that govern forest carbon, which is what turns the baseline surface from an analyst’s spreadsheet into a submission artifact. The covariate-balance SMD and the matched-control diagnostic bundle answer Verra VM0047, whose dynamic performance-benchmark logic requires that the expected forest loss attributed to a project be derived from comparable, monitored land rather than a self-selected historical rate; a logged, pre-registered matching procedure with demonstrated balance is precisely the evidence a VM0047 validation body examines, and the periodic re-estimation window the surface carries satisfies the methodology’s move away from static baselines. The leakage-belt deduction and the permanence buffer contribution map to the leakage and non-permanence provisions common to both Verra and Gold Standard land-use methodologies, which require that displaced emissions be quantified and that a share of credits be withheld against reversal risk before any issuance. The uncertainty deduction, and the fact that it is computed from the propagated variance handed over by the emission factor uncertainty mapping stage, satisfies the conservativeness principle of ISO 14064-2 and the estimation-uncertainty disclosure expectations of CSRD ESRS E1, which scrutinises land-sector removals for transparent treatment of how confident the reported figure actually is.

For debugging, treat the worst covariate SMD, the fraction of project pixels with no acceptable donor match, and the ratio of leakage-belt loss to gross avoided loss as monitored signals on every run — including the ones that pass — so that a slowly degrading donor pool or a widening jurisdictional divergence surfaces as a trend long before it breaches a hard gate. Three silent failures deserve dedicated diagnostics: a donor pool that has shrunk until nearest-neighbour matches are being drawn from ecologically distant pixels while the SMD still squeaks under the limit; a static baseline left un-regenerated past its validity window, which should raise on the run date rather than continue crediting; and a leakage belt whose own forest cover has been quietly clipped by an unrelated boundary edit, understating displacement. Validation should include a placebo test — running the matching against a pseudo-project inside the donor pool, which should yield near-zero additionality — and a sensitivity sweep on the uncertainty deduction and buffer rate against the methodology’s floors.

Conclusion #

Forest carbon baseline and additionality modelling is where a pipeline either earns or forfeits the trust of the market. The failures that have discredited so many REDD+ projects — non-comparable controls that inflated the counterfactual, ignored leakage that double-counted protection, and static baselines that kept crediting a receding threat — are all engineering failures, and all are addressable with defensible spatial statistics: covariate matching with a hard balance gate, an explicit leakage-belt deduction, a permanence buffer, an uncertainty deduction traced to a propagated envelope, and baselines rebuilt as regenerable surfaces aligned with the dynamic, jurisdictional direction of VM0047. Build the baseline this way and the credited difference becomes something an auditor can reconstruct from the pixels rather than take on faith. For a worked, task-level walkthrough of applying these techniques to an avoided-deforestation project, continue with Modeling Additionality Baselines for REDD+ Projects.

Related guides #

• Spatial Modeling & Carbon Stock Validation — the parent framework whose stock surfaces this baseline stage turns into additionality claims.
• Modeling Additionality Baselines for REDD+ Projects — the task-level guide implementing matched controls end to end.
• Emission Factor Uncertainty Mapping — the upstream envelope that supplies the uncertainty deduction applied here.
• Verra VM0047 vs Gold Standard GIS Requirements — the methodology comparison that governs which baseline approach a project may use.