"""Part 2: parameterised mantle thermal evolution and the end of geological habitability. One energy equation per planet, C dT/dt = H(t) - Q(T) [- melt heat loss, stagnant lid only], with T the mantle potential temperature, H the radiogenic heat production in the convecting mantle and Q the convective heat loss. Four formulations of Q(T) are run side by side: F1 "classical" mobile lid, Nu ~ Ra^beta with beta = 1/3 (Schubert, Turcotte & Olson 2001; Turcotte & Schubert 2014; used for exoplanets by Kite, Manga & Gaidos 2009). Calibrated on Earth's present potential temperature only; it under-predicts Earth's present heat loss (the well-known low-Urey-ratio problem). F2 "flat" mobile lid with beta = 0, i.e. heat flow almost independent of mantle temperature (Christensen 1985; Conrad & Hager 1999; Korenaga 2003, 2008). Calibrated on BOTH Earth's present potential temperature and present convective heat loss (38 TW), as in Korenaga (2008). F2h "hybrid" flat (beta = 0) while the mantle is hotter than today's Earth, classical (beta = 1/3) below, anchored to the same present-day heat loss. This is our own combination of the two published scalings (Korenaga's argument for small beta rests on melting-induced dehydration, which vanishes as melting wanes). F3 "stagnant" stagnant lid, Solomatov & Moresi (2000) scaling as implemented by Foley & Smye (2018, Astrobiology 18, 873), with their melting and degassing model and parameter values (simplified: quasi-steady lid, no crustal sequestration of heat sources, no depletion of mantle carbon; each of these simplifications LENGTHENS the degassing lifetime). End of geological habitability. The THRESHOLD is published; the mobile-lid degassing LAWS are ours: * degassing criterion: CO2 degassing per unit area below f_crit times the modern Earth's (f_crit = 0.1: Tajika 2003, 2007 as quoted by Kadoya & Tajika 2014 for the Earth's insolation; adopted for STAGNANT-lid planets by Foley & Smye 2018 and Unterborn et al. 2022). The critical rate depends on insolation (Kadoya & Tajika 2014, fig. 5) and planets below it cycle between frozen and thawed states rather than freezing for good. - stagnant lid (F3): the degassing model of Foley & Smye (2018); - mobile lid, "ridge" law (baseline, OUR construction): spreading rate x depth of the melting zone x fixed mantle carbon; ignores arc and metamorphic degassing of subducted carbon; - mobile lid, "plate-speed" law (variant, OUR construction): degassing proportional to the plate speed alone as long as any melting occurs (arc-dominated limit; the coupling used by Oosterloo et al. 2021). * melting criterion: decompression melting stops (f_crit -> 0), i.e. the adiabat no longer crosses the DRY peridotite solidus above the base of the lid (surface for a mobile lid). More permissive than the degassing criterion, but not the most permissive conceivable (hydrous melting and arc volcanism persist to lower temperatures). Thermostat (revision): the flat scaling (F2, F2h) has no feedback between mantle temperature and heat loss, so planets with a large early radiogenic budget would heat far above the peridotite liquidus (Oosterloo et al. 2021 note the same for beta = 0 and that melt transport must then cool the mantle). The mantle potential temperature of mobile-lid planets is therefore limited to T_cap (default 2200 K; zero-pressure dry liquidus 2053 K, Katz et al. 2003): above it the excess heat is assumed to be removed by magmatic (heat-pipe / magma-ocean) transport. T_cap = None switches the limit off (first-pass behaviour, clip at 4000 K). Units: time in Gyr, power in TW, temperature in K. """ import numpy as np from scipy.optimize import brentq from common import HALF_LIFE, AGE_EARTH, SEC_PER_GYR, earth_heat_initial_TW R_GAS = 8.314 T_REF = 1623.0 # K, Earth's present mantle potential temperature (1350 C; Herzberg et al. 2010) T_SURF = 273.0 # K C_EARTH = 7.0e27 / (1e12 * SEC_PER_GYR) # whole-Earth heat capacity 7e27 J/K (Korenaga 2008) in TW Gyr / K Q_CONV_EARTH = 38.0 # TW, present convective heat loss: 46 TW total (Jaupart et al. 2015) minus ~8 TW crustal radiogenic F_CC_EARTH = 0.35 # fraction of Earth's heat-producing elements stored in continental crust (~7 of 20 TW) DTSOL_DP = 120.0e-9 # K/Pa, slope of the dry solidus (Foley & Smye 2018 eq. 9, after Korenaga et al. 2002) GAMMA_AD = 1.0e-8 # K/Pa, mantle adiabat T_SOL0_HIGH = 1423.0 # K, zero-pressure dry solidus used by Foley & Smye (2018) T_SOL0_LOW = 1358.85 # K, zero-pressure dry solidus of Katz et al. (2003): 1085.7 C T_SOL0_MID = 0.5 * (T_SOL0_LOW + T_SOL0_HIGH) # 1390.9 K, baseline (mid-point of the two published values) T_LIQ0 = 2053.0 # K, zero-pressure dry peridotite liquidus of Katz et al. (2003): 1780 C T_CAP_DEFAULT = 2200.0 # K, thermostat for mobile-lid mantles (see module docstring); varied 2100-2300 K EARTH_DEGAS = 6.0e12 # mol/yr, modern Earth CO2 degassing (Marty & Tolstikhin 1998) A_EARTH_M2 = 4.0 * np.pi * 6378.1e3 ** 2 # mass-radius-gravity scalings for Earth-like composition (Valencia et al. 2006): R ~ M^0.27 EXP_R = 0.27 EXP_G = 1.0 - 2.0 * EXP_R # 0.46 EXP_A = 2.0 * EXP_R # 0.54 _H0 = earth_heat_initial_TW() _NUC = ("U238", "U235", "Th232", "K40") _H0_ARR = np.array([_H0[k] for k in _NUC]) _LAM = np.log(2.0) / np.array([HALF_LIFE[k] for k in _NUC]) def radiogenic_TW(t, fac, mass): """Heat production [TW] at age t [Gyr]. fac: (N, 4) initial concentration factors relative to the Earth for (238U, 235U, 232Th, 40K); mass: (N,) planet mass in Earth masses.""" return mass * (fac * (_H0_ARR * np.exp(-_LAM * t))).sum(axis=-1) # --------------------------------------------------------------------------------------------- # heat-loss laws # --------------------------------------------------------------------------------------------- def q_power_law(T, mass, beta, E, Q0, eta_fac=1.0): """Mobile-lid heat loss [TW]: Q = Q0 (A/A_E) (g/g_E)^beta (dT/dT_ref)^(1+beta) (eta_ref/eta(T))^beta, eta(T) = eta_fac * eta_ref * exp[(E/R)(1/T - 1/T_ref)]. For beta = 1/3 this is the classical boundary-layer scaling (independent of layer depth); for beta < 1/3 we keep the flux per unit area independent of layer depth as well (assumption).""" dT = np.maximum(T - T_SURF, 1.0) return (Q0 * mass ** EXP_A * (mass ** EXP_G / eta_fac) ** beta * (dT / (T_REF - T_SURF)) ** (1.0 + beta) * np.exp(-beta * E / R_GAS * (1.0 / T - 1.0 / T_REF))) class MobileLid: """Mobile-lid (plate-tectonic) formulations F1, F2 and F2h.""" def __init__(self, kind, E=300.0e3, f_cc=F_CC_EARTH, Q0=Q_CONV_EARTH, T_init=2000.0, beta=None, t_sol0=T_SOL0_MID, T_cap=T_CAP_DEFAULT): self.kind, self.E, self.f_cc, self.t_sol0 = kind, E, f_cc, t_sol0 self.T_cap = T_cap self.beta = {"classical": 1.0 / 3.0, "flat": 0.0, "hybrid": 0.0}[kind] if beta is None else beta self.Q0, self.T_init = Q0, T_init self.calibrated = False # heat loss [TW] def q(self, T, mass, eta_fac=1.0): if self.kind == "hybrid": return np.minimum(q_power_law(T, mass, 0.0, self.E, self.Q0, eta_fac), q_power_law(T, mass, 1.0 / 3.0, self.E, self.Q0, eta_fac)) return q_power_law(T, mass, self.beta, self.E, self.Q0, eta_fac) def prepare(self, p): """Cache the temperature-independent factors (pure speed-up; same numbers as q()).""" mass, eta = p["mass"], p["eta_fac"] p["_area"] = mass ** EXP_A p["_grav"] = mass ** EXP_G p["_pref_b"] = self.Q0 * p["_area"] * (p["_grav"] / eta) ** self.beta p["_pref_c"] = self.Q0 * p["_area"] * (p["_grav"] / eta) ** (1.0 / 3.0) p["_cap"] = C_EARTH * mass p["_h0"] = (1.0 - self.f_cc) * mass[:, None] * p["fac"] * _H0_ARR[None, :] return p def _q_fast(self, T, p): x = np.maximum(T - T_SURF, 1.0) / (T_REF - T_SURF) arr = self.E / R_GAS * (1.0 / T - 1.0 / T_REF) if self.kind == "hybrid": return np.minimum(p["_pref_b"] * x, p["_pref_c"] * x ** (4.0 / 3.0) * np.exp(-arr / 3.0)) if self.beta == 0.0: return p["_pref_b"] * x return p["_pref_b"] * x ** (1.0 + self.beta) * np.exp(-self.beta * arr) def dTdt(self, t, T, p): if "_h0" not in p: self.prepare(p) H = p["_h0"] @ np.exp(-_LAM * t) return (H - self._q_fast(T, p)) / p["_cap"] def rates(self, t, T, p): """dT/dt and both degassing laws from a single heat-loss evaluation (pure speed-up; the same numbers as dTdt(), degassing_rel() and degassing_arc_rel()).""" if "_h0" not in p: self.prepare(p) q = self._q_fast(T, p) k = (p["_h0"] @ np.exp(-_LAM * t) - q) / p["_cap"] flux_rel = q / (self._q_earth_now * p["_area"]) dT_rel = (T - T_SURF) / (T_REF - T_SURF) speed = flux_rel ** 2 / dT_rel ** 2 melt_rel = np.maximum(T - self.t_sol0, 0.0) / (T_REF - self.t_sol0) / p["_grav"] return k, speed * melt_rel, np.where(T > self.t_sol0, speed, 0.0) def degassing_rel(self, T, p): """CO2 degassing per unit area relative to the modern Earth's. Ridge degassing ~ areal spreading rate x depth of the melting zone x mantle carbon concentration (held at Earth's). Half-space cooling gives spreading rate S ~ Q^2 / (dT^2 A); the melting zone starts at P_i = (T - T_sol0) / (dTsol/dP - gamma), i.e. at depth P_i / (rho g).""" if "_h0" not in p: self.prepare(p) flux_rel = self._q_fast(T, p) / (self._q_earth_now * p["_area"]) # heat flux per area / Earth's dT_rel = (T - T_SURF) / (T_REF - T_SURF) melt_rel = np.maximum(T - self.t_sol0, 0.0) / (T_REF - self.t_sol0) / p["_grav"] return flux_rel ** 2 / dT_rel ** 2 * melt_rel def degassing_arc_rel(self, T, p): """Variant law: CO2 degassing per unit area proportional to the areal spreading (= subduction) rate alone, S ~ Q^2 / (dT^2 A), as long as the mantle still melts at all. This is the arc-dominated limit, in which the carbon carried down by the plates and returned at arcs controls the flux and the depth of the ridge melting zone does not matter.""" if "_h0" not in p: self.prepare(p) flux_rel = self._q_fast(T, p) / (self._q_earth_now * p["_area"]) dT_rel = (T - T_SURF) / (T_REF - T_SURF) return np.where(T > self.t_sol0, flux_rel ** 2 / dT_rel ** 2, 0.0) # ------------------------------------------------------------------ calibration on the Earth def calibrate(self): """classical: tune Q0 so that the Earth reaches T_REF at 4.5 Gyr from T_init (the result is insensitive to T_init). flat / hybrid: Q0 is the observed convective heat loss and T_init is tuned instead. Other beta values: tune T_init if a solution exists (heat-budget calibration), otherwise fall back to tuning Q0.""" self._q_earth_now = 1.0 cap, self.T_cap = self.T_cap, None # thermostat off while calibrating on the Earth def t_now(T_init, Q0): self.Q0 = Q0 one = {"fac": np.ones((1, 4)), "mass": np.ones(1), "eta_fac": np.ones(1)} # fresh: no stale cache T = integrate(self, np.array([T_init]), one, AGE_EARTH, dt=0.005)["T_end"] return float(T[0]) - T_REF if self.kind == "classical" or (self.kind == "power" and self._budget_fails()): self.Q0 = brentq(lambda q0: t_now(self.T_init, q0), 3.0, 120.0, xtol=1e-4) self.calibration = "T_p only" else: self.T_init = brentq(lambda ti: t_now(ti, self.Q0), 1500.0, 2600.0, xtol=1e-3) self.calibration = "T_p and heat loss" self._q_earth_now = float(self.q(np.array([T_REF]), np.ones(1))[0]) # the thermostat must not sit below the Earth-calibrated starting temperature self.T_cap = None if cap is None else max(cap, self.T_init) self.cap_raised = cap is not None and self.T_init > cap self.calibrated = True return self def _budget_fails(self): def fresh(): return {"fac": np.ones((1, 4)), "mass": np.ones(1), "eta_fac": np.ones(1)} self._q_earth_now = 1.0 cap, self.T_cap = self.T_cap, None hi = float(integrate(self, np.array([2600.0]), fresh(), AGE_EARTH, dt=0.005)["T_end"][0]) lo = float(integrate(self, np.array([1500.0]), fresh(), AGE_EARTH, dt=0.005)["T_end"][0]) self.T_cap = cap return not (lo < T_REF < hi) def earth_summary(self, f_crit=0.1): one = {"fac": np.ones((1, 4)), "mass": np.ones(1), "eta_fac": np.ones(1)} res = integrate(self, np.array([self.T_init]), one, 60.0, dt=0.005, f_crit=f_crit) H = (1.0 - self.f_cc) * float(radiogenic_TW(AGE_EARTH, one["fac"], one["mass"])[0]) Q = float(self.q(np.array([T_REF]), np.ones(1))[0]) return {"T_init_K": self.T_init, "Q_now_TW": Q, "H_mantle_now_TW": H, "urey_conv": H / Q, "cooling_K_per_Gyr": (Q - H) / C_EARTH, "lifetime_Gyr": float(res["lifetime"][0]), "lifetime_platespeed_law_Gyr": float(res["lifetime_arc"][0]), "melt_end_Gyr": float(res["melt_end"][0]), "calibration": self.calibration} class StagnantLid: """F3: Foley & Smye (2018) stagnant-lid model (their eqs. 1, 3, 4, 9-14, 24; parameter values of their table 1), simplified as described in the module docstring.""" rho, cp, k, kappa, alpha = 4000.0, 1250.0, 5.0, 1.0e-6, 3.0e-5 rho_m, L_m, rho_l = 2800.0, 6.0e5, 3300.0 c1, a_rh, c2 = 0.5, 2.5, 0.05 dphi_dP = 1.5e-10 gamma_erupt = 2.1e-8 # gamma_melt - gamma_mantle [K/Pa] mu_ref = 2.0e20 # Pa s at T_REF Rp_E, Rc_E, g_E = 6378.1e3, 3488.1e3, 9.8 C_tot_E = 1.0e22 # mol CO2 in the mantle of an Earth-mass planet def __init__(self, E=300.0e3, T_init=2000.0, t_sol0=T_SOL0_MID): self.kind, self.E, self.T_init, self.t_sol0 = "stagnant", E, T_init, t_sol0 self.f_cc = 0.0 self.calibration = "none (published parameters; viscosity 2e20 Pa s at 1623 K)" def state(self, T, p): if "_Rp" not in p: mass = p["mass"] p["_Rp"], p["_Rc"], p["_g"] = self.Rp_E * mass ** EXP_R, self.Rc_E * mass ** EXP_R, self.g_E * mass ** EXP_G p["_h0s"] = mass[:, None] * p["fac"] * _H0_ARR[None, :] * 1e12 Rp, Rc, g = p["_Rp"], p["_Rc"], p["_g"] d = Rp - Rc dT = np.maximum(T - T_SURF, 1.0) theta = self.E * dT / (R_GAS * T ** 2) mu = p["eta_fac"] * self.mu_ref * np.exp(self.E / R_GAS * (1.0 / T - 1.0 / T_REF)) Ra = self.rho * g * self.alpha * dT * d ** 3 / (self.kappa * mu) F = self.c1 * self.k * dT / d * theta ** (-4.0 / 3.0) * Ra ** (1.0 / 3.0) # W/m2 T_l = T - self.a_rh * R_GAS * T ** 2 / self.E delta = np.minimum(self.k * np.maximum(T_l - T_SURF, 1.0) / F, 0.95 * d) # quasi-steady lid P_i = np.maximum(T - self.t_sol0, 0.0) / (DTSOL_DP - GAMMA_AD) P_f = self.rho_l * g * delta melting = P_i > P_f phi = np.where(melting, 0.5 * (P_i - P_f) * self.dphi_dP, 0.0) d_m = P_i / (self.rho_l * g) v = self.c2 * self.kappa / d * (Ra / theta) ** (2.0 / 3.0) through = np.where(melting, 17.8 * np.pi * Rp * v * (d_m - delta), 0.0) # m3/s of mantle through the melt zone return {"Rp": Rp, "Rc": Rc, "delta": delta, "F": F, "P_i": P_i, "phi": phi, "through": through, "melting": melting} def dTdt(self, t, T, p): s = self.state(T, p) H = p["_h0s"] @ np.exp(-_LAM * t) A = 4.0 * np.pi * (s["Rp"] - s["delta"]) ** 2 V = 4.0 / 3.0 * np.pi * ((s["Rp"] - s["delta"]) ** 3 - s["Rc"] ** 3) dT_m = np.maximum(T - s["P_i"] * self.gamma_erupt - T_SURF, 0.0) melt_loss = s["through"] * s["phi"] * self.rho_m * (self.cp * dT_m + self.L_m) return (H - A * s["F"] - melt_loss) / (V * self.rho * self.cp) * SEC_PER_GYR def degassing_rel(self, T, p): s = self.state(T, p) v_mantle = 4.0 / 3.0 * np.pi * (s["Rp"] ** 3 - s["Rc"] ** 3) conc = self.C_tot_E * p["mass"] / v_mantle # mol/m3, Earth-like carbon per unit mass flux = s["through"] * conc * SEC_PER_GYR / 1e9 # mol/yr return (flux / (4.0 * np.pi * s["Rp"] ** 2)) / (EARTH_DEGAS / A_EARTH_M2) def calibrate(self): self.calibrated = True return self def earth_summary(self, f_crit=0.1): one = {"fac": np.ones((1, 4)), "mass": np.ones(1), "eta_fac": np.ones(1)} res = integrate(self, np.array([self.T_init]), one, 20.0, dt=0.005, f_crit=f_crit) T45 = float(integrate(self, np.array([self.T_init]), one, AGE_EARTH, dt=0.005)["T_end"][0]) s = self.state(np.array([T45]), one) Q = float(4 * np.pi * (s["Rp"][0] - s["delta"][0]) ** 2 * s["F"][0] / 1e12) H = float(radiogenic_TW(AGE_EARTH, one["fac"], one["mass"])[0]) return {"T_init_K": self.T_init, "Q_now_TW": Q, "H_mantle_now_TW": H, "urey_conv": H / Q, "cooling_K_per_Gyr": float(-self.dTdt(AGE_EARTH, np.array([T45]), one)[0]), "lifetime_Gyr": float(res["lifetime"][0]), "lifetime_platespeed_law_Gyr": float(res["lifetime"][0]), "melt_end_Gyr": float(res["melt_end"][0]), "T_at_4.5Gyr_K": T45, "calibration": self.calibration} # --------------------------------------------------------------------------------------------- # integrator (Heun / RK2, fixed step, vectorised over planets) # --------------------------------------------------------------------------------------------- def time_grid(t_max, dt, dt_early=None, t_early=1.0, dt_late=None, t_late=20.0): """Step schedule: fine steps `dt_early` up to `t_early` (the hot initial transient is stiff, especially for the stagnant lid with melt heat loss), then `dt`, then `dt_late` after `t_late`.""" edges = [0.0] parts = [] if dt_early is not None and dt_early < dt and t_max > t_early: parts.append((0.0, t_early, dt_early)) edges.append(t_early) start = edges[-1] if dt_late is not None and dt_late > dt and t_max > t_late: parts.append((start, t_late, dt)) parts.append((t_late, t_max, dt_late)) else: parts.append((start, t_max, dt)) grid = [np.array([0.0])] for a, b, h in parts: grid.append(np.linspace(a, b, int(round((b - a) / h)) + 1)[1:]) return np.concatenate(grid) def integrate(model, T0, p, t_max, dt=0.01, f_crit=0.1, record=None, dt_early=None, t_early=1.0, dt_late=None): """Integrate dT/dt from t = 0 to t_max (Heun / RK2, vectorised over planets). Returns the final temperature, the habitable lifetime (last time at which degassing per unit area >= f_crit x modern Earth; f_crit may be an array, one threshold per planet), the same for the plate-speed degassing law where the model has one (`lifetime_arc`), the time melting ends (last time degassing > 0), the peak temperature and whether the thermostat was reached, and optionally T(t) and degassing(t) at the grid times nearest to those in `record`.""" cap = getattr(model, "T_cap", None) T_hi = 4000.0 if cap is None else float(cap) T = np.minimum(np.array(T0, dtype=float), T_hi) grid = time_grid(t_max, dt, dt_early, t_early, dt_late) life = np.zeros_like(T) life_arc = np.zeros_like(T) melt_end = np.zeros_like(T) T_peak = T.copy() rec_T, rec_D, rec_t = [], [], [] rec_idx = set(int(np.argmin(np.abs(grid - x))) for x in record) if record is not None else set() want_life = getattr(model, "calibrated", False) has_arc = want_life and hasattr(model, "degassing_arc_rel") fast = want_life and hasattr(model, "rates") for i, t in enumerate(grid): k1 = None if fast: k1, D, D_arc = model.rates(t, T, p) life = np.where(D >= f_crit, t, life) melt_end = np.where(D > 0.0, t, melt_end) life_arc = np.where(D_arc >= f_crit, t, life_arc) if i in rec_idx: rec_T.append(T.copy()) rec_D.append(D.copy()) rec_t.append(t) elif want_life or i in rec_idx: D = model.degassing_rel(T, p) if want_life: life = np.where(D >= f_crit, t, life) melt_end = np.where(D > 0.0, t, melt_end) if has_arc: life_arc = np.where(model.degassing_arc_rel(T, p) >= f_crit, t, life_arc) if i in rec_idx: rec_T.append(T.copy()) rec_D.append(D.copy()) rec_t.append(t) if i == len(grid) - 1: break h = grid[i + 1] - t if k1 is None: k1 = model.dTdt(t, T, p) k2 = model.dTdt(t + h, np.clip(T + h * k1, 400.0, T_hi), p) T = np.clip(T + 0.5 * h * (k1 + k2), 400.0, T_hi) T_peak = np.maximum(T_peak, T) out = {"T_end": T, "lifetime": life, "melt_end": melt_end, "lifetime_arc": life_arc if has_arc else life, "T_peak": T_peak, "at_cap": T_peak >= T_hi - 1e-6} if record is not None: out["T_rec"] = np.array(rec_T) out["D_rec"] = np.array(rec_D) out["t_rec"] = np.array(rec_t) return out def make_model(name, E=300.0e3, f_cc=F_CC_EARTH, Q0=Q_CONV_EARTH, t_sol0=T_SOL0_MID, beta=None, T_cap=T_CAP_DEFAULT): if name == "F1_classical": return MobileLid("classical", E=E, f_cc=f_cc, t_sol0=t_sol0, T_cap=T_cap).calibrate() if name == "F2_flat": return MobileLid("flat", E=E, f_cc=f_cc, Q0=Q0, t_sol0=t_sol0, T_cap=T_cap).calibrate() if name == "F2h_hybrid": return MobileLid("hybrid", E=E, f_cc=f_cc, Q0=Q0, t_sol0=t_sol0, T_cap=T_cap).calibrate() if name == "F3_stagnant": return StagnantLid(E=E, t_sol0=t_sol0).calibrate() if name == "power": # exponent scan for the Earth twin only: no thermostat return MobileLid("power", E=E, f_cc=f_cc, Q0=Q0, t_sol0=t_sol0, beta=beta, T_cap=None).calibrate() raise ValueError(name) FORMULATIONS = ("F1_classical", "F2h_hybrid", "F2_flat", "F3_stagnant") LABELS = {"F1_classical": "F1 classical (beta = 1/3), mobile lid", "F2h_hybrid": "F2h hybrid (flat above, classical below today's T), mobile lid", "F2_flat": "F2 flat (beta = 0), mobile lid", "F3_stagnant": "F3 stagnant lid (Foley & Smye 2018)"} if __name__ == "__main__": for name in FORMULATIONS: m = make_model(name) print(name, {k: (round(v, 3) if isinstance(v, float) else v) for k, v in m.earth_summary().items()}) for beta in (0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 1 / 3): m = make_model("power", beta=beta) print(f"beta={beta:.2f}", {k: (round(v, 3) if isinstance(v, float) else v) for k, v in m.earth_summary().items()})