#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
upside-down-rocks.py
Science Journaling Club — Field Notes · Explainer · Geology

A DELIBERATELY SIMPLE physical model of the "sinkite" density inversion reported by
Rudjord & Huuse (2025), Communications Earth & Environment 6, 490,
doi:10.1038/s43247-025-02398-8.

IMPORTANT HONESTY NOTE
----------------------
This is the club's own back-of-envelope model, NOT the analysis in the paper. The
paper is an interpretation of basin-scale 3D seismic reflection data plus wireline
logs and cuttings mineralogy from hundreds of wells. It does not publish the
viscosity of Miocene ooze, and neither does anyone else, because nobody has ever
measured it. Everything below is a sanity check on whether the proposed mechanism
is even dimensionally allowed, using textbook relations an advanced high-school
student can follow. Where we had to choose a number, we say so.

WHAT THE SCRIPT COMPUTES
------------------------
Part 1 — Bulk density of two water-saturated sediments.
    For a saturated porous medium with porosity phi and grain (matrix) density
    rho_g in pore fluid of density rho_w:
            rho_bulk = phi * rho_w + (1 - phi) * rho_g
    Sand: quartz grains, rho_g = 2650 kg/m^3 (Hamilton 1976; standard value).
    Ooze: opal-A (hydrous biogenic silica from diatoms/radiolaria),
          rho_g = 2100 kg/m^3 — opal-A is markedly less dense than quartz
          because it is a hydrated, amorphous silica with structural water.
          We test 1950-2200 as a sensitivity band.
    Pore fluid: North Sea formation brine, rho_w = 1030 kg/m^3.
    We sweep realistic porosities: sand 0.30-0.50, ooze 0.45-0.75, and locate the
    window in which the YOUNGER SAND IS DENSER THAN THE OLDER OOZE — the
    gravitationally unstable arrangement the paper requires.

Part 2 — Buoyancy force per unit volume.
    f = (rho_sand - rho_ooze) * g   [N/m^3], g = 9.81 m/s^2.
    This is the driving body force of the instability. Also expressed as the
    pressure difference it generates across a 200 m thick ooze layer.

Part 3 — Sink velocity, Stokes-type.
    A sphere of radius R sinking at low Reynolds number through a fluid of
    dynamic viscosity mu, under density contrast d_rho, reaches terminal velocity
            v = 2 * d_rho * g * R^2 / (9 * mu)                 (Stokes 1851)
    We treat a km-scale sand body as a sphere of radius R = 500 m. This is a
    generous idealisation: the real bodies are irregular, the ooze is not a
    Newtonian fluid, and the sand is descending through fractures rather than
    ploughing through intact rock. Stokes gives an UPPER BOUND on speed for a
    given viscosity, so the times below are LOWER BOUNDS.
    We sweep mu over 10^6 to 10^20 Pa s (water is 1e-3; window glass at room
    temperature is ~1e18; rock salt ~1e17-1e18; the mantle ~1e21) because the
    viscosity of liquefied ooze is the single biggest unknown in the problem.
    We report time to sink 200 m (the observed vertical scale of the mounds)
    and check the Reynolds number so we know when Stokes is legitimate.

Part 4 — Liquefaction threshold.
    Terzaghi effective stress: sigma' = sigma_v - u.
    Mohr-Coulomb strength of a cohesionless sand: tau = sigma' * tan(phi_fric).
    Excess pore pressure ratio r_u = delta_u / sigma'_v0, where
            sigma'_v0 = (rho_bulk - rho_w) * g * z      (z = depth below seafloor)
    At r_u = 1 the effective stress is zero and the sand has NO shear strength:
    it is a liquid. We compute the absolute pore pressure required at several
    burial depths, and the residual strength as a function of r_u.

Part 5 — Rayleigh-Taylor dominant wavelength.
    Classical two-layer viscous RT result (Turcotte & Schubert, Geodynamics):
    for a light layer of thickness b beneath a dense layer of comparable
    thickness and viscosity, the fastest-growing wavelength is
            lambda_max = 2.568 * b
    With b = the 200-500 m ooze interval this predicts a natural spacing that we
    compare against the observed ~1 km mound scale. This is a consistency check,
    not a proof.

ASSUMPTIONS, STATED PLAINLY
---------------------------
 * Newtonian, isoviscous ooze. Almost certainly wrong; real sediments are
   visco-plastic and strain-rate dependent.
 * Single sphere, unbounded fluid. Ignores neighbouring bodies, the free
   surface, the finite ooze thickness, and the polygonal fault network that the
   paper says actually guides the sand down.
 * Densities held constant with depth; no compaction during descent.
 * No opal-A -> opal-CT diagenesis, which in reality stiffens and densifies the
   ooze and would shut the process off once it happens.
 * Liquefaction treated as an instantaneous switch rather than a cyclic process.
 * g = 9.81 m/s^2, geometry Euclidean, spherical cow duly acknowledged.

Deterministic: no RNG is used. A seed is set anyway so that any future
stochastic extension is reproducible.

Run:  python upside-down-rocks.py
"""

import math
import random
import sys

try:
    sys.stdout.reconfigure(encoding="utf-8")
except Exception:  # pragma: no cover
    pass

try:
    import numpy as np
    HAVE_NUMPY = True
except ImportError:  # pragma: no cover
    HAVE_NUMPY = False

random.seed(20250621)  # date of publication of the anchor paper; RNG unused
if HAVE_NUMPY:
    np.random.seed(20250621)

# ----------------------------------------------------------------------------
# Constants and chosen inputs
# ----------------------------------------------------------------------------
G = 9.81              # m/s^2
RHO_W = 1030.0        # kg/m^3, North Sea formation brine
RHO_QTZ = 2650.0      # kg/m^3, quartz grain density
RHO_OPAL = 2100.0     # kg/m^3, opal-A grain density (central estimate)
OPAL_LOW, OPAL_HIGH = 1950.0, 2200.0

R_BODY = 500.0        # m, radius of an idealised km-wide sand body
SINK_DISTANCE = 200.0 # m, observed vertical relief of the mounds
OOZE_THICK_MIN, OOZE_THICK_MAX = 200.0, 500.0   # m, host ooze interval
PHI_FRIC_DEG = 32.0   # degrees, friction angle of loose marine sand

SEC_PER_YEAR = 365.25 * 24 * 3600.0

RULE = "-" * 72


def bulk_density(phi, rho_grain, rho_fluid=RHO_W):
    """Saturated bulk density of a porous medium."""
    return phi * rho_fluid + (1.0 - phi) * rho_grain


def frange(a, b, step):
    """Inclusive float range, stdlib fallback."""
    out, x, n = [], a, 0
    while x <= b + 1e-9:
        out.append(round(a + n * step, 10))
        n += 1
        x = a + n * step
    return out


def header(title):
    print()
    print("=" * 72)
    print(title)
    print("=" * 72)


# ============================================================================
# PART 0 — inputs echo
# ============================================================================
header("UPSIDE-DOWN ROCKS — club model of buoyancy-driven stratigraphic inversion")
print("Anchor paper : Rudjord, J. E. & Huuse, M. (2025), Commun. Earth Environ. 6, 490")
print("               doi:10.1038/s43247-025-02398-8")
print("This script  : the Science Journaling Club's own simplified model.")
print("               Not the paper's analysis. Assumptions listed in docstring.")
print()
print("numpy available            : %s" % HAVE_NUMPY)
print("g                          : %.2f m/s^2" % G)
print("pore brine density         : %.0f kg/m^3" % RHO_W)
print("quartz grain density       : %.0f kg/m^3" % RHO_QTZ)
print("opal-A grain density       : %.0f kg/m^3 (range %.0f-%.0f)"
      % (RHO_OPAL, OPAL_LOW, OPAL_HIGH))
print("idealised body radius      : %.0f m" % R_BODY)
print("sink distance modelled     : %.0f m" % SINK_DISTANCE)
print("ooze interval thickness    : %.0f-%.0f m" % (OOZE_THICK_MIN, OOZE_THICK_MAX))
print("sand friction angle        : %.0f deg" % PHI_FRIC_DEG)

# ============================================================================
# PART 1 — Bulk densities across realistic porosity
# ============================================================================
header("PART 1 — Saturated bulk density vs porosity")
print("rho_bulk = phi*rho_fluid + (1-phi)*rho_grain")
print()
print("%-10s | %-22s | %-22s" % ("porosity", "SAND (quartz grains)", "OOZE (opal-A grains)"))
print("%-10s | %-22s | %-22s" % ("phi", "rho_bulk [kg/m^3]", "rho_bulk [kg/m^3]"))
print(RULE)

phis = frange(0.30, 0.80, 0.05)
sand_curve, ooze_curve = [], []
for phi in phis:
    rs = bulk_density(phi, RHO_QTZ)
    ro = bulk_density(phi, RHO_OPAL)
    sand_curve.append((phi, rs))
    ooze_curve.append((phi, ro))
    print("%-10.2f | %-22.1f | %-22.1f" % (phi, rs, ro))

print()
print("Realistic field values (what we adopt for the rest of the run):")
PHI_SAND = 0.40      # loose, shallow-buried marine sand
PHI_OOZE = 0.65      # biosiliceous ooze retains huge porosity in a rigid frame
RHO_SAND = bulk_density(PHI_SAND, RHO_QTZ)
RHO_OOZE = bulk_density(PHI_OOZE, RHO_OPAL)
print("  sand: phi = %.2f  ->  rho_sand = %.1f kg/m^3" % (PHI_SAND, RHO_SAND))
print("  ooze: phi = %.2f  ->  rho_ooze = %.1f kg/m^3" % (PHI_OOZE, RHO_OOZE))
print("  density contrast  d_rho = %.1f kg/m^3  (younger sand is the HEAVY one)"
      % (RHO_SAND - RHO_OOZE))
print("  ratio rho_sand/rho_ooze = %.3f" % (RHO_SAND / RHO_OOZE))

# ============================================================================
# PART 1b — The inversion window
# ============================================================================
header("PART 1b — Where is the density inversion window?")
print("For each ooze porosity, find the MAXIMUM sand porosity at which the sand")
print("is still denser than the ooze. Above that, the stack is stable and boring.")
print()
print("Crossover: phi_sand_max = 1 - (rho_ooze - rho_w) / (rho_qtz - rho_w)")
print()
print("%-12s | %-16s | %-16s | %-14s" %
      ("ooze phi", "rho_ooze", "max sand phi", "unstable?"))
print(RULE)

crossover = []
for phi_o in frange(0.45, 0.75, 0.05):
    ro = bulk_density(phi_o, RHO_OPAL)
    phi_s_max = 1.0 - (ro - RHO_W) / (RHO_QTZ - RHO_W)
    crossover.append((phi_o, ro, phi_s_max))
    flag = "YES" if phi_s_max > PHI_SAND else "no"
    print("%-12.2f | %-16.1f | %-16.3f | %-14s" % (phi_o, ro, phi_s_max, flag))

print()
print("Read that table as: at ooze porosity 0.65 the sand can be as loose as")
print("phi = %.3f and still outweigh the ooze it sits on. Real shallow-buried"
      % crossover[4][2])
print("marine sands are ~0.35-0.45. The window is wide open.")

# sensitivity to opal grain density
print()
print("Sensitivity of d_rho to the opal-A grain density we assumed")
print("(sand held at phi=%.2f, ooze at phi=%.2f):" % (PHI_SAND, PHI_OOZE))
print("%-22s | %-18s | %-18s" % ("opal grain density", "rho_ooze", "d_rho"))
print(RULE)
for rg in (OPAL_LOW, 2025.0, RHO_OPAL, 2150.0, OPAL_HIGH):
    ro = bulk_density(PHI_OOZE, rg)
    print("%-22.0f | %-18.1f | %-18.1f" % (rg, ro, RHO_SAND - ro))
print("Inversion survives the whole plausible opal-A range. It is not a knife-edge.")

# ============================================================================
# PART 2 — Buoyancy force per unit volume
# ============================================================================
header("PART 2 — Driving force")
D_RHO = RHO_SAND - RHO_OOZE
f_body = D_RHO * G
print("Buoyancy force per unit volume  f = d_rho * g")
print("  d_rho = %.1f kg/m^3" % D_RHO)
print("  f     = %.1f N/m^3  (= %.4f MPa per km of height)" % (f_body, f_body * 1000 / 1e6))
print()
for h in (50.0, 100.0, 200.0, 500.0):
    dp = f_body * h
    print("  driving pressure across %5.0f m of ooze : %8.0f Pa  = %.3f MPa"
          % (h, dp, dp / 1e6))
print()
print("For scale: %.3f MPa across the %.0f m ooze layer is about %.1f atmospheres"
      % (f_body * SINK_DISTANCE / 1e6, SINK_DISTANCE,
         f_body * SINK_DISTANCE / 101325.0))
print("of persistent, one-way push. Small. But geology has all the time it needs.")

# ============================================================================
# PART 3 — Stokes-type sink velocity and timescale
# ============================================================================
header("PART 3 — Sink velocity and time, Stokes-type drag")
print("v = 2 * d_rho * g * R^2 / (9 * mu)      R = %.0f m, d_rho = %.1f kg/m^3"
      % (R_BODY, D_RHO))
print("t = %.0f m / v ; Re = rho_ooze * v * (2R) / mu  (Stokes valid for Re << 1)"
      % SINK_DISTANCE)
print()
print("%-12s | %-14s | %-14s | %-16s | %-12s" %
      ("mu [Pa s]", "v [m/s]", "v [m/yr]", "t to sink 200 m", "Reynolds"))
print(RULE)

rows = []
for e in range(6, 21):
    mu = 10.0 ** e
    v = 2.0 * D_RHO * G * R_BODY ** 2 / (9.0 * mu)
    v_yr = v * SEC_PER_YEAR
    t_s = SINK_DISTANCE / v
    t_yr = t_s / SEC_PER_YEAR
    re = RHO_OOZE * v * (2 * R_BODY) / mu
    if t_yr < 1.0:
        tstr = "%.3g yr" % t_yr
    elif t_yr < 1e3:
        tstr = "%.1f yr" % t_yr
    elif t_yr < 1e6:
        tstr = "%.3g kyr" % (t_yr / 1e3)
    else:
        tstr = "%.3g Myr" % (t_yr / 1e6)
    rows.append((e, mu, v, v_yr, t_yr, re, tstr))
    print("1e%-10d | %-14.3e | %-14.3e | %-16s | %-12.2e"
          % (e, v, v_yr, tstr, re))

print()
print("Plausibility gate: the mounds formed inside the Miocene-Pliocene window,")
print("so the sinking must complete in roughly 1e3 to 1e7 years to be credible.")
print("Faster than ~1e3 yr and we would expect to see the process still running")
print("at the modern seabed; slower than ~1e7 yr and it never finishes.")
print()
print("%-12s | %-16s | %-30s" % ("mu [Pa s]", "t to sink 200 m", "verdict"))
print(RULE)
for e, mu, v, v_yr, t_yr, re, tstr in rows:
    if t_yr < 1e3:
        verdict = "too fast - would still be visible"
    elif t_yr <= 1e7:
        verdict = "GEOLOGICALLY PLAUSIBLE"
    else:
        verdict = "too slow - exceeds available time"
    print("1e%-10d | %-16s | %-30s" % (e, tstr, verdict))

plausible = [e for (e, mu, v, v_yr, t_yr, re, tstr) in rows if 1e3 <= t_yr <= 1e7]
print()
print("Plausible viscosity band: 1e%d to 1e%d Pa s (%d decades wide)."
      % (min(plausible), max(plausible), max(plausible) - min(plausible) + 1))
print("Context: water 1e-3, honey ~1e1, pitch ~1e8, window glass ~1e18,")
print("rock salt ~1e17-1e18, Earth's upper mantle ~1e21 Pa s.")
print("So the model needs the liquefied ooze/sand system to behave, on average,")
print("somewhere between hot pitch and rock salt. That is not an absurd ask,")
print("but it is an unmeasured one, and it is the weakest joint in the argument.")
print()
mu_check = 1e18
v_check = 2.0 * D_RHO * G * R_BODY ** 2 / (9.0 * mu_check)
print("Reynolds check at mu = 1e18 Pa s: Re = %.2e  -> deep creeping flow, Stokes OK."
      % (RHO_OOZE * v_check * 2 * R_BODY / mu_check))
print("Re first drops below 1 at mu = 1e8 Pa s and is astronomically small")
print("throughout the plausible band, so the drag law is at least self-consistent")
print("even where the geology is not. The two fastest rows (1e6, 1e7 Pa s) have")
print("Re > 1 and Stokes does not strictly apply there - they are rejected anyway.")

# body-size sensitivity: v scales as R^2
print()
print("Because v ~ R^2, body size matters enormously. At mu = 1e18 Pa s:")
print("%-14s | %-16s | %-16s" % ("radius [m]", "v [m/yr]", "t to sink 200 m"))
print(RULE)
for R in (50.0, 100.0, 250.0, 500.0, 1000.0):
    v = 2.0 * D_RHO * G * R ** 2 / (9.0 * mu_check)
    t_yr = SINK_DISTANCE / v / SEC_PER_YEAR
    print("%-14.0f | %-16.3e | %-13.4g kyr" % (R, v * SEC_PER_YEAR, t_yr / 1e3))
print("A 50 m blob is 100x slower than a 500 m one. Being huge is the whole")
print("reason these bodies could sink at all within the available time. Small")
print("load casts need much weaker sediment to do the same trick.")

# ============================================================================
# PART 4 — Liquefaction threshold
# ============================================================================
header("PART 4 — Liquefaction: the pore-pressure ratio that kills shear strength")
print("Terzaghi:      sigma' = sigma_v - u")
print("Mohr-Coulomb:  tau = sigma' * tan(phi_fric),  phi_fric = %.0f deg, no cohesion"
      % PHI_FRIC_DEG)
print("Excess pore pressure ratio  r_u = delta_u / sigma'_v0")
print("  sigma'_v0 = (rho_sand - rho_w) * g * z   (buoyant weight of overburden)")
print("  tau(r_u)  = (1 - r_u) * sigma'_v0 * tan(phi_fric)")
print("  => tau = 0 exactly at r_u = 1.00. That is the liquefaction threshold.")
print()
tanphi = math.tan(math.radians(PHI_FRIC_DEG))
print("%-10s | %-14s | %-14s | %-14s | %-14s" %
      ("depth z [m]", "sigma_v [MPa]", "u_hydro [MPa]", "sigma'_v0 [MPa]", "tau0 [MPa]"))
print(RULE)
depth_rows = []
for z in (10.0, 25.0, 50.0, 100.0, 200.0, 400.0):
    sig_v = RHO_SAND * G * z
    u_h = RHO_W * G * z
    sig_eff = sig_v - u_h
    tau0 = sig_eff * tanphi
    depth_rows.append((z, sig_v, u_h, sig_eff, tau0))
    print("%-10.0f | %-14.3f | %-14.3f | %-14.3f | %-14.3f"
          % (z, sig_v / 1e6, u_h / 1e6, sig_eff / 1e6, tau0 / 1e6))

print()
print("Extra pore pressure required to reach r_u = 1 (i.e. total liquefaction):")
print("%-10s | %-20s | %-22s | %-18s" %
      ("depth z [m]", "delta_u needed [MPa]", "total u at failure [MPa]", "u / sigma_v"))
print(RULE)
for z, sig_v, u_h, sig_eff, tau0 in depth_rows:
    print("%-10.0f | %-20.3f | %-22.3f | %-18.3f"
          % (z, sig_eff / 1e6, sig_v / 1e6, 1.000))

print()
print("Residual shear strength as pore pressure climbs (z = 50 m):")
print("%-10s | %-18s | %-22s" % ("r_u", "tau [kPa]", "% of drained strength"))
print(RULE)
z_ref = 50.0
sig_eff_ref = (RHO_SAND - RHO_W) * G * z_ref
tau_ref = sig_eff_ref * tanphi
for r_u in (0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 0.9, 0.95, 0.99, 1.0):
    tau = (1.0 - r_u) * sig_eff_ref * tanphi
    print("%-10.2f | %-18.2f | %-22.1f" % (r_u, tau / 1e3, 100.0 * tau / tau_ref))
print()
print("The headline: shear strength is LINEAR in r_u. Get the pore pressure to")
print("90% of the buoyant overburden stress and the sand keeps 10% of its")
print("strength; get to 100% and it keeps none. There is no cliff edge, no")
print("critical exotic condition - just a straight line running down to zero.")
print("Earthquake shaking does this routinely in the top tens of metres of")
print("loose saturated sand. The paper invokes exactly that trigger.")

# ============================================================================
# PART 5 — Rayleigh-Taylor dominant wavelength
# ============================================================================
header("PART 5 — Does the predicted spacing match kilometre-scale mounds?")
print("Two-layer viscous Rayleigh-Taylor, equal thickness b and equal viscosity:")
print("  lambda_max = 2.568 * b     (classical result, Turcotte & Schubert)")
print("Take b = the ooze interval thickness reported for the study area.")
print()
print("%-22s | %-22s | %-20s" % ("ooze thickness b [m]", "lambda_max [m]", "lambda_max [km]"))
print(RULE)
rt = []
for b in (150.0, 200.0, 300.0, 400.0, 500.0, 600.0):
    lam = 2.568 * b
    rt.append((b, lam))
    print("%-22.0f | %-22.0f | %-20.2f" % (b, lam, lam / 1000.0))
lam_lo = 2.568 * OOZE_THICK_MIN
lam_hi = 2.568 * OOZE_THICK_MAX
print()
print("For the reported %.0f-%.0f m ooze interval the model predicts a natural"
      % (OOZE_THICK_MIN, OOZE_THICK_MAX))
print("spacing of %.2f-%.2f km." % (lam_lo / 1000.0, lam_hi / 1000.0))
print("Observed mounds are described as 'several kilometres wide', with a")
print("characteristic scale around 1 km. Our number lands in the same decade.")
print("That is agreement at the level of 'not obviously wrong', which is all a")
print("one-line textbook formula is entitled to claim. It is NOT confirmation.")

# ============================================================================
# SUMMARY LEDGER
# ============================================================================
header("SUMMARY LEDGER — numbers quoted in the article")
print("%-42s %s" % ("sand bulk density (phi=0.40)", "%.0f kg/m^3" % RHO_SAND))
print("%-42s %s" % ("ooze bulk density (phi=0.65)", "%.0f kg/m^3" % RHO_OOZE))
print("%-42s %s" % ("density contrast d_rho", "%.0f kg/m^3" % D_RHO))
print("%-42s %s" % ("density ratio sand/ooze", "%.3f" % (RHO_SAND / RHO_OOZE)))
print("%-42s %s" % ("max sand porosity for inversion @ ooze 0.65",
                    "%.3f" % crossover[4][2]))
print("%-42s %s" % ("buoyancy force per unit volume", "%.0f N/m^3" % f_body))
print("%-42s %s" % ("driving pressure across 200 m ooze",
                    "%.3f MPa" % (f_body * SINK_DISTANCE / 1e6)))
print("%-42s %s" % ("plausible viscosity band",
                    "1e%d - 1e%d Pa s" % (min(plausible), max(plausible))))
v_ref = 2.0 * D_RHO * G * R_BODY ** 2 / (9.0 * 1e18)
print("%-42s %s" % ("sink rate at mu = 1e18 Pa s",
                    "%.2f cm/yr (= %.1f m/kyr)"
                    % (v_ref * SEC_PER_YEAR * 100, v_ref * SEC_PER_YEAR * 1000)))
print("%-42s %s" % ("time to sink 200 m at mu = 1e18 Pa s",
                    "%.1f kyr" % (SINK_DISTANCE / v_ref / SEC_PER_YEAR / 1e3)))
print("%-42s %s" % ("time to sink 200 m at mu = 1e20 Pa s",
                    "%.2f Myr" % (SINK_DISTANCE / (v_ref / 100.0) / SEC_PER_YEAR / 1e6)))
print("%-42s %s" % ("liquefaction threshold", "r_u = 1.00 (sigma' = 0)"))
print("%-42s %s" % ("delta_u needed at z = 50 m",
                    "%.3f MPa" % (sig_eff_ref / 1e6)))
print("%-42s %s" % ("drained shear strength at z = 50 m",
                    "%.1f kPa" % (tau_ref / 1e3)))
print("%-42s %s" % ("RT wavelength for b = 200-500 m",
                    "%.2f - %.2f km" % (lam_lo / 1000, lam_hi / 1000)))
print()
print("All figures in the article are drawn from exactly these numbers.")
print("End of run.")