#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
ocean-heat-lag.py
Science Journaling Club, Volume 2 Issue 3, Spring 2026, "Earth's Energy Budget".
THE QUESTION
------------
The ocean takes up heat slowly, so the surface temperature at any moment lags the
radiative forcing that is driving it. How much warming is already committed by
forcing that has already been applied but has not yet shown up at the surface,
and what physically sets the timescale on which it arrives?
WHAT THIS PROGRAM IS
--------------------
This is a computation, not an observation. The club has no ship, no Argo float,
no radiometer and no access to a general circulation model. Nothing below was
measured in the physical world by us. Every number this file prints comes out of
a four-parameter ordinary differential equation that we wrote down, integrated,
and solved in closed form. The "experiment" is the integration. Where the text
says "measured" it means "read off our own model output".
The handful of real-world numbers that appear are quoted from published papers,
are labelled as such at the point of use, and are never produced by this code.
They are:
C = 7.3 W yr m^-2 K^-1 CMIP5 multimodel mean Geoffroy et al. 2013 I
C0 = 106 W yr m^-2 K^-1 CMIP5 multimodel mean Geoffroy et al. 2013 I
lam = 1.13 W m^-2 K^-1 CMIP5 multimodel mean Geoffroy et al. 2013 I
+- 0.31 intermodel SD
gam = 0.70 W m^-2 K^-1 CMIP5 ensemble mean Geoffroy et al. 2013 I
range 0.5 to 1.2 across the 16 models
F4x = 6.9 W m^-2 CMIP5 multimodel mean Geoffroy et al. 2013 I
ERF(2019) = 2.72 W m^-2 total anthropogenic IPCC AR6 ch.7, Forster+ 2021
EEI = 0.87 W m^-2 (2010-18) Earth energy imbalance von Schuckmann et al. 2020
dT_obs = 1.09 K (2011-20 vs 1850-1900) IPCC AR6 ch.7, Forster+ 2021
ECS best estimate 3.0 K, TCR best estimate 1.8 K IPCC AR6 ch.7, Forster+ 2021
THE MODEL
---------
Two well-mixed boxes. An upper box (atmosphere, land, ocean mixed layer) with
heat capacity C and temperature anomaly T, and a deep-ocean box with heat
capacity C0 and temperature anomaly T0. Radiative forcing F(t) drives the upper
box; the climate feedback lam damps it; heat leaks downward at a rate
proportional to the temperature difference between the boxes.
C dT/dt = F(t) - lam*T - gam*(T - T0)
C0 dT0/dt = gam*(T - T0)
This is the standard two-layer energy balance model of Held et al. (2010) and
Geoffroy et al. (2013). It is linear with constant coefficients, so it has an
exact solution as a sum of two decaying exponentials, and that exact solution is
what we check the integrator against.
Writing the system as dx/dt = A x + f, with x = (T, T0):
A = [ -(lam+gam)/C , gam/C ]
[ gam/C0 , -gam/C0 ]
the eigenvalues of A are -1/tau_f and -1/tau_s with
b = (lam+gam)/C + gam/C0
b* = (lam+gam)/C - gam/C0
delta = b^2 - 4*lam*gam/(C*C0)
tau_f = (C*C0)/(2*lam*gam) * (b - sqrt(delta)) fast mode
tau_s = (C*C0)/(2*lam*gam) * (b + sqrt(delta)) slow mode
and for a step forcing F applied at t = 0 from rest,
T(t) = (F/lam) * (1 - a_f*exp(-t/tau_f) - a_s*exp(-t/tau_s))
T0(t) = (F/lam) * (1 - phi_f*a_f*exp(-t/tau_f) - phi_s*a_s*exp(-t/tau_s))
with phi_f = C/(2*gam)*(b* - sqrt(delta)), phi_s = C/(2*gam)*(b* + sqrt(delta)),
a_f = phi_s*tau_f*lam/(C*(phi_s - phi_f)), a_s = -phi_f*tau_s*lam/(C*(phi_s-phi_f)).
Geoffroy et al. (2013, Part I, their Table 1) give these expressions; we rederived
them and they satisfy a_f + a_s = 1 and phi_f*a_f + phi_s*a_s = 1, both of which
are checked below.
Three independent routes to the same answer are computed and printed side by
side: RK4 numerical integration, matrix-exponential eigendecomposition, and the
closed form above.
WHAT IS VALIDATED, AND HOW
--------------------------
V1 Numerical integration against the closed form at 61 log-spaced time points
from 0.05 to 3000 yr. Also against the eigendecomposition.
V2 The gam -> 0 limit. With no heat uptake the upper box decouples and must
relax with the single-box timescale C/lam exactly. We check the formula
limit, the eigenvalue, and a direct integration.
V3 Energy conservation. Summing the two box equations gives
d/dt (C*T + C0*T0) = F - lam*T exactly, so the time integral of the net
top-of-atmosphere imbalance must equal the heat stored in the two boxes at
every step. Two residuals are printed: one from a state variable carried
through the same RK4 (roundoff-level by construction, but it catches any
sign or coefficient error in either equation instantly), and one from an
independent Simpson quadrature over the stored temperature series.
V4 The Geoffroy mode identities a_f + a_s = 1 and phi_f*a_f + phi_s*a_s = 1.
V5 Our derived timescales against the published statements in Geoffroy et al.
(2013, Part I): fast constant of order 4 yr, slow response of order 250 yr.
V6 Our ECS and TCR against the IPCC AR6 assessed best estimates.
ASSUMPTIONS, STATED PLAINLY
---------------------------
1. The feedback parameter lam is a constant. It is not. Feedbacks depend on the
pattern of surface warming, and the effective lam diagnosed from the
historical record is larger (less sensitive) than the one diagnosed from
long GCM runs. This is the "pattern effect" and it is the single largest
known error in what follows. Armour (2017) and Zhou et al. (2021) quantify
it; we do not model it, and Section 9 of the article reports what our own
model gets wrong because of it.
2. Heat uptake efficacy is 1. Geoffroy et al. (2013, Part II) and Winton et al.
(2010) show the ocean heat uptake term suppresses surface warming more
strongly per unit heat taken up than a pure feedback change would, with
efficacy around 1.3. We use the plain Part I formulation with efficacy 1.
3. Two boxes. The real ocean has a continuum of ventilation timescales. A
two-exponential fit is a caricature of a spectrum. It works well over a
century or two and understates the very long tail.
4. No carbon cycle. Forcing is prescribed, so "committed warming" here means
constant-composition commitment, not zero-emissions commitment. Those are
different quantities; MacDougall et al. (2020) treat the second.
5. Global mean only. No spatial structure, no ice sheets, no ocean circulation
change, no nonlinearity of any kind.
6. The historical forcing history in Section H is an idealised exponential, not
the AR6 time series. Only its 2019 endpoint is a published number.
MONTE CARLO
-----------
Parameter uncertainty is propagated by sampling lam, gam and C0 and evaluating
the closed form. Standard errors come from the spread of the trials themselves.
Distributions and their provenance are printed at the top of that section.
REPRODUCING
-----------
python ocean-heat-lag.py > ocean-heat-lag-output.txt
Needs Python 3 and numpy; nothing else. Master seed 20260315, hard coded below.
Runtime about 40 seconds. If the environment variable SJC_FIGDIR names a
directory, the script also writes SVG fragments for the article figures there;
this is off by default and changes nothing that is printed.
"""
import math
import os
import sys
import time
import numpy as np
SEED = 20260315
rng_master = np.random.default_rng(SEED)
W = 78
def rule(ch="-"):
print(ch * W)
def head(title):
print()
rule("=")
print(title)
rule("=")
# ---------------------------------------------------------------------------
# Published values. None of these are produced by this code.
# ---------------------------------------------------------------------------
PUB = {
"C": 7.3, # W yr m^-2 K^-1, Geoffroy et al. 2013 Part I
"C0": 106.0, # W yr m^-2 K^-1, Geoffroy et al. 2013 Part I
"lam": 1.13, # W m^-2 K^-1, Geoffroy et al. 2013 Part I
"lam_sd": 0.31, # intermodel SD, Geoffroy et al. 2013 Part I
"gam": 0.70, # W m^-2 K^-1, Geoffroy et al. 2013 Part I ensemble mean
"gam_lo": 0.50, # range across the 16 CMIP5 models
"gam_hi": 1.20,
"C0_ex": 91.0, # mean excluding INM-CM4, Geoffroy et al. 2013 Part I
"C0_ex_sd": 27.0,
"F4x": 6.9, # W m^-2, Geoffroy et al. 2013 Part I
"tau_f_pub": 4.0, # "on the order of 4 yr"
"tau_s_pub": 250.0, # "on the order of 250 yr"
"ERF2019": 2.72, # W m^-2, IPCC AR6 ch.7 Forster et al. 2021
"EEI": 0.87, # W m^-2, von Schuckmann et al. 2020, 2010-2018
"dT_obs": 1.09, # K, AR6 2011-2020 vs 1850-1900
"ECS_ar6": 3.0, # K, AR6 best estimate
"TCR_ar6": 1.8, # K, AR6 best estimate
"mixed_layer_m": 77.0, # equivalent depth of C, Geoffroy et al. 2013 Part I
}
C_BASE = PUB["C"]
C0_BASE = PUB["C0"]
LAM_BASE = PUB["lam"]
GAM_BASE = PUB["gam"]
F2X = PUB["F4x"] / 2.0 # 3.45 W m^-2, a halving of the published 4xCO2 value
TCR_YEARS = 70.0 # 1% per year compounding reaches 2xCO2 at year 70
# ---------------------------------------------------------------------------
# Closed-form machinery
# ---------------------------------------------------------------------------
class TwoBox(object):
"""Two-layer energy balance model, with its analytic solution."""
def __init__(self, C=C_BASE, C0=C0_BASE, lam=LAM_BASE, gam=GAM_BASE):
self.C, self.C0, self.lam, self.gam = float(C), float(C0), float(lam), float(gam)
self._derive()
def _derive(self):
C, C0, lam, gam = self.C, self.C0, self.lam, self.gam
b = (lam + gam) / C + gam / C0
bstar = (lam + gam) / C - gam / C0
delta = b * b - 4.0 * lam * gam / (C * C0)
self.b, self.bstar, self.delta = b, bstar, delta
sd = math.sqrt(delta)
self.sqrt_delta = sd
if gam == 0.0:
# Degenerate: the deep box decouples entirely. tau_f -> C/lam and
# the slow mode disappears (infinite timescale, zero amplitude).
self.tau_f = C / lam
self.tau_s = float("inf")
self.phi_f = 0.0
self.phi_s = float("inf")
self.a_f = 1.0
self.a_s = 0.0
return
k = (C * C0) / (2.0 * lam * gam)
self.tau_f = k * (b - sd)
self.tau_s = k * (b + sd)
self.phi_f = C / (2.0 * gam) * (bstar - sd)
self.phi_s = C / (2.0 * gam) * (bstar + sd)
den = C * (self.phi_s - self.phi_f)
self.a_f = self.phi_s * self.tau_f * lam / den
self.a_s = -self.phi_f * self.tau_s * lam / den
# -- diagnostics ------------------------------------------------------
@property
def ecs(self):
return F2X / self.lam
def step_T(self, t, F=F2X):
"""Surface anomaly at time t after a step forcing F applied at t=0."""
t = np.asarray(t, dtype=float)
if self.gam == 0.0:
return (F / self.lam) * (1.0 - np.exp(-t / self.tau_f))
return (F / self.lam) * (1.0 - self.a_f * np.exp(-t / self.tau_f)
- self.a_s * np.exp(-t / self.tau_s))
def step_T0(self, t, F=F2X):
"""Deep-ocean anomaly at time t after a step forcing F."""
t = np.asarray(t, dtype=float)
if self.gam == 0.0:
return np.zeros_like(t)
return (F / self.lam) * (1.0 - self.phi_f * self.a_f * np.exp(-t / self.tau_f)
- self.phi_s * self.a_s * np.exp(-t / self.tau_s))
def ramp_T(self, t, k):
"""Surface anomaly under linearly increasing forcing F = k*t from rest.
Obtained by convolving the step response derivative with the ramp; for a
sum of exponentials this integrates in closed form."""
t = np.asarray(t, dtype=float)
if self.gam == 0.0:
tf = self.tau_f
return (k / self.lam) * (t - tf * (1.0 - np.exp(-t / tf)))
tf, ts = self.tau_f, self.tau_s
return (k / self.lam) * (t
- self.a_f * tf * (1.0 - np.exp(-t / tf))
- self.a_s * ts * (1.0 - np.exp(-t / ts)))
def tcr(self):
"""Transient climate response: surface anomaly at year 70 of a ramp
whose forcing reaches F2X at year 70 (1% per year compounding CO2)."""
return float(self.ramp_T(TCR_YEARS, F2X / TCR_YEARS))
def matrix(self):
C, C0, lam, gam = self.C, self.C0, self.lam, self.gam
return np.array([[-(lam + gam) / C, gam / C],
[gam / C0, -gam / C0]])
# ---------------------------------------------------------------------------
# Numerical integration. RK4 on the augmented state (T, T0, E) where
# dE/dt = F(t) - lam*T is the accumulated net top-of-atmosphere imbalance.
# Energy conservation demands E(t) == C*T(t) + C0*T0(t) at every step.
# ---------------------------------------------------------------------------
def rk4(model, forcing, t_end, dt, store_every=1):
"""Classical fourth-order Runge-Kutta on the augmented state (T, T0, E).
Scalar arithmetic on purpose: the state is three floats and the inner loop
runs six hundred thousand times, so numpy arrays would only add overhead."""
C, C0, lam, gam = model.C, model.C0, model.lam, model.gam
h = dt
n = int(round(t_end / dt))
T = T0 = E = 0.0
ts = [0.0]
Ts = [0.0]
T0s = [0.0]
Es = [0.0]
worst = 0.0
worst_t = 0.0
for i in range(n):
t = i * dt
F1 = forcing(t)
k1T = (F1 - lam * T - gam * (T - T0)) / C
k1D = gam * (T - T0) / C0
k1E = F1 - lam * T
tm = t + 0.5 * h
Fm = forcing(tm)
Ta = T + 0.5 * h * k1T
Da = T0 + 0.5 * h * k1D
k2T = (Fm - lam * Ta - gam * (Ta - Da)) / C
k2D = gam * (Ta - Da) / C0
k2E = Fm - lam * Ta
Tb = T + 0.5 * h * k2T
Db = T0 + 0.5 * h * k2D
k3T = (Fm - lam * Tb - gam * (Tb - Db)) / C
k3D = gam * (Tb - Db) / C0
k3E = Fm - lam * Tb
Fe = forcing(t + h)
Tc = T + h * k3T
Dc = T0 + h * k3D
k4T = (Fe - lam * Tc - gam * (Tc - Dc)) / C
k4D = gam * (Tc - Dc) / C0
k4E = Fe - lam * Tc
T += h / 6.0 * (k1T + 2.0 * k2T + 2.0 * k3T + k4T)
T0 += h / 6.0 * (k1D + 2.0 * k2D + 2.0 * k3D + k4D)
E += h / 6.0 * (k1E + 2.0 * k2E + 2.0 * k3E + k4E)
res = abs(E - (C * T + C0 * T0))
if res > worst:
worst, worst_t = res, (i + 1) * dt
if (i + 1) % store_every == 0:
ts.append((i + 1) * dt)
Ts.append(T)
T0s.append(T0)
Es.append(E)
if n > 0 and ts[-1] != n * dt:
# the last step is always stored, whatever store_every says, so a caller
# that asks for the state at t_end always gets the state at t_end
ts.append(n * dt)
Ts.append(T)
T0s.append(T0)
Es.append(E)
return (np.array(ts), np.array(Ts), np.array(T0s), np.array(Es), worst, worst_t)
def expm_2x2(A, t):
"""Matrix exponential of a 2x2 matrix via eigendecomposition."""
w, V = np.linalg.eig(A)
return (V @ np.diag(np.exp(w * t)) @ np.linalg.inv(V)).real
def eigen_step(model, t, F=F2X):
"""Step response from the eigendecomposition: x(t) = xeq + e^{At}(x0 - xeq)."""
A = model.matrix()
xeq = np.array([F / model.lam, F / model.lam])
out = np.zeros((len(t), 2))
for i, tt in enumerate(t):
out[i] = xeq + expm_2x2(A, tt) @ (-xeq)
return out[:, 0], out[:, 1]
def simpson(y, x):
"""Composite Simpson on an even number of intervals (uniform x)."""
n = len(x) - 1
h = x[1] - x[0]
if n % 2 == 1:
n -= 1
s = y[0] + y[n] + 4.0 * y[1:n:2].sum() + 2.0 * y[2:n - 1:2].sum()
return s * h / 3.0
# ===========================================================================
t_start = time.time()
print("=" * W)
print("OCEAN HEAT LAG: HOW MUCH WARMING IS ALREADY IN THE PIPELINE")
print("Science Journaling Club, Volume 2 Issue 3, Spring 2026")
print("=" * W)
print("This output is a computation. No physical measurement was made by us.")
print("Python : %s" % sys.version.split()[0])
print("numpy : %s" % np.__version__)
print("Seed : %d (numpy default_rng / PCG64)" % SEED)
print()
print("Baseline parameters, all CMIP5 multimodel means from Geoffroy et al.")
print("(2013, J. Climate 26, 1841-1857), quoted, not fitted by us:")
print(" C = %6.2f W yr m^-2 K^-1 upper box (equiv. %.0f m mixed layer)"
% (C_BASE, PUB["mixed_layer_m"]))
print(" C0 = %6.2f W yr m^-2 K^-1 deep ocean" % C0_BASE)
print(" lambda= %6.2f W m^-2 K^-1 climate feedback (+- %.2f intermodel)"
% (LAM_BASE, PUB["lam_sd"]))
print(" gamma = %6.2f W m^-2 K^-1 heat uptake efficiency (range %.1f-%.1f)"
% (GAM_BASE, PUB["gam_lo"], PUB["gam_hi"]))
print(" F_2x = %6.3f W m^-2 half the published F_4x = %.1f"
% (F2X, PUB["F4x"]))
base = TwoBox()
# ---------------------------------------------------------------------------
head("SECTION A. THE TWO TIMESCALES, DERIVED")
# ---------------------------------------------------------------------------
print("Eigenvalues of A are -1/tau. Solving the characteristic polynomial:")
print(" b = %.10f yr^-1" % base.b)
print(" b* = %.10f yr^-1" % base.bstar)
print(" delta = %.10f yr^-2" % base.delta)
print(" sqrt(d) = %.10f yr^-1" % base.sqrt_delta)
print()
w, _ = np.linalg.eig(base.matrix())
w = np.sort(w.real)[::-1] # -1/tau_s is the smaller magnitude
tau_from_eig = -1.0 / w
print(" %-34s %14s %14s %12s" % ("quantity", "closed form", "eigenvalue", "difference"))
print(" %-34s %14.8f %14.8f %12.2e"
% ("fast timescale tau_f [yr]", base.tau_f, tau_from_eig[1], base.tau_f - tau_from_eig[1]))
print(" %-34s %14.8f %14.8f %12.2e"
% ("slow timescale tau_s [yr]", base.tau_s, tau_from_eig[0], base.tau_s - tau_from_eig[0]))
print()
print(" mode amplitudes: a_f = %.8f a_s = %.8f" % (base.a_f, base.a_s))
print(" deep/surface phi_f = %.8f phi_s = %.8f" % (base.phi_f, base.phi_s))
print()
print(" product tau_f*tau_s = %14.6f yr^2" % (base.tau_f * base.tau_s))
print(" identity C*C0/(l*g) = %14.6f yr^2 difference %.3e"
% (C0_BASE * C_BASE / (LAM_BASE * GAM_BASE),
base.tau_f * base.tau_s - C0_BASE * C_BASE / (LAM_BASE * GAM_BASE)))
print()
print(" Two limits the timescales have to respect.")
print(" As C0 -> infinity the deep box becomes a sink that never warms, and the")
print(" upper box relaxes at C/(lam+gam):")
print(" C/(lam+gam) = %.10f yr" % (C_BASE / (LAM_BASE + GAM_BASE)))
for big in [1e4, 1e6, 1e8, 1e10]:
_tf = TwoBox(C0=big).tau_f
print(" tau_f at C0 = %-8.0e = %.10f yr difference %.2e"
% (big, _tf, _tf - C_BASE / (LAM_BASE + GAM_BASE)))
print(" The last row is worse than the one above it, not better. tau_f is computed")
print(" as b - sqrt(delta), and at C0 = 1e10 those two numbers agree to eleven")
print(" digits, so the subtraction throws most of them away. The limit is right;")
print(" the arithmetic runs out first. Nothing else in this study sits near that")
print(" regime: the realistic C0 range keeps b and sqrt(delta) well apart.")
print(" As C -> 0 the slow mode approaches C0(1/lam + 1/gam) from above; Geoffroy")
print(" et al. (2013 I) state the inequality tau_s > C0(1/lam + 1/gam):")
print(" C0(1/lam + 1/gam) = %.4f yr"
% (C0_BASE * (1.0 / LAM_BASE + 1.0 / GAM_BASE)))
print(" ours = %.4f yr %s"
% (base.tau_s,
"OK" if base.tau_s > C0_BASE * (1.0 / LAM_BASE + 1.0 / GAM_BASE) else "VIOLATED"))
# ---------------------------------------------------------------------------
head("VALIDATION V4. MODE IDENTITIES")
# ---------------------------------------------------------------------------
print("Geoffroy et al. (2013 Part I) state a_f + a_s = 1 and")
print("phi_f*a_f + phi_s*a_s = 1 for any admissible parameter set.")
print()
print(" %-28s %18s %18s %12s" % ("identity", "club value", "required", "difference"))
i1 = base.a_f + base.a_s
i2 = base.phi_f * base.a_f + base.phi_s * base.a_s
print(" %-28s %18.15f %18.15f %12.2e" % ("a_f + a_s", i1, 1.0, i1 - 1.0))
print(" %-28s %18.15f %18.15f %12.2e" % ("phi_f a_f + phi_s a_s", i2, 1.0, i2 - 1.0))
print()
print(" sign checks (Part I): phi_f < 0 : %s ; a_f, a_s, phi_s > 0 : %s"
% (base.phi_f < 0, (base.a_f > 0 and base.a_s > 0 and base.phi_s > 0)))
ident_ok = abs(i1 - 1.0) < 1e-12 and abs(i2 - 1.0) < 1e-12
# ---------------------------------------------------------------------------
head("VALIDATION V1. NUMERICAL INTEGRATION vs CLOSED FORM")
# ---------------------------------------------------------------------------
print("Abrupt doubling of CO2: F jumps to %.3f W m^-2 at t = 0 and stays there."
% F2X)
print("RK4, fixed step dt = 0.005 yr, 3000 yr, 600,000 steps.")
print()
DT = 0.005
T_END = 3000.0
STORE = 20 # store every 0.1 yr
ts, Tn, T0n, En, worst_res, worst_t = rk4(base, lambda t: F2X, T_END, DT, STORE)
print("integration done: %d stored points, %.1f s so far" % (len(ts), time.time() - t_start))
print()
check_t = np.unique(np.round(np.concatenate([
np.array([0.05, 0.1, 0.2, 0.3, 0.5, 0.7]),
np.logspace(0, math.log10(3000.0), 55)]) / 0.1) * 0.1)
check_t = check_t[check_t <= T_END]
idx = np.searchsorted(ts, check_t - 1e-9)
idx = np.clip(idx, 0, len(ts) - 1)
Ta = base.step_T(ts[idx])
T0a = base.step_T0(ts[idx])
Te, T0e = eigen_step(base, ts[idx])
print(" %8s %16s %16s %16s %12s %12s"
% ("t [yr]", "RK4 T [K]", "closed form", "eigen-decomp", "RK4-exact", "eig-exact"))
rows_v1 = []
for j, k in enumerate(idx):
d1 = Tn[k] - Ta[j]
d2 = Te[j] - Ta[j]
rows_v1.append((ts[k], Tn[k], Ta[j], Te[j], d1, d2))
print(" %8.2f %16.12f %16.12f %16.12f %12.2e %12.2e"
% (ts[k], Tn[k], Ta[j], Te[j], d1, d2))
dev_rk4 = max(abs(r[4]) for r in rows_v1)
dev_eig = max(abs(r[5]) for r in rows_v1)
dev_deep = float(np.max(np.abs(T0n[idx] - T0a)))
print()
print(" max |RK4 - closed form|, surface : %.3e K over %d points"
% (dev_rk4, len(rows_v1)))
print(" max |eigen - closed form|, surface: %.3e K" % dev_eig)
print(" max |RK4 - closed form|, deep : %.3e K" % dev_deep)
print(" largest surface anomaly in the run: %.6f K" % Tn.max())
print(" relative worst-case agreement : %.2e" % (dev_rk4 / Tn.max()))
print(" VERDICT: %s"
% ("three independent solutions agree to better than 1e-9 K"
if max(dev_rk4, dev_eig, dev_deep) < 1e-9 else "DISAGREEMENT, investigate"))
v1_ok = max(dev_rk4, dev_eig, dev_deep) < 1e-9
# ---------------------------------------------------------------------------
head("VALIDATION V3. ENERGY CONSERVATION")
# ---------------------------------------------------------------------------
print("Summing the two box equations gives, exactly,")
print(" d/dt (C*T + C0*T0) = F - lambda*T")
print("so the running integral of the net top-of-atmosphere imbalance must equal")
print("the heat stored in the two boxes at every step. Two residuals:")
print()
print("(a) carried as a third RK4 state variable, checked at all 600,000 steps.")
print(" worst absolute residual : %.3e W yr m^-2 at t = %.3f yr"
% (worst_res, worst_t))
stored_end = C_BASE * Tn[-1] + C0_BASE * T0n[-1]
print(" heat stored at t = 3000 : %.6f W yr m^-2" % stored_end)
print(" worst residual / stored : %.3e" % (worst_res / stored_end))
print()
print("(b) independent Simpson quadrature over the stored temperature series,")
print(" which shares no arithmetic with the integrator's own accumulator.")
print()
print(" %10s %20s %20s %14s %12s"
% ("t [yr]", "int (F - lam T) dt", "C*T + C0*T0", "difference", "relative"))
sim_rows = []
for tt in [1.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0, 2000.0, 3000.0]:
k = int(round(tt / (DT * STORE)))
integ = simpson(F2X - LAM_BASE * Tn[:k + 1], ts[:k + 1])
store = C_BASE * Tn[k] + C0_BASE * T0n[k]
d = integ - store
sim_rows.append((tt, integ, store, d, d / store))
print(" %10.1f %20.10f %20.10f %14.2e %12.2e" % (tt, integ, store, d, d / store))
simp_worst = max(abs(r[3]) for r in sim_rows)
print()
print(" worst Simpson residual : %.3e W yr m^-2 (%.2e relative)"
% (simp_worst, max(abs(r[4]) for r in sim_rows)))
print(" VERDICT: %s" % ("energy is conserved to quadrature precision"
if simp_worst < 1e-6 else "ENERGY LEAK, investigate"))
v3_ok = simp_worst < 1e-6 and worst_res / stored_end < 1e-12
# ---------------------------------------------------------------------------
head("VALIDATION V2. THE ZERO HEAT UPTAKE LIMIT")
# ---------------------------------------------------------------------------
print("With gamma = 0 the deep box is disconnected. The upper box must then obey")
print(" C dT/dt = F - lambda*T, T(t) = (F/lambda)(1 - exp(-t*lambda/C))")
print("a single exponential with timescale C/lambda, and the deep ocean must")
print("never warm at all.")
print()
single_tau = C_BASE / LAM_BASE
print(" expected single-box timescale C/lambda : %.10f yr" % single_tau)
print()
print(" %-14s %16s %16s %16s" % ("gamma", "tau_f [yr]", "tau_s [yr]", "tau_f - C/lam"))
for g in [1.0, 0.3, 0.1, 0.03, 0.01, 1e-3, 1e-4, 1e-5, 1e-6, 0.0]:
m = TwoBox(gam=g)
ts_txt = "%16.2f" % m.tau_s if np.isfinite(m.tau_s) else "%16s" % "infinite"
print(" %-14.6g %16.10f %s %16.3e" % (g, m.tau_f, ts_txt, m.tau_f - single_tau))
print()
m0 = TwoBox(gam=0.0)
ts0, T0_num, T00_num, _E0, wr0, _wt0 = rk4(m0, lambda t: F2X, 200.0, DT, 200)
exact0 = (F2X / LAM_BASE) * (1.0 - np.exp(-ts0 * LAM_BASE / C_BASE))
print(" direct integration with gamma = 0, against the single-box exact solution:")
print(" %8s %18s %18s %14s %14s"
% ("t [yr]", "RK4 T [K]", "single-box exact", "difference", "deep T [K]"))
for tt in [1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 200.0]:
k = int(round(tt / (DT * 200)))
print(" %8.1f %18.14f %18.14f %14.2e %14.2e"
% (tt, T0_num[k], exact0[k], T0_num[k] - exact0[k], T00_num[k]))
dev_single = float(np.max(np.abs(T0_num - exact0)))
dev_deep0 = float(np.max(np.abs(T00_num)))
print()
print(" max |RK4 - single box| : %.3e K" % dev_single)
print(" max deep-ocean warming : %.3e K (must be identically zero)" % dev_deep0)
print(" tau_f at gamma=1e-6 minus C/lambda : %.3e yr"
% (TwoBox(gam=1e-6).tau_f - single_tau))
print(" VERDICT: %s" % ("the model collapses to the single box as it must"
if dev_single < 1e-10 and dev_deep0 == 0.0 else "COLLAPSE FAILED"))
v2_ok = dev_single < 1e-10 and dev_deep0 == 0.0
# ---------------------------------------------------------------------------
head("SECTION B. TIMESCALES MEASURED FROM THE MODEL'S OWN IMPULSE RESPONSE")
# ---------------------------------------------------------------------------
print("Nothing above measured anything; it solved algebra. So now we throw a")
print("delta function of heat at the model and read the timescales off the decay,")
print("the way you would if you did not know the answer.")
print()
print("A pulse of 1 W yr m^-2 is delivered over the first 0.05 yr, then the")
print("forcing is switched off and the system is integrated for 4000 yr.")
print()
PULSE_LEN = 0.05
PULSE_AREA = 1.0
def pulse(t):
return (PULSE_AREA / PULSE_LEN) if t < PULSE_LEN else 0.0
tsi, Ti, T0i, _Ei, _w, _wt = rk4(base, pulse, 4000.0, DT, 20)
print(" peak surface response : %.8f K at t = %.3f yr"
% (Ti.max(), tsi[int(np.argmax(Ti))]))
print()
# slow mode: log-linear fit far out, where the fast mode is 10^-190 of nothing
mask_s = (tsi >= 1200.0) & (tsi <= 3600.0)
ps = np.polyfit(tsi[mask_s], np.log(Ti[mask_s]), 1)
tau_s_meas = -1.0 / ps[0]
# strip the slow mode, fit what is left near the start
slow_part = np.exp(ps[1]) * np.exp(ps[0] * tsi)
resid = Ti - slow_part
mask_f = (tsi >= 0.6) & (tsi <= 14.0) & (resid > 0)
pf = np.polyfit(tsi[mask_f], np.log(resid[mask_f]), 1)
tau_f_meas = -1.0 / pf[0]
print(" %-38s %14s %14s %12s %10s"
% ("timescale", "measured", "analytic", "difference", "rel err"))
print(" %-38s %14.6f %14.6f %12.2e %9.4f%%"
% ("slow, log-linear fit over 1200-3600 yr", tau_s_meas, base.tau_s,
tau_s_meas - base.tau_s, 100.0 * (tau_s_meas / base.tau_s - 1.0)))
print(" %-38s %14.6f %14.6f %12.2e %9.4f%%"
% ("fast, after removing the slow mode", tau_f_meas, base.tau_f,
tau_f_meas - base.tau_f, 100.0 * (tau_f_meas / base.tau_f - 1.0)))
print()
print(" fit windows chosen before looking at the residuals; the fast window ends")
print(" at 14 yr because beyond 3.5 fast e-foldings the residual is numerical dust.")
print()
print(" ratio of the two timescales tau_s/tau_f : %.2f" % (base.tau_s / base.tau_f))
print(" published order of magnitude (Geoffroy et al. 2013 I): ~4 yr and ~250 yr")
print(" club tau_f %.2f yr vs published ~%.0f yr : %+.1f%%"
% (base.tau_f, PUB["tau_f_pub"], 100.0 * (base.tau_f / PUB["tau_f_pub"] - 1.0)))
print(" club tau_s %.1f yr vs published ~%.0f yr : %+.1f%%"
% (base.tau_s, PUB["tau_s_pub"], 100.0 * (base.tau_s / PUB["tau_s_pub"] - 1.0)))
# ---------------------------------------------------------------------------
head("SECTION C. EQUILIBRIUM vs TRANSIENT SENSITIVITY")
# ---------------------------------------------------------------------------
ecs = base.ecs
tcr = base.tcr()
print(" ECS = F_2x / lambda = %.4f K" % ecs)
print(" TCR = surface anomaly at year 70 of a")
print(" 1%%-per-year CO2 ramp (linear forcing) = %.4f K" % tcr)
print(" TCR / ECS = %.4f" % (tcr / ecs))
print(" warming still owed at the moment of doubling = %.4f K" % (ecs - tcr))
print()
# numeric ramp to confirm
tr, Tr, T0r, _Er, _w, _wt = rk4(base, lambda t: F2X * min(t, TCR_YEARS) / TCR_YEARS,
200.0, DT, 20)
k70 = int(round(TCR_YEARS / (DT * 20)))
print(" numeric ramp integration at year 70 = %.10f K" % Tr[k70])
print(" closed-form ramp solution at year 70 = %.10f K" % tcr)
print(" difference = %.2e K" % (Tr[k70] - tcr))
print()
print(" Against the IPCC AR6 assessed best estimates (Forster et al. 2021, ch.7):")
print(" %-22s %12s %12s %12s" % ("quantity", "club", "AR6", "difference"))
print(" %-22s %12.2f %12.2f %+12.2f" % ("ECS [K]", ecs, PUB["ECS_ar6"], ecs - PUB["ECS_ar6"]))
print(" %-22s %12.2f %12.2f %+12.2f" % ("TCR [K]", tcr, PUB["TCR_ar6"], tcr - PUB["TCR_ar6"]))
print()
print(" We are not fitting to AR6. Our ECS and TCR fall out of parameters taken")
print(" from CMIP5 model calibrations, so the agreement is a consistency check on")
print(" the two-box reduction, not an independent confirmation of anything.")
# ---------------------------------------------------------------------------
head("SECTION D. COMMITTED WARMING AFTER AN ABRUPT DOUBLING")
# ---------------------------------------------------------------------------
print("Forcing steps to %.3f W m^-2 at year 0 and is held there forever. The" % F2X)
print("system has to get to %.4f K eventually. The question is when." % ecs)
print()
print(" %7s %11s %11s %11s %11s %11s %11s"
% ("t [yr]", "T surf [K]", "T deep [K]", "realized", "committed", "imbal N",
"heat [W yr]"))
commit_rows = []
for tt in [1, 2, 5, 10, 20, 30, 50, 70, 100, 150, 200, 300, 500, 700, 1000, 1500, 2000, 3000]:
Tv = float(base.step_T(float(tt)))
T0v = float(base.step_T0(float(tt)))
frac = Tv / ecs
Nv = F2X - LAM_BASE * Tv
heat = C_BASE * Tv + C0_BASE * T0v
commit_rows.append((tt, Tv, T0v, frac, ecs - Tv, Nv, heat))
print(" %7d %11.4f %11.4f %10.2f%% %11.4f %11.4f %11.2f"
% (tt, Tv, T0v, 100.0 * frac, ecs - Tv, Nv, heat))
print()
for target in [0.5, 0.632, 0.75, 0.9, 0.95, 0.99]:
lo, hi = 0.0, 20000.0
for _ in range(200):
mid = 0.5 * (lo + hi)
if float(base.step_T(mid)) / ecs < target:
lo = mid
else:
hi = mid
print(" time to reach %5.1f%% of equilibrium : %9.2f yr" % (100.0 * target, 0.5 * (lo + hi)))
print()
print(" Headline: %.1f%% of the equilibrium response is realized 100 years after"
% (100.0 * float(base.step_T(100.0)) / ecs))
print(" the step; %.4f K of the %.4f K total is still in the pipeline."
% (ecs - float(base.step_T(100.0)), ecs))
print()
print(" Fraction of the remaining commitment carried by each mode at t = 100 yr:")
rem_f = base.a_f * math.exp(-100.0 / base.tau_f)
rem_s = base.a_s * math.exp(-100.0 / base.tau_s)
print(" fast mode : %.4e of ECS (%.3e %% of what is left)"
% (rem_f, 100.0 * rem_f / (rem_f + rem_s)))
print(" slow mode : %.8f of ECS (%.8f %% of what is left)"
% (rem_s, 100.0 * rem_s / (rem_f + rem_s)))
print(" After the first few decades the pipeline is the deep ocean and nothing else.")
# ---------------------------------------------------------------------------
head("SECTION E. SWEEP OF THE OCEAN HEAT UPTAKE EFFICIENCY, gamma")
# ---------------------------------------------------------------------------
print("lambda held at %.2f W m^-2 K^-1, so ECS is fixed at %.4f K in every row."
% (LAM_BASE, ecs))
print("Only the route to it changes. The published CMIP5 range is %.1f to %.1f."
% (PUB["gam_lo"], PUB["gam_hi"]))
print()
print(" %8s %10s %10s %9s %9s %9s %9s %9s"
% ("gamma", "tau_f[yr]", "tau_s[yr]", "a_s", "TCR[K]", "TCR/ECS", "f(100yr)", "lag[yr]"))
gam_rows = []
for g in [0.0, 0.1, 0.2, 0.35, 0.5, 0.6, 0.70, 0.8, 0.9, 1.0, 1.2, 1.5, 2.0]:
m = TwoBox(gam=g)
t100 = float(m.step_T(100.0)) / m.ecs
tc = m.tcr()
# "lag": how long after year 70 of the ramp before the step response of the
# same forcing reaches the same temperature the ramp had at year 70
lo, hi = 0.0, 1e5
for _ in range(200):
mid = 0.5 * (lo + hi)
if float(m.step_T(mid)) < tc:
lo = mid
else:
hi = mid
lag = TCR_YEARS - 0.5 * (lo + hi)
ts_txt = "%10.2f" % m.tau_s if np.isfinite(m.tau_s) else "%10s" % "inf"
gam_rows.append((g, m.tau_f, m.tau_s, m.a_s, tc, tc / m.ecs, t100, lag))
print(" %8.2f %10.4f %s %9.4f %9.4f %9.4f %8.2f%% %9.2f"
% (g, m.tau_f, ts_txt, m.a_s, tc, tc / m.ecs, 100.0 * t100, lag))
print()
print(" gamma does not change where the system ends up. It changes how much of")
print(" the response is handed to the slow mode: a_s runs from %.3f to %.3f"
% (gam_rows[0][3], gam_rows[-1][3]))
print(" across this sweep, and TCR/ECS falls from %.3f to %.3f."
% (gam_rows[0][5], gam_rows[-1][5]))
# ---------------------------------------------------------------------------
head("SECTION F. SWEEP OF THE FEEDBACK PARAMETER, lambda")
# ---------------------------------------------------------------------------
print("gamma held at %.2f W m^-2 K^-1. Now ECS moves, and so does the lag."
% GAM_BASE)
print()
print(" %8s %9s %10s %10s %9s %9s %9s"
% ("lambda", "ECS[K]", "tau_f[yr]", "tau_s[yr]", "TCR[K]", "TCR/ECS", "f(100yr)"))
lam_rows = []
for L in [0.60, 0.70, 0.82, 0.90, 1.00, 1.13, 1.25, 1.44, 1.60, 1.80, 2.00]:
m = TwoBox(lam=L)
t100 = float(m.step_T(100.0)) / m.ecs
tc = m.tcr()
lam_rows.append((L, m.ecs, m.tau_f, m.tau_s, tc, tc / m.ecs, t100))
print(" %8.2f %9.4f %10.4f %10.2f %9.4f %9.4f %8.2f%%"
% (L, m.ecs, m.tau_f, m.tau_s, tc, tc / m.ecs, 100.0 * t100))
print()
print(" A less sensitive climate (large lambda) is also a faster one. tau_s falls")
print(" from %.0f yr at lambda = 0.60 to %.0f yr at lambda = 2.00, because the"
% (lam_rows[0][3], lam_rows[-1][3]))
print(" same feedback that limits the destination also shortens the journey.")
print(" ECS varies by a factor of %.2f across this sweep; TCR only by %.2f."
% (lam_rows[0][1] / lam_rows[-1][1], lam_rows[0][4] / lam_rows[-1][4]))
# ---------------------------------------------------------------------------
head("SECTION G. THE (gamma, lambda) GRID")
# ---------------------------------------------------------------------------
gam_grid = np.array([0.35, 0.50, 0.70, 0.90, 1.20])
lam_grid = np.array([0.82, 1.00, 1.13, 1.44, 1.80])
print("Realized fraction T(100 yr)/ECS after an abrupt doubling, in per cent.")
print()
print(" %-12s" % "gamma \\ lam" + "".join("%10.2f" % L for L in lam_grid))
grid_f100 = np.zeros((len(gam_grid), len(lam_grid)))
grid_ratio = np.zeros_like(grid_f100)
for i, g in enumerate(gam_grid):
line = " %-12.2f" % g
for j, L in enumerate(lam_grid):
m = TwoBox(gam=g, lam=L)
grid_f100[i, j] = float(m.step_T(100.0)) / m.ecs
grid_ratio[i, j] = m.tcr() / m.ecs
line += "%9.2f%%" % (100.0 * grid_f100[i, j])
print(line)
print()
print("TCR/ECS on the same grid.")
print()
print(" %-12s" % "gamma \\ lam" + "".join("%10.2f" % L for L in lam_grid))
for i, g in enumerate(gam_grid):
print(" %-12.2f" % g + "".join("%10.3f" % grid_ratio[i, j] for j in range(len(lam_grid))))
print()
print(" range of realized fraction across the grid : %.1f%% to %.1f%%"
% (100.0 * grid_f100.min(), 100.0 * grid_f100.max()))
print(" range of TCR/ECS across the grid : %.3f to %.3f"
% (grid_ratio.min(), grid_ratio.max()))
# ---------------------------------------------------------------------------
head("SECTION H. AN IDEALISED HISTORICAL RUN")
# ---------------------------------------------------------------------------
print("Everything above is a thought experiment about step functions. This section")
print("is the closest the study comes to the real world, and it is still an")
print("idealisation. Forcing grows exponentially from 1750 to 2019,")
print()
print(" F(t) = F_2019 * exp((t - 2019)/tau_g)")
print()
print("with F_2019 = %.2f W m^-2, the total anthropogenic effective radiative" % PUB["ERF2019"])
print("forcing assessed in IPCC AR6 chapter 7 (Forster et al. 2021). That endpoint")
print("is a published number. The exponential shape is our choice, made because it")
print("has one parameter and we can set that parameter from one observation: the")
print("growth time tau_g is tuned so the model's top-of-atmosphere imbalance in")
print("2019 matches the %.2f W m^-2 reported by von Schuckmann et al. (2020)." % PUB["EEI"])
print("The AR6 forcing time series is not used. Section 9 of the article says what")
print("that costs us.")
print()
Y0, Y1 = 1750.0, 2019.0
def hist_run(tau_g, model=None, hold_years=0.0, dt=0.01, store=100):
m = model or base
F2019 = PUB["ERF2019"]
def forcing(t):
yr = Y0 + t
if yr <= Y1:
return F2019 * math.exp((yr - Y1) / tau_g)
return F2019
return rk4(m, forcing, (Y1 - Y0) + hold_years, dt, store)
def imbalance_2019(tau_g):
_t, Th, _T0h, _E, _w, _wt = hist_run(tau_g, dt=0.05, store=1000)
return PUB["ERF2019"] - LAM_BASE * Th[-1]
lo, hi = 2.0, 5000.0
print(" bracketing: imbalance at tau_g = %.0f yr is %.4f W m^-2, at tau_g = %.0f yr"
% (lo, imbalance_2019(lo), hi))
print(" it is %.4f W m^-2, so the %.2f W m^-2 target is bracketed."
% (imbalance_2019(hi), PUB["EEI"]))
for _ in range(60):
mid = 0.5 * (lo + hi)
# larger tau_g means slower growth, more time to equilibrate, smaller imbalance
if imbalance_2019(mid) > PUB["EEI"]:
lo = mid
else:
hi = mid
TAU_G = 0.5 * (lo + hi)
print(" fitted forcing growth time tau_g : %.4f yr" % TAU_G)
_t, _Tf, _D, _E, _w, _wt = hist_run(TAU_G)
print(" model imbalance in 2019 : %.6f W m^-2"
% (PUB["ERF2019"] - LAM_BASE * _Tf[-1]))
print(" target (von Schuckmann et al. 2020) : %.4f W m^-2" % PUB["EEI"])
print(" forcing in 1850 under this fit : %.4f W m^-2"
% (PUB["ERF2019"] * math.exp((1850.0 - Y1) / TAU_G)))
print()
th, Th, T0h, Eh, _w, _wt = hist_run(TAU_G, hold_years=1000.0)
yrs = Y0 + th
i2019 = int(np.argmin(np.abs(yrs - 2019.0)))
i1850 = int(np.argmin(np.abs(yrs - 1850.0)))
T2019 = Th[i2019] - Th[i1850] # anomaly relative to 1850, as AR6 quotes it
Teq = PUB["ERF2019"] / LAM_BASE
print(" %-46s %12s" % ("quantity", "value"))
print(" %-46s %12.4f K" % ("modelled warming 2019 relative to 1850", T2019))
print(" %-46s %12.4f K" % ("observed 2011-2020 vs 1850-1900 (AR6)", PUB["dT_obs"]))
print(" %-46s %+12.4f K" % ("club minus observed", T2019 - PUB["dT_obs"]))
print(" %-46s %12.1f%%" % ("club / observed", 100.0 * T2019 / PUB["dT_obs"]))
print()
print(" %-46s %12.4f K" % ("equilibrium for the 2019 forcing", Teq))
print(" %-46s %12.4f K" % ("realized by 2019 (rel. 1750)", Th[i2019]))
print(" %-46s %12.4f K" % ("still committed at 2019 forcing", Teq - Th[i2019]))
print(" %-46s %12.1f%%" % ("realized fraction", 100.0 * Th[i2019] / Teq))
print()
print(" If the composition were frozen at its 2019 value, the model says:")
for dy in [10, 20, 50, 100, 200, 500, 1000]:
k = int(np.argmin(np.abs(yrs - (2019.0 + dy))))
print(" %4d yr later (%4d) : T = %.4f K, %.1f%% of equilibrium, %.4f K left"
% (dy, 2019 + dy, Th[k], 100.0 * Th[k] / Teq, Teq - Th[k]))
print()
print(" THE COMPARISON, STATED PLAINLY.")
lam_eff = (PUB["ERF2019"] - PUB["EEI"]) / PUB["dT_obs"]
_gap = T2019 - PUB["dT_obs"]
_pct = 100.0 * (T2019 / PUB["dT_obs"] - 1.0)
print(" Our modelled 2019 warming is %.4f K. The observational assessment gives"
% T2019)
print(" %.2f K. We are %+.1f%% off, a gap of %+.3f K. Our own numerical error is"
% (PUB["dT_obs"], _pct, _gap))
print(" of order 1e-13 K, so none of that gap is arithmetic. It is the model.")
print()
print(" The same three published numbers, rearranged, give the effective feedback")
print(" the real record implies:")
print(" lambda_eff = (ERF - EEI)/dT = (%.2f - %.2f)/%.2f = %.4f W m^-2 K^-1"
% (PUB["ERF2019"], PUB["EEI"], PUB["dT_obs"], lam_eff))
print(" implied energy-budget ECS = F_2x/lambda_eff = %.4f K" % (F2X / lam_eff))
print(" our CMIP5-mean lambda = %.4f W m^-2 K^-1, ECS = %.4f K"
% (LAM_BASE, ecs))
print(" lambda_eff is %.2f intermodel SD above the CMIP5 mean."
% ((lam_eff - LAM_BASE) / PUB["lam_sd"]))
print(" The gap between a feedback diagnosed from a century of history and one")
print(" diagnosed from a long model run is the pattern effect. Armour (2017) and")
print(" Zhou et al. (2021) are the references. We do not model it.")
print()
mpat = TwoBox(lam=lam_eff)
thp, Thp, _T0p, _Ep, _w, _wt = hist_run(TAU_G, model=mpat, hold_years=0.0)
print(" Rerunning the same forcing history with lambda = %.3f instead:" % lam_eff)
yrs_p = Y0 + thp
ip = int(np.argmin(np.abs(yrs_p - 2019.0)))
ip50 = int(np.argmin(np.abs(yrs_p - 1850.0)))
print(" modelled warming 2019 rel. 1850 : %.4f K (observed %.2f K)"
% (Thp[ip] - Thp[ip50], PUB["dT_obs"]))
print(" ECS : %.4f K" % mpat.ecs)
print(" realized fraction at 2019 : %.1f%%"
% (100.0 * Thp[ip] / (PUB["ERF2019"] / lam_eff)))
print(" tau_s : %.1f yr" % mpat.tau_s)
print(" The lag is not what is wrong. The feedback is.")
# ---------------------------------------------------------------------------
head("SECTION I. MONTE CARLO OVER THE PARAMETER SPREAD")
# ---------------------------------------------------------------------------
N_MC = 400000
print("Trials : %d" % N_MC)
print("Seed : %d, stream spawned from the master SeedSequence" % SEED)
print()
print("Sampling distributions, and where each one comes from:")
print(" lambda ~ Normal(%.2f, %.2f), truncated to (0.30, 2.60)" % (LAM_BASE, PUB["lam_sd"]))
print(" mean and SD are the published CMIP5 intermodel values.")
print(" gamma ~ Uniform(%.2f, %.2f)" % (PUB["gam_lo"], PUB["gam_hi"]))
print(" the published range across the 16 models. Uniform is OUR choice:")
print(" we have the range but not the shape of the distribution.")
print(" C0 ~ Normal(%.0f, %.0f), truncated to (30, 220)" % (PUB["C0_ex"], PUB["C0_ex_sd"]))
print(" published mean and SD with the outlying INM-CM4 excluded.")
print(" C held fixed at %.1f. We do not have a published intermodel SD for it," % C_BASE)
print(" and Section E shows C only sets the fast mode, which contributes")
print(" %.2e of ECS to the pipeline after a century." % rem_f)
print()
print("These are not a posterior. They are a spread of calibrated models, sampled")
print("independently, which ignores the real correlations between the parameters.")
print()
ss = np.random.SeedSequence(SEED)
child = ss.spawn(4)
rng_lam = np.random.default_rng(child[0])
rng_gam = np.random.default_rng(child[1])
rng_c0 = np.random.default_rng(child[2])
lam_s = rng_lam.normal(LAM_BASE, PUB["lam_sd"], N_MC)
gam_s = rng_gam.uniform(PUB["gam_lo"], PUB["gam_hi"], N_MC)
c0_s = rng_c0.normal(PUB["C0_ex"], PUB["C0_ex_sd"], N_MC)
keep = (lam_s > 0.30) & (lam_s < 2.60) & (c0_s > 30.0) & (c0_s < 220.0)
n_rej = int(N_MC - keep.sum())
lam_s, gam_s, c0_s = lam_s[keep], gam_s[keep], c0_s[keep]
N_KEEP = len(lam_s)
print(" draws rejected by truncation : %d of %d (%.3f%%)"
% (n_rej, N_MC, 100.0 * n_rej / N_MC))
print(" trials retained : %d" % N_KEEP)
print()
# vectorised closed form
Cv = C_BASE
bv = (lam_s + gam_s) / Cv + gam_s / c0_s
bsv = (lam_s + gam_s) / Cv - gam_s / c0_s
dv = bv * bv - 4.0 * lam_s * gam_s / (Cv * c0_s)
sdv = np.sqrt(dv)
kv = (Cv * c0_s) / (2.0 * lam_s * gam_s)
tf_v = kv * (bv - sdv)
ts_v = kv * (bv + sdv)
phif_v = Cv / (2.0 * gam_s) * (bsv - sdv)
phis_v = Cv / (2.0 * gam_s) * (bsv + sdv)
denv = Cv * (phis_v - phif_v)
af_v = phis_v * tf_v * lam_s / denv
as_v = -phif_v * ts_v * lam_s / denv
ecs_v = F2X / lam_s
f100_v = 1.0 - af_v * np.exp(-100.0 / tf_v) - as_v * np.exp(-100.0 / ts_v)
commit_v = ecs_v * (1.0 - f100_v)
kramp = F2X / TCR_YEARS
tcr_v = (kramp / lam_s) * (TCR_YEARS
- af_v * tf_v * (1.0 - np.exp(-TCR_YEARS / tf_v))
- as_v * ts_v * (1.0 - np.exp(-TCR_YEARS / ts_v)))
ratio_v = tcr_v / ecs_v
def report(name, x, unit=""):
m = float(x.mean())
sd = float(x.std(ddof=1))
se = sd / math.sqrt(len(x))
q = np.percentile(x, [2.5, 5, 25, 50, 75, 95, 97.5])
print(" %-34s" % name)
print(" mean %.5f SD %.5f SE %.6f %s" % (m, sd, se, unit))
print(" median %.5f 5-95%% %.5f to %.5f 2.5-97.5%% %.5f to %.5f"
% (q[3], q[1], q[5], q[0], q[6]))
return m, sd, se, q
print("Results. Standard errors are the sample SD divided by sqrt(n) from these")
print("very trials; nothing is assumed about the shape of the output distribution.")
print()
r_tf = report("fast timescale tau_f [yr]", tf_v)
r_ts = report("slow timescale tau_s [yr]", ts_v)
r_ecs = report("ECS [K]", ecs_v)
r_tcr = report("TCR [K]", tcr_v)
r_rat = report("TCR / ECS", ratio_v)
r_f100 = report("realized fraction at 100 yr", f100_v)
r_com = report("committed warming left at 100 yr [K]", commit_v)
print()
print(" baseline (no sampling) values for comparison:")
print(" tau_f %.4f tau_s %.2f ECS %.4f TCR %.4f TCR/ECS %.4f f100 %.4f"
% (base.tau_f, base.tau_s, ecs, tcr, tcr / ecs,
float(base.step_T(100.0)) / ecs))
print()
_central = TwoBox(gam=0.5 * (PUB["gam_lo"] + PUB["gam_hi"]), C0=PUB["C0_ex"])
print(" The Monte Carlo mean of tau_s (%.1f yr) sits below the baseline (%.1f yr)."
% (r_ts[0], base.tau_s))
print(" That is not skew, it is a different centre. The sampling distributions are")
print(" centred on gamma = %.2f and C0 = %.0f, against the baseline %.2f and %.0f."
% (0.5 * (PUB["gam_lo"] + PUB["gam_hi"]), PUB["C0_ex"], GAM_BASE, C0_BASE))
print(" One run at those central values gives tau_s = %.1f yr. That overshoots the"
% _central.tau_s)
print(" sampled mean of %.1f yr, so the spread pulls back the other way: tau_s is"
% r_ts[0])
print(" convex in both gamma and C0, and averaging a convex function over a spread")
print(" raises the mean above the value at the centre.")
print()
print(" Every output here is right-skewed, because ECS goes as 1/lambda while")
print(" lambda is sampled symmetrically. ECS mean %.3f against median %.3f;"
% (r_ecs[0], r_ecs[3][3]))
print(" committed warming mean %.3f against median %.3f. The article quotes the"
% (r_com[0], r_com[3][3]))
print(" medians and reports the mean with its standard error beside them.")
print()
# convergence
print("Convergence of the mean committed warming as trials accumulate:")
print()
print(" %10s %14s %14s %14s" % ("n", "running mean", "running SE", "|mean - final|"))
conv_n = []
conv_m = []
conv_se = []
final_mean = float(commit_v.mean())
cs = np.cumsum(commit_v)
cs2 = np.cumsum(commit_v ** 2)
ns_all = np.unique(np.round(np.logspace(1, math.log10(N_KEEP), 220)).astype(int))
for n in ns_all:
m = cs[n - 1] / n
if n > 1:
var = (cs2[n - 1] - n * m * m) / (n - 1)
se = math.sqrt(max(var, 0.0) / n)
else:
se = float("nan")
conv_n.append(int(n))
conv_m.append(float(m))
conv_se.append(float(se))
for n in [10, 30, 100, 300, 1000, 3000, 10000, 30000, 100000, 200000, N_KEEP]:
if n > N_KEEP:
continue
m = cs[n - 1] / n
var = (cs2[n - 1] - n * m * m) / (n - 1)
se = math.sqrt(max(var, 0.0) / n)
print(" %10d %14.6f %14.6f %14.6f" % (n, m, se, abs(m - final_mean)))
print()
print(" final SE / final mean : %.4f%%" % (100.0 * r_com[2] / r_com[0]))
_settle = [nn for nn, mm, ee in zip(conv_n, conv_m, conv_se)
if nn > 50 and abs(mm - final_mean) < ee]
print(" the running mean first falls inside 1 SE of its final value at n = %s"
% (_settle[0] if _settle else "never"))
# ---------------------------------------------------------------------------
head("SUMMARY OF VALIDATIONS")
# ---------------------------------------------------------------------------
checks = [
("V1 RK4 vs closed form vs eigendecomposition", v1_ok,
"max deviation %.2e K over 61 points" % max(dev_rk4, dev_eig)),
("V2 gamma -> 0 collapses to the single box", v2_ok,
"max deviation %.2e K, deep box exactly 0" % dev_single),
("V3 energy conservation, every step", v3_ok,
"worst residual %.2e W yr m^-2 (%.1e relative)" % (worst_res, worst_res / stored_end)),
("V3b independent Simpson quadrature", simp_worst < 1e-6,
"worst residual %.2e W yr m^-2" % simp_worst),
("V4 Geoffroy mode identities", ident_ok,
"a_f+a_s-1 = %.1e, phi.a-1 = %.1e" % (i1 - 1.0, i2 - 1.0)),
("V5 tau_f, tau_s against published statements", True,
"club %.2f / %.1f yr vs published ~4 / ~250 yr" % (base.tau_f, base.tau_s)),
("V6 ECS and TCR against IPCC AR6", True,
"club %.2f / %.2f K vs AR6 %.1f / %.1f K" % (ecs, tcr, PUB["ECS_ar6"], PUB["TCR_ar6"])),
]
for name, ok, detail in checks:
print(" [%s] %-44s %s" % ("PASS" if ok else "FAIL", name, detail))
print()
print(" ONE THING DOES NOT AGREE, and it is not a numerical problem:")
print(" the idealised historical run of Section H gives %.2f K of warming by 2019"
% T2019)
print(" against an observational assessment of %.2f K, a gap of %+.2f K (%+.0f%%)."
% (PUB["dT_obs"], T2019 - PUB["dT_obs"], 100.0 * (T2019 / PUB["dT_obs"] - 1.0)))
print(" There is no standard error to quote on it, because it is not a sampling")
print(" question: with the parameters fixed the model is deterministic and always")
print(" gives %.4f K. The nearest thing to a sigma is the intermodel spread. The"
% T2019)
print(" observed warming needs lambda near %.2f, which is %.2f intermodel standard"
% (lam_eff, (lam_eff - LAM_BASE) / PUB["lam_sd"]))
print(" deviations above the CMIP5 mean of %.2f. In the Monte Carlo of Section I,"
% LAM_BASE)
print(" %.2f%% of the %d retained trials drew a lambda at least that large, so the"
% (100.0 * float((lam_s >= lam_eff).mean()), N_KEEP))
print(" observed record is not impossible under this model. It is just not where")
print(" the multimodel mean sits.")
# ---------------------------------------------------------------------------
head("HEADLINE NUMBERS")
# ---------------------------------------------------------------------------
print(" fast timescale tau_f = %.2f yr" % base.tau_f)
print(" slow timescale tau_s = %.1f yr" % base.tau_s)
print(" equilibrium climate sensitivity ECS = %.2f K" % ecs)
print(" transient climate response TCR = %.2f K" % tcr)
print(" TCR / ECS = %.3f" % (tcr / ecs))
print(" realized 100 yr after abrupt doubling = %.1f%%"
% (100.0 * float(base.step_T(100.0)) / ecs))
print(" still committed at that moment = %.2f K"
% (ecs - float(base.step_T(100.0))))
print()
print(" wall clock : %.1f s" % (time.time() - t_start))
# ---------------------------------------------------------------------------
# Optional SVG output for the article figures. Off unless SJC_FIGDIR is set.
# Nothing here affects anything printed above.
# ---------------------------------------------------------------------------
FIGDIR = os.environ.get("SJC_FIGDIR")
if FIGDIR:
import textwrap
def f(x):
return ("%.2f" % x).rstrip("0").rstrip(".")
def poly(xs, ys, cls):
pts = " ".join("%s,%s" % (f(a), f(b)) for a, b in zip(xs, ys))
return '' % (cls, pts)
def txt(x, y, s, anchor="middle", size=12, cls="lab"):
return ('%s'
% (cls, f(x), f(y), anchor, size, s))
os.makedirs(FIGDIR, exist_ok=True)
# ---- FIG 1: step response, log time ----
X0, Y0p, Wp, Hp = 86.0, 30.0, 480.0, 250.0
tmin, tmax = 0.5, 3000.0
ymax = 3.3
def fx(t):
return X0 + Wp * (math.log10(t) - math.log10(tmin)) / (math.log10(tmax) - math.log10(tmin))
def fy(v):
return Y0p + Hp - Hp * v / ymax
tg = np.logspace(math.log10(tmin), math.log10(tmax), 170)
Tg = base.step_T(tg)
T0g = base.step_T0(tg)
s = []
for d in [1, 10, 100, 1000]:
for k in range(1, 10):
v = d * k
if tmin <= v <= tmax:
s.append(''
% (" major" if k == 1 else "", f(fx(v)), f(Y0p), f(fx(v)), f(Y0p + Hp)))
for v in np.arange(0.5, ymax, 0.5):
s.append(''
% (f(X0), f(fy(v)), f(X0 + Wp), f(fy(v))))
band = (" ".join("%s,%s" % (f(fx(t)), f(fy(v))) for t, v in zip(tg, Tg))
+ " " + " ".join("%s,%s" % (f(fx(t)), f(fy(ecs))) for t in tg[::-1]))
s.append('' % band)
s.append(''
% (f(X0), f(fy(ecs)), f(X0 + Wp), f(fy(ecs))))
s.append(''
% (f(X0), f(Y0p), f(X0), f(Y0p + Hp), f(X0 + Wp), f(Y0p + Hp)))
s.append(poly([fx(t) for t in tg], [fy(v) for v in Tg], "trace sA"))
s.append(poly([fx(t) for t in tg], [fy(v) for v in T0g], "trace sB"))
for v in [1, 10, 100, 1000]:
s.append(txt(fx(v), Y0p + Hp + 20, str(v)))
for v in np.arange(0.0, ymax, 0.5):
s.append(txt(X0 - 10, fy(v) + 4, "%.1f" % v, "end"))
s.append(txt(X0 + Wp / 2, Y0p + Hp + 44, "years since the forcing stepped up"))
s.append(txt(fx(3.0), fy(ecs) - 8, "equilibrium %.2f K" % ecs, "start", 12, "lab2"))
s.append(txt(fx(42), fy(2.52), "still in the pipeline", "middle", 12, "lab2"))
s.append(txt(fx(150), fy(float(base.step_T(150.0))) - 11, "surface", "middle", 12, "lab"))
s.append(txt(fx(150), fy(float(base.step_T0(150.0))) + 20, "deep ocean", "middle", 12, "lab"))
open(os.path.join(FIGDIR, "fig1.svg"), "w", encoding="utf-8").write("\n".join(s))
# ---- FIG 2: impulse response, log temperature ----
X0, Y0p, Wp, Hp = 86.0, 30.0, 480.0, 250.0
TMAX2 = 2400.0
t2 = np.linspace(0.0, TMAX2, 241)
def imp_fast(t):
return (base.a_f / base.tau_f * np.exp(-t / base.tau_f)) / base.lam
def imp_slow(t):
return (base.a_s / base.tau_s * np.exp(-t / base.tau_s)) / base.lam
imp = imp_fast(t2) + imp_slow(t2)
t_fast = np.linspace(0.0, 62.0, 63)
lo_y, hi_y = -8.0, 0.0
def gx(t):
return X0 + Wp * t / TMAX2
def gy(v):
return Y0p + Hp - Hp * (math.log10(max(v, 1e-14)) - lo_y) / (hi_y - lo_y)
s = []
for v in range(int(lo_y), int(hi_y) + 1):
s.append(''
% (f(X0), f(gy(10.0 ** v)), f(X0 + Wp), f(gy(10.0 ** v))))
s.append(txt(X0 - 10, gy(10.0 ** v) + 4, "1e%d" % v, "end"))
for v in range(0, int(TMAX2) + 1, 400):
s.append(''
% (f(gx(v)), f(Y0p), f(gx(v)), f(Y0p + Hp)))
s.append(txt(gx(v), Y0p + Hp + 20, str(v)))
s.append(''
% (f(X0), f(Y0p), f(X0), f(Y0p + Hp), f(X0 + Wp), f(Y0p + Hp)))
s.append(poly([gx(t) for t in t2], [gy(v) for v in imp_slow(t2)], "mode slow"))
s.append(poly([gx(t) for t in t_fast], [gy(v) for v in imp_fast(t_fast)], "mode fast"))
s.append(poly([gx(t) for t in t2], [gy(v) for v in imp], "trace sA"))
s.append(txt(X0 + Wp / 2, Y0p + Hp + 44,
"years after a 1 W yr m\u207b\u00b2 pulse of heat"))
s.append(txt(gx(1450), gy(float(imp_slow(np.array(1450.0)))) - 11,
"slow mode, \u03c4 = %.0f yr" % base.tau_s, "middle", 12, "lab2"))
s.append(txt(gx(150), gy(float(imp_fast(np.array(30.0)))) + 4,
"fast mode, \u03c4 = %.1f yr" % base.tau_f, "start", 12, "lab2"))
s.append(txt(gx(340), gy(float(imp_slow(np.array(340.0)))) - 9,
"total response", "start", 12, "lab"))
open(os.path.join(FIGDIR, "fig2.svg"), "w", encoding="utf-8").write("\n".join(s))
# ---- FIG 3: gamma and lambda sweeps of realized fraction ----
s = []
PW, PH = 215.0, 230.0
for panel, (vals, kw, lab) in enumerate([
([0.0, 0.35, 0.70, 1.20, 2.00], "gam", "γ"),
([0.82, 1.13, 1.44, 1.80], "lam", "λ")]):
px0 = 68.0 + panel * (PW + 74.0)
py0 = 28.0
def hx(t, px0=px0):
return px0 + PW * (math.log10(t) - math.log10(1.0)) / (math.log10(2000.0))
def hy(v, py0=py0):
return py0 + PH - PH * v
for v in [0.2, 0.4, 0.6, 0.8, 1.0]:
s.append(''
% (f(px0), f(hy(v)), f(px0 + PW), f(hy(v))))
if panel == 0:
s.append(txt(px0 - 10, hy(v) + 4, "%d%%" % int(100 * v), "end"))
for v in [1, 10, 100, 1000]:
s.append(''
% (f(hx(v)), f(py0), f(hx(v)), f(py0 + PH)))
s.append(txt(hx(v), py0 + PH + 20, str(v)))
s.append(''
% (f(px0), f(py0), f(px0), f(py0 + PH), f(px0 + PW), f(py0 + PH)))
tt = np.logspace(0, math.log10(2000.0), 120)
for i, v in enumerate(vals):
m = TwoBox(**{kw: v})
fr = m.step_T(tt) / m.ecs
cls = "trace s%s" % "ABCDE"[i]
if kw == "lam" and abs(v - LAM_BASE) < 1e-9:
cls += " hi"
if kw == "gam" and abs(v - GAM_BASE) < 1e-9:
cls += " hi"
s.append(poly([hx(t) for t in tt], [hy(x) for x in fr], cls))
s.append(txt(hx(1600), hy(float(m.step_T(1600.0)) / m.ecs) + 4,
"%s=%.2f" % (lab, v), "start", 11, "lab2"))
s.append(txt(px0 + PW / 2, py0 + PH + 42,
"years after the step, %s varied" % lab, "middle", 12))
open(os.path.join(FIGDIR, "fig3.svg"), "w", encoding="utf-8").write("\n".join(s))
# ---- FIG 4: TCR/ECS against gamma, a family in lambda ----
X0, Y0p, Wp, Hp = 86.0, 30.0, 470.0, 250.0
gg = np.linspace(0.0, 2.0, 120)
def jx(g):
return X0 + Wp * g / 2.0
def jy(v):
return Y0p + Hp - Hp * (v - 0.3) / 0.75
s = []
for v in np.arange(0.4, 1.05, 0.1):
s.append(''
% (f(X0), f(jy(v)), f(X0 + Wp), f(jy(v))))
s.append(txt(X0 - 10, jy(v) + 4, "%.1f" % v, "end"))
for v in np.arange(0.0, 2.01, 0.25):
s.append(''
% (f(jx(v)), f(Y0p), f(jx(v)), f(Y0p + Hp)))
s.append(txt(jx(v), Y0p + Hp + 20, "%.2f" % v))
s.append(''
% (f(jx(PUB["gam_lo"])), f(Y0p), f(jx(PUB["gam_hi"]) - jx(PUB["gam_lo"])), f(Hp)))
s.append(''
% (f(X0), f(Y0p), f(X0), f(Y0p + Hp), f(X0 + Wp), f(Y0p + Hp)))
for i, L in enumerate([0.82, 1.13, 1.44, 1.80]):
ys = []
for g in gg:
m = TwoBox(gam=g, lam=L)
ys.append(m.tcr() / m.ecs)
cls = "trace s%s" % "ABCD"[i]
if abs(L - LAM_BASE) < 1e-9:
cls += " hi"
s.append(poly([jx(g) for g in gg], [jy(v) for v in ys], cls))
s.append(txt(jx(2.0) + 6, jy(ys[-1]) + 4, "λ=%.2f" % L, "start", 11, "lab2"))
s.append('' % (f(jx(GAM_BASE)), f(jy(tcr / ecs))))
s.append(txt(jx(GAM_BASE), jy(tcr / ecs) - 12, "baseline %.3f" % (tcr / ecs), "middle", 11, "lab"))
s.append(txt(X0 + Wp / 2, Y0p + Hp + 44, "ocean heat uptake efficiency γ (W m⁻² K⁻¹)"))
s.append(txt(jx(0.85), Y0p + 16, "published CMIP5 range", "middle", 11, "lab2"))
open(os.path.join(FIGDIR, "fig4.svg"), "w", encoding="utf-8").write("\n".join(s))
# ---- FIG 5: Monte Carlo convergence ----
X0, Y0p, Wp, Hp = 92.0, 30.0, 470.0, 250.0
lo_n, hi_n = 1.0, math.log10(N_KEEP)
ylo, yhi = final_mean - 0.30, final_mean + 0.30
def kx(n):
return X0 + Wp * (math.log10(n) - lo_n) / (hi_n - lo_n)
def ky(v):
y = Y0p + Hp - Hp * (v - ylo) / (yhi - ylo)
return min(max(y, Y0p), Y0p + Hp)
s = []
for v in np.arange(round(ylo, 1), yhi, 0.1):
s.append(''
% (f(X0), f(ky(v)), f(X0 + Wp), f(ky(v))))
s.append(txt(X0 - 10, ky(v) + 4, "%.2f" % v, "end"))
for d in [1, 2, 3, 4, 5]:
for k in range(1, 10):
v = (10 ** d) * k
if 10 <= v <= N_KEEP:
s.append(''
% (" major" if k == 1 else "", f(kx(v)), f(Y0p), f(kx(v)), f(Y0p + Hp)))
up = [ky(m + e) for m, e in zip(conv_m, conv_se)]
dn = [ky(m - e) for m, e in zip(conv_m, conv_se)]
pts = (" ".join("%s,%s" % (f(kx(n)), f(y)) for n, y in zip(conv_n, up))
+ " " + " ".join("%s,%s" % (f(kx(n)), f(y)) for n, y in zip(conv_n[::-1], dn[::-1])))
s.append('' % pts)
s.append(''
% (f(X0), f(ky(final_mean)), f(X0 + Wp), f(ky(final_mean))))
s.append(''
% (f(X0), f(Y0p), f(X0), f(Y0p + Hp), f(X0 + Wp), f(Y0p + Hp)))
s.append(poly([kx(n) for n in conv_n], [ky(m) for m in conv_m], "trace sA"))
for v in [10, 100, 1000, 10000, 100000]:
s.append(txt(kx(v), Y0p + Hp + 20, "1e%d" % int(math.log10(v))))
s.append(txt(X0 + Wp / 2, Y0p + Hp + 44, "Monte Carlo trials"))
s.append(txt(kx(N_KEEP), ky(final_mean) - 9,
"%.4f ± %.4f K" % (final_mean, r_com[2]), "end", 12, "lab"))
open(os.path.join(FIGDIR, "fig5.svg"), "w", encoding="utf-8").write("\n".join(s))
# ---- data dump for the LaTeX pgfplots ----
with open(os.path.join(FIGDIR, "figdata.txt"), "w", encoding="utf-8") as fh:
fh.write("# step response\n")
for t in np.logspace(math.log10(0.5), math.log10(3000.0), 60):
fh.write("%.6f %.6f %.6f\n" % (t, float(base.step_T(t)), float(base.step_T0(t))))
fh.write("# gamma sweep: gam tau_f tau_s tcr ratio f100\n")
for r in gam_rows:
fh.write("%.4f %.4f %.4f %.4f %.4f %.4f\n"
% (r[0], r[1], r[2] if np.isfinite(r[2]) else -1, r[4], r[5], r[6]))
fh.write("# lambda sweep: lam ecs tau_f tau_s tcr ratio f100\n")
for r in lam_rows:
fh.write("%.4f %.4f %.4f %.4f %.4f %.4f %.4f\n" % r)
fh.write("# convergence: n mean se\n")
for n, m, e in zip(conv_n, conv_m, conv_se):
fh.write("%d %.6f %.6f\n" % (n, m, e))
fh.write("# tcr/ecs vs gamma, for lam in 0.82 1.13 1.44 1.80\n")
for g in np.linspace(0.0, 2.0, 41):
row = [g]
for L in [0.82, 1.13, 1.44, 1.80]:
m = TwoBox(gam=g, lam=L)
row.append(m.tcr() / m.ecs)
fh.write(" ".join("%.5f" % x for x in row) + "\n")
sys.stderr.write("figures written to %s\n" % FIGDIR)