""" core-leakage.py -- the Science Journaling Club's own simplified model for "The Core Is Leaking: Ruthenium-100, Hawaiian Lava, and Gold From the Center of the Earth" (Field Notes, Paper Analysis, Geochemistry). THIS IS NOT THE PAPER'S ANALYSIS. Messling et al. (2025, Nature 642, 376-380, doi:10.1038/s41586-025-09003-0) did the measuring, the chemistry and the real error budget. What follows is a back-of-the-envelope reconstruction that a high-school club can actually check line by line. Where our numbers land near theirs, that is a sanity check, not a replication. ============================================================================= PART A -- TWO-COMPONENT ISOTOPE MIXING ============================================================================= Ruthenium isotope compositions are quoted in epsilon notation: the deviation of a sample's 100Ru/101Ru from a terrestrial standard, in parts per 10,000. So eps100Ru = +0.09 means the sample is 0.09/10000 = 9 parts per MILLION rich in 100Ru. That is the whole signal. Nine parts per million. If a mantle source is a mechanical mixture of ordinary mantle and core metal, the isotope ratio of the mixture is NOT the simple average of the two end-members -- it is weighted by how much ruthenium each end-member brings to the table. The core is ~800x richer in Ru than the mantle, so a vanishingly small mass of core dominates the isotope budget. For a mass fraction f of core material: eps_mix = [ f*Cc*eps_c + (1-f)*Cm*eps_m ] / [ f*Cc + (1-f)*Cm ] Invert for f: f = Cm*(eps_mix - eps_m) / [ Cc*(eps_c - eps_mix) + Cm*(eps_mix - eps_m) ] ASSUMPTIONS (A1-A7), stated plainly: A1 Two components only. Real mantle sources are messier: recycled crust, late-veneer heterogeneity, and melt-stage fractionation all exist. A2 eps100Ru of the ambient convecting mantle is 0.00 +/- 0.02 (2s), by definition of the modern-mantle reference. A3 eps100Ru of the core is +0.25, allowed to range uniformly over +0.15 to +0.35. This is NOT measured -- nobody has a core sample. It is inferred from meteorite Mo-Zr-Ru correlations and the non-carbonaceous chondrite array. It is the single softest number in the model. A4 Ru concentrations: core 4000 ppb, primitive mantle 5.0 ppb (McDonough 2003; Palme & O'Neill 2014). Both given +/- 25% / +/- 20%. A5 The measured Hawaiian excess is eps100Ru = +0.09 +/- 0.03 (2 s.e. on the Hawaiian sample population). External 2s reproducibility on a single measurement is much worse, +/- 0.13 on reference material OREAS 684; the paper's resolution comes from repeating the measurement, not from any single heroic run. We model BOTH. A6 The erupted basalt inherits its source's core mass fraction unchanged. This is false in detail (partial melting concentrates siderophiles into the melt) but it keeps the gold estimate conservative -- a LOWER bound. A7 No isotope fractionation during melting, transport or eruption. Mass- independent 100Ru excesses are nucleosynthetic; melting cannot make them. Uncertainty is propagated by Monte Carlo (200,000 seeded draws) rather than by calculus, because f is a ratio of sums and its distribution is skewed. ============================================================================= PART B -- THE INDEPENDENT TUNGSTEN CHECK ============================================================================= The same f, pushed through the same mixing equation with tungsten's numbers, predicts a 182W deficit. Tungsten is quoted in mu notation: parts per MILLION deviation in 182W/184W. The core is only ~36x richer in W than the mantle, so W is a far blunter instrument -- but it is an INDEPENDENT one, and the two must agree or the story collapses. B1 W: core 470 ppb, primitive mantle 13 ppb (McDonough 2003). B2 mu182W of the core = -220 +/- 30 ppm relative to the silicate Earth (Hf-W systematics; Kleine & Walker 2017). Also soft. B3 Observed Hawaiian OIB mu182W runs about -5 to -12 ppm; the global OIB range reaches about -18 to -20 ppm (Mundl et al. 2017). ============================================================================= PART C -- THE GOLD LEDGER ============================================================================= How much core-derived gold reaches the crust per year at Hawaii? Au_flux = Q * rho * f * C_Au,core C1 Q = present-day Hawaiian magmatic volume flux = 0.15 km3/yr, range 0.08-0.21 km3/yr. The low end is Kilauea's measured magma supply alone (0.079 +/- 0.004 km3/yr, Dvorak & Dzurisin 1993); the high end includes the whole active shield complex plus intrusive material. C2 rho = 2900 kg/m3 (2800-3000), dense basalt. C3 C_Au,core = 500 ppb (300-800). Core gold is an inference, not a measurement. C4 The gold rides along passively with the entrained core metal. No separate gold-enrichment mechanism is invoked. C5 A plain wedding band is taken as 4.0 g of gold. C6 "All the gold ever mined" is taken as 2.2e8 kg (~220,000 t, World Gold Council stock estimate). Used only for scale. C7 A melt-enrichment VARIANT is also reported: if gold is perfectly incompatible and the melt fraction is F = 0.08, the erupted basalt scavenges gold from 1/F = 12.5x its own mass of source mantle. This is the generous end. The headline number does not use it. Nothing here is a prediction about ore deposits. Hawaiian basalt is not a gold mine and never will be; the concentrations are parts per billion. Seeded RNG: results are bit-for-bit reproducible. """ import numpy as np RNG = np.random.default_rng(2025_09_003) NDRAW = 200_000 SEP = "=" * 74 sub = lambda s: print("\n" + s + "\n" + "-" * len(s)) def pct(a, q): return float(np.percentile(a, q)) def band(a, label, unit="", sci=False): fmt = (lambda v: f"{v:.4g}") if not sci else (lambda v: f"{v:.4e}") print(f" {label:<34s} {fmt(pct(a,50))} {unit}" f" [68% {fmt(pct(a,16))} - {fmt(pct(a,84))}]" f" [95% {fmt(pct(a,2.5))} - {fmt(pct(a,97.5))}]") # ===================================================================== print(SEP) print("CORE LEAKAGE -- CLUB MIXING + MASS-FLUX MODEL") print("Anchor: Messling et al. (2025) Nature 642, 376-380") print(" doi:10.1038/s41586-025-09003-0") print("This is the club's own simplified model, NOT the paper's analysis.") print(SEP) # ------------------------------------------------------- fixed inputs EPS_MEAS, EPS_MEAS_2SE = 0.09, 0.03 # Hawaiian population mean, 2 s.e. EPS_EXT_2SD = 0.13 # single-run external reproducibility EPS_MANTLE, EPS_MANTLE_2SD = 0.00, 0.02 EPS_CORE_LO, EPS_CORE_HI = 0.15, 0.35 RU_CORE, RU_CORE_REL = 4000.0, 0.25 # ppb RU_MANT, RU_MANT_REL = 5.0, 0.20 # ppb W_CORE, W_MANT = 470.0, 13.0 # ppb MU_CORE, MU_CORE_SD = -220.0, 30.0 # ppm Q_LO, Q_MID, Q_HI = 0.08, 0.15, 0.21 # km3/yr RHO_LO, RHO_HI = 2800.0, 3000.0 # kg/m3 AU_CORE, AU_LO, AU_HI = 500.0, 300.0, 800.0 # ppb RING_G = 4.0 # grams of gold in a plain band MINED_KG = 2.2e8 # all gold ever mined, kg CHAIN_MYR = 85.0 # Hawaiian-Emperor chain age, Myr MELT_F = 0.08 # melt fraction, variant only sub("PART 0 -- THE SIGNAL, IN PLAIN UNITS") print(f" Hawaiian eps100Ru excess = +{EPS_MEAS:.2f} +/- {EPS_MEAS_2SE:.2f} (2 s.e.)") print(f" ... same thing in parts per million = +{EPS_MEAS*100:.0f} +/- {EPS_MEAS_2SE*100:.0f} ppm") print(f" Single-measurement 2s reproducibility (OREAS 684) = {EPS_EXT_2SD:.2f} eps" f" = {EPS_EXT_2SD*100:.0f} ppm") print(f" -> one measurement CANNOT see this. The population mean can.") print(f" Signal-to-noise on a single run = {EPS_MEAS/EPS_EXT_2SD:.2f}") print(f" Ru core/mantle concentration ratio = {RU_CORE/RU_MANT:.0f}x (Ru leverage)") print(f" W core/mantle concentration ratio = {W_CORE/W_MANT:.0f}x (W leverage)") # ===================================================================== sub("PART A -- MIXING FRACTION f (analytic central case)") def solve_f(eps_mix, eps_core, eps_mant, c_core, c_mant): num = c_mant * (eps_mix - eps_mant) den = c_core * (eps_core - eps_mix) + num return num / den f_central = solve_f(EPS_MEAS, 0.25, EPS_MANTLE, RU_CORE, RU_MANT) print(f" eps_mix={EPS_MEAS}, eps_core=0.25, eps_mantle=0.00," f" Cc={RU_CORE:.0f} ppb, Cm={RU_MANT:.1f} ppb") print(f" f_central = {f_central:.3e} = {f_central*100:.4f} % core by mass") print(f" = 1 part in {1/f_central:,.0f}") sub("PART A -- MONTE CARLO PROPAGATION (200,000 seeded draws)") eps_mix_d = RNG.normal(EPS_MEAS, EPS_MEAS_2SE / 2.0, NDRAW) eps_core_d = RNG.uniform(EPS_CORE_LO, EPS_CORE_HI, NDRAW) eps_mant_d = RNG.normal(EPS_MANTLE, EPS_MANTLE_2SD / 2.0, NDRAW) ru_core_d = RNG.normal(RU_CORE, RU_CORE * RU_CORE_REL / 2.0, NDRAW) ru_mant_d = RNG.normal(RU_MANT, RU_MANT * RU_MANT_REL / 2.0, NDRAW) f_d = solve_f(eps_mix_d, eps_core_d, eps_mant_d, ru_core_d, ru_mant_d) ok = (f_d > 0) & (f_d < 1) & (ru_core_d > 0) & (ru_mant_d > 0) & (eps_core_d > eps_mix_d) rejected = NDRAW - int(ok.sum()) f_ok = f_d[ok] print(f" draws kept = {ok.sum():,} of {NDRAW:,}" f" (rejected {rejected:,} = {100*rejected/NDRAW:.2f}% as unphysical)") band(f_ok, "f (mass fraction)", "", sci=True) band(f_ok * 100, "f (percent core by mass)", "%") print(f" Interpretation: the Hawaiian source is ~{pct(f_ok,50)*100:.3f}% core metal,") print(f" i.e. about 1 part in {1/pct(f_ok,50):,.0f}. The published paper puts an") print(f" upper bound of <0.25% (bulk core) or <0.3% (oxide layer); our median") print(f" sits comfortably inside that, which is the sanity check we wanted.") sub("PART A -- SENSITIVITY: what if the excess were smaller?") print(f" {'eps100Ru':>10s} {'ppm':>6s} {'f (bulk core)':>16s} {'1 part in':>14s}") for e in (0.05, 0.07, 0.09, 0.11, 0.17): fv = solve_f(e, 0.25, 0.0, RU_CORE, RU_MANT) print(f" {e:>10.2f} {e*100:>6.0f} {fv:>16.3e} {1/fv:>14,.0f}") print(" (0.05 = a deliberately conservative reading of the anomaly;") print(" 0.11 = Kilauea Iki lava lake; 0.17 = Napali Member, Kauai)") sub("PART A -- SENSITIVITY: the softest number, eps100Ru of the core") print(f" {'eps_core':>10s} {'f (bulk core)':>16s} {'percent':>10s}") for ec in (0.15, 0.20, 0.25, 0.30, 0.35): fv = solve_f(EPS_MEAS, ec, 0.0, RU_CORE, RU_MANT) print(f" {ec:>10.2f} {fv:>16.3e} {fv*100:>9.4f}%") print(" A factor of ~2.3 in the assumed core composition moves f by ~2.7x.") print(" Nobody has ever measured a piece of the core. This is the weak joint.") # ===================================================================== sub("PART B -- INDEPENDENT TUNGSTEN CROSS-CHECK") mu_core_d = RNG.normal(MU_CORE, MU_CORE_SD, len(f_ok)) mu_pred = (f_ok * W_CORE * mu_core_d) / (f_ok * W_CORE + (1 - f_ok) * W_MANT) band(mu_pred, "predicted mu182W", "ppm") print(f" Observed Hawaiian OIB mu182W ~ -5 to -12 ppm (Mundl+ 2017)") print(f" Global OIB range ~ 0 to -18 ppm") inside = float(np.mean((mu_pred > -12) & (mu_pred < -5)) * 100) print(f" Fraction of our predictions landing inside the Hawaiian window: {inside:.1f}%") print(" The Ru clock and the W clock were set independently and they agree to") print(" within a factor of ~2. That agreement is the paper's real argument --") print(" not either isotope system on its own.") sub("PART B -- INVERTED: what f would W alone demand?") print(f" {'mu182W obs':>12s} {'f from W alone':>17s} {'ratio to Ru f':>15s}") f_ru = pct(f_ok, 50) for mu in (-3.0, -5.0, -8.0, -12.0, -18.0): fw = (W_MANT * mu) / (W_CORE * (MU_CORE - mu) + W_MANT * mu) print(f" {mu:>12.1f} {fw:>17.3e} {fw/f_ru:>15.2f}") # ===================================================================== sub("PART C -- THE GOLD LEDGER (headline, conservative)") q_d = RNG.triangular(Q_LO, Q_MID, Q_HI, len(f_ok)) * 1e9 # m3/yr rho_d = RNG.uniform(RHO_LO, RHO_HI, len(f_ok)) # kg/m3 au_d = RNG.triangular(AU_LO, AU_CORE, AU_HI, len(f_ok)) * 1e-9 # kg Au / kg mass_d = q_d * rho_d # kg basalt / yr core_d = mass_d * f_ok # kg core metal / yr gold_d = core_d * au_d # kg core-derived Au / yr print(f" Hawaiian magmatic volume flux Q = {Q_MID} km3/yr (range {Q_LO}-{Q_HI})") band(mass_d, "erupted+intruded basalt mass", "kg/yr", sci=True) band(core_d, "core metal delivered", "kg/yr", sci=True) band(gold_d, "CORE-DERIVED GOLD", "kg/yr") g50 = pct(gold_d, 50) g16, g84 = pct(gold_d, 16), pct(gold_d, 84) print(f"\n Median gold flux = {g50:.1f} kg/yr" f" ({g50*1000/365.25/24/60:.2f} g/min)") ring_yr = (RING_G / 1000.0) / gold_d band(ring_yr * 365.25 * 24 * 60, "time per 4.0 g wedding band", "minutes") print(f" Median: one wedding band's worth of core gold every" f" {pct(ring_yr,50)*365.25*24*60:.1f} minutes") print(f" Equivalently {1/pct(ring_yr,50):,.0f} rings per year," f" 68% range {1/pct(ring_yr,84):,.0f}-{1/pct(ring_yr,16):,.0f}") # the "how big is that lump of iron" framing side = (pct(core_d, 50) / 7800.0) ** (1 / 3) # iron, kg/m3 print(f" The core metal itself: {pct(core_d,50):,.0f} kg/yr -- a cube of iron" f" {side:.2f} m on a side, every year.") sub("PART C -- CUMULATIVE OVER THE HAWAIIAN-EMPEROR CHAIN") cum = gold_d * CHAIN_MYR * 1e6 band(cum, f"gold over {CHAIN_MYR:.0f} Myr", "kg", sci=True) band(cum / 1000.0, f"gold over {CHAIN_MYR:.0f} Myr", "tonnes", sci=True) band(cum / MINED_KG, "as multiples of all gold ever mined", "x") print(f" (all gold ever mined taken as {MINED_KG:.1e} kg = ~{MINED_KG/1000:,.0f} t)") print(" CAVEAT: flux was NOT constant over 85 Myr -- Hawaiian output has grown") print(" several-fold along the chain. Treat this as an order of magnitude only.") sub("PART C -- MELT-ENRICHMENT VARIANT (generous end)") gold_var = gold_d / MELT_F band(gold_var, "gold if Au perfectly incompatible", "kg/yr") print(f" melt fraction F = {MELT_F}; the melt scavenges {1/MELT_F:.1f}x its own") print(" mass of source mantle. We do NOT quote this as the headline.") print(f" Variant median: {pct(gold_var,50):.0f} kg/yr" f" = one band every {(RING_G/1000)/pct(gold_var,50)*365.25*24*60:.2f} minutes") sub("PART C -- WHAT FRACTION OF HAWAIIAN BASALT GOLD IS CORE-DERIVED?") AU_MANTLE_PPB = 1.7 au_from_core_ppb = f_ok * AU_CORE * 1e0 # ppb contributed au_from_mant_ppb = (1 - f_ok) * AU_MANTLE_PPB frac = au_from_core_ppb / (au_from_core_ppb + au_from_mant_ppb) * 100 band(frac, "core share of source-rock gold", "%") print(f" Primitive-mantle Au taken as {AU_MANTLE_PPB} ppb (Palme & O'Neill 2014).") print(" So even in the leakiest rocks on the planet, most of the gold is still") print(" ordinary mantle gold. The core signal is a minority shareholder.") # ===================================================================== sub("FIGURE DATA -- mixing curve (Figure 2)") print(f" {'log10(f)':>10s} {'f':>12s} {'eps100Ru':>10s} {'mu182W':>10s}") for lg in np.arange(-5.0, -1.99, 0.25): fv = 10 ** lg e = (fv * RU_CORE * 0.25) / (fv * RU_CORE + (1 - fv) * RU_MANT) m = (fv * W_CORE * MU_CORE) / (fv * W_CORE + (1 - fv) * W_MANT) print(f" {lg:>10.2f} {fv:>12.3e} {e:>10.4f} {m:>10.2f}") sub("FIGURE DATA -- histogram of f (Figure 3), 18 bins in log10 f") h, edges = np.histogram(np.log10(f_ok), bins=18) for i in range(len(h)): print(f" bin {i:>2d} log10f {edges[i]:>7.3f} to {edges[i+1]:>7.3f}" f" n = {h[i]:>7d} frac = {h[i]/len(f_ok):.4f}") sub("FIGURE DATA -- measured samples (Figure 4 panel)") SAMPLES = [ ("Kilauea Iki lava lake, Hawai'i", 0.11, 0.04), ("Hawaiian OIB mean", 0.09, 0.03), ("Napali Member, Kaua'i", 0.17, 0.13), ("Eifel peridotite (ambient mantle)", 0.02, 0.03), ("Rhenish picrite (ambient mantle)", -0.01, 0.13), ] print(f" {'sample':<36s} {'eps100Ru':>9s} {'2s':>7s} {'ppm':>7s} {'f median':>12s}") for name, e, s in SAMPLES: fv = solve_f(e, 0.25, 0.0, RU_CORE, RU_MANT) if e > 0 else float("nan") fs = f"{fv:.2e}" if fv == fv and fv > 0 else "--" print(f" {name:<36s} {e:>9.2f} {s:>7.2f} {e*100:>7.0f} {fs:>12s}") sub("HEADLINE NUMBERS FOR THE ARTICLE") print(f" signal : +{EPS_MEAS*100:.0f} ppm excess 100Ru (eps100Ru = +0.09 +/- 0.03)") print(f" core mass fraction f : {pct(f_ok,50)*100:.3f}% " f"(68% {pct(f_ok,16)*100:.3f}-{pct(f_ok,84)*100:.3f}%; " f"95% {pct(f_ok,2.5)*100:.3f}-{pct(f_ok,97.5)*100:.3f}%)") print(f" = 1 part in : {1/pct(f_ok,50):,.0f}") print(f" predicted mu182W : {pct(mu_pred,50):.1f} ppm " f"(68% {pct(mu_pred,16):.1f} to {pct(mu_pred,84):.1f})") print(f" core metal to surface : {pct(core_d,50):,.0f} kg/yr") print(f" core-derived gold : {g50:.0f} kg/yr (68% {g16:.0f}-{g84:.0f})") print(f" wedding band interval : {pct(ring_yr,50)*365.25*24*60:.1f} minutes") print(f" over 85 Myr : {pct(cum,50)/1e6/1000:.2f} million tonnes" f" = {pct(cum/MINED_KG,50):.0f}x all gold ever mined") print(f" core share of Au in the source rock: {pct(frac,50):.1f}%") print(SEP)