# Public copy: internal compute identifiers and local environment paths were removed or generalized. Scientific values are retained. Paths may need adapting and some original run inputs are not included. """ Toy climate model: CO2 -> radiative forcing -> warming. SUPERSEDED — kept only as a reference implementation. The delivered result for this project comes from real CESM2/CAM4 aquaplanet runs on the `the compute environment` cluster (see co2_sweep_diagnostics.csv and the two figures), not from this module. Nothing in the deliverables is derived from the code below. Two independent, very cheap models: 1. `two_layer` 0-D global energy balance with a mixed layer coupled to a deep ocean (Held et al. 2010; Geoffroy et al. 2013). Gives transient warming and the TCR/ECS split. 2. `ebm_1d` 1-D diffusive energy balance model on latitude (Budyko 1969; Sellers 1969; North 1975) with an ice-albedo parameterisation. CALIBRATION CAVEAT for `ebm_1d`: as parameterised here (A_OLR=203.3, D_DIFF=0.55, T_ICE=-10 C) the planet is too warm for ice to form anywhere — a diagnostic run gives ice_fraction = 0.000 and ice_edge_lat = NaN at every CO2 level from 280 to 1120 ppm, so the ice-albedo feedback NEVER ACTIVATES. With the albedo held uniform the diffusive response is very nearly uniform in latitude too: doubling CO2 warms the tropics and the poles by the same 1.86 K. So this module as shipped does NOT demonstrate ice-albedo feedback, polar amplification, or ice-edge retreat. Producing those would require retuning (a larger A_OLR and/or smaller D_DIFF to cool the high latitudes below T_ICE); that retuning was never completed and the module was dropped in favour of the CESM runs. Both run in well under a second on one CPU core. """ import numpy as np # ---------------------------------------------------------------------------- # CO2 radiative forcing # ---------------------------------------------------------------------------- C0_PPM = 280.0 # pre-industrial CO2 F_2X = 3.71 # W m-2 forcing per doubling of CO2 (IPCC AR5/AR6) def co2_forcing(ppm, c0=C0_PPM): """Radiative forcing of CO2 [W m-2]. Myhre et al. (1998) log law.""" return (F_2X / np.log(2.0)) * np.log(np.asarray(ppm, float) / c0) # ---------------------------------------------------------------------------- # 1. Zero-dimensional two-layer ocean model # ---------------------------------------------------------------------------- # C dT/dt = F(t) - lam*T - gam*(T - Td) mixed layer / atmosphere # Cd dTd/dt = gam*(T - Td) deep ocean # Heat capacities in W yr m-2 K-1 (C = 8 ~ 60 m slab, Cd = 100 ~ 800 m). LAM = 1.13 # net climate feedback [W m-2 K-1] -> ECS = F_2X/LAM = 3.3 K GAM = 0.70 # ocean heat uptake efficiency [W m-2 K-1] C_UP, C_DEEP = 8.0, 100.0 def two_layer(ppm_series, dt=0.5, lam=LAM, gam=GAM, c_up=C_UP, c_deep=C_DEEP): """Integrate the two-layer model over a CO2 pathway (one value per year). Returns dict with year, ppm, forcing, T (surface anomaly, K), T_deep (K) and ocean heat uptake N (W m-2). """ ppm = np.asarray(ppm_series, float) nyr = ppm.size t_fine = np.arange(0.0, nyr, dt) f_fine = np.interp(t_fine, np.arange(nyr), co2_forcing(ppm)) T = np.zeros(t_fine.size) Td = np.zeros(t_fine.size) for i in range(t_fine.size - 1): # forward Euler dT = (f_fine[i] - lam * T[i] - gam * (T[i] - Td[i])) / c_up dTd = gam * (T[i] - Td[i]) / c_deep T[i + 1] = T[i] + dt * dT Td[i + 1] = Td[i] + dt * dTd keep = np.isclose(t_fine % 1.0, 0.0) return dict(year=t_fine[keep], ppm=ppm[: keep.sum()], forcing=f_fine[keep], T=T[keep], T_deep=Td[keep], N=gam * (T[keep] - Td[keep])) def ecs(lam=LAM): """Equilibrium climate sensitivity [K per CO2 doubling].""" return F_2X / lam # ---------------------------------------------------------------------------- # 2. One-dimensional diffusive EBM with ice-albedo feedback # ---------------------------------------------------------------------------- # C dT/dt = Q*S(x)*(1-alpha) - (A + B*T) + d/dx[ D (1-x^2) dT/dx ], x = sin(lat) Q_SOLAR = 341.5 # global-mean insolation S0/4 [W m-2] S2 = -0.482 # 2nd Legendre coefficient of the annual-mean insolation A_OLR = 203.3 # outgoing longwave intercept [W m-2] (at 280 ppm) B_OLR = 2.00 # outgoing longwave slope [W m-2 K-1] D_DIFF = 0.55 # meridional heat diffusivity [W m-2 K-1] ALB_OCEAN, ALB_ICE, T_ICE = 0.32, 0.62, -10.0 # ice forms below -10 C def ebm_1d(ppm=C0_PPM, n=90, years=400, dt=0.02, d_diff=D_DIFF): """Equilibrium zonal-mean surface temperature [deg C] vs latitude. Implicit-in-diffusion, explicit-in-albedo time stepping to steady state. """ x = np.linspace(-1 + 1 / n, 1 - 1 / n, n) # cell centres, sin(lat) dx = x[1] - x[0] xe = np.linspace(-1, 1, n + 1) # cell edges w = (1 - xe ** 2) / dx ** 2 # edge diffusion weights # tridiagonal Laplacian-like operator L (no flux through the poles) L = np.zeros((n, n)) for i in range(n): if i > 0: L[i, i - 1] = w[i] if i < n - 1: L[i, i + 1] = w[i + 1] L[i, i] = -(w[i] + w[i + 1]) S = 1.0 + S2 * 0.5 * (3 * x ** 2 - 1) # insolation shape A_eff = A_OLR - co2_forcing(ppm) # more CO2 -> less OLR C_slab = 5.0 # W yr m-2 K-1 M = np.eye(n) * (C_slab / dt + B_OLR) - d_diff * L # implicit operator Minv = np.linalg.inv(M) T = 15.0 - 30.0 * x ** 2 # first guess [deg C] for _ in range(int(years / dt)): alb = np.where(T <= T_ICE, ALB_ICE, ALB_OCEAN) rhs = C_slab / dt * T + Q_SOLAR * S * (1 - alb) - A_eff T = Minv @ rhs lat = np.degrees(np.arcsin(x)) area = np.full(n, 1.0 / n) # equal-area in sin(lat) ice = T <= T_ICE return dict(lat=lat, x=x, T=T, area=area, ice=ice, T_global=float((T * area).sum()), ice_fraction=float(area[ice].sum()), ice_edge_lat=(float(np.min(np.abs(lat[ice]))) if ice.any() else np.nan))