Skip to content

Interventions API

The src.interventions module models super-greenhouse gas (GHG) injection and its radiative effects on the Martian climate. Radiative forcing efficiencies are calibrated to Mars conditions (Marinova et al. 2005) — they differ from Earth IPCC values because Mars lacks water-vapour overlap bands and has a much thinner CO₂ column.

Radiative Forcing Model

The net greenhouse amplification from injected GHGs:

\[ \Delta F_\text{total} = \sum_i \eta_i \cdot C_i \]

where \(\eta_i\) is the Mars-specific radiative forcing efficiency (W m⁻² ppb⁻¹) and \(C_i\) is the concentration in ppb.

The updated greenhouse factor fed into the temperature ODE:

\[ \gamma_\text{new} = \gamma_\text{base} \cdot \left(1 + \frac{\Delta F_\text{total}}{F_\text{ref}}\right) \]

Supported Compounds

Compound Lifetime (yr) Notes
CF4 >50,000 Extremely long-lived
SF6 3,200 High GWP
C2F6 10,000 —
NF3 500 —
C3F8, CHF3, CH2F2, CH3F varies HFC/PFC family
CH4 12 Short-lived, synergistic
N2O 114 —
C3H8, CH3Cl varies —

See list_compounds() for current values loaded at runtime.

GHG Compounds Registry

src.interventions.compounds

GHG compound registry for Mars terraforming interventions.

Radiative forcing efficiencies are from Marinova et al. (2005), "Radiative-convective model of warming Mars with artificial greenhouse gases", Journal of Geophysical Research, Table 2. Mars-specific values differ from Earth (IPCC AR6) because Mars lacks water-vapour overlap bands and has a thinner CO₂ column.

Atmospheric lifetimes are from IPCC AR5/AR6 and Ravishankara et al. (1993).

CompoundProperties dataclass

Physical and radiative properties of a super-greenhouse gas.

Attributes

molecular_weight_g_mol : float Molecular weight in g/mol. atmospheric_lifetime_yr : float First-order e-folding decay lifetime (years). On Mars, photolysis and chemistry are slower than on Earth, so lifetimes are longer. rf_efficiency_W_m2_ppb : float Radiative forcing efficiency in W m⁻² per ppb (mixing ratio). Marinova 2005 Mars-specific values. gwp100 : float Global Warming Potential relative to CO₂ over 100 years (Earth-based reference, provided for familiarity). description : str Human-readable name.

Source code in package/src/interventions/compounds.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
@dataclass(frozen=True)
class CompoundProperties:
    """Physical and radiative properties of a super-greenhouse gas.

    Attributes
    ----------
    molecular_weight_g_mol : float
        Molecular weight in g/mol.
    atmospheric_lifetime_yr : float
        First-order e-folding decay lifetime (years).  On Mars, photolysis
        and chemistry are slower than on Earth, so lifetimes are longer.
    rf_efficiency_W_m2_ppb : float
        Radiative forcing efficiency in W m⁻² per ppb (mixing ratio).
        Marinova 2005 Mars-specific values.
    gwp100 : float
        Global Warming Potential relative to CO₂ over 100 years (Earth-based
        reference, provided for familiarity).
    description : str
        Human-readable name.
    """

    molecular_weight_g_mol:  float
    atmospheric_lifetime_yr: float
    rf_efficiency_W_m2_ppb:  float
    gwp100:                  float
    description:             str

get_compound(name)

Return properties for a compound by name.

Raises

KeyError If the compound is not in the registry. Available names are listed by :func:list_compounds.

Source code in package/src/interventions/compounds.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def get_compound(name: str) -> CompoundProperties:
    """Return properties for a compound by name.

    Raises
    ------
    KeyError
        If the compound is not in the registry.  Available names are
        listed by :func:`list_compounds`.
    """
    if name not in COMPOUNDS:
        available = ", ".join(sorted(COMPOUNDS))
        raise KeyError(
            f"Unknown compound '{name}'. Available: {available}"
        )
    return COMPOUNDS[name]

list_compounds()

Return sorted list of registered compound names.

Source code in package/src/interventions/compounds.py 
133
134
135
def list_compounds() -> list[str]:
    """Return sorted list of registered compound names."""
    return sorted(COMPOUNDS)

Radiative Forcing

src.interventions.forcing

Radiative forcing from super-greenhouse gases — Mars.

Converts atmospheric GHG masses to concentration (ppb), computes total radiative forcing ΔF (W/m²), and updates the planet's greenhouse factor.

All tensor operations stay on the planet's device. Python-float constants are broadcast safely by PyTorch without allocating CPU intermediates.

References

Marinova et al. (2005), "Radiative-convective model of warming Mars with artificial greenhouse gases", J. Geophys. Res., 110, E03002.

compute_concentration_ppb(ghg_state, total_atm_mass_kg)

Convert atmospheric masses to molar mixing ratios (ppb).

Uses the molar-fraction definition:

ppb_i = (M_i / MW_i) / (M_atm / MW_atm) × 1e9

where M_i is the atmospheric mass of compound i (kg), MW_i is its molecular weight (g/mol), M_atm is the total atmospheric mass (kg), and MW_atm is the mean molecular weight of the Mars atmosphere (g/mol).

Parameters

ghg_state : GHGState Current atmospheric GHG mass state. total_atm_mass_kg : torch.Tensor Total Mars atmospheric mass in kg (from mars.atmosphere.atmospheric_mass).

Returns

dict[str, torch.Tensor] Per-compound concentration in ppb, on the same device as ghg_state.

Source code in package/src/interventions/forcing.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def compute_concentration_ppb(
    ghg_state: GHGState,
    total_atm_mass_kg: torch.Tensor,
) -> dict[str, torch.Tensor]:
    """Convert atmospheric masses to molar mixing ratios (ppb).

    Uses the molar-fraction definition:

        ppb_i = (M_i / MW_i) / (M_atm / MW_atm) × 1e9

    where M_i is the atmospheric mass of compound i (kg), MW_i is its
    molecular weight (g/mol), M_atm is the total atmospheric mass (kg),
    and MW_atm is the mean molecular weight of the Mars atmosphere (g/mol).

    Parameters
    ----------
    ghg_state : GHGState
        Current atmospheric GHG mass state.
    total_atm_mass_kg : torch.Tensor
        Total Mars atmospheric mass in kg (from mars.atmosphere.atmospheric_mass).

    Returns
    -------
    dict[str, torch.Tensor]
        Per-compound concentration in ppb, on the same device as ghg_state.
    """
    ppb: dict[str, torch.Tensor] = {}
    for name in ghg_state.compounds:
        mass_kg  = ghg_state.get_mass_kg(name)                    # kg, on device
        mw_ratio = get_compound(name).molecular_weight_g_mol / _MARS_ATM_MW_G_MOL
        ppb[name] = (mass_kg / total_atm_mass_kg) / mw_ratio * 1e9
    return ppb

delta_F_total(ghg_state, total_atm_mass_kg)

Compute total radiative forcing from all tracked GHGs (W/m²).

ΔF = Σ_i η_i [W m⁻² ppb⁻¹] × C_i [ppb]

This is the linear (optically thin) approximation appropriate for trace gases at concentrations < ~1000 ppb, which covers realistic injection rates over 50–100 year horizons.

Parameters

ghg_state : GHGState Current atmospheric GHG mass state. total_atm_mass_kg : torch.Tensor Total Mars atmospheric mass in kg.

Returns

torch.Tensor Scalar total forcing (W/m²) on the GHGState's device.

Source code in package/src/interventions/forcing.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def delta_F_total(
    ghg_state: GHGState,
    total_atm_mass_kg: torch.Tensor,
) -> torch.Tensor:
    """Compute total radiative forcing from all tracked GHGs (W/m²).

    ΔF = Σ_i  η_i [W m⁻² ppb⁻¹]  ×  C_i [ppb]

    This is the linear (optically thin) approximation appropriate for trace
    gases at concentrations < ~1000 ppb, which covers realistic injection
    rates over 50–100 year horizons.

    Parameters
    ----------
    ghg_state : GHGState
        Current atmospheric GHG mass state.
    total_atm_mass_kg : torch.Tensor
        Total Mars atmospheric mass in kg.

    Returns
    -------
    torch.Tensor
        Scalar total forcing (W/m²) on the GHGState's device.
    """
    device = ghg_state.device
    dF     = torch.zeros((), dtype=TF_DTYPE, device=device)
    ppb    = compute_concentration_ppb(ghg_state, total_atm_mass_kg)

    for name, conc in ppb.items():
        eta  = get_compound(name).rf_efficiency_W_m2_ppb   # Python float
        dF   = dF + eta * conc

    return dF

delta_F_from_composition(composition, total_pressure)

Compute total GHG radiative forcing from atmosphere composition (W/m²).

Uses partial pressures (Pa) to derive mole-fraction concentrations (ppb), then applies compound-specific RF efficiencies:

ppb_i  =  P_i / P_total × 10⁹
ΔF     =  Σ_i  η_i × ppb_i

Only species present in the COMPOUNDS registry contribute — background gases (CO₂, N₂, Ar) are silently skipped.

Parameters

composition : dict[str, torch.Tensor] Species → partial pressure (Pa), as stored in mars.atmosphere.composition. total_pressure : torch.Tensor Total atmospheric surface pressure (Pa).

Source code in package/src/interventions/forcing.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def delta_F_from_composition(
    composition: dict,
    total_pressure: torch.Tensor,
) -> torch.Tensor:
    """Compute total GHG radiative forcing from atmosphere composition (W/m²).

    Uses partial pressures (Pa) to derive mole-fraction concentrations (ppb),
    then applies compound-specific RF efficiencies:

        ppb_i  =  P_i / P_total × 10⁹
        ΔF     =  Σ_i  η_i × ppb_i

    Only species present in the COMPOUNDS registry contribute — background
    gases (CO₂, N₂, Ar) are silently skipped.

    Parameters
    ----------
    composition : dict[str, torch.Tensor]
        Species → partial pressure (Pa), as stored in ``mars.atmosphere.composition``.
    total_pressure : torch.Tensor
        Total atmospheric surface pressure (Pa).
    """
    from src.interventions.compounds import COMPOUNDS, get_compound
    device = total_pressure.device
    dF = torch.zeros((), dtype=TF_DTYPE, device=device)
    for name, P_i in composition.items():
        if name in COMPOUNDS:
            ppb = P_i.to(device) / total_pressure * 1e9
            eta = get_compound(name).rf_efficiency_W_m2_ppb
            dF = dF + eta * ppb
    return dF

update_greenhouse_factor(mars, delta_F, baseline_ghf=None, baseline_olr=None)

Apply cumulative GHG radiative forcing via the energy-balance GHF formula.

Derivation

At the original (CO₂-only) equilibrium::

F_in  =  OLR_baseline  =  ε σ (T_eq / GHF_base)⁴          (1)

Injected GHGs trap ΔF extra W/m². At the new radiative equilibrium the atmosphere must emit more to compensate — T rises until::

F_in + ΔF  =  ε σ (T_eq_new / GHF_base)⁴                  (2)

Combining (1) and (2) gives the new effective greenhouse factor::

GHF_new  =  GHF_base × (1 + ΔF / F_in_base)^(1/4)          (3)

where F_in_base = ε σ (T₀ / GHF_base)⁴ is the baseline OLR (≈ mean absorbed solar flux for the initial conditions).

This formula: * Is independent of instantaneous surface temperature — no singularity when Mars is in a cold polar winter (T ≈ 149 K). * Grows monotonically with cumulative ΔF — no year-to-year oscillation. * Is stable for any physically reachable ΔF over 50–200 year horizons.

Parameters

mars : Mars The Mars instance whose greenhouse_factor will be updated. delta_F : torch.Tensor Cumulative total forcing from all GHGs (W/m²), on mars._device. baseline_ghf : torch.Tensor, optional Initial CO₂-only GHF. Callers must pass the value saved at t=0 so that ΔF is applied relative to the fixed baseline (not compounded). Defaults to the current mars.thermal.greenhouse_factor (first call only). baseline_olr : torch.Tensor, optional Pre-computed F_in_base = ε σ (T₀/GHF_base)⁴ in W/m². Callers should cache this at initialisation to avoid recomputing it at every annual update. Computed on-the-fly if not provided.

Source code in package/src/interventions/forcing.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
def update_greenhouse_factor(
    mars,
    delta_F: torch.Tensor,
    baseline_ghf: torch.Tensor | None = None,
    baseline_olr: torch.Tensor | None = None,
) -> None:
    """Apply cumulative GHG radiative forcing via the energy-balance GHF formula.

    Derivation
    ----------
    At the original (CO₂-only) equilibrium::

        F_in  =  OLR_baseline  =  ε σ (T_eq / GHF_base)⁴          (1)

    Injected GHGs trap ΔF extra W/m².  At the **new** radiative equilibrium
    the atmosphere must emit more to compensate — T rises until::

        F_in + ΔF  =  ε σ (T_eq_new / GHF_base)⁴                  (2)

    Combining (1) and (2) gives the new effective greenhouse factor::

        GHF_new  =  GHF_base × (1 + ΔF / F_in_base)^(1/4)          (3)

    where F_in_base = ε σ (T₀ / GHF_base)⁴ is the baseline OLR (≈ mean
    absorbed solar flux for the initial conditions).

    This formula:
    * Is **independent of instantaneous surface temperature** — no singularity
      when Mars is in a cold polar winter (T ≈ 149 K).
    * Grows **monotonically** with cumulative ΔF — no year-to-year oscillation.
    * Is **stable** for any physically reachable ΔF over 50–200 year horizons.

    Parameters
    ----------
    mars : Mars
        The Mars instance whose greenhouse_factor will be updated.
    delta_F : torch.Tensor
        **Cumulative** total forcing from all GHGs (W/m²), on mars._device.
    baseline_ghf : torch.Tensor, optional
        Initial CO₂-only GHF.  Callers must pass the value saved at t=0 so
        that ΔF is applied relative to the fixed baseline (not compounded).
        Defaults to the current mars.thermal.greenhouse_factor (first call only).
    baseline_olr : torch.Tensor, optional
        Pre-computed F_in_base = ε σ (T₀/GHF_base)⁴ in W/m².
        Callers should cache this at initialisation to avoid recomputing
        it at every annual update.  Computed on-the-fly if not provided.
    """
    if float(delta_F.item()) <= 0.0:
        return   # no forcing to apply

    d        = mars._device
    sb       = STEFAN_BOLTZMANN.to(d)

    GHF_base = (baseline_ghf if baseline_ghf is not None
                else mars.thermal.greenhouse_factor)

    if baseline_olr is not None:
        F_in_base = baseline_olr
    else:
        T0        = mars.thermal.surface_temperature
        F_in_base = _MARS_EMISSIVITY * sb * (T0 / GHF_base) ** 4.0

    # Equation (3): monotonic, singularity-free
    GHF_new = GHF_base * (1.0 + delta_F / F_in_base) ** 0.25

    mars.thermal.greenhouse_factor = GHF_new.clamp(min=1.0)

Intervention Controller

src.interventions.controller

InterventionController — annual GHG injection scheduler for Mars.

The controller drives terraforming experiments by injecting super-greenhouse gases on a yearly schedule. All atmospheric state lives in mars.atmosphere.composition — injected GHGs are just additional entries in the same dict alongside CO₂, N₂, Ar. There is no separate GHGState object.

Architecture

mars.atmosphere.composition = {
    "CO2": 580 Pa,   # pre-existing
    "N2":   15 Pa,   # pre-existing
    "CF4":   X Pa,   # added by inject()
    "SF6":   Y Pa,   # added by inject()
    ...
}

mars.inject(schedule) converts kg → Pa and adds directly to composition. mars.decay_ghg(dt_years) decays all COMPOUNDS-registered species in place. mars.thermal.greenhouse_factor and mars.delta_F are always current.

The controller adds only one piece of state not derivable from physics: cumulative_injected_kg — how much total mass has ever been injected per species. This is a reporting quantity, not physics state.

InterventionSnapshot dataclass

Climate state at the end of one Mars year of intervention.

All tensor fields live on the planet's device.

Source code in package/src/interventions/controller.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
@dataclass
class InterventionSnapshot:
    """Climate state at the end of one Mars year of intervention.

    All tensor fields live on the planet's device.
    """

    year:                     int                         # 1-based year counter
    time_s:                   torch.Tensor                # elapsed seconds

    # Climate state (mirrors TimeController Snapshot for compatibility)
    surface_temperature:      torch.Tensor                # K  (annual mean)
    temp_min:                 torch.Tensor                # K  (annual minimum)
    temp_max:                 torch.Tensor                # K  (annual maximum)
    surface_pressure:         torch.Tensor                # Pa (annual mean)
    ice_mass:                 torch.Tensor                # kg (end-of-year)
    solar_flux:               torch.Tensor                # W m⁻² (annual mean)
    orbital_angle:            torch.Tensor                # rad (end-of-year)
    greenhouse_factor:        torch.Tensor                # GHF at end of year

    # Intervention-specific
    delta_F:                  torch.Tensor                # total GHG forcing (W/m²)
    ghg_partial_pressure_Pa:  dict[str, torch.Tensor]     # Pa per GHG compound
    cumulative_injected_kg:   dict[str, torch.Tensor]     # total injected to date

    @property
    def time(self) -> torch.Tensor:
        """Alias for time_s — keeps compatibility with Snapshot consumers."""
        return self.time_s

time property

Alias for time_s — keeps compatibility with Snapshot consumers.

InterventionController

Annual GHG injection scheduler + climate simulation driver.

Injected gases become part of mars.atmosphere.composition — the same dict that holds CO₂, N₂, and Ar. No separate atmospheric bookkeeper exists.

Parameters

mars : Mars Configured Mars instance. injection_schedule_kg_yr : dict[str, float] {compound_name: kg_per_Mars_year} for each injected species. All compound names must be present in the COMPOUNDS registry. dt : float Sub-annual integration timestep in seconds. Default 3600 s (1 hour). accuracy : Accuracy Integration accuracy passed to the underlying TimeController. compile : bool If True, wrap physics methods with torch.compile.

Examples

from src.celestials import Mars from src.interventions import InterventionController mars = Mars() ic = InterventionController(mars, {"CF4": 1e9, "SF6": 5e8}, dt=21600) history = ic.run(n_years=50)

mars.atmosphere.composition now contains CF4 and SF6 partial pressures

print(mars.atmosphere.composition.keys()) print(f"ΔF = {mars.delta_F.item():.3f} W/m²")

Source code in package/src/interventions/controller.py 
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
class InterventionController:
    """Annual GHG injection scheduler + climate simulation driver.

    Injected gases become part of ``mars.atmosphere.composition`` — the same
    dict that holds CO₂, N₂, and Ar.  No separate atmospheric bookkeeper exists.

    Parameters
    ----------
    mars : Mars
        Configured Mars instance.
    injection_schedule_kg_yr : dict[str, float]
        ``{compound_name: kg_per_Mars_year}`` for each injected species.
        All compound names must be present in the COMPOUNDS registry.
    dt : float
        Sub-annual integration timestep in seconds.  Default 3600 s (1 hour).
    accuracy : Accuracy
        Integration accuracy passed to the underlying TimeController.
    compile : bool
        If True, wrap physics methods with torch.compile.

    Examples
    --------
    >>> from src.celestials import Mars
    >>> from src.interventions import InterventionController
    >>> mars = Mars()
    >>> ic = InterventionController(mars, {"CF4": 1e9, "SF6": 5e8}, dt=21600)
    >>> history = ic.run(n_years=50)
    >>> # mars.atmosphere.composition now contains CF4 and SF6 partial pressures
    >>> print(mars.atmosphere.composition.keys())
    >>> print(f"ΔF = {mars.delta_F.item():.3f} W/m²")
    """

    def __init__(
        self,
        mars,
        injection_schedule_kg_yr: dict[str, float],
        dt: float = 3600.0,
        accuracy: Accuracy = Accuracy.FAST,
        compile: bool = False,
    ) -> None:
        for name in injection_schedule_kg_yr:
            get_compound(name)

        self._mars     = mars
        self._device   = mars._device
        self._schedule = {k: float(v) for k, v in injection_schedule_kg_yr.items()}
        self._tc       = TimeController(mars, dt=dt, accuracy=accuracy, compile=compile)

        # Reporting-only: cumulative kg injected per compound (not physics state)
        self._cumulative_injected_kg: dict[str, torch.Tensor] = {
            k: torch.zeros((), dtype=TF_DTYPE, device=self._device)
            for k in self._schedule
        }

        # Prime baselines on Mars so delta_F is well-defined from run() start
        if self._schedule:
            mars._init_ghg()

        self._elapsed     = torch.zeros((), dtype=TF_DTYPE, device=self._device)
        self._year_s      = float(_MARS_YEAR_T.item())
        self.all_hourly: list = []   # full hourly trace accumulated across all years

    # ------------------------------------------------------------------
    # Main loop
    # ------------------------------------------------------------------

    def run(
        self,
        n_years: int,
        callback: Optional[Callable[[InterventionSnapshot], None]] = None,
    ) -> list[InterventionSnapshot]:
        """Run the intervention simulation for *n_years* Mars years.

        Per year:
            1. mars.inject(schedule)        → partial pressures updated in composition
            2. TimeController.run(1 year)   → ODE integrates with updated GHF
            3. mars.decay_ghg(1.0)          → COMPOUNDS species decay in composition
            4. InterventionSnapshot         → read from mars state

        Returns
        -------
        list[InterventionSnapshot]
            One snapshot per year (length = n_years).
        """
        history: list[InterventionSnapshot] = []

        for year in range(1, n_years + 1):
            # Step 1 — inject; GHF synced automatically inside mars.inject()
            if self._schedule:
                self._mars.inject(self._schedule)
                for name, kg in self._schedule.items():
                    self._cumulative_injected_kg[name] = (
                        self._cumulative_injected_kg[name] + kg
                    )

            # Step 2 — simulate one Mars year
            year_history = self._tc.run(duration=self._year_s)
            self.all_hourly.extend(year_history)
            self._elapsed = self._elapsed + self._year_s

            # Step 3 — decay; GHF resynced automatically inside mars.decay_ghg()
            self._mars.decay_ghg(dt_years=1.0)

            # Step 4 — snapshot from Mars (single source of truth)
            snap = self._make_snapshot(year, year_history)
            history.append(snap)

            if callback is not None:
                callback(snap)

        return history

    # ------------------------------------------------------------------
    # Private helpers
    # ------------------------------------------------------------------

    def _make_snapshot(self, year: int, year_history: list) -> InterventionSnapshot:
        def _mean(attr: str) -> torch.Tensor:
            return torch.stack([getattr(s, attr) for s in year_history]).mean()

        temps = torch.stack([s.surface_temperature for s in year_history])

        ghg_pp = {
            name: P.clone()
            for name, P in self._mars.atmosphere.composition.items()
            if name in COMPOUNDS
        }
        return InterventionSnapshot(
            year                    = year,
            time_s                  = self._elapsed.clone(),
            surface_temperature     = temps.mean(),
            temp_min                = temps.min(),
            temp_max                = temps.max(),
            surface_pressure        = _mean("surface_pressure"),
            ice_mass                = year_history[-1].ice_mass.clone(),
            solar_flux              = _mean("solar_flux"),
            orbital_angle           = year_history[-1].orbital_angle.clone(),
            greenhouse_factor       = self._mars.thermal.greenhouse_factor.clone(),
            delta_F                 = self._mars.delta_F.clone(),
            ghg_partial_pressure_Pa = ghg_pp,
            cumulative_injected_kg  = {
                k: v.clone() for k, v in self._cumulative_injected_kg.items()
            },
        )

run(n_years, callback=None)

Run the intervention simulation for n_years Mars years.

Per year
  1. mars.inject(schedule) → partial pressures updated in composition
  2. TimeController.run(1 year) → ODE integrates with updated GHF
  3. mars.decay_ghg(1.0) → COMPOUNDS species decay in composition
  4. InterventionSnapshot → read from mars state
Returns

list[InterventionSnapshot] One snapshot per year (length = n_years).

Source code in package/src/interventions/controller.py 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def run(
    self,
    n_years: int,
    callback: Optional[Callable[[InterventionSnapshot], None]] = None,
) -> list[InterventionSnapshot]:
    """Run the intervention simulation for *n_years* Mars years.

    Per year:
        1. mars.inject(schedule)        → partial pressures updated in composition
        2. TimeController.run(1 year)   → ODE integrates with updated GHF
        3. mars.decay_ghg(1.0)          → COMPOUNDS species decay in composition
        4. InterventionSnapshot         → read from mars state

    Returns
    -------
    list[InterventionSnapshot]
        One snapshot per year (length = n_years).
    """
    history: list[InterventionSnapshot] = []

    for year in range(1, n_years + 1):
        # Step 1 — inject; GHF synced automatically inside mars.inject()
        if self._schedule:
            self._mars.inject(self._schedule)
            for name, kg in self._schedule.items():
                self._cumulative_injected_kg[name] = (
                    self._cumulative_injected_kg[name] + kg
                )

        # Step 2 — simulate one Mars year
        year_history = self._tc.run(duration=self._year_s)
        self.all_hourly.extend(year_history)
        self._elapsed = self._elapsed + self._year_s

        # Step 3 — decay; GHF resynced automatically inside mars.decay_ghg()
        self._mars.decay_ghg(dt_years=1.0)

        # Step 4 — snapshot from Mars (single source of truth)
        snap = self._make_snapshot(year, year_history)
        history.append(snap)

        if callback is not None:
            callback(snap)

    return history