rscm_components/components/
carbon_cycle.rs

1//! Carbon cycle component
2//!
3//! A simple one-box carbon cycle model that tracks atmospheric CO2 concentrations
4//! and land uptake based on emissions and temperature.
5
6use crate::constants::GTC_PER_PPM;
7use ode_solvers::Vector3;
8use rscm_core::component::{
9    Component, GridType, InputState, OutputState, RequirementDefinition, RequirementType,
10    ScalarWindow,
11};
12use rscm_core::errors::RSCMResult;
13use rscm_core::ivp::{get_last_step, IVPBuilder, IVP};
14use rscm_core::state::StateValue;
15use rscm_core::timeseries::{FloatValue, Time};
16use rscm_core::ComponentIO;
17use serde::{Deserialize, Serialize};
18use std::sync::Arc;
19
20type ModelState = Vector3<FloatValue>;
21
22/// Parameters for the one-box carbon cycle component
23#[derive(Debug, Clone, Serialize, Deserialize)]
24pub struct CarbonCycleParameters {
25    /// Timescale of the box's response
26    /// unit: yr
27    pub tau: FloatValue,
28    /// Pre-industrial atmospheric CO2 concentration
29    /// unit: ppm
30    pub conc_pi: FloatValue,
31    /// Sensitivity of lifetime to changes in global-mean temperature
32    /// unit: 1 / K
33    pub alpha_temperature: FloatValue,
34}
35
36/// Solver options for the ODE integration
37#[derive(Debug, Clone, Serialize, Deserialize)]
38pub struct SolverOptions {
39    pub step_size: FloatValue,
40}
41
42/// One-box carbon cycle component
43///
44/// This component models the carbon cycle using a simple one-box model where:
45/// - CO2 emissions increase atmospheric concentrations
46/// - Land uptake removes CO2 at a rate that depends on the concentration anomaly
47/// - The uptake rate is temperature-dependent
48///
49/// The governing equations are:
50/// $$ \frac{dC}{dt} = E - \frac{C - C_0}{\tau \exp(\alpha_T \cdot T)} $$
51///
52/// Where:
53/// - $C$ is atmospheric CO2 concentration (ppm)
54/// - $E$ is emissions (GtC/yr converted to ppm/yr)
55/// - $C_0$ is pre-industrial concentration (ppm)
56/// - $\tau$ is the baseline lifetime (yr)
57/// - $\alpha_T$ is the temperature sensitivity (1/K)
58/// - $T$ is the surface temperature anomaly (K)
59#[derive(Debug, Clone, Serialize, Deserialize, ComponentIO)]
60#[component(tags = ["carbon-cycle", "simple", "stable"], category = "Carbon Cycle")]
61#[inputs(
62    emissions { name = "Emissions|CO2|Anthropogenic", unit = "GtC / yr" },
63    temperature { name = "Surface Temperature", unit = "K" },
64)]
65#[states(
66    concentration { name = "Atmospheric Concentration|CO2", unit = "ppm" },
67    cumulative_emissions { name = "Cumulative Emissions|CO2", unit = "Gt C" },
68    cumulative_uptake { name = "Cumulative Land Uptake", unit = "Gt C" },
69)]
70pub struct CarbonCycle {
71    parameters: CarbonCycleParameters,
72    solver_options: SolverOptions,
73}
74
75impl CarbonCycle {
76    /// Create a new carbon cycle component from parameters
77    pub fn from_parameters(parameters: CarbonCycleParameters) -> Self {
78        Self {
79            parameters,
80            solver_options: SolverOptions { step_size: 0.1 },
81        }
82    }
83
84    /// Set custom solver options
85    pub fn with_solver_options(self, solver_options: SolverOptions) -> Self {
86        Self {
87            parameters: self.parameters,
88            solver_options,
89        }
90    }
91}
92
93#[typetag::serde]
94impl Component for CarbonCycle {
95    fn definitions(&self) -> Vec<RequirementDefinition> {
96        Self::generated_definitions()
97    }
98
99    fn solve(
100        &self,
101        t_current: Time,
102        t_next: Time,
103        input_state: &InputState,
104    ) -> RSCMResult<OutputState> {
105        let inputs = CarbonCycleInputs::from_input_state(input_state);
106
107        let y0 = ModelState::new(
108            inputs.concentration.at_start(),
109            inputs.cumulative_uptake.at_start(),
110            inputs.cumulative_emissions.at_start(),
111        );
112
113        let solver = IVPBuilder::new(Arc::new(self.to_owned()), input_state, y0);
114
115        let mut solver = solver.to_rk4(t_current, t_next, self.solver_options.step_size);
116        solver.integrate().expect("Failed solving");
117
118        let results = get_last_step(solver.results(), t_next);
119
120        let outputs = CarbonCycleOutputs {
121            concentration: results[0],
122            cumulative_uptake: results[1],
123            cumulative_emissions: results[2],
124        };
125
126        Ok(outputs.into())
127    }
128}
129
130impl IVP<Time, ModelState> for CarbonCycle {
131    fn calculate_dy_dt(
132        &self,
133        _t: Time,
134        input_state: &InputState,
135        y: &Vector3<FloatValue>,
136        dy_dt: &mut Vector3<FloatValue>,
137    ) {
138        let inputs = CarbonCycleInputs::from_input_state(input_state);
139
140        // Inputs come from input_state (exogenous data at start of timestep)
141        let emissions = inputs.emissions.at_start();
142        let temperature = inputs.temperature.at_start();
143
144        // State variables come from the ODE state vector y
145        let conc = y[0];
146
147        // dC / dt = E - (C - C_0) / (tau * exp(alpha_temperature * temperature))
148        let lifetime =
149            self.parameters.tau * (self.parameters.alpha_temperature * temperature).exp();
150        let uptake = (conc - self.parameters.conc_pi) / lifetime; // ppm / yr
151
152        dy_dt[0] = emissions / GTC_PER_PPM - uptake; // ppm / yr
153        dy_dt[1] = uptake * GTC_PER_PPM; // GtC / yr
154        dy_dt[2] = emissions // GtC / yr
155    }
156}
rscm_components/components/carbon_cycle.rs

rscm_components/components/
carbon_cycle.rs
