Part I: Solid Earth | Chapter 3

Igneous Petrology

From magma genesis to crystalline rock: the physics and chemistry of melting

3.1 Magma Generation

Magma is generated when rocks in the mantle or crust partially melt. Three principal mechanisms lower the solidus or raise the local temperature above it: decompression melting, flux melting (addition of volatiles), and heat addition.

Solidus Depression by Water

Water dramatically lowers the solidus temperature of mantle peridotite. The dry peridotite solidus is approximately:

$$T_{\text{solidus}}^{\text{dry}}(P) \approx 1100 + 3.5 P \quad (°\text{C}, P \text{ in kbar})$$

With water present, the solidus drops by several hundred degrees. The wet peridotite solidus can be approximated as:

$$T_{\text{solidus}}^{\text{wet}}(P) \approx 1000 + 0.1 P^2 - 5P \quad (°\text{C}, P \text{ in kbar})$$

The depression $\Delta T = T^{\text{dry}} - T^{\text{wet}}$ can exceed 300°C at subduction zone pressures. This is why island arc volcanism occurs — water released from the subducting slab lowers the solidus of the overlying mantle wedge.

Derivation: Decompression Melting

At mid-ocean ridges, melting occurs because upwelling mantle follows an adiabatic path that crosses the solidus. The mantle adiabatic gradient is:

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

where $\alpha \approx 3 \times 10^{-5}$ K$^{-1}$ is thermal expansion, $g = 9.8$ m/s$^2$,$T \approx 1600$ K is potential temperature, and $c_p \approx 1200$ J kg$^{-1}$ K$^{-1}$. This gives:

$$\left(\frac{dT}{dz}\right)_{\text{adiabat}} \approx \frac{3 \times 10^{-5} \times 9.8 \times 1600}{1200} \approx 0.4 \text{ °C/km}$$

The solidus gradient is steeper ($\sim 3$ °C/km near the surface), so the adiabat eventually crosses the solidus during ascent. The depth of initial melting depends on mantle potential temperature $T_p$:

$$z_{\text{melt}} = \frac{T_p - T_{\text{solidus}}(z=0)}{\left(\frac{dT_{\text{solidus}}}{dz}\right) - \left(\frac{dT}{dz}\right)_{\text{adiabat}}}$$

For $T_p = 1350$°C and a surface solidus of ~1100°C with a solidus gradient of 3°C/km:

$$z_{\text{melt}} \approx \frac{1350 - 1100}{3 - 0.4} \approx 96 \text{ km}$$

The melt fraction increases as the mantle continues to rise above this depth. The isentropic productivity is approximately:

$$\left(\frac{\partial F}{\partial P}\right)_S \approx \frac{1}{T} \frac{\Delta S_{\text{fus}}}{c_p} \frac{dT_{\text{solidus}}}{dP} \approx 1\text{-}2 \%/\text{kbar}$$

3.2 Fractional Crystallization

Bowen's Reaction Series

N.L. Bowen (1928) established that minerals crystallize from cooling magma in a predictable sequence. The series has two branches:

Discontinuous (Ferromagnesian)

Olivine (~1800°C)
Pyroxene (~1200°C)
Amphibole (~1100°C)
Biotite (~900°C)

Continuous (Plagioclase)

Ca-rich (anorthite) (~1500°C)
Intermediate
Na-rich (albite) (~1100°C)

Both converge to: K-feldspar, Muscovite, Quartz (~700°C)

Derivation: Rayleigh Fractionation

When crystals are immediately removed from the melt (perfect fractional crystallization), the evolution of trace element concentrations follows the Rayleigh distillation law.

Consider a trace element with bulk partition coefficient $D = C_s / C_L$, where$C_s$ is the concentration in the solid and $C_L$ in the liquid. For an infinitesimal amount of crystallization $dF$ (where $F$ is melt fraction remaining):

$$C_L F = C_L(F + dF)(F + dF) + D \cdot C_L \cdot (-dF)$$

Mass balance on the trace element during crystallization of mass $dm$:

$$d(C_L F) = -D \cdot C_L \cdot dF$$$$F \, dC_L + C_L \, dF = -D \cdot C_L \, dF$$$$F \, dC_L = -(D-1) C_L \, dF$$$$\frac{dC_L}{C_L} = (D-1) \frac{dF}{F}$$

Note: since $F$ decreases as crystallization proceeds, and we lose mass $dF < 0$ from the liquid, we write $dF$ as a decrease. Integrating from $F = 1$ (initial) to $F$:

$$\int_{C_0}^{C_L} \frac{dC_L}{C_L} = (D-1) \int_1^F \frac{dF}{F}$$$$\ln\frac{C_L}{C_0} = (D-1) \ln F$$

$$\boxed{C_L = C_0 F^{(D-1)}}$$

Key behaviors:

Incompatible elements ($D \ll 1$): $C_L/C_0 = F^{(D-1)} \approx F^{-1}$ — concentration increases dramatically as $F \to 0$. Examples: Rb, Ba, Nb, Zr
Compatible elements ($D \gg 1$): $C_L/C_0 = F^{(D-1)} \to 0$ quickly — element is depleted from melt. Examples: Ni ($D \approx 10$ in olivine), Cr
Perfectly compatible ($D = 1$): $C_L = C_0$ — no fractionation

For batch (equilibrium) melting, the contrasting equation is:

$$C_L = \frac{C_0}{D + F(1 - D)}$$

3.3 Binary Phase Diagrams

Eutectic Systems: Diopside-Anorthite

The diopside (CaMgSi$_2$O$_6$) — anorthite (CaAl$_2$Si$_2$O$_8$) system is a classic eutectic. These two minerals do not form solid solutions (simplified), so the phase diagram has two liquidus curves meeting at the eutectic point.

The liquidus temperature for each component is depressed by the presence of the other, described approximately by the freezing point depression:

$$\ln X_A = -\frac{\Delta H_{\text{fus},A}}{R}\left(\frac{1}{T} - \frac{1}{T_{m,A}}\right)$$

where $X_A$ is the mole fraction, $\Delta H_{\text{fus}}$ is the enthalpy of fusion, and $T_{m,A}$ is the melting temperature of pure A. The eutectic is where both curves intersect. For Di-An:

Diopside melting point: $T_m^{\text{Di}} = 1392$°C
Anorthite melting point: $T_m^{\text{An}} = 1553$°C
Eutectic: ~1274°C at ~42 wt% An (58% Di)

Solid Solution: The Olivine System

The forsterite-fayalite system exhibits complete solid solution. Both liquidus and solidus curves span the composition range. At any temperature between liquidus and solidus, the equilibrium compositions of melt and crystal are given by the lever rule:

$$\frac{\text{mass of liquid}}{\text{mass of solid}} = \frac{X_{\text{bulk}} - X_{\text{solid}}}{X_{\text{liquid}} - X_{\text{bulk}}}$$

Forsterite ($T_m = 1890$°C) crystallizes at higher temperatures than fayalite ($T_m = 1205$°C), so the first olivine to crystallize from a melt is always more Mg-rich than the bulk composition.

Peritectic Systems

In a peritectic system, an early-formed phase reacts with the liquid to produce a new phase. The classic example is the forsterite-silica system, where enstatite forms by peritectic reaction:

$$\text{Mg}_2\text{SiO}_4 \text{ (forsterite)} + \text{SiO}_2 \text{ (melt)} \rightarrow 2\text{MgSiO}_3 \text{ (enstatite)}$$

At the peritectic point, three phases coexist (forsterite + enstatite + liquid), giving$F = 2 - 3 + 1 = 0$ (invariant at constant pressure). The peritectic temperature is ~1557°C.

3.4 IUGS Classification

The TAS Diagram

The Total Alkali-Silica (TAS) diagram is the standard IUGS classification for volcanic rocks. It plots $(\text{Na}_2\text{O} + \text{K}_2\text{O})$ vs. $\text{SiO}_2$ (both in weight percent). Key boundaries:

Ultrabasic: SiO$_2$ < 45% (picrite, komatiite)
Basic: 45-52% (basalt, basanite)
Intermediate: 52-63% (andesite, trachyandesite)
Acid: >63% (dacite, rhyolite, trachyte)

The alkali-silica dividing line between alkaline and subalkaline series is approximated by:

$$(\text{Na}_2\text{O} + \text{K}_2\text{O}) = 0.37 + 0.14 \times \text{SiO}_2 - 0.0005 \times \text{SiO}_2^2$$

CIPW Normative Mineralogy

The CIPW norm (Cross, Iddings, Pirsson, Washington) converts a bulk chemical analysis into a theoretical mineral assemblage assuming anhydrous, equilibrium crystallization. The algorithm proceeds in a fixed sequence:

Convert oxide weight percentages to molecular proportions: $n_i = \text{wt\%}_i / M_i$
Allocate accessory minerals (apatite from P$_2$O$_5$, ilmenite from TiO$_2$)
Distribute Al$_2$O$_3$: first to orthoclase (with K$_2$O), then albite (Na$_2$O), then anorthite (CaO)
Remaining CaO goes to diopside (with MgO, FeO)
Remaining MgO + FeO forms hypersthene or olivine
Excess or deficit SiO$_2$ determines: quartz (oversaturated) vs nepheline (undersaturated)

The presence of normative quartz or nepheline is a fundamental classification criterion: quartz-normative rocks are silica-oversaturated (tholeiitic basalts), while nepheline-normative rocks are silica-undersaturated (alkaline basalts).

3.5 Magma Mixing and Crustal Assimilation

Binary Magma Mixing

When two magmas of different compositions mix, the resulting hybrid composition for any element or isotope ratio lies on a straight line in element-element space (for elements) or a hyperbola in isotope-element space. For a fraction $f$ of magma B mixed with $(1-f)$ of magma A:

$$C_{\text{mix}} = f \cdot C_B + (1-f) \cdot C_A$$

For isotope ratios (e.g., $^{87}\text{Sr}/^{86}\text{Sr}$), the mixing relationship becomes hyperbolic because isotope ratios are weighted by elemental concentrations:

$$\left(\frac{^{87}\text{Sr}}{^{86}\text{Sr}}\right)_{\text{mix}} = \frac{f \cdot C_B^{\text{Sr}} \cdot R_B + (1-f) \cdot C_A^{\text{Sr}} \cdot R_A}{f \cdot C_B^{\text{Sr}} + (1-f) \cdot C_A^{\text{Sr}}}$$

where $R_A$ and $R_B$ are the isotope ratios of end-members A and B.

AFC: Assimilation-Fractional Crystallization

DePaolo (1981) developed the AFC model to describe the simultaneous assimilation of country rock and fractional crystallization. The key parameter is:

$$r = \frac{\dot{M}_a}{\dot{M}_c}$$

where $\dot{M}_a$ is the rate of assimilation and $\dot{M}_c$ is the rate of crystallization. The trace element evolution follows:

$$\frac{C_L}{C_0} = F^{-z} + \frac{r}{r-1} \frac{C_a}{z C_0}\left(1 - F^{-z}\right)$$

where $z = (r + D - 1)/(r - 1)$, $F$ is the remaining melt fraction,$C_a$ is the concentration in the assimilant, and $D$ is the bulk partition coefficient. This model explains why granites have elevated $^{87}\text{Sr}/^{86}\text{Sr}$ratios compared to their basaltic parents.

Texture and Cooling Rate

The texture of an igneous rock records its cooling history. Crystal nucleation rate$J$ and growth rate $G$ both depend on the degree of undercooling$\Delta T = T_{\text{liquidus}} - T$:

$$J \propto \exp\left(-\frac{\Delta G^*}{kT}\right) \cdot \exp\left(-\frac{E_a}{kT}\right)$$

where $\Delta G^*$ is the activation energy for nucleation and $E_a$ is the activation energy for diffusion. At small undercooling, growth dominates and few large crystals form (phaneritic texture). At large undercooling, nucleation dominates and many small crystals form (aphanitic texture). Extreme undercooling produces glass.

Phaneritic: Slow cooling in plutonic environment (granite, gabbro)
Aphanitic: Rapid cooling at surface (basalt, rhyolite)
Porphyritic: Two-stage cooling — large phenocrysts in fine groundmass
Glassy: Very rapid quenching (obsidian, pumice)
Pegmatitic: Very slow cooling with fluxing volatiles ($> 1$ cm crystals)

3.6 Volcanic Eruption Dynamics

Magma Viscosity

The viscosity of silicate melts varies over many orders of magnitude and controls eruption style. The Arrhenius relationship gives:

$$\eta = A \exp\left(\frac{E_a}{RT}\right)$$

Viscosity increases with SiO$_2$ content (more polymerized melt), decreasing temperature, and increasing crystal content. Water dramatically reduces viscosity by breaking Si-O-Si bridges. Typical values:

Basalt at 1200°C: $\eta \approx 10^1 - 10^2$ Pa·s
Andesite at 1000°C: $\eta \approx 10^3 - 10^4$ Pa·s
Rhyolite at 800°C (dry): $\eta \approx 10^8 - 10^{12}$ Pa·s
Rhyolite at 800°C (6 wt% H$_2$O): $\eta \approx 10^4 - 10^5$ Pa·s

Conduit Flow and Eruption Rate

For steady Poiseuille flow through a cylindrical conduit of radius $R_c$ and length $L$:

$$Q = \frac{\pi R_c^4}{8\eta} \frac{\Delta P}{L}$$

where $\Delta P = (\rho_{\text{rock}} - \rho_{\text{magma}})gL + P_{\text{chamber}}$ is the driving pressure. The strong dependence on conduit radius ($R_c^4$) means small changes in conduit geometry profoundly affect eruption rate. The mass eruption rate is:

$$\dot{M} = \rho_{\text{magma}} Q = \frac{\pi \rho R_c^4 \Delta P}{8\eta L}$$

Eruption styles are classified by their Volcanic Explosivity Index (VEI), which correlates with magma viscosity and volatile content. Effusive eruptions (VEI 0-1) characterize low-viscosity basalts; Plinian eruptions (VEI 5-8) involve highly viscous, volatile-rich magmas that fragment explosively.

Computational Simulations

Binary Eutectic Phase Diagram: Diopside-Anorthite

Python

script.py95 lines

import numpy as np

# Binary eutectic: Diopside - Anorthite
# Using simplified liquidus curves from thermodynamic data

R = 8.314  # J/mol/K

# Melting points and enthalpies of fusion
Tm_Di = 1392 + 273.15  # K
Tm_An = 1553 + 273.15  # K
dH_Di = 138000.0       # J/mol (enthalpy of fusion)
dH_An = 130000.0       # J/mol

# Liquidus: ln(X_A) = -(dH/R)(1/T - 1/Tm)
# Solve for T given X

def liquidus_T(X, Tm, dH):
    """Temperature on liquidus for component with mole fraction X."""
    if X <= 0 or X > 1:
        return None
    inv_T = 1.0/Tm - R * np.log(X) / dH
    return 1.0 / inv_T

print("Diopside-Anorthite Binary Eutectic Phase Diagram")
print("=" * 65)
print(f"  Tm(Di) = {Tm_Di-273.15:.0f} C, dH(Di) = {dH_Di/1000:.0f} kJ/mol")
print(f"  Tm(An) = {Tm_An-273.15:.0f} C, dH(An) = {dH_An/1000:.0f} kJ/mol")
print()

# Compute liquidus curves
print(f"{'X_An':<10} {'T_Di_liq (C)':<15} {'T_An_liq (C)':<15} {'Stable liq. (C)'}")
print("-" * 55)

eutectic_X = None
eutectic_T = 3000.0

compositions = np.arange(0.0, 1.01, 0.05)
for X_An in compositions:
    X_Di = 1.0 - X_An
    T_Di = liquidus_T(X_Di, Tm_Di, dH_Di) - 273.15 if X_Di > 0.001 else None
    T_An = liquidus_T(X_An, Tm_An, dH_An) - 273.15 if X_An > 0.001 else None

if T_Di is not None and T_An is not None:
        T_stable = min(T_Di, T_An)
        T_liq = max(T_Di, T_An)
        print(f"  {X_An:<8.2f} {T_Di:<13.0f} {T_An:<13.0f} {T_liq:<10.0f}")

# Find eutectic (where the two curves are closest)
        diff = abs(T_Di - T_An)
        if diff < abs(eutectic_T - 0):
            eutectic_T = (T_Di + T_An) / 2
            eutectic_X = X_An
    elif T_Di is not None:
        print(f"  {X_An:<8.2f} {T_Di:<13.0f} {'---':<13} {T_Di:<10.0f}")
    elif T_An is not None:
        print(f"  {X_An:<8.2f} {'---':<13} {T_An:<13.0f} {T_An:<10.0f}")

# Better eutectic search
best_diff = 1e6
for X_An in np.linspace(0.01, 0.99, 10000):
    X_Di = 1.0 - X_An
    T_Di = liquidus_T(X_Di, Tm_Di, dH_Di)
    T_An = liquidus_T(X_An, Tm_An, dH_An)
    if T_Di is not None and T_An is not None:
        diff = abs(T_Di - T_An)
        if diff < best_diff:
            best_diff = diff
            eutectic_X = X_An
            eutectic_T = (T_Di + T_An) / 2 - 273.15

print()
print(f"Eutectic point: X_An = {eutectic_X:.3f} ({eutectic_X*100:.1f}% An), T = {eutectic_T:.0f} C")
print()

# Cooling path example
print("Cooling Path: X_An = 0.70 (70% Anorthite)")
print("-" * 50)
X_bulk = 0.70
T_liq = liquidus_T(X_bulk, Tm_An, dH_An) - 273.15
print(f"  Liquidus temperature: {T_liq:.0f} C (Anorthite begins to crystallize)")
print(f"  First crystal: Pure anorthite")
print(f"  As T decreases, anorthite crystallizes, liquid moves toward eutectic")
print(f"  At eutectic ({eutectic_T:.0f} C): both Di and An crystallize simultaneously")
print(f"  Final rock: {X_bulk*100:.0f}% An + {(1-X_bulk)*100:.0f}% Di")

# Lever rule at intermediate T
T_test = (T_liq + eutectic_T) / 2
# At T_test, find X_liquid on An liquidus
# ln(X_An_liq) = -(dH_An/R)(1/T - 1/Tm_An)
T_K = T_test + 273.15
X_An_liq = np.exp(-(dH_An/R) * (1.0/T_K - 1.0/Tm_An))
frac_liquid = (X_bulk - 1.0) / (X_An_liq - 1.0) if X_An_liq != 1.0 else 1.0
frac_liquid = min(1.0, max(0.0, (1.0 - X_bulk) / (1.0 - X_An_liq))) if X_An_liq < 1 else 1.0
print(f"  At T = {T_test:.0f} C: X_liquid = {X_An_liq:.3f}, fraction liquid ~ {frac_liquid:.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Rayleigh Fractionation Curves for Trace Elements

Python

script.py79 lines

import numpy as np

# Rayleigh fractionation: C_L = C_0 * F^(D-1)
# Batch melting: C_L = C_0 / (D + F(1-D))

print("Rayleigh Fractionation: Trace Element Evolution")
print("=" * 70)

# Partition coefficients for basaltic system
elements = {
    'Ni':  {'D': 10.0, 'type': 'highly compatible'},
    'Cr':  {'D': 4.0,  'type': 'compatible'},
    'Sc':  {'D': 1.5,  'type': 'slightly compatible'},
    'Y':   {'D': 0.5,  'type': 'slightly incompatible'},
    'Zr':  {'D': 0.1,  'type': 'incompatible'},
    'Rb':  {'D': 0.01, 'type': 'highly incompatible'},
    'Nb':  {'D': 0.005,'type': 'very incompatible'},
    'Th':  {'D': 0.001,'type': 'extremely incompatible'},
}

print(f"{'Element':<10} {'D':<8} {'Type':<25} | C_L/C_0 at F =")
print(f"{'':>43} | {'0.9':<8} {'0.7':<8} {'0.5':<8} {'0.3':<8} {'0.1':<8} {'0.01'}")
print("-" * 100)

for elem, props in elements.items():
    D = props['D']
    ratios = []
    for F in [0.9, 0.7, 0.5, 0.3, 0.1, 0.01]:
        C_ratio = F**(D - 1)
        ratios.append(C_ratio)
    ratio_str = ' '.join([f'{r:<8.3f}' if r < 100 else f'{r:<8.1f}' for r in ratios])
    print(f"  {elem:<8} {D:<8.3f} {props['type']:<23} | {ratio_str}")

print()
print()
print("Comparison: Rayleigh vs Batch Melting (D = 0.01, highly incompatible)")
print("=" * 60)
D = 0.01
print(f"{'F (melt frac)':<16} {'C_L/C_0 (Rayleigh)':<22} {'C_L/C_0 (Batch)':<22} {'Ratio'}")
print("-" * 65)

for F in [1.0, 0.8, 0.6, 0.4, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005]:
    rayleigh = F**(D - 1)
    batch = 1.0 / (D + F * (1 - D))
    ratio = rayleigh / batch if batch > 0 else 0
    print(f"  {F:<14.3f} {rayleigh:<20.2f} {batch:<20.2f} {ratio:<8.2f}")

print()
print("Key insight: At low melt fractions, Rayleigh produces much higher")
print("enrichments than batch melting for incompatible elements.")

print()
print()
print("Crystal Fractionation Modeling: Basalt to Rhyolite")
print("=" * 60)
print("Major element evolution (simplified, wt%):")
print(f"{'F':<8} {'SiO2':<10} {'MgO':<10} {'FeO':<10} {'CaO':<10} {'Na2O+K2O'}")
print("-" * 50)

# Simplified major element fractionation trends
# Starting basalt composition
for F in [1.0, 0.8, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05]:
    # Approximate SiO2 enrichment during differentiation
    SiO2 = 50 + 25 * (1 - F)**1.2
    MgO = 8.0 * F**2.5
    FeO = 10.0 * F**0.8
    CaO = 11.0 * F**1.2
    alkali = 3.0 * F**(-0.5) if F > 0.05 else 13.4
    alkali = min(alkali, 15.0)

rock = ""
    if SiO2 < 52: rock = "basalt"
    elif SiO2 < 57: rock = "basaltic andesite"
    elif SiO2 < 63: rock = "andesite"
    elif SiO2 < 69: rock = "dacite"
    else: rock = "rhyolite"

print(f"  {F:<6.2f} {SiO2:<10.1f} {MgO:<10.1f} {FeO:<10.1f} {CaO:<10.1f} {alkali:<8.1f}  ({rock})")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

TAS Classification of Volcanic Rocks

Python

script.py92 lines

import numpy as np

# Total Alkali-Silica (TAS) diagram classification
# Based on Le Bas et al. (1986)

def classify_TAS(sio2, na2o_k2o):
    """Classify volcanic rock using TAS diagram."""
    alk = na2o_k2o
    si = sio2

# Simplified TAS classification boundaries
    if si < 41:
        if alk > 3: return "Foidite"
        return "Picrobasalt"
    elif si < 45:
        if alk > 5: return "Foidite"
        if alk > 3: return "Tephrite/Basanite"
        return "Picrobasalt"
    elif si < 52:
        if alk > 9.3: return "Phonotephrite"
        if alk > 7: return "Tephriphonolite"
        if alk > 5: return "Trachybasalt"
        if alk > 3.3: return "Basaltic trachyandesite"
        return "Basalt"
    elif si < 57:
        if alk > 11.5: return "Phonolite"
        if alk > 9.3: return "Tephriphonolite"
        if alk > 7: return "Trachyandesite"
        if alk > 5: return "Basaltic trachyandesite"
        return "Basaltic andesite"
    elif si < 63:
        if alk > 11.5: return "Phonolite"
        if alk > 9.3: return "Trachyte/Trachydacite"
        if alk > 7: return "Trachyandesite"
        return "Andesite"
    elif si < 69:
        if alk > 11.5: return "Phonolite"
        if alk > 7.3: return "Trachyte/Trachydacite"
        return "Dacite"
    else:
        if alk > 7.3: return "Trachyte/Trachydacite"
        return "Rhyolite"

# Alkali-subalkaline dividing line (Irvine & Baragar, 1971)
def is_alkaline(sio2, na2o_k2o):
    boundary = 0.37 + 0.14 * sio2 - 0.0005 * sio2**2
    return na2o_k2o > boundary

# Sample volcanic rock analyses
samples = [
    {"name": "MORB (Mid-Ocean Ridge)", "SiO2": 50.5, "Na2O": 2.5, "K2O": 0.2},
    {"name": "Hawaii tholeiite", "SiO2": 49.5, "Na2O": 2.3, "K2O": 0.5},
    {"name": "Hawaii alkali basalt", "SiO2": 47.0, "Na2O": 3.0, "K2O": 1.2},
    {"name": "Andes andesite", "SiO2": 58.0, "Na2O": 3.8, "K2O": 1.5},
    {"name": "Mt St Helens dacite", "SiO2": 64.5, "Na2O": 4.5, "K2O": 1.5},
    {"name": "Yellowstone rhyolite", "SiO2": 76.5, "Na2O": 3.8, "K2O": 5.0},
    {"name": "Vesuvius phonolite", "SiO2": 55.0, "Na2O": 6.5, "K2O": 6.0},
    {"name": "Etna trachybasalt", "SiO2": 48.0, "Na2O": 3.5, "K2O": 2.0},
    {"name": "Iceland picrite", "SiO2": 43.0, "Na2O": 1.5, "K2O": 0.3},
    {"name": "Cascades basalt", "SiO2": 51.0, "Na2O": 3.0, "K2O": 0.8},
    {"name": "Japan arc andesite", "SiO2": 60.0, "Na2O": 3.5, "K2O": 1.8},
    {"name": "Bishop Tuff rhyolite", "SiO2": 77.0, "Na2O": 3.5, "K2O": 5.2},
]

print("TAS Classification of Volcanic Rock Samples")
print("=" * 85)
print(f"{'Sample':<28} {'SiO2':<8} {'Alk':<8} {'TAS Name':<25} {'Series'}")
print("-" * 85)

for s in samples:
    alk = s["Na2O"] + s["K2O"]
    tas_name = classify_TAS(s["SiO2"], alk)
    series = "Alkaline" if is_alkaline(s["SiO2"], alk) else "Subalkaline"
    print(f"  {s['name']:<26} {s['SiO2']:<8.1f} {alk:<8.1f} {tas_name:<25} {series}")

print()
print()
print("Subalkaline Subdivision: Tholeiitic vs Calc-alkaline")
print("=" * 60)
print("Discriminant: FeO*/MgO ratio trends")
print(f"{'SiO2':<10} {'FeO*/MgO (Thol)':<20} {'FeO*/MgO (Calc-alk)'}")
print("-" * 50)

for sio2 in [48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70]:
    # Tholeiitic: FeO*/MgO increases then decreases (iron enrichment)
    x = (sio2 - 48) / 22.0
    thol = 1.0 + 2.5 * x - 1.5 * x**2 if x < 0.8 else 2.0 - 0.5 * x
    # Calc-alkaline: FeO*/MgO stays relatively constant then decreases
    calc = 1.0 + 0.3 * x - 0.8 * x**2
    calc = max(calc, 0.3)
    print(f"  {sio2:<8} {thol:<18.2f} {calc:<18.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Magma Generation: Decompression Melting and Solidus Curves

Python

script.py88 lines

import numpy as np

# Decompression melting at mid-ocean ridges
# Compare adiabat with dry and wet solidus

print("Decompression Melting at Mid-Ocean Ridges")
print("=" * 70)

# Parameters
g = 9.81
rho_m = 3300.0
alpha = 3e-5
cp = 1200.0

# Adiabatic gradient: dT/dz = alpha * g * T / cp
# At T = 1600 K: dT/dz = 3e-5 * 9.81 * 1600 / 1200 = 0.39 C/km

print()
print("Adiabatic gradient: dT/dz = alpha*g*T/cp")
print(f"  = {alpha} * {g} * 1600 / {cp} = {alpha*g*1600/cp:.3f} C/km")
print()

# Depth-pressure conversion: P = rho*g*z, 1 kbar = 100 MPa
# z(km) = P(kbar) * 100e6 / (rho_m * g) / 1000 = P * 30.93 km/kbar

print("Mantle Temperature Profiles vs Solidus Curves")
print("-" * 70)
print(f"{'Depth (km)':<14} {'P (kbar)':<10} {'T_adiab (C)':<14} {'T_dry_sol (C)':<16} {'T_wet_sol (C)':<16} {'Melt?'}")
print("-" * 70)

depths = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 160, 180, 200]

# Potential temperatures to test
Tp_values = [1280, 1350, 1450]  # normal, mid, hot

for Tp in Tp_values:
    print(f"  Potential temperature Tp = {Tp} C:")
    melt_started = False
    for z in depths:
        P = rho_m * g * z * 1000 / 1e8  # kbar
        # Adiabatic temperature
        T_ad = Tp + 0.4 * z  # approximate, 0.4 C/km
        # Dry peridotite solidus (simplified Hirschmann 2000)
        T_dry = 1120 + 3.2 * P + 0.001 * P**2
        # Wet peridotite solidus
        T_wet = 970 + 0.5 * P**2 / (P + 5)  # very simplified

melting = ""
        if T_ad > T_dry:
            melting = "DRY MELT"
            melt_started = True
        elif T_ad > T_wet:
            melting = "wet melt"
            melt_started = True

print(f"  {z:<12} {P:<8.1f} {T_ad:<12.0f} {T_dry:<14.0f} {T_wet:<14.0f} {melting}")

# Find intersection depth
    for z_test in range(200, 0, -1):
        P_t = rho_m * g * z_test * 1000 / 1e8
        T_a = Tp + 0.4 * z_test
        T_d = 1120 + 3.2 * P_t + 0.001 * P_t**2
        if T_a >= T_d:
            print(f"  -> Onset of dry melting at z ~ {z_test} km (P ~ {P_t:.1f} kbar)")
            break
    print()

print()
print("Melt Fraction vs Depth (for Tp = 1350 C)")
print("-" * 50)
print(f"{'Depth (km)':<14} {'F (%)':<10} {'Melt thickness'}")
print("-" * 40)

# Simplified: F increases by ~1%/kbar above solidus
total_melt_thickness = 0
Tp = 1350
for z in range(120, 0, -5):
    P = rho_m * g * z * 1000 / 1e8
    T_ad = Tp + 0.4 * z
    T_dry = 1120 + 3.2 * P + 0.001 * P**2
    if T_ad > T_dry:
        dT = T_ad - T_dry
        F = min(dT / 5.0, 25)  # approx 0.2%/C, max 25%
        total_melt_thickness += F / 100 * 5  # km of melt per 5 km column
        print(f"  {z:<12} {F:<8.1f} {total_melt_thickness:.1f} km (cumulative)")

print(f"  Total melt column: ~{total_melt_thickness:.1f} km -> ~{total_melt_thickness:.0f} km oceanic crust")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Magma Viscosity and Eruption Dynamics

Python

script.py80 lines

import numpy as np

# Magma viscosity as function of composition, temperature, and water content
# Simplified Giordano et al. (2008) model

R = 8.314

print("Magma Viscosity: Effect of Composition and Temperature")
print("=" * 70)

# Approximate Arrhenius parameters for different compositions
compositions = [
    {'name': 'Basalt (dry)', 'A': 1e-5, 'Ea': 200000, 'SiO2': 50},
    {'name': 'Basalt (2% H2O)', 'A': 1e-6, 'Ea': 160000, 'SiO2': 50},
    {'name': 'Andesite (dry)', 'A': 1e-6, 'Ea': 260000, 'SiO2': 58},
    {'name': 'Andesite (2% H2O)', 'A': 1e-7, 'Ea': 220000, 'SiO2': 58},
    {'name': 'Rhyolite (dry)', 'A': 1e-8, 'Ea': 400000, 'SiO2': 75},
    {'name': 'Rhyolite (3% H2O)', 'A': 1e-7, 'Ea': 280000, 'SiO2': 75},
    {'name': 'Rhyolite (6% H2O)', 'A': 1e-6, 'Ea': 200000, 'SiO2': 75},
]

temps = [700, 800, 900, 1000, 1100, 1200, 1300]

print(f"{'Composition':<22} ", end='')
for T in temps:
    print(f"{'T='+str(T)+'C':<12}", end='')
print()
print("-" * 106)

for comp in compositions:
    print(f"  {comp['name']:<20} ", end='')
    for T_C in temps:
        T_K = T_C + 273.15
        eta = comp['A'] * np.exp(comp['Ea'] / (R * T_K))
        if eta < 1:
            print(f"{'<1':>10}  ", end='')
        elif eta < 1e6:
            print(f"{eta:>10.0f}  ", end='')
        else:
            print(f"{eta:>10.1e}  ", end='')
    print(" Pa.s")

print()
print()
print("Conduit Flow: Eruption Rate (Poiseuille Flow)")
print("=" * 65)
print("Q = pi * R^4 * dP / (8 * eta * L)")
print()

# Conduit parameters
L_conduit = 5000.0   # conduit length (m)
dP = 10e6            # overpressure (Pa) = 10 MPa
rho_magma = 2600.0   # kg/m^3

print(f"  Conduit length: {L_conduit/1000:.0f} km, Overpressure: {dP/1e6:.0f} MPa")
print()
print(f"{'Magma type':<20} {'eta (Pa.s)':<15} {'R=1m Q(m^3/s)':<16} {'R=5m Q(m^3/s)':<16} {'R=10m Q(m^3/s)'}")
print("-" * 80)

magma_types = [
    ('Basalt', 100),
    ('Basaltic andesite', 1000),
    ('Andesite', 10000),
    ('Dacite', 1e6),
    ('Rhyolite (wet)', 1e5),
    ('Rhyolite (dry)', 1e10),
]

for name, eta in magma_types:
    Qs = []
    for R_c in [1.0, 5.0, 10.0]:
        Q = np.pi * R_c**4 * dP / (8 * eta * L_conduit)
        Qs.append(Q)
    print(f"  {name:<18} {eta:<13.0e} {Qs[0]:<14.2e} {Qs[1]:<14.2e} {Qs[2]:<10.2e}")

print()
print("Key insight: Eruption rate scales as R^4 and 1/eta")
print("A 5m conduit gives 625x more flux than 1m conduit!")
print("Water reduces rhyolite viscosity by 5 orders of magnitude")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Mineralogy & Crystallography Next: Sedimentary Processes →