Monte Carlo Uncertainty Propagation for Emission Factors #

A total-emissions figure is only as defensible as the uncertainty that ships with it. This guide is the task-level recipe under Emission Factor Uncertainty Mapping, the variance-propagation stage within the Spatial Modeling & Carbon Stock Validation framework. It consumes the calibrated factor distributions produced upstream and the per-pixel biomass error bands emitted by biomass estimation from LiDAR & SAR fusion, then converts them into a single number an auditor can reconstruct: a mean total CO2e, a 95% confidence interval, and the conservative discount a registry deducts before crediting.

The engineering problem is that emission factors and activity data each carry their own error, those errors are frequently correlated, and the analytic error-propagation formula most spreadsheets use assumes they are not. When a project multiplies a factor by an activity quantity and sums thousands of such products across a landscape, the naive first-order propagation understates the joint variance wherever factors share a source table or activity strata share a survey instrument. Monte Carlo simulation sidesteps the closed-form algebra entirely: it draws correlated samples from the input distributions, recomputes the total emissions on every draw, and reads the confidence interval straight off the empirical distribution of results. The output is not a smoother number — it is a distribution the ground-truth alignment for carbon models stage can interrogate and a verifier can rerun.

Root Cause Analysis: Why Analytic Propagation Understates the Interval #

The default tool for combining uncertainties is the first-order Taylor expansion — the “add the squares of the relative errors” rule embedded in most carbon spreadsheets and in the IPCC Approach 1 guidance. It is fast, closed-form, and correct only under assumptions that spatial emissions inventories routinely violate. Three of those violations dominate production failures.

First, the formula assumes the inputs are independent, and factors are not. When a hundred parcels draw their soil-carbon factor from the same regional lookup table, an error in that table moves every parcel in the same direction at once. Independent errors partially cancel under summation; shared errors accumulate. Ignoring a factor correlation of 0.6 across a large aggregation can understate the standard error of the total by 40–70%, manufacturing a tight interval that collapses the moment a verifier asks how the factors were sourced. Monte Carlo handles this natively because the correlation lives in the sampling step, not in the algebra.

Second, the formula assumes symmetric, near-Gaussian inputs, and many factors are not. Emission factors bounded below by zero, or derived from ratio estimators, are frequently lognormal or otherwise skewed. A symmetric error bar around a skewed factor puts the reported mean in the wrong place and mis-sizes both tails, which matters precisely because registries care about the conservative tail. Sampling from the declared distribution family — lognormal, triangular, truncated normal — preserves that shape all the way to the total.

Third, spatial autocorrelation in activity data breaks the independence of the summands themselves. Neighbouring cells share land-use classifiers, survey instruments, and acquisition geometry, so their activity errors covary even when the factors do not. Treating each cell’s contribution as an independent random variable — the implicit assumption when you sum analytic variances cell by cell — understates the aggregate variance for the same reason it does with correlated factors. A correlated draw over the activity field, or an explicit spatial correlation length, restores the honest interval.

Diagnostic Pipeline: Pre-Flight Uncertainty Metadata Validation #

Monte Carlo will happily simulate garbage. Before any draw runs, the inputs must be inspected: every factor and activity term needs a mean, a dispersion (standard deviation or a confidence interval it can be derived from), a named distribution family, and — where relevant — a correlation structure that is actually usable. The gate below rejects the failure conditions above rather than letting them produce a plausible-looking but indefensible interval. Strict validation here is the same discipline the parent emission factor uncertainty mapping stage applies to its rasters.

import numpy as np
import structlog

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

VALID_FAMILIES = {"normal", "lognormal", "truncated_normal", "triangular"}
MAX_CV = 1.0          # a CV above 1.0 signals a mis-parsed dispersion, not real spread
MIN_CV = 1e-4         # a near-zero CV usually means a missing uncertainty column


def validate_uncertainty_inputs(terms: list[dict], corr: np.ndarray) -> None:
    """Gate emission-factor / activity uncertainty metadata before simulation.

    Each term must declare: name, mean, sd (or a CI to derive it), distribution.
    `corr` is the n x n correlation matrix over the terms in order.
    Raises on any condition that would make the Monte Carlo draw meaningless.
    """
    n = len(terms)
    for i, t in enumerate(terms):
        for key in ("name", "mean", "sd", "distribution"):
            if key not in t or t[key] is None:
                raise ValueError(f"term[{i}] missing required field '{key}'")
        if t["distribution"] not in VALID_FAMILIES:
            raise ValueError(f"{t['name']}: unknown distribution '{t['distribution']}'")
        if t["mean"] <= 0 and t["distribution"] == "lognormal":
            raise ValueError(f"{t['name']}: lognormal requires a positive mean")
        cv = t["sd"] / t["mean"] if t["mean"] else np.inf
        if not (MIN_CV <= cv <= MAX_CV):
            raise ValueError(
                f"{t['name']}: coefficient of variation {cv:.3f} outside "
                f"[{MIN_CV}, {MAX_CV}] — check the dispersion column")

    # Correlation matrix must be square, symmetric, unit-diagonal, and PSD;
    # an invalid matrix silently distorts the joint draw.
    if corr.shape != (n, n):
        raise ValueError(f"correlation matrix {corr.shape} != ({n}, {n})")
    if not np.allclose(corr, corr.T, atol=1e-8):
        raise ValueError("correlation matrix is not symmetric")
    if not np.allclose(np.diag(corr), 1.0, atol=1e-8):
        raise ValueError("correlation matrix diagonal must be 1.0")
    eig_min = float(np.linalg.eigvalsh(corr).min())
    if eig_min < -1e-8:
        raise ValueError(
            f"correlation matrix is not positive semi-definite (min eig {eig_min:.2e})")

    log.info("uncertainty_inputs_validated", n_terms=n, min_eigenvalue=round(eig_min, 6))

A run that raises here never reaches the simulator. The most common trigger in practice is a coefficient of variation above one, which nearly always means a raw variance was parsed into the standard deviation column, followed by a correlation matrix that lost positive semi-definiteness after a manual edit.

Deterministic Transformation Logic: Vectorized Correlated Monte Carlo #

The simulator draws correlated samples for every input term at once, vectorized across all 10,000 iterations, then computes the total emissions per iteration in a single matrix operation. Correlation is imposed with a Gaussian copula: draw standard-normal variates through the Cholesky factor of the correlation matrix, map each column to a uniform, then push it through the inverse CDF of that term’s declared distribution. This preserves both the marginal shapes and the cross-term dependence without assuming everything is jointly Gaussian. The run is seeded for reproducibility, and a convergence gate refuses to return an interval that has not stabilized.

import numpy as np
from scipy import stats
import structlog

log = structlog.get_logger()


def _sample_marginal(term: dict, u: np.ndarray) -> np.ndarray:
    """Map uniforms u in (0,1) to a term's declared distribution via its inverse CDF."""
    mean, sd, fam = term["mean"], term["sd"], term["distribution"]
    if fam == "normal":
        return stats.norm.ppf(u, loc=mean, scale=sd)
    if fam == "lognormal":
        # Convert arithmetic mean/sd to the log-space parameters.
        sigma = np.sqrt(np.log1p((sd / mean) ** 2))
        mu = np.log(mean) - 0.5 * sigma ** 2
        return stats.lognorm.ppf(u, s=sigma, scale=np.exp(mu))
    if fam == "truncated_normal":
        a = (0.0 - mean) / sd  # truncate at zero — areas and factors are non-negative
        return stats.truncnorm.ppf(u, a=a, b=np.inf, loc=mean, scale=sd)
    if fam == "triangular":
        half = sd * np.sqrt(6.0)  # symmetric triangle whose sd matches the declared sd
        lower, upper = mean - half, mean + half
        return stats.triang.ppf(u, c=0.5, loc=lower, scale=upper - lower)
    raise ValueError(f"unsupported distribution: {fam}")


def propagate_emissions_mc(
    factor_terms: list[dict],
    activity_terms: list[dict],
    corr: np.ndarray,
    n_iter: int = 10_000,
    conservative_percentile: float = 2.5,
    seed: int = 20260714,
    mc_se_tol: float = 0.005,
) -> dict:
    """Monte Carlo propagation of factor and activity uncertainty to total CO2e.

    factor_terms[i] pairs positionally with activity_terms[i]: emissions for
    source i on one iteration are factor_i * activity_i, summed over sources.
    `corr` is the correlation matrix over the interleaved [factors, activities].
    Returns mean, sd, 95% CI, a convergence flag, and the conservative estimate.
    """
    rng = np.random.default_rng(seed)
    terms = factor_terms + activity_terms
    n = len(terms)
    if corr.shape != (n, n):
        raise ValueError("correlation matrix must cover all factor+activity terms")

    # Gaussian copula: correlated standard normals -> uniforms -> marginals.
    chol = np.linalg.cholesky(corr)                      # raises if not PSD
    z = rng.standard_normal(size=(n_iter, n)) @ chol.T   # (n_iter, n) correlated normals
    u = stats.norm.cdf(z)                                # uniforms preserving rank corr

    samples = np.empty((n_iter, n), dtype="float64")
    for j, term in enumerate(terms):
        samples[:, j] = _sample_marginal(term, u[:, j])

    n_src = len(factor_terms)
    factors = samples[:, :n_src]                         # (n_iter, n_src)
    activities = samples[:, n_src:]                      # (n_iter, n_src)
    totals = (factors * activities).sum(axis=1)          # total CO2e per iteration

    mean = float(totals.mean())
    sd = float(totals.std(ddof=1))
    p_lo, p_hi = np.percentile(totals, [2.5, 97.5])
    conservative = float(np.percentile(totals, conservative_percentile))

    # Convergence gate: the standard error of the mean must be small relative
    # to the mean, or the interval is still noisy and must not be reported.
    mc_se = sd / np.sqrt(n_iter)
    converged = bool((mc_se / mean) < mc_se_tol) if mean else False
    if not converged:
        log.warning("mc_not_converged", relative_mc_se=round(mc_se / mean, 5),
                    tolerance=mc_se_tol, n_iter=n_iter)

    result = {
        "mean_tco2e": round(mean, 2),
        "sd_tco2e": round(sd, 2),
        "ci95_lower": round(float(p_lo), 2),
        "ci95_upper": round(float(p_hi), 2),
        "conservative_tco2e": round(conservative, 2),
        "conservative_percentile": conservative_percentile,
        "n_iter": n_iter,
        "seed": seed,
        "converged": converged,
        "relative_mc_se": round(mc_se / mean, 5) if mean else None,
    }
    log.info("mc_propagation_complete", **result)
    return result

The conservative estimate is the load-bearing output. Rather than reporting the mean and subtracting a fixed percentage, the routine reads the 2.5th percentile of the simulated distribution directly, which is exactly the “lower bound of the 95% confidence interval” that Verra and IPCC conservativeness rules ask a project to credit against. Because the percentile comes from the correlated draw, it already reflects the widened interval that analytic propagation would have hidden — the deduction is neither an arbitrary haircut nor an over-tight number a verifier can reject.

Compliance Gating & Audit Trail Generation #

The simulation output becomes a submission artifact only when it is bound to the assumptions that produced it. Three gates make that binding auditable. The convergence flag must be true before any figure is quoted; a converged=False result is a draft, not a claim, and the logged relative Monte Carlo standard error tells a reviewer how far off it was. The conservative percentile is recorded explicitly so a verifier can confirm the project credited the lower bound rather than the mean — this maps directly onto the conservativeness principle in ISO 14064-3 and the uncertainty-deduction clauses of the Verra VM-series. The seed, sample count, and input covariance are serialized alongside the result so the entire interval can be regenerated bit-for-bit, which is what turns “we ran a Monte Carlo” into a reproducible control an auditor can rerun.

Map the outputs to the frameworks as follows. The 95% CI and the conservative estimate answer ISO 14064-3, which expects the reported total to be both reproducible and conservative; reading the credited figure off the correlated lower tail satisfies both at once. The relative interval width feeds Verra VM-series uncertainty deductions, where a methodology-specific tolerance decides how much of the mean a project may claim. The distribution shape and the declared input families support CSRD ESRS E1 disclosure, which scrutinizes climate figures for transparent treatment of estimation uncertainty rather than a single unqualified number. Because the correlation matrix travels with the result, the audit trail also records the modelling choice most likely to be challenged — how strongly the factors were assumed to move together — and lets a reviewer stress-test it. Those covariance assumptions belong in the versioned factor store described in versioning emission factor databases for reproducible MRV, so the exact factor revision behind a submitted figure is always recoverable.

Production Integration #

Deploy the simulator inside the MRV orchestration layer along a fixed ingest → diagnose → transform → validate → export → submit sequence:

1. Ingest. Pull the emission-factor distributions and their correlation matrix from the versioned factor database, and the activity-data quantities with their dispersions from the inventory tables, pinning the exact factor revision so the run is reproducible.
2. Diagnose. Run validate_uncertainty_inputs to confirm every term carries a mean, dispersion, distribution family, and a positive semi-definite correlation matrix; reject the batch on any breach rather than substituting a default.
3. Transform. Call propagate_emissions_mc with n_iter=10_000 and the project seed, drawing correlated factor and activity samples through the Gaussian copula and reducing to the total-emissions distribution.
4. Validate. Assert converged is True and that the relative interval width sits within the methodology tolerance; a non-converged or implausibly wide result is quarantined for review, never quoted.
5. Export. Serialize the mean, 95% CI, conservative estimate, seed, sample count, and correlation matrix into the run’s audit record so the figure can be regenerated on demand.
6. Submit. Forward the conservative estimate as the creditable volume and the full distribution as supporting evidence, feeding the baseline logic in forest carbon baseline & additionality modeling and the calibration checks in ground-truth alignment for carbon models.

Run the simulation as an idempotent task keyed on the input revision and seed, so a rerun over unchanged inputs reproduces the identical interval and a changed factor revision produces a new, traceable record. By validating uncertainty metadata before it simulates, imposing correlation in the draw rather than assuming it away, gating on convergence, and reading the conservative figure straight off the empirical tail, Monte Carlo propagation turns a fragile analytic estimate into a total-emissions figure that survives third-party verification and carries its own reproducible proof.

Related guides #

• Emission Factor Uncertainty Mapping — the parent variance stage this simulation implements for total emissions.
• Biomass Estimation from LiDAR & SAR Fusion — the upstream source of the per-pixel error bands that seed the factor distributions.
• Forest Carbon Baseline & Additionality Modeling — the downstream consumer that credits against the conservative estimate.
• Versioning Emission Factor Databases for Reproducible MRV — where the factor distributions and covariance assumptions are pinned per revision.