# On the Application of 238U, 235U, and 232Th Decay Series Disequilibria to Understanding the Processes and Timescales of Oceanic Basalt Petrogenesis

Kenneth W W Sims, Gregory J Stark, Lynne J Elkins, Mark K Reagan, Peter

B Kelemen, Janne Blichert-Toft

► To cite this version:

Kenneth W W Sims, Gregory J Stark, Lynne J Elkins, Mark K Reagan, Peter B Kelemen, et al.. On the Application of 238U, 235U, and 232Th Decay Series Disequilibria to Understanding the Processes and Timescales of Oceanic Basalt Petrogenesis. Geosciences, 2026, 16 (1), pp.34. 10.3390/geosciences16010034. [https://doi.org/10.3390/geosciences16010034](https://doi.org/10.3390/geosciences16010034). [hal-05457168](https://HAL.inist--CNRS.fr/hal-05457168v1)

HAL Id: hal-05457168

https://hal.science/hal-05457168v1

Submitted on 13 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Review

# On the Application of <sup>238</sup>U, <sup>235</sup>U, and <sup>232</sup>Th Decay Series Disequilibria to Understanding the Processes and Timescales of Oceanic Basalt Petrogenesis

Kenneth W. W. Sims <fnref n="1" />,* ID, Gregory J. Stark <fnref n="1" /> ID, Lynne J. Elkins <fnref n="2" /> ID, Mark K. Reagan <fnref n="3" /> ID, Peter B. Kelemen <fnref n="4" />
and Janne Blichert-Toft <fnref n="5" /> ID

## Abstract

Understanding how processes of magma genesis and magma differentiation control and modify the chemical composition of erupted lavas from the geochemical measurements of the latter is an under-constrained inverse problem, as there is only one known parameter—the measured composition of the erupted lava—but two unknown parameters—the chemical composition and lithology of the source before melting and how melting, crystallization, and melt-rock interactions act to alter the lava en route to the surface. In this invited contribution, we review nearly seven decades of scientific research that demonstrate the potential of U and Th decay series measurements for unraveling the complexities of oceanic magma-tism. We review the underlying nuclear theory, geochemical principles, and application of the <sup>238</sup>U, <sup>235</sup>U, and <sup>232</sup>Th decay series for (i) defining the timescales of magma genesis during decompression mantle melting, (ii) establishing the timescales of magma recharge and magma degassing, and (iii) determining the eruption ages of oceanic Quaternary volcanism.

**Keywords:** U-series disequilibria; basaltic volcanism; high-resolution age dating; young volcanism; mid-ocean ridge basalts; ocean island basalts; geochemistry

<figure><img src="image_2.png" /></figure>

Check for updates

Academic Editor: Changqian Ma

Received: 8 November 2025

Revised: 23 December 2025

Accepted: 24 December 2025

Published: 8 January 2026

Copyright: © 2026 by the authors.

Licensee MDPI, Basel, Switzerland.

This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license.

## 1. Introduction

Understanding basalt petrogenesis from geochemical measurements of erupted lavas is an under-constrained inverse problem because there is one known and two unknown parameters. The known parameter is the measured composition of the erupted lava. The two unknown parameters are (1) the chemical composition and lithology of the source before melting and (2) how processes of magma genesis and magma differentiation (i.e., melting, crystallization, and melt-rock interactions) control and modify the chemical composition of the erupted lava during the transformation from solid source rock to lava.

Fortunately, the U and Th decay series nuclides are ideally suited to the study of magma genesis (Figure 1). The unique efficacy of the <sup>238</sup>U, <sup>235</sup>U, and <sup>232</sup>Th decay series

<footnote n="1">1 Wyoming High Precision Isotope Laboratory, Department of Geology and Geophysics, University of Wyoming, Laramie, WY 82071, USA; greg.stark@uwyo.edu</footnote>

<footnote n="2">2 Department of Earth and Space Science, West Chester University, West Chester, PA 19383, USA; lelkins810@wcupa.edu</footnote>

<footnote n="3">3 Department of Earth and Environmental Sciences, University of Iowa, Iowa City, IA 52242, USA; mark-reagan@uiowa.edu</footnote>

<footnote n="4">4 Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964, USA; peterk@ldeo.columbia.edu</footnote>

<footnote n="5">5 Laboratoire de Géologie de Lyon, Ecole Normale Supérieure de Lyon, CNRS UMR 5276, 46 Allée d'Italie, 69007 Lyon, France; janne.blichert-toft@ens-lyon.fr</footnote>

<footnote>* Correspondence: ksims7@uwyo.edu; Tel.: +1-(307)-766-3306</footnote>
is that radioactive principles constrain the ratios of parent-to-daughter nuclides to be in secular equilibrium before melting, thereby eliminating one of the two typical unknowns in geochemical analysis, namely, the initial chemical composition of the unmelted source in terms of some critical elemental ratios (e.g., U/Th, Th/Ra, Ra/Pb, U/Pa, etc.). Additionally, each decay chain (Figure 2) has multiple parent-daughter ratio pairs whose chemical behaviors are starkly different during magmatic processes, and whose half-lives range from days to hundreds of thousands of years and are thus comparable to the timescales over which magmatic processes occur.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic showing the transformation processes that influence magma compositions and timescales of the magmatic processes that can be discerned by measurement of U and Th decay series disequilibria. In this schematic, Z is the depth of melting, ΔZ is the length of the melt column, W₀ is the velocity of the upwelling (solid) mantle, and ω is the melt velocity. The timescales τ₁, τ₂, τ₃, τ₄, and τ₅ are defined, respectively, as (1) the timescale of solid mantle upwelling; (2) the timescale of melt migration; (3) the delay time, or storage time, of magma before eruption; (4) the timescale of degassing; and (5) the time since the lava of interest erupted. Black arrows signify solid mantle upwelling and red arrows signify melt migration.</figcaption></figure>

The application of U and Th decay series disequilibria to studies of magma genesis was thoroughly reviewed in the Reviews in Mineralogy and Geochemistry volume on U-series in 2003 [1]. Over the past two decades, however, there have been enough significant advances in our understanding of this topic to warrant an updated review. The goal of the present review thus is to bring the community up to the current state-of-the-art on U and Th decay series disequilibria with specific focus on ²³⁸U, ²³⁵U, and ²³²Th decay series systematics as a chronometer of igneous processes (Figure 1). To make the overview as complete and clear as possible, the sections to follow this Introductory section have been organized as follows:

* Section 2 presents an overview of the theoretical formalism of the U and Th decay series' unique isotope systematics, using ²³⁸U-²³⁴U-²³⁰Th within the ²³⁸U decay series as an explicit example.

* Section 3 examines the geochemical principles and assumptions particular to the application of the ²³⁸U, ²³⁵U, and ²³²Th decay series to decompression mantle melting in oceanic settings. It focuses on the differences between time-independent melting mod-
els such as batch and tractional melting, which rely solely on elemental fractionation, and time-dependent melting models that in turn take into consideration the timescales of melting and incorporate daughter ingrowth during melting. Specific examples of how different melting models and different lithologies produce divergent results, as well as the reasoning behind interpreting and selecting any model over another, are then considered. Melt generation at both divergent plate boundaries (MORB) and intraplate hotspots (OIB) are assessed.

* Section 4 presents the application of the ²³⁸U, ²³⁵U, and ²³²Th decay series for dating eruption ages of volcanic rocks using three different approaches: (1) age limits, whereby the presence of disequilibria can place bounds on the eruption age of a lava; (2) internal mineral isochron methods, which establish the crystallization age of a volcanic rock; and (3) U-Th-Ra 'model' ages, which couple measurements of U and Th disequilibria with geological, geochemical, and isotopic constraints to establish eruption ages for a lava or a lava sequence.

* Section 5 presents the application of ²²⁶Ra-²¹⁰Pb-²¹⁰Po disequilibria to establishing the timescales of magma recharge and degassing, which is directly linked to a volcano's activity and fundamental to our understanding of physical eruption dynamics and hazard assessment. We look at three specific examples: Kilauea in the Hawaiian Archipelago, Tristan de Cunha in the South Atlantic Ocean, and Mount Erebus Volcano in Antarctica.

* Section 6 provides a historical overview of current and past analytical methodologies used to measure ²³⁸U, ²³⁵U, and ²³²Th decay series isotopes. We highlight recent improvements in mass spectrometry that have greatly influenced this field of study, as well as the alpha counting methods used for ²¹⁰Pb-²¹⁰Po measurements.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Schematic of the ²³⁸U, ²³⁵U, and ²³²Th decay series. Boxes for each radioisotope include approximate half-lives and type of decay (α and/or β⁻). Green, blue, and orange denote the isotopes within the ²³⁸U, ²³⁵U, and ²³²Th decay series, respectively. The table includes common parent-daughter pairs utilized in studying igneous processes and their associated measurement techniques. Modified from [2].</figcaption></figure>
## 2. The Theoretical Construct of Multiple Nuclides Decaying in a Series

The U and Th decay series chains (Figure 2) comprise three naturally occurring decay series: the ²³⁸U series (referred to as the uranium series, but historically referred to as the radium series), the ²³⁵U series (referred to as the actinium series), and the ²³²Th series (referred to as the thorium series). There is also a ²³⁷Np series that is no longer extant because of the relatively short half-life of parental ²³⁷Np (2.2 × 10⁶ years) relative to the age of the elements in the solar system and the age of Earth. However, this chain is regularly produced artificially in the nuclear industry, and ²³³Pa, which is the progeny of ²³⁷Np, is often used as a standard ('spike') for ²³¹Pa isotope dilution measurements.

Bateman [3] first formalized the symmetrical form of multiple nuclides decaying in a series that is the central premise of the ²³⁸U, ²³⁵U, and ²³²Th decay series. Several review texts have presented the theoretical basis of the decay series (e.g., [4]). Here, we present the simplest case in which two nuclides in a decay series, a parent-daughter pair, are linked by principles of radioactive decay and ingrowth (Figures 3 and 4).

<figure><img src="image_5.png" /><figcaption>**Figure 3.** A portion of the ²³⁸U decay series, extending to ²²²Rn. Sphere sizes represent the relative natural abundances (i.e., larger spheres represent more atoms; actual sizes would vary by many orders of magnitude). The arrow widths correspond to the magnitude of the decay constant (λ). The equations reflect radioactive secular equilibrium within the decay series, in which the activities (defined as A = λN) of parent and daughter nuclides are equal. λ = decay constant; N = number of atoms; A = activity. Modified from [2].</figcaption></figure>

Consider, for example, the ²³⁸U series: ²³⁸U → ²³⁴U → ²³⁰Th → ²²⁶Ra → ... → ²⁰⁶Pb. As will be shown in the following derivation, radioactive secular equilibrium within a decay series is defined as λ₁N₁ = λ₂N₂ = λ₃N₃ or (A₁) = (A₂) = (A₃), where λᵢ is the decay constant of isotope i; Nᵢ is the number of atoms of isotope i; and (Aᵢ) is the activity of isotope i (parenthesis around an isotope denote activity) and is defined as:

$$ (A) = - \frac{dN}{dt} = \lambda N \qquad (1) $$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Attainment of radioactive secular equilibrium after an initial fractionation between hypothetical long-lived parent (N₁; t₁/₂ = 10 h) and short-lived daughter (N₂; t₁/₂ = 1 h) isotopes. Two scenarios are shown: zero initial daughter isotope (A–C) and initial daughter excess (D–F). (A) In the zero initial daughter scenario, the number of moles of the daughter isotope increases until it is in secular equilibrium with the parent isotope, which occurs after about five half-lives of the daughter. After secular equilibrium is reached, the ratio of parent to daughter is a constant value of (λ₂ – λ₁)/λ₁. (B) In terms of activity, the daughter isotope activity increases from zero until it comes into secular equilibrium with the parent, after which the activities of the parent and daughter isotope remain equal. (C) A plot showing the activity ratio (parentheses denote activity) of N₂/N₁. The activity ratio increases from zero until secular equilibrium is reached. The trend in the number of moles and, hence, the activity of N₂ is non-linear due to the decay of N₂. The daughter excess scenario in (D–F) exhibits an inverse relationship to the zero initial daughter scenario. (D) The initial number of moles of the daughter isotope decreases until secular equilibrium with the parent isotope is reached, after which the ratio of parent to daughter is again a constant value of (λ₂ – λ₁)/λ₁. (E) The activity of the daughter isotope decreases from its initial excess until it comes into secular equilibrium with the parent, after which the activities of parent and daughter are equal. (F) The activity ratio of N₂/N₁ decreases from an initial daughter excess until secular equilibrium is reached after about five half-lives of the daughter isotope. Modified from [2].</figcaption></figure>

Activity typically has units of decays per second per gram (dps/g or Becquerels/g), decays per minute per gram (dpm/g), or decays per year per gram (dpy/g) depending on which time unit (seconds, minutes, years) is used for the decay constant. So, in a decay series, where N₁ → N₂ → N₃:

$$ -\frac{dN_1}{dt} = \lambda_1 N_1 \qquad (2) $$

$$ \frac{dN_2}{dt} = \lambda_1 N_1 - \lambda_2 N_2 \qquad (3) $$
Since N₁ undergoes radioactive decay, its abundance at time t is determined by:

$$N_1 = N_1^0 e^{-\lambda_1 t} \quad (4)$$

where the term N₁⁰ is the number of atoms at t = 0. Substituting Equation (4) into Equation (3) yields:

$$\frac{dN_2}{dt} = \lambda_1 N_1^0 e^{-\lambda_1 t} - \lambda_2 N_2 \quad (5)$$

and

$$\frac{dN_2}{dt} + \lambda_2 N_2 - \lambda_1 N_1^0 e^{-\lambda_1 t} = \quad (6)$$

Solving for N₂:

$$N_2 = \frac{\lambda_1}{\lambda_2 - \lambda_1} N_1^0 (e^{-\lambda_1 t} - e^{-\lambda_2 t}) + N_2^0 e^{-\lambda_2 t} \quad (7)$$

The first term gives the number of atoms of N₂ that formed from decay of N₁ but have not yet decayed to N₃. The second term gives the number of atoms N₂ that remain from initial N₂⁰. If there are no atoms of the original daughter N₂⁰, then Equation (7) reduces to:

$$N_2 = \frac{\lambda_1}{\lambda_2 - \lambda_1} N_1^0 (e^{-\lambda_1 t} - e^{-\lambda_2 t}) \quad (8)$$

The attainment of secular equilibrium amongst a parent-daughter pair is a function of half-life and can be estimated by exploring the limits of Equation (8), to be in place after six or more half-lives have passed. After six half-lives of daughter N₂, e⁻⁷t² → , so:

$$N_2 \approx \frac{\lambda_1}{\lambda_2 - \lambda_1} N_1^0 (e^{-\lambda_1 t}) \approx \frac{\lambda_1}{\lambda_2 - \lambda_1} N_1 \quad (9)$$

If the parent has a significantly longer half-life than the daughter, then:

$$\lambda_2 - \lambda_1 \approx \lambda_2, N_2 \approx \frac{\lambda_1}{\lambda_2} N_1, \text{ and } N_2\lambda_2 \approx N_1\lambda_1 \quad (10)$$

The attainment of secular equilibrium after an initial parent-daughter fractionation is graphically illustrated in Figure 4. In this hypothetical parent-daughter system, the parent isotope (N₁; t₁/₂ = 10 h) has a half-life 10 times longer than its shorter-lived daughter (N₂; t₁/₂ = 1 h). Two scenarios are illustrated. In the first example, the daughter isotope has been entirely removed from the parent (Figure 4A-C), so N₂⁰ = 0. For this situation, the number of daughter atoms (N₂) increases over time until equilibrium is established, after which the ratio of parent to daughter is a constant value of (λ₂ - λ₁)/λ₁, which, as shown above, is reasonably approximated as λ₂/λ₁ since λ₂ - λ₁ ≈ λ₂, or in this case, N₁/N₂ = 10. In the second example, the activity of the daughter is enriched by a factor of 300% over the parent nuclide (Figure 4D-F). Although the activity of N₂ initially is greater than that of N₁ by a factor of three, then because of the factor of ten difference in half-lives and the definition of activity (Nλ), the number of atoms of N₂ is still less than the number of atoms of N₁ by a factor of five. In this scenario, the activity of the daughter atoms (N₂) decreases over time until equilibrium and (λ₂ - λ₁)/λ₁.

Now consider the actual ²³⁸U-²³⁰Th system ²³⁸U→₉₂U→ ²³⁴U→ ²³⁰Th→ ²²⁶Ra→ ... → ²⁰⁶Pb. The ²³⁸U-²³⁴U-²³⁰Th example is ideal for illustration, as (i) it is the most commonly used decay series chronometer, and, in general terms, (ii) these parent-daughter pairs are commonly used for most geochronological applications.

$$(^{230}\text{Th})_A = (^{230}\text{Th})_{AX} + (^{230}\text{Th})_{AS} \quad (11)$$
where $(^{230}\text{Th})_A$ is the total activity, $(^{230}\text{Th})_{AX}$ is the unsupported excess activity ($\lambda N$, where parentheses denote activity), and $(^{230}\text{Th})_{AS}$ is the $^{230}\text{Th}$ supported by the decay of $^{238}\text{U}$. In a system in which geochemical fractionation of Th from U ceases, the unsupported excess, $(^{230}\text{Th})_{AX}$, decays according to simple radioactive principles, such that at time t:

$$(^{230}\text{Th})_{AX} = (^{230}\text{Th})_{AX}^{0} e^{-\lambda_{230}t} \qquad (12)$$

where $(^{230}\text{Th})_{AX}^0$ is the initial unsupported excess. Normalized to $(^{232}\text{Th})$, the total $(^{230}\text{Th}/^{232}\text{Th})$ excess is:

$$\left(\frac{^{230}\text{Th}}{^{232}\text{Th}}\right)_{AX} = \left(\frac{^{230}\text{Th}}{^{232}\text{Th}}\right)_{AX}^{0} e^{-\lambda_{230}t} \qquad (13)$$

Now let us look at the second component in the total activity, which is the supported $(^{230}\text{Th})$ from the decay of $^{238(4)}\text{U}$. From Equation (8):

$$N_2 = \frac{\lambda_1}{\lambda_2 - \lambda_1} N_1^0 (e^{-\lambda_1 t} - e^{-\lambda_2 t})$$

In terms of $^{230}\text{Th}$:

$$^{230}\text{Th}_{AS} = \frac{\lambda_{234}}{\lambda_{230} - \lambda_{234}} {}^{234}\text{U}^0 (e^{-\lambda_{234}t} - e^{-\lambda_{230}t}) \qquad (14)$$

Given that $^{238}\text{U} \rightarrow {}^{234}\text{U} \rightarrow {}^{230}\text{Th}$, if $^{234}\text{U}$ is in secular equilibrium with $^{238}\text{U}$, then $(^{234}\text{U}) = (^{238}\text{U})$ and thus $^{234}\text{U}^0\lambda_{234} = (^{238}\text{U})_A$, $\lambda_{230} - \lambda_{234} \approx \lambda_{230}$, and $e^{-\lambda_{234}t} \approx 1$. From this follows that:

$$\lambda_{230} {}^{230}\text{Th}_{AS} = \lambda_{234} {}^{234}\text{U}^0 (1 - e^{-\lambda_{230}t}) \qquad (15)$$

$$(^{230}\text{Th})_{AS} = (^{238}\text{U})(1 - e^{-\lambda_{230}t}) \qquad (16)$$

From Equation (12):

$$(^{230}\text{Th})_{AX} = (^{230}\text{Th})_{AX}^{0} e^{-\lambda_{230}t}$$

and Equation (11):

$$(^{230}\text{Th})_A = (^{230}\text{Th})_{AX} + (^{230}\text{Th})_{AS}$$

so that:

$$(^{230}\text{Th})_A = (^{230}\text{Th})_{AX} e^{-\lambda_{230}t} + (^{238}\text{U})(1 - e^{-\lambda_{230}t}) \qquad (17)$$

When normalized to $^{232}\text{Th}$:

$$\left(\frac{^{230}\text{Th}}{^{232}\text{Th}}\right)_A = \left(\frac{^{230}\text{Th}}{^{232}\text{Th}}\right)_{AX} e^{-\lambda_{230}t} + \left(\frac{^{238}\text{U}}{^{232}\text{Th}}\right) (1 - e^{-\lambda_{230}t}) \qquad (18)$$

This is a linear equation in the space of $(^{230}\text{Th}/^{232}\text{Th})$ versus $(^{238}\text{U}/^{232}\text{Th})$ and is the underlying equation used for the $^{238}\text{U}-^{230}\text{Th}$ chronometer [5]. Similar variants exist for the other U-series chronometers. It should be noted that $^{234}\text{Th}$ and $^{234}\text{Pa}$, while being important tracers in chemical oceanography [6], are essentially disregarded for dating purposes, as they are assumed to be in secular equilibrium due to their very short half-lives $(^{234}\text{Th} t_{1/2} = 24.1 \text{ d}, ^{234}\text{Pa} t_{1/2} = 1.18 \text{ m})$.
### 3. U-Th-Ra and U-Pa Series Constraints on Decompression Melting at Divergent Plate Boundaries and Intraplate Hotspots

Beneath both mid-ocean ridge spreading centers and intraplate hotspot settings, the dominant mode of melt production is adiabatic decompression of solid upwelling mantle material. This melting process induces chemical fraction between parent-daughter nuclides that results in U-series disequilibria. The timescale over which this U-series disequilibrium can be used to understand the timescale of melting is highly dependent on the half-life of the daughter nuclide in the parent-daughter pair (Figure 2; Table 1). In the following, we examine the U-series models that have been proposed to constrain the timescales of adiabatic decompression melting of solid upwelling mantle.

<table><caption>**Table 1.** Summary of important U-series radionuclides.</caption><thead><tr><th>Isotope</th><th>Abundance</th><th>Half-Life</th><th>Decay Constant</th><th>Reference</th></tr></thead><tbody><tr><td><sup>238</sup>U</td><td>99.2742%</td><td>4.4683 ± 0.0024 × 109 yr</td><td>1.5513 × 10<sup>-10</sup> yr<sup>-1</sup></td><td>[7]</td></tr><tr><td><sup>234</sup>U</td><td>0.0054%</td><td>245,620 ± 260 yr</td><td>2.82206 × 10<sup>-6</sup> yr<sup>-1</sup></td><td>[8]</td></tr><tr><td><sup>230</sup>Th</td><td>0.02%</td><td>75,584 ± 110 yr</td><td>9.1706 × 10<sup>-6</sup> yr<sup>-1</sup></td><td>[8]</td></tr><tr><td><sup>226</sup>Ra</td><td>Trace</td><td>1600 ± 7 yr</td><td>4.332 × 10<sup>-4</sup> yr<sup>-1</sup></td><td>[9]</td></tr><tr><td><sup>222</sup>Rn</td><td>Trace</td><td>3.8235 ± 0.0004 d</td><td>1.8 × 10<sup>-1</sup> d<sup>-1</sup></td><td>[9]</td></tr><tr><td><sup>210</sup>Pb</td><td>Trace</td><td>22.6 ± 0.1 yr</td><td>3.07 × 10<sup>-2</sup> yr<sup>-1</sup></td><td>[10]</td></tr><tr><td><sup>210</sup>Po</td><td>Trace</td><td>138.4 ± 0.1 d</td><td>5.008 × 10<sup>-3</sup> d<sup>-1</sup></td><td>[10]</td></tr><tr><td><sup>235</sup>U</td><td>0.7204%</td><td>7.0381 ± 0.0048 × 108 yr</td><td>9.8485 × 10<sup>-10</sup> yr<sup>-1</sup></td><td>[7]</td></tr><tr><td><sup>231</sup>Pa</td><td>100%</td><td>32,760 ± 220 yr</td><td>2.1158 × 10<sup>-5</sup> yr<sup>-1</sup></td><td>[11]</td></tr><tr><td><sup>227</sup>Ac</td><td>Trace</td><td>21.77 ± 0.02 yr</td><td>3.184 × 10<sup>-2</sup> yr<sup>-1</sup></td><td>[10]</td></tr><tr><td><sup>232</sup>Th</td><td>99.98%</td><td>1.401 ± 0.007 × 1010 yr</td><td>4.948 × 10<sup>-11</sup> yr<sup>-1</sup></td><td>[12]</td></tr><tr><td><sup>228</sup>Ra</td><td>Trace</td><td>5.75 ± 0.03 yr</td><td>1.21 × 10<sup>-1</sup> yr<sup>-1</sup></td><td>[9]</td></tr><tr><td><sup>224</sup>Ra</td><td>Trace</td><td>3.6319 ± 0.0023 d</td><td>1.9085 × 10<sup>-1</sup> d<sup>-1</sup></td><td>[10]</td></tr></tbody></table>

Two end-member processes (Figures 5 and 6) have been proposed to account for the <sup>238</sup>U-<sup>230</sup>Th disequilibria measured in oceanic basalts: (1) time-independent melting, in which net elemental fractionation of Th from U at low melt fractions (and porosities) is the driving mechanism creating disequilibria [13–15], and (2) time-dependent melting, where <sup>230</sup>Th grows in over the timescale of melt generation [16–19]. If net elemental fractionation of U/Th during partial melting is the dominant effect, then the U/Th ratio of the mantle source is best determined from the <sup>230</sup>Th/<sup>232</sup>Th activity ratio of the lavas. On the other hand, if <sup>230</sup>Th ingrowth is the fundamental mechanism generating <sup>238</sup>U-<sup>230</sup>Th disequilibria, then the U/Th source ratio is best determined from the <sup>238</sup>U/<sup>232</sup>Th activity ratio of the lavas.

In time-dependent adiabatic decompression melting models (Figures 5 and 6), in-growth of <sup>230</sup>Th produces elevated (<sup>230</sup>Th/<sup>232</sup>Th) leading to higher (<sup>230</sup>Th/<sup>238</sup>U). Therefore, when upwelling is fast (e.g., beneath Kilauea volcano, which is the central locus of the Hawaiian plume), ingrowth is minimal and net elemental fractionation of Th/U during partial melting is the dominant effect [14,15], whereas when upwelling is slow (e.g., at mid-ocean ridges), <sup>230</sup>Th ingrowth is the fundamental mechanism producing U-Th disequilibria.

#### 3.1. *Time-Dependent Melting Models*

Several models relating the isotopic effects of nuclide ingrowth in U-series disequilibria to the timescales of the melting process have been proposed (e.g., [14,17–28]; however, credit for the time-dependent models (Figures 5 and 6) should be given to (a) McKenzie [20] for dynamic melting and (b) Spiegelman and Elliott [18] for reactive porous flow melting, also called chromatographic melting.
Magma Genesis: Two end-member models

Time-independent melting

<figure><img src="image_7.png" /></figure>

$$ (1) \ (^{230}\text{Th}/^{232}\text{Th})_L = (^{238}\text{U}/^{232}\text{Th})_S $$

$$ (2) \ (^{230}\text{Th}/^{232}\text{Th})_L \text{ is better correlated with } ^{87}\text{Sr}/^{86}\text{Sr} $$

Time-dependent melting

<figure><img src="image_8.png" /></figure>

$$ (1) \ (^{238}\text{U}/^{232}\text{Th})_L = (^{238}\text{U}/^{232}\text{Th})_S $$

$$ (2) \ (^{238}\text{U}/^{232}\text{Th})_L \text{ is better correlated with } ^{87}\text{Sr}/^{86}\text{Sr} $$

**Figure 5.** Schematic showing the systematics of ²³⁸U-²³⁰Th in the two end-member models proposed to account for the ²³⁸U-²³⁰Th disequilibria measured in oceanic basalts: (1) Time-independent melting in which net elemental fractionation of Th from U at low melt fractions (and porosities) is the driving mechanism creating disequilibria [13,14,18] and (2) time-dependent melting, where ²³⁰Th grows in over the timescale of melt generation [16,18]. If net elemental fractionation of U/Th during partial melting is the dominant effect, then the U/Th ratio of the mantle source is best determined from the ²³⁰Th/²³²Th activity ratio of the lavas. On the other hand, if ²³⁰Th ingrowth is the fundamental mechanism generating ²³⁸U-²³⁰Th disequilibria, then the U/Th source ratio is best determined from the ²³⁸U/²³²Th activity ratio of the lavas, and ²³⁸U/²³²Th will be best correlated with long-lived radiogenic isotope ratios, such as ⁸⁷Sr/⁸⁶Sr.

While these two end-member models, (a) dynamic versus (b) chromatographic melting, differ mainly in their treatment of the melt extraction process (i.e., (a) equilibrium in the trapped pore space plus fractional melting for the extracted portion, versus (b) porous flow of a migrating liquid that continuously reacts and equilibrates with the surrounding solid matrix), they both incorporate the effect of radioactive ingrowth and decay of the daughter isotopes (e.g., ²³⁰Th, ²²⁶Ra, ²³¹Pa) during melt production. As such, these models predict mantle melting and melt transport rates based on U-Th-Ra and U-Pa disequilibria. These U-series mantle melting rates can, in turn, be combined with theoretical predictions of melt productivity (e.g., [29,30]) to calculate solid mantle upwelling rates.

In the following, we briefly review the two different end-member melting models incorporating the timescales of the melting process and compare these models with a global MORB database (Figures 7–11). We then show, through examples, how the U-series data measured in lavas can be compared with these models to place unique temporal constraints on many of the magmatic processes associated with adiabatic decompression and time-dependent melting.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Schematic illustrating the differences in the dynamic melting modeling of McKenzie [20] and the chromatographic melting model of Spiegelman and Elliott [18]. (A) In dynamic melting, melts (red spheres) produced in the upwelling mantle are retained in pore spaces and remain in chemical equilibrium with the solid matrix (green spheres) until a critical threshold porosity (φ) is reached. Any melts produced in excess of this critical threshold value will be instantaneously extracted into a melt channel, where they quickly move to the surface. As a result, melts do not react with the solid matrix as they travel upwards. In contrast to chromatographic melting, disequilibria are primarily created at the base of the melt column where parental uranium is extracted and is a function of the timescale of melting. (B) In chromatographic melting, chemical equilibrium is maintained through continuous interaction between melt and solid during melt transport. Different elements travel at different velocities through the melt column according to their mineral/melt partition coefficients. In this scenario, the daughter is less compatible than the parent, so it travels faster through the melt column. This difference in melt column residence times is represented by the vertical distance between parent and daughter decays (denoted by short curved arrows). Because the daughter is extracted at the surface before it decays, the daughter becomes enriched relative to the parent, which has a longer residence time in the melt column. Parent and daughter activities are denoted by C<sub>p</sub> and C<sub>d</sub>, respectively. In the subsolidus mantle, parent and daughter travel at the same velocity (the mantle upwelling rate) and maintain secular equilibrium.</figcaption></figure>

## 3.2. *Dynamic Melting*

Solutions for <sup>238</sup>U-<sup>230</sup>Th-<sup>236</sup>Ra disequilibria during dynamic melting were first derived by McKenzie [20] and have subsequently been presented in several papers pertaining to the production of U-series disequilibria in basaltic melts [17,24,25,27,28,31]. By assuming a constant melt rate (Γ, mass fraction per time) and a constant porosity (φ, volume fraction of melt), the (<sup>230</sup>Th/<sup>238</sup>U) and (<sup>226</sup>Ra/<sup>230</sup>Th) can be solved analytically and expressed as a function of φ and Γ (see McKenzie [20] or Williams and Gill [17] for the full derivation):

$$ \frac{^{230}\text{Th}}{^{238}\text{U}} = \frac{F_{\text{Th}}(K_{\text{U}} + \lambda_{\text{Th}})}{F_{\text{U}}(K_{\text{Th}} + \lambda_{\text{Th}})} \qquad (19) $$

and

$$ \frac{^{226}\text{Ra}}{^{230}\text{Th}} = \left[ \frac{F_{\text{Th}}(K_{\text{U}} + \lambda_{\text{Ra}})}{F_{\text{U}}(K_{\text{Th}} + \lambda_{\text{Ra}})} \right] \left[ 1 + \frac{\lambda_{\text{Th}}(K_{\text{U}} - K_{\text{Th}})}{(K_{\text{Th}} + \lambda_{\text{Ra}})(K_{\text{U}} + \lambda_{\text{Th}})} \right] \qquad (20) $$
where

$$F_i = \frac{\phi \rho_f}{D_i \rho_s (1 - \phi) + \phi \rho_f} \qquad (21)$$

and

$$K_i = \frac{F_{Th}(1 - D_i)}{\phi \rho_f} \Gamma_0 \qquad (22)$$

where Dᵢ is the bulk mineral/melt partition coefficient, and ρf and ρs are the densities of the melt and peridotite source (2800 and 3300 kg·m⁻³), respectively. Note that an expression similar to Equation (19) can be derived for (²³¹Pa/²³⁵U).

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Global MORB data set for which U-series nuclides have been measured [14,19,25,32–53]. (A) (²³⁰Th/²³²Th) vs. (²³⁸U/²³²Th); (B) (²²⁶Ra/²³⁰Th) vs. (²³⁰Th/²³⁸U); and (C) (²³¹Pa/²³⁵U) vs. (²³⁰Th/²³⁸U) for the global data set. The solid line in (A) and the dashed lines in (B,C) are equilines, where the activities of the parent and daughter isotopes are equal.</figcaption></figure>

In the dynamic melting model, as formulated by McKenzie [20], melts produced in the upwelling mantle are retained in the pore space and are in equilibrium with the solid residue until a critical threshold porosity is achieved, after which any produced melt in excess of this critical threshold value is extracted instantaneously such that the porosity remains constant at the threshold value (Figure 6A). Dynamic melting takes into consideration the timescale of the melting process, and hence, ingrowth produces enhanced excesses of short-lived daughter nuclides such as ²²⁶Ra and to a lesser extent ²³⁰Th (Figure 8). In dynamic melting, the produced melts are also assumed to move rapidly to the surface as they form, without reacting with surrounding rocks during transport. Because most of the parental uranium is extracted near the base of the melting column, disequilibria are created primarily at the bottom of the melt column instead of throughout its entire length. (This contrasts with the outcomes of chromatographic melting presented in the next section).
With dynamic melting, melt velocities can thus be constrained by the shortest half-life of the daughter nuclide in radioactive disequilibrium with its immediate parent.

## Dynamic Melting

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Results of dynamic melting calculations for an adiabatically decompressing peridotitic mantle showing (226Ra/230Th) and (231Pa/235U) versus (230Th/238U) gridded melting model results for the solid mantle upwelling rates (W) and residual threshold porosity values (f) listed. Model results are shown for (A) the analytical solution to dynamic melting; (B) the numerical solution assuming constant bulk solid/liquid partitioning and melting rates; (C) the full numerical solution with variable elemental partitioning and melting rates, all with T<sub>p</sub> = 1300 °C; (D) the analytical solution to dynamic melting; (E) the numerical solution assuming constant bulk solid/liquid partitioning and melting rates; and (F) the full numerical solution with variable elemental partitioning and melting rates, all with T<sub>p</sub> = 1400 °C. Global MORB data are shown for comparison (see the caption in Figure 7 for the data sources).</figcaption></figure>

Once the critical threshold porosity is reached, the porosity of the melt zone is constant (Figures 6 and 8). Trace element enrichments in the melts are intermediate between batch melting and accumulated fractional melting, depending on the value of the threshold porosity relative to the distribution coefficients of the element(s) considered [17,31]. Thus, the retained melt porosity plays an important role in slowing the extraction of an incompatible element from the solid. With a small critical porosity, dynamic melting is closer to fractional melting, in which incompatible elements are efficiently stripped from the residue, whereas with a large critical porosity, dynamic melting is more akin to batch melting, in which incompatible elements reside in chemical equilibrium with the solid.

### 3.3. Chromatographic Porous Flow Melting

The steady-state chromatographic porous flow melting model of Spiegelman and Elliot [18], in its main form, assumes that chemical equilibrium is maintained between migrating liquid and the solid matrix. The essential characteristic of chromatographic melting is that continuous melt-solid interaction between the ascending melt and mantle source results in different elements traveling through the melting column at different velocities according to their relative melt/solid partition coefficients (i.e., D<sub>i</sub> values) (Figure 6B).
For non-radioactive elements, the equilibrium porous flow melting model produces trace-element enrichments that are identical to those of equilibrium batch melting [18]. However, for the U-series isotopes that have half-lives that are comparable to melt migration timescales, chromatographic effects can have a significant influence on resulting isotopic concentrations and activity ratios in the partial melt. If the daughter isotope is more incompatible than the parent isotope, which is largely true for the typical mantle lithologies considered here (e.g., [54] and references therein), the residence time of the daughter isotope in the melting column is shorter relative to the parent nuclide, as the more incompatible daughter nuclide travels preferentially with the melt while the more compatible element travels more slowly. This residence time effect enhances the production of daughter nuclides relative to their rates of decay in the melting column, which produces relatively large excesses in daughter isotopes in the resultant liquids. The extent of this daughter nuclide enhancement depends on the half-life of the daughter nuclide, and, consequently, is relatively large for $^{226}$Ra and much smaller for $^{230}$Th and $^{231}$Pa. The extent of ingrowth is also dependent on the length of the melting column, such that longer melt columns generate greater enhancements of the daughter nuclides due to longer nuclide residence times.

The one-dimensional (1D) chromatographic porous flow model of Spiegelman and Elliott [18] has seven variable parameters (which are not all independent): the solid 1D velocity (W), the liquid 1D velocity (ω), the porosity (φ), the height of the melting column (Δz), the melt fraction (F), the melting rate (Γ), and the bulk mineral/melt partition coefficients (Dᵢ) for the trace elements being considered. In the fully expanded form of the model, most of these variables may vary with depth and, therefore, with time, though the full melt column height is fixed, and the solid velocity is typically treated as an externally imposed constant controlled by regional mantle convection and migration rates (W₀). The main parameters utilized in melting models, along with their typical ranges, are summarized in Table 2.

<table><caption>**Table 2.** Main parameters used in melting models.</caption><thead><tr><th>Parameters</th><th>Notation</th><th>Units</th><th>Range</th></tr></thead><tbody><tr><td>Degree of melting</td><td>F</td><td>-</td><td>0–0.2</td></tr><tr><td>Solid upwelling rate (1D)</td><td>W</td><td>m·yr<sup>-1</sup></td><td>0–1</td></tr><tr><td>Melting rate</td><td>Γ</td><td>kg·m<sup>-3</sup>·yr<sup>-1</sup></td><td>10<sup>-3</sup> to 10<sup>-5</sup></td></tr><tr><td>Solidus temperature</td><td>T</td><td>K</td><td>1000–1400</td></tr><tr><td>Length of melting column</td><td>Z</td><td>km</td><td>0–150</td></tr><tr><td>Plate velocity</td><td>-</td><td>m·yr<sup>-1</sup></td><td>0.01–0.1</td></tr><tr><td>Matrix porosity</td><td>φ</td><td>-</td><td>0–0.05</td></tr><tr><td>Solid density</td><td>ρ<sub>s</sub></td><td>kg·m<sup>-3</sup></td><td>3300</td></tr><tr><td>Melt density</td><td>ρ<sub>f</sub> or ρ<sub>m</sub></td><td>kg·m<sup>-3</sup></td><td>2800</td></tr><tr><td>Melt velocity (1D)</td><td>ω</td><td>m·yr<sup>-1</sup></td><td>-</td></tr><tr><td>Lithosphere thickness</td><td>h</td><td>km</td><td>0–100</td></tr><tr><td>Decay Constants</td><td>λ</td><td>yr<sup>-1</sup></td><td>>10<sup>-1</sup> to 10<sup>-11</sup></td></tr><tr><td>Bulk Partition Coefficient</td><td>D<sub>i</sub></td><td>-</td><td>0.00001 to 0.1</td></tr><tr><td>Buoyancy flux</td><td>B</td><td>Mg·s<sup>-1</sup></td><td>0.5–9</td></tr></tbody></table>

Based on the formulations of Spiegelman and Elliott [18], Spiegelman [55] developed a numerical solution for $^{238}$U-$^{230}$Th-$^{226}$Ra disequilibria generated during chromatographic equilibrium porous flow melting. While Spiegelman and Elliott [18] pointed out that the system of equations to produce U-series disequilibria during chromatographic melting could be simplified sufficiently to be solved analytically, a set of analytically solvable expressions was not developed until Sims et al. [25]. (As discussed in Sims et al. [25], however, determining this solution required correcting an error in Equation (15) from
Spiegelman and Elliott [18].) We present the analytical approximation solution below for completeness and then compare outcomes of the analytical approach to those of full numerical solutions (Figure 9).

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Results of chromatographic/equilibrium porous flow melting calculations for an adiabatically decompressing peridotitic mantle showing (A,C) ($^{226}$Ra/$^{230}$Th) and (B,D) ($^{231}$Pa/$^{235}$U) versus ($^{230}$Th/$^{238}$U) gridded melting model results for the solid mantle upwelling rates (W) and residual threshold porosity values (f) listed, all at T$_{p}$ = 1400 °C. (A,B) show results for the analytical solution after Sims et al. [25], assuming constant bulk solid/melt partitioning and melting rates, while (C,D) show outcomes for variable mineralogy and bulk partitioning and variable melting rates, as described in the text. Global MORB data are shown for comparison (see the caption in Figure 7 for the data sources).</figcaption></figure>

As originally derived in Spiegelman and Elliott [18], the effects of transport of melt through a melting column on the chemistry of radiogenic isotopes can be separated from melting by expressing the concentration of an element as:

$$c_i = \alpha_i c_{bi} \qquad (23)$$

where $c_i$ is the concentration of an element measured at the surface, $c_{bi}$ is the concentration of a stable element due to batch melting, and $\alpha_i$ is the enrichment factor due to transport. For equilibrium porous flow, because re-equilibration of solid and liquid is assumed to occur continuously, stable element concentrations in the melt are determined by the batch
melting equation. Spiegelman and Elliott [18] developed an alternative approach that assumes true fractional melting, i.e., where the solid does not re-equilibrate with the melt that passes through it (see Appendix of Spiegelman and Elliott [18]). This approach, along with additional methods, have recently been implemented by Elkins and Spiegelman [56] and will be addressed in further detail below.

For the decay chain $^{238}$U-$^{230}$Th-$^{226}$Ra, the change in $\alpha_i$ with position in a 1D melt column for each isotope can be expressed by the following equations:

$$ \frac{d\alpha_0}{d\zeta} = -\lambda_0\tau_0\alpha_0 \qquad (24) $$

$$ \frac{d\alpha_1}{d\zeta} = \lambda_1(\tau_0\alpha_0 - \tau_1\alpha_1) \qquad (25) $$

$$ \frac{d\alpha_2}{d\zeta} = \lambda_2(\tau_1\alpha_1 - \tau_2\alpha_2) \qquad (26) $$

where $\zeta$ is the dimensionless distance (z/d) along a column of length d, $\lambda_i$ are the decay constants, and $\tau_i = d/\omega_{eff}^i$ is the effective transport time of an element along a column of length d where the effective velocity $\omega_{eff}^i$ is:

$$ \omega_{eff}^i = \frac{\rho_f\phi\omega + \rho_s(1-\phi)D_iW}{\rho_f\phi + \rho_s(1-\phi)D_i} \qquad (27) $$

where $\rho_f$ and $\rho_s$ are the melt and solid densities, respectively, $\phi$ is the porosity, $\omega$ is the melt velocity, W is the solid mantle upwelling velocity, and $D_i$ is the bulk mineral/melt partition coefficient of element i.

For a one-dimensional steady-state melting column, the average melt velocity is:

$$ \omega = \frac{\Gamma d}{\rho_f |\phi|} \qquad (28) $$

where $\Gamma$ is the melting rate, which is assumed to be constant in the analytical approximation of Sims et al. [25] and can be determined using the equation:

$$ \Gamma = \frac{W\rho_s F_{max}}{d} \qquad (29) $$

where $F_{max}$ is the maximum degree of melting at the top of the melt column. Note that if element i is highly incompatible ($D_i \ll 1$), the effective velocity ($\omega_{eff}^i$) approaches the melt velocity ($\omega$). Furthermore, the difference in effective velocities between elements with different bulk partition coefficients decreases as the porosity increases.

Equations (24)–(26) comprise a coupled system of differential equations, where the decay of the parent isotope increases the concentration of its daughter isotope as the melt migrates through the column chromatographically, resulting in an enrichment ($\alpha_i$) of that daughter isotope. If both the melt velocity ($\omega$) and porosity ($\phi$) are constant throughout the melt column, this system of Equations (24)–(26) can be solved analytically to determine the enrichment factor $\alpha_i$ (Equation (23)) of the daughter isotopes due to chromatographic porous flow [25]. In this approach, the constant porosity is a local liquid-filled fraction defined with reference to melt migration and not by the overall melt fraction, which starts at zero at the solidus and increases linearly along the melt column.

Solving Equations (24)–(26), Sims et al. [25] obtained the following (where the subscripts 0, 1, and 2 refer to $^{238}$U, $^{230}$Th, and $^{226}$Ra, respectively):

$$ \alpha_0(\zeta) = \alpha_0^0 e^{-\lambda_0 \tau_0 \zeta} \qquad (30) $$
$$ \alpha_1(\zeta) = \frac{\lambda_1 \tau_0 \alpha_0^0}{\lambda_1 \tau_1 - \lambda_0 \tau_0} e^{-\lambda_0 \tau_0 \zeta} + \left[ \alpha_1^0 - \frac{\lambda_1 \tau_0 \alpha_0^0}{\lambda_1 \tau_1 - \lambda_0 \tau_0} \right] e^{-\lambda_1 \tau_1 \zeta} \quad (31) $$

$$ \alpha_2(\zeta) = \frac{\lambda_2 \lambda_1 \tau_1 \tau_0 \alpha_0^0}{(\lambda_1 \tau_1 - \lambda_0 \tau_0)(\lambda_2 \tau_2 - \lambda_0 \tau_0)} e^{-\lambda_0 \tau_0 \zeta} + \frac{\lambda_2 \tau_1}{\lambda_2 \tau_2 - \lambda_1 \tau_1} \left[ \alpha_1^0 - \frac{\lambda_1 \tau_0 \alpha_0^0}{\lambda_1 \tau_1 - \lambda_0 \tau_0} \right] e^{-\lambda_1 \tau_1 \zeta} + \left[ \alpha_2^0 - \frac{\lambda_2 \lambda_1 \tau_1 \tau_0 \alpha_0^0}{(\lambda_1 \tau_1 - \lambda_0 \tau_0)(\lambda_2 \tau_2 - \lambda_0 \tau_0)} - \frac{\lambda_2 \tau_1}{\lambda_2 \tau_2 - \lambda_1 \tau_1} \left( \alpha_1^0 - \frac{\lambda_1 \tau_0 \alpha_0^0}{\lambda_1 \tau_1 - \lambda_0 \tau_0} \right) \right] e^{-\lambda_2 \tau_2 \zeta} \quad (32) $$

where $\alpha_0^0$, $\alpha_1^0$, and $\alpha_2^0$ are constants of integration. For an initial state of secular equilibrium at the base of the melting column, all $\alpha_i^0 = 1$.

As shown in Sims et al. [25], this analytical solution can be further simplified by realizing that the parent isotopes generally are more compatible than the daughter isotopes (i.e., $D_U > D_{Th} >> D_{Ra}$ and $D_U >> D_{Pa}$) and, thus, over the length of the melt column, the change in the daughter nuclide's activity, $\alpha$, is large compared to the change in the parental $\alpha$ (i.e., $\Delta\alpha_U < \Delta\alpha_{Th} << \Delta\alpha_{Ra}$ and $\Delta\alpha_U << \Delta\alpha_{Pa}$). As such, one can hold the parental $\alpha_i$ values constant, further reducing Equations (31) and (32) to:

$$ \alpha_1(\zeta) = \alpha_1^0 e^{-\lambda_1 \tau_1 \zeta} + \frac{\tau_0}{\tau_1} \alpha_0 \left[ 1 - e^{-\lambda_1 \tau_1 \zeta} \right] \quad (33) $$

$$ \alpha_2(\zeta) = \alpha_2^0 e^{-\lambda_2 \tau_2 \zeta} + \frac{\tau_1}{\tau_2} \alpha_1 \left[ 1 - e^{-\lambda_2 \tau_2 \zeta} \right] \quad (34) $$

For fixed bulk partition coefficients, total melt fractions, and melt column lengths, the two analytical solutions above (that is, the full solution in Equations (30)–(32) versus the approximate analytical solution in Equations (33) and (34)) provide nearly identical results (Figure 9). The reason for this similarity is that the half-lives of $^{230}$Th and $^{226}$Ra are very different, resulting in $\lambda_1 \tau_1 << \lambda_2 \tau_2$ (where $\lambda\tau$ represents the ratio of the effective transport time with the time required for the concentration of the isotope to decrease by half), thus allowing the $^{238}$U-$^{230}$Th and $^{230}$Th-$^{226}$Ra systems to be treated as if they are decoupled.

The simplest way to use the chromatographic porous flow model with U-series disequilibria to obtain estimates of time-dependent parameters of melting is to use observational constraints to estimate the height of the melt column, the melt production rate, and the mineral/melt partition coefficients, and then to compare empirically measured U-series data with a forward model that calculates the porosity and melting rate. Both numerical and analytical solutions can be used to calculate the U-series activity ratios for the 1D chromatographic porous flow model of Spiegelman and Elliott [18]. These solutions show that ($^{226}$Ra/$^{230}$Th) disequilibrium is controlled mostly by the porosity of the melt region (which controls the velocity of the melt relative to the solid) for melting rates under $2 \times 10^{-2} \text{ kg} \cdot \text{m}^{-3} \cdot \text{yr}^{-1}$. In contrast, ($^{230}$Th/$^{238}$U) and ($^{231}$Pa/$^{235}$U) disequilibria are controlled more strongly by the melting rate, which in turn is related primarily to the solid mantle upwelling rate (Figures 9–11) and secondarily to the behavior of the melting lithology (see further discussion below).

## Comparison of Analytical Solutions with the Full Numerical Solution for Chromatographic Porous Flow Melting

Because the full equations for porous flow melting consider a depth-dependent porosity structure using a simplified form of Darcy's Law [18,55], whereas the analytical solutions assume a constant porosity over the length of the melting column [25], the absolute values of the calculated rate-determinant parameters can be significantly different, contingent on the initial assumptions of the porosity structure of the mantle. By assuming a constant porosity, the analytical solution may overestimate the true porosity, which most strongly influences ($^{226}$Ra/$^{230}$Th), while underestimating the solid mantle upwelling velocity, which more significantly influences ($^{230}$Th/$^{238}$U). This difference occurs because the full equations
take into consideration the details of the porosity distribution and melt velocity along the melt column, whereas the analytical approach treats these as constants (Figure 9).

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Results of reactive porous flow partial melting calculations for an adiabatically decompressing peridotitic mantle at T<sub>p</sub> = 1300 °C showing (<sup>226</sup>Ra/<sup>230</sup>Th) and (<sup>231</sup>Pa/<sup>235</sup>U) versus (<sup>230</sup>Th/<sup>238</sup>U) gridded melting model results for the solid mantle upwelling rates (W) and residual threshold porosity values (f) listed. Model results are shown for (A) equilibrium transport, (B) scaled transport with Da = 0.1, (C) disequilibrium transport with constant bulk solid/melt partition coefficients and melting rates, (D) equilibrium transport, (E) scaled transport with Da = 0.1, and (F) disequilibrium transport with variable bulk solid/melt partition coefficients and melting rates, as described in the text. Global MORB data are shown for comparison (see the caption in Figure 7 for the data sources).</figcaption></figure>

An alternative approach to analytical solutions that can be computed using relatively simple scripts or spreadsheets is incremental or discretized calculations, which can approximate non-linear solutions to calculate aggregated melt compositions along a melt column. However, incremental calculations introduce compounding errors that can only be avoided by using extremely small step sizes, which rapidly become computationally expensive. To date, some incremental approaches (e.g., [27,28]) have been used to calculate nuclide concentrations directly, but the method can be unwieldy; even small step sizes can sometimes fail to accurately track particularly rapid changes in nuclide concentration.

When variables such as porosity, melting rate, and mineralogically controlled trace element partitioning are expected to vary dramatically during the melting process, non-linear solutions to porous flow are needed to account for the expected variations in input parameters, and may produce significantly different outcomes in the melt composition (for examples of non-linear changes in input variables, consider Darcy's Law flow with compaction [18]; for melt productivity variations predicted by the MELTS family of models for decompressing peridotites, see Asimow [57]; for the partitioning effects from the sudden disappearance of a critical mineral such as garnet, which can significantly fractionate U from Th, see [54,58–63]). Hence, the solutions for melt composition by chromatographic porous flow melting in such scenarios are better approximated using iterative numerical routines, rather than simple progressive (Eulerian) time stepping methods that may am-
plify or compound small errors (e.g., [64]). Because highly incompatible trace element concentrations change very rapidly in some parts of the domain but are more stable in others in response to changes in melt fraction, partial melting models also tend to be very “stiff,” meaning that the dependent variable(s) change(s) extremely rapidly in response to small changes in the independent variable(s) (e.g., [65]). The method for accurately solving stiff equations numerically is to use extremely small step sizes to avoid overshooting the true solution, but this is computationally very expensive and is sometimes so extreme that it is not practical. Adaptive time-stepping routines, which use small step sizes only where changes are large, and implicit methods that solve for the entire melting regime as a continuous solution thus tend to be the most successful computational approaches (e.g., [66]).

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Results of reactive porous flow partial melting calculations for an adiabatically decompressing peridotitic mantle at $T_p = 1400$ °C showing $(^{226}_{230}\text{Ra})/(^{230}_{238}\text{Th})$ and $(^{231}_{235}\text{Pa})/(^{230}_{238}\text{Th})$ versus $(^{230}_{Th})/(^{238}_{238}\text{U})$ gridded melting model results for the solid mantle upwelling rates (W) and residual threshold porosity values (f) listed. Model results are shown for (A) equilibrium transport, (B) scaled transport with $Da = 0.1$, (C) disequilibrium transport with constant bulk solid/melt partition coefficients and melting rates, (D) equilibrium transport, (E) scaled transport with $Da = 0.1$, and (F) disequilibrium transport with variable bulk solid/melt partition coefficients and melting rates, as described in the text. Global MORB data are shown for comparison (see the caption in Figure 7 for the data sources).</figcaption></figure>

Spiegelman [55] developed a numerical solution to the equilibrium porous flow melting formulations by Spiegelman & Elliott [18] and produced UserCalc, an open-source, publicly accessible computational tool for conducting calculations from user input parameters. Elkins and Spiegelman [56] built on this prior work by developing an updated, open-source computational tool (pyUserCalc) that reproduces the functionality of UserCalc with added computational methods based on the Appendix of Spiegelman and Elliott [18] to calculate U-series disequilibria produced in partial melts by disequilibrium (fractional) porous flow. Their method also permits near-fractional melting using a scaled reactivity rate, expressed as a Damköhler number (where the Damköhler number is the ratio of a
chemical reaction rate—in this case, the mass chemical exchange rate between solid and liquid—to a physical transport rate, here, the solid mantle upwelling rate).

Recent results for chromatographic melting by equilibrium, disequilibrium, and partial equilibrium reactive porous flow, using a chemical equilibration rate scaled by a Damköhler number, illustrate that partial melting scenarios for an average or depleted mantle peridotite, using experimentally expected mineralogy, thermodynamically predicted melt productivity, and changes in porosity with depth, are unable to generate as large a range in ($^{230}$Th/$^{238}$U) as the analytical solution predicts for mid-ocean ridge environments [25,67,68], even for a sufficiently hot mantle to place the peridotite solidus and thus the onset of melting within the garnet peridotite stability field (e.g., a mantle potential temperature ($T_p$) ≥ 1400 °C). Specifically, predicted final melt compositions have low ($^{230}$Th/$^{238}$U) close to 1.0 for most peridotite melting scenarios. This is partly due to the low initial melt productivities predicted by thermodynamic modeling, such that most or all melting occurs in the spinel stability field at typical mid-ocean ridge mantle temperatures. In chemical equilibrium transport scenarios, signatures from small amounts of melting in the garnet stability field may then be overprinted by subsequent progressive melting that is garnet-free. Higher temperatures hypothesized for hotspot locations produce a broader range of outcomes because of enhanced melting in the presence of residual garnet, but the overall mismatch with analytical solutions is noteworthy and suggests that minor amounts of other garnet-bearing lithologies may play a role on a global scale (see further discussion of heterogeneous melting effects below) [67,68].

Alternatively, Krein et al. [69] considered partial melting of relatively cold ($T_p$ = 1200 °C) lherzolitic mantle with stable plagioclase in the shallow melting regime and suggested that it is in fact possible to generate relatively high ($^{230}$Th/$^{238}$U) in partial melts, without requiring deep melting in the presence of garnet. While these outcomes are reproducible by multiple melting models, plagioclase is unlikely to be stable in adiabatic regimes except at mid-ocean ridges with rapid spreading rates, where the lithosphere is very thin, suggesting such a scenario cannot fully explain all global basalt compositions, and non-peridotitic source lithologies are still likely necessary in some settings [68].

To better illustrate the model differences explored above, Figures 10 and 11 show a series of outcomes comparing model results for two adiabatic melting scenarios, each considering a lherzolitic mantle undergoing solid decompression along a moderately warm adiabat that crosses the lherzolite solidus. The first scenario considers a moderate $T_p$ of 1300 °C, for which melting initiates in the spinel peridotite stability field, while the other uses a higher $T_p$ of 1400 °C with a solidus depth and onset of melting within the garnet peridotite stability field. At each potential temperature, the outcomes for constant elemental partition coefficients and a constant melting rate were computed first, to directly compare the simplified analytical approach with more complex numerical solutions that permit changes to porosities with depth and thus variable liquid flow rates. A second set of calculations for variable elemental partitioning and melting rates with depth using only numerical models was then considered, using the methods of Lambart [70], Elkins et al. [67], and Elkins and Lambart [68], where melt fractions with depth are determined using the Melt-PX parameterization developed from high-pressure melting experiments [71] and accompanying solid mineral proportions are computed using the thermodynamic pMELTS model at a given pressure and liquid fraction [72]. The fixed mineral/melt partition coefficients used to determine bulk elemental partitioning for these examples are provided in Table 3.
<table><caption>**Table 3.** Default mineral/melt partition coefficients (Di) used to calculate bulk rock partition coefficients in model calculations for this study. Modified from Elkins and Lambart [68].</caption><thead><tr><th>Lithology</th><th>Phase</th><th>D<sub>U</sub></th><th>D<sub>Th</sub></th><th>D<sub>Pa</sub>*</th><th>D<sub>Ra</sub>*</th><th>Reference</th></tr></thead><tbody><tr><td>Garnet Peridotite</td><td>Garnet</td><td>0.038</td><td>0.017</td><td>0.00001</td><td>0.00001</td><td>RD 1097-5 experiment [73]</td></tr><tr><td rowspan="5">KG1 Pyroxenite</td><td>Clinopyroxene</td><td>0.0030</td><td>0.0040</td><td>0.00001</td><td>0.00001</td><td>RD 1097-5 experiment [73]</td></tr><tr><td>Olivine</td><td>0.00005</td><td>0.00047</td><td>0.00001</td><td>0.00001</td><td>RD 1097-5 experiment [73]</td></tr><tr><td>Orthopyroxene</td><td>0.0078</td><td>0.0086</td><td>0.00001</td><td>0.00001</td><td>TM0500-3 experiment [73]</td></tr><tr><td>Pigeonite</td><td>0.00050</td><td>0.00018</td><td>0.00001</td><td>0.00001</td><td>W&amp;P 122-1 experiment [63]</td></tr><tr><td>Spinel</td><td>0.012</td><td>0.0024</td><td>0.00001</td><td>0.00001</td><td>Lunar basalt [74]</td></tr><tr><td rowspan="6">Spinel Peridotite</td><td>Cpx</td><td>0.0080</td><td>0.0070</td><td>0.00001</td><td>0.00001</td><td>TM 1094-9 experiment [75]</td></tr><tr><td>Olivine</td><td>0.00005</td><td>0.00047</td><td>0.00001</td><td>0.00001</td><td>RD 1097-5 experiment [73]</td></tr><tr><td>Orthopyroxene</td><td>0.0024</td><td>0.0027</td><td>0.00001</td><td>0.00001</td><td>RD 1097-2 experiment [73]</td></tr><tr><td>Plagioclase</td><td>0.0006</td><td>0.0034</td><td>0.00001</td><td>0.02000</td><td>D<sub>U</sub>, D<sub>Th</sub> calculated after [54];<br>D<sub>Ra</sub> from [76]</td></tr><tr><td>Pigeonite</td><td>0.0005</td><td>0.00018</td><td>0.00001</td><td>0.00001</td><td>W&amp;P 122-1 experiment [63]</td></tr><tr><td>Spinel</td><td>0.012</td><td>0.0024</td><td>0.00001</td><td>0.00001</td><td>Lunar basalt [74]</td></tr><tr><td rowspan="7">G2 Pyroxenite</td><td>Garnet</td><td>0.02405</td><td>0.00415</td><td>0.00001</td><td>0.00001</td><td>A343 experiment [77]</td></tr><tr><td>Clinopyroxene</td><td>0.0041</td><td>0.0032</td><td>0.00001</td><td>0.00001</td><td>A343 experiment [77]</td></tr><tr><td>Olivine</td><td>0.00005</td><td>0.00047</td><td>0.00001</td><td>0.00001</td><td>RD 1097-5 experiment [73]</td></tr><tr><td>Orthopyroxene</td><td>0.0078</td><td>0.0086</td><td>0.00001</td><td>0.00001</td><td>TM0500-3 experiment [73]</td></tr><tr><td>Plagioclase</td><td>0.0006</td><td>0.0034</td><td>0.00001</td><td>0.02000</td><td>D<sub>U</sub>, D<sub>Th</sub> calculated after [54];<br>D<sub>Ra</sub> from [76]</td></tr><tr><td>Pigeonite</td><td>0.0096</td><td>0.010</td><td>0.00001</td><td>0.00001</td><td>Experiment 18 [78]</td></tr><tr><td>Spinel</td><td>0.046</td><td>0.016</td><td>0.00001</td><td>0.00001</td><td>Maximum measured [79]</td></tr><tr><td>MIX1G, m7-16</td><td>Garnet</td><td>0.013</td><td>0.0032</td><td>0.00001</td><td>0.00001</td><td>Experimental results [79]</td></tr><tr><td rowspan="7">Pyroxenite</td><td>Clinopyroxene</td><td>0.017</td><td>0.015</td><td>0.00001</td><td>0.00001</td><td>Experimental results [79]</td></tr><tr><td>Olivine</td><td>0.00005</td><td>0.00047</td><td>0.00001</td><td>0.00001</td><td>RD 1097-5 experiment [73]</td></tr><tr><td>Spinel</td><td>0.046</td><td>0.016</td><td>0.00001</td><td>0.00001</td><td>Maximum measured [79]</td></tr><tr><td>Orthopyroxene</td><td>0.0078</td><td>0.0086</td><td>0.00001</td><td>0.00001</td><td>TM0500-3 experiment [73]</td></tr><tr><td>Plagioclase</td><td>0.0006</td><td>0.0034</td><td>0.00001</td><td>0.02000</td><td>D<sub>U</sub>, D<sub>Th</sub> calculated after [54];<br>D<sub>Ra</sub> from [76]</td></tr><tr><td>Pigeonite</td><td>0.0096</td><td>0.010</td><td>0.00001</td><td>0.00001</td><td>Experiment 18 [78]</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

At $T_p = 1300$ °C, the model results in Figure 10 illustrate that regardless of whether input variables are simplified or held constant, very little variation in $(^{230}\text{Th}/^{238}\text{U})$ can be achieved by peridotitic melting. For chemical disequilibrium transport (i.e., fractional melting), constant porosity, partitioning, and melting rates overall generate lower $(^{226}\text{Ra}/^{230}\text{Th})$ and $(^{231}\text{Pa}/^{235}\text{U})$ than full numerical solutions to porous flow melting, but for equilibrium porous flow scenarios, there is little difference between the two outcomes at relatively low mantle temperatures. As explained above, more significant differences in outcomes are observed at higher temperatures where $T_p = 1400$ °C, with particularly broad ranges in $(^{230}\text{Th}/^{238}\text{U})$ for equilibrium melting with constant input parameters, and again particularly low $(^{226}\text{Ra}/^{230}\text{Th})$ and $(^{231}\text{Pa}/^{235}\text{U})$ during disequilibrium transport with constant inputs (Figure 11).

### 3.4. Dynamic Versus Chromatographic Melting: Comparison of the Two Types of Time-Dependent Models

Dynamic melting models provide results that are similar in some ways to those obtained from the chromatographic porous flow melting model in that the extent of $(^{226}\text{Ra}/^{230}\text{Th})$ disequilibrium is controlled mainly by the porosity of the melt region, and $(^{230}\text{Th}/^{238}\text{U})$ and $(^{231}\text{Pa}/^{235}\text{U})$ disequilibria are controlled mainly by the melting rate (Figure 8). In all cases, the values of these parameters are controlled by the bulk partition coefficients ($D_i$) used for the elements U, Th, Ra, and Pa, which in turn depend on residual mineralogy, which changes with pressure in the mantle.

In the traditional formulation of McKenzie [20], dynamic melting differs from chromatographic melting in that the produced melts are assumed to move instantly to the

<footnote>* By convention, $D_{Pa}$ and $D_{Ra}$ are set equal to very small values ($1 \times 10^{-5}$) for most mantle minerals, except $D_{Ra}$ in plagioclase.</footnote>
surface as they form, and, consequently, they do not react with the solid during transport or experience further radioactive decay after extraction. Because most of the incompatible trace elements, including uranium, are extracted in the earliest stages of partial melting, disequilibria are dominantly created at the base of the melt column instead of throughout the entire length of the melt column, as occurs during chromatographic melting. Unlike chromatographic porous flow, where melt velocities are determined explicitly, with dynamic melting, melt velocities are only constrained by the need to preserve the shortest half-life of the nuclide in radioactive disequilibrium with its immediate parent. Generally, the half-life of ²²⁶Ra is used in this context. A few studies have argued that observed ²¹⁰Pb excesses place even shorter limits on the timescale of melt extraction, requiring more rapid melt ascent velocities. However, as discussed below, (²¹⁰Pb/²²⁶Ra) fractionation is probably produced by ²²²Rn accumulation or depletion (e.g., [53,80]).

In other words, in traditional dynamic melting [20], the timescale of melt migration is not explicitly determined or considered. Therefore, for a given set of partition coefficients, the dynamic melting model requires threshold or escape porosities that are lower than the maximum porosity at the top of the melt zone predicted by a chromatographic melting model. With dynamic melting, the inferred melting rates are also slightly lower. This stems from the shorter residence time of the parent nuclide in the melting column obtained with dynamic melting for a given melting rate. To achieve the same residence time for the parent nuclide (i.e., same time for in-growth of the daughter nuclide), the melting rate must be slower. To place quantitative constraints on magma extraction rates, later formulations for dynamic melting, such as the incremental model of Stracke et al. [28], explicitly calculated the outcomes of different ascent velocities for extracted magmas and then compared the model predictions to actual data.

In the dynamic melting formulations of McKenzie [20], the liquid trapped in pore spaces below the critical or threshold porosity is also assumed to reach chemical equilibrium with the solid residue, regardless of the time elapsed, while any melt fraction beyond the threshold porosity is extracted instantaneously and does not react with solids in the rest of the melting column. Beyond the critical threshold porosity, the porosity is assumed to be constant throughout the melting zone. The model thus makes two simultaneous but opposing assumptions about melt-rock chemical equilibration and magma extraction, while neglecting compaction effects on fluid velocity. While this formulation produces partial melts with incompatible element enrichments that are intermediate between equilibrium and accumulated fractional melts [17,31], and the extraction rates of incompatible elements are partly controlled by the porosity, mechanistically the underlying assumptions about chemical equilibration and transport do not effectively simulate significant aspects of the chemical behavior of a system experiencing coupled two-phase (solid and liquid) flow.

In comparison, porous flow models [18] may range from equilibrium porous flow, which has an infinite Damköhler number (i.e., melt-solid reaction is very fast relative to melt migration velocities), to scaled porous flow with a moderate Damköhler number (which calculates a continuous rate-limited or incomplete chemical exchange between solid and liquid), to a “zero-Damköhler number” model, i.e., truly time-dependent fractional melting where no chemical exchange beyond initial melting occurs because of the comparatively rapid magma transport rate [56]. The dynamic melting formulations with threshold equilibrium porosities are closest to the moderate Damköhler number scenario, but computed outcomes between the two are distinct [68].

For a full comparison with the chromatographic modeling results explored above, we here present equivalent results for adiabatic decompression partial melting of a mantle lherzolite, using both a simplified analytical solution (after [31]) and numerical modeling outcomes for dynamic melting (after [24,68,81,82]; Figure 8). The outcomes of the
numerical approach are comparable to those for incremental models with very small step sizes [27,28,68]. Consistent with prior results, it is observed that dynamic melting outcomes in Figure 8 generally produce higher ($^{230}$Th/$^{238}$U) in partial melts for a given garnet peridotite melting scenario than the disequilibrium or scaled equilibrium porous flow models in Figure 8, particularly for higher T$_{P}$ = 1400 °C with stable residual garnet. This outcome likely occurs because (1) the initial dynamic melting fraction below the critical threshold reaches perfect chemical equilibrium; (2) progressive melting in the spinel stability field has even less of an overprinting effect during dynamic melting; and (3) the porosity is fixed throughout the melt regime, which affects magma transport rates [67,68,82]. Because of the simplified mechanisms assumed for chemical equilibration in the dynamic melting formulation, however, we suggest that the disequilibrium porous flow transport model results shown in Figures 10 and 11 are more realistic in simulating the outcomes of magma generation with rapid magma transport and restricted or limited chemical exchange, even though the more restricted range of $^{238}$U-$^{230}$Th outcomes produced by the porous flow models place more stringent constraints on the process of mantle peridotite melting.

Finally, despite the observed differences in melting outcomes between the analytical and numerical approaches to both chromatographic and dynamic melting, differences which may be important in some geologic settings, the simplified approaches that permit analytical solutions may also be particularly useful for placing additional constraints on the melting regime. This is because analytical formulas are amenable to inverse modeling methods, where observed lava compositions are used to directly calculate the required starting conditions, rather than comparing systematic sets of model predictions to data sets using sensitivity analyses for specific variables, such as mantle upwelling rate (W), residual porosity (φ), or degree of melting (F). The inverse method is similar to approaches such as those of Sims et al. [26], who estimated bulk mineral/melt partition coefficients for U-series nuclides by inverting simplified model formulas that assume constant melting rates and partitioning behavior. Because the chromatographic model considers the timescale of melt migration and the analytical solution assumes a constant porosity, the velocity of the melt can in fact be determined explicitly (Equation (28)). With a numerical approach to chromatographic melting, on the other hand, melt velocity varies as a function of the porosity distribution of the melt column and thus cannot be directly determined from the final melt compositions. In the latter case, it is instead more useful to calculate melt transport times by integrating the melt velocity over the length of the melt column.

### 3.5. *Sources of Uncertainty in Melting Models*

There are four main sources of uncertainty in the interpretation of U-series measurements for modeling melting processes, namely, (1) modification of lava compositions by post-eruptive contamination; (2) uncertainties in sample ages and magma storage times; (3) uncertainties in measured U, Th, Pa, and Ra partition coefficients; and (4) unknown modification of magma compositions by assimilation, particularly lithospheric melt-rock-brine reactions.

#### 3.5.1. Uncertainties Introduced by the Modification of Lava Compositions by Post-Eruptive Contamination

Secondary alteration or post-eruptive contamination has the potential to significantly fractionate parent-daughter nuclide pairs and perturb U-series. Contamination can be particularly problematic in submarine lavas, as ocean waters contain relatively high concentrations of U [6] due to uranium's high solubility as a hexavalent cation in ocean waters. However, seawater alteration of submarine basalts is reasonably easy to discern because the ($^{234}$U/$^{238}$U) in ocean waters ubiquitously is 1.14 owing to uranium's long residence time [83–86]. Thus, any basalt that has been altered by seawater will have ($^{234}$U/$^{238}$U)
greater than unity (Figure 12). The threshold of what the petrological community considers a tolerable excess of $^{234}$U as an indication of the absence of alteration has changed considerably over time with the improved ability to measure abundance sensitivity, given that down mass tailing of $^{235}$U on $^{234}$U can increase the measured ($^{234}$U/$^{238}$U). Hence, it remains uncertain whether (a) many of the studies showing $^{234}$U excesses in Figure 12 are an indication of imprecise ($^{234}$U/$^{238}$U) measurements and are viably interpretable in terms of petrological processes, or (b) if those samples have been altered and should be excluded from petrological modeling. A classic example of how even minor $^{234}$U excesses can indicate perturbation of ($^{230}$Th/$^{238}$U) in depleted basalts via alteration is the Kolbeinsey ridge (Figure 13), where there is a clear correlation between ($^{230}$Th/$^{238}$U) and ($^{234}$U/$^{238}$U) [50].

Because of the long half-life of $^{234}$U ($t_{1/2}$ = 245 kyr), for young subaerial samples, ($^{234}$U/$^{238}$U) = 1 is a necessary but insufficient condition for ruling out secondary alteration. In this case, other indices such as K$_2$O/P$_2$O$_5$ > 1 can be used to indicate the presence or absence of chemical alteration.

### 3.5.2. Uncertainties in Sample Ages and Magma Storage Times

Because of the relatively short half-lives of the U-series daughter isotopes, as samples age either on the surface or during magma storage, the parent/daughter activity ratios will decay toward their equilibrium values. If the sample age is old and unknown, or if the magma storage interval is significant and not accounted for, one would infer parent/daughter fractionations during the melting process that are too small. Therefore, it is important to explicitly know a sample's age, and to have some geologic constraints on magma storage time when evaluating melting processes using U-series disequilibria. This issue is especially significant for $^{226}$Ra, which has a short half-life (1600 ± 7 yr) relative to the expected timescales of melting, melt migration, and magma storage processes, and, since $^{226}$Ra in partial melts is largely controlled by the porosity of the melt zone, uncertainties in sample ages and magma storage times translate into large uncertainties in the calculated porosity. However, because of the differences in $^{226}$Ra, $^{230}$Th, and $^{231}$Pa half-lives, if a lava sample has a primary (that is, melting-derived) $^{226}$Ra excess relative to $^{230}$Th, any age effects are small enough to be inconsequential for measured $^{238}$U-$^{230}$Th or $^{235}$U-$^{231}$Pa disequilibria. Similarly, if a sample has a $^{210}$Pb ($t_{1/2}$ = 22.6 ± 0.1 yr) excess or deficit relative to $^{226}$Ra, then age affects are not likely to be an issue for $^{230}$Th, $^{231}$Pa, or $^{226}$Ra.

An important nuclide that could be used to address uncertainties associated with sample age is $^{227}$Ac, a daughter product of $^{231}$Pa with a half-life of 21.77 yr, similar to $^{210}$Pb. Unlike $^{210}$Pb, which has intermediate volatile species such as $^{222}$Rn in the $^{226}$Ra-$^{210}$Pb decay scheme, which could be affected by gas exsolution and either loss or accumulation in the magma (Figure 2), $^{227}$Ac is a direct decay product of $^{231}$Pa and is not volatile. However, because of the difficult analytical methods required to measure $^{227}$Ac [87], only a few studies [44,88,89] have reported $^{231}$Pa-$^{227}$Ac disequilibria in erupted magmas that were age-constrained (i.e., zero age at the time of collection), and, so far, all analyzed magmas have ($^{227}$Ac/$^{231}$Pa) of unity, suggesting that either (1) the timescale of magma transport and storage is at least 100 years or (2) there was no initial fractionation of ($^{227}$Ac/$^{231}$Pa).
<figure><img src="image_16.png" /><figcaption>**Figure 12.** ($^{234}$U/$^{238}$U) vs. ($^{230}$Th/$^{238}$U) for the global MORB data set (see the caption in Figure 7 for the data sources [14,19,25,32–53]), including (A) the Mid-Atlantic Ridge; (B) East Pacific Rise; and (C) Juan de Fuca Ridge, Gorda Ridge, Southwest Indian Ridge, and Southeast Indian Ridge. Error bars are 2σ analytical uncertainties. Note that many samples show ($^{234}$U/$^{238}$U) out of equilibrium, indicative of seawater contamination.</figcaption></figure>

### 3.5.3. Uncertainties in Measured U, Th, Ra Partition Coefficients

Porosities and melting rates inferred from partial melting predictions using both chromatographic and dynamic melting models are highly dependent on the absolute and relative values chosen for the bulk mineral/melt U, Th, Pa, and Ra partition coefficients. The period of time during which the parent element resides in the melt column is critical to the outcomes of time-dependent ingrowth models; therefore, these models are especially sensitive to the $D_i$ values chosen for the parent elements U and Th. For the combined $^{238}$U-$^{230}$Th-$^{226}$Ra disequilibria, the value chosen for $D_{Th}$ is also particularly important, as $^{230}$Th is both a parent and daughter isotope, and therefore, its D value affects both ($^{230}$Th/$^{238}$U) and ($^{226}$Ra/$^{230}$Th) disequilibria. For example, with the chromatographic equilibrium porous flow model, decreasing the value of $D_{Th}$ by a factor of 3 (while holding $D_U$ and $D_{Ra}$ constant) changes the reference porosity required to achieve a given combination of ($^{230}$Th/$^{238}$U) and
$(^{226}_{230}\text{Ra}/^{230}\text{Th})$ by up to a factor of 3 for low melting rates, and by an order of magnitude for high melting rates (e.g., [25]).

<figure><img src="image_17.png" /></figure>

**Figure 13.** Uranium isotope calculated mixing line between seawater [($^{234}$U/$^{238}$U) = 1.14] and unaltered local Kolbeinsey melt compositions. This correlation based on high-precision mass spectrometric analyses demonstrates the need to have ($^{234}$U/$^{238}$U) in equilibrium within the uncertainty of modern mass spectrometric measurements.

The bulk melt-rock partitioning behavior of U, Th, Ra, and Pa during melting may primarily be controlled by the presence and modes of some of the major mantle mineral phases (i.e., olivine, pyroxene, and garnet), and thus by both the lithology and bulk composition of the source rock and by the pressure and temperature of melting (e.g., [54] and references therein). Clinopyroxene and garnet typically host the highest concentrations of U and Th in peridotitic rocks (e.g., [54] and references therein), and specific mineral/melt partition coefficients in turn depend on the chemical compositions of these minerals, particularly their octahedral site radii, which are controlled by bulk rock composition and temperature and pressure conditions (e.g., [79]). A compilation of experimentally determined D<sub>U</sub> and D<sub>Th</sub> values in typical mantle minerals illustrates that estimated partition coefficients can vary by an order of magnitude due to experimental methods, measurement uncertainty, and differences in mineral composition and pressure (see Table 3 and references therein). Such large variations in clinopyroxene and garnet D<sub>U</sub> and D<sub>Th</sub> in turn can produce large variations in model outcomes (e.g., [25,50]), and, therefore, also in the calculated or inferred porosities and melting rates for both chromatographic and dynamic melting models (Figures 8–11).

An additional uncertainty is the role of grain boundaries as an important repository for highly incompatible elements (e.g., [90–92]). These could play a role similar to that of retained melt, increasing bulk distribution coefficients and thus influencing estimates of porosity and transport rates required to produce a given set of disequilibria. At equilibrium, grain boundaries in mantle peridotites could contain on the order of 10,000 times higher concentrations of Th, and 50,000 times more U, than the crystalline minerals ([91], Figure 2 and accompanying text). For typical peridotites with a grain size of 1 to 0.1 mm and grain boundary widths ~1 nm, assuming cubic grains, grain boundaries could comprise ~0.6 to
6 ppmv (perhaps 0.5 to 5 ppmw) of the rock, with a melt/(grain boundary) distribution coefficient ~10⁴ higher than distribution coefficients for the crystalline minerals, increasing bulk distribution coefficient values by ~0.005 to 0.05.

As a result of the relatively high sensitivity of model outcomes to comparatively small changes in partitioning behavior, solid mantle upwelling rates and porosities can only be estimated from measured (²³⁰Th/²³⁸U), (²²⁶Ra/²³⁰Th), and (²³¹Pa/²³⁵U) disequilibria (e.g., [25]). That said, while the implied uncertainty in model outcomes means that precise upwelling rates and reference porosities cannot be inferred quantitatively from time-dependent melting models, sensitivity analyses of the impacts of varying the upwelling rate and reference porosity on (²³⁰Th/²³⁸U), (²²⁶Ra/²³⁰Th), and (²³¹Pa/²³⁵U) consistently produce qualitatively and even quantitatively different outcomes for distinct melting scenarios, making it possible nevertheless to reliably interpret the results of melting models (e.g., [68]).

### 3.5.4. Modification of Magma Compositions by Assimilation

Magmas that pool in or flow through a partially crystalline reservoir with continuous or episodic recharge may also assimilate or undergo chemical exchange and reaction with crystallized material from prior magma injection episodes (e.g., [93,94]). Continuous mineral precipitation and dissolution, as well as chemical or diffusive exchange, are very complex to model accurately due to the open nature of the system and lack of necessary constraints on the minerals present or the timing of magma recharge; however, such scenarios could potentially impact or modify quite significantly U-series disequilibria (e.g., [95,96]). For example, continuous and ongoing chemical exchange between a specific group of minerals and a surrounding magma could modify or buffer the disequilibrium between a parent-daughter isotope pair, making it more difficult to draw unique conclusions about magmatic age and storage times or mineralogy, and potentially obscuring the underlying melting signal (e.g., [96]). However, other evidence for melt-rock interactions and mineral-liquid exchange can help to constrain potential model scenarios. For example, while (²²⁶Ra/²³⁰Th) and (²³⁰Th/²³⁸U) can be modeled in isolation, constraints derived from other evidence for cumulate mineral-liquid exchange limit potential model scenarios when other trace elements are considered. Similarly, ²¹⁰Pb excesses could be modeled by preferential melt-cumulate rock exchange by diffusive processes involving plagioclase or sulfide minerals, as Pb is expected to behave more compatibly than Ra in both plagioclase and sulfide minerals; but again, this particular scenario would be expected to result in correlations between geochemical indices of plagioclase or sulfide melting and (²¹⁰Pb/²²⁶Ra) disequilibria and are generally not observed [53].

## 3.6. Mid-Ocean Ridge Basalts (MORB)

²³⁸U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹Pa disequilibria have been measured in MORB sample suites by high-precision mass spectrometric methods across the full range of ocean ridge spreading rates from fast to ultraslow (Figure 7), radically transforming our understanding of MORB petrogenesis. The following discusses a few important observations drawn from individual regional and case studies for mid-ocean ridge settings and presents an updated model for global MORB genesis with broader implications for the dynamics of mantle melt generation and oceanic crustal construction.

### 3.6.1. Global Correlation of (²³⁰Th/²³⁸U) with Ridge Axial Depth

An early first-order, seminal observation of U-series fractionation in MORB was that (²³⁰Th/²³⁸U) disequilibria vary systematically with ridge axial seafloor depth [34], which, in turn, principally is controlled by oceanic crustal thickness [97]. This correlation with (²³⁰Th/²³⁸U) was attributed to the fact that the depth of the peridotite solidus is controlled
by the mantle potential temperature, which has been predicted as a controlling factor in producing varying crustal thicknesses and, by inference, basaltic geochemical compositions [97]. Using time-dependent melting models, Bourdon et al. [34] explained the broad global negative correlation between axial ridge depth and ($^{230}$Th/$^{238}$U) by inferring mantle temperature control over melt column depth (i.e., the depth to the peridotite solidus). Because the solidus depth controls the length of the melt column, this depth (whether inferred or an explicit model parameter) influences both the duration of melting during upwelling and the overall degree and quantity of melt produced in the garnet peridotite stability field, where $D_U > D_{Th}$, producing ($^{230}$Th/$^{238}$U) > 1 in the magma [58–61, 63, 73, 75, 98, 99].

An important location to see this systematic relationship between ($^{230}$Th/$^{238}$U) and seafloor axial depth is the Arctic or North Atlantic mid-ocean ridge system [50,100] (Figure 14). The North Atlantic data set nearly spans the complete range of both axial depth (250–4131 m) and ($^{230}$Th/$^{238}$U) (0.955–1.296) for MORB globally. The high ($^{230}$Th/$^{238}$U) values measured in MORB from the Kolbeinsey [50], Mohns, and Knipovich Ridges [100] reflect melting in the presence of residual garnet. In contrast, basalts from the Eastern Volcanic Zone (EVZ) of the Gakkel Ridge erupt through very thin crust [101], and major element data (e.g., high MgO and Na₂) suggest that they are the result of smaller degrees of melting [102,103]. The overall low ($^{230}$Th/$^{238}$U) measured in EVZ Gakkel MORB suggests a relatively shallow melt regime dominated by spinel peridotite [100,101]. The negative trend of axial depth versus ($^{230}$Th/$^{238}$U) initially observed and modeled by Bourdon et al. [34] is also considerably stronger for Arctic MORB than for the global data set. Elliott and Spiegelman [104] measured a weakly correlated $r^2$ coefficient of 0.58 for all global ridge data, while the Arctic ridges when taken alone [50,100] have ($^{230}$Th/$^{238}$U) values that are more strongly correlated with axial depth ($r^2$ = 0.72).

The scatter observed in the global data set likely represents the influence and importance of several additional variables on melting, such as solidus depth, source heterogeneity, spreading rate, and the efficiency of melt extraction. The Kolbeinsey Ridge is compositionally more uniform and likely representative of elevated mantle temperature effects. The Kolbeinsey Ridge has MORB with notably high ($^{230}$Th/$^{238}$U) and shallow depths that are related to elevated crustal thickness adjacent to the Iceland hotspot. These data may indicate enhanced melting in the garnet peridotite stability field due to higher local mantle temperatures [50]. Basalts from ridges farther to the north (e.g., the Mohns and Knipovich Ridges), on the other hand, exhibit radiogenic isotope and trace element compositions indicative of greater source heterogeneity, suggesting that other lithologic sources of garnet such as pyroxenite may have increased their ($^{230}$Th/$^{238}$U) compositions and also the crustal thickness [100,105] (see below for further discussion of heterogeneity effects).

### 3.6.2. Two-Porosity Modeling

Two important observations indicate that MORB requires melt transport under heterogeneous conditions, including at least two porosities that generate distinct magma flow rates under equilibrium porous flow conditions: (1) trace element ratios indicate that MORB are mixtures of small-degree and large-degree melts, as demonstrated by the relationship between ($^{238}$U/$^{230}$Th) versus $\alpha_{Sm/Nd}$ (Figure 15; [14]), and, even more definitively, (2) a well-defined negative correlation exists between ($^{230}$Th/$^{238}$U) and ($^{226}$Ra/$^{230}$Th) for demonstrably young basalts (<375 years based on the width of the axial spreading trough (AST) and spreading rate) derived from an isotopically homogenous source at 9–10° N East Pacific Rise (Figure 16; [26]). This interpretation, that the magma source and transport regime is polybaric, was first proposed by Kelemen et al. [106] based on observations from ophiolite complexes that show melt channels and evidence for melt-rock reactions
plus the observed inverse correlation between (²²⁶Ra/²³⁰Th) and (²³⁰Th/²³⁸U) in basalts of unknown age from Juan de Fuca and Gorda ridges.

<figure><img src="image_18.png" /></figure>

**Figure 14.** Global (²³⁰Th/²³⁸U) versus Axial Depth. Diagram of measured (²³⁰Th/²³⁸U) versus axial ridge depth showing data from this study (Mohns, Knipovich, and Gakkel Ridges) and the global ridge system [26,34,41,47,50,51,107–110]. This data set is restricted to age-constrained unaltered (²³⁴U/²³⁸U = 1) oceanic samples with (²²⁶Ra/²³⁰Th) and/or (²³¹Pa/²³⁵U) out of equilibrium or stratigraphic age constraints. The curves show results of peridotite reactive porous flow melting from UserCalc [55] after calculations by Bourdon et al. [34] for various representative melting rates (M) and residual porosities (φ) as labeled. See Elkins et al. [100] for full details.

Early multiple-porosity interpretations (Figures 17 and 18) invoked 1D chromatographic melting with two distinct fluid porosity regimes and suggested a modeling approach to calculate the outcomes of partial melting in such a scenario [23,111,112]. In the Jull et al. [112] model, which was developed to explain the negatively correlated (²³⁰Th/²³⁸U) versus (²²⁶Ra/²³⁰Th) in age-constrained lavas from the 9–10° N East Pacific Rise, the end-member lavas (enriched MORB [E-MORB]) with high (²³⁰Th/²³⁸U) and low (²²⁶Ra/²³⁰Th) derived from deep in the garnet stability field and were transported to the surface on a timescale that allowed the (²²⁶Ra/²³⁰Th) to decay back to equilibrium. Lavas (depleted MORB [D-MORB]) with low (²³⁰Th/²³⁸U) and high (²²⁶Ra/²³⁰Th), on the other hand, derived from the shallow part of the melt column, where the youngest melt component formed via continued decompression melting in the spinel peridotite stability field, thus preserving the short-lived ²²⁶Ra signal.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** $^{230}$Th excess versus alpha Sm/Nd for Hawaiian samples and axial 9–10° N, Juan de Fuca, Garrett transform, and Siqueiros transform samples [14,25,26,41]. Modified from Sims et al. [25]. Alpha Sm/Nd is a measure of the Sm/Nd fractionation during MORB petrogenesis and is the Sm/Nd concentration ratio measured in the lava normalized to a model Sm/Nd source ratio inferred from the measured $^{143}$Nd/$^{144}$Nd of the lavas (for details see [14,26]).</figcaption></figure>

Because the lavas from this correlated sample suite preserve relatively uniform ra-
diogenic isotope signatures ($^{143}$Nd/$^{144}$Nd varies from 0.512378 ± 11 to 0.512400 ± 10,
$^{87}$Sr/$^{86}$Sr varies from 0.70244 ± 5 to 0.70257 ± 6, and $^{176}$Hf/$^{177}$Hf varies from 0.283169 ± 5
to 0.283197 ± 5), source heterogeneity cannot be the only factor influencing $^{238}$U-$^{230}$Th-
$^{226}$Ra systematics in mid-ocean ridge settings, and magma transport from a progressive
melting regime must also be a major control. As is typical of models with multiple vari-
able parameters, both chromatographic and dynamic melting (Figure 18) can be used
to successfully predict the inverse correlation at the 9–10° N East Pacific Rise for a ho-
mogeneous initial source [28,112], provided that the former models include channels of
focused melt flow, and the latter models include efficient melt extraction throughout a 2D,
“triangular” melting region. Evidence from ophiolites demonstrates the presence of focused
flow of MORB-like melt in dunite dissolution channels, surrounded by highly depleted
harzburgites that are very far from equilibrium with MORB-like melts (e.g., [106,113]),
providing observational support for chromatographic models that include conduits for
focused melt transport.

Recent outcomes of systematic modeling of global MORB by porous flow melting by
Elkins and Lambart [68] have also demonstrated that no single melting scenario developed
to date can generate all the observed global data, even accounting for lithologic hetero-
geneity (see below), and that a range of transport mechanisms and degrees of melt-rock
equilibration are necessary on a global scale. Other studies (e.g., [114,115]) likewise have
predicted from fluid dynamics and two-phase flow constraints that magma flow through a
porous network must be spatially heterogeneous, with varying extents of flow in localized,
high-flux “channels” surrounded by diffuse, low-flux matrices.
### 3.6.3. Effects of Source Heterogeneity on U-Series Disequilibria in MORB

Because $^{230}$Th excesses in MORB support melting in the presence of residual garnet, and because there is a continuum of radiogenic isotopic compositions and trace element abundances that are often correlated with $^{238}$U-$^{230}$Th, decompression melting of a heterogeneous mantle can explain much of the short-range and global compositional variability observed in MORB (e.g., [50,51,79,100,116]).

<figure><img src="image_20.png" /><figcaption>**Figure 16.** $^{230}$Th/$^{238}$U versus $^{226}$Ra/$^{230}$Th for axial 9-10°N samples showing inverse correlation. Modified from Sims et al. [26].</figcaption></figure>

A key petrologic characteristic of decompression melting of a heterogeneous mantle is that different lithologies have distinct melting behavior, including different solidus depths and different melt productivities (e.g., [43,71,117–126]). Recycled crustal material, which may be the principal carrier of isotopically enriched signatures in the mantle (e.g., [117]), likely originates as gabbroic and basaltic oceanic crust that metamorphoses to become eclogite under upper mantle conditions. In addition to the preservation of these eclogites during mantle convection, processes such as physical or mechanical stirring of blobs and veins, and deep partial melting followed by melt-rock reaction with ambient peridotites, can produce a range of broadly pyroxenitic compositions with melting behavior distinct from that of peridotites (e.g., [71,117]). In general, pyroxenites have deeper solidi and narrower melting intervals than peridotite rocks along the same mantle adiabat, making them more fertile and magmatically productive than the dominant mantle peridotites (e.g., [71,117]). Due to their deeper solidi, partial melts of pyroxenite may also be extracted from a physically larger source region than partial melts of peridotite alone. Current estimates suggest that, overall, the convecting mantle contains approximately 5% of such recycled and/or hybrid, pyroxenite lithologies on average, but locally the quantity of pyroxenite in a zone of partially melting mantle may be significantly higher (e.g., [117]).

Because of their overall higher melt productivity, pyroxenite veins in peridotite are also expected to contribute disproportionately to an aggregated total melt relative to their absolute abundance, though the degree to which this is the case depends significantly on the type of pyroxenite (e.g., [51,71,88,116–120,127] and references therein). The disproportionate contribution of partial melts from pyroxenites to the overall magma is also
mitigated by latent heat exchange as the two lithologies melt at different depths and in distinct proportions, though the difference in quantity from each lithology is still significant (e.g., [68,70,117,121,122]). For example, for a decompressing mantle source region containing 5% of an eclogite-like (experimental mixture G2; [77,121]) pyroxenite and 95% average peridotite, estimates suggest that 14–24% of the total melt produced by mass will derive from the pyroxenite.

<figure><img src="image_21.png" /></figure>

**Figure 17.** Results of one- and two-porosity partial melting calculations for ($^{230}$Th/$^{238}$U) and ($^{226}$Ra/$^{230}$Th) for progressive adiabatic decompression melting of a two-layer (spinel- and garnet-bearing) peridotitic mantle as a function of melting depth. After Jull et al. [112] using a constant melting rate and threshold or reference porosity of 0.001, a solid mantle upwelling rate of 1 cm/yr., and bulk solid/melt partition coefficients, all after Jull et al. [112] except where otherwise indicated. (A) Results for chromatographic/equilibrium porous flow melting with a single porous network for liquid transport using the numerical solution to porous flow with variable porosity with depth; after Spiegelman & Elliott [18] and Spiegelman [55]. (B) Results for dynamic melting, considering a solid mantle upwelling rate of 1 cm/yr. and extracted liquid transport rate of 500 cm/yr., as well as a second scenario for upwelling rates of 10 cm/yr. and instantaneous liquid transport (green lines), which replicates dynamic melting conditions; after Stracke et al. [28]. (C) Results for an example two-porosity melting with a high-porosity channel that forms at a depth of 70 km, as described in the text of Jull et al. [112], where $f_1$ indicates the porosities with depth for the low-porosity region, $f_2$ shows the porosities with depth for the high-porosity region, and resulting liquid activity ratios are shown for both liquid zones (high-porosity results are in purple). The model scenario shown assumes that 1% of the total melt flux (S) remains in the low-porosity zone, with 99% entering the high-porosity channel, and a horizontal extent of channelization (C) of 1%.

Beyond the overall liquid mass proportion, because pyroxenite lithologies, perhaps ultimately derived from recycled crustal material, are relatively enriched in incompatible elements, whereas much of the ambient peridotitic rock are trace element-depleted residues of prior partial melting (the depleted MORB mantle source, “DMM”), the pyroxenite partial melts are further expected to contribute an outsized proportion of incompatible elements to resulting mixtures (e.g., [43,119,120,128]).

A systematic treatment of bi-lithologic mantle melting beneath mid-ocean ridges (Figures 19 and 20) shows that, consequently, if a mixed pyroxenite-peridotite source partially melts, then the radiogenic isotopic and trace element composition of a lava formed by melting of that source will be highly biased from that of the average mantle or the dominant peridotite, and is instead a combined, unevenly weighted function of (1) the pyroxenite type and its melting behavior; (2) the proportion of pyroxenite in the initial
unmelted solid; and (3) the combined degrees of melting of both lithologies. Uranium-series disequilibria are further influenced by both the magma transport regime and the effects of mantle temperature on mineralogy in each lithology [68]. Essentially, the isotopic compositions of the erupted basalts thus track differences in the relative proportions of mafic lithologies (i.e., pyroxenite and/or eclogite) and ambient peridotite in their melting mantle source, where lavas with more radiogenic $^{87}$Sr/$^{86}$Sr and $^{206}$Pb/$^{204}$Pb, less radiogenic $^{143}$Nd/$^{144}$Nd and $^{176}$Hf/$^{177}$Hf, and lower U/Th ratios record a melt source containing higher quantities of pyroxenite rocks.

<figure><img src="image_22.png" /></figure>

**Figure 18.** Effects of end-member melting models—dynamic melting [20], chromatographic melting [18], two-porosity/polybaric porous flow [112], and dynamic melting using variable upwelling and melt extraction velocities [28]. The last two models both successfully explain the ($^{226}$Ra/$^{230}$Th) and $^{230}$Th/$^{238}$U) observed in the 9–10° N EPR data [26]. However, only the two-porosity/polybaric porous flow model [112] explains all of the major and trace element data observed in the 9–10° N EPR data and also the melt channels as well as evidence for melt-rock reactions observed in ophiolite complexes by Kelemen et al. [106]. The numbers within the melt columns on the left correspond to the trends in $^{226}$Ra xs vs. $^{230}$Th xs shown on the right. Red dashed lines depict the garnet-spinel phase transition and black dashed lines are equilines.
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

**Figure 19.** Partial melting calculation results after Elkins & Lambart [68], showing (A) ($^{226}$Ra/$^{230}$Th) and (B) ($^{231}$Pa/$^{235}$U) vs. ($^{230}$Th/$^{238}$U) in computed partial melts where mantle Tp = 1300 °C, (C) ($^{226}$Ra/$^{230}$Th), and (D) ($^{231}$Pa/$^{235}$U) vs. ($^{230}$Th/$^{238}$U) compositions for Tp = 1400 °C. Partial melts of G2 pyroxenite are shown in red, and those of coexisting peridotite in blue for a final melting pressure of 5 kbar. Calculated results are shown for G2 pyroxenite abundances in the solid source of 1 to 50% and a range of magma transport models (equilibrium, disequilibrium, and scaled disequilibrium transport with a Damköhler number of 0.1), all for a representative set of results using a solid mantle upwelling of 3 cm/yr. and a residual maximum porosity of 0.8%. Mass-balanced binary melt mixtures (purple symbols) have been computed assuming that all coexisting partial melts are extracted and homogenized for each melting scenario. Global MORB data are shown for comparison (see the caption for Figure 13 for the data sources).

Several studies (e.g., [50,51,105,109]) have found that MORB with relatively enriched isotopic compositions (i.e., relatively radiogenic $^{87}$Sr/$^{86}$Sr, $^{206}$Pb/$^{204}$Pb, $^{207}$Pb/$^{204}$Pb, and $^{208}$Pb/$^{204}$Pb, and unradiogenic $^{143}$Nd/$^{144}$Nd and $^{176}$Hf/$^{177}$Hf) exhibit both (a) elevated ($^{230}$Th/$^{238}$U) and (b) low, age-constrained ($^{231}$Pa/$^{235}$U) and ($^{226}$Ra/$^{230}$Th), when compared to the overall global data set. These are properties that isotopically enriched MORB share with some highly enriched ocean island basalts (e.g., [50,51,105,109]) despite the expected passive upwelling setting for mid-ocean ridges, and in some cases (e.g., the North Atlantic) far lower solid mantle upwelling rates than expected for hotspots. Studies of enriched MORB, as well as systematic modeling of heterogeneous mantle melting using
time-dependent melting models [67,68], have demonstrated that the high (230Th/238U) accompanied by low (226Ra/230Th) and (231Pa/235U) observed in some MORB cannot be produced by peridotite melting, whereas lithologies with higher fusibility and a greater stability range for garnet (i.e., garnet-bearing pyroxenites; e.g., [71,117] and references therein) can produce the observed patterns. Essentially, rapid melting of peridotite alone cannot achieve the same signatures, even with accelerated solid upwelling rates and elevated mantle temperatures, thereby necessitating the presence of at least moderate quantities of mafic (pyroxenite) rocks that melt more rapidly in these settings.

<figure><img src="image_25.png" /></figure>

**Figure 20.** Calculated porous flow partial melting results after Elkins & Lambart [68] for pure peridotite and peridotite–pyroxenite weighted melt mixtures and computed axial depths, along with peridotite melt modeling results for (230Th/238U) and axial depth after Bourdon et al. [34] for comparison. (A) Global MORB data (symbols and data sources are the same as in Figure 7). The melting modeling results from Elkins & Lambart [68] were computed assuming a solid mantle upwelling rate of 3 cm/yr., a residual maximum porosity of 0.8%, and a mantle potential temperatures of (B) 1300 °C and (C) 1400 °C. Symbol colors indicate mass-balanced mixing compositions of the partial melts of a peridotite coexisting with a series of pyroxenite compositions, as well as a pure peridotite source, as labeled. Results are shown for equilibrium porous flow (filled symbols), disequilibrium porous flow (open symbols), and disequilibrium porous flow scaled with a Damköhler number of 0.1 (shaded symbols). Additional results for pure peridotite melting at Tp = 1200 °C are
included in panel (B) for comparison with outcomes using alternative U and Th partition coefficients [69,129] illustrated with hexagon symbols. Crustal thicknesses were determined given the total quantity of melt produced over a 1D melting column for the two melting lithologies, assuming a 15% density contrast between mantle and oceanic crust and full melt extraction to form the crust. The seafloor depth was then computed for a given crustal thickness using a simple isostatic model with 200 km compensation depth, dynamic support of the surrounding plate, and a typical ocean density profile as a reference column, with 7 km crustal thickness, 3 km of ocean water, and underlying layers and densities after Lin et al. [130].

### 3.6.4. U-Series Techniques for Deciphering the Origins of MORB

Below, we outline how U-series disequilibria are a useful tool for deciphering the heterogeneous mantle and melting origins of MORB by exploring a series of case studies in greater detail.

Uranium-series disequilibria in normal MORB (N-MORB) and E-MORB from the East Pacific Rise have helped to assess the likely origins of E-MORB overall. Waters et al. [51] observed that E-MORB and N-MORB from the 9–10° N EPR region form a continuum of compositions that are best explained by mixing of melts derived from a heterogeneous mantle source with at least two major components that are both elementally and isotopically distinct. All young N-MORB and E-MORB collected from within ~5 km of the AST along 9–10° N EPR have significant $^{230}$Th excesses [26,107], most probably indicating the presence of garnet in the mantle source residue and recent fractionation of U from Th compared to the half-life of $^{230}$Th. These observations are at odds with studies of E-MORB that have called upon a two-stage process involving (1) ancient cryptic metasomatism and (2) melting of that cryptically metasomatized spinel peridotite beneath the ridge axis to generate E-MORB compositions [131]. However, it is important to note that metasomatism, recycling of oceanic crust (MORB), and recycling of ocean island basalts (OIB) that produce fertile, garnet-rich lithologies in the resulting source rocks are possible candidates for the origin of the mantle source enrichment of E-MORB. However, the large $^{230}$Th excesses observed in EPR E-MORB are difficult to reconcile with a model that involves large degrees of melting of ancient cryptically metasomatized spinel peridotite as tested by Donnelly et al. [131].

Waters et al. [51] used time-dependent melt modeling to show that mixing of deeply generated garnet pyroxenite partial melts with shallowly generated, accumulated pyroxenite-peridotite melt mixtures in a two-porosity melting and melt transport regime results in excellent matches to both the incompatible element abundances and ratios of E-MORB and the trajectory of the EPR mixing array. In particular, the models produce both elevated highly incompatible element abundances and a shallow MREE/HREE slope, a common feature of E-MORB previously attributed to melting of spinel peridotite [131]. Waters et al. [51] thus proposed that inferred mixing trends among trace elements and long-lived radiogenic isotopes observed in Pacific MORB reflect mixing of melts during melt extraction and crustal storage, rather than large-degree melting of cryptically metasomatized spinel peridotite as proposed by Donnelly et al. [131]. Waters et al. [51] further suggested a model for the genesis of E-MORB by low-degree melting of an enriched, garnet-bearing source that is melting deep and in the deep corners and base of a triangular MORB melting regime.

Taken together, the results of the two-porosity melting mixing model and the non-systematic U-series ages of E-MORB and N-MORB also suggest that E-MORB and N-MORB may be produced concurrently and occupy axial melt chambers (AMC) at the same time. However, along-axis variations in the aspect ratios of axial melt lenses may inhibit along-axis mixing (e.g., [132]) and hence the ability of N-MORB and E-MORB melts to mix. Because even small increases in the proportion of garnet pyroxenite versus peridotite
partial melt contributions can result in the eruption of E-MORB and not N-MORB, the less robust nature of magmatism at 9°30′ N EPR and the western limb of the small overlapping spreading center at 9°37′ N compared with 9°50′ N EPR may decrease melt sill thickness and inhibit along-axis mixing of geochemically variable melts. Thus, along-axis differences in the physical properties of the AMC owing to variations in magma supply may be directly related to the greater abundance of E-MORB erupted near 9°30′ N and 9°37′ N compared with 9°50′ N EPR [51].

At intermediate spreading rates, the compositions of basalts from the Southeast Indian Ridge (SEIR) vary along the ridge axis in both ($^{230}$Th/$^{238}$U) and ($^{226}$Ra/$^{230}$Th) [109]. The data exhibit an inverse correlation between ($^{226}$Ra/$^{230}$Th) and ($^{230}$Th/$^{238}$U) similar to that of the global array [26], but ($^{230}$Th/$^{238}$U) does not correlate in a straightforward way with axial depth or inferred crustal thickness, suggesting additional controls over melt composition and magma generation beyond mantle temperature variations. According to Russo et al. [109], time-dependent modeling suggests that the observed systematics in ($^{226}$Ra/$^{230}$Th) and ($^{230}$Th/$^{238}$U), axial depth, crustal thickness, and Th/U ratios are best explained by an along-axis temperature gradient accompanied by variations in lithologic composition along the axis (e.g., the presence of pyroxenite veins). Specifically, Russo et al. [109] determined that melt supply and crustal thickness decrease systematically from west to east along the ridge axis, but the melting rate required by U-series disequilibria only correlates with the melt supply in the western and central parts of the study area. In the east, the melting rate required to explain the observed data must be high despite the lower overall melt supply to the ridge axis, which suggests the increased presence of fusible pyroxenite veins in an overall cooler regime.

Uranium-series disequilibria have been useful for exploring the role of heterogeneous melting not only at fast and intermediate spreading rates, but also at slow spreading rates, as explored above for the slow- to ultraslow-spreading Arctic ridge system. Along the Kolbeinsey, Mohns, Knipovich, and Gakkel Ridges, slow spreading is distinct from underlying mantle temperature, which is elevated near Iceland, such that spreading and inferred passive upwelling rates can be isolated from mantle temperature effects on the peridotite solidus depth [50,100]. Radiogenic isotopes and trace element compositions also indicate that basalts from the Mohns and Knipovich Ridges are products of melting a heterogeneous source mantle with compositional characteristics significantly different from the mantle beneath the Kolbeinsey Ridge (e.g., [100,133,134]). For U-series disequilibria, the high-temperature but likely peridotite-dominated system beneath the Kolbeinsey Ridge produces elevated ($^{230}$Th/$^{238}$U) and ($^{231}$Pa/$^{235}$U), while lavas from the Mohns and Knipovich Ridges have relatively high ($^{230}$Th/$^{238}$U) but low age-constrained ($^{231}$Pa/$^{235}$U) and ($^{226}$Ra/$^{230}$Th), best explained by melting a pyroxenite-bearing mantle source [50,100].

In another slow-spreading ridge example, MORB from five second-order segments between 5 and 11° S along the southern Mid-Atlantic Ridge (MAR) exhibit a wide range of $^{230}$Th excesses from 0 to 32%, as well as heterogeneous trace element and radiogenic isotope compositions [135] and axial depths. Notably, however, ($^{230}$Th/$^{238}$U) does not correlate with axial depth, indicating that it is not controlled by regional changes in mantle temperature, much like observations along the SEIR [109]. As seen elsewhere, mixing of partial melts from a heterogeneous mantle source is necessary to explain the combination of distinct source compositions and melting rates required to generate U-series disequilibria.

A follow-up study of basalt compositions along the Northern Kolbeinsey Ridge (NKR) and nearby Jan Mayen Island hotspot provides an additional, useful example of the method for exploring magma origins at slow spreading but anomalously high magma production rates. Elkins et al. [105,136] evaluated $^{238}$U-$^{230}$Th-$^{226}$Ra disequilibria, radiogenic isotope variations, and trace element compositions for this volcanically active region adjacent
to the Jan Mayen Fracture Zone and observed high (²³⁰Th/²³⁸U) accompanied by low (²²⁶Ra/²³⁰Th) in age-constrained basalts from both the NKR and Jan Mayen Island, similar to Mohns and Knipovich Ridge lavas farther to the north. The U-series compositions of lavas for the region correlate with radiogenic isotopes, indicating a source heterogeneity effect. Trace element compositions of Jan Mayen Island and NKR lavas are distinct, however, suggesting slightly different underlying mantle compositions that contain different types or quantities of enriched, pyroxenitic material [136]. In particular, the combination of high (²³⁰Th/²³⁸U) with low (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) observed in NKR, Jan Mayen, Mohns, and Knipovich Ridge basalts requires faster melting in the presence of residual garnet than would occur in peridotite. To achieve the observed disequilibria, the presence of a garnet-bearing and highly fusible rock type is necessary, supporting the presence of fertile, isotopically enriched lithologies in the underlying mantle regime [100,105].

The case studies described here demonstrate that the combination of (1) U-series disequilibrium measurements in fresh, age-constrained oceanic lavas and (2) thoughtful use of time-dependent melt models can help to distinguish between unique scenarios for mantle melt generation, such as variations in solid mantle upwelling rate, diverse sampling of melt components from a multiple porosity decompressing regime, and the presence of heterogeneous source lithologies with distinct melting behavior and melt productivities. Specifically, time-dependent modeling of partial melt generation for a range of lithologic types and upwelling rates illustrates distinct differences between the outcomes of rapid melting due to fast mantle upwelling and rapid melting due to the presence of highly fusible material, making it possible to differentiate between those two scenarios [67,68,100,105]. This technique has successfully been applied to investigations of MORB generation in a variety of settings and across all spreading rates. Lastly, it is important to note that the efficiency of focused versus diffuse flow is also a critical consideration [112,137,138].

### 3.7. Ocean Island Basalts (OIB)

On a global basis, ocean island basalts (OIB), like MORB, are characterized by (²³⁰Th/²³⁸U) activity ratios greater than one [1, 14, 15, 25, 139–141]. These ²³⁰Th excesses indicate that the magmas were likely generated by partial melting in the presence of garnet, which is not surprising given the thickness of the lithosphere in intraplate settings. Also, like MORB, decompression melting due to adiabatic mantle upwelling generates OIB; as such, the preferred melting models for investigating magma origins are similar. The main difference is that hotspot settings have higher mantle potential temperatures, buoyancy fluxes, and melting rates, and are argued to be the surface manifestation of deep mantle plumes [25, 98, 142–148]. Because those buoyancy fluxes and melting rates also vary systematically from the center to the edges of upwelling mantle plumes, one of the most important outcomes from U-series measurements in OIB has been the use of ²³⁸U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹Pa disequilibria to affirm the existence of mantle plumes, to determine the melt productivity of a lava's underlying mantle source, and to understand and map the upwelling structure of individual mantle plumes [25, 139–141, 149].

## Using U-Series Disequilibria to Evaluate Upwelling in Mantle Plumes

### Hawaiian Lava Origins: Evidence for a fast-upwelling plume with structure

The Hawaiian Archipelago is one of the most widely studied (and perhaps the best understood) examples of ocean island volcanism. The islands of Hawaii are at the southeastern end of a linear chain of volcanic islands, the 6000 km long Hawaiian-Emperor ridge. These islands are inferred to be the surface expression of a mantle plume having sat beneath the westward-moving Pacific plate for ~75 Myr. Based on the volume and composition of these basalts, as well as geophysical observations such as the geoid height and bathymetric
swell, it is generally thought that the northwest motion of the Pacific plate is moving the main island of Hawaii off the hotspot center and that the new submarine volcano Kama'ehuakanaloa seamount (formerly known as Loihi) to the southeast represents the initial stages of melting produced at the leading edge of the plume [14,25,139].

Numerous studies on lavas from the Hawaiian Islands have shown that over the lifespan of a single volcano, the major and trace element chemistry of the erupted basalts vary significantly and systematically (e.g., [150,151]). This variation is thought to be the result of movement of the Pacific plate over a stationary mantle plume, and is typified by four distinct phases of volcano evolution: (1) a pre-shield stage, during which both tholeiitic and alkaline lavas are erupted; (2) a tholeiitic shield-building stage, which represents about 98% of the volcano's total volume; (3) a post-caldera or post-shield stage during which the lavas become more alkalic in composition and erupt as small parasitic cones, forming a thin veneer over the shield-stage tholeiites; and (4) a late or rejuvenated alkalic volcanism stage, which occurs after an ill-defined period of quiescence. While the time required for the entire evolution of one volcano (ca. 2–5 Ma) is much longer than the time-span that can be investigated using <sup>238</sup>U-<sup>230</sup>Th disequilibria, it is possible to obtain from several different volcanoes a suite of historic or young radiocarbon dated lavas that range in composition from primitive tholeiite to primitive alkali basalt and collectively encompass all four stages of Hawaiian volcano evolution.

Past U-series work on Hawaii has included both early alpha counting studies [152–155] and more recent mass spectrometric studies using TIMS, SIMS, and MC-ICPMS [14,25,44,139,156,157]. Most U-series studies have focused on the active shield-stage tholeiitic volcanoes, Kilauea and Mauna Loa [14,25,139,152–159]. However, <sup>238</sup>U-<sup>230</sup>Th-<sup>226</sup>Ra and <sup>235</sup>U-<sup>231</sup>Pa disequilibria have also been measured in young tholeiitic and alkaline lavas from the preshield-stage submarine volcano Kama'ehuakanaloa [14,25,160] and in alkaline lavas from the post-shield stage volcanoes, Hualalai, and Mauna Kea, and the rejuvenated stage volcano Haleakala [14,25,44]. <sup>238</sup>U-<sup>230</sup>Th disequilibria further have been measured in young alkaline submarine flows from the South Arch [14].

Early U- and Th-series measurements that used radioactive counting techniques demonstrated the existence of significant radioactive disequilibria in lavas from Kilauea, suggesting that chemical fractionation between <sup>238</sup>U-<sup>230</sup>Th and <sup>230</sup>Th-<sup>226</sup>Ra occurred during the formation of Kilauea magmas on timescales < 375 kyr and 8 kyr, respectively. However, the origin of this fractionation, whether melting or secondary processes such as assimilation of altered oceanic crust, was highly uncertain [152]. Subsequent mass spectrometric studies with higher-precision measurements have consistently demonstrated that the <sup>230</sup>Th excesses in tholeiitic lavas from Kama'ehuakanaloa, Kilauea, and Mauna Loa are small (2 to 6%) [14,25,139] and constant [157] compared to, e.g., the large <sup>230</sup>Th excesses seen in MORB (up to 25%), even though MORB are also mostly tholeiites. Note that Cohen et al. [161] showed much larger extents of <sup>230</sup>Th excesses for Mauna Loa, but as shown by Sims et al. [25], U-series systematics in these Mauna Loa lavas have been changed by secondary alteration.

In contrast, young alkalic lavas from the post-shield stage volcanoes, Hualalai and Mauna Kea, and the active rejuvenated stage volcano, Haleakala, have much larger <sup>230</sup>Th excesses than Hawaiian tholeiites, ranging from 15 to 30% [14,25,44]. The observation that these <sup>230</sup>Th excesses vary systematically across the range of compositions observed for Hawaiian lavas, with the smallest excesses observed in the tholeiites and the largest excesses seen in basanitic lavas from Haleakala, demonstrates that (<sup>230</sup>Th/<sup>238</sup>U) in Hawaiian lavas is correlated with the overall extent of melting in the mantle source inferred from the major and trace element chemistry of the lavas [14]. This correlation provided definitive evidence that U-Th fractionation is related to the melting process, at a time when the
origins of $(^{230}_{Th}/^{238}U)$ disequilibria in magmatic systems, particularly in MORB, were debated. Additionally, this conclusion was contemporaneous with early partitioning studies [58,59,61] showing that $(^{230}_{Th}/^{238}U)$ excesses can be produced by melting in the presence of garnet. And, unlike MORB, whose initial depth of melting has been debated (e.g., [162]), the conclusion that OIB are generated by melting of garnet-bearing mantle source rocks makes sense because the thickness of the oceanic lithosphere beneath ocean islands and the likely higher mantle temperature both require higher average depths of melting.

In contrast to $(^{230}_{Th}/^{238}U)$, $(^{226}_{Ra}/^{230}_{Th})$ disequilibria in Hawaiian shield stage tholeiites are much larger, with $^{226}_{Ra}$ excesses ranging from 7 to 15% [25,139,155,156,159]. For alkalic lavas from Hualalai and Haleakala, $^{226}_{Ra}$ excesses are even greater, up to 30% [25,44]. In Hawaiian lavas, $(^{231}_{Pa}/^{235}_{U})$ suggests even greater extents of parent-daughter fractionation, with $^{231}_{Pa}$ excesses ranging from 10 to 15% for tholeiitic lavas from Kilauea, Mauna Loa, and Kama'ehuakanaloa, and up to 78% for the basanitic lavas from Haleakalā [25,163]. Like $(^{230}_{Th}/^{238}U)$, the magnitude of the $^{235}_{U}$-$^{231}_{Pa}$ disequilibria correlates with the major and trace element chemistry of the lavas, and, by inference, the fraction of melting in their mantle source and/or the isotopic enrichment of the source.

While it can be shown that $^{238}_{U}$-$^{230}_{Th}$ and $^{235}_{U}$-$^{231}_{Pa}$ disequilibria in Hawaiian lavas are mainly sensitive to melt fraction because of the long half-lives of $^{230}_{Th}$ and $^{231}_{Pa}$ [14,15,25], interpretation of $^{226}_{Ra}$, because of its much shorter half-life, requires dynamic and chromatographic melting models that consider the timescales of melt generation and melt extraction (e.g., [18,20]). Sims et al. [25] showed that $^{238}_{U}$-$^{230}_{Th}$-$^{226}_{Ra}$ and $^{235}_{U}$-$^{231}_{Pa}$ disequilibria in Hawaiian lavas coupled with forward time-dependent models (dynamic and porous flow) can provide model-constrained estimates of porosity, melt transport times, and solid mantle upwelling rate within the melting Hawaiian plume. As discussed earlier, during ingrowth melting processes $^{238}_{U}$-$^{230}_{Th}$-$^{226}_{Ra}$ and $^{235}_{U}$-$^{231}_{Pa}$ disequilibria are uniquely sensitive to melting rates and porosities.

This variability in the solid mantle upwelling structure of the Hawaiian plume from the center to its leading and trailing edges was critical to further solidify the existence of mantle plumes (Figure 21). While there had previously been considerable evidence for the existence of a deep-rooted mantle plume beneath Hawaii, few studies constrained the nature of the melting region and the dynamics of the plume from both geophysical and geochemical observations. Watson and McKenzie [146] used an axisymmetric plume model to calculate the melt production rate and bulk major and trace element chemistry of magmas, but they did not account for the spatial and compositional variability seen in the chemistry of erupted lavas, nor did they consider the possible effects of the overriding Pacific plate on underlying plume dynamics. Richardson and McKenzie [24] used the plume models developed by Watson and McKenzie [146] to calculate the effects of plume upwelling on U-series disequilibria in partial melts and infer the melting rate (which, in turn, is related to the solid mantle upwelling rate) and porosity of the melting region. However, because there was very little U-series data from volcanoes other than Kilauea and Mauna Loa, they could not correlate spatial variations in the U-series disequilibria along the plume track with the dynamics of the plume. With a larger data set available, Sims et al. [25] mapped the spatial variation of U-series data in young lavas from the center of the inferred Hawaiian plume to its fringes, demonstrating that solid mantle upwelling rates are highest in the plume center and lowest at its leading and trailing edges (Figure 21). The calculated U-series upwelling rates are comparable to upwelling velocities calculated from the axisymmetric plume models of Watson and McKenzie [146] and Hauri et al. [98], but the U-series upwelling rates also suggest an asymmetry not predicted by prior plume models. In particular, the interpretation that solid mantle upwelling is slightly greater
beneath Kama'ehuakanaloa (the leading edge) than Hualalai and Haleakala on the trailing side appears more consistent with the plume model of Ribe and Christensen [164], which considers the influence of the overriding Pacific Plate on the plume's upwelling dynamics. Later studies [44,139,165] support the interpretations of Sims et al. [25] with further measurements of Kilauea, Kama'ehuakanaloa, and Haleakala.

<figure><img src="image_26.png" /></figure>

**Figure 21.** Solid mantle upwelling rate calculated from $^{238}$U-$^{230}$Th-$^{226}$Ra and $^{235}$U-$^{231}$Pa disequilibria versus the radial position of the volcanoes relative to the center of the Hawaiian plume. The center of the plume (inset) is located so that the active volcanoes Loihi, Kilauea, and Mauna Loa are positioned appropriately for their model ages and eruption rates [166]. Also shown are the model curves for solid mantle upwelling from Watson and McKenzie [146], Hauri et al. [98], Ribe and Christenson [164] and Jull and Ribe [167]. The model of Watson and McKenzie [146] assumes a constant viscosity across the plume, whereas the Hauri et al. model incorporates a temperature-dependent viscosity. The upwelling rates from the plume model of Jull and Ribe [167] consider the effect of motion of the Pacific plate and depletion buoyancy on the dynamics of the plume, whereas the calculations of Hauri et al. [98] and Watson and McKenzie [146] do not. The $^{230}$Th and $^{231}$Pa data show a clear asymmetry in upwelling velocity, which is evident in the results of Jull and Ribe [167]. These estimates will need to be refined with more accurate values for distribution coefficients and better models for melt productivity, but they nevertheless illustrate that U-series disequilibrium isotopic data may provide a much needed constraint on the dynamics of plumes. The data from Pietruszka et al., [165] for Kama'ehuakanaloa shown on the plot are not differentiated by alkalinity as they were in previous studies. The images on the right show the plume's (a) potential temperature; (b) depletion (i.e., removal of denser minerals from the residue); and (c) melting rate along the plume center from the reference model of Ribe and Christenson [164]. See map inset for the location of the A-A' section line in (a-c). Modified from [25,164,167].

Subsequent work by Kokfelt et al. [140] and Bourdon et al. [149] showed similar variations in the activities of ($^{230}$Th/$^{238}$U) and ($^{231}$Pa/$^{235}$U) when plotted as a function of the distance from the centers of the Azores and Iceland hotspots. However, these later studies do not consider the type of variations in upwelling that result from plume-lithospheric interaction that has been suggested for Hawaii [1,25].
## U-series disequilibria to understand global hotspot buoyancy flux

Chabaux and Allègre [168] showed a broad correlation between (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰Th), which they argued was due to variations in the degree of mantle melting. However, this correlation can also be interpreted as a product of variations in hot spot buoyancy flux as estimated by Sleep [1, 35, 149, 169].

Bourdon et al. [149] demonstrated an apparent link between the buoyancy flux (B) of a hotspot (defined as the product of the speed, the cross-sectional area, and the density anomaly of the upwelling) and the upwelling velocity determined from time-dependent melting models when compared with measured (²³⁰Th/²³⁸U) and (²³¹Pa/²³⁵U). From the data and calculations of Bourdon et al. [149], the range of buoyancy fluxes (B) found in mantle plumes is 0.5 to 1 Mg s⁻¹ for hotspots like Afar and Iceland and up to between 3.5 and 9 Mg s⁻¹ for Hawaii, depending on how the buoyancy flux is calculated (Ribe [170] or Sleep [169], respectively). The range in (²³⁰Th/²³⁸U) varies from ~1.4 (Afar) to ~1.08 (Hawaii) and the range of ²³¹Pa/²³⁵U varies from ~1.4 (Afar) to ~1.08 (Hawaii) [149].

However, it should be noted that the Samoan data presented in Sims and Hart [171] and Sims et al. [42] show that for age-constrained samples with ²²⁶Ra excesses, (²³⁰Th/²³⁸U) = 1.02–1.13 and (²³¹Pa/²³⁵U) = 1.05–1.33 are both low given the hotspot’s estimated buoyancy flux [169,170]. These Samoan data change the apparent form of the Bourdon et al. [149] correlation, making it much more hyperbolic than previously estimated and implying a higher excess potential temperature anomaly in the plume-generating boundary layer. Alternatively, the buoyancy flux for Samoa may be significantly underestimated, as Samoa sits adjacent to a complicated tectonic regime, namely the Tonga subduction zone [42,171].

Two important shortcomings of this approach, namely making global comparisons of U-series disequilibria in hotspot lavas with calculated hotspot buoyancy flux, are (1) that because they do not integrate across the whole plume [169,170] the calculated buoyancy fluxes could be significantly biased; and (2) these generalized comparisons do not consider the local influence of the motion of the overriding plate [25,170]. For example, in hotspots like Hawaii, where the motion of the overriding Pacific plate influences the upwelling rate [1, 25, 169], the extent of U-series disequilibria in young lavas varies from Kilauea to Haleakala, with (²³⁰Th/²³⁸U) ranging from 1.02 (Kilauea) to 1.3 (Haleakala) and (²³¹Pa/²³⁵U) from 1.3 (Kilauea) to 2.2 (Haleakala). This spans the entire range of (²³⁰Th/²³⁸U) for all plumes across the entire global correlation shown in Bourdon et al. [149]. While it is important to note that this range in Hawaiian (²³⁰Th/²³⁸U) disequilibria is considered in the error bars shown for Hawaii, such simple 1D calculations are necessarily only first-order approximations and fraught with uncertainty. A full understanding of how U-series measurements are related to hotspot buoyancy flux requires integrating those measurements across a whole plume section and then incorporating the influence of the motion of the overriding plate on the underlying upwelling plume [164], and also accounting for complicating factors such as source heterogeneity and melt transport times.

## 3.8. Global Correlations of U-Th Disequilibria with Radiogenic Isotopes in OIB and MORB

An important advance in our understanding of the Earth's mantle was the recognition that mid-ocean ridge and ocean island basalts have distinct radiogenic isotopic ratios that result from partial melting of different mantle sources, which mix to form a mantle array, consisting of a depleted end-member and an enriched end-member [172]. As the science of mantle isotope dynamics developed, both through the addition of more isotope systems and the measurement of many more samples across the ocean basins, the number of required enriched end-members has increased significantly, leading to multiple proposed mixing end-members (DMM, HIMU, EM1 and EM2), primarily distinguished by differences in the
three Pb isotope systems ($^{206}$Pb/$^{204}$Pb; $^{207}$Pb/$^{204}$Pb and $^{208}$Pb/$^{204}$Pb) and producing what is now known as the mantle tetrahedron (e.g., [173–175] and references therein). In the decades since, extensive measurements of radiogenic isotope and trace element compositions in hotspot lavas have led to a very large, multi-dimensional global data set, leading to an expanded determination of end-member compositions as greater statistical distributions have been achieved. Likewise, updated models for magma generation, transport, and mixing in a complex mantle regime have produced a more nuanced understanding of what information oceanic lavas record (e.g., [176,177]), and more sophisticated statistical treatments have explored what structures are contained in the multi-dimensional global data set [178].

While radiogenic isotope compositions in lavas have been pivotal in defining our understanding of global mantle heterogeneity, they only track ancient mantle processes. Radiogenic isotope tracers of Sr, Nd, Hf, and Pb have parent half-lives of 0.7–106 Gyr and respond to changes in parent/daughter ratios (Rb/Sr, Sm/Nd, Lu/Hf, U/Pb, Th/Pb) over timescales comparable to mantle convection times. In contrast, the $^{238}$U-$^{230}$Th isotopic system is intrinsically different in two important ways: first, the half-life of daughter isotope $^{230}$Th is only 75 kyr; secondly, the abundance of $^{230}$Th relative to parent $^{238}$U is determined by the radioactive principles governing U decay-series systematics (Section 2 above). As a result, after ~375 kyr, the Th isotopic composition of a mantle source for partial melts directly reflects its $^{238}$U/$^{232}$Th source ratio. Thus, comparison of ($^{230}$Th/$^{232}$Th), ($^{238}$U/$^{232}$Th), and ($^{230}$Th/$^{238}$U) with longer-lived radiogenic isotopes (e.g., $^{87}$Sr/$^{86}$Sr) and their record of longer-lived processes can provide fundamental insights into both the nature and evolution of Earth's mantle sources and the modes of melt generation and magma transport.

Some of the earliest comparative studies of U-series disequilibria and radiogenic isotopes in oceanic settings were made by Allègre and Condomines [13], who compared a small number of U-series measurements to $^{87}$Sr/$^{86}$Sr isotopes in lavas from a limited number of MORB and OIB, using counting techniques. They qualitatively argued that there is a much better inverse linear correlation between ($^{230}$Th/$^{232}$Th) and $^{87}$Sr/$^{86}$Sr than between ($^{238}$U/$^{232}$Th) and $^{87}$Sr/$^{86}$Sr. They concluded that ($^{230}$Th/$^{232}$Th) in a lava better represents the original mantle source ($^{238}$U/$^{232}$Th) composition compared to ($^{238}$U/$^{232}$Th), which they posited was fractionated by processes in which Th is more incompatible than U. As such, they suggested that ($^{238}$U/$^{230}$Th) fractionation was due to time-independent melting, and that all fractionation on a U-Th isochron diagram was due to a horizontal displacement without any change in ($^{230}$Th/$^{232}$Th).

Later, using non-linear statistical regression techniques and a much larger oceanic basalt database that included EM2 and HIMU (but not EM1), Sims and Hart [171] examined the relationship between ($^{230}$Th/$^{232}$Th), ($^{238}$U/$^{232}$Th), and the long-lived isotope ratios of Sr, Nd, and Pb, all measured by mass spectrometric methods (Figure 22). Sims and Hart [171] quantitatively demonstrated that: (1) the relationships between ($^{230}$Th/$^{232}$Th), ($^{238}$U/$^{232}$Th), and the long-lived isotope ratios of $^{87}$Sr/$^{86}$Sr and $^{143}$Nd/$^{144}$Nd can be approximated by two-component mixing; (2) ($^{230}$Th/$^{232}$Th) is somewhat better correlated with $^{87}$Sr/$^{86}$Sr and $^{143}$Nd/$^{144}$Nd than $^{238}$U/$^{232}$Th; and (3) the global basalt arrays of ($^{230}$Th/$^{232}$Th) and ($^{238}$U/$^{232}$Th) compared with Pb isotopic compositions require four mixing components, which correspond to the end-members that define the mantle tetrahedron (DMM, HIMU, EM1 and EM2).

### 3.9. Average Depths of Mantle Melting from ($^{231}$Pa/$^{235}$U) and ($^{230}$Th/$^{238}$U) in MORB and OIB

Correlations between ($^{231}$Pa/$^{235}$U) and ($^{230}$Th/$^{238}$U) in MORB and OIB can provide important constraints on the mineralogy and melting behavior of the melting mantle
source. While there have been relatively few U-series studies incorporating $^{235}$U-$^{231}$Pa systematics, considerable information can be gained by incorporating measurements of $^{235}$U-$^{231}$Pa [25,35,141,149,163,179]. Pickett and Murrell [163] observed that in plots of ($^{231}$Pa/$^{235}$U) versus ($^{230}$Th/$^{238}$U), ocean island basalts define a linear trend that is clearly distinct from the MORB trend. These two trends, which persist with a larger database, are best explained by variable clinopyroxene/garnet ratios in the mantle source, which, in turn, has been interpreted to reflect variation in the initial pressure of melting [25,35,141,149]. Using the partitioning data of Landwehr et al. [63] as a function of pressure with simple batch melting models, the distinction between MORB and OIB (Figure 23) can be explained by a variable clinopyroxene/garnet ratio in the mantle source, which for a peridotitic mantle is determined by the initial pressure of melting [1,35,42]. The OIB data lie along melting trends that are consistent with melting at greater pressures (well below the garnet-spinel transition) with a lower clinopyroxene/garnet ratio, while MORB can be produced by a shallower average depth of melting and are characterized by larger $^{231}$Pa excess at a given ($^{230}$Th/$^{238}$U).

<figure><img src="image_27.png" /><figcaption>**Figure 22.** $^{230}$Th/$^{232}$Th versus $^{87}$Sr/$^{86}$Sr diagram for mid-ocean ridge, ocean island, and continental basalts compiled from the literature [14,25,26,32–53,89]. This new data set based on mass spectrometry data shows a hyperbolic rather than a linear trend. If this is interpreted as reflecting a mixing line, it indicates that the Th/Sr ratio in the enriched components (low $^{230}$Th/$^{232}$Th ratios) is greater than the Th/Sr ratio in the depleted component.</figcaption></figure>

In a follow-up assessment of the correlations observed in plots of ($^{231}$Pa/$^{235}$U) versus ($^{230}$Th/$^{238}$U) in OIB, and to assess the mechanisms that control $^{238}$U-$^{230}$Th and $^{235}$U-$^{231}$Pa in lavas, Prytulak and Elliot [141] measured $^{238}$U-$^{230}$Th and $^{235}$U-$^{231}$Pa in lavas from Pico volcano in the Azores. They argued that the young mafic lavas of Pico are products of smaller degrees of mantle melting than MORB, but that due to Pico's proximity to the Mid-Atlantic Ridge and the thin overlying oceanic lithosphere, their degree of melting was large compared to many other OIB. Albeit tightly clustered, data from the basaltic Pico lavas have surprisingly high extents of disequilibria (($^{230}$Th/$^{238}$U) = 1.18–1.27 and
(²³¹Pa/²³⁵U) = 1.477–1.504) compared to both young Hawaiian lavas, which as described above include tholeiitic basalts, alkali olivine basalts, and basanites (with (²³⁰Th/²³⁸U) = 1.02–1.31 and (²³¹Pa/²³⁵U) = 1.10–1.78; [25]), and Samoa alkaline basalts ((²³⁰Th/²³⁸U) = 1.03–1.16 and (²³¹Pa/²³⁵U) = 1.05–1.33; [42]). Thus, from these plots it was reasonable to conclude that Pico lavas originate from a mantle source with a low clinopyroxene/garnet ratio and an overall deeper pressure of melting [141]. However, simple batch melting models predict an extent of melting that is unreasonably small (~0.3%). Prytulak and Elliott [141] instead argued that only a fraction of the disequilibria observed in Pico lavas can be produced by net elemental fractionation, and that ingrowth of the daughter nuclides during melting is required to produce the observed ²³⁸U-²³⁰Th and ²³⁵U-²³¹Pa disequilibria.

<figure><img src="image_28.png" /></figure>

**Figure 23.** (²³⁰Th/²³⁸U) versus (A) (²³¹Pa/²³⁵U) and (B) (²²⁶Ra/²³⁰Th) for a compilation of OIB and MORB data compared with simple batch melting models. The D<sub>U</sub>, D<sub>Th</sub>, and D<sub>Pa</sub> values and respective color schemes for the different lines are the same as shown in Figure 20, except for D<sub>Ra</sub>, which is held constant and assumed to be 4E-5 in all cases. Also shown in (A), are the recent partitioning data of Landwehr et al. [63] as a function of pressure (light grey curves). As discussed in Bourdon and Sims [1], the distinction between MORB and OIB can be explained by a variable cpx/garnet ratio in the mantle source, determined by the initial pressure of melting. The OIB data lie along melting trends that are consistent with melting at a pressure above the garnet–spinel transition and having a low cpx/garnet ratio. In contrast, MORB can be produced by a shallower average depth of melting
and are characterized by larger $^{231}$Pa excesses at a given $^{230}$Th/$^{238}$U. In (B), the $^{226}$Ra excesses for all OIB are smaller than MORB and also tend to be positively correlated. For MORB, the only data set whose age are known to be young relative to the half-life of $^{226}$Ra are from 9 to 10° N EPR [26]. These MORB samples show an inverse correlation between $^{230}$Th excess, $^{226}$Ra excess (as well as major and trace element compositions) which was interpreted as mixing of two types of melt: one having a high-$^{230}$Th component equilibrated with deep, undepleted garnet peridotites (>70 kms); the other having a high-$^{226}$Ra component preserving characteristics indicative of equilibration with highly depleted, residual harzburgite in the uppermost spinel-facies mantle (<70 kms) [112].

### 3.10. A Final Comment on Time-Independent Melting Versus Time-Dependent Melting in MORB and OIB

As discussed above, there are two major competing mechanisms for explaining an excess of daughter over parent isotope in a U-series decay chain: (1) time-independent melting models, which rely entirely on elemental fractionation of a parent-daughter nuclide pair; and (2) time-dependent melting models in which 'ingrowth' of daughter nuclides provides an alternative method to producing U-series excesses [18,20]. In all cases, MORB, because of their low upwelling rates during melting, require daughter nuclide ingrowth to produce $^{238}$U-$^{230}$Th, $^{235}$U-$^{231}$Pa, and $^{230}$Th-$^{226}$Ra disequilibria. However, in OIB, there are cases, such as the main stage of Hawaii, where the inferred upwelling rate is very fast compared to daughter nuclide ingrowth of $^{230}$Th and $^{231}$Pa and where it is highly likely that elemental fraction of Th/U and Pa/U is the dominant mechanism for $^{238}$U-$^{230}$Th and $^{235}$U-$^{231}$Pa disequilibria. In contrast, there are also hotspots, such as Pico in the Azores, where rates of upwelling are small enough that ingrowth of daughter nuclides is required to produce the observed $^{238}$U-$^{230}$Th and $^{235}$U-$^{231}$Pa disequilibria [141]. In yet other places, such as Samoa, the Galapagos, and Pitcairn, the generation of disequilibria appears to lie along a continuum, where both ingrowth and elemental fractionation of U/Th and U/Pa contribute to measured ($^{238}$U/$^{230}$Th) and ($^{235}$U/$^{231}$Pa) in lavas. As such, it is critical that, before interpreting any U-series disequilibria in lavas, one examines the data in the context of other major and trace element information, as well as the geological and tectonic setting [15]. In contrast, because of the much shorter half-life of $^{226}$Ra, ingrowth is necessary to produce $^{230}$Th-$^{226}$Ra disequilibria in all cases.

## 4. Principles and Applications of U and Th Decay Series for Dating Eruptions of Unknown Age Lavas

Because of their wide-ranging half-lives, U and Th decay series nuclides offer an important set of tools for dating volcanic activity on timescales from weeks up to ~375 kyr. The fundamental premise of the U- and Th-series geochronometers is that chemical processes fractionate parent-daughter pairs within a decay series. This parent-daughter system returns to a state of 'secular' equilibrium (i.e., $\lambda_1 N_1 = \lambda_2 N_2$) over a period directly proportional to the daughter isotope's half-life (Figure 24). The typical range of disequilibria depends on the half-life of the daughter isotope and is limited by measurement uncertainties, but in most cases the determinable age limit on these techniques is about five times the half-life of the daughter isotope.

The simplest geochemical scenario is when the chemical fractionation of the parent-daughter ratio is instantaneous. That is, the event occurs over a short timespan relative to the half-life of the parent-daughter isotope system being considered. Other, more complex situations include continuous or episodic fractionation; deconvolving these types of multifaceted scenarios introduces additional variables and uncertainty. In many geological circumstances, the elemental fractionation creating the disequilibrium between parent-daughter pairs can reasonably be considered to be instantaneous relative to the half-life of the daughter isotope (e.g., $^{230}$Th, $t_{1/2}$ = ~75 kyr). In contrast, in shorter-lived
parent-daughter systems, the chemical fractionation process might be occurring over a
timeframe that is either similar to (e.g., <sup>226</sup>Ra, t<sub>1/2</sub> = 1.6 kyr), or even longer than, five times
the isotope's half-life (e.g., <sup>228</sup>Ra, t<sub>1/2</sub> = 5.77 a, or <sup>222</sup>Rn, t<sub>1/2</sub> = 3.85 d).

<figure><img src="image_29.png" /></figure>

**Figure 24.** Schematic of decay of <sup>230</sup>Th, <sup>226</sup>Ra, and <sup>210</sup>Pb in the melt and ingrowth in the mantle residue after a fractionation event occurring at t<sub>1</sub>. Solid trendlines are for an excess in the daughter isotope at t<sub>1</sub>, while dashed trendlines are for a deficit in the daughter isotope at t<sub>1</sub>. Inset shows the short-lived excess of <sup>210</sup>Pb. Approximate valid range of U-Th and Th-Ra dating methods are shown at the bottom in kyr.

Elemental fractionation between a parent-daughter activity ratio is required to start a chronometer. Then, once elemental fractionation for the parent-daughter system ceases, that chronometer/activity ratio returns to equilibrium at a rate proportional to the various progenies' half-lives. The half-life of the longest-lived daughter of a decay chain is the ultimate determinant of when full equilibrium is attained throughout the decay chain. Thus, this final state of equilibrium occurs after ~1.5 Myr for the <sup>238</sup>U-series, (based on five times the half-life of <sup>234</sup>U); ~200 kyr for the <sup>235</sup>U-series (based on five times the half-life of <sup>231</sup>Pa); and ~35 yr for the <sup>232</sup>Th-series (based on five times the half-life of <sup>228</sup>Ra). Once in equilibrium, the decay series chronometers cease to provide temporal information. In this regard, it is important to note that once the whole decay chain is in radioactive equilibrium, the final daughter Pb isotope in-grows at a rate proportional to the half-life of its ultimate U or Th parent. As such, the initial fractionations that occurred among the decay chain nuclides have no subsequent effect on the results of <sup>238</sup>U-<sup>206</sup>Pb, <sup>235</sup>U-<sup>207</sup>Pb, and <sup>232</sup>Th-<sup>208</sup>Pb chronometers when applied to timescales on the orders of millions to billions of years (i.e., deep time). Note, however, that significant initial disequilibrium between <sup>238</sup>U, <sup>234</sup>U, and <sup>230</sup>Th or <sup>235</sup>U and <sup>231</sup>Pa can measurably affect ages for geochronometers involving <sup>206</sup>Pb or <sup>207</sup>Pb near the lower limit of their useful age ranges (e.g., [180]).
Igneous processes efficiently fractionate daughter nuclides from longer-lived parents in the U and Th decay series. This elemental fractionation, coupled with the U and Th decay series systematics, provides a wide array of U and Th decay series chronometers with meaningful age resolution for igneous systems (Figures 2 and 24).

Because of the short half-lives of the daughter isotopes, it is reasonable to assume that, prior to partial melting, U-series isotopes in a lava's mantle or crustal source were in secular equilibrium (i.e., the activities of the parent and daughter isotopes were equal). Partial melting chemically fractionates parent-daughter elements, creating isotopic disequilibria. These disequilibria subsequently decay back to equilibrium at a rate proportional to the half-lives of the daughter isotopes. Later magmatic processes such as crystallization, assimilation, and magma degassing can also induce parent-daughter elemental fractionation, which can augment or overprint the melting-induced disequilibria.

Importantly, regardless of how the parent-daughter disequilibrium is created during a lava's genesis and evolution, and in the absence of ongoing secondary chemical alteration, once the lava has erupted and solidified, fractionation of the parent-daughter nuclides ceases (Figure 24).

The most common parent-daughter pairs used for determining the ages of volcanic rocks are within the $^{238}$U decay chain, namely $^{238}$U-$^{230}$Th, $^{230}$Th-$^{226}$Ra, $^{226}$Ra-$^{210}$Pb, and $^{210}$Pb-$^{210}$Pb. A few studies have used $^{235}$U-$^{231}$Pa, and even fewer have used $^{232}$Th-$^{228}$Ra or $^{231}$Pa-$^{227}$Ac (see Sims et al. [88] as the only study that considers all of these parent-daughter nuclide pairs in igneous systems). A summary of pertinent U-series nuclides, including relative atomic abundances, half-lives, and decay constants, is provided in Table 1.

Dating of crystallization or eruption of igneous rocks using measurements of disequilibria among parent-daughter pairs in the U and Th decay series utilizes four different methods: (1) the presence or absence of U-series disequilibria to establish age limits based on the half-lives of the daughter isotopes; (2) an internal isochron method, commonly obtained on mineral separates; (3) an external isochron method, or U- and Th- series model ages, which measure the change in the excesses or deficits of the daughter isotopes relative to an assumed or known initial value; and (4) in the case of continuously recharging and degassing magma systems, modeling the parent-daughter disequilibrium using multiple U and Th decay series progeny.

It is important to note that U-series dating provides the best age constraints for dating young submarine basalts. Accurate dating of young (<375 kyr) ocean floor volcanism is an important but problematic issue, with numerous implications for topics ranging from understanding crustal construction at mid-ocean ridges to understanding how the various geological, hydrothermal, and biological aspects of the ridge system relate to one another. However, the dating of young oceanic basalts is not trivial; the half-lives of isotopic systems traditionally used for dating basalts, such as $^{147}$Sm-$^{143}$Nd, $^{87}$Rb-$^{87}$Sr, $^{238}$U-$^{206}$Pb, and $^{176}$Lu-$^{176}$Hf, are much too long; K(Ar)-Ar methods suffer from the low abundance of K in most basalts, K uptake from seawater, and/or open system behavior of Ar [181]; and cosmogenic nuclides are not applicable because of cosmic ray attenuation in the water column. Because of the importance of dating young mid-ocean ridge volcanism to the understanding of oceanic crustal construction and a lack of other viable chronometric options, there was considerable early use of U-series measurements to date oceanic basalts (e.g., [182]). Most notable and transformative was the development of mass spectrometric methods for measuring Th isotopes in basalts [183–185]. Because of these concerted efforts, U-series dating of ocean floor volcanism is now well established thanks to an important body of work that includes [26,40,42,48,51–53,80,89,107,108,110,182,183,185–192].
## 4.1. Establishing Lava Eruption Age Limits Based on the Presence or Absence of U-Series Disequilibria

The presence or absence of U-series disequilibria can place ‘first-order’ limits on the age of a lava sample, assuming that the disequilibrium was magmatic in origin and not a result of secondary processes.

In the case of mantle melting such as MORB or OIB, the timescale of mantle circulation is on the order of hundreds of millions to billions of years and it can therefore safely be assumed that the mantle source starts with U-series nuclides in a state of radioactive equilibrium. During petrogenetic processes, most significantly during partial melting, U-series disequilibrium is created by fractionating parent-daughter nuclide pairs. In the absence of secondary processes (e.g., post-eruptive alteration), any disequilibria generated by magmatic processes are “locked in” once the lava has erupted and solidified. The activity ratio then acts as a “stop-watch”, with parent and daughter returning to secular equilibrium after about five half-lives of the daughter isotope in question. Consequently, disequilibria between a parent-daughter pair provides a maximum eruption age. Since $^{238}\text{U-}^{230}\text{Th-}^{226}\text{Ra-}^{210}\text{Pb}$ have half-lives ranging from ~20 yr to 75 kyr, these techniques are appropriate for dating MORB erupted from 0.1–375 kyr (see Table 4 for list of age limits).

Most, but not all, young lavas of known age exhibit disequilibria amongst parent-daughter ratios. Thus, the absence of disequilibria can imply a lower age limit. However, caution with this approach is warranted, because a lava may never have had disequilibria when it erupted (e.g., [193]). Nonetheless, when multiple U-series age determinants are measured on numerous samples of similar composition from the same volcanic field, comparing the presence or absence of different disequilibria can be used to reasonably infer a sample’s age range. For example, suppose one lava sample from a volcano has significant $^{238}\text{U-}^{230}\text{Th}$, $^{235}\text{U-}^{231}\text{Pa}$, and $^{230}\text{Th-}^{226}\text{Ra}$ disequilibria; in that case, $^{230}\text{Th-}^{226}\text{Ra}$ isotopic systematics are most constraining, because $^{226}\text{Ra}$ has the shortest half-life, and require the lava to be less than 8 kyr, providing it has remained a closed system since eruption. Suppose another sample of similar composition from the same field has $^{238}\text{U-}^{230}\text{Th}$ and $^{235}\text{U-}^{231}\text{Pa}$ disequilibria, but no $^{230}\text{Th-}^{226}\text{Ra}$ disequilibrium; in that case, the sample is < 160 kyr and is ‘likely’ older than 8 kyr. Finally, all samples with no $^{238}\text{U-}^{230}\text{Th}$ disequilibria are ‘likely’ > 375 kyr, and so forth (Table 4). Examples where multiple nuclides are used in concert to constrain sample ages include Goldstien et al. [185], Sims et al. [42,80], and Waters et al. [52].

<table><caption>**Table 4.** Age limits of U-series radionuclides.</caption><thead><tr><td>Nuclide Pair</td><td>Age Constraint</td></tr></thead><tbody><tr><td><sup>210</sup>Po/<sup>210</sup>Pb ≠ 1</td><td>&lt;2 yr</td></tr><tr><td><sup>210</sup>Pb/<sup>226</sup>Ra ≠ 1</td><td>&lt;100 yr</td></tr><tr><td><sup>226</sup>Ra/<sup>230</sup>Th ≠ 1</td><td>&lt;8 kyr</td></tr><tr><td><sup>231</sup>Pa/<sup>235</sup>U ≠ 1</td><td>&lt;160 kyr</td></tr><tr><td><sup>230</sup>Th/<sup>238</sup>U ≠ 1</td><td>&lt;375 kyr</td></tr></tbody></table>

Additionally, an important and common usage of U-series age constraints is to apply age limits from short-lived nuclides to interpret longer-lived disequilibria in terms of magmatic processes. For example, when a lava of unknown age has significant $^{226}\text{Ra-}^{230}\text{Th}$ disequilibria, it is reasonable to assume that the $^{230}\text{Th-}^{238}\text{U}$ disequilibria has not been significantly changed by decay and can be interpreted in terms of melting and melt transport processes (e.g., [26,42,50,107,112,182,194]). Similarly, if a solidified lava has significant $^{210}\text{Pb-}^{226}\text{Ra}$ disequilibria, then measurements of $^{226}\text{Ra}$-
$^{230}$Th disequilibria in that lava can also be reasonably interpreted in terms of magmatic processes [42,52,80,100,186,195–197].

## 4.2. Internal Mineral Isochron Method

Internal isochrons make use of the differential retention of certain radionuclides by different minerals (e.g., Ra is preferentially incorporated into feldspar). This technique requires phenocryst-rich samples and is therefore not suitable for MORB glasses. Even when samples contain abundant phenocrysts, the resolution of this method is limited by uncertainties in the phenocryst ages (i.e., have all the phenocrysts resided in the crustal storage time for a significant period relative to the half-life of the daughter nuclide?).

Nonetheless, internal isochrons (typically mineral–mineral or mineral–glass) provide a powerful method to determine the age of crystallization in volcanic rocks (Figure 25). This technique was first applied to U-series dating by Allègre [5] and has since been applied to numerous volcanic settings [13,159,194,198–204].

<figure><img src="image_30.png" /></figure>

**Figure 25.** Schematic of internal isochrons in ($^{230}$Th/$^{232}$Th) versus ($^{238}$U/$^{232}$Th) space. Isochrons shown for a zero-age lava and for lavas with ages of one and two half-lives of $^{230}$Th. Arrows show directions of Th and U enrichment and subsequent decay toward the equiline, where the activities of $^{230}$Th and $^{238}$U are in equilibrium.

The U-series isochron method is akin to the isochron method applied to long-lived radiogenic systems such as $^{87}$Rb/$^{86}$Sr-$^{87}$Sr/$^{86}$Sr (see geochronological chapters in this volume), with one crucial difference: the daughter isotope is radioactive. Nonetheless, and akin to other isotope systems such as $^{87}$Rb-$^{87}$Sr, the U-series Bateman equation of ingrowth and decay (see Equation (8)) is linear when normalized to a stable element and plotted in the appropriate coordinates.

For the $^{238}$U-$^{230}$Th system, $^{232}$Th, which has a half-life of 14.01 Gyr, is effectively stable over the half-life of $^{230}$Th (<375 Kyr) and is used as a normalizing isotope. For the $^{230}$Th-$^{226}$Ra system, there is no stable, normalizing Ra isotope for $^{226}$Ra. Albeit imperfect, stable Ba is typically used as a Ra proxy. $^{235}$U-$^{231}$Pa also lacks a stable, normalizing isotope for $^{231}$Pa, and an internal mineral isochron method to the $^{235}$U-$^{231}$Pa system has yet to be used. Lead concentrations can be used as a normalizing factor for a $^{210}$Pb-$^{226}$Ra geochronometer [205]. However, internal isochrons using this technique are commonly non-linear because the timescales of magma differentiation are long with respect to the half-life of $^{210}$Pb, and thus, this technique has not been widely applied to dating lavas. In the following subsections,
we examine the application of the isochron method to determine the ages of crystallization for both $(^{238}\text{U})/(^{232}\text{Th})-(^{230}\text{Th})/(^{232}\text{Th})$ and $(^{230}\text{Th})/\text{Ba}-(^{226}\text{Ra})/\text{Ba}$.

### 4.2.1. $(^{238}\text{U})/(^{232}\text{Th})-(^{230}\text{Th})/(^{232}\text{Th})$

As discussed above, the Bateman equation for the $^{238}\text{U}$-$^{230}\text{Th}$ system, when normalized to $^{232}\text{Th}$, is linear in the space of $(^{230}\text{Th})/(^{232}\text{Th})$ versus $(^{238}\text{U})/(^{232}\text{Th})$ (see Equation (20); Figure 25). In these coordinates, b (the intercept with the equiline) is given by the expression $(^{230}\text{Th})_{\text{XS}} \cdot e^{-\lambda^{230}\text{t}}$ and provides an important measure of the long-term $(^{238}\text{U})/(^{232}\text{Th})$ ratio in the lava's source and, by extrapolation, the extent of initial $(^{230}\text{Th})/(^{238}\text{U})$ disequilibrium, and m (the slope of a linear isochron) is given by the relationship $(1 - e^{-\lambda^{230}\text{t}})$ and provides the measure of time since crystallization, as defined by:

$$t = \frac{\ln(1 - m)}{-\lambda} \qquad (35)$$

The assumptions for the $(^{238}\text{U})/(^{232}\text{Th})-(^{230}\text{Th})/(^{232}\text{Th})$ isochron method are the same as for most other isotope systems, namely that the minerals are cogenetic and coeval, have uniform $(^{238}\text{U})/(^{230}\text{Th})$ at a given time (i.e., the array is linear), and have remained closed to migration of U and Th since crystallization. However, because the timescales over which magma chamber processes operate are potentially similar to the half-life of $^{230}\text{Th}$, there are two additional concerns that must be considered, namely, (1) did the minerals grow quickly or slowly? and (2) was the magma residence time long or short relative to the half-life of $^{230}\text{Th}$ or to diffusional equilibration of U and Th between magmas and crystals?

Depending on the sequence of petrological events (e.g., melting, mixing, crystallization, assimilation, etc.) relative to the half-life of the daughter isotope being considered, mineral isochrons can be simple and easily interpreted in terms of eruption ages (Figure 25). However, nuclide behavior also can be complex and non-systematic. In other words, was crystallization a short, one-time event that occurred soon before or during the eruption, or did the crystals grow episodically over timescales on the order of tens of thousands of years? Furthermore, there is the possibility of intermittent recharge to a magma chamber and, thus, the mixing of crystals of different age and magma residence time (Figure 26). As such, because the isochron is constructed from mineral phenocrysts, which may have resided in a crustal magma chamber for a significant period, particularly for more evolved magmas, the calculated age can be a combination of both eruption age and magma chamber residence time. The best way to resolve this issue is to separate and measure just the groundmass phases, as these minerals grew syn-eruptively and thus provide an explicit eruption age. Although the effective separation of such groundmass phases is tedious and time-consuming [204], results for such careful sample preparation can provide both an eruption age and constrain the timescale of crystallization for the older phenocryst phases (Figure 27). In such cases, combining the U-series groundmass isochron age with cosmogenic $^{3}\text{He}$ concentrations (Figure 28) provide a unique and self-consistent quantification of both the eruption age and the erosion rate [204].

There are numerous examples of $^{238}\text{U}$-$^{230}\text{Th}$ mineral isochrons in arc and continental systems (e.g., [200]), but only a few in oceanic settings [108]. This dearth of $^{238}\text{U}$-$^{230}\text{Th}$ isochrons for MORB is because most MORB samples are phenocryst-poor, with the exception being some glasses that are plagioclase-phyric.

### 4.2.2. $(^{230}\text{Th})/\text{Ba}-(^{226}\text{Ra})/\text{Ba}$ Isochrons

Radium is an alkali-earth metal and, therefore, phenocrysts such as plagioclase, K-feldspar, and clinopyroxene are ideal target minerals for examining crystallization and eruption ages in young lavas using the $^{230}\text{Th}$-$^{226}\text{Ra}$ disequilibria chronometer. Use of
mineral-glass $^{226}$Ra-$^{230}$Th disequilibria requires knowledge of the initial $^{226}$Ra present in the minerals at the time of crystallization. Because no longer-lived or stable isotopes of Ra exists, most studies have used Ba concentrations (another alkali-earth metal) as an analog for the Ra initially present in the crystals (e.g., [206–208]).

<figure><img src="image_31.png" /></figure>

**Figure 26.** Three simplified, potential ($^{238}$U/$^{232}$Th)-($^{230}$Th/$^{232}$Th) isochron scenarios. (A) Polymer-miner-alic isochron representing the lava eruption age. In this scenario, a melt is extracted from the mantle, quickly crystallizes several phases relative to the half-life of $^{230}$Th, and then erupts soon thereafter as a porphyritic lava. The $^{230}$Th in the different solidified mineral phases of the lava (represented by the different colored shapes) then decays toward equilibrium for one half-life. (B) Poly-miner-alic isochron representing mixing of two populations of minerals of different ages, producing an isochron of intermediate age, which is the average of the two. This situation is often the case where magmas of like composition (e.g., basalts in a mid-ocean ridge or ocean island setting) are mixing or picking up crystals of similar composition from older cumulate piles. (C) Extraction of melt from the mantle, which rapidly crystallizes several phases relative to the half-life of $^{230}$Th and then resides in the magma chamber for about ten thousand years, followed by a second episode of crystallization, and then erupts. In this scenario, a good example would be when the second batch of minerals is the fine-grained eruptive groundmass.

The use of Ba as a stable analog for Ra requires that Ra and Ba behave identically during crystal growth, which typically is not the case. Theoretical work predicts a significant difference between $D_{Ra}$ and $D_{Ba}$, with calculated crystal/melt partition coefficients for Ra
significantly lower than those of Ba in plagioclase and pyroxene [54,209]. Not accounting for this difference will produce significant inaccuracies in calculated ages. Furthermore, Ra is preferentially excluded from these minerals, leading to lower ($^{226}$Ra)/Ba ratios in crystals than in the magmas from which they crystallized. This difference in partitioning behavior will produce the observed non-collinearity between multiple cogenetic phases, even in the absence of open-system behavior. As a result, with the Ra-Ba chronometer, it is difficult to distinguish the effects of initial Ra-Ba fractionation from decreases in ($^{226}$Ra)/Ba due to the aging of crystals.

<figure><img src="image_32.png" /></figure>

**Figure 27.** U-Th mineral isochrons for Bluewater flow calculated using (1) all groundmass and phenocryst phases and (2) groundmass phases only. These isochrons are determined using a non-linear, maximum-likelihood method to estimate the slope and intercept of the best-fit line and Monte Carlo methods to estimate 95% confidence limits [210]. The error bars on the data points represent two-sigma analytical uncertainties. Larger uncertainties for the olivine and plagioclase phenocrysts are due to their much lower concentrations. Note the significant difference in ($^{238}$U/$^{232}$Th) between the plagioclase phenocrysts and plagioclase groundmass suggesting that either their magma sources were distinct or inclusions (mineral or melt) are causing the difference. The U and Th concentrations of the olivine phenocrysts are also higher than expected from experimental partitioning [211], suggesting again that mineral or melt inclusions are controlling the U and Th budget. These observations demonstrate the importance of using the lava's groundmass phases for determining its eruption age.

Cooper et al. [159] developed a technique to correct $^{226}$Ra-$^{230}$Th crystal ages for the large difference in $D_{Ra}$ and $D_{Ba}$ in plagioclase and then applied this method to minerals separated from lavas of the 1955 eruption of Kilauea Volcano, Hawaii. Because the compositions of plagioclase and pyroxene were in equilibrium with the liquid in which they had erupted, they argued that their calculated crystallization age of 1000 ± 350 years. (or given the lava's age, 950 ± 350 years before the eruption) provides a minimum measurement of the residence time of the magma in the crustal reservoir. From this crystallization age, they then calculated a maximum average magmatic cooling rate of 0.1 °C/yr by combining the average crystal age with the temperature interval between saturation and eruption. Subsequent studies have applied this method, or some other variant of the application of the Ra
chronometer, to MORB and OIB from a variety of locations, including the 1996 North Gorda Ridge eruption [108] and Iceland [203,212]. These studies demonstrated that $^{226}$Ra-$^{230}$Th dating of crystals can be an effective tool for examining magmatic processes that occurred within the past ~10 kyr, although a correction for impurities in the bulk separates must be applied and the effects of differential partitioning behavior of Ra and Ba considered and accounted for [88,213].

<figure><img src="image_33.png" /></figure>

**Figure 28.** $^{3}$He exposure age as a function of erosion rate for the three sites sampled for exposure-age dating on the Bluewater flow. Two-sigma analytical uncertainties are shown by the dashed lines. Errors in production rate, scaling, and attenuation length will shift the curves uniformly up or down, but not change the relative distance between them. Given that the exposure age must equal the eruption age at all sites, different curves indicate that different erosion rates must have affected each site. The U-Th isochron age with one sigma (68% confidence level) and two-sigma (95% confidence level) are shown. The extrapolation of the U-Th ages to the $^{3}$He age curves provides an estimate of the erosion rates for the different samples (indicated by arrows at left of figure). Minimum erosion rate for BW-1 is calculated using the geomorphologic estimate of Dunbar and Phillips [214] who infer at least 10 cm of erosion. The maximum erosion rate of ZB-93–1 is determined by extrapolation of the upper age limit (two-sigma) for the U-Th groundmass-mineral isochron Modified from Sims et al [204].

Additionally, it is also important to note that because the timescales of crystallization are comparable to the half-life of $^{226}$Ra, the calculated ages represent an integrated growth rate of a crystal or crystal population. Continuous crystal growth over the timescales of the half-life $^{226}$Ra is an important issue that was not considered in these past studies, but which can significantly impact the calculated crystallization age. An illustrative example of how measurements of only $^{226}$Ra/$^{230}$Th in mineral phases and host glasses can lead to erroneous results was demonstrated by Sims et al. [88] on recently erupted lava bombs from Mount Erebus volcano (Figure 29). The calculated $^{226}$Ra-$^{230}$Th isochron ages for anorthoclase and host phonolite glass separated from the lava bombs gave incorrect crystallization ages ranging from many hundred years to a few thousand years, depending on $D_{Ra}/D_{Ba}$ [88,89]. While these ages may seem reasonable, measurements of ($^{228}$Ra/$^{232}$Th) in these same crystals were out of equilibrium, which is highly problematic for the $^{226}$Ra-$^{230}$Th isochron ages given that $^{228}$Ra has a half-life of only 5.77 years. To resolve this discrepancy,
Sims et al. [88] measured all the relevant nuclides from the ²³⁸U, ²³⁵U, and ²³²Th decay series in the anorthoclase megacrysts and glasses of recently erupted, known-age lava bombs and then created a finite-addition numerical model which incorporates the effects of magma recharge, anorthoclase crystal growth, and radioactive decay in both the magma and the anorthoclase crystals. Their model results show that by incorporating nuclide in-growth and decay during continuous crystallization they can successfully reproduce the measured ²²⁶Ra/²³⁰Th, ²¹⁰Pb/²²⁶Ra, ²²⁸Ra/²³²Th, and ²²⁸Th/²³²Th in both the anorthoclase megacrysts and phonolite melt (Figure 29E), demonstrating that when the timescale of crystallization is comparable to the half-life of ²²⁶Ra, the simple ²³⁰Th-²²⁶Ra isochron techniques typically used in most U-series studies provide incorrect ages. Thus, caution is warranted when interpreting ²²⁶Ra-²³⁰Th ages in the absence of other age constraints, either within the decay series or with other dating methods.

<figure><img src="image_34.png" /></figure>

**Figure 29.** Crystal ages for Mount Erebus anorthoclase megacrysts [88]. (A) Backscattered electron image (BSE) and Ca element map of a single anorthoclase crystal. The darkest areas in the BSE image are the anorthoclase host crystals. The abundant irregular areas, slightly brighter than the anorthoclase, are melt inclusions. Other included phases are pyroxene (light grey phase on the Ca map), apatite (second brightest in BSE, highlighted in red on the Ca map), and magnetite (brightest in BSE, dark on Ca map). The size of the apatite crystals has been slightly increased to improve visibility. (B) Ra/Ba isochrons. ²²⁶Ra/Ba E-3 versus ²³⁰Th/Ba E-3 for anorthoclase and glass separates plotted on a conventional Ba-normalized ²³⁰Th-²²⁶Ra isochron diagram. (C) Calculated Ra age versus D<sub>Ra</sub>/D<sub>Ba</sub> showing how the age of crystallization changes as a function of D<sub>Ra</sub>/D<sub>Ba</sub>. Calculated D<sub>Ra</sub>/D<sub>Ba</sub> values for Or<sub>18</sub> and An<sub>17</sub> are taken from Blundy & Wood [54]. The value of D<sub>Ra</sub>/D<sub>Ba</sub> of 0.1 is a lower limit (because of potential ²²⁶Ra decay). (D) Time-series data for ²²⁸Ra over the 34-year historical record from 1972 to 2005 in years since 2012 from Sims et al. [88] and Reagan et al. [215]. (²²⁸Ra/²³²Th)* was determined using ²²⁸Th/²³²Th measured by alpha spectroscopy as a proxy. The half-life of ²²⁸Th is 1.9 years, so for this 2005 sample (measured in 2008) the (²²⁸Ra/²³²Th)* inferred from (²²⁸Th/²³²Th) is a minimum value as the system had not yet attained full equilibrium. (E) For the Erebus lavas the equilibrium (²¹⁰Pb/²²⁶Ra) and (²²⁷Ac/²³¹Pa) suggest that the magma residence time is long compared with the half-lives of ²¹⁰Pb (t<sub>1/2</sub> = 22.6 years) and ²²⁷Ac (t<sub>1/2</sub> = 21.77 years). However, the observation that (²²⁸Ra/²³⁰Th) (t<sub>1/2</sub>=5.77 years) is out of equilibrium in the anorthoclase suggests that the Erebus magmatic differentiation and crystallization processes are continuous and continuing. With this observation in mind, an open-system, finite-element,
continuous-crystallization model was developed that incorporates ingrowth and decay of the different nuclides in the continuously growing anorthoclase crystals and associated phonolitic melt. Because the Ba and Th concentrations are measured in both the anorthoclase crystals and the melt, they can constrain the anorthoclase/phonolite melt partitioning of Ba and Th. It is found that the only way to replicate the observed Ba and Th concentrations and $^{238}$U, $^{235}$U, and $^{232}$Th decay series data is to incorporate magma recharge into the shallow system. To successfully model the present data set the recharge rate has to exceed the crystallization rate, which implies that the shallow magma reservoir within Mount Erebus is growing. This has important implications for hazard assessment as the modeled duration of reservoir growth is 2000 years, which coincides, within error, with the ages of the two youngest lavas from the flanks of Mount Erebus (the Northwest and Upper Ice Tower Ridge flows; [216]). Modified after Sims et al. [88].

### 4.3. U-Th-Ra 'Model' Ages

If the initial extent of disequilibrium at the time of eruption is known, then the difference between the initial activity ratios and the measured activity ratios can be used to determine the lava's eruption age. These "model ages" improve the resolution of the U-series ages by roughly an order of magnitude for each parent-daughter pair; by tens of thousands of years for $^{230}$Th model ages, and by hundreds to thousands of years for $^{226}$Ra model ages.

It is important to note that such U-series ages are 'model ages', in that they assume an unvarying source (both spatially and temporally), and represent the time since partial melting and are, therefore, the integration of a lava's magma chamber residence time and its eruption age. When determining U-series model ages, the following conditions must be met and demonstrated:

1. U-series disequilibria were initially created by primary magmatic processes and the lavas being dated have not undergone significant secondary alteration. The lavas have remained a closed system with respect to Th/U, Ra/Th, and Pa/U after eruption.

2. After melt generation and transport, the magma has not resided in a magma chamber for periods of time that are significant relative to the half-lives of $^{230}$Th, $^{231}$Pa, and $^{226}$Ra. If it has, then the determined age will be the sum of the integrated crustal residence time and the age of the lava, as the disequilibria is assumed to have been initiated by a melting process, which is when the radiometric clock was started. Furthermore, if a magma chamber is filled repeatedly and slowly with respect to the half-life of the geochronometer nuclide, then the age of the erupted lava will be the sum of the integrated magma residence time and the time since eruption,

3. The initial extent of disequilibria in all lavas under consideration is the same, which requires the sources of all lavas to have constant Th/U, Ra/Th, and Pa/U. This further requires that lavas were generated by similar degrees of melting at similar depths of melting from a relatively homogenous mantle source that had not recently lost partial melts. It also requires that all melts have similar transport rates (e.g., [26,52,53,107]).

Thus, if after eruption the lava remains a closed system, disequilibria between parent and daughter will decay back to equilibrium at a rate proportional to the half-life of the daughter isotope, and this decay can be used for dating lavas of unknown age (Figures 4 and 24), given that the above conditions are met.

If the above requirements are met, then U-Th and Th-Ra model ages of unknown age lavas can be determined from the relationships:

$$t = -\frac{1}{\lambda_{230}} \ln \left[ \frac{\left( \frac{^{230}\text{Th}}{^{232}\text{Th}} \right)_{\text{m}} - \left( \frac{^{238}\text{U}}{^{232}\text{Th}} \right)_{\text{m}}}{\left( \frac{^{230}\text{Th}}{^{232}\text{Th}} \right)_{0} - \left( \frac{^{238}\text{U}}{^{232}\text{Th}} \right)_{\text{m}}} \right] \quad (36)$$
and

$$t = -\frac{1}{\lambda_{226}} \ln \left[ \frac{\left( \frac{^{226}\text{Ra}}{^{230}\text{Th}} \right)_{\text{m}} - 1}{\left( \frac{^{226}\text{Th}}{^{233}\text{Th}} \right)_{0} - 1} \right] \quad (37)$$

where the subscript m represents the sample's measured activity ratio, the subscript 0 is its initial activity ratio upon eruption, and $\lambda_{230}$ and $\lambda_{226}$ are the decay constants for $^{230}\text{Th}$ and $^{226}\text{Ra}$, respectively. In Equation (36), $^{232}\text{Th}$ is used as a stable analog for $^{230}\text{Th}$ [5,13]. The different form of Equation (37) is used because no stable analog exists for $^{226}\text{Ra}$. The equation for determining $^{231}\text{Pa}$ model ages is of a similar form to the equation for $^{226}\text{Ra}$.

In these equations, the one additional unknown, besides time, is the initial activity ratio [i.e., $(^{230}\text{Th}/^{232}\text{Th})_0$, $(^{226}\text{Ra}/^{230}\text{Th})_0$, and $(^{231}\text{Pa}/^{235}\text{U})_0$. U-series ages are then determined by estimating this initial activity ratio from data for young, “known-age” lavas and comparing this with the sample’s measured activity ratio.

Probably the best examples of model age dating are for MORB [26,52,53,107,182–186]. The most well-studied ridge in the mid-ocean ridge system with regard to U-series dating is the 9–10° N segment of the East Pacific Rise [26,52,53,107,184–186] because this particular segment of the EPR was a focus of the US National Science Foundation Ridge and Ridge Y2K programs. In studies of samples from this region, $^{238}\text{U}-^{230}\text{Th}-^{226}\text{Ra}$ dating evolved to include both larger sample suites of well-located and geologically well-constrained samples, as well as measurements of $^{226}\text{Ra}$, which provided finer time resolution. These spatial and temporal geological constraints for the axial lavas were also essential to the application of $(^{226}\text{Ra}/^{230}\text{Th})$ for dating off-axis samples. All lavas collected within the axial valley trough at 9-10°N EPR were geologically constrained to be less than 375 years old, such that $(^{230}\text{Th}/^{232}\text{Th})_0$, $(^{226}\text{Ra}/^{230}\text{Th})_0$, and $(^{231}\text{Pa}/^{235}\text{U})_0$ were all known explicitly for the magmatic system [26].

A fundamental observation of Sims et al. [107] was that many off-axis samples have higher U-Th and Th-Ra disequilibria, and thus younger model ages, than would be predicted for their off-axis locations and the time-integrated spreading rate of the ridge (Figure 30). Collectively, the resulting understanding from the U-series chronology and enhanced observations and mapping has led to a fundamental paradigm shift for models of mid-ocean ridge construction. The distribution of young lavas requires that a significant component of oceanic crustal construction at the fast-spreading EPR involves multiple volcanic and tectonic processes acting in concert to form a complex patchwork of lava ages [26,31,52,53,107,192,217], whereas previously it was assumed that nearly all volcanic activity was limited to a very narrow axial summit trough. One possible mechanism to explain such extensive off-axis volcanism is that bending and unbending stresses play a key role in triggering volcanism on the flanks of the EPR, by simultaneously opening tensile cracks near the surface and increasing the pore pressure of melt bodies trapped in the lower crust [217]. Similar systematics were also suggested for samples from the Southwest Indian Ocean Ridge (SWIR; Figure 31) where young ages were inferred for lavas far off-axis and high on fault scarps [110]. These important observations provide evidence that crustal accretion of young volcanic rocks is not confined to a narrow central spreading axis as had long been the standard paradigm.
<figure><img src="image_35.png" /></figure>

**Figure 30.** (A) ($^{238}$U/$^{232}$Th) versus ($^{230}$Th/$^{232}$Th) for the 9°50′ N off-axis samples compared to the 9–10° N EPR axial samples [26,107]. Errors represent two sigma analytical uncertainties. Solid line represents the zero-age axial trend used for calculating “trend line” ages; dashed line is the equiline. For most of the off-axis samples ($^{230}$Th/$^{232}$Th) and ($^{238}$U/$^{232}$Th) are similar to the axial samples from 9°50′ N EPR; the three samples that lie below the range of axial values are interpreted to be older. The inset provides regional context and shows the 9–10° N EPR axial samples [26,107] compared with the Siqueiros Transform samples [26,41]. Also shown are the near-axial basalts from 9°37′ N EPR [183]. (B) $^{230}$Th/$^{238}$U versus $^{226}$Ra/$^{230}$Th for 9°50′ N EPR off-axis samples compared with the 9°50′ N EPR axial samples [26,107]. Error bars represent two sigma analytical uncertainties. Because of the significant difference in the half-lives of $^{226}$Ra and $^{230}$Th, ($^{226}$Ra/$^{230}$Th) decays back to equilibrium before ($^{230}$Th/$^{238}$U) (note decay curve on plot). The correlation observed among the zero-age samples (shown by solid line) is used as the trend line from which Th-Ra ages are calculated (see text for details). Isochrons of equal age are plotted. Note that several off-axis samples have ($^{226}$Ra/$^{230}$Th) in equilibrium, indicating that their ages are older (>8 ka) than can be determined through $^{230}$Th-$^{226}$Ra disequilibria. Inset provides regional context and shows ($^{230}$Th/$^{238}$U) versus ($^{226}$Ra/$^{230}$Th) for the 9°50′ N EPR off-axis samples from this study compared with the 9–10° N EPR axial samples [26,107], three “near-axis” 9°30′ N EPR samples, including one E-type from close to the 9°37′ N DEVAL [188], and the Siqueiros transform samples [41]. (C) 9°50′ N EPR samples ($^{230}$Th/$^{238}$U) versus their distance from the center of the ridge axis, or Axial Summit Trough. Decay curves show how ($^{230}$Th/$^{238}$U) decreases as a function of distance from the AST. Initial values for decay curves come from the average, red, and high and low values of ($^{230}$Th/$^{238}$U) for the 9°50′ N EPR axial samples as reported in Sims et al. [26]. (D) 9°50′ N EPR samples ($^{226}$Ra/$^{230}$Th) versus their location relative to the AST. Decay curves show how ($^{226}$Ra/$^{230}$Th) decreases as a function of distance from the AST. The initial values for these decay curves come from the average, red, and high and low values of ($^{226}$Ra/$^{230}$Th) for the 9°50′ N EPR axial samples [26]. (E) U-Th and Th-Ra trend-line model ages for the 9°50′ N off-axis samples (Figure 9; [107]) compared with ages calculated from the half-spreading rate (5.5 cm/yr) as determined from paleomagnetic data [218]. Note that a majority of the measured off-axis basalt samples from the
9°50′ N section of the EPR have U-Th and Th-Ra disequilibria that are larger, and model ages that are younger, than would be predicted from their off-axis distance and the time-integrated spreading rate. There are, however, a few off-axis samples with U-Th model ages that are consistent with their spreading rate ages.

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

**Figure 31.** Measured U-series disequilibria in SWIR MORB as a function of distance from the spreading axis. The figure on the left is a bathymetric image with colors corresponding to elevation in meters below sea level, along with sample locations. (a) Variable ($^{230}$Th/$^{238}$U) for both spreading axis lavas (green triangles) and rift valley lavas (red circles) with respect to the spreading axis indicate no correlation between disequilibrium and crustal age. Decay curves bracketing the spreading axis track the amount of disequilibrium in a basalt assumed to have an initial excess of ~27% as a function of distance from the axis, given a spreading rate of 7.1 mm/yr. Based on their measured U-Th disequilibrium and current locations, five rift valley lavas could not have erupted within the spreading axis. For comparison, dashed decay curve represents spreading rate of 55 mm/yr (9° N EPR). (b) As in (a), solid and dashed curves contrast different spreading rate decay curves for ($^{226}$Ra/$^{230}$Th) disequilibrium. Three of the five lavas with U-Th disequilibrium have excess $^{226}$Ra constraining their ages to <8 ka. Alternatively, at 7.1 mm/yr, an age < 8 kyr is equivalent to 55 m of movement since emplacement.

4.4. * Dating Young Volcanic Rocks with* $^{226}$Ra-$^{210}$Pb

While there have been numerous measurements of (²¹⁰Pb/²²⁶Ra) in young vol-
canic arc lavas [187,205,219–222], there are relatively few such measurements in MORB
(e.g., [50,53,80,100,105,187]) and in oceanic intra plate settings (e.g., [42,50,53,223,224]).

For most MORB lavas, $^{226}$Ra-$^{210}$Pb disequilibria have been used as a chronometer to establish whether lavas are greater or less than ~100 years old based on the half-life of $^{210}$Pb ($t_{1/2}$ = 22 years) and the presence or absence of disequilibria [42,50,53,187,225]. However, as recently shown, for some lavas it is also possible to determine model eruption ages on a decadal scale using $^{226}$Ra-$^{210}$Pb disequilibria [80].
If the required source homogeneity criteria are met, then $^{226}$Ra-$^{210}$Pb model ages can be determined from the relationship:

$$T = \frac{1}{\lambda_{210}} \ln \frac{(^{210}\text{Pb}/^{226}\text{Ra})_m - 1}{(^{210}\text{Pb}/^{226}\text{Ra})_o - 1} \quad (38)$$

where $(^{210}\text{Pb}/^{226}\text{Ra})_m$ is the measured ratio, $(^{210}\text{Pb}/^{226}\text{Ra})_o$ is the calculated initial ratio at the time of eruption of the lavas, and $\lambda_{210}$ is the decay constant of $^{210}\text{Pb}$. Again, calculating $(^{210}\text{Pb}/^{226}\text{Ra})$ model eruption ages requires an assessment of the initial extent of $^{210}\text{Pb}$ excess at the time of eruption of the samples in question.

Sims et al. [80] used $(^{210}\text{Pb}/^{226}\text{Ra})$ to calculate model submarine lava ages from the ABE vent site in the Lau basin with decadal scale resolution. A crucial criterion for using $(^{210}\text{Pb}/^{226}\text{Ra})$ to calculate model ages at this location was that the ABE vent site lavas had essentially constant $(^{226}\text{Ra}/^{230}\text{Th})$ and $(^{230}\text{Th}/^{238}\text{U})$ disequilibria (Figure 32). Given this constraint, Sims et al. [80] calculated the time of $(^{210}\text{Pb}/^{226}\text{Ra})$ decay since eruption, obtaining model eruption ages of 0 to 100 years, with an average eruption age for the vent field of roughly 53 ± 29 years at the time of analysis in 2013.

<figure><img src="image_38.png" /><figcaption>Figure 32. (A) Map of the ABE vent field showing sample locations and model eruption ages. Image is DLS-120 Side Scan backscatter image. Vent sites, faults, and approximate location of faults and lava flow fronts from Ferrini et al. [226]. (B) $(^{210}\text{Pb}/^{226}\text{Ra})$ versus $(^{226}\text{Ra}/^{230}\text{Th})$ showing ABE vent field lavas and age trajectories calculated assuming $(^{210}\text{Ra}/^{226}\text{Ra})_0 = 1.82$ or 2.0 and $(^{226}\text{Ra}/^{230}\text{Th})_0 = 3.02$ or 3.55 and then interpolating a model zero-age isochron between the two end-members. The majority of samples have $(^{210}\text{Pb}/^{226}\text{Ra}) > 1$, indicating that they erupted within the past 100 years. RC-086 likely erupted more than 100 years ago; if this sample had $(^{226}\text{Ra}/^{230}\text{Th}) = 3.05-3.63$, then 1800–2500 years have elapsed since eruption. (C) Sample photos taken by ROV Jason 2 arranged in order of increasing model age from left to right. Modified from Sims et al. [80].</figcaption></figure>

## 4.5. Dating Young Volcanic Rocks with $^{210}\text{Po}$-$^{210}\text{Pb}$

In subaerial lavas and shallow submarine volcanoes, $^{210}\text{Po}$ ($t_{1/2}$ = 138.4 days) is posited to degas entirely at magmatic temperatures, whereas its grandparent nuclide, $^{210}\text{Pb}$
(t₁/₂ = 22.6 years), remains almost entirely in the melt (e.g., [227] and references therein). The intermediate radionuclide in the decay chain, ²¹⁰Bi, has a relatively short half-life (t₁/₂ = 5.012 days) and is thus assumed to be in equilibrium with its parent, ²¹⁰Pb. Therefore, by assuming that ²¹⁰Po is completely degassed and thus has an initial concentration of zero at the time of eruption, a series of measurements of the ²¹⁰Po activity can be made in the weeks and months following eruption and an ingrowth curve can be fitted to the data to determine the maximum eruption age (Figure 33).

<figure><img src="image_39.png" /><figcaption>**Figure 33.** (A) Extrapolated best-fit ²¹⁰Po ingrowth curve assuming 100% Po degassing. Calculated maximum eruption age shown with 2-sigma analytical error. (B,C) Bathymetric images showing dome growth occurring sometime between March 1999 and April 2005. Modified after Sims et al. [42].</figcaption></figure>

The ²¹⁰Po activity as a function of time is described by the Bateman equation:

$$ (^{210}\text{Po}) = (^{210}\text{Po})_0 e^{-\lambda t} + (^{210}\text{Pb})(1 - e^{-\lambda t}) \quad (39) $$

Assuming $(^{210}\text{Po})_0 = 0$, this simplifies to:

$$ (^{210}\text{Po}) = (^{210}\text{Pb})(1 - e^{-\lambda t}) \quad (40) $$

The assumption that initial ²¹⁰Po activity is zero generally holds true for subaerial basaltic eruptions and for submarine eruptions at less than 3 km depth [228]. Lavas erupted at greater depths are subject to higher hydrostatic pressures, which may limit the degassing efficiency. Rapid eruption and quenching of lavas may also reduce the degassing efficiency in differentiated magmas [202,221]. Nevertheless, lavas that actively degas until eruption most commonly lose >90% of their Po [195,220,229–231], which means that calculated eruption ages assuming ²¹⁰Po = 0 may typically be too old by up to about 3 weeks.

Validation of the ²¹⁰Po dating technique, and the assumption that Po is completely degassed during eruption, has been completed in a few important studies, indicating that ²¹⁰Po can be a powerful tool to precisely date recent, unobserved volcanic eruptions. Two seminal examples include: