This protocol is a simplified two-box model of the global carbon cycle. It simulates how volcanic degassing of CO₂ is partitioned between the atmosphere and the ocean when no long-term sink (weathering) is present.
The key idea is:
- Without weathering, CO₂ accumulates in the system indefinitely.
- The atmosphere may temporarily hold more CO₂ than its 'equilibrium share' if exchange with the ocean is slow.
- This shows why silicate weathering, the long-term sink, is essential for stabilizing climate.
We represent two reservoirs:
1. Atmosphere (A, in ppm) – the concentration of CO₂ in the air.
2. Ocean (O, in GtCO₂) – the ocean inventory of dissolved CO₂.
Parameters:
- f_atm_eq: the equilibrium fraction of the total carbon that resides in the atmosphere (e.g. 0.05 = 5%).
- tau_ex (τ): the characteristic exchange time (yr) between atmosphere and ocean.
- D: volcanic degassing input (MtCO₂/yr), constant until a cutoff time.
- A*: the target atmospheric CO₂ (ppm) that we want to investigate.
- t_off: the cutoff time when degassing is switched off (optional).
Input adds CO₂ to the system:
D_ppm = D / 7800
Exchange moves CO₂ between A and O toward equilibrium:
F_ex = (A - f * C_tot) / tau
Without weathering:
C_tot = A + O/7.8 (in ppm-equivalents)
1. Choose parameters:
- A small or large f.
- A fast (tau = 50 yr) or slow (tau = 5000 yr) exchange.
- A volcanic input (e.g. 300 Mt/yr).
- A target CO₂ level (e.g. 500 ppm or 1000 ppm).
2. Run the model and read:
- t_hit: time when the atmosphere first reaches the target.
- t_off: time when degassing is switched off.
- A_eq_final: the final atmospheric concentration after relaxation.
- t_to_95eq, t_to_99eq: time after cutoff needed to reach 95% and 99% of equilibrium.
3. Interpret the graphs:
- Before cutoff: atmosphere rises as volcanic input dominates.
- After cutoff: atmosphere relaxes toward equilibrium with the ocean.
Example 1: Fast exchange, small atmosphere
- f = 0.02, tau = 50 yr, D = 300 Mt/yr, A* = 500 ppm.
- Atmosphere reaches 500 ppm in about 650 kyr; equilibrium after cutoff is 500 ppm.
- Lesson: with fast exchange, the atmosphere never overshoots much.
Example 2: Slow exchange, small atmosphere
- f = 0.02, tau = 5000 yr, D = 300 Mt/yr, A* = 500 ppm.
- Atmosphere reaches 500 ppm in ~120 kyr, much faster than equilibrium theory predicts.
- Equilibrium after cutoff is higher than 500 ppm unless cutoff is tuned.
- Lesson: with slow exchange, the atmosphere gets 'bloated' with CO₂.
1. Compare two scenarios: tau = 50 yr and tau = 5000 yr, same f and D.
- Question: Which system reaches 500 ppm faster? Why?
2. Equilibrium vs. overshoot: f=0.05, D=300 Mt/yr, cutoff at 100 kyr.
- Question: What is the final equilibrium A_eq? Is it the same as at cutoff?
3. Importance of weathering: Run with input left on forever.
- Question: Why does CO₂ keep rising without bound? What process is missing in this model?
This model shows that without weathering, volcanic CO₂ would accumulate endlessly. The atmosphere can fill up quickly—sometimes much faster than expected—if exchange with the ocean is slow. Only silicate weathering provides the long-term negative feedback that stabilizes climate on million-year timescales.
The Walker thermostat describes the long-term climate regulation mechanism driven by the feedback between volcanic CO₂ outgassing and silicate weathering (Walker, Hays & Kasting, 1981). Over millions of years, this negative feedback stabilizes Earth’s surface temperature. This exercise uses the interactive Python/Colab model of the Walker thermostat. The goal is to understand how volcanic outgassing, CO2 concentration, and surface temperature interact in long-term climate regulation.
How does volcanic CO₂ outgassing affect Earth’s long-term climate?
By changing the volcanic input flux (F), how much must atmospheric CO₂ (P/P₀) shift to restore balance with silicate weathering?
How strong is the climate stabilization feedback from weathering?
By adjusting the CO₂ sensitivity parameter (β) and the temperature sensitivity constant (T*), can the system keep Earth’s temperature anomaly (ΔT) close to zero?
How does climate sensitivity influence equilibrium?
If climate sensitivity for CO₂ doubling (S₂x) is high vs. low, how does this affect the final equilibrium P/P₀ and ΔT?
What is the equilibrium outcome for different scenarios?
For a given combination of F, β, T*, and S₂x, what values of CO₂ (P/P₀) and temperature anomaly (ΔT) emerge as the stable equilibrium?
Weathering flux:
W / W0 = (P / P0)^beta * exp(ΔT / Tstar)
describes how weathering (W) changes compared to today’s reference level (W₀).
The first part, (P/P₀)^β, means weathering speeds up if there is more CO₂ in the air. The higher β is, the stronger this effect.
The second part, exp(ΔT/T*), means weathering also speeds up if the planet gets warmer. The smaller T* is, the more quickly temperature changes affect weathering.
So in plain words:
👉 Weathering today (W) depends on two things — how much CO₂ is in the atmosphere (P/P₀), and how warm the surface is (ΔT). More CO₂ and higher temperature both make weathering faster.
F = normalized volcanic outgassing (F_out / W0). If F = 1, volcanic CO2 release equals present-day reference. If F > 1, more CO2 is released. If F < 1, less CO2 is released.
beta = exponent that controls how strongly weathering depends on CO2. Typical values are between 0.2 and 0.6. Larger beta means stronger sensitivity. If beta is small (e.g., 0.2), the weathering rate increases only weakly when CO₂ rises; the feedback is weaker. If beta is larger (e.g., 0.6), the weathering rate increases much more strongly with the same rise in CO₂; the feedback is stronger. In practice, larger beta means the system “notices” CO₂ changes more strongly, so chemical weathering reacts faster, drawing down CO₂ more efficiently and stabilizing climate more tightly.
Tstar = temperature scaling constant in Kelvin. It represents how strongly weathering reacts to changes in surface temperature. Typical values are 10–20 K. If T* is small (closer to 10 K), the exponential grows quickly: even a modest warming produces a large increase in weathering rate. This means the climate feedback is strong. If T* is large (closer to 20 K or more), the exponential grows more slowly: the same warming produces only a small increase in weathering. The feedback is weaker. In practice, T* controls the thermal sensitivity of the thermostat: it tells us how many degrees of warming are needed to accelerate weathering significantly. A smaller T* makes the system react more sharply to temperature shifts, while a larger T* makes it more sluggish.
β (beta) and T* work like two knobs on the Earth’s climate thermostat.
β (beta) = CO₂ sensitivity.
It controls how strongly silicate weathering responds to the amount of CO₂ in the air.
If β is large, weathering reacts strongly when CO₂ rises, removing it quickly and stabilizing climate.
If β is small, weathering responds weakly, so CO₂ can build up more easily.
T* = temperature sensitivity.
It sets how strongly weathering responds to surface warming.
A small T* means even a small warming makes weathering speed up a lot (strong feedback).
A large T* means weathering reacts slowly to warming (weak feedback).
Together:
High β + low T* = very stable climate (fast reaction to both CO₂ and warming).
Low β + high T* = sluggish climate (weak reaction, bigger swings in CO₂ and temperature).
S2x = climate sensitivity for doubling CO2, measured in Kelvin. It is the temperature increase expected if CO2 doubles. It tells us how much the global mean surface temperature would rise if the atmospheric CO₂ concentration doubled relative to a reference state. For example: If S₂x = 2 K, then doubling CO₂ causes a 2-degree warming. If S₂x = 4.5 K, the same doubling produces a much stronger 4.5-degree warming. Typical values considered realistic are between 2 and 4.5 K.
Lambda = climate sensitivity factor in K per natural logarithm unit of CO2. Calculated as S2x / ln(2). Lambda (Λ) is the climate sensitivity factor: it tells us how many degrees the Earth’s surface temperature changes for a given proportional change in atmospheric CO₂. It is derived directly from S2x (the warming expected for a doubling of CO₂). In practice, if Λ is larger, the climate reacts more strongly to changes in CO₂, so even small shifts in CO₂ produce noticeable temperature changes. If Λ is smaller, the climate response is weaker.
In pills:
S₂x is the warming knob. It tells us how much the planet heats up if CO₂ doubles.
Λ is general: warming for any proportional CO₂ change.
β and T* are the cooling knobs. They control how strongly chemical weathering reacts (to CO₂ and to temperature) to pull CO₂ back down.
So the balance is:
S₂x pushes the climate warmer,
β and T* push back by removing CO₂ faster, helping to cool and stabilize the system.
P/P0 = relative atmospheric CO2 at equilibrium, compared to present-day reference. If P/P0 = 1, CO2 is the same as today. If P/P0 > 1, more CO2. If P/P0 < 1, less CO2. Within the Walker thermostat framework—realistic values of volcanic outgassing (F), weathering sensitivity, and climate sensitivity—plausible equilibrium P/P₀ typically stays within about 0.7 to 1.5: P/P₀ ≈ 0.8 to 1 corresponds to cooling trends: examples include moderate decreases in CO₂ due to reduced volcanism or stronger weathering feedbacks. P/P₀ ≈ 1 to 1.3–1.5 corresponds to warming trends: such as modest increases in CO₂ from elevated volcanic flux or weaker weathering. These ranges are illustrative, not absolute. In extreme model parameter combinations, P/P₀ might venture outside these bands—but for educational purposes, framing a target window of 0.7 to 1.5 helps students gauge reasonable equilibrium outcomes without needing external data.
ΔT = surface temperature anomaly in Kelvin, relative to present. If ΔT = 0, the global mean surface temperature is the same as today. If ΔT > 0, the climate is warmer than today by that many degrees. For example, ΔT = +2 means 2 K warmer. If ΔT < 0, the climate is cooler than today. For example, ΔT = –3 means 3 K colder. In the Walker thermostat model, ΔT is not an arbitrary number: it is the result of the balance between volcanic outgassing and silicate weathering, mediated by CO₂. Larger CO₂ levels (P/P₀ > 1) push ΔT positive, while lower CO₂ levels (P/P₀ < 1) drive ΔT negative. Typical values in classroom exercises are modest, often in the range –3 K to +3 K, which correspond to realistic climate shifts rather than extreme snowball or runaway conditions.
Exercise 1 — Volcanic Outgassing and the Walker Thermostat
In this exercise we test what happens if volcanic CO₂ release (F) is smaller or larger than today, and how the Walker thermostat stabilizes the climate.
Fix the parameters:
β = 0.3, T* = 13.7 K, S₂x = 3.0 K.
From this we calculate Λ = 4.3 K (climate sensitivity factor).
At equilibrium:
(P/P₀)_eq = F^(1 / (β + Λ/T*))
ΔT_eq = Λ · ln(P/P₀)
Task:
Compute equilibrium CO₂ (P/P₀) and temperature anomaly (ΔT) for:
F = 0.8 (weaker volcanism)
F = 1.0 (present volcanism)
F = 1.2 (stronger volcanism)
Solutions:
F = 0.8 → (P/P₀) ≈ 0.70, ΔT ≈ −1.6 K (cooler climate)
F = 1.0 → (P/P₀) = 1.00, ΔT = 0 K (reference climate)
F = 1.2 → (P/P₀) ≈ 1.34, ΔT ≈ +1.3 K (warmer climate)
What to look for:
When outgassing decreases (F < 1), the system settles to lower CO₂ and cooler climate.
When outgassing increases (F > 1), CO₂ and temperature rise.
In both cases, the changes in ΔT are moderated by the thermostat: the Earth does not swing as wildly as the forcing itself, because weathering reacts and damps the effect.
Exercise 2 — Climate sensitivity
Keep F = 1.2, beta = 0.3, Tstar = 13.7.
Compare results for:
S2x = 2.0
S2x = 3.0
S2x = 4.5
Results:
For S₂x = 2.0 → Lambda = 2.9 K → (P/P0) ≈ 1.43, ΔT ≈ +1.0 K
For S₂x = 3.0 → Lambda = 4.3 K → (P/P0) ≈ 1.34, ΔT ≈ +1.3 K
For S₂x = 4.5 → Lambda = 6.5 K → (P/P0) ≈ 1.27, ΔT ≈ +1.5 K
Explain how stronger climate sensitivity modifies equilibrium ΔT. As climate sensitivity (S₂x) increases, the equilibrium warming ΔT gets larger, even though the CO₂ increase (P/P₀) required is smaller. Stronger sensitivity means the same CO₂ perturbation produces more warming, so the thermostat can balance volcanic input with a lower P/P₀.
😱😱😱😱
Use ChatGPT to generate the notebook. Open ChatGPT and paste the following prompt exactly.
Create a single Jupyter/Colab cell that defines a minimal interactive Walker thermostat model (Walker, Hays & Kasting, 1981) with ipywidgets. Use only numpy and matplotlib (no seaborn). Provide two separate plots: (1) equilibrium ΔT (K) vs F, and (2) equilibrium P/P0 vs F, each with a vertical guide at F=1 and a marker at the current state. Implement equations in ASCII:
W/W0 = (P/P0)^beta * exp(ΔT / Tstar)
ΔT = Lambda * ln(P/P0), Lambda = S2x / ln(2)
At steady state: F = W/W0 ⇒ (P/P0) = F^(1/(beta + Lambda/Tstar)), ΔT = Lambda * ln(P/P0)
Expose interactive sliders for:
F in [0.5, 1.5], step 0.01
beta in [0.2, 0.6], step 0.01
Tstar in [8, 25] K, step 0.1
S2x in [1.5, 6.0] K, step 0.1
F_min, F_max to set the plotted F-range.
Print a brief numeric summary each time (beta, Tstar, S2x, Lambda; equilibrium P/P0 and ΔT; a one-line interpretation based on F≶1). Enforce the physical condition beta + Lambda/Tstar > 0 with a clear error message if violated. Do NOT set custom colors or styles. Keep the cell self-contained.
Step 2 — Copy the code ChatGPT returns and move to Google Colab. Open https://colab.research.google.com, create a new notebook, paste the generated code into the first cell, and run it with Shift+Enter.
Step 3 — If sliders do not appear, enable widgets once in Colab. Run the following in a new cell, then restart the runtime from the “Runtime” menu and re-run the notebook from the top.
%pip install ipywidgets==8
from google.colab import output
output.enable_custom_widget_manager()
This guide explains, step by step, every input and output in the single‑cell Colab model for vacuum impacts and the companion crater‑size estimate. Language is simple on purpose. The model is teaching‑level: no atmosphere, no fragmentation, no explicit layering or topography. Use results as order‑of‑magnitude estimates (often within a factor of ~2 under reasonable conditions).
The projectile is a uniform sphere. The target is treated as an airless surface. Energy is purely kinetic:
E = 1/2 * m * v^2
Crater size follows a compact gravity‑dominated scaling with an empirical prefactor K that we calibrate to Chicxulub.
Each target body provides a surface gravity g (m/s^2) and a representative bulk target density rho_t (g/cm^3). These affect only the crater estimate, not the energy.
Step 1 — Choose the target body.
Pick Earth, Moon, Mars, Mercury, Venus, or Ceres. The notebook loads g and rho_t for that body. These values set the gravity context for the crater estimate.
Step 2 — Set the impactor density rho_i (g/cm^3).
This is the bulk density of the projectile. Typical: ~1–2 (porous/icy), ~2.5–3.5 (stony), ~5–8 (metal‑rich). Higher rho_i means more mass at the same size, so more energy and a larger crater.
Step 3 — Set the impactor diameter D_imp (km).
The model assumes a sphere. It computes R = 0.5 * D_imp, V = (4/3) * pi * R^3, and m = rho_i * V. Mass grows with the cube of size. Doubling D_imp makes mass ~8 times larger.
Step 4 — Set the impact speed v (km/s).
This is the speed at the surface just before impact. Energy depends on v^2. Keep v physically plausible for the chosen body.
Step 5 — Set the impact angle theta (degrees from the horizontal).
90° is vertical. 45° is most probable. The crater law includes (sin(theta))^(1/3), so shallow angles reduce crater size. Very shallow angles (<~15°) can produce elongated craters not captured by this simple law.
Step 6 — Optional: override the target density rho_t (g/cm^3).
If you know the local target is lighter or denser than the default, enter a new value. If you leave it as is, the default for the selected target is used.
Step 7 — Set the scaling constant K (dimensionless).
K normalizes the crater‑diameter law and absorbs both transient‑to‑final conversion and geology effects. We calibrate K to Chicxulub. A practical default is K = 13.0, which reproduces a ~190 km final crater for D_imp ≈ 12 km, v ≈ 20 km/s, theta ≈ 60° on an Earth‑like shallow‑marine carbonate–evaporite target. Angle and speed variations move the fitted K within ~12–15 in typical cases.
Step 8 — Set the plot ranges (v_min, v_max in km/s; D_min, D_max in km).
These ranges only define the horizontal axis of the two plots. They do not change the computed outputs. A marker shows the exact point for your current parameters.
Mass (kg):
m = rho_i * (4/3) * pi * (0.5 * D_imp)^3
Kinetic Energy (J):
E = 1/2 * m * v^2
TNT Equivalents:
E_kt = E / 4.184e12
E_Mt = E / 4.184e15
Hiroshima equivalents (15 kt each): N_Hiroshima = E / (15 * 4.184e12)
Momentum and Specific Energy:
p = m * v
q = E / m
Estimated Final Crater Diameter (km):
D_crater = K * (rho_i / rho_t)^(1/3) * (sin(theta))^(1/3) * g^(-0.22) * v^(0.44) * D_imp^(0.78)
Energy vs. speed: the curve shows the v^2 growth of E (reported in Mt TNT). The marker is your chosen speed. Density and size are held at your values.
Crater size vs. impactor size: the curve follows approximately D_imp^0.78. The marker is your chosen D_imp. Speed, angle, densities, gravity and K are held at your values.
Keep units straight. v is in km/s and D_imp is in km in the crater law. Internally the code converts units as needed.
Angle is from the horizontal. If you think in “impact angle from vertical”, replace theta by (90° − your angle).
Do not over‑interpret a single number. The real world is messy. Expect a factor‑of‑two uncertainty on crater diameter.
For very large craters (>~300 km), basin physics becomes important. The single‑law approximation degrades; use with caution.
Target Earth; rho_t ≈ 2.7 g/cm^3. Set rho_i = 3.0 g/cm^3, D_imp = 12 km, v = 20 km/s, theta = 60° from horizontal, K = 13.0. The model returns a final crater diameter of about 190 km by construction, matching Chicxulub‑scale constraints.
This exercise introduces a minimal model for the rise of atmospheric oxygen. The system is expressed in units of PAL (Present Atmospheric Level, where 1 PAL ≈ 21% O2 by volume).
The governing equation is:
dx/dt = P_photo − k_total * x
where:
- x(t) = O2 in PAL
- P_photo = ε_photo * P0 (oxygen input from photosynthesis, PAL/Myr)
- k_total = k_atm + k_fe (total sink rate, 1/Myr)
Analytic solution (for k_total > 0):
x(t) = x_eq + (x0 − x_eq) * exp(−k_total * t)
with equilibrium value:
x_eq = P_photo / k_total
τ = 1 / k_total
P0 (PAL/Myr): baseline photosynthetic oxygen source.
ε_photo (0–2): efficiency multiplier for photosynthesis.
k_atm (1/Myr): atmospheric sink proportional to O2.
k_fe (1/Myr): geological sink proportional to O2.
k_total: sum of k_atm and k_fe.
x0 (PAL): initial O2 concentration.
years_total: total time of the simulation (Myr).
If P_photo = 0 → no oxygen rise. If P_photo > 0 and k_total > 0 → oxygen increases towards equilibrium x_eq. The timescale τ = 1/k_total determines how quickly equilibrium is approached.
Example: P0 = 0.2, ε_photo = 1.0 → P_photo = 0.2 PAL/Myr. If k_atm = 0.02 and k_fe = 0.08, then k_total = 0.10. Equilibrium: x_eq = 0.2 / 0.10 = 2 PAL (≈ 42% O2). Timescale: τ = 10 Myr.
1. Set ε_photo = 0. What happens to oxygen? Why?
2. Increase k_fe while keeping k_atm constant. How does equilibrium x_eq change?
3. Decrease k_total (e.g. both k_atm and k_fe). How does the timescale τ change?
Students can generate the Colab notebook by pasting the following prompt into ChatGPT:
Create a single Jupyter/Colab cell that implements the simplified OxygenRise model of atmospheric oxygen:
Model:
dx/dt = P_photo − k_total * x
where:
- x(t) = O2 in PAL (Present Atmospheric Level, 1 PAL = 21% O2 by volume)
- P_photo = eps_photo * P0 (oxygen input from photosynthesis, in PAL/Myr)
- k_total = k_atm + k_fe (total sink rate, 1/Myr)
Analytic solution:
x(t) = x_eq + (x0 − x_eq) * exp(−k_total * t), where x_eq = P_photo / k_total if k_total > 0
Interactive sliders:
- years_total in [10, 6000]
- dt in [0.01, 2.0]
- P0 in [0.0, 1.0]
- eps_photo in [0.0, 2.0]
- k_atm in [0.0, 0.5]
- k_fe in [0.0, 0.5]
- x0 in [0.0, 2.0]
- logy (checkbox for log scale on O2)
Outputs:
1. Numeric summary: values of P_photo, k_total, x_eq, τ
2. Plot 1: O2 (% of PAL) vs time, with equilibrium line and markers
3. Plot 2: Fluxes (P_photo vs sink = k_total·x) vs time
Use only numpy, matplotlib, and ipywidgets. Keep the cell self-contained.
This simplified model illustrates how populations migrate between two regions when a climate-controlled corridor opens and closes. The model couples logistic population growth with a climate threshold condition that allows or prevents migration.
O(t), D(t): Population sizes in the Origin and Destination regions (individuals).
O0, D0: Initial populations at time t=0 (individuals).
rO, rD: Intrinsic growth rates (per year, 1/yr). This is the maximum fractional growth rate when population density is very low. For example, rO = 0.02 means a 2% increase per year at small population size.
KO, KD: Carrying capacities (individuals). These are the maximum sustainable populations supported by each region’s environment.
m_base: Baseline migration rate (fraction of Origin per year if corridor is open). For example, m_base = 0.01 means that each year 1% of the Origin population migrates to the Destination if the corridor is open.
M(t): Migration flow (individuals per timestep Δt).
C(t): Climate index (dimensionless). Controls whether the corridor is open.
Defined as: C(t) = C0 + A * sin(2πt/P + φ) + shock(t).
C0: Baseline climate index (dimensionless).
A: Amplitude of the climate oscillation (dimensionless).
P: Period of the oscillation (years). This sets the cycle length.
φ (phi_deg): Phase shift (degrees). Determines the horizontal shift of the sine wave along the time axis.
T(t): Threshold function (dimensionless). The corridor opens only if C(t) ≥ T(t).
Defined as: T(t) = T0 + trend * t.
T0: Baseline threshold value (dimensionless).
trend: Drift of the threshold per year (1/yr). Represents long-term changes, such as tectonic uplift or progressive aridity.
shock(t): Gaussian perturbation applied to the climate index C(t).
shock_year: Center year of the perturbation (years).
shock_width: Width of the Gaussian (years).
shock_amp: Amplitude of the perturbation (dimensionless). Positive values make it easier for the corridor to open; negative values make it harder.
years: Total duration of the simulation (years).
dt: Timestep of the simulation (years).
- Population trajectories O(t) and D(t).
- Climate index C(t) versus threshold T(t), with shaded intervals showing when the corridor is open.
- Migration flow M(t)/dt (individuals per year).
Students can generate the Colab notebook by pasting the following prompt into ChatGPT:
Create a single Jupyter/Colab cell that implements the Climate-Window Migration Model with logistic population growth. Define the following parameters and provide interactive sliders using ipywidgets:
- O0, D0 (initial populations)
- rO, rD (growth rates, 1/yr)
- KO, KD (carrying capacities)
- m_base (baseline migration rate)
- C0, A, P, phi_deg (climate index sinusoid)
- T0, trend (threshold function)
- shock_year, shock_width, shock_amp (shock event)
- years, dt (time and step)
Outputs must include:
1. Population trajectories O(t), D(t).
2. Climate index vs threshold, with shaded regions showing corridor open times.
3. Migration flow (individuals per year).
Use only numpy, matplotlib, and ipywidgets. Keep the cell self-contained.