"""Synthesis: probability that a past visit / industrial episode would be known today. Combines the measured ingredients (s01-s06) with a small number of stated assumptions, and propagates the assumptions by Monte Carlo. Every assumed range is listed in ASSUMPTIONS and written to synthesis_summary.json. Revision after review (the numbered points are the referee's): * every installation probability is split into a SURFACE part (assumed survival and discovery of remains at or near today's land surface; no dataset behind it) and a DEEP-TIME part (PBDB probe x E x p_recognise; the only part tied to data). Both are written out, with the equivalent look-back time of the surface term and deep-only values by era (point 3) * exclusion counts are marginalised over the Monte Carlo draws (solve E[exp(-N p)] = 0.05) instead of 3 / median, and P(no known trace | one installation per stellar passage) is computed for the sibling fly-by study's passage counts (point 2) * human-scale footprint: variants with recognition applied per inspection and with a lowland-concentration factor (point 8) * carbon: the measured Cenozoic contribution and the assumed pre-Cenozoic contribution are reported separately; fossil carbon emitted to date is 500 GtC, not 700 (point 9) * fast markers: DM-SED coverage is a lower bound; the CENOGRID spacing coverage is added as a second floor (point 1) * A_PATCH and SCALE_UP are listed as assumptions and their effect is tabulated (point 9) * Moon comparison: like-for-like preservation ratio (point 7) Outputs: synthesis_funnel.csv, synthesis_funnel_sensitivity.csv, synthesis_installations.csv, synthesis_exclusion.csv, synthesis_human_scale_variants.csv, synthesis_carbon.csv, synthesis_fastmarkers.csv, synthesis_sampling_floors.csv, synthesis_headline.csv, synthesis_summary.json """ import os, json import numpy as np, pandas as pd from common import DATA, SEED, save_json, load_json, A_LAND from s04_cenogrid import amp, mass_from_amp, MASSES, TAU_REC NMC = 20000 ERAS = [(0.0, 3.0, "0-3 Ma"), (3.0, 66.0, "3-66 Ma"), (66.0, 252.0, "66-252 Ma"), (252.0, 539.0, "252-539 Ma"), (539.0, 2400.0, "539-2400 Ma")] WINDOWS = {"last_500Myr": 500.0, "last_2400Myr": 2400.0} # windows matching the sibling fly-by study's passage counts (systems passing within 1 pc of the Sun) FLYBY = {"last_2330Myr": (2330.0, 21000.0, 11000.0, 59000.0), "last_470Myr": (470.0, 4100.0, 2600.0, 8000.0)} AREAS = [1.0, 100.0, 1e4] TAUS = [3e-4, 1e-2, 1.0] H0 = 1e-2 # Myr, value of h used in s03 A_PATCH = 100.0 # km2, ASSUMED size of one independently preserved depositional patch M_NONMARINE = 0.5 # Sadler exponent for terrestrial sections (Macrostrat units: 0.49) M_MARINE = 0.38 # Sadler exponent for marine units (Macrostrat units: 0.38) SCALE_UP = 8.0 # ASSUMED: continental crust of the globe / Macrostrat North America column area ASSUMPTIONS = { "E_eyes": "log-uniform 3-300: close inspections by anyone who would report an artefact layer, per PBDB terrestrial collection", "p_recognise": "uniform 0.1-1: an inspected artefact-bearing bed is recognised and published as artificial; also absorbs how densely artefacts fill the footprint and whether they survive burial in recognisable form", "h_collection": "log-uniform 3-100 kyr: time window within which an inspected bed would hold the artefacts (time-averaging, reworking, surface residence of durable objects)", "t_surface": "log-uniform 0.03-1 Myr: e-folding time for remains to stay discoverable at or near the land surface (no dataset)", "p_surface": "found if still at the surface: U(0.1,0.7) for 1 km2, U(0.5,1) for 100 km2, U(0.8,1) for 1e4 km2, 1 for 1e6 km2 (no dataset)", "recognition_model": "baseline applies p_recognise once per event (a common-mode factor); a variant applies it per inspection. They differ only for the 1e6 km2 footprint, where several inspections are expected", "placement": "baseline places footprints uniformly on land; a variant multiplies the expected number of inspections by a lowland-concentration factor c_low, log-uniform 1-5 (settlements and preserved terrestrial sediment both favour depositional lowlands)", "F_res_Mesozoic": "log-uniform 0.14-0.8: fraction of Mesozoic time covered by d13C stratigraphy at <=50 kyr spacing (DM-SED gives 0.14 as a floor)", "F_res_Paleozoic": "log-uniform 0.05-0.5 (DM-SED: 0.20)", "F_res_Precambrian": "log-uniform 3e-4 to 1e-2 (DM-SED: 0.008 Neoproterozoic, ~0 older)", "noise_factor": "detection threshold relative to the Paleocene CENOGRID curve: U(1,2) Mesozoic, U(1.5,4) Paleozoic, U(2,6) Precambrian (assumed; no injection test was run on pre-Cenozoic data)", "p_attribution": "log-uniform 0.003-0.3: a catalogued carbon-cycle event is recognised as industrial (no data constrain this)", "ice_core": "a pulse of >= 100 GtC equivalent at any time in the last 0.8 Myr is taken as certainly registered by ice-core CO2 (assumed)", "slow_release": "releases lasting >= 100 kyr are judged against Myr-scale variability: detection threshold doubled (assumed)", "land_area": "148.9e6 km2 at all times", "A_PATCH": "100 km2: size of one independently preserved depositional patch. Drives 'P(some part preserved)' in the funnel (model output, not a measurement); sensitivity in synthesis_funnel_sensitivity.csv", "SCALE_UP": "8: continental crust of the globe / Macrostrat North America column area. Drives the 'regions worldwide holding the layer' counts for fast markers (model output)", "sadler_extrapolation": "completeness at 300 yr - 10 kyr is extrapolated from unit durations of Myr with exponents 0.5 (non-marine) and 0.38 (marine), 3-5 decades outside the fitted range", } def era_mean(x, y, lo, hi, w=None): s = (x >= lo) & (x < hi) if s.sum() == 0: return 0.0 if w is None: return float(np.mean(y[s])) return float(np.average(y[s], weights=w[s])) def logu(rng, a, b, n): return np.exp(rng.uniform(np.log(a), np.log(b), n)) def q595(x): return [float(v) for v in np.percentile(x, [5, 50, 95])] def n95_marginal(p, cl=0.05): """Number of independent events N such that E_p[exp(-N p)] = cl: the 95% upper limit on the number of past events given that none is known, with the systematic uncertainty in p (common to all events) marginalised over the Monte Carlo draws.""" lo, hi = 0.0, 1.0 while np.mean(np.exp(-hi * p)) > cl: hi *= 2.0 if hi > 1e15: return float("inf") for _ in range(200): mid = 0.5 * (lo + hi) if np.mean(np.exp(-mid * p)) > cl: lo = mid else: hi = mid return 0.5 * (lo + hi) def main(): rng = np.random.default_rng(SEED) cov = pd.read_csv(os.path.join(DATA, "macrostrat_coverage_1myr.csv")) bins = pd.read_csv(os.path.join(DATA, "macrostrat_bins.csv")) ms = load_json("macrostrat_summary.json") mp = load_json("mapsample_summary.json") mo = pd.read_csv(os.path.join(DATA, "mapsample_outcrop_by_age.csv")) sl = pd.read_csv(os.path.join(DATA, "pbdb_psampled_slices.csv")) rec = pd.read_csv(os.path.join(DATA, "cenogrid_recovery.csv")) dm = pd.read_csv(os.path.join(DATA, "dmsed_coverage.csv")) sig = load_json("signals_summary.json") cg = load_json("cenogrid_summary.json") out = {"assumptions": ASSUMPTIONS, "n_mc": NMC} age = cov["age_ma"].values c_glob = 1.0 / mp["ratio_usca_to_global"] # North America has more sedimentary cover than the globe out["na_to_global_correction"] = c_glob # unit duration by age (time-weighted over the geological bins) def unit_dur(lo, hi): w = np.minimum(bins["age_base_ma"], hi) - np.maximum(bins["age_top_ma"], lo) ok = (w > 0) & bins["na_unit_dur_med_myr"].notna() return float(np.average(bins.loc[ok, "na_unit_dur_med_myr"], weights=w[ok])) unw = ms["unit_statistics_weighted_vs_unweighted"] # ------------------------------------------------------------------ funnel (central values; MODEL OUTPUT) frows = []; srows = [] mgrid = mo["age_ma"].values; mw = np.gradient(mgrid) for lo, hi, lab in ERAS: f_nm = era_mean(age, cov["north_america_nonmarine_farea"].values, lo, hi) * c_glob f_sed = era_mean(age, cov["north_america_sed_farea"].values, lo, hi) * c_glob nm_share = f_nm / f_sed if f_sed > 0 else 0.0 Tu = unit_dur(lo, hi) f_slice = era_mean(mgrid, mo["global_sed_outcrop_frac_per_myr_slice"].values, max(lo, 0.0), hi, mw) ss = sl[(sl["age_ma"] >= lo) & (sl["age_ma"] < hi)] for tau in TAUS: win = tau + H0 C = min(1.0, (win / Tu) ** M_NONMARINE) for A in AREAS: F2 = f_nm * C p_any = 1.0 - np.exp(-F2 * max(1.0, A / A_PATCH)) key = f"A{int(A)}_tau{tau:g}" p_pb = float(np.average(ss[f"Pterr_{key}"], weights=ss["weight_myr"])) if len(ss) else 0.0 nov = float(np.average(ss[f"Nov_terr_tau{tau:g}"], weights=ss["weight_myr"])) if len(ss) else 0.0 frows.append({"era": lab, "area_km2": A, "occupation_myr": tau, "f_nonmarine_sediment_of_that_age_survives": f_nm, "unit_duration_myr": Tu, "completeness_at_window": C, "f_window_strata_survive": F2, "P_some_part_preserved_anywhere": p_any, "f_horizon_is_mapped_bedrock_surface": f_slice * nm_share * win, "expected_outcrop_of_horizon_km2": A * f_slice * nm_share * win, "n_terr_collections_worldwide_in_window": nov, "P_pbdb_collection_inside_footprint_and_window": p_pb}) # sensitivity of the modelled preservation probability (10,000 km2, 10 kyr) to its assumptions win = 1e-2 + H0 ukey = {"0-3 Ma": "0-66 Ma", "3-66 Ma": "0-66 Ma", "66-252 Ma": "66-252 Ma", "252-539 Ma": "252-539 Ma", "539-2400 Ma": "539-2400 Ma"}[lab] Tu_unw = unw[ukey].get("na_unit_dur_median_myr_unweighted", Tu) for patch in (10.0, 100.0, 1000.0): for m_sad in (0.38, 0.5, 0.6): for tu_lab, tu in (("time-weighted", Tu), ("unweighted", Tu_unw)): Cs = min(1.0, (win / tu) ** m_sad) srows.append({"era": lab, "a_patch_km2": patch, "sadler_exponent": m_sad, "unit_duration_kind": tu_lab, "unit_duration_myr": tu, "completeness": Cs, "P_some_part_preserved_10000km2_10kyr": 1.0 - np.exp(-f_nm * Cs * max(1.0, 1e4 / patch))}) funnel = pd.DataFrame(frows) funnel.to_csv(os.path.join(DATA, "synthesis_funnel.csv"), index=False, float_format="%.4g") fsens = pd.DataFrame(srows) fsens.to_csv(os.path.join(DATA, "synthesis_funnel_sensitivity.csv"), index=False, float_format="%.4g") out["preservation_sensitivity_10000km2_10kyr"] = { lab: {"baseline": float(funnel[(funnel["era"] == lab) & (funnel["area_km2"] == 1e4) & np.isclose(funnel["occupation_myr"], 1e-2)]["P_some_part_preserved_anywhere"].iloc[0]), "min": float(fsens[fsens["era"] == lab]["P_some_part_preserved_10000km2_10kyr"].min()), "max": float(fsens[fsens["era"] == lab]["P_some_part_preserved_10000km2_10kyr"].max())} for _, _, lab in ERAS} # ------------------------------------------------------------------ installations (Monte Carlo) # NOTE: the order of the random draws is kept exactly as in the first pass so that the combined # medians are unchanged; new draws (c_low) are taken after all first-pass draws. E = logu(rng, 3, 300, NMC); prec = rng.uniform(0.1, 1.0, NMC); h = logu(rng, 3e-3, 0.1, NMC) tsurf = logu(rng, 0.03, 1.0, NMC) psurf = {1.0: rng.uniform(0.1, 0.7, NMC), 100.0: rng.uniform(0.5, 1.0, NMC), 1e4: rng.uniform(0.8, 1.0, NMC), 1e6: np.ones(NMC)} t = sl["age_ma"].values; w = sl["weight_myr"].values allwin = dict(WINDOWS); allwin.update({k: v[0] for k, v in FLYBY.items()}) pre = (t >= 3.0) def deep_matrix(A, tau, per_inspection=False, c_low=None): """P(found and recognised through buried strata | age t), one row per draw.""" sc = (tau + h) / (tau + H0) if A <= 1e4: p0 = sl[f"Pterr_A{int(A)}_tau{tau:g}"].values hit = p0[None, :] * (sc * E)[:, None] if c_low is not None: hit = hit * c_low[:, None] if per_inspection: return 1.0 - np.exp(-hit * prec[:, None]) return np.minimum(1.0, hit) * prec[:, None] nov = sl[f"Nov_terr_tau{tau:g}"].values # dispersed footprint: Poisson over collections nhit = nov[None, :] * (sc * E * A / A_LAND)[:, None] if c_low is not None: nhit = nhit * c_low[:, None] if per_inspection: return 1.0 - np.exp(-nhit * prec[:, None]) return (1.0 - np.exp(-nhit)) * prec[:, None] def combine(A, deep): """Window averages of total / deep / surface parts, one value per draw.""" recent = psurf[A][:, None] * np.exp(-t[None, :] / tsurf[:, None]) res = {} for wn, T in allwin.items(): s = t < T d_int = (w[s][None, :] * deep[:, s] * (1 - recent[:, s])).sum(1) r_int = psurf[A] * tsurf * (1 - np.exp(-T / tsurf)) res[wn] = {"total": (d_int + r_int) / T, "deep": d_int / T, "surface": r_int / T} return res irows = []; xrows = [] inst_draws = {} for A in AREAS + [1e6]: for tau in TAUS: deep = deep_matrix(A, tau) res = combine(A, deep) inst_draws[(A, tau)] = res if A <= 1e4: p0 = sl[f"Pterr_A{int(A)}_tau{tau:g}"].values p0fm = sl[f"Pterr_fm_A{int(A)}_tau{tau:g}"].values else: p0 = p0fm = None nov = sl[f"Nov_terr_tau{tau:g}"].values # central (median-parameter) values by era Em, pm, hm, tm, psm = 30.0, 0.55, np.sqrt(3e-3 * 0.1), np.sqrt(0.03), float(np.median(psurf[A])) scm = (tau + hm) / (tau + H0) deepm = (np.minimum(1.0, p0 * scm * Em) if p0 is not None else (1.0 - np.exp(-nov * scm * Em * A / A_LAND))) * pm recm = psm * np.exp(-t / tm) tot = 1 - (1 - deepm) * (1 - recm) r = {"area_km2": A, "occupation_myr": tau} for lo, hi, lab in ERAS: r[f"P_known_central_{lab}"] = era_mean(t, tot, lo, hi, w) r[f"P_deep_only_central_{lab}"] = era_mean(t, deepm, lo, hi, w) for wn in WINDOWS: q = np.percentile(res[wn]["total"], [5, 16, 50, 84, 95]) r[f"P_known_{wn}_p05"], r[f"P_known_{wn}_p16"], r[f"P_known_{wn}_median"], r[f"P_known_{wn}_p84"], r[f"P_known_{wn}_p95"] = q r[f"P_known_{wn}_mean"] = float(res[wn]["total"].mean()) for part in ("deep", "surface"): qq = q595(res[wn][part]) r[f"P_{part}_part_{wn}_p05"], r[f"P_{part}_part_{wn}_median"], r[f"P_{part}_part_{wn}_p95"] = qq r[f"deep_time_share_{wn}_median"] = float(np.median(res[wn]["deep"] / res[wn]["total"])) r[f"P_deep_part_{wn}_max_over_draws"] = float(res[wn]["deep"].max()) # conditional values the Fermi model can use directly # (a) recent: P(known | it happened within the last 1 Myr) -- assumed surface term only rec1 = psurf[A] * tsurf * (1 - np.exp(-1.0 / tsurf)) / 1.0 r["P_known_if_within_last_1Myr_p05"], r["P_known_if_within_last_1Myr_median"], r["P_known_if_within_last_1Myr_p95"] = q595(rec1) teff = psurf[A] * tsurf r["surface_lookback_equivalent_myr_p05"], r["surface_lookback_equivalent_myr_median"], r["surface_lookback_equivalent_myr_p95"] = q595(teff) # (b) old: P(known | age between 3 Ma and T), buried strata only -- PBDB probe x E x p_recognise for Tn, T in (("3-500Ma", 500.0), ("3-2400Ma", 2400.0)): s = pre & (t < T) dold = (w[s][None, :] * deep[:, s]).sum(1) / (T - 3.0) r[f"P_deep_only_{Tn}_p05"], r[f"P_deep_only_{Tn}_median"], r[f"P_deep_only_{Tn}_p95"] = q595(dold) r[f"P_deep_only_{Tn}_max_over_draws"] = float(dold.max()) if p0 is not None: # same-formation collections counted once (lower bracket of the PBDB probe) r["ratio_formation_dedup_3_500Myr"] = float(np.sum(w[(t > 3) & (t < 500)] * p0fm[(t > 3) & (t < 500)]) / max(np.sum(w[(t > 3) & (t < 500)] * p0[(t > 3) & (t < 500)]), 1e-300)) # the measured quantity itself: share of 3-500 Ma land space-time inside a terrestrial collection r["pbdb_probe_3_500Ma"] = float(np.sum(w[(t >= 3) & (t < 500)] * p0[(t >= 3) & (t < 500)]) / np.sum(w[(t >= 3) & (t < 500)])) # extreme corner of the priors: E x p_recognise = 300, h = 100 kyr scx = (tau + 0.1) / (tau + H0) dx = np.minimum(1.0, p0 * scx * 300.0) r["P_deep_extreme_corner_last_500Myr"] = float(np.sum(w[t < 500] * dx[t < 500]) / 500.0) irows.append(r) # ---- exclusion counts: none known today -> how many could there have been? ------------- for wn, T in WINDOWS.items(): p = res[wn]["total"] xrows.append({"scenario": f"{A:g} km2, {tau * 1e6:g} yr", "area_km2": A, "occupation_myr": tau, "window": wn, "P_known_median": float(np.median(p)), "P_known_mean": float(p.mean()), "N95_marginalised": n95_marginal(p), "N95_naive_3_over_median": 3.0 / float(np.median(p)), "N95_at_p95": 3.0 / float(np.percentile(p, 95)), "N95_at_p05": 3.0 / float(np.percentile(p, 5)), "N95_marginalised_deep_part_only": n95_marginal(res[wn]["deep"])}) for wn, (T, n_mid, n_lo, n_hi) in FLYBY.items(): p = res[wn]["total"] rr = {"scenario": f"{A:g} km2, {tau * 1e6:g} yr", "area_km2": A, "occupation_myr": tau, "window": wn, "P_known_median": float(np.median(p)), "P_known_mean": float(p.mean()), "N95_marginalised": n95_marginal(p), "N95_naive_3_over_median": 3.0 / float(np.median(p))} for nm_, n_ in (("central", n_mid), ("low", n_lo), ("high", n_hi)): rr[f"n_passages_{nm_}"] = n_ rr[f"P_no_known_trace_one_per_passage_{nm_}"] = float(np.mean(np.exp(-n_ * p))) rr["expected_known_one_per_passage_at_median_p"] = float(n_mid * np.median(p)) rr["max_fraction_of_passages_95"] = float(min(1.0, n95_marginal(p) / n_mid)) xrows.append(rr) inst = pd.DataFrame(irows) inst.to_csv(os.path.join(DATA, "synthesis_installations.csv"), index=False, float_format="%.4g") excl = pd.DataFrame(xrows) excl.to_csv(os.path.join(DATA, "synthesis_exclusion.csv"), index=False, float_format="%.4g") # ------------------------------------------------------------------ global carbon episodes def curve(row): return np.array([row[f"P_flag_{M}GtC"] for M in MASSES], float) lm = np.log10(MASSES) rec3 = pd.read_csv(os.path.join(DATA, "cenogrid_recovery_k3.csv")) curves = {"5sigma_singular_outlier": {r["interval"]: curve(r) for _, r in rec.iterrows()}, "3sigma_among_largest_excursions": {r["interval"]: curve(r) for _, r in rec3.iterrows()}} def S(m_eff, which, ep): """Flagged fraction for an equivalent pulse of m_eff GtC, interpolated on the injection-recovery curve.""" m_eff = np.asarray(m_eff, float) return np.interp(np.log10(np.maximum(m_eff, 1.0)), lm, ep[which], left=0.0, right=ep[which][-1]) ceno_parts = [("Quaternary", 0.8, 2.58), ("Pliocene", 2.58, 5.33), ("Miocene", 5.33, 23.03), ("Oligocene", 23.03, 33.9), ("Eocene", 33.9, 56.0), ("Paleocene", 56.0, 66.0)] F_M = logu(rng, 0.14, 0.8, NMC); F_P = logu(rng, 0.05, 0.5, NMC); F_PC = logu(rng, 3e-4, 1e-2, NMC) n_M = rng.uniform(1, 2, NMC); n_P = rng.uniform(1.5, 4, NMC); n_PC = rng.uniform(2, 6, NMC) p_attr = logu(rng, 0.003, 0.3, NMC) cases = (("500 GtC (fossil carbon emitted to date)", 500.0), ("700 GtC (fossil + land-use carbon to date; upper bound)", 700.0), ("2000 GtC", 2000.0), ("5000 GtC (all conventional fossil fuel)", 5000.0)) crows = [] for crit, ep in curves.items(): for label, M in cases: for D_yr in (300.0, 1e4, 1e6): x = D_yr / (TAU_REC * 1e3); g = (1 - np.exp(-x)) / x m_eff = float(mass_from_amp(float(amp(M)) * g)) # slow (>= 100 kyr) releases are judged against Myr-scale variability: threshold doubled slow = 2.0 if D_yr >= 1e5 else 1.0 # 0-0.8 Ma: ice-core CO2 (assumed to register any pulse >= 100 GtC equivalent) p_ice = 1.0 if m_eff >= 100 else 0.0 p_cz_deep = sum((b_ - a_) * float(S(m_eff / slow, nm, ep)) for nm, a_, b_ in ceno_parts) / (66.0 - 0.8) p_cz = (0.8 * p_ice + (66.0 - 0.8) * p_cz_deep) / 66.0 p_mz = F_M * S(m_eff / (n_M * slow), "Paleocene", ep) p_pz = F_P * S(m_eff / (n_P * slow), "Paleocene", ep) p_pc = F_PC * S(m_eff / (n_PC * slow), "Paleocene", ep) pre500 = (186 * p_mz + 248 * p_pz) / 500.0 pre2400 = (186 * p_mz + 287 * p_pz + 1861 * p_pc) / 2400.0 ev500 = 66 * p_cz / 500.0 + pre500 ev2400 = 66 * p_cz / 2400.0 + pre2400 r = {"criterion": crit, "case": label, "total_GtC": M, "release_duration_yr": D_yr, "equivalent_pulse_GtC": m_eff, "peak_d13C_permil": float(amp(M)) * g, "recorded_peak_after_mixing_permil": float(amp(M)) * g * cg["template_peak_fraction"], "P_event_0-66Ma": p_cz, "P_event_0.8-66Ma_deep_sea_only": p_cz_deep, "P_event_66-252Ma_median": float(np.median(p_mz)), "P_event_252-539Ma_median": float(np.median(p_pz)), "P_event_539-2400Ma_median": float(np.median(p_pc)), "measured_cenozoic_contribution_500Myr": 66 * p_cz / 500.0, "measured_cenozoic_contribution_2400Myr": 66 * p_cz / 2400.0} for nm, arr in (("500Myr", pre500), ("2400Myr", pre2400)): q = q595(arr) r[f"assumed_preCenozoic_contribution_{nm}_p05"], r[f"assumed_preCenozoic_contribution_{nm}_median"], r[f"assumed_preCenozoic_contribution_{nm}_p95"] = q for nm, arr in (("last_500Myr", ev500), ("last_2400Myr", ev2400)): q = np.percentile(arr, [5, 50, 95]) r[f"P_event_{nm}_p05"], r[f"P_event_{nm}_median"], r[f"P_event_{nm}_p95"] = q q = np.percentile(arr * p_attr, [5, 50, 95]) r[f"P_known_{nm}_p05"], r[f"P_known_{nm}_median"], r[f"P_known_{nm}_p95"] = q r[f"N95_marginalised_catalogued_{nm}"] = n95_marginal(arr) r[f"N95_marginalised_recognised_{nm}"] = n95_marginal(arr * p_attr) crows.append(r) carb = pd.DataFrame(crows) carb.to_csv(os.path.join(DATA, "synthesis_carbon.csv"), index=False, float_format="%.4g") # ------------------------------------------------------------------ human-scale footprint: modelling variants # (drawn after every first-pass draw so that earlier results are unchanged) c_low = logu(rng, 1.0, 5.0, NMC) vrows = [] for A, tau in ((1e6, 3e-4), (1e6, 1e-2), (1e4, 1e-2)): for vlab, kw in (("baseline: recognition once per event, uniform placement", {}), ("recognition per inspection", {"per_inspection": True}), ("lowland concentration x1-5", {"c_low": c_low}), ("per inspection + lowland concentration", {"per_inspection": True, "c_low": c_low})): res = combine(A, deep_matrix(A, tau, **kw)) r = {"area_km2": A, "occupation_myr": tau, "variant": vlab} for wn in WINDOWS: p = res[wn]["total"] r[f"P_known_{wn}_p05"], r[f"P_known_{wn}_median"], r[f"P_known_{wn}_p95"] = q595(p) r[f"N95_marginalised_{wn}"] = n95_marginal(p) vrows.append(r) var = pd.DataFrame(vrows) var.to_csv(os.path.join(DATA, "synthesis_human_scale_variants.csv"), index=False, float_format="%.4g") # ------------------------------------------------------------------ fast markers confined to the episode layer mrows = [] na_cols = ms["n_columns"]["north_america"] for D_myr, sp, nm in ((3e-4, 0.1, "300 yr"), (1e-2, 3.33, "10 kyr"), (1.0, 333.0, "1 Myr")): for lo, hi, lab in ERAS[1:]: ncol = era_mean(age, cov["north_america_marine_fcols"].values, lo, hi) * na_cols Tu = unit_dur(lo, hi) C = min(1.0, (D_myr / Tu) ** M_MARINE) n_reg = ncol * C * SCALE_UP * c_glob # DM-SED coverage for this era (time-weighted over its bins) -- a LOWER BOUND def cover(setname): d = dm[(dm["set"] == setname) & (np.isclose(dm["max_spacing_kyr"], sp))] d = d[~d["bin"].isin(["2.58-500 Ma", "2.58-2400 Ma", "Mesozoic", "Paleozoic"])] ww = np.minimum(d["age_base_ma"], hi) - np.maximum(d["age_top_ma"], lo) ok = ww > 0 return float(np.average(d.loc[ok, "frac_time_ge1_site"], weights=ww[ok])) mrows.append({"episode_duration": nm, "era": lab, "marine_columns_NA": ncol, "unit_duration_myr": Tu, "completeness": C, "regions_worldwide_holding_the_layer_MODEL": n_reg, "P_preserved_somewhere_MODEL": 1 - np.exp(-n_reg), # same with the steepest measured Sadler slope (deep sea, 0.57) instead of 0.38 "regions_worldwide_MODEL_sadler0.57": ncol * min(1.0, (D_myr / Tu) ** 0.57) * SCALE_UP * c_glob, "P_preserved_somewhere_MODEL_sadler0.57": 1 - np.exp(-ncol * min(1.0, (D_myr / Tu) ** 0.57) * SCALE_UP * c_glob), "required_sample_spacing_kyr": sp, "dmsed_floor_frac_time_sampled_any_element": cover("all"), "dmsed_floor_frac_time_sampled_Pb": cover("pb"), "dmsed_floor_frac_time_sampled_carbon": cover("carbon")}) fast = pd.DataFrame(mrows) fast.to_csv(os.path.join(DATA, "synthesis_fastmarkers.csv"), index=False, float_format="%.4g") # floors on "share of 2.58-500 Ma time with at least one section sampled at spacing <= X": # DM-SED alone, CENOGRID alone, and the two combined (CENOGRID for 2.58-66 Ma, DM-SED for 66-500 Ma; # where DM-SED also covers Cenozoic time the larger of the two is taken, which can only undercount). sc_ = pd.DataFrame(cg["spacing_coverage"]) frow = [] for X in (0.1, 1.0, 3.33, 10.0, 50.0): d_all = dm[(dm["set"] == "all") & np.isclose(dm["max_spacing_kyr"], X)].set_index("bin")["frac_time_ge1_site"] dm_tot = float(d_all["2.58-500 Ma"]) * (500 - 2.58) dm_cz = float(d_all["Neogene"]) * (23.03 - 2.58) + float(d_all["Paleogene"]) * (66 - 23.03) cgx = sc_[np.isclose(sc_["max_spacing_kyr"], X)].iloc[0] cg_cz = float(cgx["covered_myr_2.58-66Ma"]) comb = (dm_tot - dm_cz) + max(dm_cz, cg_cz) frow.append({"max_spacing_kyr": X, "dmsed_only_share_2.58-500Ma": dm_tot / (500 - 2.58), "cenogrid_only_share_2.58-500Ma": cg_cz / (500 - 2.58), "combined_floor_share_2.58-500Ma": comb / (500 - 2.58), "cenogrid_share_of_2.58-66Ma": float(cgx["frac_time_2.58-66Ma"]), "cenogrid_share_of_0-2.58Ma": float(cgx["frac_time_0-2.58Ma"]), "dmsed_share_Mesozoic": float(d_all["Mesozoic"]), "dmsed_share_Paleozoic": float(d_all["Paleozoic"])}) floors = pd.DataFrame(frow) floors.to_csv(os.path.join(DATA, "synthesis_sampling_floors.csv"), index=False, float_format="%.4g") dmb = pd.read_csv(os.path.join(DATA, "dmsed_bins.csv")).set_index("bin") out["sampling_floors"] = { "table": frow, "expected_existing_samples_in_random_300yr_window": { "dmsed_2.58-500Ma": float(dmb.loc["2.58-500 Ma", "expected_samples_in_window_300yr"]), "cenogrid_0-66Ma": float(cg["samples_per_myr_0_66Ma"] * (3e-4 + 1e-3))}, "note": ("Lower bounds from two compilations. DM-SED holds no Cenozoic carbonate d13C; ocean-drilling isotope records " "other than the CENOGRID splice, core-scanning data, varved and ice-core records are not counted. Even where " "sampling is dense, bioturbation (about 7 kyr in the deep sea) smears a 300-yr layer.")} # ------------------------------------------------------------------ Moon versus Earth, like for like (preservation only) pres = {r["era"]: r["P_some_part_preserved_anywhere"] for r in out_records(funnel, 1e4, 1e-2)} out["moon_vs_earth_preservation"] = { "earth_P_some_part_of_10000km2_site_preserved_MODEL": pres, "lunar_retention_of_100m_plus_features_over_1Gyr": "about 1 (Fassett & Thomson 2014 scaling; s06)", "preservation_ratio_range": [float(1.0 / max(pres[k] for k in pres if k != "0-3 Ma")), float(1.0 / min(pres[k] for k in pres if k != "0-3 Ma"))], "note": "Preservation only. The first-pass '3-4 orders of magnitude' compared lunar preservation with terrestrial P(known), which also contains examination."} # ------------------------------------------------------------------ headline table hrows = [] for A, lab in ((1.0, "Surface installation 1 km2, 10 kyr"), (100.0, "Surface installation 100 km2, 10 kyr"), (1e4, "Surface installation 10,000 km2, 10 kyr"), (1e6, "Dispersed built footprint 1e6 km2 (human-scale), 300 yr")): tau = 3e-4 if A == 1e6 else 1e-2 r = inst[(inst["area_km2"] == A) & (np.isclose(inst["occupation_myr"], tau))].iloc[0] hrows.append({"scenario": lab, "P_known_500Myr": r["P_known_last_500Myr_median"], "lo_500": r["P_known_last_500Myr_p05"], "hi_500": r["P_known_last_500Myr_p95"], "P_known_2400Myr": r["P_known_last_2400Myr_median"], "lo_2400": r["P_known_last_2400Myr_p05"], "hi_2400": r["P_known_last_2400Myr_p95"], "deep_part_500Myr": r["P_deep_part_last_500Myr_median"], "deep_lo_500": r["P_deep_part_last_500Myr_p05"], "deep_hi_500": r["P_deep_part_last_500Myr_p95"], "deep_part_2400Myr": r["P_deep_part_last_2400Myr_median"], "deep_lo_2400": r["P_deep_part_last_2400Myr_p05"], "deep_hi_2400": r["P_deep_part_last_2400Myr_p95"], "surface_part_500Myr": r["P_surface_part_last_500Myr_median"], "surface_part_2400Myr": r["P_surface_part_last_2400Myr_median"], "P_known_if_within_last_1Myr": r["P_known_if_within_last_1Myr_median"], "surface_lookback_equivalent_myr": r["surface_lookback_equivalent_myr_median"], "P_event_500Myr": np.nan, "P_event_2400Myr": np.nan, "basis": "deep part = PBDB probe x E x p_recognise (data x assumptions); surface part = assumed"}) for M, D in ((500.0, 300.0), (700.0, 300.0), (2000.0, 300.0), (5000.0, 300.0), (5000.0, 1e4), (5000.0, 1e6)): r = carb[(carb["total_GtC"] == M) & (carb["release_duration_yr"] == D) & (carb["criterion"] == "5sigma_singular_outlier")].iloc[0] r3 = carb[(carb["total_GtC"] == M) & (carb["release_duration_yr"] == D) & (carb["criterion"] == "3sigma_among_largest_excursions")].iloc[0] hrows.append({"scenario": f"Global fossil-carbon episode {int(M)} GtC in {int(D)} yr", "P_known_500Myr": r["P_known_last_500Myr_median"], "P_event3sigma_500Myr": r3["P_event_last_500Myr_median"], "P_event3sigma_2400Myr": r3["P_event_last_2400Myr_median"], "lo_500": r["P_known_last_500Myr_p05"], "hi_500": r["P_known_last_500Myr_p95"], "P_known_2400Myr": r["P_known_last_2400Myr_median"], "lo_2400": r["P_known_last_2400Myr_p05"], "hi_2400": r["P_known_last_2400Myr_p95"], "P_event_500Myr": r["P_event_last_500Myr_median"], "P_event_2400Myr": r["P_event_last_2400Myr_median"], "P_event_measured_cenozoic_part_500Myr": r["measured_cenozoic_contribution_500Myr"], "P_event_measured_cenozoic_part_2400Myr": r["measured_cenozoic_contribution_2400Myr"], "P_event_assumed_preCenozoic_part_500Myr_p05": r["assumed_preCenozoic_contribution_500Myr_p05"], "P_event_assumed_preCenozoic_part_500Myr_p95": r["assumed_preCenozoic_contribution_500Myr_p95"], "basis": "Cenozoic part: CENOGRID injection-recovery (measured). Pre-Cenozoic part: assumed coverage and noise factors. P_known = P_event x assumed attribution"}) pu = sig["pu244_exclusion"] c0, c1 = pu["conditional_exclusion_window_myr"]; tent = pu["tentative_exclusion_window_myr"] hrows.append({"scenario": "Global 244Pu fallout >= 10% of 20th-century weapons tests (CONDITIONAL on crust uptake 0.17 for a pulse)", "P_known_500Myr": (c1 - c0) / 500.0, "lo_500": 0.0, "hi_500": tent / 500.0, "P_known_2400Myr": (c1 - c0) / 2400.0, "lo_2400": 0.0, "hi_2400": tent / 2400.0, "P_event_500Myr": np.nan, "P_event_2400Myr": np.nan, "basis": "published deep-sea 244Pu counts cover 1.34-9 Myr (Crust-3) and weakly 0.5-25 Myr (a second crust); zero if a fallout pulse is not taken up like interstellar dust"}) head = pd.DataFrame(hrows) head.to_csv(os.path.join(DATA, "synthesis_headline.csv"), index=False, float_format="%.3g") out["headline"] = hrows out["funnel_10000km2_10kyr"] = out_records(funnel, 1e4, 1e-2) out["unit_duration_myr_by_era"] = {lab: unit_dur(lo, hi) for lo, hi, lab in ERAS} out["exclusion"] = xrows out["human_scale_variants"] = vrows hv = var[(var["area_km2"] == 1e6) & np.isclose(var["occupation_myr"], 3e-4)] out["human_scale_envelope_500Myr"] = {"median_min": float(hv["P_known_last_500Myr_median"].min()), "median_max": float(hv["P_known_last_500Myr_median"].max()), "p05_min": float(hv["P_known_last_500Myr_p05"].min()), "p95_max": float(hv["P_known_last_500Myr_p95"].max())} out["human_scale_envelope_2400Myr"] = {"median_min": float(hv["P_known_last_2400Myr_median"].min()), "median_max": float(hv["P_known_last_2400Myr_median"].max()), "p05_min": float(hv["P_known_last_2400Myr_p05"].min()), "p95_max": float(hv["P_known_last_2400Myr_p95"].max())} save_json(out, "synthesis_summary.json") pd.set_option("display.width", 280); pd.set_option("display.max_columns", 60) print(funnel[(np.isclose(funnel["occupation_myr"], 1e-2))].to_string()) print(json.dumps(out["preservation_sensitivity_10000km2_10kyr"], indent=1)) print(inst.T.to_string()) print(excl.T.to_string()) print(var.T.to_string()) print(carb.T.to_string()) print(fast.to_string()) print(floors.T.to_string()) print(json.dumps(out["sampling_floors"]["expected_existing_samples_in_random_300yr_window"], indent=1)) print(json.dumps(out["moon_vs_earth_preservation"], indent=1)) print(head.T.to_string()) def out_records(funnel, A, tau): return funnel[(funnel["area_km2"] == A) & np.isclose(funnel["occupation_myr"], tau)].to_dict("records") if __name__ == "__main__": main()