#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ lightning-tree.py ================= Science Journaling Club — companion analysis for the Field Note "The Tree That Wants to Be Struck by Lightning". Anchor paper ------------ Gora, E. M., et al. (2025). How some tropical trees benefit from being struck by lightning: evidence for Dipteryx oleifera and other large-statured trees. New Phytologist 246, 1554-1566. doi:10.1111/nph.70062 WHAT THIS SCRIPT IS ------------------- This is the CLUB'S OWN simplified model. It is NOT the analysis published in the paper, and none of the numbers below were taken from the paper's own modelling. The paper reports field observations (93 tracked strikes, liana and neighbour mortality, a 14-fold lifetime-fecundity estimate). We reproduce none of that machinery. What we do here is ask, with high-school circuit physics and a crude Monte Carlo, whether those observations are even *physically reasonable*, and where the simple physics breaks. THE THREE CALCULATIONS ---------------------- PART A — Circuit physics, two limiting cases. We model a struck Dipteryx trunk and a liana draped on it as two resistive paths from crown to ground. Wood is treated as a homogeneous ohmic cylinder, R = rho * L / A, which is a serious simplification (see ASSUMPTIONS). Limit A1 "shared voltage": the two paths are a genuine parallel pair across the same two nodes (crown, ground). Current divides by conductance. Limit A2 "shared current": the flash attaches to the liana directly, or the liana bridges a gap with no thick alternative beside it, so the liana carries a current set upstream rather than by its own resistance. These two limits give wildly different answers. That gap is the point of the article: the observed 78%-of-lianas-dead / tree-untouched asymmetry can only live near limit A2. Limit A1 contains an exact cancellation (derived and verified numerically in PART B) in which stem thickness drops out entirely. PART C — Heating. Energy deposited = (action integral) x R, spread over the mass in the path. Temperature rise from E = m c dT, then latent heat for the tissue water. We report time-to-threshold rather than a final temperature, because the liana's "final temperature" in limit A2 is physically absurd (it disintegrates in the first microseconds). PART D — Fitness. A 300-year, one-year-timestep Monte Carlo. Strike hazard scales with height squared and crown area; per-strike effects (liana knock-down, canopy gap) follow the paper's reported magnitudes. We run two different contrasts, because they answer two different questions: D1 "tolerance alone" - identical trees, identical strike rate, one survives strikes and one does not. This is the contrast the paper's 14-fold figure refers to ("the ability to survive lightning increases lifetime fecundity 14-fold"). D2 "whole strategy" - a tall, wide-crowned tolerant tree against a shorter, narrower intolerant one that is struck less often. This is the evolutionary trade-off version, and it should give a SMALLER number, because the avoider buys something real with its shorter stature. ASSUMPTIONS (stated plainly; every one of them is arguable) ---------------------------------------------------------- 1. Wood is ohmic, homogeneous, isotropic and temperature-independent. False on all four counts. Real green wood is anisotropic, its resistivity falls steeply with temperature and rises catastrophically once sap boils, and much lightning current travels on a surface flashover channel rather than through tissue at all. 2. Trunk conducting cross-section = 60% of basal area (outer sapwood plus the moist cambial zone); heartwood is assumed non-participating. 3. Resistivities: Dipteryx trunk 120 ohm.m (deliberately at the LOW end, because the paper's hypothesis is that this species is unusually conductive); liana stem 180 ohm.m. NOTE this reverses the general pattern reported by Gora & Yanoviak (2015), who found lianas LESS resistive than temperate trees; we are treating Dipteryx as the outlier. PART E sweeps this ratio from 0.25 to 4 precisely because we do not trust it. 4. Lightning current: a "continuing current" phase of 200 A for 115 ms, which is the standard thermal-damage phase used by Gora et al. (2017), plus a 30 kA first return stroke reported separately. The return stroke is short and mostly flashes over; the continuing current does the cooking. 5. Grounding resistance is ignored (ohms, versus kilohms-to-megohms of stem). 6. Fitness model: seed output proportional to crown area, liana load and canopy gap; no seed predation, no masting, no density dependence, no genetics, no pollinators. It is a cartoon with the right shape. Run: python lightning-tree.py Deps: numpy (2.x). Falls back to stdlib if numpy is unavailable. Seed: fixed at 20250326 (the paper's publication date) for reproducibility. """ from __future__ import annotations import math try: import numpy as np HAVE_NUMPY = True except Exception: # pragma: no cover HAVE_NUMPY = False SEED = 20250326 # ============================================================================ # helpers # ============================================================================ def banner(title: str) -> None: print() print("=" * 78) print(title) print("=" * 78) def rule(title: str = "") -> None: print() if title: print("-- " + title + " " + "-" * max(0, 74 - len(title))) else: print("-" * 78) def eng(x: float, unit: str = "", sig: int = 4) -> str: """Human-readable engineering notation.""" if x == 0: return "0 " + unit prefixes = [(1e12, "T"), (1e9, "G"), (1e6, "M"), (1e3, "k"), (1.0, ""), (1e-3, "m"), (1e-6, "u"), (1e-9, "n")] a = abs(x) for scale, p in prefixes: if a >= scale: return f"{x/scale:.{sig}g} {p}{unit}" return f"{x:.{sig}g} {unit}" # ============================================================================ # PART 0 — parameters # ============================================================================ banner("PART 0 - MODEL PARAMETERS (club's own choices; see docstring)") # --- geometry ------------------------------------------------------------- H_TREE = 40.0 # m, canopy height of a large Dipteryx oleifera L_TRUNK = 42.0 # m, crown-to-root current path (height + root path) DBH = 1.00 # m, trunk diameter at breast height SAPWOOD_FRAC = 0.60 # fraction of basal area that conducts D_LIANA = 0.040 # m, liana stem diameter (4 cm: a big one) L_LIANA = 45.0 # m, liana stem path length (it wanders) # --- material ------------------------------------------------------------- RHO_E_TRUNK = 120.0 # ohm.m, Dipteryx sapwood (LOW end, deliberately) RHO_E_LIANA = 180.0 # ohm.m, liana stem RHO_D_TRUNK = 900.0 # kg/m3, fresh dense tropical hardwood RHO_D_LIANA = 700.0 # kg/m3, fresh liana stem C_TRUNK = 2000.0 # J/(kg.K), wet wood, ~45% moisture C_LIANA = 3000.0 # J/(kg.K), wetter liana tissue, ~60% moisture W_TRUNK = 0.45 # water mass fraction W_LIANA = 0.60 # water mass fraction L_VAP = 2.26e6 # J/kg, latent heat of vaporisation of water T_AMBIENT = 25.0 # degC T_CAMBIUM_DEATH = 60.0 # degC, lethal for living cambium T_BOIL = 100.0 # degC # --- lightning ------------------------------------------------------------ I_CONT = 200.0 # A, continuing-current phase T_CONT = 0.115 # s, its duration I_PEAK = 30e3 # A, first return-stroke peak S_RETURN = 5.5e4 # A^2.s, median first-stroke action integral A_TRUNK = SAPWOOD_FRAC * math.pi * (DBH / 2.0) ** 2 A_LIANA = math.pi * (D_LIANA / 2.0) ** 2 print(f" trunk: H={H_TREE:.0f} m, path L={L_TRUNK:.0f} m, DBH={DBH*100:.0f} cm") print(f" conducting area A = {A_TRUNK:.4f} m^2 ({SAPWOOD_FRAC:.0%} of basal area)") print(f" rho_e = {RHO_E_TRUNK:.0f} ohm.m, rho_d = {RHO_D_TRUNK:.0f} kg/m^3, c = {C_TRUNK:.0f} J/kg/K") print(f" liana: d={D_LIANA*100:.0f} cm, path L={L_LIANA:.0f} m") print(f" conducting area A = {A_LIANA:.3e} m^2") print(f" rho_e = {RHO_E_LIANA:.0f} ohm.m, rho_d = {RHO_D_LIANA:.0f} kg/m^3, c = {C_LIANA:.0f} J/kg/K") print(f" flash: continuing current {I_CONT:.0f} A for {T_CONT*1e3:.0f} ms; " f"return stroke peak {I_PEAK/1e3:.0f} kA") print(f" area ratio A_trunk / A_liana = {A_TRUNK/A_LIANA:,.0f}") # ============================================================================ # PART A — resistances and the two limiting cases # ============================================================================ banner("PART A - THE CIRCUIT: TWO RESISTORS, TWO DIFFERENT STORIES") R_TRUNK = RHO_E_TRUNK * L_TRUNK / A_TRUNK R_LIANA = RHO_E_LIANA * L_LIANA / A_LIANA G_TRUNK = 1.0 / R_TRUNK G_LIANA = 1.0 / R_LIANA M_TRUNK = RHO_D_TRUNK * L_TRUNK * A_TRUNK M_LIANA = RHO_D_LIANA * L_LIANA * A_LIANA rule("resistance and mass of each path") print(f" R_trunk = rho L / A = {RHO_E_TRUNK:.0f} x {L_TRUNK:.0f} / {A_TRUNK:.4f} " f"= {eng(R_TRUNK,'ohm')}") print(f" R_liana = rho L / A = {RHO_E_LIANA:.0f} x {L_LIANA:.0f} / {A_LIANA:.3e} " f"= {eng(R_LIANA,'ohm')}") print(f" R_liana / R_trunk = {R_LIANA/R_TRUNK:,.0f}") print(f" m_trunk = {M_TRUNK:,.0f} kg ({M_TRUNK/1000:.1f} tonnes of wood in the path)") print(f" m_liana = {M_LIANA:,.1f} kg") print(f" m_trunk / m_liana = {M_TRUNK/M_LIANA:,.0f}") # --- Limit A1: shared voltage (true parallel divider) ---------------------- rule("LIMIT A1 - shared voltage: the two paths span the same crown-to-ground gap") G_TOT = G_TRUNK + G_LIANA f_liana = G_LIANA / G_TOT I_l_A1 = I_CONT * f_liana I_t_A1 = I_CONT * (1.0 - f_liana) P_t_A1 = I_t_A1 ** 2 * R_TRUNK P_l_A1 = I_l_A1 ** 2 * R_LIANA pm_t_A1 = P_t_A1 / M_TRUNK pm_l_A1 = P_l_A1 / M_LIANA print(f" conductances: G_trunk = {eng(G_TRUNK,'S')}, G_liana = {eng(G_LIANA,'S')}") print(f" current fraction taken by the liana = {f_liana:.6f} ({f_liana*100:.4f} %)") print(f" I_trunk = {I_t_A1:.3f} A I_liana = {I_l_A1:.4f} A") print(f" P_trunk = I^2 R = {eng(P_t_A1,'W')} P_liana = {eng(P_l_A1,'W')}") print(f" power per kg: trunk {pm_t_A1:,.0f} W/kg liana {pm_l_A1:,.0f} W/kg") print(f" >>> liana / trunk heating rate per kg = {pm_l_A1/pm_t_A1:.3f}") dT_t_A1 = pm_t_A1 * T_CONT / C_TRUNK dT_l_A1 = pm_l_A1 * T_CONT / C_LIANA print(f" temperature rise over the {T_CONT*1e3:.0f} ms continuing current:") print(f" trunk dT = {dT_t_A1:.2f} K -> {T_AMBIENT + dT_t_A1:.1f} degC") print(f" liana dT = {dT_l_A1:.2f} K -> {T_AMBIENT + dT_l_A1:.1f} degC") print(" VERDICT: in this limit NOTHING is damaged. The model fails to") print(" reproduce the paper's observation. That failure is diagnostic.") # --- Limit A2: shared current -------------------------------------------- rule("LIMIT A2 - shared current: the liana carries the full continuing current") P_t_A2 = I_CONT ** 2 * R_TRUNK P_l_A2 = I_CONT ** 2 * R_LIANA pm_t_A2 = P_t_A2 / M_TRUNK pm_l_A2 = P_l_A2 / M_LIANA print(f" same I = {I_CONT:.0f} A through each path, separately:") print(f" P_trunk = {eng(P_t_A2,'W')} -> {pm_t_A2:,.0f} W/kg") print(f" P_liana = {eng(P_l_A2,'W')} -> {eng(pm_l_A2,'W/kg')}") ratio_A2 = pm_l_A2 / pm_t_A2 print(f" >>> liana / trunk heating rate per kg = {ratio_A2:,.0f}") # closed form check closed = (RHO_E_LIANA / RHO_E_TRUNK) * (A_TRUNK / A_LIANA) ** 2 * (RHO_D_TRUNK / RHO_D_LIANA) print(f" closed form (rho_l/rho_t)(A_t/A_l)^2(rho_dt/rho_dl) = {closed:,.0f} [check]") print(f" bracket spanned by the two limits: {pm_l_A1/pm_t_A1:.2f}x ... {ratio_A2:,.0f}x") print(f" that is a factor of {ratio_A2/(pm_l_A1/pm_t_A1):,.0f} of pure modelling ignorance.") # ============================================================================ # PART B — the cancellation that ruined our afternoon # ============================================================================ banner("PART B - WHY LIMIT A1 IS BLIND TO THICKNESS (exact cancellation)") print(""" For a path of length L, area A, resistivity rho_e, density rho_d holding a voltage V across its ends: R = rho_e L / A P = V^2 / R = V^2 A / (rho_e L) m = rho_d L A P/m = V^2 / (rho_e rho_d L^2) <-- A has cancelled exactly Heating per kilogram in a shared-voltage divider depends ONLY on material and path length. Making a stem thinner gives it less current AND less mass to spread that current over, in precisely equal measure. The 4 cm vine and the 1 m trunk are, to this model, the same object.""") V_shared = I_CONT / G_TOT pm_t_id = V_shared ** 2 / (RHO_E_TRUNK * RHO_D_TRUNK * L_TRUNK ** 2) pm_l_id = V_shared ** 2 / (RHO_E_LIANA * RHO_D_LIANA * L_LIANA ** 2) print(f"\n numerical check at V = {eng(V_shared,'V')}:") print(f" identity trunk {pm_t_id:,.1f} W/kg liana {pm_l_id:,.1f} W/kg") print(f" simulation trunk {pm_t_A1:,.1f} W/kg liana {pm_l_A1:,.1f} W/kg") print(f" predicted ratio = (rho_t rho_dt L_t^2)/(rho_l rho_dl L_l^2) = " f"{(RHO_E_TRUNK*RHO_D_TRUNK*L_TRUNK**2)/(RHO_E_LIANA*RHO_D_LIANA*L_LIANA**2):.4f}") print(f" simulated ratio = {pm_l_A1/pm_t_A1:.4f} [agrees]") print("\n A doubling of liana diameter changes this ratio by exactly 0.000.") for d in (0.01, 0.02, 0.04, 0.08, 0.16): a = math.pi * (d / 2) ** 2 r = RHO_E_LIANA * L_LIANA / a g = 1 / r ff = g / (G_TRUNK + g) il = I_CONT * ff mm = RHO_D_LIANA * L_LIANA * a pmv = (il ** 2 * r) / mm print(f" liana d = {d*100:5.1f} cm -> I = {il:9.4f} A, m = {mm:8.2f} kg, " f"P/m = {pmv:,.1f} W/kg") # ============================================================================ # PART C — heating, boiling, and time-to-destruction # ============================================================================ banner("PART C - HEATING: ENERGY IN, TEMPERATURE OUT, TIME TO FAILURE") E_cambium_t = C_TRUNK * (T_CAMBIUM_DEATH - T_AMBIENT) E_cambium_l = C_LIANA * (T_CAMBIUM_DEATH - T_AMBIENT) E_boil_t = C_TRUNK * (T_BOIL - T_AMBIENT) + W_TRUNK * L_VAP E_boil_l = C_LIANA * (T_BOIL - T_AMBIENT) + W_LIANA * L_VAP rule("per-kilogram energy budgets of the two tissues") print(f" to kill cambium (reach {T_CAMBIUM_DEATH:.0f} degC):") print(f" trunk {E_cambium_t:,.0f} J/kg liana {E_cambium_l:,.0f} J/kg") print(f" to reach {T_BOIL:.0f} degC AND flash all tissue water to steam:") print(f" trunk {E_boil_t:,.0f} J/kg liana {E_boil_l:,.0f} J/kg") rule("time to cross each threshold, at the heating rates computed above") cases = [ ("trunk, shared-voltage limit", pm_t_A1, E_cambium_t, E_boil_t), ("liana, shared-voltage limit", pm_l_A1, E_cambium_l, E_boil_l), ("trunk, full 200 A", pm_t_A2, E_cambium_t, E_boil_t), ("liana, full 200 A", pm_l_A2, E_cambium_l, E_boil_l), ] print(f" {'case':32s} {'t to 60 degC':>16s} {'t to steam-out':>18s}") for name, pm, ec, eb in cases: tc = ec / pm tb = eb / pm print(f" {name:32s} {eng(tc,'s'):>16s} {eng(tb,'s'):>18s}") print(f"\n the continuing current lasts {T_CONT*1e3:.0f} ms = {T_CONT:.3f} s.") t_boil_liana_A2 = E_boil_l / pm_l_A2 print(f" >>> at full current the liana reaches total steam-out in " f"{eng(t_boil_liana_A2,'s')},") print(f" i.e. in the first {t_boil_liana_A2/T_CONT*100:.4f}% of the flash.") t_boil_trunk_A2 = E_boil_t / pm_t_A2 print(f" >>> the trunk at the SAME current would need {t_boil_trunk_A2:.1f} s,") print(f" which is {t_boil_trunk_A2/T_CONT:,.0f}x longer than the flash lasts.") rule("how much current does each actually need, to be destroyed?") def threshold_current(R, m, E_per_kg, dt): """Smallest steady current that delivers E_per_kg to mass m in time dt.""" return math.sqrt(E_per_kg * m / (R * dt)) I_crit_liana = threshold_current(R_LIANA, M_LIANA, E_boil_l, T_CONT) I_crit_trunk = threshold_current(R_TRUNK, M_TRUNK, E_boil_t, T_CONT) I_crit_liana_camb = threshold_current(R_LIANA, M_LIANA, E_cambium_l, T_CONT) print(f" liana, full steam-out : I_crit = {I_crit_liana:8.2f} A " f"({I_crit_liana/I_CONT*100:.1f}% of the continuing current)") print(f" liana, cambium death : I_crit = {I_crit_liana_camb:8.2f} A " f"({I_crit_liana_camb/I_CONT*100:.1f}% of it)") print(f" trunk, full steam-out : I_crit = {I_crit_trunk:8.0f} A " f"({I_crit_trunk/I_CONT:,.0f}x the continuing current)") print(f" >>> safety ratio in current, trunk vs liana = {I_crit_trunk/I_crit_liana:,.0f}x") print(f" >>> safety ratio in energy , trunk vs liana = " f"{(I_crit_trunk/I_crit_liana)**2:,.0f}x") rule("the return stroke, for scale") E_rs_trunk = S_RETURN * R_TRUNK print(f" action integral of a median first return stroke S = {S_RETURN:,.0f} A^2.s") print(f" if the trunk conducted ALL of it: E = S x R = {eng(E_rs_trunk,'J')}") print(f" -> dT = {E_rs_trunk/(M_TRUNK*C_TRUNK):.1f} K spread over {M_TRUNK/1000:.1f} tonnes") print(" (In reality most of the return stroke flashes over the wet bark surface.") print(" Our model has no flashover in it; this is its largest single omission.)") # ============================================================================ # PART D — fitness model # ============================================================================ banner("PART D - FITNESS: 300-YEAR MONTE CARLO, TOLERATE vs AVOID") # strategies ----------------------------------------------------------------- # D1: identical morphology, differing only in whether a strike kills. # `sublethal` = fractional loss of seed output in the year after a non-fatal # strike, decaying with SUBLETHAL_TAU. The paper reports 5.7x more leaf loss # in struck non-Dipteryx than in struck Dipteryx, hence 0.08 vs 0.46. TOLERANT = dict(name="tolerant (Dipteryx-like)", H=40.0, crown_r=12.0, p_die=0.05, sublethal=0.08) INTOLERANT = dict(name="intolerant, same shape", H=40.0, crown_r=12.0, p_die=0.64, sublethal=0.46) # D2: the evolutionary trade-off. The avoider is shorter and narrower, so it # is struck less often, but it also intercepts less light. PRONE = TOLERANT AVOID = dict(name="lightning-avoiding neighbour", H=36.0, crown_r=10.6, p_die=0.64, sublethal=0.46) SUBLETHAL_TAU = 8.0 # yr, recovery time of a scorched crown LAMBDA_PRONE = 1.0 / 56.0 # paper: mean interval between strikes ~56 yr YEARS = 300 N_TREES = 20000 MATURITY = 60 # yr, first reproduction RAMP = 40 # yr, to full reproductive capacity H_BACKGROUND = 0.0060 # /yr background mortality of a canopy tree LIANA_MORT_COEF = 0.90 # lianas raise mortality by up to this fraction LIANA_FEC_COEF = 0.70 # lianas cut seed output by up to this fraction LIANA_R = 0.050 # /yr logistic regrowth rate of liana load LIANA_KNOCKDOWN = 0.22 # paper: 78% of lianas killed per strike GAP_BONUS = 0.50 # seed bonus immediately after a strike clears 2.1 Mg GAP_TAU = 22.0 # yr, decay of that bonus as neighbours refill def strike_rate(H, crown_r): """Hazard scales as height^2 x crown area, normalised to the prone tree.""" A_ref = math.pi * PRONE["crown_r"] ** 2 return LAMBDA_PRONE * (H / PRONE["H"]) ** 2 * (math.pi * crown_r ** 2) / A_ref lam_p = strike_rate(PRONE["H"], PRONE["crown_r"]) lam_a = strike_rate(AVOID["H"], AVOID["crown_r"]) print(f" strike hazard, prone = {lam_p:.6f} /yr (1 per {1/lam_p:.0f} yr)") print(f" strike hazard, avoiding = {lam_a:.6f} /yr (1 per {1/lam_a:.0f} yr)") print(f" the prone morphology is struck {(lam_p/lam_a - 1)*100:.0f}% more often") print(f" (paper reports 49-68% more often for D. oleifera proportions)") print(f" crown area ratio = {(PRONE['crown_r']/AVOID['crown_r'])**2:.3f}") def simulate(strategy, n=N_TREES, years=YEARS, rng=None): """Vectorised cohort simulation. Returns dict of summary arrays.""" lam = strike_rate(strategy["H"], strategy["crown_r"]) p_die = strategy["p_die"] crown_area = math.pi * strategy["crown_r"] ** 2 crown_rel = crown_area / (math.pi * PRONE["crown_r"] ** 2) alive = np.ones(n, dtype=bool) liana = np.full(n, 0.05) # liana load, 0..1 gap = np.zeros(n) # canopy-gap bonus, 0..1 seeds = np.zeros(n) n_strikes = np.zeros(n, dtype=np.int64) death_year = np.full(n, years, dtype=np.int64) surv_curve = np.zeros(years + 1) seed_curve = np.zeros(years + 1) liana_curve = np.zeros(years + 1) surv_curve[0] = 1.0 for t in range(1, years + 1): # --- liana dynamics (logistic regrowth toward full infestation) ---- liana = liana + LIANA_R * liana * (1.0 - liana) gap = gap * math.exp(-1.0 / GAP_TAU) # --- lightning ---------------------------------------------------- struck = alive & (rng.random(n) < lam) n_strikes += struck liana = np.where(struck, liana * LIANA_KNOCKDOWN, liana) gap = np.where(struck, 1.0, gap) killed_by_bolt = struck & (rng.random(n) < p_die) # --- background + liana-driven mortality --------------------------- h = H_BACKGROUND * (1.0 + LIANA_MORT_COEF * liana) died = alive & ((rng.random(n) < h) | killed_by_bolt) death_year = np.where(died & alive, t, death_year) alive = alive & ~died # --- reproduction -------------------------------------------------- if t > MATURITY: size_f = min(1.0, (t - MATURITY) / RAMP) out = (crown_rel * size_f * (1.0 - LIANA_FEC_COEF * liana) * (1.0 + GAP_BONUS * gap)) seeds += np.where(alive, np.maximum(out, 0.0), 0.0) seed_curve[t] = np.mean(np.where(alive, np.maximum(out, 0.0), 0.0)) surv_curve[t] = alive.mean() liana_curve[t] = liana[alive].mean() if alive.any() else np.nan return dict(seeds=seeds, strikes=n_strikes, death_year=death_year, surv=surv_curve, seedrate=seed_curve, lianabar=liana_curve, lam=lam) if HAVE_NUMPY: res_t = simulate(TOLERANT, rng=np.random.default_rng(SEED)) res_i = simulate(INTOLERANT, rng=np.random.default_rng(SEED + 1)) res_p = res_t res_a = simulate(AVOID, rng=np.random.default_rng(SEED + 2)) mt, mi = res_t["seeds"].mean(), res_i["seeds"].mean() mp, ma = res_p["seeds"].mean(), res_a["seeds"].mean() RATIO = mt / mi # D1 - the paper's contrast RATIO_TRADEOFF = mp / ma # D2 - the trade-off contrast rule("D1 - TOLERANCE ALONE (identical trees, 20,000 each)") print(f" {'':34s} {'tolerant':>14s} {'intolerant':>14s}") print(f" {'mean lifetime seed output (a.u.)':34s} {mt:>14.2f} {mi:>14.2f}") print(f" {'median age at death (yr)':34s} " f"{np.median(res_t['death_year']):>14.0f} {np.median(res_i['death_year']):>14.0f}") print(f" {'fraction alive at 300 yr':34s} " f"{res_t['surv'][-1]:>14.3f} {res_i['surv'][-1]:>14.3f}") print(f" {'mean strikes received':34s} " f"{res_t['strikes'].mean():>14.2f} {res_i['strikes'].mean():>14.2f}") print(f" {'mean liana load, live trees':34s} " f"{np.nanmean(res_t['lianabar'][100:]):>14.3f} " f"{np.nanmean(res_i['lianabar'][100:]):>14.3f}") print() print(f" >>> LIFETIME SEED OUTPUT RATIO (tolerant / intolerant) = {RATIO:.2f}x") print(f" paper's reported figure: ~14x") print(f" discrepancy factor: {max(RATIO,14)/min(RATIO,14):.2f} " f"({'below' if RATIO < 14 else 'above'})") rule("D2 - THE TRADE-OFF (tall+tolerant vs shorter+intolerant)") print(f" mean lifetime seed output: prone {mp:.2f} avoiding {ma:.2f}") print(f" >>> RATIO = {RATIO_TRADEOFF:.2f}x") print(" Smaller, as it should be: the avoider genuinely buys lower strike") print(" risk with its shorter stature. Growing tall is not free.") # decomposition ---------------------------------------------------------- rule("where the D1 advantage comes from (one factor off at a time)") def variant(a_strat=None, b_strat=None, n=8000, **over): global LIANA_FEC_COEF, GAP_BONUS, LIANA_MORT_COEF keep = (LIANA_FEC_COEF, GAP_BONUS, LIANA_MORT_COEF) LIANA_FEC_COEF = over.get("fec", keep[0]) GAP_BONUS = over.get("gap", keep[1]) LIANA_MORT_COEF = over.get("mort", keep[2]) A = a_strat or TOLERANT B = b_strat or INTOLERANT if "p_die_b" in over: B = dict(B, p_die=over["p_die_b"]) p = simulate(A, n=n, rng=np.random.default_rng(SEED + 7)) b = simulate(B, n=n, rng=np.random.default_rng(SEED + 8)) LIANA_FEC_COEF, GAP_BONUS, LIANA_MORT_COEF = keep return p["seeds"].mean() / b["seeds"].mean() print(f" full model (8,000 trees) {variant():6.2f}x") print(f" no liana effect on seed output (fec=0) {variant(fec=0.0):6.2f}x") print(f" no canopy-gap bonus (gap=0) {variant(gap=0.0):6.2f}x") print(f" no liana effect on mortality (mort=0) {variant(mort=0.0):6.2f}x") print(f" none of the three " f"{variant(fec=0.0, gap=0.0, mort=0.0):6.2f}x") print(" (what survives the last line is pure survival.)") # sensitivity ------------------------------------------------------------ rule("sensitivity of D1: is 14x inside our plausible range?") print(f" {'p(die | struck)':>16s} {'liana fec.':>12s} {'gap bonus':>12s} {'ratio':>10s}") lo = hi = None grid = [] for pd in (0.56, 0.64, 0.78): for fec in (0.40, 0.70, 0.85): for gp in (0.25, 0.50, 0.80): r = variant(fec=fec, gap=gp, p_die_b=pd) grid.append(r) lo = r if lo is None else min(lo, r) hi = r if hi is None else max(hi, r) print(f" {pd:>16.2f} {fec:>12.2f} {gp:>12.2f} {r:>10.2f}x") print(f"\n range across that grid: {lo:.2f}x to {hi:.2f}x " f"(median {sorted(grid)[len(grid)//2]:.2f}x)") print(f" the paper's 14x {'IS' if lo <= 14 <= hi else 'IS NOT'} inside it.") # ------------------------------------------------------------------ figs banner("FIGURE DATA (copied by hand into the article's inline SVGs)") print("\n[FIG2] heating rate per kilogram, W/kg, log scale") fig2 = [ ("trunk, shared-voltage", pm_t_A1), ("liana, shared-voltage", pm_l_A1), ("liana at 5% of current", (0.05 * I_CONT) ** 2 * R_LIANA / M_LIANA), ("liana at 25% of current", (0.25 * I_CONT) ** 2 * R_LIANA / M_LIANA), ("liana, full 200 A", pm_l_A2), ] for nm, v in fig2: print(f" {nm:26s} {v:>16,.1f} log10 = {math.log10(v):5.2f}") print(f" destruction threshold (liana) {E_boil_l/T_CONT:>12,.0f} W/kg " f"log10 = {math.log10(E_boil_l/T_CONT):5.2f}") print(f" destruction threshold (trunk) {E_boil_t/T_CONT:>12,.0f} W/kg " f"log10 = {math.log10(E_boil_t/T_CONT):5.2f}") print("\n[FIG3] temperature vs time at full 200 A (deg C), log-time") for tt in (1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.115, 1.0, 10.0, 100.0): Tl = T_AMBIENT + min(pm_l_A2 * tt, E_boil_l) / C_LIANA Tt = T_AMBIENT + min(pm_t_A2 * tt, E_boil_t) / C_TRUNK Tl = min(Tl, 100.0 + 1e-9) if pm_l_A2 * tt > C_LIANA * 75 else Tl Tt = min(Tt, 100.0 + 1e-9) if pm_t_A2 * tt > C_TRUNK * 75 else Tt print(f" t = {eng(tt,'s'):>10s} liana {Tl:9.1f} C trunk {Tt:9.1f} C") print("\n[FIG4] survivorship and mean liana load, D1 (every 20 yr)") print(f" {'yr':>5s} {'S_tol':>9s} {'S_intol':>9s} {'L_tol':>9s} {'L_intol':>9s}") for t in range(0, YEARS + 1, 20): print(f" {t:>5d} {res_t['surv'][t]:>9.3f} {res_i['surv'][t]:>9.3f} " f"{res_t['lianabar'][t]:>9.3f} {res_i['lianabar'][t]:>9.3f}") print("\n[FIG4b] cumulative mean seed output, D1 (every 20 yr)") cum_t = np.cumsum(res_t["seedrate"]) cum_i = np.cumsum(res_i["seedrate"]) for t in range(60, YEARS + 1, 20): print(f" {t:>5d} tolerant {cum_t[t]:>8.2f} intolerant {cum_i[t]:>8.2f} " f"ratio {cum_t[t]/max(cum_i[t],1e-9):>6.2f}") else: # pragma: no cover RATIO = float("nan") RATIO_TRADEOFF = float("nan") print(" numpy unavailable - fitness Monte Carlo skipped.") # ============================================================================ # PART E — sensitivity of the circuit model # ============================================================================ banner("PART E - HOW MUCH DO OUR RESISTIVITY GUESSES MATTER?") print(f" {'rho_liana/rho_trunk':>22s} {'A1 ratio':>12s} {'A2 ratio':>14s}") for k in (0.25, 0.5, 1.0, 1.5, 2.0, 4.0): rl = RHO_E_TRUNK * k Rl = rl * L_LIANA / A_LIANA Gl = 1 / Rl ff = Gl / (G_TRUNK + Gl) a1 = ((I_CONT * ff) ** 2 * Rl / M_LIANA) / pm_t_A1 a2 = (I_CONT ** 2 * Rl / M_LIANA) / pm_t_A2 print(f" {k:>22.2f} {a1:>12.3f} {a2:>14,.0f}") print(" The A2 column moves by a factor of 16 across the whole sweep;") print(" the A1 column never leaves the neighbourhood of 1. The qualitative") print(" conclusion - that thickness only matters in the shared-current") print(" limit - is robust to everything we do not know about resistivity.") # ============================================================================ # SUMMARY # ============================================================================ banner("KEY NUMBERS USED IN THE ARTICLE") rows = [ ("trunk conducting cross-section", f"{A_TRUNK:.4f} m^2"), ("liana conducting cross-section", f"{A_LIANA:.3e} m^2"), ("area ratio", f"{A_TRUNK/A_LIANA:,.0f}x"), ("R_trunk", eng(R_TRUNK, "ohm")), ("R_liana", eng(R_LIANA, "ohm")), ("mass in trunk path", f"{M_TRUNK/1000:.1f} t"), ("mass in liana path", f"{M_LIANA:.1f} kg"), ("liana share of current (shared-V)", f"{f_liana*100:.4f} %"), ("heating ratio, shared-voltage limit", f"{pm_l_A1/pm_t_A1:.2f}x"), ("heating ratio, shared-current limit", f"{ratio_A2:,.0f}x"), ("trunk dT over 115 ms (shared-V)", f"{dT_t_A1:.2f} K"), ("liana dT over 115 ms (shared-V)", f"{dT_l_A1:.2f} K"), ("liana time to steam-out at 200 A", eng(t_boil_liana_A2, "s")), ("trunk time to steam-out at 200 A", f"{t_boil_trunk_A2:.0f} s"), ("current that destroys the liana", f"{I_crit_liana:.1f} A"), ("current that destroys the trunk", f"{I_crit_trunk:,.0f} A"), ("trunk safety margin in energy", f"{(I_crit_trunk/I_crit_liana)**2:,.0f}x"), ("fitness ratio D1, our Monte Carlo", f"{RATIO:.2f}x"), ("fitness ratio D2 (trade-off)", f"{RATIO_TRADEOFF:.2f}x"), ("fitness ratio, the paper", "~14x"), ] for k, v in rows: print(f" {k:<38s} {v:>20s}") print() print(" Reminder: every number above is the club's own model output, not a") print(" result from Gora et al. (2025). The paper measured outcomes in a") print(" forest; we only checked that the physics does not forbid them.") print()