""" zombie-volcano.py -- the Science Journaling Club's own simplified model for "A Volcano That Has Been Dead for 250,000 Years Is Still Breathing" (Field Note, Paper Analysis, Volcanology). THIS IS NOT THE PAPER'S ANALYSIS. Liu et al. (2025, PNAS 122(18) e2420996122, doi:10.1073/pnas.2420996122) did the seismic tomography, the anisotropy inversion and the petrophysical modelling of Uturuncu's interior. Nothing in this file touches their data. What follows is a deformation-side exercise the club ran on its own, on synthetic ground-motion data, to answer a narrower question: what can the SHAPE of the surface deformation alone tell you, and what can it not tell you? ============================================================================= PART A -- THE MOGI POINT SOURCE ============================================================================= Mogi (1958) solved the surface displacement produced by a small pressurised sphere buried in an elastic half-space. For a source at depth d below a flat free surface, undergoing a volume change dV, the surface displacements at radial distance r from the point directly above the source are uz(r) = (1 - nu) * dV * d / ( pi * (r^2 + d^2)^(3/2) ) [vertical] ur(r) = (1 - nu) * dV * r / ( pi * (r^2 + d^2)^(3/2) ) [radial] with nu Poisson's ratio (0.25 here, the standard crustal value). Three properties of this solution do all the work in this article: P1 uz(r) has the SAME SIGN at every r. The factor d/(r^2+d^2)^(3/2) is strictly positive for d > 0 and all finite r. So sign(uz) = sign(dV), everywhere, with no exceptions. An inflating point source lifts the whole surface; a deflating one lowers the whole surface. NEITHER can produce a ring of subsidence around a central uplift. Any observed change of sign in the vertical field therefore requires at least two sources of opposite sign. This is a theorem about the model, not a fitting result, and the fits below only measure how badly a single source fails. P2 ur(r) is maximised at r = d / sqrt(2) = 0.7071 d. Differentiate r*(r^2+d^2)^(-3/2): the stationary point is where d^2 = 2 r^2. So the radius of maximum horizontal motion reads the source depth directly. P3 uz(r)/uz(0) = d^3 / (r^2+d^2)^(3/2), so the half-width, the radius at which the uplift has fallen to half its peak, is r_half = d * sqrt(2^(2/3) - 1) = 0.7664 d. A wide uplift bump means a deep source. This is why a 70 km wide uplift cannot come from a shallow chamber, no matter how hard anyone squints at it. ASSUMPTIONS AND LIMITS (A1-A8), stated plainly: A1 Elastic, homogeneous, isotropic half-space. The real crust under the Altiplano-Puna is layered, hot, and partly molten, and a viscoelastic or temperature-dependent rheology changes the inferred depths by tens of percent and the inferred volumes by more. A2 Point source. A real magma body 200 km across is not a point. The Mogi solution is only accurate when the source radius is much less than its depth, which is badly violated by the Altiplano-Puna Magma Body. A3 Flat free surface. Uturuncu's summit stands 6,008 m above sea level and about 1.7 km above the surrounding altiplano. Topography matters and is ignored here. A4 Axial symmetry. We model a 1-D radial profile, not the full 2-D field that a real InSAR inversion fits. A5 The "observations" are SYNTHETIC. We build a sombrero-shaped profile that resembles the published Uturuncu pattern (about 10 mm/yr of central uplift over a dome a few tens of kilometres wide, ringed by roughly 1 to 2 mm/yr of subsidence, the whole field more than 100 km across; Pritchard & Simons 2002; Fialko & Pearse 2012; Henderson & Pritchard 2017), then add seeded Gaussian noise at 0.20 mm/yr, a plausible rate uncertainty for a multi-year stacked InSAR velocity field. We are NOT inverting real interferograms. The exercise tests what the shape implies, using data whose true answer we know. A6 Sources are co-axial: both sit on the same vertical line under the summit. Real deformation sources can be offset laterally. A7 Steady rates. Uturuncu's uplift rate has in fact slowed, from about 10 mm/yr in the 1990s to nearer 3 mm/yr recently, which by itself argues against a steadily filling reservoir. A8 dV is the volume change of the SOURCE, not the volume of material added to it. For a compressible fluid these differ by a large factor, and that difference is the whole point of Part D. ============================================================================= PART B -- THE INVERSION ============================================================================= Given a radial profile of observed vertical rate, we fit: One source: free depth d, free volume rate dV. Two sources: free depths d1 < d2, free volume rates dV1, dV2. For a FIXED set of depths the problem is linear in the volume rates, because uz is proportional to dV. So we grid-search the depths (which is the hard, non-linear part) and solve the volume rates exactly by least squares at every grid node. No optimiser, no starting guess, no local minimum to fall into. Misfit is reported as root-mean-square residual in mm/yr, and as variance reduction, VR = 1 - sum(res^2)/sum(obs^2), in percent. ============================================================================= PART C -- WHY THE SOMBRERO IS DIAGNOSTIC ============================================================================= Three separate quantities are reported: C1 The RMS misfit of the best single source against the best pair. C2 The sign test: the best single source's predicted displacement in the moat annulus, compared with the observed subsidence there. By P1 the prediction cannot be negative if the fitted dV is positive. C3 The moat-to-peak amplitude ratio, which a single source fixes at zero by construction and which the data put near 0.17. ============================================================================= PART D -- THE ERUPTION ARITHMETIC ============================================================================= Plain division. Take the fitted inflation rate of the shallow source and ask how long it would take to accumulate volumes the size of real eruptions, and compare it with Uturuncu's own long-term magma output rate. D1 Uturuncu's edifice volume is about 85 km3, erupted between roughly 890 ka and 250 ka (Sparks et al. 2008; Muir et al. 2015). That gives a long-run average magma output rate of about 1.3e5 m3/yr. D2 Reference eruption sizes, dense-rock equivalent: VEI 4 = 0.1 km3, VEI 5 = 1 km3, VEI 6 = 10 km3. The Altiplano-Puna Volcanic Complex ignimbrites are 500 to 2,500 km3 (Pastos Grandes, Atana). D3 Density of dacitic magma taken as 2,400 kg/m3. Density of supercritical H2O-CO2 fluid in the 2 to 15 km depth range taken as 200 to 400 kg/m3 (used only to show the mass contrast, an order-of-magnitude point). D4 The imaged melt fraction in the Altiplano-Puna Magma Body maxes out at about 25% (Liu et al. 2025). Crystal mushes lock up rheologically and stop being eruptible somewhere near 50% crystals, so extracting 1 km3 of eruptible melt from 25% mush means reorganising about 4 km3 of mush. That factor of 4 is a rule of thumb, not a measurement. Seeded RNG. Results are bit-for-bit reproducible on any machine. """ import math import random import numpy as np SEED = 20250428 # the paper's publication date, 28 April 2025 RNG = np.random.default_rng(SEED) random.seed(SEED) NU = 0.25 # Poisson's ratio SEP = "=" * 74 def sub(s): print("\n" + s + "\n" + "-" * len(s)) # ====================================================================== # PART A -- the Mogi point source # ====================================================================== def mogi_uz(r, d, dV, nu=NU): """Vertical surface displacement above a Mogi source. Metres.""" r = np.asarray(r, dtype=float) return (1.0 - nu) * dV * d / (math.pi * (r * r + d * d) ** 1.5) def mogi_ur(r, d, dV, nu=NU): """Radial (horizontal, outward) surface displacement. Metres.""" r = np.asarray(r, dtype=float) return (1.0 - nu) * dV * r / (math.pi * (r * r + d * d) ** 1.5) def unit_uz(r, d, nu=NU): """uz per unit volume change: the column of the design matrix.""" return mogi_uz(r, d, 1.0, nu) print(SEP) print("ZOMBIE VOLCANO -- CLUB MOGI-SOURCE DEFORMATION MODEL") print("Anchor: Liu, Kendall, Zhang, Blundy, Pritchard, Hudson & MacQueen") print(" (2025) PNAS 122(18), e2420996122") print(" doi:10.1073/pnas.2420996122") print("This is the club's own simplified model on SYNTHETIC data,") print("NOT the paper's analysis. The paper is a seismology paper.") print(SEP) sub("PART A0 -- THE THREE PROPERTIES OF A MOGI SOURCE") print(" uz(r) = (1-nu) dV d / (pi (r^2+d^2)^1.5) [vertical]") print(" ur(r) = (1-nu) dV r / (pi (r^2+d^2)^1.5) [radial]") print(f" nu = {NU}") print() print(" P1 sign(uz) = sign(dV) at EVERY radius. A single point source") print(" cannot change the sign of the vertical displacement anywhere.") print(" Demonstration, one inflating source at d = 20 km, dV = +1e7 m3:") _r_demo = np.array([0.0, 10e3, 30e3, 60e3, 100e3, 200e3, 500e3]) _uz_demo = mogi_uz(_r_demo, 20e3, 1e7) * 1000.0 print(" r (km) : " + " ".join(f"{v/1e3:>8.0f}" for v in _r_demo)) print(" uz (mm): " + " ".join(f"{v:>8.4f}" for v in _uz_demo)) print(f" minimum over that range = {_uz_demo.min():.6f} mm " f"(sign never flips)") print() r_ur_peak = 20e3 / math.sqrt(2.0) print(" P2 ur is maximal at r = d/sqrt(2) = 0.7071 d.") print(f" For d = 20 km that is r = {r_ur_peak/1e3:.2f} km.") _rr = np.linspace(1.0, 60e3, 60001) _num_peak = _rr[np.argmax(mogi_ur(_rr, 20e3, 1e7))] print(f" Numerical check by brute-force scan: {_num_peak/1e3:.2f} km.") print() half_coeff = math.sqrt(2.0 ** (2.0 / 3.0) - 1.0) print(f" P3 uz half-width r_half = d * sqrt(2^(2/3) - 1) = {half_coeff:.4f} d.") print(f" An uplift bump 70 km across (r_half = 35 km) needs") print(f" d = 35 / {half_coeff:.4f} = {35.0/half_coeff:.1f} km. Deep.") print(" Ratio uz/ur = d/r at every radius, so a single pair of") print(" vertical and horizontal measurements already fixes the depth.") # ====================================================================== # PART B -- build the synthetic sombrero # ====================================================================== sub("PART B -- THE SYNTHETIC SOMBRERO (our fabricated 'observations')") R_MAX = 120e3 R_STEP = 2e3 r = np.arange(0.0, R_MAX + R_STEP, R_STEP) # 0 to 120 km, 2 km spacing # "Truth" model. Chosen so the profile resembles the published Uturuncu # pattern; NOT taken from any paper's inversion. TRUE_D1 = 20e3 # shallow inflating source, m TRUE_D2 = 80e3 # deep deflating source, m TRUE_UZ1_0 = 0.014 # its contribution to peak uplift, m/yr TRUE_UZ2_0 = -0.004 # its contribution at r = 0, m/yr # convert those peak amplitudes back to volume rates TRUE_DV1 = TRUE_UZ1_0 * math.pi * TRUE_D1 ** 2 / (1.0 - NU) TRUE_DV2 = TRUE_UZ2_0 * math.pi * TRUE_D2 ** 2 / (1.0 - NU) signal = mogi_uz(r, TRUE_D1, TRUE_DV1) + mogi_uz(r, TRUE_D2, TRUE_DV2) NOISE_SD = 0.20e-3 # 0.20 mm/yr, a plausible stacked-InSAR rate error noise = RNG.normal(0.0, NOISE_SD, size=r.size) obs = signal + noise obs_mm = obs * 1000.0 sig_mm = signal * 1000.0 print(f" radial profile: {r.size} samples, 0 to {R_MAX/1e3:.0f} km, " f"{R_STEP/1e3:.0f} km spacing") print(f" seeded noise : Gaussian, sd = {NOISE_SD*1000:.2f} mm/yr, " f"seed = {SEED}") print() print(" TRUE (hidden) model used to manufacture the data:") print(f" source 1 depth {TRUE_D1/1e3:>6.1f} km " f"dV = {TRUE_DV1:+.4e} m3/yr (inflating)") print(f" source 2 depth {TRUE_D2/1e3:>6.1f} km " f"dV = {TRUE_DV2:+.4e} m3/yr (deflating)") print(f" |dV2|/dV1 = {abs(TRUE_DV2)/TRUE_DV1:.2f}") print(f" net dV = {TRUE_DV1 + TRUE_DV2:+.4e} m3/yr") print() peak_mm = obs_mm[0] moat_idx = int(np.argmin(obs_mm)) moat_r = r[moat_idx] / 1e3 moat_mm = obs_mm[moat_idx] # zero crossing, linear interpolation on the noiseless signal cross = None for i in range(len(sig_mm) - 1): if sig_mm[i] > 0.0 >= sig_mm[i + 1]: t = sig_mm[i] / (sig_mm[i] - sig_mm[i + 1]) cross = (r[i] + t * R_STEP) / 1e3 break print(" Shape of the resulting 'observed' profile:") print(f" peak uplift at r = 0 {peak_mm:+.2f} mm/yr") print(f" zero crossing {cross:.1f} km") print(f" deepest subsidence {moat_mm:+.2f} mm/yr " f"at r = {moat_r:.0f} km") print(f" moat-to-peak amplitude ratio {abs(moat_mm)/peak_mm:.3f}") print(f" still subsiding at r = 120 km {obs_mm[-1]:+.2f} mm/yr") print() print(" For comparison, the published Uturuncu pattern: roughly 10 mm/yr of") print(" central uplift in the 1990s over a dome tens of km wide, ringed by") print(" 1 to 4 mm/yr of subsidence reaching past 70 km radius (Pritchard &") print(" Simons 2002; Fialko & Pearse 2012; Henderson & Pritchard 2017).") print(" Our synthetic is a caricature of that, built so we know the answer.") print() print(" Profile, every 10 km (mm/yr):") print(" r(km) : " + " ".join(f"{r[i]/1e3:>7.0f}" for i in range(0, r.size, 5))) print(" uz : " + " ".join(f"{obs_mm[i]:>7.2f}" for i in range(0, r.size, 5))) # ====================================================================== # PART C -- fit one source, then two # ====================================================================== def fit_n_sources(depths, r_, obs_): """Least-squares volume rates for FIXED depths. Returns (dV, pred, rms).""" G = np.column_stack([unit_uz(r_, d) for d in depths]) dV, *_ = np.linalg.lstsq(G, obs_, rcond=None) pred = G @ dV res = obs_ - pred rms = float(np.sqrt(np.mean(res ** 2))) return dV, pred, rms def var_reduction(obs_, pred_): return 100.0 * (1.0 - np.sum((obs_ - pred_) ** 2) / np.sum(obs_ ** 2)) sub("PART C1 -- BEST SINGLE MOGI SOURCE") d_grid = np.arange(2e3, 200e3 + 1.0, 250.0) # 2 to 200 km, 250 m steps best1 = None curve1 = [] for d in d_grid: dV, pred, rms = fit_n_sources([d], r, obs) curve1.append(rms * 1000.0) if best1 is None or rms < best1[2]: best1 = (d, dV[0], rms, pred) curve1 = np.array(curve1) d1b, dV1b, rms1, pred1 = best1 vr1 = var_reduction(obs, pred1) pred1_mm = pred1 * 1000.0 print(f" grid: {d_grid.size} depths from {d_grid[0]/1e3:.0f} to " f"{d_grid[-1]/1e3:.0f} km, volume rate solved exactly at each") print() print(f" best-fit depth {d1b/1e3:.2f} km") print(f" best-fit volume rate {dV1b:+.4e} m3/yr " f"({'inflating' if dV1b > 0 else 'deflating'})") print(f" RMS misfit {rms1*1000:.3f} mm/yr") print(f" variance reduction {vr1:.2f} %") print(f" noise floor (the best any model could do) " f"{NOISE_SD*1000:.2f} mm/yr") print(f" misfit / noise floor {rms1/NOISE_SD:.1f}x") print() print(f" predicted uz at r = 0 {pred1_mm[0]:+.2f} mm/yr " f"(observed {obs_mm[0]:+.2f})") print(f" predicted uz minimum {pred1_mm.min():+.4f} mm/yr " f"(observed {obs_mm.min():+.2f})") print(f" sign changes in prediction: {int(np.sum(np.diff(np.sign(pred1_mm)) != 0))}") print(f" sign changes in data : " f"{int(np.sum(np.diff(np.sign(sig_mm)) != 0))}") print() print(" THE MOAT TEST. Take the annulus 40 to 100 km, where the data") print(" subside, and compare:") moat_mask = (r >= 40e3) & (r <= 100e3) print(f" observed mean uz there {obs_mm[moat_mask].mean():+.3f} mm/yr") print(f" single-source prediction {pred1_mm[moat_mask].mean():+.3f} mm/yr") print(f" mean signed error " f"{(pred1_mm[moat_mask] - obs_mm[moat_mask]).mean():+.3f} mm/yr") print(f" RMS error in the moat " f"{np.sqrt(np.mean((pred1_mm[moat_mask]-obs_mm[moat_mask])**2)):.3f} mm/yr") print() print(" Note the single source cannot get the sign right in the moat") print(" AND the sign right at the centre. It is not being out-fitted by a") print(" better model; it is barred by P1 from making the shape at all.") sub("PART C1b -- WHAT THE SINGLE-SOURCE MISFIT LOOKS LIKE ACROSS DEPTH") print(" There is no depth that rescues a single source. RMS misfit versus") print(" assumed depth, sampled every 10 km:") print(" depth(km) : " + " ".join( f"{d_grid[i]/1e3:>6.0f}" for i in range(32, d_grid.size, 40))) print(" RMS(mm/yr): " + " ".join( f"{curve1[i]:>6.3f}" for i in range(32, d_grid.size, 40))) print(f" flattest value anywhere on the grid: {curve1.min():.3f} mm/yr") print(f" worst value anywhere on the grid : {curve1.max():.3f} mm/yr") print(" Even the best depth leaves a misfit several times the noise.") sub("PART C2 -- BEST PAIR OF MOGI SOURCES") # coarse grid over (d1, d2), d1 < d2 g1 = np.arange(5e3, 60e3 + 1.0, 1e3) g2 = np.arange(30e3, 200e3 + 1.0, 2e3) best2 = None for a in g1: for b in g2: if b <= a + 5e3: continue dV, pred, rms = fit_n_sources([a, b], r, obs) if best2 is None or rms < best2[0]: best2 = (rms, a, b, dV, pred) # local refinement, 50 m steps rms2, a0, b0, _, _ = best2 for a in np.arange(a0 - 1.5e3, a0 + 1.5e3 + 1.0, 50.0): for b in np.arange(b0 - 4e3, b0 + 4e3 + 1.0, 50.0): if a <= 1e3 or b <= a + 2e3: continue dV, pred, rms = fit_n_sources([a, b], r, obs) if rms < best2[0]: best2 = (rms, a, b, dV, pred) rms2, d2a, d2b, dV2, pred2 = best2 vr2 = var_reduction(obs, pred2) pred2_mm = pred2 * 1000.0 print(f" coarse grid {g1.size} x {g2.size} depth pairs, then refined at 50 m") print() print(f" shallow source depth {d2a/1e3:>7.2f} km " f"dV = {dV2[0]:+.4e} m3/yr " f"({'inflating' if dV2[0] > 0 else 'deflating'})") print(f" deep source depth {d2b/1e3:>7.2f} km " f"dV = {dV2[1]:+.4e} m3/yr " f"({'inflating' if dV2[1] > 0 else 'deflating'})") print(f" |dV_deep| / dV_shallow {abs(dV2[1])/abs(dV2[0]):.2f}") print(f" net volume rate {dV2.sum():+.4e} m3/yr") print() print(f" RMS misfit {rms2*1000:.3f} mm/yr") print(f" variance reduction {vr2:.2f} %") print(f" misfit / noise floor {rms2/NOISE_SD:.2f}x") print() print(" RECOVERY OF THE HIDDEN TRUTH:") print(f" depth 1 true {TRUE_D1/1e3:>6.2f} km recovered {d2a/1e3:>6.2f} km " f"error {100*(d2a-TRUE_D1)/TRUE_D1:+.1f} %") print(f" depth 2 true {TRUE_D2/1e3:>6.2f} km recovered {d2b/1e3:>6.2f} km " f"error {100*(d2b-TRUE_D2)/TRUE_D2:+.1f} %") print(f" dV 1 true {TRUE_DV1:+.3e} recovered {dV2[0]:+.3e} " f"error {100*(dV2[0]-TRUE_DV1)/abs(TRUE_DV1):+.1f} %") print(f" dV 2 true {TRUE_DV2:+.3e} recovered {dV2[1]:+.3e} " f"error {100*(dV2[1]-TRUE_DV2)/abs(TRUE_DV2):+.1f} %") sub("PART C3 -- HEAD TO HEAD, AND WHY THE MOAT IS THE DIAGNOSTIC") print(f" {'model':<22s}{'params':>8s}{'RMS mm/yr':>12s}{'VR %':>10s}" f"{'sign flips':>12s}") print(" " + "-" * 64) print(f" {'one Mogi source':<22s}{2:>8d}{rms1*1000:>12.3f}{vr1:>10.2f}" f"{int(np.sum(np.diff(np.sign(pred1_mm)) != 0)):>12d}") print(f" {'two Mogi sources':<22s}{4:>8d}{rms2*1000:>12.3f}{vr2:>10.2f}" f"{int(np.sum(np.diff(np.sign(pred2_mm)) != 0)):>12d}") print(f" {'the data themselves':<22s}{'-':>8s}{'-':>12s}{'-':>10s}" f"{int(np.sum(np.diff(np.sign(sig_mm)) != 0)):>12d}") print() print(f" misfit ratio, one source / two sources : {rms1/rms2:.2f}x") print(f" residual variance removed by the 2nd source: " f"{100*(1 - (rms2/rms1)**2):.2f} %") print() print(" The sombrero is diagnostic for a reason that has nothing to do with") print(" statistics. Adding a second source does not merely fit better; a") print(" single source is mathematically incapable of producing the shape.") print(" Two numbers make that concrete:") print(f" moat depth in the data {moat_mm:+.2f} mm/yr") print(f" deepest a single source can go, given it fits the centre: " f"{pred1_mm.min():+.4f} mm/yr") print(" The second of those can never be negative while the first is. The") print(" gap is not a fitting residual. It is a structural impossibility.") sub("PART C4 -- DEPTH READING FROM THE RADIAL DISPLACEMENT") ur_shallow = mogi_ur(r, d2a, dV2[0]) ur_deep = mogi_ur(r, d2b, dV2[1]) ur_tot_mm = (ur_shallow + ur_deep) * 1000.0 uz_tot_mm = pred2_mm print(" Horizontal motion is the other half of the measurement, and it") print(" carries the depth in a way the vertical does not.") print() print(f" shallow source: ur peaks at r = d/sqrt(2) = {d2a/math.sqrt(2)/1e3:.2f} km") print(f" deep source : ur peaks at r = d/sqrt(2) = {d2b/math.sqrt(2)/1e3:.2f} km") print() i_ur = int(np.argmax(ur_tot_mm)) print(f" combined radial field peaks at r = {r[i_ur]/1e3:.0f} km, " f"{ur_tot_mm[i_ur]:+.2f} mm/yr") print() print(" Vertical and radial, every 20 km (mm/yr):") print(" r(km) : " + " ".join(f"{r[i]/1e3:>7.0f}" for i in range(0, r.size, 10))) print(" uz : " + " ".join(f"{uz_tot_mm[i]:>7.2f}" for i in range(0, r.size, 10))) print(" ur : " + " ".join(f"{ur_tot_mm[i]:>7.2f}" for i in range(0, r.size, 10))) print() print(" For a single source uz/ur = d/r exactly. Measuring both at one") print(" radius fixes the depth with no inversion at all:") for rr in (10e3, 20e3, 40e3): uz_s = float(mogi_uz(rr, d2a, dV2[0])) ur_s = float(mogi_ur(rr, d2a, dV2[0])) print(f" r = {rr/1e3:>3.0f} km : uz/ur = {uz_s/ur_s:.4f} " f"-> d = {rr*uz_s/ur_s/1e3:.2f} km") # ====================================================================== # PART D -- the eruption arithmetic # ====================================================================== sub("PART D -- WHAT THE INFLATION RATE WOULD HAVE TO BUILD") DV_INFLATE = float(dV2[0]) # our fitted shallow-source volume rate LAST_ERUPT_KA = 250.0 # Liu et al. 2025: 250 +/- 5 ka EDIFICE_KM3 = 85.0 # Sparks et al. 2008 EDIFICE_SPAN_YR = 640e3 # about 890 ka to 250 ka UNREST_YEARS = 33.0 # 1992 InSAR discovery to 2025 RHO_MAGMA = 2400.0 # kg/m3, dacite RHO_FLUID_LO, RHO_FLUID_HI = 200.0, 400.0 # kg/m3, supercritical H2O-CO2 MELT_FRAC = 0.25 # Liu et al. 2025, APMB maximum MUSH_FACTOR = 1.0 / MELT_FRAC # mush to reorganise per unit of melt q_long = EDIFICE_KM3 * 1e9 / EDIFICE_SPAN_YR print(f" fitted shallow inflation rate dV = {DV_INFLATE:.4e} m3/yr") print(f" = {DV_INFLATE/1e9:.5f} km3/yr") print(f" = {DV_INFLATE/3.15576e7:.3f} m3/s") print(f" published estimate for Uturuncu: about 1 m3/s, 1e-2 km3/yr") print(f" (Sparks et al. 2008, from 1992-2006 InSAR). Same ballpark,") print(f" which is the only claim we make for our number.") print() print(f" Uturuncu's own long-run magma output rate:") print(f" edifice {EDIFICE_KM3:.0f} km3 over {EDIFICE_SPAN_YR/1e3:.0f} kyr " f"= {q_long:.3e} m3/yr") print(f" RATIO, inflation rate / long-run magma output rate: " f"{DV_INFLATE/q_long:.0f}x") print() print(" That ratio is the first thing that should make you suspicious. If") print(" the volume change were magma arriving, Uturuncu would currently be") print(f" supplied {DV_INFLATE/q_long:.0f} times faster than it was during the " f"640,000 years") print(" when it was actually building a mountain.") sub("PART D2 -- THE REDUCTIO") cum_since_erupt = DV_INFLATE * LAST_ERUPT_KA * 1e3 cum_since_unrest = DV_INFLATE * UNREST_YEARS print(f" If the present rate had run steadily since the last eruption") print(f" ({LAST_ERUPT_KA:.0f} ka), the accumulated volume would be:") print(f" {cum_since_erupt:.4e} m3 = {cum_since_erupt/1e9:,.0f} km3") print(f" = {cum_since_erupt/1e9/EDIFICE_KM3:.0f}x the entire volcano") print(f" = {cum_since_erupt/1e9/2500:.1f}x the Atana ignimbrite (2,500 km3),") print(" the largest eruption the Altiplano-Puna has ever produced.") print() print(" Nothing like that is underground. So the rate is not steady, and") print(" whatever the ground is doing it has not been doing it for long.") print() print(f" Since the uplift was first measured ({UNREST_YEARS:.0f} years):") print(f" accumulated volume change {cum_since_unrest:.4e} m3 " f"= {cum_since_unrest/1e9:.3f} km3") print(f" that is {cum_since_unrest/1e9/0.1:.1f}x a VEI 4, " f"{cum_since_unrest/1e9/1.0:.2f}x a VEI 5, " f"{cum_since_unrest/1e9/10.0:.3f}x a VEI 6") sub("PART D2b -- THE NET BUDGET, WHICH POINTS THE WRONG WAY") net_dv = float(dV2.sum()) print(" The fit does not just find a shallow source filling up. It finds a") print(" deep source emptying out, and the deep one is the larger of the two.") print(f" shallow source {dV2[0]:+.4e} m3/yr") print(f" deep source {dV2[1]:+.4e} m3/yr") print(f" NET {net_dv:+.4e} m3/yr " f"= {net_dv/1e9:+.4f} km3/yr") print(f" Over the {UNREST_YEARS:.0f} years of measured unrest the net volume " f"change is") print(f" {net_dv*UNREST_YEARS/1e9:+.2f} km3. Negative. Taken at face value the " f"system as a whole") print(" is losing volume while its top inflates, which is exactly what you") print(" expect if something mobile is moving UP out of a deep reservoir") print(" rather than new material arriving from below.") print() print(" No eruptible body is being assembled on that budget. It is being") print(" disassembled and redistributed upward.") sub("PART D3 -- CENTURIES NEEDED, THE LADDER") targets = [ ("VEI 4, a modest dome or flow", 0.1), ("VEI 5, 1991 Pinatubo scale", 1.0), ("VEI 6, 1883 Krakatau scale", 10.0), ("Uturuncu's whole edifice", 85.0), ("Pastos Grandes ignimbrite", 500.0), ("Atana ignimbrite", 2500.0), ] print(f" {'target':<30s}{'km3':>9s}{'yr, naive':>12s}" f"{'yr, mush':>11s}{'yr, long-run':>14s}") print(" " + "-" * 76) for name, v in targets: vol = v * 1e9 naive = vol / DV_INFLATE mushy = vol * MUSH_FACTOR / DV_INFLATE longrun = vol / q_long print(f" {name:<30s}{v:>9.1f}{naive:>12,.0f}{mushy:>11,.0f}{longrun:>14,.0f}") print() print(" 'naive' : every cubic metre of volume change is eruptible magma.") print(f" 'mush' : the magma body is {MELT_FRAC*100:.0f}% melt, so " f"{MUSH_FACTOR:.0f} m3 of mush must be") print(" reorganised per m3 of eruptible melt extracted.") print(" 'long-run' : magma arrives at the rate that actually built the") print(" volcano over the last 890,000 years.") print() v5_centuries = 1e9 * MUSH_FACTOR / DV_INFLATE / 100.0 print(f" Headline: assembling one cubic kilometre of eruptible magma at the") print(f" observed inflation rate, allowing for the {MELT_FRAC*100:.0f}% melt " f"fraction the") print(f" paper images, takes {v5_centuries:.1f} centuries. At the rate that " f"actually built") print(f" the mountain it takes {1e9/q_long/100:.0f} centuries.") sub("PART D4 -- MASS, NOT VOLUME: WHY FLUID IS THE CHEAP EXPLANATION") m_magma = DV_INFLATE * RHO_MAGMA m_fl_lo = DV_INFLATE * RHO_FLUID_LO m_fl_hi = DV_INFLATE * RHO_FLUID_HI print(f" The same {DV_INFLATE:.3e} m3/yr of volume change costs:") print(f" as dacitic magma at {RHO_MAGMA:.0f} kg/m3 " f"{m_magma:.3e} kg/yr = {m_magma/1e9:.2f} Mt/yr") print(f" as supercritical fluid at {RHO_FLUID_LO:.0f} kg/m3 " f"{m_fl_lo:.3e} kg/yr = {m_fl_lo/1e9:.2f} Mt/yr") print(f" as supercritical fluid at {RHO_FLUID_HI:.0f} kg/m3 " f"{m_fl_hi:.3e} kg/yr = {m_fl_hi/1e9:.2f} Mt/yr") print(f" mass ratio magma / fluid : {m_magma/m_fl_hi:.0f}x to " f"{m_magma/m_fl_lo:.0f}x") print() print(" And that is still an overstatement of the fluid's cost, because a") print(" compressible fluid can change the volume of a pore network by") print(" expanding in place, with no new mass arriving at all. Volume change") print(" is not mass change. The Mogi model cannot tell the two apart, and") print(" that single blind spot is what forty years of deformation-only") print(" interpretation at Uturuncu kept walking into.") # ====================================================================== # PART E -- numbers the article and the interactive page quote # ====================================================================== sub("PART E -- HEADLINE NUMBERS FOR THE ARTICLE AND THE INTERACTIVE PAGE") print(f" synthetic peak uplift {peak_mm:+.2f} mm/yr") print(f" synthetic moat minimum {moat_mm:+.2f} mm/yr " f"at {moat_r:.0f} km") print(f" zero crossing {cross:.1f} km") print(f" moat/peak ratio {abs(moat_mm)/peak_mm:.3f}") print() print(f" ONE SOURCE depth {d1b/1e3:.1f} km") print(f" ONE SOURCE dV {dV1b:+.3e} m3/yr") print(f" ONE SOURCE RMS {rms1*1000:.2f} mm/yr") print(f" ONE SOURCE variance reduction {vr1:.1f} %") print(f" ONE SOURCE min predicted uz {pred1_mm.min():+.4f} mm/yr") print() print(f" TWO SOURCE shallow depth {d2a/1e3:.1f} km") print(f" TWO SOURCE shallow dV {dV2[0]:+.3e} m3/yr") print(f" TWO SOURCE deep depth {d2b/1e3:.1f} km") print(f" TWO SOURCE deep dV {dV2[1]:+.3e} m3/yr") print(f" TWO SOURCE RMS {rms2*1000:.2f} mm/yr") print(f" TWO SOURCE variance reduction {vr2:.2f} %") print(f" TWO SOURCE |deep|/shallow {abs(dV2[1])/abs(dV2[0]):.2f}") print(f" TWO SOURCE net dV {dV2.sum():+.3e} m3/yr") print(f" net over 33 yr {dV2.sum()*UNREST_YEARS/1e9:+.2f} km3") print() print(f" misfit ratio 1-source : 2-source {rms1/rms2:.1f}x") print(f" inflation / long-run magma supply {DV_INFLATE/q_long:.0f}x") print(f" volume if steady since 250 ka {cum_since_erupt/1e9:,.0f} km3 " f"= {cum_since_erupt/1e9/EDIFICE_KM3:.0f}x the edifice") print(f" 33 years of inflation {cum_since_unrest/1e9:.2f} km3") print(f" centuries for 1 km3 eruptible {v5_centuries:.1f}") print(f" magma:fluid mass ratio {m_magma/m_fl_hi:.0f}x to " f"{m_magma/m_fl_lo:.0f}x") print(SEP) # ---------------------------------------------------------------------- # Data dumps for the article's hand-built SVG figures. # ---------------------------------------------------------------------- sub("APPENDIX -- PROFILE TABLE FOR THE FIGURES (mm/yr)") print(f" {'r(km)':>7s}{'observed':>11s}{'1-source':>11s}" f"{'2-source':>11s}{'radial':>10s}") for i in range(0, r.size, 2): print(f" {r[i]/1e3:>7.0f}{obs_mm[i]:>11.3f}{pred1_mm[i]:>11.3f}" f"{pred2_mm[i]:>11.3f}{ur_tot_mm[i]:>10.3f}") sub("APPENDIX -- SINGLE-SOURCE RMS VERSUS ASSUMED DEPTH (mm/yr)") for i in range(0, d_grid.size, 20): print(f" {d_grid[i]/1e3:>7.1f} km {curve1[i]:>8.3f}") print(SEP)