"""Signal models: magnitudes of an Anthropocene-like episode and how long each marker lasts. Everything here is arithmetic on published magnitudes (sources in SOURCES); nothing is fitted. Outputs: signals_carbon.csv, signals_radionuclides.csv, signals_pu244.json, signals_moon.csv, signals_summary.json """ import os, json import numpy as np, pandas as pd from common import DATA, save_json, A_EARTH, A_LAND from s04_cenogrid import amp, mass_from_amp, template, M0, DSRC, TAU_REC # peak of the recorded (bioturbated) response relative to the instantaneous shift: 0.82 for 100-kyr recovery, 7-kyr mixing PEAK_FRAC = float(template(np.linspace(0, 100, 2001)).max()) SOURCES = { "carbon_to_date": "Friedlingstein et al. 2025 (Global Carbon Budget 2024, ESSD 17, 965): cumulative 1850-2023 fossil CO2 emissions 490 +- 25 GtC, land-use change 225 +- 65 GtC, total 710 +- 70 GtC. Only the fossil part (about 500 GtC by 2025) adds 13C-depleted carbon to the exchangeable ocean-atmosphere-biosphere pool; land-use carbon moves within that pool (referee point 9)", "pu_fallout": "UNSCEAR 2000 Annex C: ~11 PBq of 239+240Pu dispersed by atmospheric weapons tests; 240/239 atom ratio 0.18", "pu244_ratio": "Steier et al. 2013 (NIM B 294, 160): global-fallout 244Pu/239Pu = (5.7 +- 1.0)e-5; a Baltic sediment core gives 7.9e-5", "pu244_ism": "Wallner et al. 2021 (Science 372, 742) as tabulated by Wang et al. 2021 (arXiv:2105.05178): interstellar 244Pu flux 980 +- 180 atoms/cm2/Myr over 0-9 Myr; crust layers 0-1.34, 1.34-4.57, 4.57-9.0 Myr with 34, 141, 6.3 counts; uptake 0.17", "pu244_2015": "Wallner et al. 2015 (Nat. Commun. 6, 5956): crust 0.5-25 Myr and sediment 0.5-2.2 Myr, at most one count per layer", "u236": "Sakaguchi et al. 2009: ~900 kg 236U in global fallout (a coral-based estimate gives 1060 kg)", "i129": "pre-nuclear marine 129I/127I = 1.5e-12 (Fehn and others)", "moon_diffusion": "Fassett & Thomson 2014 (JGR Planets 119, 2255): topographic diffusivity 5.5 m2/Myr on the maria over 3 Gyr", "moon_boulders": "Basilevsky et al. 2013, 2015: half of metre-sized boulders destroyed in tens of Myr, 90-99% in 200-300 Myr", "moon_gardening": "Speyerer et al. 2016 (Nature 538, 215): top centimetres of regolith churned >100 times faster than earlier models (the exact timescale was not verified here and is not used in any number); Apollo deep cores undisturbed below a few decimetres for ~0.5 Gyr", } LN2 = np.log(2.0) YEAR = 3.156e7 def carbon_table(): rows = [] for label, M in (("fossil carbon emitted to date", 500.0), ("fossil + land-use carbon to date (upper bound for the isotope effect)", 700.0), ("RCP8.5 to 2100-like", 2000.0), ("all conventional fossil fuel", 5000.0), ("PETM-size", 10000.0)): for D_yr in (300.0, 1e4, 1e5, 1e6): x = D_yr / (TAU_REC * 1e3) g = (1 - np.exp(-x)) / x # peak response of a first-order reservoir to a boxcar release a_pulse = float(amp(M)); a_eff = a_pulse * g rows.append({"case": label, "total_GtC": M, "release_duration_yr": D_yr, "mean_rate_GtC_per_yr": M / D_yr, "pulse_amplitude_permil": a_pulse, "duration_factor": g, "peak_amplitude_permil": a_eff, "recorded_peak_after_mixing_permil": a_eff * PEAK_FRAC, "equivalent_pulse_GtC": float(mass_from_amp(a_eff)), "years_of_2023_emissions": M / 11.0}) return pd.DataFrame(rows) def radionuclides(): """Time for which an Anthropocene-scale inventory stays above a measurement floor.""" pu239_atoms = 6.6e15 * 24110 * YEAR / LN2 # 6.6 PBq of 239Pu pu239_cm2 = pu239_atoms / (A_EARTH * 1e10) pu244_cm2_lo, pu244_cm2_hi = pu239_cm2 * 5.7e-5, pu239_cm2 * 7.94e-5 ocean_u_kg = 3.3e-9 * 1.37e21 # 3.3 ppb U in 1.37e21 kg sea water r236 = 1000.0 / ocean_u_kg * 238 / 236 # atom ratio after ocean mixing (1000 kg 236U) rows = [ {"nuclide": "240Pu", "half_life_yr": 6561, "anthropocene_level": 0.18 * pu239_cm2, "unit": "atoms/cm2 (global-mean fallout)", "floor": 1e4, "floor_basis": "AMS, ~1e6 atoms in a 100 cm2 sample"}, {"nuclide": "239Pu", "half_life_yr": 24110, "anthropocene_level": pu239_cm2, "unit": "atoms/cm2 (global-mean fallout)", "floor": 1e4, "floor_basis": "AMS, ~1e6 atoms in a 100 cm2 sample"}, {"nuclide": "129I", "half_life_yr": 15.7e6, "anthropocene_level": 4e-11, "unit": "129I/127I, ocean-mixed (assumes ~25x pre-nuclear)", "floor": 1e-13, "floor_basis": "AMS limit ~1e-14; in-situ fissiogenic 129I in old sediment ~1e-13"}, {"nuclide": "236U", "half_life_yr": 23.42e6, "anthropocene_level": r236, "unit": "236U/238U, ocean-mixed (1000 kg fallout)", "floor": 1e-12, "floor_basis": "in-situ nucleogenic 236U/238U in sediment 1e-14 to 1e-12"}, {"nuclide": "244Pu", "half_life_yr": 81.3e6, "anthropocene_level": 0.5 * (pu244_cm2_lo + pu244_cm2_hi), "unit": "atoms/cm2 (global-mean fallout)", "floor": 500.0, "floor_basis": "interstellar 244Pu accumulated in a 1-Myr layer (245-980 atoms/cm2/Myr)"}, ] for r in rows: r["factor_above_floor"] = r["anthropocene_level"] / r["floor"] r["window_yr"] = r["half_life_yr"] * np.log2(r["factor_above_floor"]) # the floor is uncertain by a factor ~10 either way: window +- 3.3 half-lives r["window_lo_yr"] = max(0.0, r["window_yr"] - 3.32 * r["half_life_yr"]) r["window_hi_yr"] = r["window_yr"] + 3.32 * r["half_life_yr"] out = pd.DataFrame(rows) meta = {"pu239_atoms_total": pu239_atoms, "pu239_mass_kg": pu239_atoms / 6.022e23 * 0.239, "pu239_atoms_per_cm2_global_mean": pu239_cm2, "pu244_atoms_per_cm2_global_mean": [pu244_cm2_lo, pu244_cm2_hi], "pu244_mass_g_global": [pu244_cm2_lo * A_EARTH * 1e10 / 6.022e23 * 244, pu244_cm2_hi * A_EARTH * 1e10 / 6.022e23 * 244], "u236_u238_ocean_mixed": r236} return out, meta def pu244_exclusion(meta): """What the published deep-sea 244Pu measurements exclude, and on what condition. Revised after review. (1) The 244Pu counted in the old layers is itself an unexplained excess, so a technogenic pulse is bounded by the WHOLE signal: its expected counts must not exceed counts + 2 sigma. The first pass used 3 sigma of the count error as if the counts were a known background, which gave a bound about three times too strong. (2) The bound assumes that a past fallout pulse ends up in an open-ocean ferromanganese crust with the same uptake (0.17) that Wallner et al. 2021 adopt for the steady interstellar dust influx. The one case where that can be tested, the 20th-century fallout in the top layer of the same crust, does not support it (see modern_layer_test), so the exclusion is CONDITIONAL on that uptake. """ spike_lo, spike_hi = meta["pu244_atoms_per_cm2_global_mean"] res = {"layers": []} # Tabulated 'interstellar' flux Phi (Wang et al. 2021, Table 2). Their Eq. 7 puts the deposition per cm2 # of Earth at Phi/4 (cross-section / surface). Both readings are carried; the weaker bound (Phi taken # as a terrestrial deposition) is the one quoted. layers = [("Crust-3/A 0-1.34 Myr", 0.0, 1.34, 930.0, 34.0, 17.0), ("Crust-3/B 1.34-4.57 Myr", 1.34, 4.57, 1990.0, 141.0, 19.0), ("Crust-3/C 4.57-9.0 Myr", 4.57, 9.0, 270.0, 6.3, 4.3)] for name, t0, t1, phi, counts, err in layers: T = t1 - t0 decay = 2 ** (-(0.5 * (t0 + t1)) / 81.3) r = {"layer": name, "counts": counts, "count_err": err} for conv, div in (("flux_in_space", 4.0), ("terrestrial", 1.0)): nat = phi / div * T # natural atoms/cm2 deposited in the layer for sp, lab in ((spike_lo, "lo"), (spike_hi, "hi")): ratio = sp * decay / nat r[f"bomb_spike_over_natural_{conv}_{lab}"] = ratio r[f"expected_counts_from_bomb_scale_pulse_{conv}_{lab}"] = counts * ratio # largest pulse compatible with the layer: its counts may not exceed counts + 2 sigma r[f"max_allowed_fraction_of_bomb_fallout_{conv}_{lab}"] = (counts + 2 * err) / (counts * ratio) res["layers"].append(r) def rng_of(layer, conv): v = [res["layers"][layer][f"max_allowed_fraction_of_bomb_fallout_{conv}_{lab}"] for lab in ("lo", "hi")] return [float(min(v)), float(max(v))] res["max_fraction_layerB_terrestrial"] = rng_of(1, "terrestrial") res["max_fraction_layerB_flux_in_space"] = rng_of(1, "flux_in_space") res["max_fraction_layerC_terrestrial"] = rng_of(2, "terrestrial") res["max_fraction_layerC_flux_in_space"] = rng_of(2, "flux_in_space") res["quoted_limit_fraction_of_bomb_fallout"] = float(max(res["max_fraction_layerB_terrestrial"] + res["max_fraction_layerC_terrestrial"])) # Modern-layer test: the same arithmetic applied to the top layer, which did receive bomb fallout. A = res["layers"][0] pred = sorted(A[k] for k in A if k.startswith("expected_counts_from_bomb_scale_pulse")) # The published interstellar residual is 34 +- 17 counts AFTER subtracting the anthropogenic part. If the # subtracted anthropogenic 244Pu were N counts, the Poisson error on the residual would be >= sqrt(34 + N), # so N <= 17^2 - 34. n_anthro_max = A["count_err"] ** 2 - A["counts"] res["modern_layer_test"] = { "predicted_bomb_counts_if_uptake_0p17": [float(pred[0]), float(pred[-1])], "published_residual_counts": [A["counts"], A["count_err"]], "max_anthropogenic_counts_compatible_with_published_error": float(n_anthro_max), "implied_modern_uptake_relative_to_model": [float(n_anthro_max / pred[-1]), float(n_anthro_max / pred[0])], "note": ("Sixty years after the tests, the crust has taken up at most this fraction of what the model assumes for a " "pulse (bound from the published error bar; the true anthropogenic count was not available to this study). " "Plausible reasons: bomb plutonium is still largely in the water column and is scavenged preferentially at " "ocean margins. Whether a pulse is eventually incorporated at 0.17 over 1e3-1e4 yr is untested.")} lim = res["quoted_limit_fraction_of_bomb_fallout"] lo_u, hi_u = res["modern_layer_test"]["implied_modern_uptake_relative_to_model"] res["limit_if_pulse_uptake_no_better_than_modern_layer"] = [float(lim / hi_u), float(lim / lo_u)] res["bomb_244Pu_mass_g"] = meta["pu244_mass_g_global"] res["conditional_exclusion_window_myr"] = [1.34, 9.0] res["tentative_exclusion_window_myr"] = 25.0 res["note"] = ("Layers B and C of Crust-3 (1.34-9.0 Myr; Wallner et al. 2021) hold 141 and 6.3 counts of 244Pu. A pulse equal to " "the 20th-century weapons-test fallout, IF incorporated at the interstellar uptake of 0.17, would have added the " "counts listed. The 0.5-25 Myr statement rests on a different crust (Wallner et al. 2015; 0, 1 and 1 counts in " "three layers) and is weak.") return res def moon(): kappa = 5.5 # m2/Myr, Fassett & Thomson 2014 (km-scale craters; smaller features degrade more slowly per unit kappa) d_eq = 150.0 # m, crater equilibrium diameter on ~3.5 Gyr maria (100-300 m) rows = [] for L in (1.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 1e4): t_eq = 3500.0 * (L / d_eq) ** 2 # scaled from crater equilibrium on the maria t_diff = 0.07 * L ** 2 / kappa # relief of scale L erased when kappa*t ~ 0.07 L^2 rows.append({"horizontal_scale_m": L, "erasure_time_myr_equilibrium_scaling": min(t_eq, 3500.0), "erasure_time_myr_diffusion": min(t_diff, 3500.0)}) return pd.DataFrame(rows) def coal_crosscheck(): """Order-of-magnitude side estimate (assumptions, not data): ancient land surface laid bare by coal mining.""" cum_coal_gt = 350.0 # ASSUMED cumulative world production (not verified against a primary source; plausible range 300-450 Gt) density = 1.35 # t/m3 seam_m = 3.0 # assumed mean mined thickness (1.5-10 m) area_km2 = cum_coal_gt * 1e9 / density / seam_m / 1e6 seam_kyr = 20.0 # time represented by a mined seam (1-2 m of coal ~ 10-20 m of peat at ~1 mm/yr) # space-time fraction of 3-500 Ma land surface laid bare, to set beside the PBDB probe at 1 km2 # (probability that a 1 km2 x 20 kyr cell holds a terrestrial fossil collection) st_frac = area_km2 * (seam_kyr / 1e3) / (A_LAND * 497.0) try: with open(os.path.join(DATA, "pbdb_summary.json")) as f: probe1 = json.load(f)["time_averages"]["Pterr_A1_tau0.01"]["3_500Myr"] except Exception: probe1 = float("nan") return {"assumed_cumulative_coal_Gt": cum_coal_gt, "assumed_mean_seam_m": seam_m, "bedding_area_exposed_km2": area_km2, "fraction_of_land_area": area_km2 / A_LAND, "assumed_time_per_seam_kyr": seam_kyr, "spacetime_fraction_of_3_500Ma_land_laid_bare": st_frac, "pbdb_probe_1km2_20kyr_3_500Ma": probe1, "ratio_coal_to_pbdb_probe_at_1km2": st_frac / probe1, "note": "assumption-level estimate used only to argue that the 'eyes factor' E is of order 10, not 1000"} def main(): c = carbon_table(); c.to_csv(os.path.join(DATA, "signals_carbon.csv"), index=False, float_format="%.4g") rn, meta = radionuclides(); rn.to_csv(os.path.join(DATA, "signals_radionuclides.csv"), index=False, float_format="%.4g") pu = pu244_exclusion(meta) m = moon(); m.to_csv(os.path.join(DATA, "signals_moon.csv"), index=False, float_format="%.4g") out = {"sources": SOURCES, "fallout_meta": meta, "pu244_exclusion": pu, "coal_crosscheck": coal_crosscheck(), "M0_GtC": M0, "d13c_source_permil": DSRC, "tau_rec_kyr": TAU_REC} save_json(out, "signals_summary.json") pd.set_option("display.width", 250) print(c.to_string()); print(rn.to_string()); print(json.dumps(pu, indent=1)); print(m.to_string()) print(json.dumps(out["coal_crosscheck"], indent=1)); print(json.dumps(meta, indent=1)) if __name__ == "__main__": main()