Global scale modeling of melting and isotopic evolution of Earth's mantle

H. J. van Heck, J. H. Davies, T. Elliott, and D. Porcelli School of Earth and Ocean Sciences, Cardiff University, Main Building, Park Place, Cardiff, CF10 3AT, Wales, UK Department of Earth Sciences, Utrecht University, Heidelberglaan 2, 3584 CS Utrecht, the Netherlands Department of Earth Sciences, University of Bristol, Wills Memorial Building, Queen’s Road, Bristol BS8 1RJ, UK Department of Earth Sciences, University of Oxford, South Parks Road, Oxford OX1 3AN, UK

convection while tracking chemical information throughout the convective mantle.
We have implemented a new way to track both bulk compositions and concentrations of trace elements in a finite element mantle convection code. Our approach is to track bulk compositions and trace element abundances via particles. One value on each particle represents bulk composition, and can be interpreted as the basalt component. 10 In our model, chemical fractionation of bulk composition and trace elements happens at self-consistent, evolving melting zones. Melting is defined via a composition-dependent solidus, such that the amount of melt generated depends on pressure, temperature and bulk composition of each particle. A novel aspect is that we do not move particles that undergo melting; instead we transfer the chemical information carried by the particle 15 to other particles. Molten material is instantaneously transported to the surface layer, thereby increasing the basalt component carried by the particles close to the surface, and decreasing the basalt component in the residue.
The model is set to explore a number of radiogenic isotopic systems but as an example here the trace elements we choose to follow are the Pb isotopes and their radioac-In the Earth's mantle, chemical heterogeneities in bulk composition and trace element concentration and isotope composition are continuously created by melting. Oceanic crust is produced by partial melting at oceanic spreading centres where most mantle melting occurs, and also where most chemical heterogeneity is generated. This heterogeneous material is brought into the deeper mantle via subduction of oceanic 15 lithosphere. Here it mixes. To a lesser extent, melting also happens on continents and beneath oceanic lithosphere to create ocean island basalts (OIB), providing a second mechanism for creating heterogeneity. This makes melting a first order feature to be implemented in thermo-chemical convection codes. In addition to this continuous generation of heterogeneities, chemically distinct material might have survived for bil- 20 lions of years, originating much earlier in Earth history e.g. linked to core formation processes, mantle magma oceans, or asteroid bombardment.
Numerical mantle convection codes have been developed that are capable of tracking chemical heterogeneities (for an overview see Tackley, 2007). Although different techniques each have advantages (level-set Samuel and Evonuk, 2010,  In this paper, we will deal with global scale convection, melting, and the tracking of heterogeneities resulting from melting. Christensen and Hofmann (1994) were the first to demonstrate a method to track the evolution of recycled oceanic crust and its influence on the chemistry of the mantle. After that study was published, many followed a similar approach, either tracking preset heterogeneities (e.g. Davies, 2002;5 Zhong and Hager, 2003;Nakagawa and Tackley, 2004), or having the chemical heterogeneities emerge via melting during the calculation at fixed melting zones (Walzer and Hendel, 1999;Davies, 2002;Huang and Davies, 2007a, b, c), moving melting zones that follow force-balanced plates or imposed plate motions (Brandenburg and van Keken, 2007a;Brandenburg and Van Keken, 2007b;Brandenburg et al., 2008), or freely moving melting location via a melting phase diagram (De Smet et al., 1998;Van Thienen et al., 2004;Xie and Tackley, 2004a, b;Nakagawa et al., 2009Nakagawa et al., , 2010. In addition to tracking of bulk compositions, it is also possible to track the distribution and evolution of trace element abundances. Of particular interest are both parent and daughter isotopes of radiogenic systems. Since the radiogenic parents decay in a very 15 predictable manner (following known decay constants) their relative abundances compared to their daughter isotopes can be used as clocks. Different elements behave differently during partial melting (via different partition coefficients), so when the techniques of tracking and melting are implemented the system of segregation and decay can be tracked. We follow the well studied Uranium-Thorium-Lead system (U-Th-Pb), 20 as has previously been done in a similar way in: Christensen and Hofmann (1994); Xie and Tackley (2004a); Brandenburg et al. (2008).
For purely thermal convection, numerous analytical solutions exist, which can be used to benchmark numerical codes. (For example Rayleigh-Taylor instabilities; Nu-Ra scalings; corner flow.) In this way, numerical codes can be checked for accuracy 25 for simple setups and then used to study more complicated scenarios. For thermochemical convection on the other hand, very few analytical solutions exist. This makes the benchmarking of thermo-chemical convection difficult. There is one semi-analytical theory available for the evolution of Pb isotopes in the mantle (Rudge, 2006 use this theory and compare its predictions to the outcome of newly performed numerical calculations.
In this paper, we present the details of a newly implemented method of melt and chemical heterogeneity tracking in a global scale convection code. We will show a good fit of our model data to the analytical prediction of Pb-isochron ages as function of 5 melting age distribution functions. Through this we validate our implementation.

Methods
In this section we first describe the numerical code and calculation setup before presenting more details about using particles to track bulk and trace element compositions. Then we move on to describe the implementation of melting, including changes in man-10 tle bulk compositions and fractionation of trace elements. We conclude this section by providing some details about the initial setup of the trace elements, as done for the calculation presented.

Fluid flow convection
We use TERRA, a parallel, well-established, benchmarked, spherical convection code 15 (Baumgardner, 1985;Bunge and Baumgardner, 1995;Yang and Baumgardner, 2000;Stegman et al., 2002Stegman et al., , 2003Köstler, 2011;Davies et al., 2013). The grid covers the full 3-D spherical shell. At each radial layer the grid is a regular subdivision of an icosahedron (Baumgardner and Frederickson, 1985 Assuming incompressibility and the Boussinesq approximation the equations for convective flow in the mantle can be expressed non-dimensionally as: where length is non-dimensionalised by D the depth of the mantle; time is nondimensionalised by D 2 κ −1 (κ the thermal diffusivity), and temperature by T the tem- , where H is the heat generation rate per unit volume and c p is the specific heat at constant pressure. Material movement is driven by buoyancy forces resulting from horizontal differences in density as expressed in the momentum equation (Eq. 2). Temperature is advected with the material flow, diffuses 15 and is produced internally as described by Eq. (3); Eq. (3) describes the advection of temperature with the material flow, and its diffusion, and the heat produced internally; while Eq. (1) is the continuity equation which ensures conservation of mass. Equations (1) and (2)  To mimic most of Earth's evolution we ran the calculation over a time corresponding to several billions of years. The calculation started at a time 3.6 Ga and ran forward until present day. By doing so we skip the first billion year of the Earth's evolution. This is 5 done because in Earth's early evolution, mantle temperatures were most likely higher than they are at present, leading to lower viscosity and higher vigour in the convection. Numerically this type of convection would be harder to solve for accurately. Moreover, the style of convection and so whether or not it can be treated using the same sort of model is also uncertain. We note that the viscosity is not temperature-dependent in the 10 actual simple case presented here, and so the starting time could be older.

Particles and bulk composition
We use particles for tracking chemical information, and thus for dealing with melting. A schematic illustration of how this works is given in Fig. 1. The particles are advected using a second order Runge-Kutta advection scheme. Each particle is linked to an 15 array to carry the particle information. Three values are used to indicate the particles position; one to indicate the mass that the particle represents; one to track at what time the particle last produced or received melt; one to track the bulk composition; and one for the abundance of each isotope which is tracked in the system. Using the setup presented in this paper that means that each particle is linked to an array with 13 num-20 bers. The mass that a particle represents is attributed at initialisation as the volume of the node, multiplied by the density, divided by the number of particles attributed to that node. After initialisation, the mass a particle represents is only changed when the particle is split, or merged with another particle (i.e. for numerical, not physical reasons. See Sect. 2.4 for more details.). Bulk composition is tracked by a single number 25 (C value) for each particle. The C value can vary between 0 and 1. Zero we interpret to mean effectively completely infertile (in basalt), and 1 to mean completely fertile (in basalt). We use a simple set-up in the actual test that we present, described in the sec-9559 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tion above. This includes on initialisation giving all particles a value of 0.6, representing a homogeneous partly fertile mantle.

Splitting and merging
In order to keep proper coverage of particles throughout the mantle, we sometimes need to split and merge particles. Particles were split whenever less than a threshold 5 number of particles are present in a grid cell. The threshold for the calculations presented here was three. All particles in such cells are split in half, creating two new particles each representing half the original mass. One particle copies the position of the old particle while the other is placed in the same grid cell, mirrored over the z axis.
Each new particle receives half of the atoms of each isotope the old particle was car-10 rying, while bulk composition and melting time are simply copied.
Merging of particles is needed sometimes as well. The threshold number of particles required in a grid cell before merging was undertaken can again be varied and was set at 35 in the calculations presented here. When 2 particles merge, their mass and isotopes are simply summed. Which 2 particles are merged is determined by their 15 location in the array allocated to the node, so effectively they are selected ad random. The position of the new particles is the average location determined by weighting the old positions according to the particle masses. The bulk composition is also the mass weighted average of the two original compositions. The melting time of the new particle is copied from one of the old particles picked at random. We note that these rules 20 for splitting and merging particles conserves the global bulk composition, mass and isotope abundances.

Trace elements
We track the evolution of different trace elements through the domain, focussing on several isotopes of Uranium (U), Thorium (Th) and Lead (Pb). 238  and 207 Pb/ 204 Pb with time, but the slope of the correlation between these two ratios defines the time of U-Pb fractionation. In reality, it is unlikely the mantle comprises discrete reservoirs, fractionated at the same time and more plausibly a mixture of sources fractionated at variable times. In this case the slope of an array of 206 Pb/ 204 Pb and 207 Pb/ 204 Pb ratios for mantle-derived samples still carries age information about man-10 tle evolution (e.g. Allègre et al., 1980), but interpretation of such "pseudo-isochrons" (see Rudge, 2006) is more complex. One potassium isotope ( 40 K) is also tracked since, next to U and Th, 40 K is the isotope that generates the bulk of the mantle's internal heating.
For all radioactive parents we first estimate the total amount of each isotope at present day. That is the amount in all reservoirs combined, i.e. the total budget for the Earth. After that is done, we estimate the amount at the start of the calculation (3.6 Ga ago) by adding the amount that has decayed since then, via standard exponential decay (Eq. 4): 20 where X can be either U, Th, or K; X s is the abundance at the start of the calculation; X pd the abundance at present day; λ is the decay constant, and t is the time between the time at the start and present day.

Melting
We now describe the melting algorithm starting with an overview. The melting algorithm 25 is implemented on particles whose temperature exceeds their solidus. The amount of melt and its content of trace elements is calculated for each such particle. This information is then passed to near surface particles conserving energy, mass, bulk composition and atoms of trace elements. Note that, in contrast to many other implementations, the melting particles are not moved as part of the melting event. We next describe the general assumptions underlying our algorithm, and then describe the specific choices and 5 give more details. Our melting relationships follow from 3 assumptions. The first assumption is that the proportion of fusible (or basaltic) material in a particle can be represented by a compositional parameter C. Following a melting event, all or part of this fusible component is removed from the melt-producing particle, thereby depleting it. The degree of melting, 10 F , is given directly by the change in composition C as follows, where, C n is the new bulk composition and C o is the previous bulk composition, as recorded at the particle. The second assumption is that the solidus is dependent on the compositional pa-15 rameter, C. i.e. T s = f (C), where T s is the solidus temperature, and f (C) is a function of C. Physically the function f (C) must be monotonic, i.e. the solidus temperature must increase steadily as the composition becomes more depleted. Since f (C) is monotonic, its inverse function, f −1 = g(T s ) (also monotonic) exists. The function g gives the composition, C, as a function of the solidus temperature C = g(T s ). 20 The third assumption is that, following the melting event the temperature of the particle will be the temperature of the solidus for its new composition; this is achieved by changing the composition, not the temperature. Melt can be extracted until the residual composition is so refractory (C = 0), that no further melting occurs.
Using these assumptions, the degree of melting is explicitly calculated as follows. 25 First, at each time step, temperature (T ) is interpolated from the grid to the particles using linear interpolation (using barycentric-based finite element shape functions). Then using the temperature difference between the composition-dependent solidus and the 9562 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | actual temperature on the particle, we calculate the new composition of the residue C n (assuming it is in thermal equilibrium with its new solidus (T ); and only when the temperature exceeds the solidus.) Then using Eq. (5) we can calculate the degree of melting, F , using the new compo-5 sition. For numerical reasons we set a threshold for F of at least 0.0001; degrees of melting lower than this are ignored. For this work we make the simplifying assumption that the functions relating solidus temperature to composition and the inverse are linear; i.e. functions f and g are linear. This simplification is justified for this work since our goal is to demonstrate and then test the method with a simple semi-analytical model. We also assume that the solidus temperature is a function of pressure, and make the reasonable, simplifying assumption for this work, that it is a linear relationship (see Fig. 2). As a result we get that where, C n is the new bulk composition; T the temperature at the particle; T m−1 the melt-15 ing temperature at the surface for material of composition of C = 1; z the depth of the particle; dT m /dz the slope of the solidus; and T Comp the compositionally dependent temperature difference between the solidii for material of composition zero and material of composition one. The total amount of melt is calculated by multiplying the mass of the particle with 20 the degree of melting. This melt (basalt, C = 1) is extracted from the melt producing particle and brought to one or more particles close to the surface. We consider the particles at the surface, directly above the melt producing particle, and keep updating their composition towards pure basalt until all melt is stored. We do this starting at the surface layer of the grid and continue downwards to underlying layers if required such that a layer of pure basalt forms and all the melt is accommodated. Since not all particles represent the same amount of mass, we are careful to ensure that all melt is stored and mass conservation is obeyed. Note that the melting/residue particle keeps the same mass. This reflects the fact that the melting column subsides (compacts) in the following way: Since basalt is removed, and the particle's mass is conserved, 5 the removed material is implicitly replaced by depleted material (of composition C = 0, because fusible component is linearly linked to degree of melting, Eq. 5). The basalt is stored higher in the column, where it takes the space of fully depleted material. Any intervening layers are unchanged. Note, after melting that while the temperature of the melt producing particle in the 10 residue is unchanged (and hence energy is conserved in the melting event), it is in thermal equilibrium with its new solidus since its composition has changed appropriately. The implementation presented here does not take the effect of latent heat or thermal advection by melt migration into account. Neglecting this will only have an effect on the thermal evolution of calculations, and then only ones with massive magmatism. The 15 effect on the chemical evolution, considered here, will be minimal. This aspect could be added to the model if future applications required it. We then bring the trace isotopes with the melt to the surface. For this we assume simple batch melting partitioning, 20 where A m−i is the abundance of isotope i [number of moles] that is moved to the melt; A s−i the abundance that was present in the solid before melting; F is the degree of melting (Eq. 5); and D i the isotope (and element) specific partition coefficient. Note that the right hand side starts with a multiplication with F , which is not needed when elemental fractionation is described in terms of concentrations. Since our approach We define a melting age for each particle. This is the most recent time that a particle changes its bulk composition due to melting, either by producing or receiving melt. The particles track their melting age, saving this time as one of their attributes.

Initialisation of chemistry
The bulk composition is initialised with a C value (bulk composition) of 0.6 for each 5 particle. The initialisation of trace elements is done in terms of concentrations, once attributed to the particles these concentrations are translated to abundances via the masses of the particles. At initialisation we homogeneously depleted the top 30 % of the mantle in all trace elements to 98 %. This was done as an end-member model (e.g. Armstrong, 1968) to account for the removal of heat-forming elements to the continen-10 tal crust before 3.6 Ga. We thus implicitly assume that fluxes to the continental crust are balanced by fluxes back to the mantle, although we do not explicitly model this process or its potentially heterogeneous distribution.  Table 2 and expressed algebraically as:   have decayed between the formation of the Earth and the time the calculation starts (3.6 Ga). These values are calculated using Eq. (4), with T E , the age of the Earth, and the decay constants as listed in Table 3.

Figure 3 shows snapshots of temperature distribution and bulk composition in the do-
15 main. Figure 4 shows snapshots of the temperature, melt production, bulk composition and melting age (i.e. time since melting), at the end of the calculation. The relationships between these parameters is clearly visible; in our model, mantle material melts at focussed regions of high temperature close to the surface (plumes), where the bulk composition gets altered. flow (subduction zones). Segregated material has reached the core-mantle boundary (CMB) within 500 million years. (Seen in time series of bulk compositional field and values of the domain rms-velocity which is above 1 cm yr −1 ; not shown.) The snapshots are representative of the structures as they develop over time; note both the amount of melt produced over time (Fig. 5c), and the increase in the average melting age over 5 time (Fig. 6) are steady. Figure 7 shows the radial distribution of Pb isotopes and Pb isotope ratios radially. The more basaltic rich surface layer shows up as an increase for example in 204 Pb, while just beneath we see a decrease which goes with the thin underlying residual layer. Deeper in the mantle the figure shows that in this case it is relatively well-mixed with limited variation.

Melting diagnostics
As shown in Fig. 5a our method conserves bulk composition (average value stays constant over time). The figure also shows that the surface average bulk composition becomes more basaltic than the global average, as expected. Figure 5b shows that the total number of particles present in the domain stays roughly constant at a total 15 of 1.2 billion, although particles are continuously merged and split (Fig. 5b). The fact that this method conserves bulk composition shows that the splitting and merging is not affecting the average composition. Figure 5c shows total melt production over time.
As shown, there is limited variation in the amount of melt produced as function of time.
Melt production never stops.

Pb-pseudo-isochrons vs. melting age distributions
Following Rudge (2006) we can compare our findings with his analytical solution linking pseudo-isochron ages based on Pb-isotope distributions to the distribution of melting ages in the mantle. Using; We can plot the Pb-isotope ratios carried by the particles, for example in the top layer of the model, at different times. Following Rudge (2006) we fit a geometric mean 5 regression line (also known as the reduced major axis regression line, as in Fig. 8) to these data at the different times; i.e. evaluating β of Eq. (16). Then, using Eq. (16) we can evaluate the pseudo-isochron age τ ddi , for each of these different times. Similarly, following Rudge (2006) a pseudo-isochron age can be obtained by looking at the distribution of melting ages in the mantle as follows; where q m (τ) is the probability density function of the particle melting ages; f (τ) is an arbitrary given function e.g. (e λ 238Tm − 1) 2 ; τ ddi is the pseudo-isochron age;T m is a random 15 variable which gives the distributions of the parcel ages that have undergone melting. τ s is the starting age of the model, 3.6 Ga. Rudge's theory is based on a statistical box-model and in particular assumes: 1 strong mixing. 2 heavy averaging. Figure 6 shows the build up of q m (τ) over time in our model. Since melt production is fairly constant, the histograms show a steady increase in melting age. 20 Figure 9b shows the time evolution of the obtained pseudo-isochron ages based on Pb isotopes (τ ddi from Eqs. 16 and 17; black and red curves in Fig. 9b), and particle melting ages (blue line in Fig. 9b Fig. 9b; showing Pb isotopes on the surface, our results were subsampled and used only data on Pb isotopes from particles that had undergone melting at least once (had a melting age), and were at the surface layer of the model. For the Pb isotopes sampled at the melt, our results were subsampled and used only data on Pb isotopes from melt that was produced in that time step (i.e. not taking the information from 5 a particle but from the material that moved to the surface due to melting). In this case, the pseudo-isochron has a value of 0 for the first 500 million year of the calculated time since until that time only unfractionated material is sampled.

Discussion
We have implemented tracking of radioactively decaying isotope systems in a numeri-10 cal model of mantle convective flows. Through this integration we have a tool that will allow experiments on the evolution of the Earth's mantle providing stronger constrains on both the spatial and temporal evolution. The key process is fractionation of both bulk composition and trace elements on melting. The algorithm that deals with this process (Sect. 2.6), is built on using particles to track the advection of chemical concentration 15 (bulk composition) and abundances (trace elements) with fluid flow, and move information between those particles upon melting. An advantage of this method of moving melt (i.e. via information not particles) is that we can consider any degree of melting, not just the quanta of melting that must be considered in algorithms that move particles. This allows us to consider much smaller degrees of melting. This is important for example 20 for incompatible elements in the residue.
By deliberately keeping the system simple we can compare and test the results of our numerical experiments to a quasi-analytical solution (Sect. 3, Rudge, 2006). This solution links the melting time distribution of the whole mantle to the pseudo-isochrons that can be measured in lead isotopes sampled only at the surface of the domain. the end of the modeled time the misfit is around 2 %. We note that our Eq. (17) from Rudge (2006) assumes (1) a well mixed planet, and (2) that the number of melting events that are averaged before sampling (N) is large (heavy averaging), a generalisation from Rudge et al. (2005). As regards the mixing, we note the homogeneous distribution of Pb isotopes, both radially (Fig. 7), and laterally (Fig. 4d) supports strong 5 mixing. As regards averaging we note that Rudge (2006) suggests that the dependence of the pseudo-isochron age on N is fairly weak. The good match is achieved for both sampling the surface and sampling the melt. When sampling the melt, N = 1, while particles sampled at the surface carry the signature of a collection of multiple melting events (larger N). The results presented here also support that the dependence of the pseudo-isochron age on N is weak. The good match gives us confidence in the method and therefore opens the opportunity to extract information about the interior distributions of chemical heterogeneity from surface observations. The sampling location can be important as can be seen in Fig. 9b, where we show the evolution of the pseudo-isochron based on random sampling across the surface 15 and sampling of the melt just after fractionation. Early on in the calculation, the difference between our subset sampled at the surface and the melting location is substantial mainly because the material sampled at melting locations has not gone through melting and fractionation before, whereas the surface contains fractionated material immediately after the calculation starts. On the long time scale (> BY) the difference 20 between the two isochrons is very small, again supporting the idea of a strongly mixed reservoir.
In our model setup, due to the convective pattern resulting from limited variation in viscosity, melting predominantly happens at the top of circular upwellings, i.e. plume like structures. Most melting in the terrestrial mantle (and chemical fractionation) hap- 25 pens at elongate plate boundaries, oceanic spreading ridges, and the downwelling counterpart, subduction zones. Although different in detail, the global scale distribution and evolution of both melting ages and isotope patterns seems largely unaffected. This is also shown by the limited difference observed in the pseudo-isochrons based on samples taken at random across the surface vs. those taken at the location of melt production.
Although the pseudo-isochrons we find via the melting ages and lead isotopes are consistent, in absolute value they do deviate from the one observed in lead isotopes in nature. As mentioned already the simulation case presented is intentionally simple 5 to allow a direct comparison with the analytical solution of Rudge (2006), as a result it is not Earth-like in every respect. In particular, the vigour of convection (mean velocity 1.5 cm yr −1 ) is much lower than Earth (current surface velocity 5 cm yr −1 RMS).
Also the model vigour is constant, while on Earth it is expected to be more vigorous in the hotter past. The combined effect is that the number of overturns in the simulated 10 case will be many times less than for Earth (Huang and Davies, 2007a). Therefore the number of passages through melting zones will also be much lower in this simulation than Earth. Since this model case is neither Earth-like in its vigour nor its melting a difference between pseudo-isochrons ages is not a surprise. Since more melting would remove more of the older heterogeneities and therefore reduce the pseudo-chron age 15 it might be expected therefore that more realistic models will have the potential to reconcile these differences. Future work is planned to investigate this. Future implementations will be extended with routines to allow trace elements to move to/from both continent and atmosphere reservoirs (for noble gasses), and extensions on how chemical structures affect the flow field. The good comparison to analyt-20 ical theory presented in this work, gives confidence that the current implementation is a good basis from which to include more complex and Earth-like processes into future numerical experiments. By doing so we can shift the focus from comparing numerical experiments to analytical solutions, to comparing them to observations. 25 We presented and tested a new implementation for tracking bulk chemistry and trace element abundance in a global mantle convection code that includes melting. A notable feature of the melting routine is that we transport information between particles, rather than move particles. We showed that our implementation is robust in the sense that (1) it conserves composition, (2) it conserves trace element abundance, (3) it matches the Rudge (2006) quasi-analytical solution for the prediction of isochron ages based on the distribution of melting age and pseudo-isochron ages based on lead isotopes at Introduction  Comput. Phys., 9, 207-215, 1995. 9557 Christensen, U. R. and Hofmann, A. W.: Segregation of subducted oceanic crust in the convecting mantle, J. Geophys. Res.-Sol. Ea., 99, 19867-19884, 1994. 9555, 9556 Davies, D. R., Davies, J. H., Hassan, O., Morgan, K., and Nithiarasu, P.: Investigations 5 into the applicability of adaptive finite element methods to two-dimensional infinite Prandtl number thermal and thermochemical convection, Geochem. Geophy. Geosy., 8, Q05010, doi:10.1029/2006GC001470, 2007 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Nakagawa, T., Tackley, P. J., Deschamps, F., and Connolly, J. A.: Incorporating self-consistently calculated mineral physics into thermochemical mantle convection simulations in a 3-D spherical shell and its influence on seismic anomalies in Earth's mantle, Geochem. Geophy. Geosy., 10, Q03004, doi:10.1029/2008GC002280, 2009 The influence of MORB and 5 harzburgite composition on thermo-chemical mantle convection in a 3-D spherical shell with self-consistently calculated mineral physics, Earth Planet. Sc. Lett., 296, 403-412, 2010. 9556 Oldham, D. and Davies, J. Figure 4. Snapshots taken at the end of the calculation of (a) temperature 50 km below the surface; (b) time since melting (in billion years) 50 km below the surface; (c) basalt fraction at the surface; (d) 206 Pb/ 204 Pb molar ratio at a depth of 1300 km. Note that melting in our model happens at the top of regions of central upwellings (plumes). Also; as the lead isotope figure (d) shows, the mid mantle seems to be fairly homogenous on the scale modelled here. For the melting time (b), basalt fraction (c) and lead isotope ratio (d) the values were linearly interpolated from the particles to the grid before plotting. The height of the bins shows the summed particle mass in that bin. (right) Pseudo-isochron ages as determined via the Pb isotopes sampled at the surface (red), the melting ages distribution (blue), and Pb isotopes sampled at the melt (black). The horizontal axis shows time since the start of the calculation, and the vertical axis the pseudo-isochron age.