L03: reconstruct earth heat budget, marginal stability, Nu-Ra scaling

Author: timlichtenberg

book/03_heat_energy/figures/earth_heat_budget.avif

@@ -103,7 +103,7 @@ The total heat flowing out of Earth's interior is approximately **47 TW** {cite:
 :width: 450px
 :align: center
 
-Earth's surface heat-flux budget (~47 TW), decomposed into radiogenic (~20 TW), primordial / secular cooling (~15 TW), core cooling and solidification (~10 TW), and lunar tidal dissipation (~0.1 TW). Tidal contribution is too small to render at this scale. After {cite:p}`DaviesDavies2010`. Custom plot.
+Earth's surface heat-flux budget (~47 TW), decomposed into radiogenic (~20 TW), primordial / secular cooling (~15 TW), core cooling and solidification (~10 TW), and lunar tidal dissipation (~0.1 TW). Tidal contribution is too small to render at this scale. After {cite:p}`DaviesDavies2010`. Generated by `scripts/figures/L03_heat_energy/fig_earth_heat_budget.py`.
 ```
 
 Roughly 40% of Earth's current heat loss is powered by ongoing radioactive decay, with the rest being the slow release of primordial heat stored since formation. This means Earth is still cooling: its interior temperature is gradually decreasing over geological time. The relative importance of these sources varies strongly across the solar system: small bodies like the Moon have lost most of their primordial heat long ago, while a tidally heated moon like Io has a heat budget dominated by tidal dissipation rather than radiogenic or primordial heat.
@@ -291,7 +291,7 @@ $$
 :width: 600px
 :align: center
 
-Marginal stability curve for Rayleigh-Bénard convection between rigid horizontal boundaries (schematic, normalised to the canonical critical values). Convection grows for parameter combinations above the curve; conduction is the only steady transport mode below it. The minimum of the curve is the critical Rayleigh number $\mathrm{Ra}_c\approx 1708$, attained at non-dimensional wavenumber $k_c\approx 3.12$ {cite:p}`Turcotte2002`. Custom plot.
+Marginal stability curve for Rayleigh-Bénard convection between rigid horizontal boundaries (schematic, normalised to the canonical critical values). Convection grows for parameter combinations above the curve; conduction is the only steady transport mode below it. The minimum of the curve is the critical Rayleigh number $\mathrm{Ra}_c\approx 1708$, attained at non-dimensional wavenumber $k_c\approx 3.12$ {cite:p}`Turcotte2002`. Generated by `scripts/figures/L03_heat_energy/fig_marginal_stability.py`.
 ```
 
 For Earth's mantle, typical parameters are $\alpha \sim 2 \times 10^{-5}$ K$^{-1}$, $\rho \sim 4000$ kg m$^{-3}$, $g \sim 10$ m s$^{-2}$, $\Delta T \sim 2500$ K, $d \sim 3 \times 10^6$ m, $\kappa \sim 10^{-6}$ m$^2$ s$^{-1}$, and $\eta \sim 10^{21}$ Pa s. This gives:
@@ -331,7 +331,7 @@ This means convection transports heat roughly 200 times more efficiently than co
 :width: 600px
 :align: center
 
-Schematic Nusselt-Rayleigh scaling. Below the critical Rayleigh number $\mathrm{Ra}_c\approx 1708$ heat is conducted only ($\mathrm{Nu}=1$); above it the boundary-layer scaling $\mathrm{Nu}\propto\mathrm{Ra}^{1/3}$ applies (Eq. {eq}`eq:nu-ra-scaling`). Earth's mantle, with $\mathrm{Ra}\sim 10^7\text{--}10^8$, transports heat $\sim 200$ to $500$ times more efficiently than conduction alone under this idealised isoviscous scaling. Custom plot.
+Schematic Nusselt-Rayleigh scaling. Below the critical Rayleigh number $\mathrm{Ra}_c\approx 1708$ heat is conducted only ($\mathrm{Nu}=1$); above it the boundary-layer scaling $\mathrm{Nu}\propto\mathrm{Ra}^{1/3}$ applies (Eq. {eq}`eq:nu-ra-scaling`). Earth's mantle, with $\mathrm{Ra}\sim 10^7\text{--}10^8$, transports heat $\sim 200$ to $500$ times more efficiently than conduction alone under this idealised isoviscous scaling. Generated by `scripts/figures/L03_heat_energy/fig_nu_ra_scaling.py`.
 ```
 
 

book/03_heat_energy/figures/marginal_stability.avif

@@ -0,0 +1,26 @@
+{
+  "purpose": "Fig. fig:earth-heat-budget: Earth surface heat budget pie",
+  "total_TW": 47.0,
+  "components": [
+    {
+      "label": "Radiogenic\n(mantle + crust)",
+      "TW": 20.0
+    },
+    {
+      "label": "Primordial\n(secular cooling)",
+      "TW": 15.0
+    },
+    {
+      "label": "Core cooling\nand solidification",
+      "TW": 10.0
+    },
+    {
+      "label": "Tidal (lunar)",
+      "TW": 0.1
+    }
+  ],
+  "source": "Davies, J. H. & Davies, D. R. (2010), Solid Earth 1, 5-24",
+  "doi": "10.5194/se-1-5-2010",
+  "citation_key": "DaviesDavies2010",
+  "license_note": "Component figures are conventional review values."
+}

book/03_heat_energy/figures/nu_ra_scaling.avif

@@ -0,0 +1,87 @@
+"""Generate Fig. (`fig:earth-heat-budget`).
+
+Pie chart of Earth's surface heat-flux budget, ~47 TW total,
+decomposed into radiogenic, primordial / secular cooling, core
+cooling and solidification, and lunar tidal dissipation.
+
+Caption / figure id : `fig:earth-heat-budget`
+Markdown source     : book/03_heat_energy/heat_energy.md
+Citation key        : DaviesDavies2010
+
+Component values are static numerics from the Davies & Davies 2010
+review and are tabulated in `data/earth_heat_budget.json`.
+"""
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+
+from scripts.figures._shared.style import apply_style, save_figure
+
+
+REPO_ROOT = Path(__file__).resolve().parents[3]
+DATA_DIR = Path(__file__).resolve().parent / "data"
+META = DATA_DIR / "earth_heat_budget.json"
+OUT_AVIF = REPO_ROOT / "book/03_heat_energy/figures/earth_heat_budget.avif"
+
+COMPONENTS = [
+    ("Radiogenic\n(mantle + crust)", 20.0, "#d62728"),
+    ("Primordial\n(secular cooling)", 15.0, "#ff7f0e"),
+    ("Core cooling\nand solidification", 10.0, "#1f77b4"),
+    ("Tidal (lunar)", 0.1, "#2ca02c"),
+]
+
+
+def write_metadata() -> None:
+    DATA_DIR.mkdir(parents=True, exist_ok=True)
+    META.write_text(json.dumps({
+        "purpose": "Fig. fig:earth-heat-budget: Earth surface heat budget pie",
+        "total_TW": 47.0,
+        "components": [{"label": l, "TW": v} for l, v, _ in COMPONENTS],
+        "source": "Davies, J. H. & Davies, D. R. (2010), Solid Earth 1, 5-24",
+        "doi": "10.5194/se-1-5-2010",
+        "citation_key": "DaviesDavies2010",
+        "license_note": "Component figures are conventional review values.",
+    }, indent=2))
+
+
+def make_plot() -> Path:
+    apply_style()
+    write_metadata()
+    labels, sizes, colors = zip(*[(l, v, c) for l, v, c in COMPONENTS])
+
+    fig, ax = plt.subplots(figsize=(7, 5.5))
+    wedges, texts = ax.pie(
+        sizes, colors=colors,
+        labels=[f"{l}\n~{v:.1f} TW" for l, v in zip(labels, sizes)],
+        startangle=90, counterclock=False,
+        wedgeprops={"edgecolor": "white", "linewidth": 1.5},
+        textprops={"fontsize": 10},
+        labeldistance=1.1,
+    )
+
+    # Inside-wedge values for the larger components
+    for w, val in zip(wedges, sizes):
+        if val >= 5:
+            ang = (w.theta1 + w.theta2) / 2
+            import numpy as np
+            x = 0.6 * np.cos(np.radians(ang))
+            y = 0.6 * np.sin(np.radians(ang))
+            ax.text(x, y, f"{val:.0f} TW", ha="center", va="center",
+                    fontsize=11, color="white", fontweight="bold")
+
+    ax.set_title("Earth's surface heat-flux budget (~47 TW)", fontsize=12, pad=18)
+    fig.tight_layout()
+    return save_figure(fig, OUT_AVIF, avif_quality=80)
+
+
+def main() -> None:
+    out = make_plot()
+    print(f"  meta : {META}")
+    print(f"  plot : {out}")
+
+
+if __name__ == "__main__":
+    main()

book/03_heat_energy/heat_energy.md

@@ -0,0 +1,91 @@
+"""Generate Fig. (`fig:marginal-stability`).
+
+Schematic Rayleigh-Bénard marginal-stability curve normalised to the
+canonical critical values: Ra_c ≈ 1708 at non-dimensional wavenumber
+k_c ≈ 3.117 between rigid horizontal boundaries.
+
+Caption / figure id : `fig:marginal-stability`
+Markdown source     : book/03_heat_energy/heat_energy.md
+Citation key        : Chandrasekhar1961
+
+The exact stability curve solves a transcendental boundary-value
+problem; here we plot the standard analytic approximation:
+
+    Ra(k) = (k^2 + pi^2)^3 / k^2
+
+This is the *free-slip* (stress-free) result with Ra_c = 657.5 at
+k_c = pi/sqrt(2) ≈ 2.221. The classical *rigid-rigid* boundary case
+(Ra_c = 1707.762, k_c = 3.117) does not have a closed-form curve
+but is well approximated for visualisation by scaling the free-slip
+curve so its minimum sits at the rigid values, which is what we do.
+"""
+from __future__ import annotations
+
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from scripts.figures._shared.style import apply_style, save_figure
+
+
+REPO_ROOT = Path(__file__).resolve().parents[3]
+OUT_AVIF = REPO_ROOT / "book/03_heat_energy/figures/marginal_stability.avif"
+
+# Rigid-rigid critical values (Chandrasekhar 1961, Table II)
+RA_C = 1707.762
+K_C = 3.117
+
+
+def make_plot() -> Path:
+    apply_style()
+
+    # Free-slip stability curve (analytic) shifted to land on the
+    # rigid-rigid critical values for visualisation.
+    k = np.linspace(0.4, 8.0, 400)
+    pi = np.pi
+    free_slip = (k ** 2 + pi ** 2) ** 3 / k ** 2
+    free_slip_min = (4.0 / 27.0) * (3.0 * pi ** 2) ** 3 / (3 * pi ** 2) ** 2  # value at k = pi/sqrt(2)
+    # Scale and shift in k so the minimum hits (K_C, RA_C)
+    scale_k = K_C / (pi / np.sqrt(2))
+    k_shifted = k * scale_k
+    Ra_shifted = free_slip * (RA_C / free_slip_min)
+
+    fig, ax = plt.subplots(figsize=(6.5, 4.4))
+    ax.plot(k_shifted, Ra_shifted, color="#1f77b4", lw=2.0)
+
+    # Fill convection-on region
+    ax.fill_between(k_shifted, Ra_shifted, 1e6, color="#ff7f0e", alpha=0.13,
+                    label="Convection")
+    ax.fill_between(k_shifted, 0, Ra_shifted, color="#1f77b4", alpha=0.10,
+                    label="Conduction only")
+
+    # Mark critical point
+    ax.plot(K_C, RA_C, "o", color="#d62728", markeredgecolor="black",
+            markersize=8, zorder=5)
+    ax.annotate(f"$\\mathrm{{Ra}}_c \\approx {RA_C:.0f}$,\n"
+                f"$k_c \\approx {K_C:.2f}$",
+                xy=(K_C, RA_C), xytext=(K_C + 0.5, RA_C * 0.45),
+                fontsize=10,
+                arrowprops=dict(arrowstyle="->", lw=0.8, color="dimgray"))
+
+    ax.set_yscale("log")
+    ax.set_xlabel(r"Non-dimensional horizontal wavenumber $k$")
+    ax.set_ylabel(r"Rayleigh number $\mathrm{Ra}$")
+    ax.set_xlim(0, 8)
+    ax.set_ylim(5e2, 5e5)
+    ax.set_title("Marginal stability curve (rigid horizontal boundaries)")
+    ax.grid(which="both", linestyle=":", alpha=0.3)
+    ax.legend(loc="upper right", frameon=False)
+
+    fig.tight_layout()
+    return save_figure(fig, OUT_AVIF, avif_quality=80)
+
+
+def main() -> None:
+    out = make_plot()
+    print(f"  plot : {out}")
+
+
+if __name__ == "__main__":
+    main()

scripts/figures/L03_heat_energy/data/earth_heat_budget.json

@@ -0,0 +1,77 @@
+"""Generate Fig. (`fig:nu-ra-scaling`).
+
+Schematic Nusselt-Rayleigh scaling: Nu = 1 below Ra_c (pure
+conduction), Nu ~ Ra^{1/3} above (boundary-layer scaling, after the
+asymptotic boundary-layer theory of Howard 1966 and Malkus 1954).
+The shaded band marks Earth's mantle (Ra ~ 1e7 to 1e8).
+
+Caption / figure id : `fig:nu-ra-scaling`
+Markdown source     : book/03_heat_energy/heat_energy.md
+Citation key        : (textbook scaling; cited as Turcotte2002 in body)
+"""
+from __future__ import annotations
+
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from scripts.figures._shared.style import apply_style, save_figure
+
+
+REPO_ROOT = Path(__file__).resolve().parents[3]
+OUT_AVIF = REPO_ROOT / "book/03_heat_energy/figures/nu_ra_scaling.avif"
+
+RA_C = 1708.0
+
+
+def make_plot() -> Path:
+    apply_style()
+    Ra = np.logspace(2, 11, 400)
+
+    # Conduction below Ra_c, Nu ~ Ra^(1/3) above (with the canonical
+    # prefactor that gives Nu = 1 at Ra = Ra_c, i.e. continuous at the
+    # bifurcation).
+    Nu = np.ones_like(Ra)
+    above = Ra > RA_C
+    Nu[above] = (Ra[above] / RA_C) ** (1.0 / 3.0)
+
+    fig, ax = plt.subplots(figsize=(7, 4.4))
+    ax.plot(Ra, Nu, color="#1f77b4", lw=2.0,
+            label=r"$\mathrm{Nu}(\mathrm{Ra})$")
+
+    ax.axvline(RA_C, color="#d62728", lw=1.0, linestyle="--",
+               label=fr"$\mathrm{{Ra}}_c = {RA_C:.0f}$")
+
+    # Earth-mantle band
+    ax.axvspan(1e7, 1e8, color="#ff7f0e", alpha=0.18,
+               label=r"Earth's mantle (Ra $\sim 10^7$-$10^8$)")
+
+    # Asymptotic slope label
+    ax.text(1e9, (1e9 / RA_C) ** (1.0 / 3.0) * 1.4,
+            r"$\mathrm{Nu} \propto \mathrm{Ra}^{1/3}$",
+            fontsize=11, color="#1f77b4")
+    ax.text(20, 0.7, r"$\mathrm{Nu} = 1$ (conduction)",
+            fontsize=10, color="#1f77b4")
+
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+    ax.set_xlim(1e2, 1e11)
+    ax.set_ylim(0.5, 1e3)
+    ax.set_xlabel("Rayleigh number $\\mathrm{Ra}$")
+    ax.set_ylabel("Nusselt number $\\mathrm{Nu}$")
+    ax.set_title("Heat-transport efficiency: $\\mathrm{Nu}$ vs $\\mathrm{Ra}$")
+    ax.grid(which="both", linestyle=":", alpha=0.3)
+    ax.legend(loc="lower right", frameon=False)
+
+    fig.tight_layout()
+    return save_figure(fig, OUT_AVIF, avif_quality=80)
+
+
+def main() -> None:
+    out = make_plot()
+    print(f"  plot : {out}")
+
+
+if __name__ == "__main__":
+    main()