Investigating Influences on the Pb Pseudo-Isochron Using Three-Dimensional Mantle Convection Models With a Continental Reservoir

For mid-ocean ridge basalts and ocean island basalts, measurements of Pb isotope ratios show broad linear correlations with a certain degree of scatter. In 207 Pb/ 204 Pb— 206 Pb/ 204 Pb space, the best fit line defines a pseudo-isochron age ( τ Pb ) of ∼1.9 Gyr. Previous modeling suggests a relative change in the behaviors of U and Pb between 2.25 and 2.5 Ga, resulting in net recycling of HIMU (


Introduction
In the field of mantle studies, the coupled U-Th-Pb isotope systems provide valuable information on crustal recycling and mantle stirring.An integral record of U-Th-Pb fractionation during geodynamic processes is provided by measurements of the radiogenic Pb isotope ratios (Gast et al., 1964;Ito et al., 1987;Tatsumoto, 1966;White, 1985).The isotopes 238 U, 235 U, and 232 Th decay into different isotopes of Pb ( 206 Pb, 207 Pb, and 208 Pb respectively) at different rates dictated by the parent's half-life, providing powerful means with which to investigate global scale geodynamic and geochemical processes.
Variations in the Pb isotope ratios measured in mantle derived basalts (Figure 1) are driven by long term differences of U/Pb and Th/Pb in their sources.While melting at mid-ocean ridges is likely to play a role in separating U and Th from Pb, the effectiveness of this process depends on the ratio of partition coefficient (D = [element] • We present numerical geodynamic models with new parameterizations for U recycling and preferential Pb removal from subducted oceanic crust • A combination of these processes provides a good fit to both Pb pseudo-isochron and observed scatter of Pb isotope ratios measured in oceanic basalts • Our models do not require long term accumulation of subducted oceanic crust to fit geochemical constraints

Supporting Information:
Supporting Information may be found in the online version of this article. 10.1029/2021GC010309 2 of 24 et al., 1994;Peucker-Ehrenbrink et al., 1994).The former occurs due to the fluid mobility of U in the U +6 oxidation state, allowing oxidized U to be transported from the continents into the oceans where it is subsequently incorporated into altered oceanic crust via hydrothermal addition (Collerson & Kamber, 1999;Michard & Albarede, 1985).Lead and Th do not exhibit the same behavior, so under an oxidising atmosphere oceanic crust is expected to become enriched in U relative to Pb and Th (Galer & O'Nions, 1985).During subduction Pb, U and Th are removed from subducted oceanic crust to different extents at different depths.Pb is preferentially removed at shallow depths, causing it to be sequestered to the lower continental crust or lithospheric mantle (Kellogg et al., 2007;Kramers & Tolstikhin, 1997), or incorporated into arc magmas (Kelley et al., 2005;Miller et al., 1994).As the subducting slab reaches greater depths, U is preferentially removed over Pb but much of this may be re-incorporated into the upper mantle (Elliott et al., 1999;Kelley et al., 2005).
The net result of these fractionation processes is that subducted oceanic crust can acquire a high U/Pb (and Th/Pb), which with time will give rise to strongly radiogenic Pb isotope signatures.The ratios 206 Pb/ 204 Pb and 207 Pb/ 204 Pb contain information on the time integrated 238 U/ 204 Pb (μ) and 235 U/ 204 Pb while the combination of 208 Pb/ 204 Pb and 206 Pb/ 204 Pb informs on the time-integrated 232 Th/ 238 U (κ m ).Oceanic basalts with the most radiogenic Pb isotope ratios, such as lavas from Tubuai, have been dubbed HIMU (Zindler & Hart, 1986) as they reflect mantle sources with long-term elevated μ.
Lead isotope ratios measured in MORBs and OIBs show a roughly linear trend in both 207 Pb/ 204 Pb-206 Pb/ 204 Pb and 208 Pb/ 204 Pb-206 Pb/ 204 Pb space (Figure 1) (Gast et al., 1964;Tatsumoto, 1966).The slope in 206 Pb/ 204 Pb-207 Pb/ 204 Pb space (Figure 1a) potentially has age significance, dating the average time since the last U/Pb fractionation.The linear regression through this data field can be thought of as a pseudo-isochron; "pseudo" as rather than dating a single melting event this represents the sum of multiple melting events during plate creation, subduction and stirring back into the mantle.(1) The left-hand side is the gradient of the regression line in 207 Pb/ 204 Pb-206 Pb/ 204 Pb space, and λ 235 & λ 238 are the decay constants for 235 U and 238 U respectively.We use a 235 U/ 238 U ratio of 137.88, a standardised value (Steiger & Jäger, 1977) used widely in previous literature.Despite the different source regions of MORBs and OIBs, their Pb isotope composition plots on a similar array (Figure 1), with their combined τ Pb being 1.96 Ga (Figure 1a).
Another constraint on the mantle Pb isotope composition is the scatter observed in the Pb isotope arrays 207 Pb / 204 Pb-206 Pb/ 204 Pb and 208 Pb/ 204 Pb-206 Pb/ 204 Pb.This relies on the assumption that the scatter is the product of geologic processes and not just an artifact of the analytical technique.Studies quantifying the difference in Pb isotope ratios measured using different analysis techniques (Weis et al. (2006)) find marginal (±1%-2%) differences, however, the error associated with the measurements is smaller than the magnitude of the scatter.This indicates that at least some of the scatter is indeed real.We also find only a limited difference between the scatter in Pb isotopic ratios collected via different techniques from a compilation of published measurements (Figure S1 in Supporting Information S1).It is therefore reasonable to quantify the scatter observed in MORBs and OIBs in order to compare against modeled data.We quantify the scatter (d) as the average orthogonal distance away from the pseudo-isochron line, and for model data this measure is weighted by the amount of 204Pb (Figure 1).For the data in Figure 1 In order to better understand what can be inferred from the distribution of Pb isotope ratios in mantle derived basalts, statistical and numerical models of mantle processes have been employed.If observations such as τ Pb and the spread of Pb isotope ratios can be reproduced in models which reasonably parameterize terrestrial geochemical processes then they may help us to better understand the timescales of mantle mixing and key events in Earth's history which have shaped the mantle's geochemistry.
Early box modeling approaches to the problem yielded unrealistically high τ Pb (Allégre et al., 1980;Armstrong & Hein, 1973), a finding that was replicated in numerical mantle convection simulations (Christensen & Hofmann, 1994).This was attributed to the rapid ingrowth of 207 Pb early in Earth's history, thus causing early 10.1029/2021GC010309 3 of 24 differentiated material in the model to develop highly radiogenic 207 Pb/ 204 Pb, which strongly influences τ Pb .To overcome this Christensen and Hofmann (1994) chose to begin their simulations at 3.6 Ga with a uniform distribution of trace elements, under the assumption that high temperatures of the early Earth would efficiently homogenize any differentiation that had occurred.
Later numerical investigations overcame the need to initialize models from 3.6 Ga by invoking a change in the relative behavior of Pb and U at some point in Earth's history.Such a change may have been brought on due to the onset of subduction or oxidation of Earth's atmosphere and ocean.Xie and Tackley (2004) found that their simulations provided a good match to the τ Pb of oceanic basalts if the production of HIMU in recycled crust was prevented until 2.5 Ga.Around this time the atmosphere is thought to have become sufficiently rich in oxygen to allow the recycling of U from the continents into the mantle (Lyons et al., 2014).In Xie and Tackley (2004), fractionation of U and Th from Pb followed the method Christensen and Hofmann (1994) in which both magmatic and non-magmatic fractionation were incorporated into a single process, namely an unrealistically large difference between the partition coefficients D U and D Pb .A weakness of this parameterization is that it neglects the role of continental crust, which acts as a geochemical reservoir separate from the mantle.
The models presented by Brandenburg et al. (2008), guided by the results of statistical box models (Kellogg et al., 2002(Kellogg et al., , 2007)), parameterized the net effect of U recycling from the continental crust and preferential removal of Pb from subducted oceanic crust into a single process.A proportion of each element is removed from the melt to the continental crust, controlled by their corresponding extraction coefficients.In their preferred case a relative change in behavior between U, Th and Pb, likened to a change in subduction conditions, is applied at 2.25 Gyr by changing the relative removal rates of elements via extraction coefficients.With this two-stage process (pre and post change in subduction conditions) they manage to produce a slope in 207 Pb/ 204 Pb-206 Pb/ 204 Pb space similar to oceanic basalts, as well as the full range of observed Pb isotope ratios, highlighting the importance of parameterizing processes such as depth dependent extraction rates for different isotopes.
While previous models have included parameterizations that effectively model relative changes in non-magmatic fractionation processes affecting the U-Th-Pb system (Brandenburg & van Keken, 2007a;Brandenburg et al., 2008; 10.1029/2021GC010309 4 of 24 Xie & Tackley, 2004), none include a mechanism for explicitly modeling either U recycling from the continental crust or preferential removal of Pb from subducted oceanic crust.This makes it difficult to disentangle the relative effects of each process.Doing so is important as U recycling is a redox controlled process affected by the composition of Earth's atmosphere while preferential removal of Pb is related to subduction.It is unlikely that the timing of these disconnected processes is well correlated so it is sensible to separate them out if we are to better understand how tectonics influences the distribution of Pb in the mantle.Additionally, previous numer ical modeling has exclusively been conducted in 2D geometry whereas 3D geometry is preferable in order to better represent the stirring and processing efficiency of the mantle, which is key in eliminating old heterogeneity.
Here we present 3D mantle convection simulations that include separate mechanisms for recycling U from a continental reservoir into the mantle, the preferential removal of Pb, and melting fractionation.In doing so we will be able to unpick the way in which different fractionation processes affect the distribution of Pb isotopes in the mantle.We first present simulations to determine the effect of pure melt fractionation on the modeled Pb isotope ratios and a simulation with similar setup to that of a preferred case of Xie and Tackley (2004) to replicate their results.Subsequent cases feature a new parameterization for the recycling of U from the continental reservoir as well as a separate process for preferentially removing Pb from oceanic crust.In our analysis of the models we will use both the 207 Pb and 208 Pb constraints in parallel.By including communication between the mantle and continental reservoirs, we can assess how successfully we can reproduce the full characteristics of mantle Pb isotope systematics, which have long been puzzling and expressed in terms of two paradoxes.The first Pb paradox (or the future paradox) concerns the almost ubiquitous observation that Pb isotope compositions of MORBs and marine sediments plot to more radiogenic values than the array defined by meteorites in 207 Pb/ 204 Pb-206 Pb/ 204 Pb space (Allégre, 1982;Sinha & Tilton, 1973).If Earth's bulk composition is indeed similar to that of such meteorites, there must exist an un-radiogenic complement 'hidden' somewhere that we cannot observe it.The second paradox arises due to the mismatch between the upper mantle Th/U ratio (∼2.6-2.7)inferred from the composition of MORBs and the time-integrated Th/U ratio given by 208 Pb/ 204 Pb and 206 Pb/ 204 Pb ratios (∼3.2) (Galer & O'Nions, 1985;O'Nions & McKenzie, 1993).

∇ ⋅ 𝐮𝐮 =0
(2) Variables and parameters in these equations are the fluid velocity u, viscosity η, pressure P, thermal expansivity α, density ρ, acceleration due to gravity g, average mantle temperature T av , temperature Thermal conductivity Note.Note that the reference viscosity is equal to the upper mantle viscosity.1) and by the isothermal CMB.An initial thermal condition is generated from a random temperature field which is run forward for 5 Gyr.Each simulation is run from 3.6 Ga to present day as in previous studies (Brandenburg & van Keken, 2007b;Christensen & Hofmann, 1994;van Heck et al., 2016;Xie & Tackley, 2004).This avoids modeling early Earth conditions, during which it is likely that the mantle would have had a significantly lower viscosity due to extremely high temperatures.

Particles
Trace element abundances, bulk composition (C), mass and melting age (time since last melted) are stored on active tracer particles (van Heck et al., 2016) which are advected through the grid in the mantle flow.
C represents the fusible component of mantle material and can have a value from 0.0 to 1.0, with 0.0 representing completely depleted material (harzburgitic) and 1.0 representing completely enriched material (basaltic).Half of the particles are initialized with a composition of C = 0.25, 3/8 with C = 0.0, and 1/8 with C = 1.0 with each composition being evenly distributed throughout the mantle giving an average mantle composition of C = 0.25.The composition and depth of particles affects the local density, with completely basaltic material being 4% denser than completely depleted material in the upper mantle and 3% denser in the lower mantle (Ono et al., 2001).Across the CMB the buoyancy ratio (B) = 0.33, calculated as the ratio of chemical to thermal density contributions using Where Δρ b is the density difference between material with a basaltic composition and material with an average mantle composition (C = 0.25) in the lower mantle (101 kg m −3 ), α is the thermal expansion coefficient, ρ 0 is the mantle reference density and ΔT is the temperature difference across the mantle.The "basalt barrier" in the mantle transition zone is caused by the delayed transition to dense, lower-mantle mineral phases in cool subducted oceanic crust from 660 to 750 km (Irifune & Ringwood, 1993).We model this by making basalt 5% more buoyant in the transition zone compared to harzburgitic compositions (Davies, 2008).We include a simplified parameterization of the olivine system phase transitions at 410 and 660 km depth which do not vary with bulk composition (Table 2).We note that this implementation is not entirely self consistent with our bulk composition parameterization but does offer the ability to capture some of the known behaviors associated with these phase changes (Price et al., 2019;Wolstencroft & Davies, 2011).

Initializing Chemistry
Trace elements are initially distributed in equal ratios across all particles, with abundances calculated from each particle's mass.This implies that any fractionation that had taken place before this time has been efficiently mixed (Christensen & Hofmann, 1994).Initial concentrations can be found in Table 3.The present day value for 238 U is calculated from an estimate of its current concentration in bulk silicate Earth (BSE).From this, the concentrations of 235 U and 232 Th are estimated from their respective present day molar ratios to 238 U.The 204 Pb isotope is stable and so its present day global abundance, calculated by its inferred modern molar ratio to 238 U, is equal to its starting abundance.Initial values for radiogenic lead isotopes 206 Pb, 207 Pb, and 208 Pb are estimated from initial abundance on accretion taken from the Note.Density difference is relative to PREM (Dziewonski & Anderson, 1981).Canyon Diablo meteorite reference (Tatsumoto et al., 1973), plus ingrowth from decay of their parent isotopes: i Pbs= 204 Pb × i∕204 PbCD +PD (7) Where i = 206/207/208.i Pb s is the abundance of a radiogenic lead isotope at the start of the calculation (t s = 3.6 Ga), i/204 Pb CD is its ratio to 204 Pb at the formation of Earth (t 0 = 4.56 Ga), and P D difference in abundance of the parent isotope between t 0 and t s .These abundances are calculated from the decay equation: Where j = 235/238/232 and j P t represents the abundance of parental isotopes 235 U, 238 U or 232 Th at time t, j P pd is the present day parental abundance, Δt is the difference in time between t and present day, and λ j is the corresponding decay constant for the parent isotope.

Melting
Our melting method follows that of van Heck et al. (2016).The solidus, dependent on depth z and composition C is defined by Where T meltsurf (1,200 K) is the melting temperature of C = 1 material at the surface, T meltslope (2.5 K km −1 ) is the slope of the solidus, and T meltcomp (500 K) is the temperature difference between the solidi of C = 0 and C = 1 material (varies linearly).The temperature of each particle is linearly interpolated from the grid and if its solidus is exceed then melting will occur.The new composition (C n ) of a melting particle is reduced so that Under the assumption that the mass of the volume represented by the particle is in thermal equilibrium with the new solidus (T).Note that for a melting particle the composition value will always decrease, that is, the particle becomes more depleted in its basaltic component.This depletion may occur up to C = 0, at which point the particle is so refractory in composition that it can no longer melt.While this melting process does conserve energy, it neglects the effects of latent heat and thermal advection due to melt movement (Nakagawa & Tackley, 2012;Xie & Tackley, 2004).The material depleted from a particle forms a "melt package" to move chemical information.The degree of melting F is then given by Where C 0 is bulk composition of the particle prior to melting.Trace elements are also lost during melting events, with the amount being removed being given by Where i refers to the isotope being moved, A m−i is the number of atoms of each isotope that is removed from the melting particle, F is the degree of melting (Equation 11), A s−i the number of atoms of each isotope on the particle before melting occurs and D i is the partition coefficient of each isotope (listed in Table 4).The melting process therefore fractionates elements which have different partition coefficients.Each melt package instantaneously migrates, transferring its basaltic component and trace elements to particles which are not yet fully enriched in the surface cell vertically above the cell of the melting particle (van Heck et al., 2016).If all particles in the surface cell are completely enriched then the particles in the next radial cell below will be filled.This process leads to an enriched, basaltic crust forming at the surface, underlain by a depleted, harzburgitic layer.We note that our melting process does not distinguish between intrusive and extrusive volcanism, which would be important for Note.Note that D U = 0.007, D Th = 0.007 for all cases.Dashes indicate that a particular process is absent in this model and times indicate the time a parameter or process begins.Cases with U-recycling style U a and U b reserve 1/2 and 1/3 respectively of all the continental U at the indicated time for recycling into the mantle.Cases with Pb-removal styles Pb a and Pb b set the fraction of Pb removed from melt to the continental reservoir to be f pbc = 1/4 and f pbc = 1/10 respectively.Note that for cases XIE007 and XIE010, we show D Pb before 2.4 Ga, followed by D Pb after 2.4 Ga.

Table 4
Parameter Settings for Cases

Continental Crust
As well as fractionation via differences in partitions coefficients, we investigate the effects of recycling U from the continental crust to the mantle and preferential removal of Pb from oceanic crust.To model these processes we implement a "continental reservoir" which sits external to the mantle and acts as a vessel to store number of moles of trace elements which would be isolated from the mantle.Unlike in Brandenburg et al. (2008), where their continental reservoir is continuously populated by extraction of trace elements after melting, our continental reservoir is initialized with 1/3 of the global budget of U, Th, and Pb at the start of the calculation (3.6 Ga).This implies rapid early extraction of the continental crust prior to the start of the simulation (Armstrong, 1968).
Initially the trace elements are completely un-fractionated (BSE) ratios but the processes of U recycling and Pb removal from melts cause the U/Pb of the continental crust to decrease with time.This differs from Brandenburg et al. (2008) where a low U/Pb for the continental crust arises due to differences in the extraction coefficients for different elements.
For simulations which include recycling of U (Table 1), a fraction of the 235 U and 238 U in the continental reservoir at either 2.4 Ga or 0.6 Ga is reserved to be recycled into the mantle.The purpose of this is to maintain some control over the total amount of U that is recycled over the course of the simulation.Note that this reserved portion of the continental reservoir still experiences radioactive decay.Recycling begins either at 2.4 or 0.6 Ga, to coincide with the great oxygenation event (GOE) or widespread ocean oxidation (Lyons et al., 2014) respectively.The amount of a uranium isotope to be recycled at each time step (U rts ) is given by Where U cr is the remaining amount of a U isotope in the continental reservoir that has been set aside for recycling, Δt is the time step, t tot is the total time the calculation will run for, and t is the current time through the calculation.
U cr is then updated The U to be recycled in a given time step is evenly distributed across all the particles in the top most (surface) radial layer of the model, approximating the hydrothermal addition of U to oceanic crust (Collerson & Kamber, 1999;Michard & Albarede, 1985).
In modeling the preferential removal of Pb from subducted slabs we assume that all melt will become oceanic crust and eventually be subducted.Under this assumption, removing Pb from the melt has the same effect as removing it from subducted oceanic crust and means we do not have to identify downwelling regions.Thus a set fraction of all Pb from the produced melt is added to the continental reservoir in each step: Where Pb cr is the abundance of a given Pb isotope in the continental reservoir, Pb m is the abundance of the Pb isotope in a melt package and f pbc is a fixed fraction representing the amount of the Pb isotope that is removed from the melt to the continental reservoir.This process therefore removes Pb from the convecting mantle (Kelley et al., 2005;Kellogg et al., 2007;Kramers & Tolstikhin, 1997).We begin preferential Pb removal at 3.0 Ga, which is a commonly estimated time for the onset of plate tectonics (Laurent et al., 2014;Shirey & Richardson, 2011;Tang et al., 2016).

Parameter Space
For comparability between simulations we use the same partition coefficients for U (D U = 0.007) and Th (D Th = 0.007), as those used by (Xie & Tackley, 2004).The first set of models revisits previous modeling approaches (Christensen & Hofmann, 1994;Xie & Tackley, 2004), using just differences in D Pb to induce 10.1029/2021GC010309 8 of 24 fractionation.We then introduce model processes for U recycling and Pb removal from subducted oceanic crust (see Table 4).
In case MELT, we set D Pb = 0.010.There is a great deal of uncertainty in D Pb measurements, but using this value allows for minor fractionation to occur between Pb and U and Th.Hi-DPb uses the same exaggerated D Pb (0.025) relative to D U and D Th used in previous studies (Xie & Tackley, 2004) in order to draw some level of comparison.Case XIE007 has a similar setup to a simulation in Xie and Tackley (2004), where initial D Pb = D U = 0.007 before 2.4 Ga, implicitly assuming any U-Pb fractionation occurring before this time to have been homogenised.After 2.4 Ga D Pb = 0.025 as in Hi-DPb.XIE010 also includes a change in D Pb = 0.025 at 2.4 Ga but prior to this D Pb = 0.010, allowing for some fractionation to take place.
Subsequent cases have D Pb = 0.010 for the entirety of the simulation.Uranium recycling from the continental reservoir is introduced in CONTU, to coincide with the GOE (Lyons et al., 2014), with 1/2 of the U available in continental reservoir at 2.4 Ga being reserved for recycling.This is paired with Pb removal in UPb2.4,where 25% of the Pb in the melt is extracted to the continental reservoir from 3.0 Ga.The onset of U recycling is delayed in case UPb0.6 to 0.6 Ga, until when the oceans are likely to have been fully oxygenated (Lyons et al., 2014).In case LO-UPb we reduce the amount of U recycling by reserving 1/3 of the U available in the continental reservoir at 2.4 Ga for recycling.The rate of Pb removal is also reduced by extracting 10% of the Pb from the melt to the continental reservoir.The range of Pb removal rates investigated reflects that there is uncertainty over the efficiency of Pb removal relative to U removal from subducted oceanic crust (Kelley et al., 2005).
Case SCALE investigates how an arguably more realistic rate of mantle processing (Huang & Davies, 2007) affects our results, using the same setup of fractionation processes as LO-UPb (Table 4).A temporal scaling factor of 13.7 is used, calculated using the method of Huang and Davies (2007) (see Text S1 in Supporting Information S1 for details).To get a similar total removal of Pb from melt to the continental reservoir as in our other simulations, in this simulation this transfer per time-step is reduced by the same factor of 1/13.7.

Temperature Evolution
As the trace element composition of particles does not influence the dynamics of these simulations, the thermal and bulk chemistry evolution is identical for all unscaled cases (Figure 2).Although the dynamics of SCALE will be different to the unscaled cases, many of the features are common.Large hot plumes (Figure 2a) develop in the mantle in all cases, a behavior which is known to occur in models with a free slip surface (Davies, 2005).
The linear cool downwellings are generally more mobile than the plumes, which are relatively stable especially at the base.Figure 3a shows the volume averaged temperature decrease from 2010 K to around 1550 K over the course of the calculation for unscaled cases.This leaves the present day mantle cooling rate (averaged over the last 0.5 Gyr) at 60 K Gyr −1 , slightly lower than estimates of Earth's current mantle cooling rate of 73 K Gyr −1 (Labrosse & Jaupart, 2007).Surface heat flow (Figure 3b) decreases from 54 TW to 39 TW over the same period.This is in line with the current best estimates of Earth's surface heat flow, which is around 39 TW, excluding energy lost from radioactive decay in the continents (Davies & Davies, 2010).It should be noted that while the surface heat flow of the model is a good match with Earth's, the mean surface velocity is considerably lower at 0.73 cm yr −1 compared to ≈5.00 cm yr −1 for Earth (Müller et al., 2008).

Bulk Chemistry Evolution
Figure 2 shows snapshots of the temperature anomaly and bulk composition of the mantle after the final time step of the calculation for all unscaled cases.Melting is concentrated in the heads of plumes (Figure 2a) as these are the hottest regions of the model.The type of melting that we see in these models (described in 2.4) has elements of both MORB and OIB type melting.The plumes transport deep material to the surface where it may melt, as in OIBs, but the melting zone (depth to which particles may melt) is shallow (<135 km) and produces a basaltic crust underlain by a depleted residue layer, like MORB melting.Basaltic material at the surface is pushed laterally as the large upwellings interact with the surface.Where surface material collides spindly regions of downwelling form, akin to subduction zones on Earth.As the melting that we observe in these convection models 10.1029/2021GC010309 9 of 24 is not truly representative of MORB or OIB melting, we shall compare our results against the global average τ Pb and d of MORBs and OIBs (Figure 1).
The delayed phase transition in basaltic material from 660 to 720 km creates a partial barrier to subducted material, preventing some of it from reaching the lower mantle.As a result the transition zone is slightly enriched in basaltic material compared to the ambient mantle.The amount of mantle that has been processed (% of particles that have melted) increases steadily for the first 1.7 Gyr of calculation.Remelting of particles and lower melting rates cause a slower increase in % processed later on in the calculation.At present day 67% of particles have undergone at least one melt event for unscaled cases.Early in the simulation, the large volumes of subducted  basalt reaching the base of the mantle allows small piles to form here.The piles are short lived, quickly being heated up sufficiently to overcome their negative chemical buoyancy and become entrained into mantle upwellings.As the calculation progresses the mantle cools and less basalt is produced, leaving just small accumulations of basaltic material at the base of plumes.The mantle stirs the enriched basaltic, depleted, and unmelted components efficiently, stretching out heterogeneity into thin strands aligned with the direction of flow.Strong stirring coupled with decreasing melting rates leads to a decrease in the wavelength of basaltic accumulations at the surface and within the mantle over time.The Pb isotope signature, being a product of melting and other shallow processes, will to some extent reflect the time integrated signature of processing so provides an indirect way of examining the behavior of the mantle.

Understanding Pb Isotope Ratio Outputs
Figure 4a shows sample output data in 207 Pb/ 204 Pb-206 Pb/ 204 Pb space from case Hi-DPb (fixed partition coefficients).Although this is not a good reflection of the distribution of Pb in modern oceanic basalts (Figure 1), it illustrates a range of the observed model features so is provided to help explain the origin of such features.Each point is the ratio of Pb isotopes carried in a single "melt package"-the information being transferred from a melting particle to particle(s) vertically above at the surface.Generally, there are four areas of the plot that represent melt compositions with distinct melting histories (labeled i-iv).
A high density of melt packages plot in close proximity to a single, central point, labeled "i".This point roughly coincides with the BSE Pb isotope ratio for the present day as given by the input parameters.The melt recorded here is either from particles that have melted for the first time or particles that are remelting after having been recently melted and so have not had enough time to accumulate a significantly different Pb isotope composition.
Melts with Pb isotope compositions less radiogenic that BSE are bounded by a straight edge to the right (high 206 Pb/ 204 Pb, labeled "ii") and an arc to the left (low 206 Pb/ 204 Pb, labeled "iii").The straight edge represents the upper bound of 206 Pb/ 204 Pb, set by variable degrees of melting experienced by BSE-like compositions at the beginning of the calculation.The arc (iii) is for compositions that were following the bulk mantle composition and were depleted at various times, at which point they were almost completely stripped of their U complement.Melts that plot in between these two extents (ii and iii) represent a mixtures of remelts of depleted residues and BSE-like compositions.
Melt packages which have more radiogenic Pb isotope ratios than present day BSE (labeled "iv") are interpreted to be remelts of melts.These are particles that had melt added to them, so they have relatively more U than Pb and so a high μ.
Comparing the isotope ratios of melts in Figure 4a to the MORB data in Figure 1a shows that there is a much greater range of Pb isotope ratios in the model data, with both more and less radiogenic values being recorded.
What is not immediately clear, though, is how much Pb each melt package represents.Given that melt packages with unradiogenic Pb isotope ratios are derived from residues, they will have low concentrations of incompatible Pb. Figure 4b shows the same data but contoured by abundance of 204 Pb (in moles).Points around the BSE value carry significantly more Pb than those further away with more extreme Pb isotope ratios (note the logarithmic color scale).This is significant because the regression is weighted by 204 Pb abundance, so is strongly influenced by moderate Pb isotope ratios.

Fractionation Using Partition Coefficients: MELT, Hi-DPb, XIE007, XIE010
When fractionation of U and Th from Pb is solely controlled by partition coefficients, the mantle has only short wavelength Pb ratio anomalies which are randomly distributed (Figures 5a-5d).These anomalies vary in magnitude depending on the degree of fractionation between U and Pb in each case.In Of these cases, only XIE007 can produce low τ Pb similar to that of the observations, similar to previous modeling (Xie & Tackley, 2004).The range of Pb isotope ratios that are observed in melts is affected by the strength of fractionation (due to the differences between partition coefficients) and the timing of fractionation, with older and stronger heterogeneity contributing to greater scatter (Figures 6 and 7a-7d).None of these simulations, however, can match the scatter observed in oceanic basalts, with d 207 and d 208 being at least a factor of 3 lower than the observations in each case.
In Figure 8 we plot radial averages of μ, 206 Pb/ 204 Pb, and 204 Pb.Only case XIE010 is plotted of the cases which exclude a continental reservoir as each has similar results.The radial average of both μ and 206 Pb/ 204 Pb is constant throughout almost all of the mantle for case XIE010 (Figure 8) while other cases display some radial structure.Particles at the surface generally exhibit high μ, 206 Pb/ 204 Pb and 204 Pb concentrations (Figure 8) while particles in layers down to around 150 km depth display significantly lower values.This trend is attributed to the melting process.Slight increases in μ, 206 Pb/ 204 Pb and 204 Pb concentrations just above the CMB (Figure 8) are a product of chemical density differences, causing particles with a basaltic bulk composition to have a longer than average residence time at the CMB.

Recycling Uranium and Sequestering Lead: CONTU, UPb2.4, UPb0.6, LO-UPb
Including either U recycling or both Pb removal and U recycling changes the shape of the scatter in 207 Pb/ 204 Pb-206 Pb/ 204 Pb space relative to cases where fractionation is controlled by partition coefficients (Figures 6e-6i).There is no pronounced "pinch" near the average composition (as seen in Figures 6b  and 6d) and the Pb isotope ratios are commonly significantly more radiogenic than cases with fractionation only due to differences in partition coefficients.Additionally, when processes for both non-magmatic fractionation are included, d 207 and d 208 are larger than when fractionation is controlled by partition coefficients (Table 5), better matching the observations (Figure 1).The mantle in all unscaled cases develops long wavelength Pb isotope heterogeneities (Figures 5e-5h) due to the strong fractionation offered by non-magmatic processes.
The total amount of U recycled in simulation CONTU (1.7 × 10 15 mol 235 U, 6.9 × 10 16 mol 238 U) is equivalent to ∼20% of the present day global budget of 238 U. Solely recycling U succeeds in producing a τ Pb value that is similar to that of natural samples (Figure 6e), but does not produce an Earth-like distribution in 208 Pb/ 204 Pb-206 Pb/ 204 Pb (Figure 7e), instead plotting to high values of 206 Pb/ 204 Pb.
Lead is extracted at an average rate of ∼3.2 × 10 6 mol yr −1 204 Pb in UPb2.4 (Figure 9), equating to approximately 22% of the global budget being removed from the mantle to the continental crust since 3.0 Ga.There is a greater range of Pb isotope compositions compared to CONTU (Figures 6 and 7f) and the melts tend not to cluster so much around a narrow range of values.τ Pb = 1.88 Gyr, a slight increase from CONTU but similar to the global observed (Figure 1a).There is strong enrichment of U at the surface and in the transition zone relative to 204 Pb (Figure 8a) whilst the melting layer, just below the surface, is more depleted.In 208 Pb/ 204 Pb-206 Pb/ 204 Pb space (Figure 7f) there is a high degree of scatter and the gradient of the regression line through the data (0.65) is low compared to oceanic basalts (1.01).Similar to CONTU, κ m = 2.9, putting it within the range of what can be considered reasonable.The amount of 238 U recycled into the mantle in UPb0.6 corresponds to 18% of the global present day budget.Delaying recycling from 2.4 Ga in UPb2.4 to 0.6 Ga leads there to be a less radiogenic Pb isotope signature overall (Figure 6g).d 207 is largely unchanged from UPb2.4; however, d 208 is reduced by an order of magnitude (Table 5) and is much lower than the observed.The gradient in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space is closer to that of oceanic basalts compared to UPb2.4.UPb0.6 has the strongest U enrichment compared to 204 Pb at the surface of any of the cases (Figure 8a).The very low μ in the melting layers slowly increases with depth and is consistently lower than UPb2.4 until just above the CMB.
With less 238 U being recycled in LO-UPb compared to UPb2.4 (Figure 9), approximately 14% of the present day global budget of 238 U is recycled into the mantle.Lead is extracted at a rate equivalent to approximately 2.5 × 10 6 mol yr −1 of 204 Pb (Figure 9).τ Pb is marginally smaller (1.82 Gyr) than UPb2.4 (Figure 6h) and both d 207 and d 208 are slightly reduced (Table 5), bringing them closer to the observed values.The gradient of the regression line in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space (0.92) is greater than that in case UPb2.4 (Figure 7h) and similar to the observed.In this case κ m of melts is 3.2, falling in between that of UPb2.4 and UPb0.6.

Processing Rate: SCALE
In SCALE, the fraction of particles that have experienced at least one melting event is 97%, compared to 67% for unscaled cases.The physical length scale of Pb isotope anomalies in the mantle is more similar to cases with fractionation controlled by partition coefficients (Figure 5i).The Pb isotope composition of melts is more radiogenic than for case LO-UPb (case with reduced U recycling and Pb removal rates, Figures 6 and 7i), which has similar U recycling and Pb removal rates.6i), but is still similar to that of mantle derived basalts (Figure 1a).The gradient of the regression line in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space is also similar to that of case LO-UPb and the average κ m of melts is slightly increased to 3.3.Radial averages of μ, 206 Pb/ 204 Pb, and 204 Pb concentrations show little variation with depth (Figure 8) but 204 Pb concentrations are slightly lower than those of LO-UPb.

Fractionation via Differences in Partition Coefficients: MELT, Hi-DPb, XIE007, XIE010
As in previous steady state statistical (Allégre et al., 1980) and numerical (Christensen & Hofmann, 1994;Xie & Tackley, 2004)   In each of the cases MELT, Hi-DPb, XIE007 and XIE010, average Pb isotope ratios are less radiogenic (Figures 6a-6d) than the Pb isotope composition measured in oceanic basalts (Figure 1).When including the processes related to the formation (and erosion) of continental crust in CONTU, the majority of melts are more radiogenic than the geochron (Figure 6e) with the continents forming an unradiogenic complement.Despite this, CONTU still does not produce significant amounts of the most radiogenic Pb isotope compositions measured in OIBs ( 206 Pb/ 204 Pb = 21, 207 Pb/ 204 Pb = 15.75, 208 Pb/ 204 Pb = 40).Unlike in Xie and Tackley (2004), which has melt fractionation starting at 4.5 Ga, the low τ Pb of 1.76 Gyr is attained when melt fractionation only begins 3.6 Ga.However, the Pb isotope distribution in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space does not replicate the observations, with significant excess 206 Pb ingrowth compared to 208 Pb (Figure 7e).This is due to particles receiving U from the continents but not Th and so an excess of 207 Pb and 206 Pb develops relative to 208 Pb.Recycling Th with U from the continental reservoir to the mantle is considered very unlikely as Th is not mobilized by oxidative weathering, in contrast to U. Consequently, no significant Th would be added by seafloor alteration to the mafic crust during seafloor alteration in the same way U 6+ is (Hart & Staudigel, 1982, 1989).This does not prevent subduction in nature of Th-enriched sediments; however, from assessing fluxes of U from altered mafic crust and subducting sediment, Elliott et al. (1999) concluded that U would be preferentially returned to the mantle relative to Th.
In cases where Pb is extracted to the continent in addition to U recycling, the Pb isotope signatures of melts are fairly evenly distributed across the data field (Figures 6f and 7f), rather than in previous cases where a large proportion of melts have Pb isotope compositions that cluster around a narrow range of values (e.g., Figures 6a-6e).
The combination of U recycling to particles at the surface of the model and Pb removal from the melt allows radiogenic Pb isotope ratios to grow rapidly, even for particles which haven't necessarily undergone a melting event.Some of the melt produced replicates the most radiogenic Pb isotope compositions that are observed in OIBs (Figure 6f) and, as in CONTU, almost all of the melt is more radiogenic than the geochron.Combining Pb removal with U recycling somewhat remedies the unrealistic distribution of Pb ratios in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space observed in CONTU (Figure 7f); however, the gradient of the regression line is notably shallower than that of the data (Figure 1b) due to an excess of 206 Pb.The 208 Pb/ 204 Pb-206 Pb/ 204 Pb gradient may be tuned to be more Earth-like values by recycling less U (as in LO-UPb).The other end member, where Pb removal is active and we prevent U recycling (see Figure S2 in Supporting Information S1), results in a gradient in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space of 1.1, but d 208 is an order of magnitude lower than the observed, regardless of whether 25% or 10% of Pb is removed to the continental reservoir.When only Pb removal is considered, U and Th are only weakly fractionated during melting, so κ m remains BSE-like (3.9) as in cases where fractionation is controlled by partition coefficients.However, this scenario also cannot reproduce the low τ Pb of oceanic basalts (Figure S2 in Supporting Information S1).
Less U is recycled in case UPb0.6 than UPb2.4,but as this is over a much shorter time period the receiving particles at the surface develop extremely high μ, averaging 43.9 (Figure 8a).Consequently, radiogenic Pb ingrowth on these particles is rapid, and so despite having less time in which to accumulate a radiogenic Pb isotope signature, physically large wavelength variations still develop (Figure 5g).Rapid, late U recycling rates also causes κ m to become unrealistically high (κ m = 3.7) compared to what is measured in MORBs.In our model the onset of U recycling is sudden and the recycling rate decreases with time; however, in reality this is not likely to have been the case.The amount of U being recycled into the mantle may have increased steadily over time between the GOE at 10.1029/2021GC010309 17 of 24 2.4-2.1 Ga (Holland, 1985;Lyons et al., 2014) and full oceanic oxygenation at 0.6 Ga, or may have increased in a step wise fashion (Partin et al., 2013).
Our simplified end member cases UPb2.4 and UPb0.6 do not capture the exact nature of the temporal change that is expected from evidence such as the bimodal Th/U of igneous arc rocks (Liu et al., 2019).An alternative approximation in which the amount of U recycling increases from 2.4 Ga to just after 0.6 Ga may bring the gradient and scatter in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space closer to the observed data values.
Our estimate for the onset of U recycling is fairly well bounded as it is linked to the GOE and full ocean oxygenation at 2.4 Ga and 0.6 Ga respectively (Lyons et al., 2014).However, the history of plate tectonics on Earth is not so well constrained and is still widely debated.For example, the onset of modern style plate tectonics is simultaneously argued to have begun at 700 Ma (Stern, 2005;Stern et al., 2016) and during the Archean (van Kranendonk et al., 2007).We have taken 3.0 Ga as an estimate for the onset of plate tectonics (Laurent et al., 2014;Shirey & Richardson, 2011;Tang et al., 2016), and begin removing Pb from melts from this time as this process is currently intrinsically linked to subduction.If our simulations are sensitive to the time at which preferential Pb removal is initialized then this could provide insight into the time scales over which subduction as we currently understand it has been effective.To this end we set up simulations identical to case LO-UPb varying the time at which Pb is preferentially removed from melts, which we consider to be varying the time of the onset of plate tectonics (Figures S3, S4 in Supporting Information S1).As might be expected, earlier fractionation, which generates ancient heterogeneity, results in older present day τ Pb and greater scatter in the 207 Pb/ 204 Pb-206 Pb/ 204 array (Figure S3 in Supporting Information S1).The 208 Pb/ 204 Pb-206 Pb/ 204 Pb array is only weakly affected by varying the initiation time of plate tectonics between 2.4 Ga and 3.6 Ga (Figure S4 in Supporting Information S1).The "sweet spot" which best matches τ Pb and the observed scatter is when preferential Pb removal begins at 3.0-3.2Ga.This does not necessarily mean that the Pb isotope signature of MORBs and OIBs requires modern day plate tectonics to have initiated in this time window, but does suggest that some process which strongly fractionates U from Pb may have initiated around this time.This could have been a change in the style of continent extraction (Dhuime et al., 2012) or the initiation of a crustal recycling process that would not be classified as modern style plate tectonics (Baes et al., 2020;Simon et al., 2007).Caution must be applied though, as τ Pb is also weakly sensitive to the timing of U recycling (see cases UPb2.4,UPb0.6), and so the U recycling timing and rates must be well constrained.This can be achieved through reconciling the gradient and scatter of the 208 Pb/ 204 Pb-206 Pb/ 204 Pb array (Figures 7f and 7g) and ensuring that the final U concentrations in the continental crust are in line with predictions (Rudnick & Gao, 2013).
The melting zone in case LO-UPb has an average of ∼10 (Figure 8a).This is in agreement with estimates of μ for the upper mantle (Zartman & Haines, 1988).Below the melting zone the radial average μ varies from 11 to 15. Basaltic material, which generally also has a high μ signature, initially segregates to the CMB due to being cold and dense, having formed at the surface, and its intrinsic chemical density in the lower mantle.After a short time the subducted material becomes buoyant enough to be entrained into upwellings and is stirred into the mantle.By the present day there are no significant basaltic accumulations near the CMB (Figure 2b) as the accumulations are less efficiently replenished later in the simulations when mantle temperatures are lower.We therefore see material with a high μ signature distributed throughout the mantle rather than in large, dense piles at the CMB and melts containing a wide range of Pb isotopic signatures.This aligns with recent numerical modeling which argues for fragmented and dispersed recycled domains, rather than large and deep ones, in order to account for Earth's 40 Ar budget (Tucker et al., 2022).The effect of the intrinsic chemical density contrast is to generate a compositional gradient in the lowermost mantle (Figure 10), as preferred by Albarède and van der Hilst (2002), rather than distinct chemical stratification.The full suite of Pb isotope compositions, including 10.1029/2021GC010309 18 of 24 HIMU, can be accounted for without the need for a discrete, long-lived reservoir of subducted oceanic crust as in, for example, Christensen and Hofmann (1994) and Brandenburg and van Keken (2007a).Instead, in addition to the fractionation on melting incorporated in all these models, we implement the widely accepted processes of U recycling from the continental crust and preferential removal of Pb from subducted oceanic crust.The combination of these processes also provides the closest match to the scatter of the observations (d 207 ).We note, however, that some researchers have argued for alternative mechanisms of producing HIMU-like basalts.A lithospheric source of HIMU (Homrighausen et al., 2018;Weiss et al., 2016) is not readily related to the processes modeled in our study, but a proposed model of recycling U in subducted carbonate (Castillo, 2015) is similar to our CONTU scenario.
Despite the convective vigor in our models (Ra ≈ 2 × 10 7 Pa s) being lower than what is expected of the Earth, the similarity between the results of cases LO-UPb and SCALE suggests that further processing by melting would not significantly affect our results.Efficient stirring and a high rate of mantle processing in SCALE results in large scale chemical homogeneity in the mantle (Figures 5i and 8).However, at the finer scale the mantle remains very much heterogeneous (Figures 6 and 7i), suggesting that processing by melting alone cannot eliminate mantle heterogeneity.
Results for case SCALE show mantle 204 Pb concentrations that are slightly lower than in case LO-UPb (Figure 8c).This is due to our simplistic implementation of Pb removal where an arbitrary fraction is removed from the melt.Despite this fraction being decreased by a factor equal to the scaling factor, more 204 Pb is removed in SCALE than in LO-UPb.Consequently, the melts in case SCALE have more radiogenic Pb isotope ratios than in LO-UPb and oceanic basalts.This may also be influenced by the efficient stirring of the upper mantle in case SCALE, which prevents the melting layer from developing an unradiogenic Pb isotope composition as happens in un-scaled cases (Figure 8b).Material sampled from the melting layer therefore has a similar composition to the rest of the mantle.The more radiogenic material has an old Pb isotope signature, and hence τ Pb is slightly older in case SCALE than in LO-UPb.With carefully chosen parameters the same amount of Pb could be removed from the mantle in cases LO-UPb and SCALE.
It is worth noting that while our model can provide a good match to global scale averages, it falls short of being able to reconcile the Pb isotopic variations observed between different ocean basins (Hofmann, 2003).This is to be expected, not least because our convection simulations do not have an Earth-like configuration of plates at the surface.While it would be possible to run mantle circulation simulations driven by plate motion reconstructions with the same non-magmatic fractionation processes that we have presented, the temporal extent of these reconstructions (1 Ga (Merdith et al., 2021),) would limit their use.However, in fitting the global scatter we implicitly reproduce the full Pb isotopic range observed in the ocean basins, including the distinctive compositions found in the Southern Hemisphere, often dubbed DUPAL (Hart, 1984).This is achieved without the need for sediment input, which is often argued for in explaining DUPAL compositions (Rehkämper and Hofmann (1997)).While intriguing, caution is advised with such an interpretation, again due to the lack of an Earth-like tectonic regime and also as we do not track Nd and Sr isotope ratios.

Distribution of Melt Ages
Statistical box modeling by Rudge et al. (2005) and Rudge (2006) approached the problem of the observed τ Pb by relating it to the remelting rate of the mantle.Assuming a well mixed mantle, the distribution of melting ages (time since a particle last melted) may be interpreted as a probability density and can be used to calculate the pseudo-isochron age (τ PDF ) (Rudge, 2006).For our model HI-DPb, there is a good correlation between τ PDF calculated from the distribution of melting ages (Figure 11a) and τ Pb calculated from the Pb isotope distribution throughout the calculation (Figure 11b).The greatest misfit occurs in the first ∼500 Myr, before the mantle is well mixed.This is a similar result to that presented in van Heck et al. (2016), in which a similar model setup was used.The dynamics of HI-DPb and, for example, LO-UPb, are identical but due to the additional processes of U recycling and Pb removal there is a poor match between τ PDF and τ Pb for case LO-UPb (Figure 11b).Recycling of continental U and extraction of Pb from melts involves communication with an external reservoir (the models of Rudge et al. (2005) and Rudge (2006) are closed systems), which means Pb isotope heterogeneity is no longer purely a product of melting, and Pb isotope heterogeneity is generated at many length scales (Figure 5g).Consequently, for this class of models the relationship between melt age distribution and τ Pb does not apply.

Highly Radiogenic Pb
As seen in visualisations (Figure 5) and Pb isotope ratio plots (Figures 6 and 7) some particles in the system develop extreme radiogenic Pb isotope ratios in cases which include U recycling and Pb loss to the continents.Particles which have melted multiple times become strongly depleted in Pb, causing them to develop high μ.This may also be compounded by the same Pb depleted particles receiving recycled continental U. Rapid radiogenic Pb enrichment relative to 204 Pb follows, leading to the extremely radiogenic Pb isotope compositions.For example, in case LO-UPb 17.7% of melt parcels have 206 Pb/ 204 Pb > 25, which is the upper end of what is measured in OIBs.However, the concentration of Pb on these particles is low and represents just 2.3% of the total Pb in the melt.In case SCALE this is <0.05% of the total Pb in the melt.It is also worth recognizing that our models do not incorporate the shearing and diffusion processes (Kellogg & Turcotte, 1987) that may help to equilibrate strongly radiogenic compositions.
Although the most radiogenic Pb ratios produced by our models are at odds with the observations, the length scale of heterogeneity mixing involved to generate OIBs is poorly known.For example, in the statistical upper mantle assemblage (SUMA) model of Meibom and Anderson (2003), various degrees of mixing accounts for differences in the range of isotope ratios observed in different eruptive environments.As an averaging process, mixing of melts has the effect of reducing the spread of Pb isotope ratios, in the case of SUMA bringing model ratios closer to those of MORBs.Recent Nd isotope analysis of lower crustal pyroxene and plagioclase, obtained from drilling of the Mid-Atlantic Ridge, has shown that higher degrees of chemical heterogeneity exist than have been inferred from sampling erupted MORBs, even on the crystal scale (Lambart et al., 2019).Similarly, sulphide inclusions within MORBs have been shown to display Pb isotope compositions outside the range observed in whole rocks (Burton et al., 2012).We approximate magma mixing by averaging a number of melts (N) produced close to one another, similar to the Rudge (2006).When N = 50, the fraction of melt packages with extremely high 206 Pb/ 204 Pb in case LO-UPb falls to 4.6% (from 17.7%) while the fraction of Pb that is carried on these particles decreases to 0.8% (from 2.3%).Rudge (2006) showed that τ Pb has a weak dependence on N. In our models we generally find mixing to cause τ Pb to increase.For N = 50, case LO-UPb τ Pb increases from 1.82 to 2.02 Gyr (Figure S5 in Supporting Information S1).At the same time, d 207 is reduced from 0.0313 to 0.0291 and d 208 from 0.2339 to 0.1522 (Figures S5, S6 in Supporting Information S1), bringing the scatter more in line with observations.

Composition of the Continental Reservoir
Lead fluxes between the mantle and continental crust are poorly constrained but the relative distribution of Pb between these two reservoirs is better known.For example, Rudnick and Gao (2013) estimate that 43% of Earth's Pb budget resides in the continental crust, though this comes with large uncertainties.This can be used as a constraint to ensure that a reasonable amount of Pb ends up in the continental reservoir, not forgetting that we have omitted the core as a reservoir which is likely to contain unradiogenic Pb incorporated as sulphides (Kramers & Tolstikhin, 1997;Maltese & Mezger, 2020) or metallic melts (Ballhaus & Ellis, 1996).Additionally, Pb stripped from subducted oceanic crust may not end up in either the continental crust or the core, but in the lithospheric mantle, isolated from the convecting mantle (Halla, 2005;King et al., 2007).For convenience, we have called the reservoir into which all Pb removed from melts are transferred, the continent reservoir but acknowledge that other reservoirs may also account for a portion of the Pb removed from the mantle.Nonetheless, we can make some first order observations on the composition of our model continental reservoir compared to Earth's continental crust.Compared to the value of 43% for the continental crust, the fraction of Pb in the continent at present day is too low in CONTU but higher in cases UPb2.4,UPb0.6, LO-UPb and SCALE (Table 6).
Recycling U and sequestering Pb causes the continental reservoir to develop unradiogenic ratios of 206 Pb/ 204 Pb and 207 Pb/ 204 Pb (Table 6).This is complementary to the radiogenic compositions found in recycled material.Previous proposals that the continental crust could host unradiogenic Pb have suggested that it exists in the lower crust due to emplacement of Pb that is preferentially stripped from subducted oceanic lithosphere (Kramers & Tolstikhin, 1997;Zartman & Haines, 1988), something that our model cannot directly determine as we only derive bulk continental crust values.Bulk continental crust 206 Pb/ 204 Pb and 207 Pb/ 204 Pb for all cases in Table 6 are lower than the estimates of Rudnick and Goldstein (1990) for the Pb isotope composition of the bulk continental crust, and more similar to estimates of the Pb isotope composition of the lower continental crust (Rudnick & Goldstein, 1990;Zartman & Haines, 1988).For each model case the bulk continental crust 208 Pb/ 204 Pb is less radiogenic than estimates for the upper continental crust (Zartman & Haines, 1988).
Cases CONTU, UPb2.4 and UPb0.6 each have a bulk continental κ m (Table 6) which is significantly higher than estimated values of ∼4.55-5.0(Paul et al., 2003;Wedepohl, 1995;Zartman & Haines, 1988).Reduced U recycling in cases LO-UPb and SCALE decreases the κ m of the continental reservoir (to κ m = 5.84), though it remains slightly higher than the estimates.Chemical modeling by Kramers and Tolstikhin (1997) suggests μ of 4.16 and 10.2 for the lower and upper crust, respectively.Separately, the bulk continental crust μ has been estimated by Allègre and Lewin (1989) to be 9.58 ± 1.In each of the cases we present the continental reservoir develops a bulk μ which is much smaller than these estimates, smaller even than estimates for solely the lower continental crust.
Of the cases presented, LO-UPb best reflects geochemical estimates of the Pb, μ, and κ m composition of the continental crust.The biggest discrepancy is the μ of the continental reservoir, which is just a third of the estimated value (Allègre & Lewin, 1989).To some extent this may reflect our omission of the core as a geochemical reservoir, which is likely to comprise unradiogenic Pb due to early core formation (Maltese & Mezger, 2020;Wood & Halliday, 2005).Such early Pb extraction has also been invoked as a solution of the first Pb paradox (Hart & Gaetani, 2006;Kramers & Tolstikhin, 1997).

Limitations
The models presented here neither have Earth-like plates geometries at the surface (as discussed in Section 4.2) nor rheology which is dependent on temperature or composition, therefore the details of actual mantle flow may be expected to be different.We might expect Earth-like plates to lead to more effective mixing (van Keken et al., 2002) due to toroidal flow of strike-slip boundaries (Ferrachat & Ricard, 1998), but that could be countered by the fact that subducting slab curtains can reduce mixing (Barry et al., 2017).More realistic rheology would inhibit the mixing of subducted material until it has heated up and become sufficiently weak (Zhong et al., 2000) which could affect the Pb isotope ratio distribution in the mantle.Plumes would become narrower, which could limit how much material they could transport to the surface, however their velocity would also increase which would counter this to some extent.While our models do include an element of secular cooling (Figure 3a), we have simplified the CMB boundary condition by imposing a fixed temperature for the 3.6 Gyr of modeled time.A more Earth-like setup would be to include a cooling CMB which would cause hotter mantle temperatures in the past leading to more rapid processing.Our models do not include the process of diffusive homogenization (Kellogg & Turcotte, 1987) which can act when anomalies are sufficiently thinned by shearing (Xie & Tackley, 2004).It is unclear whether the oceanic crust can be sufficiently thinned for this to play a significant role for U, Th and Pb, which might be expected to have relatively low diffusion coefficients.We have neither included nor required a primitive layer at the CMB.
It is currently unclear whether the class of model presented here will be capable of successfully modeling the signatures of other isotopic systems.More sophisticated models will be required to address all of these questions.

Conclusions
Using 3-D mantle convection simulations we have demonstrated the important role that recycling of continental U and preferential removal of Pb from oceanic crust have on generating the characteristic range of Pb isotope ratios observed in oceanic basalts today.
When U -Pb fractionation is solely controlled by the partition coefficients, a relative change in behavior between U and Pb is required at around 2.4 Ga in order to reproduce the low τ Pb of oceanic basalts, as in Xie and Tackley (2004).However, the scatter achieved does not satisfy the full range of Pb isotope ratios that are observed.
If U recycling acts alone then the low τ Pb observed in MORBs and OIBs can be reproduced, but excess ingrowth of 206 Pb relative to 208 Pb causes an irreconcilable mismatch between the modeled 208 Pb/ 204 Pb-206 Pb/ 204 Pb array of melts and oceanic basalts.Likewise, simulations in which preferential Pb removal acts alone cannot reproduce the low κ m or τ Pb that is observed in nature.Simulations which include both of the above processes provide the best matches to multiple constraints, including τ Pb , the gradient of the 208 Pb/ 204 Pb-206 Pb/ 204 Pb array, the scatter d 207 and d 208 and the average κ m of melts.Unlike much previous work, our models do not require deep, large-scale, long-lived reservoirs of subducted oceanic crust in order to reconcile geochemical observations.Instead, material with a strongly radiogenic Pb isotope signature occurs at many length scales throughout the mantle.Rather than chemical stratification, the intrinsic density contrast generate a chemical gradient toward the base of the mantle.
The scatter and gradient of the 208 Pb/ 204 Pb-206 Pb/ 204 Pb array are sensitive to the timing of U recycling, so have the potential for helping to constrain this.For our preferred U recycling parameters our model estimates that a style of subduction which preferentially removed Pb from subducted oceanic crust over Th and U started at 3.0-3.2Ga.Future work should not neglect the core as a reservoir for unradiogenic Pb, this would help to achieve a more Earth-like composition of the continental crust.This research was partially funded by the NERC funded consortium, "Mantle volatiles: processes, reservoirs and fluxes" (Grant No. NE/M000397/1) of the Deep Volatiles Programme, the NERC large grant "Mantle Circulation Constrained (MC2): A multidisciplinary 4D Earth framework for understanding mantle upwellings" (Grant No. NE/T012633/1), and also by the School of Earth and Environmental Sciences, Cardiff University.Numerical calculations were undertaken at: (a) ARCHER2, the UK's national high-performance supercomputer; (b) HAWK, part of Supercomputing Wales, the national high-performance supercomputing system for Wales.Graphs were produced using the Matplotlib package for Python (Hunter, 2007) and we make use of the scientific color map, batlow, for select figures (Crameri et al., 2020).Visualizations were produced using Paraview (Ahrens et al., 2005).This is Cardiff CRediT Contribution 2.

Figure 2 .
Figure 2. Volume slices taken after the final time step of unscaled cases.(a) Colored by temperature anomaly with iso-surface for dT ≥ +600 K, where dT is the difference in temperature from the radial average.View is clipped at 25 km depth as surface temperature is fixed at 300 K (b) Same volume slice colored by bulk composition (C) White on the color scale represents ambient mantle composition of C = 0.25.The purple colors indicate harzburgitic material while green colors show basaltic material.Surfaces are drawn for regions with C ≥ 0.9.

Figure 3 .
Figure 3. (a) Volume averaged temperature of the mantle, for unscaled cases with constant internal heating rate, over time.(b) Surface heat flow (solid line), core-mantle boundary heat flow (dashed line) and internal heating (dotted line) over time.

Figure 4 .
Figure 4. (a) Scatter of Pb isotope ratios of melt packages from the final time step of case Hi-DPb.Regression line calculated using the orthogonal distance method, weighted by the abundance of 204 Pb.The mid point of the data is taken as the median value in 207 Pb/ 204 Pb and 206 Pb/ 204 Pb.See text for description of annotations i-iv.(b) Data from (a) contoured by 204 Pb abundance within each cell of a 100 × 100 grid.

Figure 6 .
Figure 6. 207Pb/ 204 Pb-206 Pb/ 204 Pb of melt packages in cases (a) MELT, (b) Hi-DPb, (c) XIE007, (d) XIE010, (e) CONTU, (f) UPb2.4,(g) UPb0.6,(h) LO-UPb, and (i) SCALE.Axes are divided into a 100 by 100 grid for color-coded contouring by 204 Pb absolute abundance.Dark dashed line is the pseudo-isochron for each case and light dashed line is pseudo-isochron for MORB data as in Figure 1a.The red dashed line is the 4.55 Gyr geochron for the initial isotopes used in these calculations.d value is the average orthogonal distance of each melt package from the pseudo-isochron.

Figure 7 .
Figure 7. 208 Pb/ 204 Pb-206 Pb/ 204 Pb of melt packages in cases (a) MELT, (b) Hi-DPb, (c) XIE007, (d) XIE010, (e) CONTU, (f) UPb2.4,(g) UPb0.6,(h) LO-UPb, and (i) SCALE.Axes are divided into a 100 by 100 grid for color-coded contouring by 204 Pb absolute abundance.Gradient is the gradient of the regression line through the data (dark dashed line).The value κ m is mean 232 Th/ 238 U of all melt produced at the final time step.Light dashed line is the regression through the mid-ocean ridge basalts data as in Figure 1b.

Figure 9 .
Figure 9. Flux of 238 U into the mantle (left panel) and 204 Pb out of the mantle (right panel) for cases UPb2.4,UPb0.6 and LO-UPb.Note that the 204 Pb flux out of the mantle is identical for cases UPb2.4 and UPb0.6.

Figure 10 .
Figure 10.Radial average of the bulk composition at present day for case LO-UPb.

Figure 11 .
Figure 11.(a) Distribution of melting ages (time since last melting) at present day for particles which have undergone melting in case HI-DPb.Note that the melting age distribution is identical for all cases except SCALE due to similar dynamics.(b) Comparison of τ PDF , calculated from the distribution of melting ages (blue circles) and τ Pb calculated from the Pb isotope composition of the melt in cases HI-DPb (orange circles) and LO-UPb (green circles) at different times in the simulation. μ T, time t, thermal diffusivity κ, radiogenic heat production H, specific heat at constant pressure C p , and bulk composition C. Other model parameters are listed in Table1.The simulations presented have been conducted on a spherical mesh with 65 layers consisting of over 10 million nodes, giving a radial and average lateral resolution of 45 km.We use a simple 2 layer vertical viscosity profile with a × 30 viscosity jump at 660 km (van Keken & Ballentine, 1998) so

Table 1 Model Parameters Geochemistry, Geophysics, Geosystems
of the upper mantle and lower mantle are 3 × 10 22 Pa s and 9 × 10 23 Pa s respectively.The lack of a viscous lithosphere gives the models a mobile surface to approximate the mobility of plate tectonics.Both surface and core-mantle boundary (CMB) are free-slip, impermeable and isothermal.The model is internally heated homogeneously at a constant rate (Table

Table 2
Olivine Phase Change Parameters for an Assumed Composition With 67% Olivine PANTON ET AL. 10.1029/2021GC010309 6 of 24 207Pb/ 204 Pb-206 Pb/ 204 Pb space, there is central point of highest Pb density at 206 Pb/ 204 Pb = 17.5 and 207 Pb/ 204 Pb = 15.43 in each of these four cases (Figures6a-6d).This point roughly falls on the 4.55 Gyr geochron, and approximately equal proportions of melt compositions plot with more and less radiogenic values.The similarity of κ m (mean 232 Th/ 238 U of melts in final time step ≈ 3.9) and the gradient in 208 Pb/ 204 Pb-206 Pb/ 204 Pb space (Figures7a-7d) is because U and Th are only weakly fractionated.
Compared to case LO-UPb the scatter in both 207 Pb/ 204 Pb-206 Pb/ 204 Pb and 208 Pb/ 204 Pb- 206Pb/ 204 Pb space is reduced.τ Pb is increased (2.01 Gyr) relative to LO-UPb (Figure modeling, our simulations MELT and Hi-DPb cannot reproduce the τ Pb observed in oceanic basalts.Increasing fractionation does, however, have a strong effect on the scatter observed in the melt (Figures 6, 7a and 7b) because of the relatively high U/Pb of the enriched basaltic crust and relatively low U/Pb of the depleted residue.Pb Is the Pseudo-Isochron Age, grad 208 Is the Gradient of the Regression Line for Pb Isotopes Plotted in 208 Pb/ 204 Pb-206 Pb/ 204 Pb Space, d 207 and d 208 Is the Average Distance Each Point Plots Away From the Regression Line in 207 Pb/ 204 Pb-206 Pb/ 204 Pb and 208 Pb/ 204 Pb-206 Pb/ 204 Pb Space Respectively (Weighted by 204 Pb Abundance) and κ m Is the Average 232 Th/ 238 U Measured in the Melt.
developing.Consequently, scatter orthogonal to the pseudo-isochron is diminished compared to Hi-DPb, and so the low d 207 and d 208 do not match the observations.The early fractionated material in XIE010 can accumulate a wider range of Pb isotope signatures (compared to BSE) before being re-sampled, resulting in more scatter compared to XIE007, but unsurprisingly this also increases τ Pb to 2.59 Gyr.Nonetheless, no cases in which fractionation is solely controlled by partition coefficient can reproduce both τ Pb and the observed scatter.Figure 8.(a) Plots of layer-averaged μ with depth for cases XIE010, CONTU, UPb2.4,UPb0.6, LO-UPb, SCALE.(b) Layer-averaged 206 Pb/ 204 Pb with depth.(c) Layer-averaged 204 Pb concentration with depth.

Table 5
Model Results for Some Key Characteristics of the Pb Isotope Composition and Distribution in Melts With Observed Values for Comparison.