A coupled geochemical-geodynamical approach for predicting mantle melting in space and time

13 Geodynamical simulations underpin our understanding of upper-mantle processes, but 14 their predictions require validation against observational data. Widely used geophysi15 cal datasets provide limited constraints on dynamical processes into the geological past, 16 whereas under-exploited geochemical observations from volcanic lavas at Earth’s surface 17 constitute a valuable record of mantle processes back in time. Here, we describe a new 18 peridotite-melting parameterization, BDD21, that can predict the incompatible-element 19 concentrations of melts within geodynamical simulations, thereby providing a means to 20 validate these simulations against geochemical datasets. Here, BDD21’s functionality is 21 illustrated using the Fluidity computational modelling framework, although it is designed 22 so that it can be integrated with other geodynamical software. To validate our melting 23 parameterization and coupled geochemical-geodynamical approach, we develop 2-D single24 phase flow simulations of melting associated with passive upwelling beneath mid-oceanic 25 ridges and edge-driven convection adjacent to lithospheric steps. We find that melt vol26 umes and compositions calculated for mid-oceanic ridges at a range of mantle temper27 atures and plate-spreading rates closely match those observed at present-day ridges. Our 28 lithospheric-step simulations predict spatial and temporal melting trends that are con29 sistent with those recorded at intra-plate volcanic provinces in similar geologic settings. 30 Taken together, these results suggest that our coupled geochemical-geodynamical approach 31 can accurately predict a suite of present-day geochemical observations. Since our results 32 are sensitive to small changes in upper-mantle thermal and compositional structure, this 33 novel approach provides a means to improve our understanding of the mantle’s thermo34 chemical structure and flow regime into the geological past. 35 Plain Language Summary 36 Earth’s mantle is a ∼ 3000 km thick layer of hot rock that lies between Earth’s 37 crust and core. Its slow, creeping, convective motion over billions of years has been in38 tegral to Earth’s thermal, chemical, tectonic and geological evolution. However, an in39 ability to reproduce observational constraints derived from the composition of volcanic 40 lavas at Earth’s surface limits our capacity to validate models of mantle convection back 41 in time. Here, we present a new framework that can predict the volume and composi42 tion of melts generated within the mantle. These predictions compare favourably with 43 those recorded by igneous rocks at Earth’s surface in two geologic settings: mid-oceanic 44 ridges, where plates move apart to drive decompression melting, and lithospheric steps, 45 where instabilities associated with changes in the thickness of Earth’s rigid outermost 46 shell generate volcanism far from plate boundaries. The approach and tools presented 47 here will allow scientists to better understand the mantle’s past structure, dynamics, evo48 lution and impact on Earth’s surface. 49

during melting, melt productivity as the mantle decompresses is defined as 152 where P , T , X, S, ∆S, C P represent pressure, temperature, melt fraction, entropy, en-153 tropy of fusion and specific heat capacity, respectively; the coefficient of thermal expan-154 sion and density are denoted by α and ρ, respectively; and subscripts s and f refer to 155 the solid and fluid phases, respectively (McKenzie, 1984). This equation describes adi-156 abatic decompression melting, and changes in temperature are tracked during melting 159 simultaneously with Equation 1. Here, we adapt these equations to allow each particle 160 to decompress along its respective temperature gradient calculated within each Fluidity 161 simulation. We assume that 163 where the derivative super-script signals that this temperature gradient is sourced from 198 where˙ II is the second invariant of the strain-rate tensor, n is the stress exponent that 199 depends on the deformation mechanism, E and V are activation energy and volume, A 200 is a pre-factor, P = ρgz is the lithostatic pressure, R is the gas constant and T is tem- 201 perature. Note that the rheological parameters, detailed in Table A1,  tal core, the oceanic lithosphere is treated as a thermal-boundary layer with a temper-232 ature profile that adheres to a half-space cooling model. 233 We define viscosity (µ) through a diffusion creep rheology law 234 µ = A × exp E * + ρ 0 gz V * R(T + ψz) .
(8) 235 Viscosity is temperature-and depth-dependent. At 660 km depth, µ reaches a value of Our mid-oceanic-ridge simulations are analysed once they achieve steady state and 240 particles are present throughout the melt region. In the lithospheric step case, we are 241 interested in the spatio-temporal melting trends as these simulations do not approach 242 a steady-state melting regime. Therefore, our lithospheric step simulations are initial-243 ized with particles distributed across the entire computational domain. Particles initial-244 ized at supra-solidus conditions are assigned initial melt fractions assuming that they 245 underwent instantaneous melting at that depth.

246
-7-manuscript submitted to Geochemistry, Geophysics, Geosystems   272 where X is the melt fraction, and whereD andP are the bulk distribution coefficients 273 for a given element within the solid assemblage and the melting assemblage, respectively 274 (Shaw, 1970).

275
To calculate c l and c s for a given particle, these equations are numerically integrated 276 using the LSODA algorithm (Hindmarsh, 1983;Petzold, 1983) from the solidus, where 277 X = 0, to the particle's present location, where X = X . The path between these lim-278 its for each particle, X(P, T ), is dictated by the flow field of our geodynamical simula-279 tions (Section 2).

280
Together,D andP determine how readily an incompatible element will partition 281 into the melt.D andP are the sum of the distribution coefficients for each mineral phase, 282 D mnl , scaled by their modal abundance in the solid and liquid phase, respectively. There-283 fore, i.e., if all mineral phases sum to 1 − X; Figure 2). Here, to parameterize F mnl (P, X) 314 for olivine, clinopyroxene, garnet and spinel, we perform linear regressions through these 315 data at constant pressure:

317
The empirically determined constants, a mnl and b mnl , from each linear regression at each 318 pressure are then combined to parameterize F mnl (P, X) via a second-order polynomial 319 regression:

321
At depths between ∼ 60-90 km, the stable aluminous phase changes from spinel 322 to garnet through the following reaction: When performing polynomial regressions from 1-7 GPa, the proportion of garnet present 324 at the solidus (i.e., F gnt when X = 0), approaches zero at ∼ 2.7 GPa. The pressure 325 range over which the spinel-garnet transition occurs is controversial: thermodynamical 326 models typically place this transition at lower pressures than it occurs at within peri- pecially heavy rare earth elements, is strongly dependent on the presence of garnet, we 329 allow the depth of the spinel-garnet-transition zone to be a user-defined variable within 330 our parameterization. We therefore calculate the modal mineralogy constants (a 0−2 and 331 b 0−2 ) for spinel-and garnet-bearing peridotite separately using experiments between 1-332 1.5 GPa and 4-7 GPa, respectively. All experiments from 1-7 GPa are used to calibrate 333 F cpx (P, X) since clinopyroxene is not precipitated or consumed at this phase boundary. 334 We assume that mantle modal mineralogy varies linearly between the mineral assemblages 335 estimated to be present at the top and bottom of the spinel-garnet-transition zone.

336
The proportion of orthopyroxene present in the mantle (F opx ) does not vary lin-337 early as a function of melt fraction, since it strongly depends on whether clinopyroxene 338 is a stable phase ( Figure 2). For simplicity, we assume that F opx accounts for the total 339 missing fraction when all other mineral contributions are combined: tion of pressure, temperature and composition using a lattice strain equation .

353
The compatibility of an element within a mineral is principally controlled by the valency sites within mineral lattices to expand and contract . As such, to use this  is 200× greater than that of Ce . The concentration of water in depleted 387 mantle is set to 100 ppm (Salters & Stracke, 2004). Before applying BDD21 to 2-D geodynamical simulations, we first describe how man-390 tle mineralogy and melt composition respond to isentropic melting in a simple 1-D case.

391
A particle of primitive mantle with T p = 1325 • C ascends vertically from 3.5 GPa to 392 the surface. As this particle rises, it intersects the solidus at ∼ 100 km depth and be-    tated, respectively. Therefore, heavy rare earth elements become increasingly incompat-428 ible in the solid phase and preferentially enter the melt phase until they are exhausted.

429
Since the integrated melt composition (C l ) is a weighted sum of all instantaneous (c l ) 430 generated thus-far, the intersection between the isopleth that defines v z = R s /3 and the anhydrous peri-455 dotite solidus. We assume that all melts generated beyond these limits freeze within the 456 lithosphere. Thus, the thickness of crust produced at the ridge axis (t c ) is calculated by

458
where M and A correspond to the the melting rate and melting region over which melt 459 reaches the ridge axis, respectively(Forsyth, 1993). We assume that the density of oceanic   in plagioclase and magma chambers are constantly replenished by primitive melts.

537
One way to avoid the pitfalls of conducting a fractional crystallisation correction 538 is to employ incompatible element ratios. The geochemical ratio 10 4 (Sm+Gd)/Ti is thought 539 to be largely unaffected by fractional crystallisation and is observed to be almost con-  face, which increases effective viscosity, further decreasing upwelling rates and, thus, caus-569 ing melting to cease at greater depths.

570
In our simulations, these effects are particularly pronounced at R s < 2.1 cm yr −1 .

571
When R s ≥ 2.1 cm yr −1 , conductive cooling is reduced and maximum melt fraction 572 and melting rates, in addition to crustal thickness, increase significantly (Figures 6a-c).

573
However, as R s increases beyond 2.1 cm yr −1 we observe a narrowing in the melt focus-

574
ing distance, x f , and a gradual reduction in crustal thickness estimates. This pattern

641
We note that in our simulations we do not consider the relative densities of fertile 642 and depleted mantle. Fertile mantle is denser than depleted mantle (Jordan, 1978 To ensure that our approach is also valid in these low-melt-productivity environments, 654 we construct and analyse a series of coupled geochemical-geodynamical simulations that 655 are applicable for such intra-plate settings.  Both geochemical modeling and mantle temperature estimates derived from seismic to-708 mographic models suggest that the lithosphere beneath these volcanic provinces is ∼ 50- However, the use of geochemical constraints remains in its infancy and will doubtless be-  We use BDD21 to predict primary magmatic compositions. It is important to re-764 member that when we compare these predictions to observations at mid-oceanic ridges, 765 we are comparing to basaltic glasses that represent the vestigial melts of complex mag- be treated with caution, it is encouraging that our mid-oceanic ridge simulations repli-774 cate the broad trends found in two very different geochemical systems -Na 2 O and 10 4 (Sm+Gd)/Ti.

775
In future, given that each incompatible element behaves differently during fractional crys-776 tallisation, it may be informative to use our predictive primary melt compositions as a 777 starting point to investigate these fractional crystallisation processes. 778 We compare our simulations to mid-oceanic ridges that have potential tempera-779 tures that are estimated to be between 1250-1650 • C using techniques that convert from  (Figure 9d).

815
The increasing disparity between predicted and observed crustal thicknesses as po- In Section 4.2, we predict that crustal thicknesses decrease as a function of spread-841 ing rate. To obtain these values, we make use of the observation of  842 that the x f can be approximated through a relationship between upwelling rate, spread-843 ing rate and depth of the anhydrous solidus. Thereby, our predicted decrease occurs be-    where R cpx = R 1 + R 2 P . (2)

29
Our geochemical melting model is constrained using experiments conducted on MM3 and 30 KR4003 peridotites (Baker & Stolper, 1994;Walter, 1998;Falloon et al., 1999). X cpx−out  Figure S1). Here, we exploit the M cpx and X cpx−out values recorded within these experiments to approximate M cpx , R 1 and R 2 as 0.18, 0.94 and −0.1, respectively. To provide a good fit for melt fraction as a function of temperature, two additional constants, 36 β 2 and B 1 , are also updated to 1.2 and 1520 • C, respectively ( Figure S2).

37
The decompression melting model parameterised by Katz et al. (2003) assumes that no 38 heat is lost during melting. However, it is necessary to modify these equations in this 39 case since within our geodynamic model we allow heat diffusion to occur. We replace the 40 adiabatic gradient term with the actual temperature gradient that is experienced by each 41 particle in the geodynamic model using Equations 3-5 in the main text. This replacement 42 requires the following thermodynamic assumption:

44
where P , T and S are pressure, temperature and entropy, respectively; and α, ρ and C P 45 denote thermal expansivity, density and heat capacity, respectively. coefficients can either be assumed to be constant using the values listed in Table S1, or 56 they can vary as a function of pressure (P ), temperature (T ) and mineral chemistry. D spl 57 is assumed to be constant as a function of P and X for all elements. In the mantle, r 3+ 0(ol) = 0.72 × 10 −10 .
Where χ Al min is the modal proportion of Al in a mineral, which in this case is olivine (ol), 79 and can be calculated by respectively. All other constants are listed in Table S2.