On the Application of ²³⁸U, ²³⁵U, and ²³²Th Decay Series Disequilibria to Understanding the Processes and Timescales of Oceanic Basalt Petrogenesis

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 ²³⁸U, ²³⁵U, and ²³²Th decay series 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.

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 has 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 models 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.

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 x 10⁶ years) relative to the age of the elements in the solar system and the age of the 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 [2] 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., [3]).

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 (Figure 3 and Figure 4). Consider, for example, the ²³⁸U series:${}^{92}{}_{238}\mathrm{U}⇥\mapsto {}^{92}{}_{234}\mathrm{U}⟶{}^{90}{}_{230}\mathrm{T}\mathrm{h}⟶{}^{88}{}_{226}\mathrm{R}\mathrm{a}⟶\dots ⟶{}^{82}{}_{206}\mathrm{P}\mathrm{b}.$ As will be shown in the following derivation, radioactive secular equilibrium within a decay series is defined as:${\mathsf{\lambda }}_1{\mathrm{N}}_1={\mathsf{\lambda }}_2{\mathrm{N}}_2={\mathsf{\lambda }}_3{\mathrm{N}}_3$ or (A₁) = (A₂) = (A₃), where ${\mathsf{\lambda }}_i$ is the decay constant of isotope $i$, ${\mathrm{N}}_i$ is the number of atoms of isotope $i$, and (${\mathrm{A}}_i$) is the activity of isotope $i$ (parenthesis around an isotope denote activity), and is defined as:

$$(A)=-{\displaystyle \frac{\mathrm{d}\mathrm{N}}{\mathrm{d}\mathrm{t}}}=\mathsf{\lambda }\mathrm{N}$$

(1)

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 ${\mathrm{N}}_1⟶{\mathrm{N}}_2⟶{\mathrm{N}}_3:$

$$-{\displaystyle \frac{\mathrm{d}{\mathrm{N}}_1}{\mathrm{d}\mathrm{t}}}={\mathsf{\lambda }}_1{\mathrm{N}}_1$$

(2)

$${\displaystyle \frac{\mathrm{d}{\mathrm{N}}_2}{\mathrm{d}\mathrm{t}}}={\mathsf{\lambda }}_1{\mathrm{N}}_1-{\mathsf{\lambda }}_2{\mathrm{N}}_2$$

(3)

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 β^-). Table includes common parent-daughter pairs utilized in studying igneous processes and their associated measurement techniques. Modified from [4].

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 β^-). Table includes common parent-daughter pairs utilized in studying igneous processes and their associated measurement techniques. Modified from [4].

Since N₁ undergoes radioactive decay, its abundance at time t is given by:

$${\mathrm{N}}_1={\mathrm{N}}_1^0{\mathrm{e}}^{-{\mathsf{\lambda }}_1\mathrm{t}}$$

(4)

where the term ${\mathrm{N}}_1^0$is the number of atoms at t = 0. Substituting Equation 4 into Equation 3 yields:

$${\displaystyle \frac{\mathrm{d}{\mathrm{N}}_2}{\mathrm{d}\mathrm{t}}}={\mathsf{\lambda }}_1{\mathrm{N}}_1^0{\mathrm{e}}^{-{\mathsf{\lambda }}_1\mathrm{t}}-{\mathsf{\lambda }}_2{\mathrm{N}}_2$$

(5)

and

$${\displaystyle \frac{\mathrm{d}{\mathrm{N}}_2}{\mathrm{d}\mathrm{t}}}+{\mathsf{\lambda }}_2{\mathrm{N}}_2-{\mathsf{\lambda }}_1{\mathrm{N}}_1^0{\mathrm{e}}^{-{\mathsf{\lambda }}_1\mathrm{t}}=Ø$$

(6)

Solving for ${\mathrm{N}}_2$:

$${\mathrm{N}}_2={\displaystyle \frac{{\mathsf{\lambda }}_1}{{\mathsf{\lambda }}_2-{\mathsf{\lambda }}_1}}{\mathrm{N}}_1^0\left({\mathrm{e}}^{-{\mathsf{\lambda }}_1\mathrm{t}}-{\mathrm{e}}^{-{\mathsf{\lambda }}_2\mathrm{t}}\right)+{\mathrm{N}}_2^0{\mathrm{e}}^{-{\mathsf{\lambda }}_2\mathrm{t}}$$

(7)

The first term gives the number of atoms ${\mathrm{N}}_2$ that formed from decay of ${\mathrm{N}}_1$, but have not yet decayed to N₃. The second term gives the number of atoms ${\mathrm{N}}_2$ that remain from initial ${\mathrm{N}}_2^0$. If there are no atoms of the original daughter ${\mathrm{N}}_2^0$, then this reduces to:

$${\mathrm{N}}_2={\displaystyle \frac{{\mathsf{\lambda }}_1}{{\mathsf{\lambda }}_2-{\mathsf{\lambda }}_1}}{\mathrm{N}}_1^0\left({\mathrm{e}}^{-{\mathsf{\lambda }}_1\mathrm{t}}-{\mathrm{e}}^{-{\mathsf{\lambda }}_2\mathrm{t}}\right)$$

(8)

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). Arrow widths correspond to the magnitude of the decay constant (λ). 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 [4].

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). Arrow widths correspond to the magnitude of the decay constant (λ). 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 [4].

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 ${\mathrm{N}}_2$:${\mathrm{e}}^{-{\mathsf{\lambda }}_2\mathrm{t}}⟶Ø,\mathrm{s}\mathrm{o}$:

$${\mathrm{N}}_2\approx {\displaystyle \frac{{\mathsf{\lambda }}_1}{{\mathsf{\lambda }}_2-{\mathsf{\lambda }}_1}}{\mathrm{N}}_1^0\left({\mathrm{e}}^{-{\mathsf{\lambda }}_1\mathrm{t}}\right)\approx {\displaystyle \frac{{\mathsf{\lambda }}_1}{{\mathsf{\lambda }}_2-{\mathsf{\lambda }}_1}}{\mathrm{N}}_1$$

(9)

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

$${\mathsf{\lambda }}_2-{\mathsf{\lambda }}_1\approx {\mathsf{\lambda }}_2,{\mathrm{N}}_2\approx {\displaystyle \frac{{\mathsf{\lambda }}_1}{{\mathsf{\lambda }}_2}}{\mathrm{N}}_1,\mathrm{and}{\mathrm{N}}_2{\mathsf{\lambda }}_2\approx {\mathrm{N}}_1{\mathsf{\lambda }}_1$$

(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_1/2 = 10 hours) has a half-life ten times longer than its shorter-lived daughter (N₂; t_1/2 = 1 hour). Two scenarios are illustrated. In the first example, the daughter isotope has been entirely removed from the parent (Figure 4A-C), so ${\mathrm{N}}_2^0$ = 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 (${\mathsf{\lambda }}_2$ – ${\mathsf{\lambda }}_1$)/${\mathsf{\lambda }}_1$, which, as shown above, is reasonably approximated as ${\mathsf{\lambda }}_2$/${\mathsf{\lambda }}_1$ since ${\mathsf{\lambda }}_2-{\mathsf{\lambda }}_1\approx {\mathsf{\lambda }}_2$ 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 (${\mathsf{\lambda }}_2$ – ${\mathsf{\lambda }}_1$)/${\mathsf{\lambda }}_1$.

Now consider the actual ²³⁸U-²³⁰Th system${}^{92}{}_{238}\mathrm{U}⇥\mapsto {}^{92}{}_{234}\mathrm{U}⟶{}^{90}{}_{230}\mathrm{T}\mathrm{h}⟶{}^{88}{}_{226}\mathrm{R}\mathrm{a}⟶\dots ⟶{}^{82}{}_{206}\mathrm{P}\mathrm{b}$. 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}{\mathrm{T}\mathrm{h})}_{\mathrm{A}}=({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}+({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{S}}$$

(11)

where (${}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}}$ is the total activity, $({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}$ is the unsupported excess activity ($\mathsf{\lambda }\mathrm{N}$, where parentheses denote activity) and $({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{S}}$ is the ²³⁰Th supported by the decay of ²³⁸U. In a system in which geochemical fractionation of Th from U ceases, the unsupported excess, (${}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}$, decays according to simple radioactive principles, such that at time t:

$$({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}=({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}^0{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}$$

(12)

where $({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}^0$ is the initial unsupported excess. Normalized to (²³²Th), the total (²³⁰Th/²³²Th) excess is:

$${\left({\displaystyle \frac{{}^{230}\mathrm{T}\mathrm{h}}{{}^{232}\mathrm{T}\mathrm{h}}}\right)}_{\mathrm{A}\mathrm{X}}={\left({\displaystyle \frac{{}^{230}\mathrm{T}\mathrm{h}}{{}^{232}\mathrm{T}\mathrm{h}}}\right)}_{\mathrm{A}\mathrm{X}}^0{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}$$

(13)

Figure 4. Attainment of radioactive secular equilibrium after an initial fractionation between hypothetical long-lived parent (N₁; t_1/2 = 10 hours) and short-lived daughter (N₂; t_1/2 = 1 hour) 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 [4].

Figure 4. Attainment of radioactive secular equilibrium after an initial fractionation between hypothetical long-lived parent (N₁; t_1/2 = 10 hours) and short-lived daughter (N₂; t_1/2 = 1 hour) 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 [4].

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

$${\mathrm{N}}_2={\displaystyle \frac{{\mathsf{\lambda }}_1}{{\mathsf{\lambda }}_2-{\mathsf{\lambda }}_1}}{\mathrm{N}}_1^0\left({\mathrm{e}}^{-{\mathsf{\lambda }}_1\mathrm{t}}-{\mathrm{e}}^{-{\mathsf{\lambda }}_2\mathrm{t}}\right)$$

In terms of ${}^{230}\mathrm{T}\mathrm{h}$:

$${{}^{230}\mathrm{T}\mathrm{h}}_{\mathrm{A}\mathrm{S}}={\displaystyle \frac{{\mathsf{\lambda }}_{234}}{{\mathsf{\lambda }}_{230}-{\mathsf{\lambda }}_{234}}}{{}^{234}\mathrm{U}}^0\left({\mathrm{e}}^{-{\mathsf{\lambda }}_{234}\mathrm{t}}-{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}\right)$$

(14)

Given that${}^{238}\mathrm{U}⟶{}^{234}\mathrm{U}⟶{}^{230}\mathrm{T}\mathrm{h}$, if ²³⁴U is in secular equilibrium with ²³⁸U, then

$\left({}^{234}\mathrm{U}\right)=\left({}^{238}\mathrm{U}\right)$ and thus ${}^{234}{\mathrm{U}}^0{\mathsf{\lambda }}_{234}={\left({}^{238}\mathrm{U}\right)}_{\mathrm{A}}$, ${\mathsf{\lambda }}_{230}-{\mathsf{\lambda }}_{234}\approx {\mathsf{\lambda }}_{230},$ and ${\mathrm{e}}^{-{\mathsf{\lambda }}_{234}\mathrm{t}}\approx 1$. From this follows that

$${\mathsf{\lambda }}_{230}{{}^{230}\mathrm{T}\mathrm{h}}_{\mathrm{A}\mathrm{S}}={\mathsf{\lambda }}_{234}{{}^{234}\mathrm{U}}^0\left(1-{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}\right)$$

(15)

$${\left({}^{230}\mathrm{T}\mathrm{h}\right)}_{\mathrm{A}\mathrm{S}}=\left({}^{238}\mathrm{U}\right)\left(1-{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}\right)$$

(16)

From Equation 12:

$$({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}=({}^{230}{\mathrm{T}\mathrm{h})}_{\mathrm{A}\mathrm{X}}^0{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}$$

and Equation 11:

$${\left({}^{230}\mathrm{T}\mathrm{h}\right)}_{\mathrm{A}}={\left({}^{230}\mathrm{T}\mathrm{h}\right)}_{\mathrm{A}\mathrm{X}}+{\left({}^{230}\mathrm{T}\mathrm{h}\right)}_{\mathrm{A}\mathrm{S}}$$

so that:

$${\left({}^{230}\mathrm{T}\mathrm{h}\right)}_{\mathrm{A}}={\left({}^{230}\mathrm{T}\mathrm{h}\right)}_{\mathrm{A}\mathrm{X}}{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}+\left({}^{238}\mathrm{U}\right)\left(1-{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}\right)$$

(17)

When normalized to ${}^{232}\mathrm{T}\mathrm{h}:$

$${\left({\displaystyle \frac{{}^{230}\mathrm{T}\mathrm{h}}{{}^{232}\mathrm{T}\mathrm{h}}}\right)}_{\mathrm{A}}={\left({\displaystyle \frac{{}^{230}\mathrm{T}\mathrm{h}}{{}^{232}\mathrm{T}\mathrm{h}}}\right)}_{\mathrm{A}\mathrm{X}}{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}+\left({\displaystyle \frac{{}^{238}\mathrm{U}}{{}^{232}\mathrm{T}\mathrm{h}}}\right)\left(1-{\mathrm{e}}^{-{\mathsf{\lambda }}_{230}\mathrm{t}}\right)$$

(18)

This is a linear equation in the space of (²³⁰Th/²³²Th) versus (²³⁸U/²³²Th) and is the underlying equation used for the ²³⁸U-²³⁰Th chronometer [5]. Similar variants exist for the other U-series chronometers. It should be noted that ²³⁴Th and ²³⁴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 (²³⁴Th t_1/2 = 24.1 d, ²³⁴Pa t_1/2 = 1.18 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.

Two end-member processes (Figure 5 and Figure 6) have been 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 [7,8,9]; and 2) and time-dependent melting, where ²³⁰Th grows in over the timescale of melt generation [10,11,12,13]. If net elemental fractionation of U/Th during partial melting is the dominant effect, then the U/Th 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.

In time-dependent adiabatic decompression melting models (Figure 5 and Figure 6), ingrowth of ²³⁰Th produces elevated (²³⁰Th/²³²Th) leading to higher (²³⁰Th/²³⁸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 [8,9]; whereas, when upwelling is slow (e.g., at mid-ocean ridges), ²³⁰Th ingrowth is the fundamental mechanism producing U-Th disequilibria.

Figure 6. Schematic illustrating the differences in the dynamic melting modeling of McKenzie [14] and the chromatographic melting model of Spiegelman and Elliott [12]. (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_p and C_d, respectively. In the subsolidus mantle, parent and daughter travel at the same velocity (the mantle upwelling rate) and maintain secular equilibrium.

Figure 6. Schematic illustrating the differences in the dynamic melting modeling of McKenzie [14] and the chromatographic melting model of Spiegelman and Elliott [12]. (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_p and C_d, respectively. In the subsolidus mantle, parent and daughter travel at the same velocity (the mantle upwelling rate) and maintain secular equilibrium.

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. [8,11,12,13,14,15,16,17,18,19,20,21,22]; however, credit for the time-dependent models (Figure 5 and Figure 6) should be given to (a) McKenzie [14] for dynamic melting and (b) Spiegelman and Elliott [12] for reactive porous flow melting, also called chromatographic melting.

While these two end-member models, dynamic (a) versus chromatographic (b) 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., [23,24]) 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 (Figure 7, Figure 8, Figure 9, Figure 10 and Figure 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.

3.2. Dynamic Melting

Solutions for ²³⁸U-²³⁰Th-²³⁶Ra disequilibria during dynamic melting were first derived by McKenzie [14] and have subsequently been presented in several papers pertaining to the production of U-series disequilibria in basaltic melts [11,18,19,21,22,47]. By assuming a constant melt rate (Γ, mass fraction per time), and a constant porosity (φ, volume fraction of melt), the (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰Th) can be solved analytically and expressed as a function of ϕ and Γ (see McKenzie [14] or Williams and Gill [11] for the full derivation):

$${\displaystyle \frac{{}^{230}\mathrm{T}\mathrm{h}}{{}^{238}\mathrm{U}}}={\displaystyle \frac{{\mathrm{F}}_{\mathrm{T}\mathrm{h}}({\mathrm{K}}_{\mathrm{U}}+{\mathsf{\lambda }}_{\mathrm{T}\mathrm{h}})}{{\mathrm{F}}_{\mathrm{U}}({\mathrm{K}}_{\mathrm{T}\mathrm{h}}+{\mathsf{\lambda }}_{\mathrm{T}\mathrm{h}})}}$$

(19)

and

$${\displaystyle \frac{{}^{226}\mathrm{R}\mathrm{a}}{{}^{230}\mathrm{T}\mathrm{h}}}=\left[{\displaystyle \frac{{\mathrm{F}}_{\mathrm{T}\mathrm{h}}\left({\mathrm{K}}_{\mathrm{U}}+{\mathsf{\lambda }}_{\mathrm{R}\mathrm{a}}\right)}{{\mathrm{F}}_{\mathrm{U}}\left({\mathrm{K}}_{\mathrm{T}\mathrm{h}}+{\mathsf{\lambda }}_{\mathrm{R}\mathrm{a}}\right)}}\right]\left[1+{\displaystyle \frac{{\mathsf{\lambda }}_{\mathrm{T}\mathrm{h}}\left({\mathrm{K}}_{\mathrm{U}}+{\mathrm{K}}_{\mathrm{U}}\right)}{\left({\mathrm{K}}_{\mathrm{T}\mathrm{h}}+{\mathsf{\lambda }}_{\mathrm{R}\mathrm{a}}\right)\left({\mathrm{K}}_{\mathrm{U}}+{\mathsf{\lambda }}_{\mathrm{T}\mathrm{h}}\right)}}\right]$$

(20)

where

$${\mathrm{F}}_{\mathrm{i}}={\displaystyle \frac{{\mathsf{\varphi }\mathsf{\rho }}_{\mathrm{f}}}{{\mathrm{D}}_{\mathrm{i}}{\mathsf{\rho }}_{\mathrm{s}}\left(1-\right)+{\mathsf{\varphi }}_{\mathrm{f}}}}$$

(21)

and

$${\mathrm{K}}_{\mathrm{i}}={\displaystyle \frac{{\mathrm{F}}_{\mathrm{T}\mathrm{h}}\left(1-{\mathrm{D}}_{\mathrm{i}}\right)}{{\mathsf{\varphi }\mathsf{\rho }}_{\mathrm{f}}}}{\mathsf{\Gamma }}_0$$

(22)

and D_i 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 8. Results of dynamic melting calculations for an adiabatically decompressing peridotitic mantle showing (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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) 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_p = 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_p = 1400 ºC. Global MORB data are shown for comparison. After Figure 7.

Figure 8. Results of dynamic melting calculations for an adiabatically decompressing peridotitic mantle showing (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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) 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_p = 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_p = 1400 ºC. Global MORB data are shown for comparison. After Figure 7.

In the dynamic melting model, as formulated by McKenzie [14], 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.

Once the critical threshold porosity is reached, the porosity of the melt zone is constant (Figure 6 and Figure 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 [11,47]. 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 similar to batch melting, in which incompatible elements reside in chemical equilibrium with the solid.

Table 1. Summary of important U-series radionuclides.

Isotope	Abundance	Half-life	Decay Constant	Reference
238U	99.2742%	4.4683 ± 0.0024 x 109 yr	1.5513 x 10–10 yr–1	[48]
234U	0.0054%	245,620 ± 260 yr	2.82206 x 10–6 yr–1	[49]
230Th	0.02%	75,584 ± 110 yr	9.1706 × 10–6 yr–1	[49]
226Ra	Trace	1,600 ± 7 yr	4.332 x 10–4 yr–1	[50]
222Rn	Trace	3.8235 ± 0.0004 d	1.8 x 10–1 d–1	[50]
210Pb	Trace	22.6 ± 0.1 yr	3.07 x 10-2 yr–1	[51]
210Po	Trace	138.4 ± 0.1 d	5.008 x 10–3 d–1	[51]
235U	0.7204%	7.0381 ± 0.0048 x 108 yr	9.8485 x 10–10 yr–1	[48]
231Pa	100%	32,760 ± 220 yr	2.1158 x 10–5 yr–1	[52]
227Ac	Trace	21.77 ± 0.02 yr	3.184 x 10–2 yr–1	[51]
232Th	99.98%	1.401 ± 0.007 x 1010 yr	4.948 x 10–11 yr–1	[53]
228Ra	Trace	5.75 ± 0.03 yr	1.21 x 10–1 yr–1	[50]
224Ra	Trace	3.6319 ± 0.0023 d	1.9085 x 10–1 d–1	[51]

Table 1. Summary of important U-series radionuclides.

Isotope	Abundance	Half-life	Decay Constant	Reference
238U	99.2742%	4.4683 ± 0.0024 x 109 yr	1.5513 x 10–10 yr–1	[48]
234U	0.0054%	245,620 ± 260 yr	2.82206 x 10–6 yr–1	[49]
230Th	0.02%	75,584 ± 110 yr	9.1706 × 10–6 yr–1	[49]
226Ra	Trace	1,600 ± 7 yr	4.332 x 10–4 yr–1	[50]
222Rn	Trace	3.8235 ± 0.0004 d	1.8 x 10–1 d–1	[50]
210Pb	Trace	22.6 ± 0.1 yr	3.07 x 10-2 yr–1	[51]
210Po	Trace	138.4 ± 0.1 d	5.008 x 10–3 d–1	[51]
235U	0.7204%	7.0381 ± 0.0048 x 108 yr	9.8485 x 10–10 yr–1	[48]
231Pa	100%	32,760 ± 220 yr	2.1158 x 10–5 yr–1	[52]
227Ac	Trace	21.77 ± 0.02 yr	3.184 x 10–2 yr–1	[51]
232Th	99.98%	1.401 ± 0.007 x 1010 yr	4.948 x 10–11 yr–1	[53]
228Ra	Trace	5.75 ± 0.03 yr	1.21 x 10–1 yr–1	[50]
224Ra	Trace	3.6319 ± 0.0023 d	1.9085 x 10–1 d–1	[51]

3.3. Chromatographic Porous Flow Melting

The steady-state chromatographic porous flow melting model of Spiegelman and Elliot [12], 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_i 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 [12]. 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 ²²⁶Ra and much smaller for ²³⁰Th and ²³¹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.

Table 2. Main parameters used in melting models.

Parameters	Notation	Units	Range
Degree of melting	F	---	0-0.2
Solid upwelling rate (1D)	W	m·yr-1	0-1
Melting rate	Γ	kg·m-3·yr-1	10−3 to 10−5
Solidus temperature	T	K	1000-1400
Length of melting column	Z	km	0-150
Plate velocity	---	m·yr-1	0.01-0.1
Matrix porosity	φ	---	0-0.05
Solid density	ρs	kg·m−3	3300
Melt density	ρf or ρm	kg·m−3	2800
Melt velocity (1D)	ω	m·yr-1
Lithosphere thickness	h	km	0-100
Decay Constants	λ	yr-1	>10-1 to 10-11
Bulk Partition Coefficient	Di	---	0.00001 to 0.1
Buoyancy flux	B	Mg·s-1	0.5-9

Table 2. Main parameters used in melting models.

Parameters	Notation	Units	Range
Degree of melting	F	---	0-0.2
Solid upwelling rate (1D)	W	m·yr-1	0-1
Melting rate	Γ	kg·m-3·yr-1	10−3 to 10−5
Solidus temperature	T	K	1000-1400
Length of melting column	Z	km	0-150
Plate velocity	---	m·yr-1	0.01-0.1
Matrix porosity	φ	---	0-0.05
Solid density	ρs	kg·m−3	3300
Melt density	ρf or ρm	kg·m−3	2800
Melt velocity (1D)	ω	m·yr-1
Lithosphere thickness	h	km	0-100
Decay Constants	λ	yr-1	>10-1 to 10-11
Bulk Partition Coefficient	Di	---	0.00001 to 0.1
Buoyancy flux	B	Mg·s-1	0.5-9

The one-dimensional (1D) chromatographic porous flow model of Spiegelman and Elliott [12] has seven variable parameters (which are not all independent): the solid 1D velocity (W), the liquid 1D velocity ($\omega$), the porosity (ϕ), the height of the melting column (∆z), the melt fraction (F), the melting rate (Γ) and the bulk mineral/melt partition coefficients (D_i) 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₀).

Based on the formulations of Spiegelman and Elliott [12], Spiegelman [55] developed a numerical solution for ²³⁸U-²³⁰Th-²²⁶Ra disequilibria generated during chromatographic equilibrium porous flow melting. While Spiegelman and Elliott [12] 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. [19]. (As discussed in Sims et al. [19], however, determining this solution required correcting an error in equation 15 from Spiegelman and Elliott [12].) 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).

As originally derived in Spiegelman and Elliott [12], 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:

$${\mathrm{c}}_{\mathrm{i}}={\mathsf{\alpha }}_{\mathrm{i}}{\mathrm{c}}_{\mathrm{b}\mathrm{i}}$$

(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 given by the batch melting equation. Spiegelman and Elliott [12] 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 Elliot [12]). This approach, along with additional methods, have recently been implemented by Elkins and Spiegelman [56] and will be addressed in further detail below.

Figure 9. Results of chromatographic/equilibrium porous flow melting calculations for an adiabatically decompressing peridotitic mantle showing (A, C) (²²⁶Ra/²³⁰Th) and (B, D) (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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) and (B) show results for the analytical solution after Sims et al. (1999) assuming constant bulk solid/melt partitioning and melting rates, while (C) and (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. After Figure 7.

Figure 9. Results of chromatographic/equilibrium porous flow melting calculations for an adiabatically decompressing peridotitic mantle showing (A, C) (²²⁶Ra/²³⁰Th) and (B, D) (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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) and (B) show results for the analytical solution after Sims et al. (1999) assuming constant bulk solid/melt partitioning and melting rates, while (C) and (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. After Figure 7.

For the decay chain ²³⁸U-²³⁰Th-²²⁶Ra, the change in ${\alpha }_i$ with position in a 1D melt column for each isotope can be expressed by the following equations:

$${\displaystyle \frac{\mathrm{d}{\mathsf{\alpha }}_0}{\mathrm{d}\mathsf{\zeta }}}={-\mathsf{\lambda }}_0{\mathsf{\tau }}_0{\mathsf{\alpha }}_0$$

(24)

$${\displaystyle \frac{\mathrm{d}{\mathsf{\alpha }}_1}{\mathrm{d}\mathsf{\zeta }}}={\mathsf{\lambda }}_1({\mathsf{\tau }}_0{\mathsf{\alpha }}_0-{\mathsf{\tau }}_1{\mathsf{\alpha }}_1)$$

(25)

$${\displaystyle \frac{\mathrm{d}{\mathsf{\alpha }}_2}{\mathrm{d}\mathsf{\zeta }}}={\mathsf{\lambda }}_2({\mathsf{\tau }}_1{\mathsf{\alpha }}_1-{\mathsf{\tau }}_2{\mathsf{\alpha }}_2)$$

(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 }_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{i}}={\displaystyle \frac{{\mathsf{\rho }}_{\mathrm{f}}\omega +{\mathsf{\rho }}_{\mathrm{s}}\left(1-\right){\mathrm{D}}_{\mathrm{i}}\mathrm{W}}{{\mathsf{\rho }}_{\mathrm{f}}+{\mathsf{\rho }}_{\mathrm{s}}\left(1-\right){\mathrm{D}}_{\mathrm{i}}}}$$

(27)

where ${\rho }_f$ and ${\rho }_s$ are the melt and solid densities, respectively, ϕ 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 ={\displaystyle \frac{\mathsf{\Gamma }\mathrm{d}}{{\mathsf{\rho }}_{\mathrm{f}}}}$$

(28)

where $\mathsf{\Gamma }$ is the melting rate, which is assumed to be constant in the analytical approximation of Sims et al. (1999), and can be determined using the equation:

$$\mathsf{\Gamma }={\displaystyle \frac{\mathrm{W}{\mathsf{\rho }}_{\mathrm{s}}{\mathrm{F}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{d}}}$$

(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$ << 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 () 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 [19]. 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. [19] obtained the following (where the subscripts 0, 1, and 2 refer to ²³⁸U, ²³⁰Th, and ²²⁶Ra, respectively):

$${\mathsf{\alpha }}_0\left(\mathsf{\zeta }\right)={\mathsf{\alpha }}_0^0{\mathrm{e}}^{-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0\mathsf{\zeta }}$$

(30)

$${\mathsf{\alpha }}_1\left(\mathsf{\zeta }\right)={\displaystyle \frac{{{\mathsf{\lambda }}_1{\mathsf{\tau }}_0\mathsf{\alpha }}_0^0}{{\mathsf{\lambda }}_1{\mathsf{\tau }}_1-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0}}{\mathrm{e}}^{-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0\mathsf{\zeta }}+\left[{\mathsf{\alpha }}_1^0-{\displaystyle \frac{{{\mathsf{\lambda }}_1{\mathsf{\tau }}_0\mathsf{\alpha }}_0^0}{{\mathsf{\lambda }}_1{\mathsf{\tau }}_1-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0}}\right]{\mathrm{e}}^{-{\mathsf{\lambda }}_1{\mathsf{\tau }}_1\mathsf{\zeta }}$$

(31)

$$\begin{gathered} {\mathsf{\alpha }}_2\left(\mathsf{\zeta }\right)={\displaystyle \frac{{\mathsf{\lambda }}_2{{\mathsf{\lambda }}_1{{\mathsf{\tau }}_1\mathsf{\tau }}_0\mathsf{\alpha }}_0^0}{{(\mathsf{\lambda }}_1{\mathsf{\tau }}_1-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0){(\mathsf{\lambda }}_2{\mathsf{\tau }}_2-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0)}}{\mathrm{e}}^{-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0\mathsf{\zeta }}+{\displaystyle \frac{{\mathsf{\lambda }}_2{\mathsf{\tau }}_1}{{\mathsf{\lambda }}_2{\mathsf{\tau }}_2-{\mathsf{\lambda }}_1{\mathsf{\tau }}_1}}\left[{\mathsf{\alpha }}_1^0-{\displaystyle \frac{{{\mathsf{\lambda }}_1{\mathsf{\tau }}_0\mathsf{\alpha }}_0^0}{{\mathsf{\lambda }}_1{\mathsf{\tau }}_1-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0}}\right]{\mathrm{e}}^{-{\mathsf{\lambda }}_1{\mathsf{\tau }}_1\mathsf{\zeta }} \\ +\left[{\mathsf{\alpha }}_2^0-{\displaystyle \frac{{\mathsf{\lambda }}_2{{\mathsf{\lambda }}_1{{\mathsf{\tau }}_1\mathsf{\tau }}_0\mathsf{\alpha }}_0^0}{{(\mathsf{\lambda }}_1{\mathsf{\tau }}_1-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0){(\mathsf{\lambda }}_2{\mathsf{\tau }}_2-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0)}} \\ -{\displaystyle \frac{{\mathsf{\lambda }}_2{\mathsf{\tau }}_1}{{\mathsf{\lambda }}_2{\mathsf{\tau }}_2-{\mathsf{\lambda }}_1{\mathsf{\tau }}_1}}\left({\mathsf{\alpha }}_1^0-{\displaystyle \frac{{{\mathsf{\lambda }}_1{\mathsf{\tau }}_0\mathsf{\alpha }}_0^0}{{\mathsf{\lambda }}_1{\mathsf{\tau }}_1-{\mathsf{\lambda }}_0{\mathsf{\tau }}_0}}\right)\right]{\mathrm{e}}^{-{\mathsf{\lambda }}_2{\mathsf{\tau }}_2\mathsf{\zeta }} \end{gathered}$$

(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. [19], 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., ${∆\alpha }_{\mathrm{U}}$< ${∆\alpha }_{\mathrm{T}\mathrm{h}}$<< ${∆\alpha }_{\mathrm{R}\mathrm{a}}$ and ${∆\alpha }_{\mathrm{U}}$<< ${∆\alpha }_{\mathrm{P}\mathrm{a}}$). As such, one can hold the parental ${\alpha }_i$ values constant, further reducing Equations 31-32 to:

$${\mathsf{\alpha }}_1\left(\mathsf{\zeta }\right)={\mathsf{\alpha }}_1^0{\mathrm{e}}^{-{\mathsf{\lambda }}_1{\mathsf{\tau }}_1\mathsf{\zeta }}+{\displaystyle \frac{{\mathsf{\tau }}_0}{{\mathsf{\tau }}_1}}{\mathsf{\alpha }}_0\left[1-{\mathrm{e}}^{-{\mathsf{\lambda }}_1{\mathsf{\tau }}_1\mathsf{\zeta }}\right]$$

(33)

$${\mathsf{\alpha }}_2\left(\mathsf{\zeta }\right)={\mathsf{\alpha }}_2^0{\mathrm{e}}^{-{\mathsf{\lambda }}_2{\mathsf{\tau }}_2\mathsf{\zeta }}+{\displaystyle \frac{{\mathsf{\tau }}_1}{{\mathsf{\tau }}_2}}{\mathsf{\alpha }}_1\left[1-{\mathrm{e}}^{-{\mathsf{\lambda }}_2{\mathsf{\tau }}_2\mathsf{\zeta }}\right]$$

(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-34) provide nearly identical results (Figure 9). The reason for this similarity is that the half-lives of ²³⁰Th and ²²⁶Ra are very different, resulting in λ₁τ₁ << λ₂τ₂ (where λτ 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 ²³⁸U-²³⁰Th and ²³⁰Th-²²⁶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 which 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 [12]. These solutions show that (²²⁶Ra/²³⁰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 2x10^-2 kg·m^-3·yr^-1. In contrast, (²³⁰Th/²³⁸U) and (²³¹Pa/²³⁵U) disequilibria are controlled more strongly by the melting rate, which in turn is related primarily to the solid mantle upwelling rate (Figure 9, Figure 10 and Figure 11) and secondarily to the behavior of the melting lithology (see further discussion below).

3.1.1. 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 [12,55], whereas the analytical solutions assume a constant porosity over the length of the melting column [19], 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 (²²⁶Ra/²³⁰Th), while underestimating the solid mantle upwelling velocity, which more significantly influences (²³⁰Th/²³⁸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).

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., [21,22]) 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 for Darcy’s Law flow with compaction [12]; 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,59,60,61,62,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 amplify 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 10. Results of reactive porous flow partial melting calculations for an adiabatically decompressing peridotitic mantle at T_p = 1300 ºC showing (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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. After Figure 7.

Figure 10. Results of reactive porous flow partial melting calculations for an adiabatically decompressing peridotitic mantle at T_p = 1300 ºC showing (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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. After Figure 7.

Spiegelman [55] developed a numerical solution to the equilibrium porous flow melting formulations by Spiegelman & Elliott [12] 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 [12] 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).

Figure 11. Results of reactive porous flow partial melting calculations for an adiabatically decompressing peridotitic mantle at T_p = 1400 ºC showing (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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. After Figure 7.

Figure 11. Results of reactive porous flow partial melting calculations for an adiabatically decompressing peridotitic mantle at T_p = 1400 ºC showing (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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. After Figure 7.

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 (²³⁰Th/²³⁸U) as the analytical solution predicts for mid-ocean ridge environments [19,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 (²³⁰Th/²³⁸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 (²³⁰Th/²³⁸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, Figure 10 and Figure 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.

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 (²³⁰Th/²³⁸U) can be achieved by peridotitic melting. For chemical disequilibrium transport (i.e., fractional melting), constant porosity, partitioning, and melting rates overall generate lower (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵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 (²³⁰Th/²³⁸U) for equilibrium melting with constant input parameters, and again particularly low (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵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 (²²⁶Ra/²³⁰Th) disequilibrium is controlled mainly by the porosity of the melt region, and (²³⁰Th/²³⁸U) and (²³¹Pa/²³⁵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 [14], dynamic melting differs from chromatographic melting in that the produced melts are assumed to move instantly to the 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. [46,81]).

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].

Lithology	Phase	DU	DTh	DPa *	DRa *	Reference
Garnet Peridotite	Garnet	0.038	0.017	0.00001	0.00001	RD 1097-5 experiment, [73]
KG1 Pyroxenite	Clinopyroxene	0.0030	0.0040	0.00001	0.00001	RD 1097-5 experiment, [73]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Orthopyroxene	0.0078	0.0086	0.00001	0.00001	TM0500-3 experiment, [73]
	Pigeonite	0.00050	0.00018	0.00001	0.00001	W&P 122-1 experiment, [74]
	Spinel	0.012	0.0024	0.00001	0.00001	Lunar basalt, [75]
Spinel Peridotite	Cpx	0.0080	0.0070	0.00001	0.00001	TM 1094-9 experiment, [76]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Orthopyroxene	0.0024	0.0027	0.00001	0.00001	RD 1097-2 experiment, [73]
	Plagioclase	0.0006	0.0034	0.00001	0.02000	DU, DTh calculated after [54]; DRa from [77]
	Pigeonite	0.0005	0.00018	0.00001	0.00001	W&P 122-1 experiment, [74]
	Spinel	0.012	0.0024	0.00001	0.00001	Lunar basalt, [75]
G2 Pyroxenite	Garnet	0.02405	0.00415	0.00001	0.00001	A343 experiment, [78]
	Clinopyroxene	0.0041	0.0032	0.00001	0.00001	A343 experiment, [78]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Orthopyroxene	0.0078	0.0086	0.00001	0.00001	TM0500-3 experiment, [73]
	Plagioclase	0.0006	0.0034	0.00001	0.02000	DU, DTh calculated after [54]; DRa from [77]
	Pigeonite	0.0096	0.010	0.00001	0.00001	Experiment 18, [79]
	Spinel	0.046	0.016	0.00001	0.00001	Maximum measured, [80]
MIX1G, m7-16	Garnet	0.013	0.0032	0.00001	0.00001	Experimental results, [80]
Pyroxenite	Clinopyroxene	0.017	0.015	0.00001	0.00001	Experimental results, [80]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Spinel	0.046	0.016	0.00001	0.00001	Maximum measured, [80]
	Orthopyroxene	0.0078	0.0086	0.00001	0.00001	TM0500-3 experiment, [73]
	Plagioclase	0.0006	0.0034	0.00001	0.02000	DU, DTh calculated after [54]; DRa from [77]
	Pigeonite	0.0096	0.010	0.00001	0.00001	Experiment 18, [79]

* By convention, D_Pa and D_Ra are set equal to very small values (1 x 10^-5) for most mantle minerals, except D_Ra in plagioclase.

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].

Lithology	Phase	DU	DTh	DPa *	DRa *	Reference
Garnet Peridotite	Garnet	0.038	0.017	0.00001	0.00001	RD 1097-5 experiment, [73]
KG1 Pyroxenite	Clinopyroxene	0.0030	0.0040	0.00001	0.00001	RD 1097-5 experiment, [73]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Orthopyroxene	0.0078	0.0086	0.00001	0.00001	TM0500-3 experiment, [73]
	Pigeonite	0.00050	0.00018	0.00001	0.00001	W&P 122-1 experiment, [74]
	Spinel	0.012	0.0024	0.00001	0.00001	Lunar basalt, [75]
Spinel Peridotite	Cpx	0.0080	0.0070	0.00001	0.00001	TM 1094-9 experiment, [76]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Orthopyroxene	0.0024	0.0027	0.00001	0.00001	RD 1097-2 experiment, [73]
	Plagioclase	0.0006	0.0034	0.00001	0.02000	DU, DTh calculated after [54]; DRa from [77]
	Pigeonite	0.0005	0.00018	0.00001	0.00001	W&P 122-1 experiment, [74]
	Spinel	0.012	0.0024	0.00001	0.00001	Lunar basalt, [75]
G2 Pyroxenite	Garnet	0.02405	0.00415	0.00001	0.00001	A343 experiment, [78]
	Clinopyroxene	0.0041	0.0032	0.00001	0.00001	A343 experiment, [78]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Orthopyroxene	0.0078	0.0086	0.00001	0.00001	TM0500-3 experiment, [73]
	Plagioclase	0.0006	0.0034	0.00001	0.02000	DU, DTh calculated after [54]; DRa from [77]
	Pigeonite	0.0096	0.010	0.00001	0.00001	Experiment 18, [79]
	Spinel	0.046	0.016	0.00001	0.00001	Maximum measured, [80]
MIX1G, m7-16	Garnet	0.013	0.0032	0.00001	0.00001	Experimental results, [80]
Pyroxenite	Clinopyroxene	0.017	0.015	0.00001	0.00001	Experimental results, [80]
	Olivine	0.00005	0.00047	0.00001	0.00001	RD 1097-5 experiment, [73]
	Spinel	0.046	0.016	0.00001	0.00001	Maximum measured, [80]
	Orthopyroxene	0.0078	0.0086	0.00001	0.00001	TM0500-3 experiment, [73]
	Plagioclase	0.0006	0.0034	0.00001	0.02000	DU, DTh calculated after [54]; DRa from [77]
	Pigeonite	0.0096	0.010	0.00001	0.00001	Experiment 18, [79]

* By convention, D_Pa and D_Ra are set equal to very small values (1 x 10^-5) for most mantle minerals, except D_Ra in plagioclase.

In other words, in traditional dynamic melting [14], 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. [22], 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 [14], 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 [11,47], 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 [12] 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 calculate 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 [47]) and numerical modeling outcomes for dynamic melting (after [18,68,82,83]; Figure 8). The outcomes of the numerical approach are comparable to those for incremental models with very small step sizes [21,22,68]. Consistent with prior results, it is observed that dynamic melting outcomes in Figure 8 generally produce higher (²³⁰Th/²³⁸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,83]. 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 Figure 10 and Figure 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 ²³⁸U-²³⁰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. [20] 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 (²³⁴U/²³⁸U) in ocean waters ubiquitously is 1.14 owing to uranium’s long residence time [84,85,86,87]. Thus, any basalt that has been altered by seawater will have (²³⁴U/²³⁸U) greater than unity (Figure 12). The threshold of what the petrological community considers a tolerable excess of ²³⁴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 ²³⁵U on ²³⁴U can increase the measured (²³⁴U/²³⁸U). Hence, it remains uncertain whether (a) many of the studies showing ²³⁴U excesses in Figure 12 are an indication of imprecise (²³⁴U/²³⁸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 ²³⁴U excesses can indicate perturbation of (²³⁰Th/²³⁸U) in depleted basalts via alteration is the Kolbeinsey ridge (Figure 13) where there is a clear correlation between (²³⁰Th/²³⁸U) and (²³⁴U/²³⁸U) [43].

Because of the long half-life of ²³⁴U (t_1/2 = 245 kyr), for young subaerial samples, (²³⁴U/²³⁸U) of one is a necessary but insufficient condition for ruling out secondary alteration. In this case, other indices such as K₂O/P₂O₅ >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 ²²⁶Ra, which has a short half-life (1,600 ± 7 yr) relative to the expected timescales of melting, melt migration, and magma storage processes; and, since ²²⁶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 ²²⁶Ra, ²³⁰Th, and ²³¹Pa half-lives, if a lava sample has a primary (that is, melting-derived) ²²⁶Ra excess relative to ²³⁰Th, any age effects are small enough to be inconsequential for measured ²³⁸U-²³⁰Th or ²³⁵U-²³¹Pa disequilibria. Similarly, if a sample has a ²¹⁰Pb (t^½ = 22.6 ± 0.1 yr) excess or deficit relative to ²²⁶Ra, then age affects are not likely to be an issue for ²³⁰Th, ²³¹Pa, or ²²⁶Ra.

Figure 12. ²³⁴U/²³⁸U and ²³⁰Th/²³⁸U for the global MORB data set (Figure 7) showing the range of uncertainty. Note that many samples show ²³⁴U/²³⁸U indicative of seawater contamination.

Figure 12. ²³⁴U/²³⁸U and ²³⁰Th/²³⁸U for the global MORB data set (Figure 7) showing the range of uncertainty. Note that many samples show ²³⁴U/²³⁸U indicative of seawater contamination.

An important nuclide that could be used to address uncertainties associated with sample age is ²²⁷Ac, a daughter product of ²³¹Pa with a half-life of 21.77 yr, similar to ²¹⁰Pb. Unlike ²¹⁰Pb, which has intermediate volatile species such as ²²²Rn in the ²²⁶Ra-²¹⁰Pb decay scheme, which could be affected by gas exsolution and either loss or accumulation in the magma (Figure 2), ²²⁷Ac is a direct decay product of ²³¹Pa and is not volatile. However, because of the difficult analytical methods required to measure ²²⁷Ac [88], only a few studies [37,89,90] have reported ²³¹Pa-²²⁷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 (²²⁷Ac/²³¹Pa) of unity, suggesting that either i) the timescale of magma transport and storage is at least 100 years or ii) there was no initial fractionation of (²²⁷Ac/²³¹Pa).

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 (Supplemental Table 1). 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 ²³⁸U-²³⁰Th-²²⁶Ra disequilibria, the value chosen for D_Th is also particularly important, as ²³⁰Th is both a parent and daughter isotope and therefore its D value affects both (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰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 (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰Th) by up to a factor of 3 for low melting rates, and by an order of magnitude for high melting rates (e.g., [19]).

Figure 13. Uranium isotope calculated mixing line between seawater [(²³⁴U/²³⁸U) = 1.14] and unaltered local Kolbeinsey melt composition. This correlation based on high-precision mass spectrometric analyses demonstrates the need to have (²³⁴U/²³⁸U) in equilibrium within the uncertainty of modern mass spectrometric measurements.

Figure 13. Uranium isotope calculated mixing line between seawater [(²³⁴U/²³⁸U) = 1.14] and unaltered local Kolbeinsey melt composition. This correlation based on high-precision mass spectrometric analyses demonstrates the need to have (²³⁴U/²³⁸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., [80]). A compilation of experimentally determined D_U and D_Th 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_U and D_Th in turn can produce large variations in model outcomes (e.g., [19,43]), and, therefore, also in the calculated or inferred porosities and melting rates for both chromatographic and dynamic melting models (Figure 8, Figure 9, Figure 10 and Figure 11).

An additional uncertainty is the role of grain boundaries as an important repository for highly incompatible elements (e.g., [91,92,93]). 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 ([93], 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., [19]). 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., [94,95]). 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., [96,97]). 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., [97]). 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 modelled 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 modelled 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 [46].

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 [27], which, in turn, principally is controlled by oceanic crustal thickness [98]. This correlation with (²³⁰Th/²³⁸U) was attributed to the fact that the depth of the peridotite solidus is controlled by mantle potential temperature, which has been predicted as a controlling factor in producing varying crustal thicknesses and, by inference, basaltic geochemical compositions [98]. Using time-dependent melting models, Bourdon et al. [27] explained the broad global negative correlation between axial ridge depth and (²³⁰Th/²³⁸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 (²³⁰Th/²³⁸U) > 1 in the magma [58,59,60,61,63,73,76,99,100].

An important location to see this systematic relationship between (²³⁰Th/²³⁸U) and seafloor axial depth is the Arctic or North Atlantic mid-ocean ridge system [43,101] (Figure 14). The North Atlantic data set nearly spans the complete range of both axial depth (250–4131 m) and (²³⁰Th/²³⁸U) (0.955–1.296) for MORB globally. The high (²³⁰Th/²³⁸U) values measured in MORB from the Kolbeinsey [43], Mohns, and Knipovich Ridges [101] 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 [102], and major element data (e.g., high MgO and Na₈) suggest that they are the result of smaller degrees of melting [103,104]. The overall low (²³⁰Th/²³⁸U) measured in EVZ Gakkel MORB suggests a relatively shallow melt regime dominated by spinel peridotite [101]. The negative trend of axial depth versus (²³⁰Th/²³⁸U) initially observed and modeled by Bourdon et al. [27] is also considerably stronger for Arctic MORB than for the global data set. Elliott and Spiegelman [105] measured a weakly correlated r² coefficient of 0.58 for all global ridge data, while the Arctic ridges when taken alone [43,101] have (²³⁰Th/²³⁸U) values that are more strongly correlated with axial depth (r² = 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 (²³⁰Th/²³⁸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 [43]. 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 (²³⁰Th/²³⁸U) compositions and also the crustal thickness [101,106] (see below for further discussion of heterogeneity effects).

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 [20,26,34,40,43,44,107,108,109,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. [27] for various representative melting rates (M) and residual porosities (/) as labeled. See Elkins et al. [101] for full details.

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 [20,26,34,40,43,44,107,108,109,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. [27] for various representative melting rates (M) and residual porosities (/) as labeled. See Elkins et al. [101] for full details.

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 (²³⁸U/²³⁰Th) versus α_Sm/Nd (Figure 15; [8]); and, even more definitively, 2) a well-defined negative correlation exists between (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰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; [20]). This interpretation, that the magma source and transport regime is polybaric, was first proposed by Kelemen et al. [111] 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.

Early multiple-porosity interpretations (Figure 17 and Figure 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 [17,113,114]. In the Jull et al. [114] 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 15. ²³⁰Th excess versus alpha Sm/Nd for Hawaiian samples and axial 9-10°N and Siqueiros transform samples. Modified from Sims et al. [112]. 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 ¹⁴³Nd/¹⁴⁴Nd of the lavas (for details see [8,20]).

Figure 15. ²³⁰Th excess versus alpha Sm/Nd for Hawaiian samples and axial 9-10°N and Siqueiros transform samples. Modified from Sims et al. [112]. 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 ¹⁴³Nd/¹⁴⁴Nd of the lavas (for details see [8,20]).

Because the lavas from this correlated sample suite preserve relatively uniform radiogenic isotope signatures (¹⁴³Nd/¹⁴⁴Nd varies from 0.512378 ± 11 to 0.512400 ± 10, ⁸⁷Sr/⁸⁶Sr varies from 0.70244 ± 5 to 0.70257± 6, and ¹⁷⁶Hf/¹⁷⁷Hf varies from 0.283169 ± 5 to 0.283197 ± 5), source heterogeneity cannot be the only factor influencing ²³⁸U-²³⁰Th-²²⁶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 variable 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 homogeneous initial source [22,114], 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., [111,115], 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 heterogeneity (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., [116,117]) 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.

Figure 16. ²³⁰Th/²³⁸U versus ²²⁶Ra/²³⁰Th for axial 9-10°N samples showing inverse correlation. Modified from Sims et al. [20].

Figure 16. ²³⁰Th/²³⁸U versus ²²⁶Ra/²³⁰Th for axial 9-10°N samples showing inverse correlation. Modified from Sims et al. [20].

3.6.3. Effects of Source Heterogeneity on U-series Disequilibria in MORB

Because ²³⁰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 ²³⁸U-²³⁰Th, decompression melting of a heterogeneous mantle can explain much of the short-range and global compositional variability observed in MORB (e.g., [43,44,80,101,118]).

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., [36,71,119,120,121,122,123,124,125,126,127,128]). Recycled crustal material, which may be the principal carrier of isotopically enriched signatures in the mantle (e.g., [119]), 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,119]). 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,119]). 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., [119]).

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., [36,44,118,119,120,121,122,129], 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,119,123,124]). For example, for a decompressing mantle source region containing 5% of an eclogite-like (experimental mixture G2; [123,130]) pyroxenite and 95% average peridotite, estimates suggest that 14-24% of the total melt produced by mass will derive from the pyroxenite.

Figure 17. Results of one- and two-porosity partial melting calculations for (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰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. [114] 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. [114] 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 [12] 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. [22]. (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. [114], where f₁ indicates the porosities with depth for the low-porosity region, f₂ 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%.

Figure 17. Results of one- and two-porosity partial melting calculations for (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰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. [114] 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. [114] 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 [12] 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. [22]. (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. [114], where f₁ indicates the porosities with depth for the low-porosity region, f₂ 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., [36,121,122,131].

A systematic treatment of bi-lithologic mantle melting beneath mid-ocean ridges (Figure 19 and Figure 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 un-melted 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 ⁸⁷Sr/⁸⁶Sr and ²⁰⁶Pb/²⁰⁴Pb, less radiogenic ¹⁴³Nd/¹⁴⁴Nd and ¹⁷⁶Hf/¹⁷⁷Hf, and lower U/Th ratios record a melt source containing higher quantities of pyroxenite rocks.

Figure 18. Effects of end-member melting models – dynamic melting [10], chromatographic melting [12], two-porosity/polybaric porous flow [114], and dynamic melting using variable upwelling and melt extraction velocities [22]. The last two models both successfully explain the (²²⁶Ra/^230Th) and ²³⁰Th/²³⁸U) observed in the 9-10°N EPR data [20]. However, only the two-porosity/polybaric porous flow model [114] 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. [111].

Figure 18. Effects of end-member melting models – dynamic melting [10], chromatographic melting [12], two-porosity/polybaric porous flow [114], and dynamic melting using variable upwelling and melt extraction velocities [22]. The last two models both successfully explain the (²²⁶Ra/^230Th) and ²³⁰Th/²³⁸U) observed in the 9-10°N EPR data [20]. However, only the two-porosity/polybaric porous flow model [114] 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. [111].

Several studies (e.g., [44,101,106,109]) have found that MORB with relatively enriched isotopic compositions (i.e., relatively radiogenic ⁸⁷Sr/⁸⁶Sr, ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb, and unradiogenic ¹⁴³Nd/¹⁴⁴Nd and ¹⁷⁶Hf/¹⁷⁷Hf) exhibit both (a) elevated (²³⁰Th/²³⁸U) and (b) low, age-constrained (²³¹Pa/²³⁵U) and (²²⁶Ra/²³⁰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., [44,101,106,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 (²³⁰Th/²³⁸U) accompanied by low (²²⁶Ra/²³⁰Th) and (²³¹Pa/²³⁵U) observed in some MORB cannot be produced by peridotite melting, whereas lithologies with higher fusilibility and a greater stability range for garnet (i.e., garnet-bearing pyroxenites; e.g., [71,119], 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 19. Partial melting calculation results after Elkins & Lambart [68], showing (A) (²²⁶Ra/²³⁰Th) and (B) (²³¹Pa/²³⁵U) vs. (²³⁰Th/²³⁸U) in computed partial melts where mantle Tp = 1300 ºC, (C) (²²⁶Ra/²³⁰Th), and (D) (²³¹Pa/²³⁵U) vs. (²³⁰Th/²³⁸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 Dahmkö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. After Figure 13.

Figure 19. Partial melting calculation results after Elkins & Lambart [68], showing (A) (²²⁶Ra/²³⁰Th) and (B) (²³¹Pa/²³⁵U) vs. (²³⁰Th/²³⁸U) in computed partial melts where mantle Tp = 1300 ºC, (C) (²²⁶Ra/²³⁰Th), and (D) (²³¹Pa/²³⁵U) vs. (²³⁰Th/²³⁸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 Dahmkö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. After Figure 13.

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. [44] 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 ²³⁰Th excesses [20,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 ²³⁰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 [134]. 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 ²³⁰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. [134].

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 global MORB data (Figure 7) and peridotite melt modeling results for (²³⁰Th/²³⁸U) and axial depth after Bourdon et al. [27] for comparison. 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,132] 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. [133].

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 global MORB data (Figure 7) and peridotite melt modeling results for (²³⁰Th/²³⁸U) and axial depth after Bourdon et al. [27] for comparison. 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,132] 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. [133].

Waters et al. [44] 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 [134]. Waters et al. [44] 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. [134]. Waters et al. [44] 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. [135]) 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 [44].

At intermediate spreading rates, the compositions of basalts from the Southeast Indian Ridge (SEIR) vary along the ridge axis in both (²³⁰Th/²³⁸U) and (²²⁶Ra/²³⁰Th) [109]. The data exhibit an inverse correlation between (²²⁶Ra/²³⁰Th) and (²³⁰Th/²³⁸U) similar to that of the global array [20], but (²³⁰Th/²³⁸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 (²²⁶Ra/²³⁰Th) and (²³⁰Th/²³⁸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 [43,101]. 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., [101,136,137]). For U-series disequilibria, the high-temperature but likely peridotite-dominated system beneath the Kolbeinsey Ridge produces elevated (²³⁰Th/²³⁸U) and (²³¹Pa/²³⁵U), while lavas from the Mohns and Knipovich Ridges have relatively high (²³⁰Th/²³⁸U) but low age-constrained (²³¹Pa/²³⁵U) and (²²⁶Ra/²³⁰Th), best explained by melting a pyroxenite-bearing mantle source [43,101].

In another slow-spreading ridge example, MORB from five second-order segments between 5-11º S along the southern Mid-Atlantic Ridge (MAR) exhibit a wide range of ²³⁰Th excesses from 0 to 32%, as well as heterogeneous trace element and radiogenic isotope compositions [138] and axial depths. Notably, however, (²³⁰Th/²³⁸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. [106,139] evaluated ²³⁸U-²³⁰Th-²²⁶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 [139]. 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 [101,106].

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,101,106]. 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 [114,140,141].

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 [8,9,19,142,143,144,145]. 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 [19,99,146,147,148,149,150,151,152]. 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 [19,142,143,145,153].

3.7.1. 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 about 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 [8,19,154].

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., [155,156]). This variation is believed 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 ²³⁸U-²³⁰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 [157,158,159,160] and more recent mass spectrometric studies using TIMS, SIMS, and MC-ICPMS [8,19,37,154,161,162]. Most U-series studies have focused on the active shield-stage tholeiitic volcanoes, Kilauea and Mauna Loa [8,19,154,157,158,159,160,161,162,163,164]. However, ²³⁸U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹Pa disequilibria have also been measured in young tholeiitic and alkaline lavas from the preshield-stage submarine volcano Kamaʻehuakanaloa [8,19,154] and in alkaline lavas from the post-shield stage volcanoes, Hualalai, and Mauna Kea, and the rejuvenated stage volcano Haleakala [8,19,37]. ²³⁸U-²³⁰Th disequilibria further have been measured in young alkaline submarine flows from the South Arch [8].

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 ²³⁸U-²³⁰Th and ²³⁰Th-²²⁶Ra occurred during the formation of Kilauea magmas on timescales less than 375 ka and 8 ka, respectively. However, the origin of this fractionation, whether melting or secondary processes such as assimilation of altered oceanic crust, was highly uncertain [157]. Subsequent mass spectrometric studies with higher-precision measurements have consistently demonstrated that the ²³⁰Th excesses in tholeiitic lavas from Kamaʻehuakanaloa, Kilauea, and Mauna Loa are small (2 to 6%) [8,19,165] and constant [162] compared to, e.g., the large ²³⁰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 ²³⁰Th excesses for Mauna Loa, but as shown by Sims et al. [19], 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 ²³⁰Th excesses than Hawaiian tholeiites, ranging from 15 to 30% [8,19,37]. The observation that these ²³⁰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 (²³⁰Th/²³⁸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 [8]. This correlation provided definitive evidence that U-Th fractionation is related to the melting process, at a time when the origins of (²³⁰Th/²³⁸U) disequilibria in magmatic systems, particularly in MORB, were debated. Additionally, this conclusion was contemporaneous with early partitioning studies [58,59,61] showing that (²³⁰Th/²³⁸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., [166]), 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 (²³⁰Th/²³⁸U), (²²⁶Ra/²³⁰Th) disequilibria in Hawaiian shield stage tholeiites are much larger, with ²²⁶Ra excesses ranging from 7 to 15% [19,160,161,164,165]. For alkalic lavas from Hualalai and Haleakala, ²²⁶Ra excesses are even greater, up to 30% [19,37]. In Hawaiian lavas, (²³¹Pa/²³⁵U) suggests even greater extents of parent-daughter fractionation, with ²³¹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ā [19,167]. Like (²³⁰Th/²³⁸U), the magnitude of the ²³⁵U-²³¹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 ²³⁸U-²³⁰Th and ²³⁵U-²³¹Pa disequilibria in Hawaiian lavas are mainly sensitive to melt fraction because of the long half-lives of ²³⁰Th and ²³¹Pa [8,9,19], interpretation of ²²⁶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., [12,14]). Sims et al. [19] showed that ²³⁸U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹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 ²³⁸U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹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 [150] 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 [18] used the plume models developed by Watson and McKenzie [150] 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. [19] 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 [150] and Hauri et al. [99], 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 [168], which considers the influence of the overriding Pacific Plate on the plume’s upwelling dynamics. Later studies [37,165,169] support the interpretations of Sims et al. [19] with further measurements of Kilauea, Kamaʻehuakanaloa, and Haleakala.

Subsequent work by Kokfelt et al. [143] and Bourdon et al. [153] showed similar variations in the activities of (²³⁰Th/²³⁸U) and (²³¹Pa/²³⁵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 [19,144].

U-series disequilibria to understand global hotspot buoyancy flux

Chabaux and Allègre [170] 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 [28,144,153,171].

Bourdon et al. [153] 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. [153], the range of buoyancy fluxes (B) found in mantle plumes is 0.5 to 1 Mg s^-1 for hotspots like Afar and Iceland and up to between 3.5 and 9 Mg s^-1 for Hawaii, depending on how the buoyancy flux is calculated (Ribe [172] or Sleep [171], 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) [153].

Figure 21. Solid mantle upwelling rate calculated from ²³⁸U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹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 [173]. Also shown are the model curves for solid mantle upwelling from Watson and McKenzie [150], Hauri et al. [99], and Jull and Ribe [174]. The model of Watson and McKenzie [150] 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 [174] take into account the effect of motion of the Pacific plate and depletion buoyancy on the dynamics of the plume, whereas the calculations of Hauri et al. [99] and Watson and McKenzie [150] do not. The ²³⁰Th and ²³¹Pa data show a clear asymmetry in upwelling velocity, which is evident in the results of Jull and Ribe [174]. 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.

Figure 21. Solid mantle upwelling rate calculated from ²³⁸U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹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 [173]. Also shown are the model curves for solid mantle upwelling from Watson and McKenzie [150], Hauri et al. [99], and Jull and Ribe [174]. The model of Watson and McKenzie [150] 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 [174] take into account the effect of motion of the Pacific plate and depletion buoyancy on the dynamics of the plume, whereas the calculations of Hauri et al. [99] and Watson and McKenzie [150] do not. The ²³⁰Th and ²³¹Pa data show a clear asymmetry in upwelling velocity, which is evident in the results of Jull and Ribe [174]. 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.

However, it should be noted that the Samoan data presented in Sims and Hart [175] and Sims et al. [35] 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 [171,172]. These Samoan data change the apparent form of the Bourdon et al. [153] correlation, making it much more hyberbolic 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 (the Tonga subduction) [35,175].

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 [171,172] 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 [19,172]. For example, in hotspots like Hawaii where the motion of the overriding Pacific plate influences the upwelling rate [19,144,171], 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. [153]. 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 [168], 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 [176]. 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 (²⁰⁶Pb/²⁰⁴Pb; ²⁰⁷Pb/²⁰⁴Pb and ²⁰⁸Pb/²⁰⁴Pb) and producing what is now known as the mantle tetrahedron (e.g. [177,178,179], 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., [180,181]), and more sophisticated statistical treatments have explored what structures are contained in the multi-dimensional global data set [182].

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 ²³⁸U- ²³⁰Th isotopic system is intrinsically different in two important ways: first, the half-life of daughter isotope ²³⁰Th is only 75 ka; secondly, the abundance of ²³⁰Th relative to parent ²³⁸U is determined by the radioactive principles governing U decay-series systematics (Section 2 above). As a result, after about 375 ka, the Th isotopic composition of a mantle source for partial melts directly reflects its ²³⁸U/²³²Th source ratio. Thus, comparison of (²³⁰Th/²³²Th), (²³⁸U/²³²Th), and (²³⁰Th/²³⁸U) with longer-lived radiogenic isotopes (e.g., ⁸⁷Sr/⁸⁶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 [7], who compared a small number of U-series measurements to ⁸⁷Sr/⁸⁶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 (²³⁰Th/²³²Th) and ⁸⁷Sr/⁸⁶Sr than between (²³⁸U/²³²Th) and ⁸⁷Sr/⁸⁶Sr. They concluded that (²³⁰Th/²³²Th) in a lava better represents the original mantle source (²³⁸U/²³²Th) composition compared to (²³⁸U/²³²Th), which they posited was fractionated by processes in which Th is more incompatible than U. As such, they suggested that (²³⁸U/²³⁰Th) fractionation was due to time-independent melting, and that all fractionation on a U-Th isochron diagram was due to a horizonal displacement without any change in (²³⁰Th/²³²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 [175] examined the relationship between (²³⁰Th/²³²Th), (²³⁸U/²³²Th), and the long-lived isotope ratios of Sr, Nd, and Pb, all measured by mass spectrometric methods (Figure 22). Sims and Hart [175] quantitively demonstrated that: 1) the relationships between (²³⁰Th/²³²Th), (²³⁸U/²³²Th), and the long-lived isotope ratios of ⁸⁷Sr/⁸⁶Sr and ¹⁴³Nd/¹⁴⁴Nd can be approximated by two-component mixing; 2) (²³⁰Th/²³²Th) is somewhat better correlated with ⁸⁷Sr/⁸⁶Sr and ¹⁴³Nd/¹⁴⁴Nd than ²³⁸U/²³²Th, and 3) the global basalt arrays of (²³⁰Th/²³²Th) and (²³⁸U/²³²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).

Figure 22. ²³⁰Th/²³²Th versus ⁸⁷Sr/⁸⁶Sr diagram for mid-ocean ridge, ocean island, and continental basalts compiled from the literature [8,13,19,20,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,90]. 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 ²³⁰Th/²³²Th ratios) is greater than the Th/Sr ratio in the depleted component.

Figure 22. ²³⁰Th/²³²Th versus ⁸⁷Sr/⁸⁶Sr diagram for mid-ocean ridge, ocean island, and continental basalts compiled from the literature [8,13,19,20,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,90]. 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 ²³⁰Th/²³²Th ratios) is greater than the Th/Sr ratio in the depleted component.

3.9. Average Depths of Mantle Melting from (²³¹Pa/²³⁵U) and (²³⁰Th/²³⁸U) in MORB and OIB

Correlations between (²³¹Pa/²³⁵U) and (²³⁰Th/²³⁸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 ²³⁵U-²³¹Pa systematics, considerable information can be gained by incorporating measurements of ²³⁵U-²³¹Pa [19,28,145,153,167,183]. Pickett and Murrell [167] observed that in plots of (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸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 [19,28,145,153]. 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 [28,35,144]. 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 ²³¹Pa excess at a given (²³⁰Th/²³⁸U).

In a follow-up assessment of the correlations observed in plots of (²³¹Pa/²³⁵U) versus (²³⁰Th/²³⁸U) in OIB, and to assess the mechanisms that control ²³⁸U-²³⁰Th and ²³⁵U-²³¹Pa in lavas, Prytulak and Elliot [145] measured ²³⁸U-²³⁰Th and ²³⁵U-²³¹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 ((²³⁰Th/²³⁸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; [19]), and Samoa alkaline basalts ((²³⁰Th/²³⁸U) = 1.03-1.16 and (²³¹Pa/²³⁵U) = 1.05-1.33; [35]). 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 [145]. However, simple batch melting models predict an extent of melting that is unreasonably small (~0.3%). Prytulak and Elliott [145] 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.

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 [12,14]. In all cases, MORB, because of their low upwelling rates during melting, require daughter nuclide ingrowth to produce ²³⁸U-²³⁰Th, ²³⁵U-²³¹Pa, and ²³⁰Th-²²⁶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 ²³⁰Th and ²³¹Pa and where it is highly likely that elemental fraction of Th/U and Pa/U is the dominant mechanism for ²³⁸U-²³⁰Th and ²³⁵U-²³¹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 ²³⁸U-²³⁰Th and ²³⁵U-²³¹Pa disequilibria [145]. 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 (²³⁸U/²³⁰Th) and (²³⁵U/²³¹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 [9]. In contrast, because of the much shorter half-life of ²²⁶Ra, ingrowth is necessary to produce ²³⁰Th-²²⁶Ra disequilibria in all cases.

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_U, D_Th, and D_Pa values and respective color schemes for the different lines are the same as given in Figure 20, except for D_Ra, 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 [144], 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 ²³¹Pa excesses at a given ²³⁰Th/²³⁸U. In (B), the ²²⁶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 ²²⁶Ra are from 9-10°N EPR [20]. These MORB samples show an inverse correlation between ²³⁰Th excess, ²²⁶Ra excess (as well as major and trace element compositions) which was interpreted as mixing of two types of melt: one having a high-²³⁰Th component equilibrated with deep, undepleted garnet peridotites (> 70 kms); the other having a high-²²⁶Ra component preserving characteristics indicative of equilibration with highly depleted, residual harzburgite in the uppermost spinel-facies mantle (< 70 kms) [114].

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_U, D_Th, and D_Pa values and respective color schemes for the different lines are the same as given in Figure 20, except for D_Ra, 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 [144], 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 ²³¹Pa excesses at a given ²³⁰Th/²³⁸U. In (B), the ²²⁶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 ²²⁶Ra are from 9-10°N EPR [20]. These MORB samples show an inverse correlation between ²³⁰Th excess, ²²⁶Ra excess (as well as major and trace element compositions) which was interpreted as mixing of two types of melt: one having a high-²³⁰Th component equilibrated with deep, undepleted garnet peridotites (> 70 kms); the other having a high-²²⁶Ra component preserving characteristics indicative of equilibration with highly depleted, residual harzburgite in the uppermost spinel-facies mantle (< 70 kms) [114].

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 about three hundred and seventy-five thousand years. 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. ${\mathsf{\lambda }}_1{\mathrm{N}}_1={\mathsf{\lambda }}_2{\mathrm{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. ²³⁰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., ²²⁶Ra, t_1/2 = 1.6 kyr), or even longer than, five times the isotope’s half-life (e.g., ²²⁸Ra, t_1/2 = 5.77 a, or ²²²Rn, t_1/2 = 3.85 d).

Figure 24. Schematic of decay of ²³⁰Th, ²²⁶Ra, and ²¹⁰Pb in the melt and ingrowth in the mantle residue after a fractionation event occurring at t₁. Inset shows the short-lived excess of ²¹⁰Pb. Approximate valid range of U-Th and Th-Ra dating methods are given at the bottom in kyr.

Figure 24. Schematic of decay of ²³⁰Th, ²²⁶Ra, and ²¹⁰Pb in the melt and ingrowth in the mantle residue after a fractionation event occurring at t₁. Inset shows the short-lived excess of ²¹⁰Pb. Approximate valid range of U-Th and Th-Ra dating methods are given 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 Ma for the ²³⁸U-series, (based on five times the half-life of ²³⁴U); ~200 ka for the ²³⁵U-series (based on five times the half-life of ²³¹Pa); and ~35 a for the ²³²Th-series (based on five times the half-life of ²²⁸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 ²³⁸U-²⁰⁶Pb, ²³⁵U-²⁰⁷Pb, and ²³²Th-²⁰⁸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 ²³⁸U, ²³⁴U, and ²³⁰Th or ²³⁵U and ²³¹Pa can measurably affect ages for geochronometers involving ²⁰⁶Pb or ²⁰⁷Pb near the lower limit of their useful age ranges (e.g., [184]).

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 (Figure 2 and Figure 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 ²³⁸U decay chain, namely ²³⁸U-²³⁰Th, ²³⁰Th-²²⁶Ra, ²²⁶Ra-²¹⁰Pb, and ²¹⁰Po-²¹⁰Pb. A few studies have used ²³⁵U-²³¹Pa, and even fewer have used ²³²Th-²²⁸Ra or ²³¹Pa-²²⁷Ac (see Sims et al. [89] 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: (i) the presence or absence of U-series disequilibria to establish age limits based on the half-lives of the daughter isotopes; (ii) an internal isochron method, commonly obtained on mineral separates; (iii) 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 iv) 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 ¹⁴⁷Sm-¹⁴³Nd, ⁸⁷Rb-⁸⁷Sr, ²³⁸U-²⁰⁶Pb, and ¹⁷⁶Lu-¹⁷⁶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 [185]; 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., [186]). Most notable and transformative was the development of mass spectrometric methods for measuring Th isotopes in basalts [187,188,189]. 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 [20,33,35,41,44,45,46,81,90,107,108,110,186,187,189,190,191,192,193,194,195,196].

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 ²³⁸U-²³⁰Th-²²⁶Ra-²¹⁰Pb have half-lives ranging from ~20 years to 75 ka, these techniques are appropriate for dating MORB erupted from 0.1-375 ka (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. [197]). 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 ²³⁸U-²³⁰Th, ²³⁵U-²³¹Pa, and ²³⁰Th-²²⁶Ra disequilibria; in that case, ²³⁰Th-²²⁶Ra isotopic systematics are most constraining, because ²²⁶Ra has the shortest half-life, and require the lava to be less than 8 ka, providing it has remained a closed system since eruption. Suppose another sample of similar composition from the same field has ²³⁸U-²³⁰Th and ²³⁵U-²³¹Pa disequilibria, but no ²³⁰Th-²²⁶Ra disequilibrium; in that case, the sample is less than 160 ka and is ‘likely’ older than 8 ka. Finally, all samples with no ²³⁸U-²³⁰Th disequilibria are ‘likely’ >375ka, and so forth (Table 4). Examples where multiple nuclides are used in concert to constrain sample ages include Goldstien et al. [189], Sims et al. [35,81], and Waters et al. [45].

Table 4. Age limits of U-series radionuclides.

Nuclide Pair	Age Constraint
210Po/210Pb ≠ 1	< 2 yr
210Pb/226Ra ≠ 1	< 100 yr
226Ra/230Th ≠ 1	< 8 kyr
231Pa/235U ≠ 1	< 160 kyr
230Th/238U ≠ 1	< 375 kyr

Table 4. Age limits of U-series radionuclides.

Nuclide Pair	Age Constraint
210Po/210Pb ≠ 1	< 2 yr
210Pb/226Ra ≠ 1	< 100 yr
226Ra/230Th ≠ 1	< 8 kyr
231Pa/235U ≠ 1	< 160 kyr
230Th/238U ≠ 1	< 375 kyr

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 ²²⁶Ra-²³⁰Th disequilibria, it is reasonable to assume that the ²³⁰Th-²³⁸U disequilibria has not been significantly changed by decay and can be interpreted in terms of melting and melt transport processes (e.g. [20,35,43,107,114,186,198]). Similarly, if a solidified lava has significant ²¹⁰Pb-²²⁶Ra disequilibria, then measurements of ²²⁶Ra-²³⁰Th disequilibria in that lava can also be reasonably interpreted in terms of magmatic processes [35,45,81,101,190,199,200,201].

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 [7,164,198,202,203,204,205,206,207,208].

Figure 25. Schematic of internal isochrons in (²³⁰Th/²³²Th) versus (²³⁸U/²³²Th) space. Isochrons shown for a zero-age lava and for lavas with ages of one and two half-lives of ²³⁰Th. Arrows show directions of Th and U enrichment and subsequent decay toward the equiline, where the activities of ²³⁰Th and ²³⁸U are in equilibrium.

Figure 25. Schematic of internal isochrons in (²³⁰Th/²³²Th) versus (²³⁸U/²³²Th) space. Isochrons shown for a zero-age lava and for lavas with ages of one and two half-lives of ²³⁰Th. Arrows show directions of Th and U enrichment and subsequent decay toward the equiline, where the activities of ²³⁰Th and ²³⁸U are in equilibrium.

The U-series isochron method is akin to the isochron method applied to long-lived radiogenic systems such as ⁸⁷Rb/⁸⁶Sr-⁸⁷Sr/⁸⁶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 ⁸⁷Rb-⁸⁷Sr, the U-series Bateman equation of ingrowth and decay (see Eq. 8) is linear when normalized to a stable element and plotted in the appropriate coordinates.

For the ²³⁸U-²³⁰Th system, ²³²Th, which has a half-life of 14.01 Gyr, is effectively stable over the half-life of ²³⁰Th (<375 Kyr) and is used as a normalizing isotope. For the ²³⁰Th-²²⁶Ra system, there is no stable, normalizing Ra isotope for ²²⁶Ra. Albeit imperfect, stable Ba is typically used as a Ra proxy. ²³⁵U-²³¹Pa also lacks a stable, normalizing isotope for ²³¹Pa, and an internal mineral isochron method to the ²³⁵U-²³¹Pa system has yet to be used. Lead concentrations can be used as a normalizing factor for a ²¹⁰Pb-²²⁶Ra geochronometer [209]. 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 ²¹⁰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 (²³⁸U)/(²³²Th)-(²³⁰Th)/(²³²Th) and (²³⁰Th)/Ba-(²²⁶Ra)/Ba.

4.2.1.(238. U)/(²³²Th) - (²³⁰Th)/(²³²Th)

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

$$\mathrm{t}={\displaystyle \frac{\mathrm{l}\mathrm{n}\left(1-\mathrm{m}\right)}{-\mathsf{\lambda }}}$$

(35)

The assumptions for the (²³⁸U)/(²³²Th)-(²³⁰Th)/(²³²Th) isochron method are the same as for most other isotope systems, namely that the minerals are cogenetic and coeval, have uniform (²³⁸U/²³⁰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 ²³⁰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 ²³⁰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 [208], 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 ³He concentrations (Figure 28) provide a unique and self-consistent quantification of both the eruption age and the erosion rate [208].

There are numerous examples of ²³⁸U-²³⁰Th mineral isochrons in arc and continental systems (e.g. [204]), but only a few in oceanic settings [108]. This dearth of ²³⁸U-²³⁰Th isochrons for MORB is because most MORB samples are phenocryst-poor, with the exception being some glasses that are plagioclase-phyric.

Figure 26. Three simplified, potential (²³⁸U/²³²Th)-(²³⁰Th/²³²Th) isochron scenarios. (A) Polymineralic 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 ²³⁰Th, and then erupts soon thereafter as a porphyritic lava. The ²³⁰Th in the different solidified mineral phases of the lava then decays toward equilibrium for one half-life. (B) Polymineralic 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 ²³⁰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.

Figure 26. Three simplified, potential (²³⁸U/²³²Th)-(²³⁰Th/²³²Th) isochron scenarios. (A) Polymineralic 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 ²³⁰Th, and then erupts soon thereafter as a porphyritic lava. The ²³⁰Th in the different solidified mineral phases of the lava then decays toward equilibrium for one half-life. (B) Polymineralic 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 ²³⁰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.

4.2.2.(230. Th)/Ba - (²²⁶Ra)/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 ²³⁰Th-²²⁶Ra disequilibria chronometer. Use of mineral-glass ²²⁶Ra-²³⁰Th disequilibria requires knowledge of the initial ²²⁶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., [210,211,212]).

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,213]. Not accounting for this difference will produce significant inaccuracies in calculated ages. Furthermore, Ra is preferentially excluded from these minerals, leading to lower (²²⁶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 (²²⁶Ra)/Ba due to the aging of crystals.

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 nonlinear, maximum-likelihood method to estimate the slope and intercept of the best-fit line and Monte Carlo methods to estimate 95% confidence limits [214]. 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 (²³⁸U/²³²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 [215], 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.

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 nonlinear, maximum-likelihood method to estimate the slope and intercept of the best-fit line and Monte Carlo methods to estimate 95% confidence limits [214]. 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 (²³⁸U/²³²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 [215], 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. [164] developed a technique to correct ²²⁶Ra-²³⁰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 1,000 ± 350 yrs. (or given the lava’s age, 950 ± 350 yrs. before 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 [207,216]. These studies demonstrated that ²²⁶Ra-²³⁰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 [89,217].

Figure 28. ³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 dotted 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 ³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 [218] 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.

Figure 28. ³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 dotted 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 ³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 [218] 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.

Additionally, it is also important to note that because the timescales of crystallization are comparable to the half-life of ²²⁶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 ²²⁶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 ²²⁶Ra/²³⁰Th in mineral phases and host glasses can lead to erroneous results was demonstrated by Sims et al. [89] on recently erupted lava bombs from Mount Erebus volcano (Figure 29). The calculated ²²⁶Ra-²³⁰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 [89,90]. While these ages may seem reasonable, measurements of (²²⁸Ra/²³²Th) in these same crystals were out of equilibrium, which is highly problematic for the ²²⁶Ra-²³⁰Th isochron ages given that ²²⁸Ra has a half-life of only 5.77 years. To resolve this discrepancy, Sims et al. [89] 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 29. Crystal ages for Mount Erebus anorthoclase megacrysts [89]. (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 Ca map), apatite (second brightest in BSE, highlighted in red on 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_Ra/D_Ba showing how the age of crystallization changes as a function of D_Ra/D_Ba. Calculated D_Ra/D_Ba values for Or₁₈ and An₁₇ are taken from Blundy & Wood [54]. The value of D_Ra/D_Ba 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.Time is since 2012. (²²⁸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. Reagan et al. (1992) data are shown for comparison. (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_1/2 = 22· 6 years) and ²²⁷Ac (t_1/2 = 21·77 years) (not shown). However, the observation that (²²⁸Ra/²³⁰Th) (t1_/2= 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 ²³⁸U, ²³⁵U, and ²³²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; [219]). Modified after Sims et al. [89].

Figure 29. Crystal ages for Mount Erebus anorthoclase megacrysts [89]. (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 Ca map), apatite (second brightest in BSE, highlighted in red on 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_Ra/D_Ba showing how the age of crystallization changes as a function of D_Ra/D_Ba. Calculated D_Ra/D_Ba values for Or₁₈ and An₁₇ are taken from Blundy & Wood [54]. The value of D_Ra/D_Ba 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.Time is since 2012. (²²⁸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. Reagan et al. (1992) data are shown for comparison. (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_1/2 = 22· 6 years) and ²²⁷Ac (t_1/2 = 21·77 years) (not shown). However, the observation that (²²⁸Ra/²³⁰Th) (t1_/2= 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 ²³⁸U, ²³⁵U, and ²³²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; [219]). Modified after Sims et al. [89].

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 ²³⁰Th model ages, and by hundreds to thousands of years for ²²⁶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:

• 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.
• 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 ²³⁰Th, ²³¹Pa, and ²²⁶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,
• 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., [20,45,46,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 (Figure 4 and Figure 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:

$$\mathrm{t}=-{\displaystyle \frac{1}{{\mathsf{\lambda }}_{230}}}\mathrm{l}\mathrm{n}\left[{\displaystyle \frac{{\left(\frac{{}^{230}\mathrm{T}\mathrm{h}}{{}^{232}\mathrm{T}\mathrm{h}}\right)}_{\mathrm{m}}-{\left(\frac{{}^{238}\mathrm{U}}{{}^{232}\mathrm{T}\mathrm{h}}\right)}_{\mathrm{m}}}{{\left(\frac{{}^{230}\mathrm{T}\mathrm{h}}{{}^{232}\mathrm{T}\mathrm{h}}\right)}_0-{\left(\frac{{}^{238}\mathrm{U}}{{}^{232}\mathrm{T}\mathrm{h}}\right)}_{\mathrm{m}}}}\right]$$

(36)

and

$$\mathrm{t}=-{\displaystyle \frac{1}{{\mathsf{\lambda }}_{226}}}\mathrm{l}\mathrm{n}\left[{\displaystyle \frac{{\left(\frac{{}^{226}\mathrm{R}\mathrm{a}}{{}^{230}\mathrm{T}\mathrm{h}}\right)}_{\mathrm{m}}-1}{{\left(\frac{{}^{226}\mathrm{T}\mathrm{h}}{{}^{233}\mathrm{T}\mathrm{h}}\right)}_0-1}}\right]$$

(37)

where the subscript m represents the sample’s measured activity ratio, the subscript 0 is its initial activity ratio upon eruption, and λ₂₃₀ and λ₂₂₆ are the decay constants for ²³⁰Th and ²²⁶Ra, respectively. In Equation 36, ²³²Th is used as a stable analog for ²³⁰Th [5,7]. The different form of Equation 37 is used because no stable analog exists for ²²⁶Ra. The equation for determining ²³¹Pa model ages is of a similar form to the equation for ²²⁶Ra.

In these equations, the one additional unknown, besides time, is the initial activity ratio [i.e. (²³⁰Th/²³²Th)₀, (²²⁶Ra/²³⁰Th)₀, and (²³¹Pa/²³⁵U)₀]. 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 [20,45,46,107,186,187,188,189,190,195,196]. 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 [20,45,46,107,188,189,190,195,196] 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, ²³⁸U-²³⁰Th-²²⁶Ra dating evolved to include both larger sample suites of well-located and geologically well-constrained samples, as well as measurements of ²²⁶Ra, which provided finer time resolution. These spatial and temporal geological constraints for the axial lavas were also essential to the application of (²²⁶Ra/²³⁰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 (²³⁰Th/²³²Th)₀, (²²⁶Ra/²³⁰Th)₀, and (²³¹Pa/²³⁵U)₀ were all known explicitly for the magmatic system [20].

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 [20,45,46,107,195,196,220], 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 [220]. Similar systematics were also suggested for samples from the Southweast 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.

4.4. Dating Young Volcanic Rocks with ²²⁶Ra-²¹⁰Pb

While there have been numerous measurements of (²¹⁰Pb/²²⁶Ra) in young volcanic arc lavas [191,209,221,222,223,224], there are relatively few such measurements in MORB (e.g. [43,46,81,101,106,191]) and in oceanic intra plate settings (e.g. [35,163,225,226]).

Figure 30. (A) (²³⁸U/²³²Th) versus (²³⁰Th/²³²Th) for the 9°50’N off-axis samples compared to the 9-10°N EPR axial samples [20,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 (²³⁰Th/²³²Th) and (²³⁸U/²³²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. Inset provides regional context and shows the 9-10°N EPR axial samples [20,107] compared with the Siqueiros Transform samples [20,34]. Also shown are the near-axial basalts from 9°37’N EPR [187]. (B) ²³⁰Th/²³⁸U versus ²²⁶Ra/²³⁰Th for 9°50’N EPR off-axis samples compared with the 9°50’N EPR axial samples [20,107]. Error bars represent two sigma analytical uncertainties. Because of the significant difference in the half-lives of ²²⁶Ra and ²³⁰Th, (²²⁶Ra/²³⁰Th) decays back to equilibrium before (²³⁰Th/²³⁸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 (²²⁶Ra/²³⁰Th) in equilibrium, indicating that their ages are older (> 8 ka) than can be determined through ²³⁰Th-²²⁶Ra disequilibria. Inset provides regional context and shows (²³⁰Th/²³⁸U) versus (²²⁶Ra/²³⁰Th) for the 9°50’N EPR off-axis samples from this study compared with the 9-10°N EPR axial samples [20,107], three “near-axis” 9°30N EPR samples, including one E-type from close to the 9°37’N DEVAL [192], and the Siqueiros transform samples [34]. (C) 9°50’N EPR samples (²³⁰Th/²³⁸U) versus their distance from the center of the ridge axis, or Axial Summit Trough. Decay curves show how (²³⁰Th/²³⁸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 (²³⁰Th/²³⁸U) for the 9°50’N EPR axial samples as reported in Sims et al. [20]. (D) 9°50’N EPR samples (²²⁶Ra/²³⁰Th) versus their location relative to the AST. Decay curves show how (²²⁶Ra/²³⁰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 (²²⁶Ra/²³⁰Th) for the 9°50’N EPR axial samples [20]. (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 [227]. 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 30. (A) (²³⁸U/²³²Th) versus (²³⁰Th/²³²Th) for the 9°50’N off-axis samples compared to the 9-10°N EPR axial samples [20,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 (²³⁰Th/²³²Th) and (²³⁸U/²³²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. Inset provides regional context and shows the 9-10°N EPR axial samples [20,107] compared with the Siqueiros Transform samples [20,34]. Also shown are the near-axial basalts from 9°37’N EPR [187]. (B) ²³⁰Th/²³⁸U versus ²²⁶Ra/²³⁰Th for 9°50’N EPR off-axis samples compared with the 9°50’N EPR axial samples [20,107]. Error bars represent two sigma analytical uncertainties. Because of the significant difference in the half-lives of ²²⁶Ra and ²³⁰Th, (²²⁶Ra/²³⁰Th) decays back to equilibrium before (²³⁰Th/²³⁸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 (²²⁶Ra/²³⁰Th) in equilibrium, indicating that their ages are older (> 8 ka) than can be determined through ²³⁰Th-²²⁶Ra disequilibria. Inset provides regional context and shows (²³⁰Th/²³⁸U) versus (²²⁶Ra/²³⁰Th) for the 9°50’N EPR off-axis samples from this study compared with the 9-10°N EPR axial samples [20,107], three “near-axis” 9°30N EPR samples, including one E-type from close to the 9°37’N DEVAL [192], and the Siqueiros transform samples [34]. (C) 9°50’N EPR samples (²³⁰Th/²³⁸U) versus their distance from the center of the ridge axis, or Axial Summit Trough. Decay curves show how (²³⁰Th/²³⁸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 (²³⁰Th/²³⁸U) for the 9°50’N EPR axial samples as reported in Sims et al. [20]. (D) 9°50’N EPR samples (²²⁶Ra/²³⁰Th) versus their location relative to the AST. Decay curves show how (²²⁶Ra/²³⁰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 (²²⁶Ra/²³⁰Th) for the 9°50’N EPR axial samples [20]. (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 [227]. 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.

For most MORB lavas, ²²⁶Ra-²¹⁰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 ²¹⁰Pb (t_1/2 = 22 years) and the presence or absence of disequilibria [35,43,46,191,228]. However, as recently shown, for some lavas it is also possible to determine model eruption ages on a decadal scale using ²²⁶Ra-²¹⁰Pb disequilibria [81].

If the required source homogeneity criteria are met, then ²²⁶Ra-²¹⁰Pb model ages can be determined from the relationship:

$$T={\displaystyle \frac{1}{{\lambda }_{210}}}\ln {\displaystyle \frac{{\left({}^{210}Pb/{}^{226}Ra\right)}_m-1}{{\left({}^{210}Pb/{}^{226}Ra\right)}_o-1}}$$

(38)

where (²¹⁰Pb/²²⁶Ra)_m is the measured ratio, (²¹⁰Pb/²²⁶Ra)_o is the calculated initial ratio at the time of eruption of the lavas, and λ₂₁₀ is the decay constant of ²¹⁰Pb. Again, calculating (²¹⁰Pb/²²⁶Ra) model eruption ages requires an assessment of the initial extent of ²¹⁰Pb excess at the time of eruption of the samples in question.

Sims et al. [81] used (²¹⁰Pb/²²⁶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 (²¹⁰Pb/²²⁶Ra) to calculate model ages at this location was that the ABE vent site lavas had essentially constant (²²⁶Ra/²³⁰Th) and (²³⁰Th/²³⁸U) disequilibria (Figure 32). Given this constraint, Sims et al. [81] calculated the time of (²¹⁰Pb/²²⁶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.

4.4. Dating Young Volcanic Rocks with ²¹⁰Po-²¹⁰Pb

In subaerial lavas and shallow submarine volcanoes, ²¹⁰Po (t_1/2 = 138.4 days) is posited to degas entirely at magmatic temperatures, whereas its grandparent nuclide, ²¹⁰Pb (t_1/2 = 22.6 years), remains almost entirely in the melt (e.g. [229], and references therein). The intermediate radionuclide in the decay chain, ²¹⁰Bi, has a relatively short half-life (t_1/2 = 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). The ²¹⁰Po activity as a function of time is described by the Bateman equation:

$$\left({}^{210}\mathrm{P}\mathrm{o}\right)={\left({}^{210}\mathrm{P}\mathrm{o}\right)}_0{\mathrm{e}}^{-\mathsf{\lambda }\mathrm{t}}+\left({}^{210}\mathrm{P}\mathrm{b}\right)\left(1-{\mathrm{e}}^{-\mathsf{\lambda }\mathrm{t}}\right)$$

(39)

Assuming ${\left({}^{210}\mathrm{P}\mathrm{o}\right)}_0=0$, this simplifies to:

$$\left({}^{210}\mathrm{P}\mathrm{o}\right)=\left({}^{210}\mathrm{P}\mathrm{b}\right)\left(1-{\mathrm{e}}^{-\mathsf{\lambda }\mathrm{t}}\right)$$

(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 [230]. 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 [206,223]. Nevertheless, lavas that actively degas until eruption most commonly lose >90% of their Po [199,222,231,232,233], which means that calculated eruption ages assuming ²¹⁰Po = 0 may typically be too old by up to about 3 weeks.

Figure 31. Measured U-series disequilibria in SWIR MORB as a function of distance from the spreading axis. (a) Variable (²³⁰Th/²³⁸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 (²²⁶Ra/²³⁰Th) disequilibrium. Three of the five lavas with U-Th disequilibrium have excess ²²⁶Ra constraining their ages to < 8 ka. Alternatively, at 7.1 mm/yr, an age < 8 ka is equivalent to 55 m of movement since emplacement.

Figure 31. Measured U-series disequilibria in SWIR MORB as a function of distance from the spreading axis. (a) Variable (²³⁰Th/²³⁸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 (²²⁶Ra/²³⁰Th) disequilibrium. Three of the five lavas with U-Th disequilibrium have excess ²²⁶Ra constraining their ages to < 8 ka. Alternatively, at 7.1 mm/yr, an age < 8 ka is equivalent to 55 m of movement since emplacement.

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:

• In a study of lavas from Vailulu’u and Malumalu in Samoa, Sims et al. [35] demonstrate that the calculated ²¹⁰Po age, which provides an eruption date of 8 November, 2004 (+62 days/-67 days), falls within a window of bathymetric imaging showing new dome growth occurring sometime between March 1999 and April 2005 (Figure 33).
• In a study of the EPR, Tolstoy et al. [234] found a 2+ year buildup of seismic activity preceding a sudden drop in seismic activity on 22 January, 2006. While this event was interpreted as a singular eruption along the East Pacific Rise, the ²¹⁰Po dates in lava samples from this event show a wider range of eruption dates, from June 2005 to January 2006 [234]. The apparent discrepancy between the calculated ages and seismic observations is controversial; however, a simple interpretation that reconciles this discrepancy is that the ²¹⁰Po ages may reflect the buildup of seismic activity representing numerous smaller, distant eruptive events rather than precursors to a singular eruptive event.

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. [235]. (b) (²¹⁰Pb/²²⁶Ra) versus (²²⁶Ra/²³⁰Th) showing ABE vent field lavas and age trajectories calculated assuming (²¹⁰Ra/²²⁶Ra)₀ = 1.82 or 2.0 and (²²⁶Ra/²³⁰Th)₀ = 3.02 or 3.55 and then interpolating a model zero-age isochron between the two end-members. The majority of samples have (²¹⁰Pb/²²⁶Ra) > 1, indicating that they erupted within the past 100 years. RC-086 likely erupted more than 100 years ago; if this sample had (²²⁶Ra/²³⁰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.

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. [235]. (b) (²¹⁰Pb/²²⁶Ra) versus (²²⁶Ra/²³⁰Th) showing ABE vent field lavas and age trajectories calculated assuming (²¹⁰Ra/²²⁶Ra)₀ = 1.82 or 2.0 and (²²⁶Ra/²³⁰Th)₀ = 3.02 or 3.55 and then interpolating a model zero-age isochron between the two end-members. The majority of samples have (²¹⁰Pb/²²⁶Ra) > 1, indicating that they erupted within the past 100 years. RC-086 likely erupted more than 100 years ago; if this sample had (²²⁶Ra/²³⁰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.

5. Timescales of Magma Degassing

Magma recharge and degassing are directly linked to volcanic activity. Determining the timescales of these shallow-level processes is fundamental to our understanding of physical eruption dynamics and hazard assessment.

At many quiescent but potentially dangerous volcanoes, persistent degassing of volatile constituents dissolved in magma after melting is a common manifestation of volcanic activity. Placing constraints on magma residence times in shallow degassing reservoirs and/or feeding systems is critical to understanding magma dynamics beneath these volcanoes. Measurement of radioactive disequilibria between short-lived ²³⁸U-series isotopes in volcanic gases and lavas can provide constraints on the timescales of shallow-level magma dynamics (e.g. magma recharge and degassing). While there are a few studies that have measured ²²²Rn, a noble gas, or its daughter nuclides, ²¹⁰Po, ²¹⁰Pb, and ²¹⁰Bi in gas emissions from actively degassing open-conduit systems [229,230,232], most inferences about magma degassing timescales have come from measurements of ratios between ²¹⁰Po, ²¹⁰Pb, and ²²⁶Ra in lavas.

Figure 33. (A) Extrapolated best-fit ²¹⁰Po ingrowth curve assuming 100% Po degassing. Calculated maximum eruption age given with 2-sigma analytical error. (B and C) Bathymetric images showing dome growth occurring sometime between March 1999 and April 2005. Modified after Sims et al. [35].

Figure 33. (A) Extrapolated best-fit ²¹⁰Po ingrowth curve assuming 100% Po degassing. Calculated maximum eruption age given with 2-sigma analytical error. (B and C) Bathymetric images showing dome growth occurring sometime between March 1999 and April 2005. Modified after Sims et al. [35].

Measurements of (²¹⁰Pb)/(²²⁶Ra) in erupted lavas are ideal for constraining the processes and timescales of magma degassing because a) ²²²Rn, which is an intermediary nuclide between ²²⁶Ra and ²¹⁰Pb in the ²³⁸U-decay series (Figure 2 and Figure 34a), is a noble gas that is highly volatile; and b) the 22.3 year half-life of ²¹⁰Pb is appropriate for monitoring gas fluxes in magmas over a timeframe from several decades to years before eruption. Degassing processes generally have a negligible effect on the concentrations and activities of ²¹⁰Pb and ²²⁶Ra themselves, as Ra is not volatile and only about 1% of Pb typically degasses from magmas [229,230,232]. In contrast, Rn strongly partitions from magma into gas [199,236]. Hence, if magma degassing endures for years to decades, (²¹⁰Pb/²²⁶Ra) values are lowered as a function of the half-life of ²¹⁰Pb (Figure 34a and b) and the ²²²Rn degassing efficiency [221,225]. Similarly, if ²²²Rn concentrates in the gas phase of a magma that has stalled somewhere, then ²¹⁰Pb generated by decay of ²²²Rn can lead to a ²¹⁰Pb excess (Figure 34a and b) in the magma (e.g. [209,222,236]). However, because of the large differences in the half-lives of ²²²Rn (3.8 days) and ²¹⁰Pb, the erupted lavas must derive from a small part of the whole magma system, and the ratio of the magma that is degassing to the amount of magma that accumulates the ²²²Rn and erupts must be on the order of tens to thousands [222,225]. Lastly, if stored magma is undersaturated in volatiles or has maintained a steady pressure and temperature for a hundred years or more, (²¹⁰Pb)/(²²⁶Ra) may attain a state of radioactive equilibrium [89].

Like ²²²Rn, ²¹⁰Po, the immediate daughter of ²¹⁰Pb, is volatile at magmatic temperatures. With its short half-life of 138 days, ²¹⁰Po is useful for constraining degassing of magmas in the timeframe from months to weeks before eruption. Most lavas erupt with magmatic (²¹⁰Po/²¹⁰Pb) values near zero, indicating that they were degassing Po and other volatiles until eruption [199,222,231,232,233]. There are exceptions to this rule, however, including at Mount Erebus, which is discussed below. Prior work has attributed elevated (²¹⁰Po/²¹⁰Pb) values to (a) less efficient degassing resulting from a rapid rise of magma to the surface, (b) late-stage magmatic stagnation, or (c) infusion of Po from streaming gasses after magmas reach the near-surface environment (e.g., [206,238]).

Figure 34. (A) Schematic illustration showing how accumulation or degassing of ²²²Rn can lead to either (²¹⁰Pb) excesses relative to (²²⁶Ra) or (²¹⁰Pb) deficits relative to (²²⁶Ra). (B) Global compilation of (²¹⁰Pb/²²⁶Ra) versus (²²⁶Ra/²³⁰Th) for arc and oceanic samples. Data from Gauthier & Condomines [221], Turner et al. [209], Rubin et al. [191], Reagan et al. [222,223]; Berlo & Turner [224]. Note both ²¹⁰Pb excesses and deficits relative to ²²⁶Ra.

Figure 34. (A) Schematic illustration showing how accumulation or degassing of ²²²Rn can lead to either (²¹⁰Pb) excesses relative to (²²⁶Ra) or (²¹⁰Pb) deficits relative to (²²⁶Ra). (B) Global compilation of (²¹⁰Pb/²²⁶Ra) versus (²²⁶Ra/²³⁰Th) for arc and oceanic samples. Data from Gauthier & Condomines [221], Turner et al. [209], Rubin et al. [191], Reagan et al. [222,223]; Berlo & Turner [224]. Note both ²¹⁰Pb excesses and deficits relative to ²²⁶Ra.

In the following, we briefly discuss several examples using (²¹⁰Pb/²²⁶Ra) and, where available, (²¹⁰Po/²¹⁰Pb) to understand magma degassing processes in mid-ocean ridge (9°03N EPR and the Lau Basin ABE vent site) and ocean island/plume settings (Hawaii, Tristan de Cunha, Mount Erebus Volcano). These examples are chosen to illustrate situations where (²¹⁰Pb /²²⁶Ra) is greater than one, less than one, and at unity and clearly demonstrate that the relationships between crystallization, recharge, and volatile behavior are complex, but tractable, through the measurement of (²¹⁰Pb/²²⁶Ra) coupled with other petrological, geochemical, and isotopic constraints.

5.1. Magma Degassing in Ocean Ridge Settings

While there is a substantial data set for (²¹⁰Pb/²²⁶Ra) in young arc volcanic lavas (Figure 34b), there are relatively few studies using (²¹⁰Pb/²²⁶Ra) to understand magma degassing in MORB. Early measurements of (²¹⁰Pb/²²⁶Ra) disequilibria in MORB reported predominantly small ²¹⁰Pb deficits, with ²¹⁰Pb/²²⁶Ra as low as 0.86 [191,239], which were argued to be a result of melting (Figure 34b and Figure 35b), as Ra is significantly more incompatible in crystals than Pb during mantle melting [54].

However, small ²¹⁰Pb excesses (~10%) were measured in young MORB (Figure 34b and Figure 35b) from the Axial Seamount on the Juan de Fuca Ridge [191,240]. Radon is more incompatible than Pb during mantle melting making ²¹⁰Pb excesses more difficult to explain by simple mantle melting. However, because Pb is expected to behave more compatibly than Ra in both plagioclase and sulfide minerals, the ²¹⁰Pb excesses could have been generated by preferential melting of plagioclase or sulfide minerals in the mantle or crust (perhaps from melt-rock reactions). However, these minerals would need to be young and have pre-existing ²¹⁰Pb excesses. However, if this were the case, then one would expect correlations between geochemical indices of plagioclase or sulfide melting and (²¹⁰Pb/²²⁶Ra) disequilibria, which are not observed in MORB with ²¹⁰Pb excesses.

Figure 35. (A) (²²⁶Ra/²³⁰Th) versus (²¹⁰Pb/²²⁶Ra) for lavas from the overlapping spreading center at 9°03 N EPR. Note the significant excesses of ²¹⁰Pb relative to ²²⁶Ra for the basalts and basaltic andesites. (B) Axis-perpendicular schematic cross-section of the east limb of the 9103′N OSC, illustrating how Waters et al. [46] envision the formation of ²¹⁰Pb excesses in basalt beneath a dacite magma body and the eruption of an off-axis pillow ridge from an off-axis magma body. (C) The inset illustration is an enlarged view of the sub-axial magma bodies and corresponds to the dashed box above. The illustration is not to scale. Additionally, while the illustration suggests the presence of a dacitic magma body overlying the sub-axial magma sill, the location of dacitic magma is not explicitly known. Detailed images and information on the seismically determined distribution and geometry of subsurface melt can be found in Kent et al. [241] and Singh et al. [242].

Figure 35. (A) (²²⁶Ra/²³⁰Th) versus (²¹⁰Pb/²²⁶Ra) for lavas from the overlapping spreading center at 9°03 N EPR. Note the significant excesses of ²¹⁰Pb relative to ²²⁶Ra for the basalts and basaltic andesites. (B) Axis-perpendicular schematic cross-section of the east limb of the 9103′N OSC, illustrating how Waters et al. [46] envision the formation of ²¹⁰Pb excesses in basalt beneath a dacite magma body and the eruption of an off-axis pillow ridge from an off-axis magma body. (C) The inset illustration is an enlarged view of the sub-axial magma bodies and corresponds to the dashed box above. The illustration is not to scale. Additionally, while the illustration suggests the presence of a dacitic magma body overlying the sub-axial magma sill, the location of dacitic magma is not explicitly known. Detailed images and information on the seismically determined distribution and geometry of subsurface melt can be found in Kent et al. [241] and Singh et al. [242].

More recently, larger excesses of ²¹⁰Pb have been observed at two spreading-ridge locations, with (²¹⁰Pb/²²⁶Ra) as high as 1.69 at 9°03N EPR (Figure 35; [44]) and 1.82 at the Lau Basin ABE vent site (Figure 32; [81]). In both cases, the large ²¹⁰Pb excesses are posited to be caused by shallow magma chamber accumulation of ²²²Rn produced by deeper magmatic degassing. In general, buoyancy and viscosity relationships (i.e. Stokes' Law) suggest that gas bubbles will rise fast enough through basaltic magma (but not through more viscous, SiO₂-rich magmas) to allow differential motion of Rn-bearing gas and melt phases over timescales similar to the half-life of ²²²Rn. However, it is worth noting that, as suggested for Pinatubo [243], Rn-bearing gas may rise even more rapidly through pipe vesicles and fractured viscous magma.

In the case of 9°03’ N EPR (Figure 35), there is a wide range of erupted lava compositions that includes viscous, dacitic lavas and an imaged melt lens [241]. Waters et al. [46] envisioned a scenario in which a relatively large volume of basaltic magma residing within the lower to middle crust is continuously degassing. Radon-bearing, CO₂-rich gas bubbles rise and accumulate within a smaller, shallow reservoir of basalt and basaltic andesite beneath a larger volume of more viscous, dacite magma. ²²²Rn accumulates in the resulting bubble-rich basaltic froth and decays to ²¹⁰Pb. Density instability due to accumulation of this buoyant basaltic froth beneath dacitic magma may also facilitate magma mixing between basaltic and dacitic magma end-members, resulting in eruptions of mafic to intermediate lavas with large ²¹⁰Pb excesses.

In the case of the ABE lavas (Figure 32), which erupted off-axis from the main Lau Spreading center, Sims et al. [81] similarly posit that their ²¹⁰Pb excesses were created by a low-density, high-viscosity, rheological magma barrier which prevented volatile phase migration and caused a build-up of water and ²²²Rn in an underlying magma layer. This hypothesis is supported by two observations about the ABE lavas: 1) they have geochemical signatures suggesting they have been significantly influenced by subduction processes and therefore are relatively water-rich, and 2) they are relatively SiO₂-rich with high viscosities, which would cause a rheological barrier.

5.2. Magma Degassing in Ocean Island Plume Settings

5.2.1. Kilauea

Kilauea volcano on the island of Hawaii is one of the Earth’s most active volcanoes. During the period from 1983 to 2018, Kilauea erupted nearly continuously, predominantly from the Pu‘u‘ O‘o vent, with brief eruptions from other vents and fissure systems. Gas emissions at the summit caldera occurred over the entire period [244]. Girard et al. [163] analyzed (²¹⁰Pb/²²⁶Ra) ratios in lavas erupted over this period up through 2008. With the exception of lavas emitted when vents first opened, lavas from Pu‘u ‘O‘o, the nearby Kupahianaha vent, and Halema’uma’u within the summit caldera had (²¹⁰Pb)/(²²⁶Ra) = 0.80±0.11. These ²¹⁰Pb deficits were attributed to loss of Rn partitioned into CO₂-rich gasses from the point of volatile saturation at about 30 km depth to the surface over an average of 8 years. These data show that lavas erupted steadily at Kilauea rise to the surface after volatile saturation without significantly stalling within its magma chamber system. In contrast, more crystalline lavas that erupted immediately after the vent openings at Pu’u’ O’o and Napau, which is up-rift of Pu’u’ O’o, had near-equilibrium (²¹⁰Pb/²²⁶Ra) values, implying that these lavas had been stored in the rift system for several decades or more before being pushed out by newly intruding magmas.

Explosions in March 2008 opened a vent along the lower east wall of Halema’uma’u, ejecting basalts and altered lithic fragments. These explosions were followed by effusion of lavas that filled the crater. This marked the first eruption in recorded history when lavas simultaneously erupted from the summit and the East Rift Zone at Kilauea. Tephra samples with altered rock fragments and Pele’s hair had ²¹⁰Pb excesses over ²²⁶Ra resulting from ²²²Rn decay as well as minor Pb condensation and infusion from magmatic gasses (Figure 36). Lava samples erupted during this episode have (²¹⁰Pb/²²⁶Ra) values ranging from 0.8 to 1, suggesting they represented both stored and new lavas [163].

Several samples erupted between 2005 and 2008 were also serially analyzed for (²¹⁰Po). This serial analysis allows fitting of a radiogenic ingrowth curve and thus determination of both (²¹⁰Po) and (²¹⁰Pb) at the time of eruption. All lavas erupted from Pu‘u ‘O‘o have near-zero initial (²¹⁰Po)/(²¹⁰Pb) values (Figure 36). Several samples from Halema’uma’u had similarly low values, including one with an equilibrium (²¹⁰Pb)/(²²⁶Ra) value. Thus, newly arrived magmas as well as some of the shallowly stored magmas must have degassed until eruption. Pele’s hair and tephra samples erupted in 2008 from Halema’uma’u, in contrast, had significant excesses of (²¹⁰Po) over (²¹⁰Pb) with measured (²¹⁰Po)/(²¹⁰Pb) values as high as 6.5 soon after eruption. This is despite acid leaching of the Kilauea samples to remove ²¹⁰Po-bearing sublimates from sample surfaces. These ²¹⁰Po excesses were created by the infusion of Po from magmatic gasses into particles making up the fragmental deposits. A presence of sublimated Po on the surfaces of tephra particles was clearly shown by the high (²¹⁰Po) values measured in leachates, particularly from the Pele’s hair sample (Figure 36). This illustrates a rarely mentioned, but noteworthy, hazard associated with erupting, as well as passively degassing, volcanoes, namely the presence of short-lived radionuclides, including those that decay by energetic alpha emissions, on ash and smoke particles in volcanic plumes.

Figure 36. Measured (²¹⁰Po) values in decays per minute per gram (circles with 2σ error bars) and ingrowth curves for two samples erupted from Kilauea plotted against the number of days between the eruption and the date of analysis. Also plotted are (²²⁶Ra) values for these samples. Samples have been leached with 0.5 M HCl to remove polonium from surfaces. (A) Lava sample erupted in 2005. Partitioning of polonium into the gas phase led to an initial (²¹⁰Po) value near 0. The best fit ingrowth curve gives an initial (²¹⁰Pb) value of 0.204 dpm/g, which is below the measured (²²⁶Ra) value of 0.244 dpm/g. Deficits of ²¹⁰Pb with respect to ²²⁶Ra in Kilauea lavas are likely caused by persistent ²²²Rn loss during magma degassing before eruption [163]. For this sample, 6 years of persistent Rn-degassing would be needed. (B) Pele’s hair erupted in 2008. The leached sample had excess ²¹⁰Po, which decayed towards an initial (²¹⁰Pb) value above that of (²²⁶Ra) suggesting that the Pele’s hair fibers contained Pb absorbed from the gas phase. The leachate for this sample had (²¹⁰Po) = 952, illustrating the high (²¹⁰Po) concentrations in sublimated particles generated from cooling magmatic gasses [163].

Figure 36. Measured (²¹⁰Po) values in decays per minute per gram (circles with 2σ error bars) and ingrowth curves for two samples erupted from Kilauea plotted against the number of days between the eruption and the date of analysis. Also plotted are (²²⁶Ra) values for these samples. Samples have been leached with 0.5 M HCl to remove polonium from surfaces. (A) Lava sample erupted in 2005. Partitioning of polonium into the gas phase led to an initial (²¹⁰Po) value near 0. The best fit ingrowth curve gives an initial (²¹⁰Pb) value of 0.204 dpm/g, which is below the measured (²²⁶Ra) value of 0.244 dpm/g. Deficits of ²¹⁰Pb with respect to ²²⁶Ra in Kilauea lavas are likely caused by persistent ²²²Rn loss during magma degassing before eruption [163]. For this sample, 6 years of persistent Rn-degassing would be needed. (B) Pele’s hair erupted in 2008. The leached sample had excess ²¹⁰Po, which decayed towards an initial (²¹⁰Pb) value above that of (²²⁶Ra) suggesting that the Pele’s hair fibers contained Pb absorbed from the gas phase. The leachate for this sample had (²¹⁰Po) = 952, illustrating the high (²¹⁰Po) concentrations in sublimated particles generated from cooling magmatic gasses [163].

5.2.2. Tristan de Cunha

Tristan de Cunha island is a remote oceanic volcanic center in the South Atlantic Ocean at the west end of the Walvis Ridge, which has been attributed to mantle melting in a long-lived mantle plume. On July 19-20, 2004, a period of intense seismic activity was recorded by newly installed Comprehensive Test Ban Treaty Organization/International Monitoring System seismometers on Tristan da Cunha, indicating that there had been an eruption on a rift zone flanking the nearby Nightingale Island. A few days later, pumice began to wash up on local beaches. Samples of the pumice were collected by the British Geological Survey, and one was analyzed for (²¹⁰Pb/²²⁶Ra), as reported by Reagan et al. [223]. The pumice was phonolite, a highly differentiated lava produced by crystal fractionation from an alkaline basalt. ²¹⁰Pb was enriched over ²²⁶Ra by about 50% in the phonolite, despite accompanying evidence for degassing of volatile elements such as Po, and thus Rn, associated with the eruption. Instead of gas accumulation, the ²¹⁰Pb excess in this case was attributed to the extensive fractionation of feldspars and amphiboles required to produce phonolite from basalt, as both of these phases preferentially incorporate Ra over Pb into their crystal structures. Modeling of these and other data showed that the original basaltic magma was low in volume and had intruded into the rift zone no more than two centuries before eruption. The significant ²¹⁰Pb excess in the pumice showed that the crystal fractionation to generate the phonolite was likely ongoing until the time of eruption. Serial analysis of ²¹⁰Po in the same sample (Figure 37) showed that it had an initial (²¹⁰Po)/(²¹⁰Pb) value of 0.15 upon eruption, illustrating that the pumice in question was indeed ejected during the July 2004 eruption [223].

Figure 37. Measured (²¹⁰Po) and ingrowth curves for an acid-leached sample of phonolite pumice collected from a beach on Tristan da Cunha in September 2004. Also plotted are (²²⁶Ra) values for this sample. The illustrated ²¹⁰Po deficit with respect to ²¹⁰Pb shows that this pumice erupted within a few months of the date of collection. The most likely origin of the pumice was a subsea eruption associated with an earthquake swarm offshore of Nightingale Island on July 29 and 30, 2004. The (²¹⁰Pb) excess over (²²⁶Ra) calculated from the (²¹⁰Po) ingrowth curve was attributed to rapid crystal fractionation of Ra-bearing minerals from the magma over a period of several decades to a century before the eruption [223].

Figure 37. Measured (²¹⁰Po) and ingrowth curves for an acid-leached sample of phonolite pumice collected from a beach on Tristan da Cunha in September 2004. Also plotted are (²²⁶Ra) values for this sample. The illustrated ²¹⁰Po deficit with respect to ²¹⁰Pb shows that this pumice erupted within a few months of the date of collection. The most likely origin of the pumice was a subsea eruption associated with an earthquake swarm offshore of Nightingale Island on July 29 and 30, 2004. The (²¹⁰Pb) excess over (²²⁶Ra) calculated from the (²¹⁰Po) ingrowth curve was attributed to rapid crystal fractionation of Ra-bearing minerals from the magma over a period of several decades to a century before the eruption [223].

5.2.3. Mount Erebus Volcano

Mount Erebus Volcano has a HIMU radiogenic Sr, Nd, Hf, and Pb isotopic signature [245] and a thermal structure suggestive of deep mantle upwelling [246,247]. Mount Erebus also has the most persistent, continuously degassing lava lake on Earth [89,90,245,248]. The Erebus lava lake degasses both by quiescent emissions and strombolian eruptions with distinct gas compositions [249]. These originate from different depths within the magma conduit, with the quiescent emissions coming from shallow passive degassing, and deeper eruptions deriving from deep within the system ([90], and references therein).

Measurements of U and Th decay series nuclides in young Erebus lava bombs (collected at the time of eruption) show (²¹⁰Pb/²²⁶Ra) of unity and are thus in radioactive equilibrium [89]. This observation, that ²¹⁰Pb is in equilibrium with ²²⁶Ra, is rather surprising as one would expect that removal of ²²²Rn by magma degassing would have a significant impact on the (²¹⁰Pb/²²⁶Ra) of the magma. To explain the measured (²¹⁰Pb/²²⁶Ra) of unity in the young, known-age Erebus phonolite bombs, Sims et al. [89] posited that magma degassing for the phonolitic magma occurred deeper and hundreds of years before eruption and therefore (²¹⁰Pb) and (²²⁶Ra) had reached a state of radioactive equilibrium. In this scenario, the anorthoclase-rich phonolite circulating in the upper part of the magmatic system is made of old, degassed magma that does not communicate with gas emissions derived from deeper parts of the magma chamber. These ongoing emissions are rising to the surface and causing the periodic strombolian eruptions (see e.g., Sims et al. [90], and references therein for an overview of degassing processes on Erebus volcano).

This scenario is consistent with a) the observation that the melt inclusions from anorthoclase megacrysts are highly degassed compared to melt inclusions in olivine from parental basalts from around the flanks of the volcano [250,251,252]; b) calculated magma residence times, using measurements of all relevant ²³⁸U, ²³⁵U, and ²³²Th decay nuclides from phonolite glass and anorthoclase megacrysts in known-age lava bombs, are on the order of a thousand years (Figure 29); and c) for the 2005 lava bomb, whose eruption date (16 December, 2005) is known explicitly, serial analysis showed (²¹⁰Po)/(²¹⁰Pb) = 0.28, which is significantly higher than observed in most zero-aged lavas, indicating that ²¹⁰Po was not completely degassed from the lava at the time of eruption. This (²¹⁰Po)/(²¹⁰Pb) value was modeled to show that the shallow magma/lava lake system, on average, degasses about 1% of its Po per day [89].

6. Brief Overview of Techniques for Measuring U-Series Nuclides in Oceanic Volcanic Rocks

The U and Th decay series isotope systematics (see section 2) make the measurement of U and Th decay series nuclides quite challenging for two reasons:

• Because the equilibrium abundances of the U and Th decay series nuclides are proportional to their respective half-lives (a=λn, where a=activity, λ= ln(2)/half-life, and n=number of atoms), the short-lived progeny of ²³²Th, ²³⁸U, and ²³⁵U occur in very low concentrations. For example, in most oceanic basalts, ²³²Th and ²³⁸U are generally found in weight concentrations of ng/g to ug/g, ²³⁴U and ²³⁰Th in concentrations of pg/g, ²²⁶Ra in concentrations of fg/g, and ²²⁸Ra and ²¹⁰Pb in concentrations of ag/g.
• Isotope ratios involving long-lived and short-lived nuclides are often quite extreme as the equilibrium atomic ratios of the various elements are set by the ratio of their decay constants; as such the measured ratios are typically on the order of 10^-5 for ²³⁴U/²³⁸U, 10^-6 for ²³⁰Th/²³²Th, 10^-3 for ²²⁶Ra/²²⁸Ra, and 10^-8 for ²¹⁰Pb/²⁰⁸Pb.

Early measurements of U decay series nuclides relied exclusively on alpha- and gamma-spectrometric methods that counted individual decay events [186,199,231,233,253]. These counting methods provided seminal and cutting-edge U and Th-decay series data and are still used today for short-lived nuclides such as ²¹⁰Pb, ²¹⁰Po, ²²⁸Th, and ²²⁷Ac [37,45,81,88,89,163,206,254,255]. However, for the longer-lived U and Th decay series isotopes ²³²Th, ²³⁸U, ²³⁰Th, and ²²⁶Ra, these decay counting techniques have very large uncertainties due to their low count rates. Thus, mass spectrometric methods that count atoms rather than decay events, such as TIMS, SIMS, and MC-ICPMS, significantly improve sensitivity and precision for the measurement of the longer-lived isotopes such as ²³⁸U-²³⁴U-²³⁰Th-²²⁶Ra and ²³⁵U-²³¹Pa, while simultaneously decreasing analysis times significantly [35,167,210,256,257,258,259,260,261,262,263,264]. Readers are referred to the above papers for details on chemical separations and purification methods and counting and mass spectrometric techniques. For a one-to-one comparison of samples measured by both alpha and gamma versus mass spectrometry, see Sims et al. [262].

Finally, quality assurance needs careful attention, particularly when calculating U-series disequilibria, as these ages involve multiple measurements and therefore propagation of many and varied sources of analytical uncertainty. This includes determining the uncertainties of each individual factor leading to the analytical measurement. Such uncertainty propagation can become complicated when there is covariance between measurement errors. Several papers explicitly address the proper handling of errors associated with U and Th decay series chronological measurements (e.g., [265], and references therein). It is important to note that reference materials, run as unknowns, are essential for establishing both the precision and accuracy of the measurements. One crucial quality of U and Th decay series measurements is that the equilibrium decay principles establish the relative abundances of the parent-daughter pairs. Determining both the precision and accuracy of these abundances requires use of unaltered quality assurance materials (e.g., well characterized international rock standards) that are old relative to the half-lives of the isotope system of interest (e.g. [262,263,265,266]).

7. Summary

This paper demonstrates the potential of U and Th decay series for unraveling the temporal complexities of oceanic magmatism. While there admittedly are remaining uncertainties limiting the results that are obtainable from U and Th decay series measurements, there are several robust results.

• Defining the timescales of magma genesis during decompression mantle melting. The formation and transport of basaltic magma in the mantle during adiabatic decompression can be characterized in terms of the velocity of the upwelling solid mantle, the rate of melting, the melt velocity, melt transport processes (time scale, focused and diffuse, closed or open system), and the delay time, or storage time, of magma before eruption. While the velocity of upwelling mantle and the corresponding melting rate can be estimated from models coupled with thermodynamic data and geophysical observations (e.g., [23,24,150,267]), the melt velocity and magma storage times are difficult to estimate from first principles and require measurement. Fortunately, because the daughter isotopes ²³⁰Th and ²²⁶Ra of the ²³⁸U decay series and the daughter isotope ²³¹Pa of the ²³⁵U decay series have half-lives that bracket the timescales over which melting and melt extraction are thought to occur, measurements of (²³⁰Th/²³⁸U), (²²⁶Ra/²³⁰Th), and (²³¹Pa/²³⁵U) ratios in young basaltic lavas provide unique information on the timescales and modes of melt generation and transport.
• Determining the eruption ages of oceanic Quaternary volcanism. Without lava ages it is not possible to relate volcanic activity to tectonics processes, establish the timing and duration of changes in lava compositions, and link hydrothermal and biological activities to fluctuations in submarine and subaerial oceanic volcanism. Thus, age constraints are not a detail, they are necessary. Without them, our understanding of oceanic systems will remain limited. Dating young submarine lavas has proven to be notoriously difficult – long-lived radiogenic isotope systems such as Rb-Sr, Sm-Nd, Lu-Hf, and U-Th-Pb have half-lives that are too long and hence lack the temporal resolution for precise ages; cosmogenic age dating on the ocean floor is not possible; potassium-argon dating techniques suffer from a myriad of problems when applied to young oceanic basalts; and sedimentation rates are known to vary widely over short spatial scales. Thus, the visual appearance of lavas provides only a subjective, semi-quantitative estimate of age [268]. Fortunately, ²³⁸U and ²³²Th decay series disequilibria can provide robust constraints on the eruption ages of young oceanic basalts. The application of the ²³⁸U, ²³⁵U, and ²³²Th for dating eruption ages of volcanic rocks uses three different approaches: 1) age limits, whereby the absence or 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 estimate the eruption age of an individual lava or a sequence of lavas.
• Establishing the timescales of magma recharge and degassing. Processes operating within magma reservoirs affect both the physical state and the chemical composition of magmas within them. Measurements of (²¹⁰Po), (²¹⁰Pb), and (²¹⁰Bi) in gas emissions from actively degassing open-conduit systems and (²¹⁰Pb)/(²²⁶Ra) and (²²⁸Ra/²²⁸Th) in young known-age erupted lavas provide otherwise unobtainable constraints on the timescales of magma recharge and degassing. These timescales are directly linked to a volcano’s activity and are thus fundamental to our understanding of physical eruption dynamics and hazard assessment.