Part II: Geophysics | Chapter 3

Earth's Interior

PREM, seismic discontinuities, and the thermal structure of our planet

2.1 The Preliminary Reference Earth Model (PREM)

The Preliminary Reference Earth Model (PREM) by Dziewonski and Anderson (1981) is the standard 1D radial model of Earth's interior. It specifies density $\rho(r)$, P-wave velocity $V_P(r)$, S-wave velocity $V_S(r)$, and attenuation as functions of radius, constrained by body-wave travel times, surface wave dispersion, and normal mode frequencies.

Derivation 1: The Adams-Williamson Equation

The density distribution within Earth can be estimated from seismic velocities using the Adams-Williamson equation. Start with hydrostatic equilibrium:

$$\frac{dP}{dr} = -\rho(r) g(r)$$

The adiabatic bulk modulus relates pressure changes to density changes:

$$K_S = \rho \frac{dP}{d\rho}\bigg|_S = \rho \phi$$

where the seismic parameter $\phi$ is directly measurable:

$$\phi = V_P^2 - \frac{4}{3}V_S^2 = \frac{K_S}{\rho}$$

Combining hydrostatic equilibrium with the adiabatic equation of state:

$$\frac{d\rho}{dr} = -\frac{\rho(r) g(r)}{\phi(r)} = -\frac{\rho(r) g(r)}{V_P^2 - \frac{4}{3}V_S^2}$$

This is the Adams-Williamson equation. Starting from the surface with known density and integrating inward, it yields the density profile. Deviations from the Adams-Williamson prediction indicate compositional changes or phase transitions.

Historical Context

The determination of Earth's internal structure is one of the great achievements of geophysics. Key milestones include: Wiechert (1897) proposed an iron core; Oldham (1906) identified the core seismologically; Gutenberg (1914) placed the core-mantle boundary at 2900 km; Lehmann (1936) discovered the inner core; Birch (1952) used the Adams-Williamson equation to constrain mantle composition; and Dziewonski and Anderson (1981) published PREM.

2.2 Seismic Discontinuities

Earth's interior contains several major discontinuities where seismic velocities change abruptly, indicating changes in composition or mineral phase:

Moho (Mohorovicic, 1909): Crust-mantle boundary at ~7 km (oceanic) to ~35-70 km (continental). $V_P$ jumps from ~6.5 to ~8.1 km/s.
410 km discontinuity: Olivine to wadsleyite phase transition. $V_P$ increases by ~5%.
660 km discontinuity: Ringwoodite to bridgmanite + ferropericlase. Marks the base of the transition zone.
D'' layer (~2700 km): A ~200 km thick heterogeneous region above the CMB with possible post-perovskite phase.
Core-Mantle Boundary (2891 km): $V_P$ drops from ~13.7 to ~8.1 km/s; $V_S$ drops to zero.
Inner Core Boundary (5150 km): $V_S$ reappears as liquid iron solidifies.

Derivation 2: Reflection Coefficient at a Discontinuity

When a seismic wave encounters a sharp velocity contrast, energy is reflected and transmitted. For normal incidence of a P-wave on an interface between media with impedances$Z_1 = \rho_1 V_{P1}$ and $Z_2 = \rho_2 V_{P2}$:

$$R = \frac{Z_2 - Z_1}{Z_2 + Z_1} = \frac{\rho_2 V_{P2} - \rho_1 V_{P1}}{\rho_2 V_{P2} + \rho_1 V_{P1}}$$

The transmission coefficient is:

$$T = \frac{2 Z_1}{Z_2 + Z_1}$$

The energy reflection coefficient is $R^2$, and energy conservation requires$R^2 + \frac{Z_2}{Z_1}T^2 = 1$. At the CMB, the large impedance contrast produces strong reflections (PcP, ScS phases), while at the Moho the contrast generates prominent reflected phases used in crustal studies.

2.3 Geothermal Gradient and Heat Flow

Derivation 3: Conductive Geotherm with Heat Production

In the continental crust, heat is produced by radioactive decay of U, Th, and K. The steady-state heat conduction equation with internal heat production $A$ is:

$$k \frac{d^2T}{dz^2} + A(z) = 0$$

where $k$ is thermal conductivity and $z$ is depth (positive downward). If heat production decreases exponentially with depth, $A(z) = A_0 e^{-z/h_r}$ where$h_r$ is the characteristic depth (~10 km):

$$T(z) = T_0 + \frac{q_m}{k}z + \frac{A_0 h_r^2}{k}\left(1 - e^{-z/h_r}\right)$$

where $T_0$ is surface temperature and $q_m$ is the mantle heat flow (heat flux from below the radioactive layer). The surface heat flow is:

$$q_0 = q_m + A_0 h_r$$

This is the Lachenbruch linear heat flow relation (1968): surface heat flow equals mantle heat flow plus the integrated crustal radiogenic contribution. Typical values are $q_0 = 60$-90 mW/m$^2$ for continents and $q_m \approx 20$-30 mW/m$^2$.

Derivation 4: Adiabatic Temperature Gradient in the Mantle

Below the lithosphere, heat transport is dominated by convection, and the temperature follows an adiabatic gradient. From thermodynamics:

$$\left(\frac{dT}{dz}\right)_{\text{ad}} = \frac{\alpha g T}{C_P}$$

where $\alpha$ is the thermal expansion coefficient, $g$ is gravitational acceleration, $T$ is absolute temperature, and $C_P$ is specific heat at constant pressure. For the upper mantle with $\alpha \approx 3 \times 10^{-5}$ K$^{-1}$,$g \approx 10$ m/s$^2$, $T \approx 1600$ K, and $C_P \approx 1200$ J/(kg$\cdot$K):

$$\left(\frac{dT}{dz}\right)_{\text{ad}} = \frac{3 \times 10^{-5} \times 10 \times 1600}{1200} \approx 0.4 \text{ K/km}$$

The adiabatic gradient (~0.3-0.5 K/km) is much smaller than the conductive gradient in the lithosphere (~15-25 K/km), which is why the mantle geotherm is relatively flat.

2.4 Normal Modes and Earth's Free Oscillations

Derivation 5: Eigenfrequencies of a Homogeneous Sphere

After a great earthquake, the entire Earth rings like a bell at discrete frequencies. These normal modes are standing waves described by spherical harmonics. For a homogeneous, isotropic sphere of radius $R$, the eigenfrequencies of spheroidal modes$_nS_l$ satisfy:

$$\omega_{nl}^2 = \frac{V_S^2}{R^2}\left[l(l+1) + \text{correction terms}\right]$$

where $n$ is the radial order and $l$ is the angular degree. The fundamental radial mode $_0S_0$ has a period of about 20.5 minutes, corresponding to the entire Earth expanding and contracting. The football mode $_0S_2$ has a period of about 54 minutes.

Normal mode frequencies provide the most precise constraints on Earth's radial structure. The splitting of degenerate modes (due to rotation, ellipticity, and 3D heterogeneity) constrains lateral variations in velocity and density.

Hugo Benioff first observed free oscillations after the 1960 Chilean earthquake (M$_w$ 9.5), confirming theoretical predictions by Lamb (1882) and Jeans (1923). Freeman Gilbert and George Backus developed the theoretical framework for interpreting normal mode data in the 1960s, which ultimately led to the PREM model.

Toroidal and Spheroidal Modes

Normal modes are classified as either spheroidal ($_nS_l$) or toroidal ($_nT_l$). Spheroidal modes involve radial motions (compression and dilation) and are sensitive to both P-wave and S-wave velocity structure. Toroidal modes involve purely tangential motions (torsional oscillations) and are sensitive only to S-wave structure. The fundamental modes ($n = 0$) have the simplest radial structure, while higher overtones ($n > 0$) have more complex patterns with $n$ radial nodes.

The frequency of each mode depends on Earth's density and elastic structure. By measuring thousands of mode frequencies with precision of $\sim 10^{-4}$ mHz, the radial profiles of velocity and density are determined to high accuracy. Mode coupling (interaction between modes of different type and degree due to rotation and lateral heterogeneity) provides additional constraints on 3D structure.

The attenuation of normal modes, measured by the quality factor $Q$, constrains the anelastic properties of the interior. The asthenosphere has very low $Q$(~80-120), indicating significant dissipation, likely due to partial melt. The inner core has anomalously low $Q$ for shear waves (~100), which is not fully understood but may relate to grain boundary sliding or the presence of melt along grain boundaries.

2.5 Composition and Phase Transitions

The major compositional layers of Earth are constrained by matching seismic velocities and densities to laboratory measurements at high pressure and temperature. Francis Birch (1952) established that the mantle is composed of ultramafic silicates (olivine, pyroxene, garnet) while the core is dominantly iron with ~10% light elements (O, Si, S, C, H).

The 410 km and 660 km discontinuities correspond to well-characterized phase transitions in the olivine system. At 410 km (~14 GPa), olivine ($\alpha$-phase) transforms to wadsleyite ($\beta$-phase). At 520 km, wadsleyite transforms to ringwoodite ($\gamma$-phase). At 660 km (~24 GPa), ringwoodite decomposes to bridgmanite (MgSiO$_3$ perovskite) plus ferropericlase.

The 410 km transition has a positive Clapeyron slope ($dP/dT > 0$), so it occurs at shallower depth in hot regions. The 660 km transition has a negative Clapeyron slope ($dP/dT < 0$), which can impede mantle convection across this boundary. Whether mantle convection is layered or whole-mantle remains debated, though tomographic images of subducted slabs penetrating the 660 suggest at least partial whole-mantle flow.

The Post-Perovskite Phase

Murakami et al. (2004) discovered that bridgmanite (MgSiO$_3$ perovskite) transforms to a post-perovskite phase at pressures above ~125 GPa, corresponding to conditions in the D$''$ layer. This phase transition has a steep positive Clapeyron slope (~8 MPa/K), meaning it occurs at shallower depth in cold regions (like the base of subducted slabs). The double-crossing of this phase boundary by a mantle geotherm may explain seismic observations of paired discontinuities in D$''$.

Inner Core Anisotropy

Seismic waves traveling along Earth's rotation axis through the inner core arrive about 3-4 seconds faster than waves traveling in the equatorial plane. This ~3% anisotropy suggests preferred alignment of hexagonal close-packed iron crystals, possibly due to convective or solidification textures. The inner core is also heterogeneous, with the western hemisphere showing a distinct seismic signature from the eastern hemisphere, and evidence for a distinct innermost inner core with different anisotropy. The inner core is growing at a rate of ~1 mm/yr as the outer core slowly solidifies, releasing latent heat that helps power the geodynamo.

2.6 Earth's Magnetic Field and the Geodynamo

Earth's magnetic field is generated by convection of liquid iron in the outer core, acting as a self-sustaining dynamo. The field is approximately dipolar at the surface, with current dipole moment $\sim 8 \times 10^{22}$ A$\cdot$m$^2$. The magnetic Reynolds number governs whether fluid motions can sustain a magnetic field:

$$Rm = \frac{vL}{\eta_m}$$

where $v$ is the fluid velocity, $L$ is the length scale, and $\eta_m$ is the magnetic diffusivity (~2 m$^2$/s for liquid iron). For the outer core with$v \sim 5 \times 10^{-4}$ m/s and $L \sim 2000$ km, $Rm \sim 500$, well above the critical value of ~10-50 needed for dynamo action.

The geomagnetic field reverses polarity at irregular intervals, with a mean reversal rate of ~4-5 per Myr during the current epoch. During a reversal, the dipole field weakens and the non-dipole components become relatively stronger. The process takes ~1000-10,000 years. Long intervals without reversals are called superchrons: the Cretaceous Normal Superchron lasted from 124-84 Ma (~40 Myr without reversal).

Joseph Larmor (1919) first proposed the dynamo mechanism. Walter Elsasser (1946) and Edward Bullard (1949) developed the mathematical theory. Gary Glatzmaier and Paul Roberts (1995) produced the first self-consistent numerical geodynamo simulation that exhibited spontaneous reversals.

Computational Lab: Earth's Interior

PREM Model, Geotherms, and Seismic Discontinuities

Python

script.py344 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# EARTH'S INTERIOR SIMULATION
# ============================================================
print("EARTH'S INTERIOR SIMULATION")
print("=" * 70)

# ============================================================
# 1. Simplified PREM-like Model
# ============================================================
print()
print("1. Simplified PREM Model: Velocity and Density vs Depth")
print("-" * 70)

R_earth = 6371.0  # km

# Define layers: (name, top_depth_km, bottom_depth_km, rho, Vp, Vs)
layers = [
    ("Upper Crust",    0,    15,  2.6, 5.8, 3.2),
    ("Lower Crust",   15,    35,  2.9, 6.8, 3.9),
    ("Uppermost Mantle", 35, 210, 3.4, 8.1, 4.5),
    ("Transition Zone", 210, 410, 3.5, 8.6, 4.7),
    ("TZ (wadsleyite)", 410, 520, 3.7, 9.5, 5.2),
    ("TZ (ringwoodite)", 520, 660, 3.9, 10.3, 5.6),
    ("Lower Mantle",  660, 2740,  4.9, 13.0, 7.0),
    ("D'' layer",    2740, 2891,  5.5, 13.7, 7.3),
    ("Outer Core",   2891, 5150, 10.5, 10.0, 0.0),
    ("Inner Core",   5150, 6371, 13.0, 11.2, 3.6),
]

print(f"  {'Layer':22} {'Depth (km)':16} {'rho (g/cc)':12} {'Vp (km/s)':12} {'Vs (km/s)':12} {'Z (kg/m2s)'}")
print("  " + "-" * 80)

for name, ztop, zbot, rho, vp, vs in layers:
    Z = rho * 1000 * vp * 1000  # impedance in SI
    print(f"  {name:20} {ztop:5}-{zbot:8} {rho:10.1f} {vp:10.1f} {vs:10.1f} {Z:12.2e}")

# ============================================================
# 2. Adams-Williamson Density Integration
# ============================================================
print()
print()
print("2. Adams-Williamson Density Integration (Mantle)")
print("-" * 60)

# Simple integration through the mantle
depths = np.linspace(35, 2891, 200)
rho_surf = 3.38  # density at top of mantle (g/cc -> kg/m^3 * 1000)
g_avg = 10.0  # m/s^2 (approximately constant in mantle)

# Approximate Vp and Vs profiles (simplified)
def vp_profile(z):
    if z < 410:
        return 8.0 + z * 0.004
    elif z < 660:
        return 9.5 + (z - 410) * 0.003
    else:
        return 10.75 + (z - 660) * 0.0013

def vs_profile(z):
    if z < 410:
        return 4.5 + z * 0.002
    elif z < 660:
        return 5.0 + (z - 660) * 0.002
    else:
        return 5.95 + (z - 660) * 0.0006

rho_aw = [rho_surf]
dz = (depths[1] - depths[0]) * 1000  # meters

print(f"  {'Depth (km)':14} {'Vp (km/s)':12} {'Vs (km/s)':12} {'phi (km2/s2)':14} {'rho (kg/m3)'}")
print("  " + "-" * 60)

for i in range(1, len(depths)):
    z = depths[i]
    vp = vp_profile(z)
    vs = vs_profile(z)
    phi = (vp**2 - 4.0/3.0 * vs**2)  # km^2/s^2
    phi_si = phi * 1e6  # m^2/s^2
    drho = -rho_aw[-1] * 1000 * g_avg / phi_si * dz / 1000  # kg/m^3 -> g/cc
    rho_aw.append(rho_aw[-1] + drho)

if i % 25 == 0:
        print(f"  {z:12.0f} {vp:10.2f} {vs:10.2f} {phi:12.2f} {rho_aw[-1]*1000:10.0f}")

# ============================================================
# 3. Geotherm Calculation
# ============================================================
print()
print()
print("3. Continental Geotherm")
print("-" * 60)

k_cond = 3.0  # W/(m*K)
T0 = 15.0     # surface temperature (C)
A0 = 2.5e-6   # surface heat production (W/m^3)
hr = 10000.0  # characteristic depth (m)
qm = 0.025    # mantle heat flow (W/m^2)
q0 = qm + A0 * hr
print(f"  Surface heat flow q0 = {q0*1000:.1f} mW/m^2")
print(f"  Mantle heat flow qm = {qm*1000:.1f} mW/m^2")
print(f"  Radiogenic contribution = {A0*hr*1000:.1f} mW/m^2")
print()

print(f"  {'Depth (km)':14} {'T_cond (C)':14} {'T_adiabat (C)':16} {'Gradient (C/km)'}")
print("  " + "-" * 55)

alpha_th = 3e-5
Cp = 1200.0
T_mantle = 1350.0  # potential temperature

for z_km in [0, 5, 10, 15, 20, 30, 40, 50, 75, 100, 150, 200, 300, 410, 660]:
    z = z_km * 1000
    # Conductive geotherm (valid in lithosphere)
    T_cond = T0 + qm/k_cond * z + A0*hr**2/k_cond * (1 - np.exp(-z/hr))
    # Adiabatic geotherm (valid in convecting mantle)
    T_ad = T_mantle + 273.15
    g_local = 10.0
    T_ad = (T_mantle + 273.15) * np.exp(alpha_th * g_local * z / Cp) - 273.15
    grad = alpha_th * g_local * (T_ad + 273.15) / Cp * 1000  # C/km
    print(f"  {z_km:12} {T_cond:12.1f} {T_ad:14.1f} {grad:10.2f}")

# ============================================================
# 4. Reflection Coefficients at Major Discontinuities
# ============================================================
print()
print()
print("4. Reflection Coefficients at Major Boundaries")
print("-" * 60)

boundaries = [
    ("Moho",       2.9, 6.8, 3.4, 8.1),
    ("410 km",     3.5, 8.9, 3.7, 9.5),
    ("660 km",     3.9, 10.3, 4.4, 10.75),
    ("CMB",        5.5, 13.7, 10.0, 8.1),
    ("ICB",       12.1, 10.3, 13.0, 11.2),
]

print(f"  {'Boundary':14} {'Z1 (MRayl)':14} {'Z2 (MRayl)':14} {'R':10} {'R^2 (%)'}")
print("  " + "-" * 55)

for name, rho1, vp1, rho2, vp2 in boundaries:
    Z1 = rho1 * vp1  # (g/cc)(km/s) = 10^6 kg/(m^2 s) = MRayl
    Z2 = rho2 * vp2
    R = (Z2 - Z1) / (Z2 + Z1)
    print(f"  {name:12} {Z1:12.2f} {Z2:12.2f} {R:8.4f} {R**2*100:8.2f}")

# ============================================================
# 5. Normal Mode Periods
# ============================================================
print()
print()
print("5. Earth Normal Mode Periods (observed)")
print("-" * 55)

modes = [
    ("0S0", "Breathing", 20.5),
    ("0S2", "Football", 54.7),
    ("0S3", "l=3 spheroidal", 35.6),
    ("1S2", "1st overtone", 26.9),
    ("0T2", "l=2 toroidal", 44.2),
    ("0S29", "Rayleigh-like", 5.7),
]

print(f"  {'Mode':10} {'Description':20} {'Period (min)':14} {'Freq (mHz)'}")
print("  " + "-" * 50)
for mode, desc, period in modes:
    freq = 1.0 / (period * 60) * 1000
    print(f"  {mode:8} {desc:18} {period:12.1f} {freq:8.3f}")

# ============================================================
# PLOTS
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# Plot 1: PREM velocity profile
ax = axes[0, 0]
for name, ztop, zbot, rho, vp, vs in layers:
    ax.plot([ztop, zbot], [vp, vp], 'b-', linewidth=2)
    ax.plot([ztop, zbot], [vs, vs], 'r-', linewidth=2)
    ax.plot([ztop, zbot], [rho, rho], 'g-', linewidth=2)

ax.plot([], [], 'b-', linewidth=2, label='Vp (km/s)')
ax.plot([], [], 'r-', linewidth=2, label='Vs (km/s)')
ax.plot([], [], 'g-', linewidth=2, label='rho (g/cc)')

for z_disc in [35, 410, 660, 2891, 5150]:
    ax.axvline(x=z_disc, color='gray', linestyle=':', alpha=0.3)

ax.set_xlabel('Depth (km)')
ax.set_ylabel('Velocity (km/s) or Density (g/cc)')
ax.set_title('Simplified PREM Model')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 6371)

# Plot 2: Geotherm
ax = axes[0, 1]
z_crust = np.linspace(0, 200, 200)
T_cond_arr = T0 + qm/k_cond * z_crust*1000 + A0*hr**2/k_cond * (1 - np.exp(-z_crust*1000/hr))

z_mantle = np.linspace(0, 660, 200)
T_ad_arr = (T_mantle + 273.15) * np.exp(alpha_th * 10.0 * z_mantle*1000 / Cp) - 273.15

ax.plot(T_cond_arr, z_crust, 'r-', linewidth=2, label='Conductive geotherm')
ax.plot(T_ad_arr, z_mantle, 'b-', linewidth=2, label='Mantle adiabat')

# Solidus (approximate)
z_sol = np.linspace(0, 660, 100)
T_sol = 1100 + z_sol * 3.5  # very simplified solidus
ax.plot(T_sol, z_sol, 'k--', linewidth=1.5, label='Dry solidus (approx)')

ax.set_xlabel('Temperature (C)')
ax.set_ylabel('Depth (km)')
ax.set_title('Continental Geotherm')
ax.invert_yaxis()
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 3500)
ax.set_ylim(660, 0)

# Plot 3: Adams-Williamson density
ax = axes[1, 0]
ax.plot(np.array(rho_aw)*1000, depths, 'b-', linewidth=2, label='Adams-Williamson')

# PREM reference points
prem_z = [35, 210, 410, 660, 2740, 2891]
prem_rho = [3400, 3500, 3700, 4400, 5500, 5500]
ax.plot(prem_rho, prem_z, 'ro', markersize=8, label='PREM reference')

ax.set_xlabel('Density (kg/m^3)')
ax.set_ylabel('Depth (km)')
ax.set_title('Density Profile: Adams-Williamson')
ax.invert_yaxis()
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 4: Pressure and gravity profiles
ax = axes[1, 1]
z_full = np.linspace(0, 6371, 500)
g_profile = np.zeros_like(z_full)
P_profile = np.zeros_like(z_full)

# Simplified gravity (roughly PREM-like)
for i, z in enumerate(z_full):
    r = R_earth - z
    if r > 0:
        # Approximate g(r): increases slightly in mantle, peaks near CMB
        if z < 2891:
            g_profile[i] = 9.81 + z * 0.0003
        elif z < 5150:
            g_profile[i] = 10.68 - (z - 2891) * 0.0025
        else:
            g_profile[i] = 5.05 - (z - 5150) * 0.0041

# Pressure integration
dz_p = (z_full[1] - z_full[0]) * 1000
for i in range(1, len(z_full)):
    z = z_full[i]
    if z < 35:
        rho_p = 2800
    elif z < 660:
        rho_p = 3500 + z * 0.5
    elif z < 2891:
        rho_p = 4500 + (z-660) * 0.5
    elif z < 5150:
        rho_p = 10500
    else:
        rho_p = 13000
    P_profile[i] = P_profile[i-1] + rho_p * g_profile[i] * dz_p

P_GPa = P_profile / 1e9

ax2 = ax.twiny()
ax.plot(g_profile, z_full, 'b-', linewidth=2, label='g (m/s^2)')
ax2.plot(P_GPa, z_full, 'r-', linewidth=2, label='P (GPa)')

ax.set_xlabel('Gravity (m/s^2)', color='blue')
ax2.set_xlabel('Pressure (GPa)', color='red')
ax.set_ylabel('Depth (km)')
ax.set_title('Gravity and Pressure vs Depth')
ax.invert_yaxis()
ax.tick_params(axis='x', labelcolor='blue')
ax2.tick_params(axis='x', labelcolor='red')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Magnetic Field Properties")
print("-" * 55)

Rm_core = 500  # magnetic Reynolds number
v_core = 5e-4  # m/s
L_core = 2000e3  # m
eta_mag = 2.0  # m^2/s

print(f"  Magnetic Reynolds number: Rm = {Rm_core}")
print(f"  Core flow velocity: {v_core*1000:.1f} mm/s")
print(f"  Magnetic diffusivity: {eta_mag} m^2/s")
print(f"  Ohmic decay time: {L_core**2/eta_mag/(3.156e7*1000):.0f} kyr")
print()
print("  Geomagnetic Reversal Statistics:")
print("  " + "-" * 45)

reversal_periods = [
    ("Present epoch", 4.5, "reversals/Myr"),
    ("Cretaceous superchron", 0.0, "84-124 Ma"),
    ("Late Jurassic", 3.0, "reversals/Myr"),
    ("Permo-Triassic superchron", 0.0, "250-310 Ma"),
]

print(f"  {'Period':30} {'Rate':18} {'Notes'}")
print("  " + "-" * 55)
for name, rate, notes in reversal_periods:
    print(f"  {name:28} {rate:16.1f} {notes}")

# Dipole field calculation
M_dipole = 8e22  # A.m^2
mu_0 = 4e-7 * np.pi
R_E = 6.371e6

print()
print("  Magnetic field strength vs latitude:")
print(f"  {'Latitude':12} {'B_r (uT)':14} {'B_theta (uT)':16} {'|B| (uT)'}")
print("  " + "-" * 50)

for lat_m in [0, 15, 30, 45, 60, 75, 90]:
    theta_m = np.radians(90 - lat_m)
    Br = mu_0 * M_dipole * 2 * np.cos(theta_m) / (4 * np.pi * R_E**3) * 1e6
    Bt = mu_0 * M_dipole * np.sin(theta_m) / (4 * np.pi * R_E**3) * 1e6
    B_total = np.sqrt(Br**2 + Bt**2)
    print(f"  {lat_m:10} {Br:12.1f} {Bt:14.1f} {B_total:10.1f}")

print()
print()
print("[Chart saved: PREM profile, geotherm, density, gravity/pressure]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Phase Transitions and Mineral Physics at High Pressure

Python

script.py72 lines

import numpy as np

# ============================================================
# MINERAL PHYSICS AT HIGH PRESSURE
# ============================================================
print("MINERAL PHYSICS: HIGH PRESSURE PHASE TRANSITIONS")
print("=" * 70)

# ============================================================
# 1. Birch-Murnaghan Equation of State
# ============================================================
print()
print("1. Birch-Murnaghan EOS: Volume vs Pressure")
print("-" * 55)

# Third-order Birch-Murnaghan for MgSiO3 bridgmanite
K0 = 261.0  # bulk modulus (GPa)
K0p = 4.0   # pressure derivative
V0 = 24.45  # molar volume (cm^3/mol)

print(f"  K0 = {K0} GPa, K0' = {K0p}, V0 = {V0} cm^3/mol")
print()
print(f"  {'P (GPa)':12} {'V/V0':10} {'rho/rho0':12} {'K (GPa)':12} {'Depth (km)'}")
print("  " + "-" * 55)

for P in [0, 25, 50, 75, 100, 125, 136]:
    # Iteratively solve BM3 for V/V0
    x = 1.0  # V/V0 guess
    for iteration in range(100):
        f = 0.5 * ((x)**(-2.0/3.0) - 1)
        P_calc = 3*K0*f*(1 + 2*f)**(5.0/2.0) * (1 + 1.5*(K0p - 4)*f)
        dPdx = -K0 * x**(-5.0/3.0) * (1 + 2*f)**(5.0/2.0) * (1 + 1.5*(K0p-4)*f)
        x = x - (P_calc - P) / dPdx
        x = max(x, 0.5)
        if abs(P_calc - P) < 0.01:
            break
    rho_ratio = 1.0 / x
    K_eff = K0 * (1 + 2*f)**(5.0/2.0) * (1 + (3*K0p - 5)*f)
    depth = P / 0.033  # rough conversion
    print(f"  {P:10} {x:8.4f} {rho_ratio:10.4f} {K_eff:10.1f} {depth:8.0f}")

# ============================================================
# 2. Clapeyron Slopes of Major Transitions
# ============================================================
print()
print()
print("2. Major Phase Transitions in the Mantle")
print("-" * 65)

transitions = [
    ("Olivine -> Wadsleyite", 410, 13.5, 2.5, "+", "Exothermic"),
    ("Wadsleyite -> Ringwoodite", 520, 18.0, 5.0, "+", "Exothermic"),
    ("Ringw -> Brid + Fp", 660, 23.5, -2.5, "-", "Endothermic"),
    ("Post-perovskite", 2700, 125, 8.0, "+", "Exothermic"),
]

print(f"  {'Transition':28} {'Depth (km)':14} {'P (GPa)':10} {'dP/dT (MPa/K)':16} {'Type'}")
print("  " + "-" * 75)

for name, depth, P, slope, sign, thtype in transitions:
    print(f"  {name:26} {depth:12} {P:8.1f} {sign}{abs(slope):14.1f} {thtype}")

# Effect on convection
print()
print("  Effect on mantle convection:")
print("    410 km (+slope): assists downgoing slabs (deepens in cold)")
print("    660 km (-slope): resists slab penetration (shallows in cold)")
print("    Post-perovskite: may enhance D'' instabilities")

print()
print("Complete: mineral physics calculations finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Seismology Next: Geodynamics →