Soil carbon modeling

This page reproduces the two statistical models in Bob Carpenter's soil-carbon case study. Both models describe carbon moving between two soil pools and leaving the system as evolved CO2. They differ only in whether the reported CO2 means are treated as direct observations of that process or as noisy summaries of latent sample-level values.

The short time series below is a build fixture rather than the case study's experimental data. The process equations, initial condition, priors, and two observation models are preserved. Open either comparison to see the complete generated Stan program next to the exact Julia source evaluated by this docs build; readers do not have to run the model to reveal the Stan.

The two-pool process model

Let and denote carbon in the two soil pools. Carbon leaves pool 1 at rate and pool 2 at rate , while fractions and feed those losses back into the other pool:

The mixing fraction divides the initial total carbon between the pools. At each observation time, cumulative evolved CO2 is the initial total minus the carbon still present:

In StanBlocks, the ODE right-hand side is a typed, Stan-only @deffun. ode_rk45 uses Stan's current variadic interface, so the physical parameters are passed directly instead of being packed into the legacy theta, x_r, and x_i arrays. The solver returns a time-by-state matrix; selecting its two columns and converting them with to_vector gives the process prediction.

The reusable soil_dynamics submodel owns both the positive rate/feedback parameters and the ODE solve. A parent model writes eCO2_hat ~ soil_dynamics(; ...): this is SLIC composition, not a probability statement. It inlines the submodel hygienically and binds its returned prediction to eCO2_hat.

Direct residual model

The first model compares the reported mean at each time directly with the ODE prediction. Its positive sigma consequently absorbs both measurement uncertainty and process/model discrepancy. Call-time keyword arguments become Stan data, while names first introduced on the left of ~ become parameters or transformed quantities. StanBlocks infers their shapes from the ODE result and likelihood, and emits the declared lower bounds into Stan.

using StanBlocks

@deffun @stanonly begin
    two_pool_feedback(
        t::real, carbon::vector[ny],
        k1::real, k2::real, alpha21::real, alpha12::real,
    )::vector[ny] = begin
        dcarbon::vector[ny]
        dcarbon[1] = -k1 * carbon[1] + alpha12 * k2 * carbon[2]
        dcarbon[2] = alpha21 * k1 * carbon[1] - k2 * carbon[2]
        dcarbon
    end
end

soil_data = (;
    totalC_t0=7.7,
    t0=0.0,
    ts=[1.0, 2.0, 3.0, 4.0],
    eCO2mean=[0.3, 0.8, 1.5, 2.1],
)
soil_measurement_data = Base.merge(
    soil_data,
    (; eCO2sd=[0.10, 0.10, 0.15, 0.20]),
)

soil_dynamics = @slic begin
    k1 ~ std_normal(; lower=0)
    k2 ~ std_normal(; lower=0)
    alpha21 ~ std_normal(; lower=0)
    alpha12 ~ std_normal(; lower=0)
    gamma ~ beta(10, 1)

    carbon0 = append_row(
        rep_vector(gamma * totalC_t0, 1),
        (1 - gamma) * totalC_t0,
    )
    trajectory = ode_rk45(
        two_pool_feedback, carbon0, t0, to_array_1d(ts),
        k1, k2, alpha21, alpha12,
    )
    return totalC_t0 -
        to_vector(trajectory[:, 1]) -
        to_vector(trajectory[:, 2])
end

soil_measurement = @slic begin
    eCO2_hat ~ soil_dynamics(; totalC_t0, t0, ts)
    sigma ~ cauchy(0, 1; lower=0)
    eCO2mean ~ normal(eCO2_hat, sigma)
end
soil_measurement_posterior = soil_measurement(; soil_data...)

functions {
vector two_pool_feedback(
    real t,
    vector carbon,
    real k1,
    real k2,
    real alpha21,
    real alpha12
) {
    int ny = dims(carbon)[1];
    vector[ny] dcarbon;
    dcarbon[1] = (((-k1) * carbon[1]) + (alpha12 * k2 * carbon[2]));
    dcarbon[2] = ((alpha21 * k1 * carbon[1]) - (k2 * carbon[2]));
    return dcarbon;
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    return to_vector(normal_rng(a, b));
}
}
data {
    real totalC_t0;
    int ts_n;
    real t0;
    vector[ts_n] ts;
    int eCO2mean_n;
    vector[eCO2mean_n] eCO2mean;
}
transformed data {
}
parameters {
    real<lower=0> eCO2_hat_k1;
    real<lower=0> eCO2_hat_k2;
    real<lower=0> eCO2_hat_alpha21;
    real<lower=0> eCO2_hat_alpha12;
    real<lower=0, upper=1> eCO2_hat_gamma;
    real<lower=0> sigma;
}
transformed parameters {
    vector[(1 + 1)] eCO2_hat_carbon0 = append_row(rep_vector((eCO2_hat_gamma * totalC_t0), 1), ((1 - eCO2_hat_gamma) * totalC_t0));
    array[ts_n] vector[(1 + 1)] eCO2_hat_trajectory = ode_rk45(
        two_pool_feedback,
        eCO2_hat_carbon0,
        t0,
        to_array_1d(ts),
        eCO2_hat_k1,
        eCO2_hat_k2,
        eCO2_hat_alpha21,
        eCO2_hat_alpha12
    );
    vector[ts_n] eCO2_hat = ((totalC_t0 - to_vector(eCO2_hat_trajectory[:, 1])) - to_vector(eCO2_hat_trajectory[:, 2]));
}
model {
    eCO2_hat_k1 ~ std_normal();
    eCO2_hat_k2 ~ std_normal();
    eCO2_hat_alpha21 ~ std_normal();
    eCO2_hat_alpha12 ~ std_normal();
    eCO2_hat_gamma ~ beta(10, 1);
    sigma ~ cauchy(0, 1);
    eCO2mean ~ normal(eCO2_hat, sigma);
}
generated quantities {
    vector[eCO2mean_n] eCO2mean_likelihood = normal_lpdfs(eCO2mean, eCO2_hat, sigma);
    vector[eCO2mean_n] eCO2mean_gen = normal_vector_rng(eCO2mean_n, eCO2_hat, sigma);
}

Separate measurement error

The extension keeps the process model unchanged and inserts a positive latent eCO2 value at every observation time. The process residual scale sigma describes variation around the ODE trajectory; the known eCO2sd values then describe uncertainty in the reported means. This is the distinction the direct model cannot make.

Base.merge constructs the second data set from the common process data plus the measurement standard deviations. The latent vector's length is tied explicitly to ts, and lower = 0 preserves the original model's physical support. Most importantly, the model reuses the same soil_dynamics submodel rather than copying the rates, priors, initial condition, and ODE. The generated-Stan view shows the full inlined program despite that source-level composition.

soil_latent = @slic begin
    eCO2_hat ~ soil_dynamics(; totalC_t0, t0, ts)
    sigma ~ cauchy(0, 1; lower=0)
    eCO2 ~ normal(eCO2_hat, sigma; n=length(ts), lower=0)
    eCO2mean ~ normal(eCO2, eCO2sd)
end
soil_latent_posterior = soil_latent(; soil_measurement_data...)

functions {
vector two_pool_feedback(
    real t,
    vector carbon,
    real k1,
    real k2,
    real alpha21,
    real alpha12
) {
    int ny = dims(carbon)[1];
    vector[ny] dcarbon;
    dcarbon[1] = (((-k1) * carbon[1]) + (alpha12 * k2 * carbon[2]));
    dcarbon[2] = ((alpha21 * k1 * carbon[1]) - (k2 * carbon[2]));
    return dcarbon;
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    vector scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    vector b
) {
    int n = anontok__1;
    return to_vector(normal_rng(a, b));
}
}
data {
    real totalC_t0;
    int ts_n;
    real t0;
    vector[ts_n] ts;
    int eCO2mean_n;
    vector[eCO2mean_n] eCO2mean;
    int eCO2sd_n;
    vector[eCO2sd_n] eCO2sd;
}
transformed data {
}
parameters {
    real<lower=0> eCO2_hat_k1;
    real<lower=0> eCO2_hat_k2;
    real<lower=0> eCO2_hat_alpha21;
    real<lower=0> eCO2_hat_alpha12;
    real<lower=0, upper=1> eCO2_hat_gamma;
    real<lower=0> sigma;
    vector<lower=0>[num_elements(ts)] eCO2;
}
transformed parameters {
    vector[(1 + 1)] eCO2_hat_carbon0 = append_row(rep_vector((eCO2_hat_gamma * totalC_t0), 1), ((1 - eCO2_hat_gamma) * totalC_t0));
    array[ts_n] vector[(1 + 1)] eCO2_hat_trajectory = ode_rk45(
        two_pool_feedback,
        eCO2_hat_carbon0,
        t0,
        to_array_1d(ts),
        eCO2_hat_k1,
        eCO2_hat_k2,
        eCO2_hat_alpha21,
        eCO2_hat_alpha12
    );
    vector[ts_n] eCO2_hat = ((totalC_t0 - to_vector(eCO2_hat_trajectory[:, 1])) - to_vector(eCO2_hat_trajectory[:, 2]));
}
model {
    eCO2_hat_k1 ~ std_normal();
    eCO2_hat_k2 ~ std_normal();
    eCO2_hat_alpha21 ~ std_normal();
    eCO2_hat_alpha12 ~ std_normal();
    eCO2_hat_gamma ~ beta(10, 1);
    sigma ~ cauchy(0, 1);
    eCO2 ~ normal(eCO2_hat, sigma);
    eCO2mean ~ normal(eCO2, eCO2sd);
}
generated quantities {
    vector[eCO2mean_n] eCO2mean_likelihood = normal_lpdfs(eCO2mean, eCO2, eCO2sd);
    vector[eCO2mean_n] eCO2mean_gen = normal_vector_rng(eCO2mean_n, eCO2, eCO2sd);
}

These pages reproduce the model definitions rather than the original case study's fitting and plotting workflow. The source case study remains the place to compare posterior predictive trajectories and diagnose how the two error models allocate uncertainty.