Melting conditions and mantle source composition from probabilis-2 tic joint inversion of major and rare earth element concentrations 3

17 The chemical composition of erupted basalts provides a record of the 18 thermo-chemical state of their source region and the melting conditions that 19 lead to their formation. Here we present the first probabilistic inversion 20 framework capable of inverting both trace and major element data of mafic 21 volcanic rocks to constrain mantle potential temperature, depth of melting, 22 and major and trace element source composition. The inversion strategy is 23 based on the combination of i) a two-phase multi-component reactive trans24 port model, ii) a thermodynamic solver for the evolution of major elements 25 and mineral/liquid phases, (iii) a disequilibrium model of trace element par26 titioning and iv) an adaptive Markov chain Monte Carlo algorithm. The 27 mechanical and chemical evolution of melt and solid residue are therefore 28 modelled in an internallyand thermodynamically-consistent manner. 29 We illustrate the inversion approach and its sensitivity to relevant model 30 parameters with a series of numerical experiments with increasing level of 31 Email address: oliveira.bravo.b@gmail.com (B. Oliveira) Preprint submitted to Geochimica et Cosmochimica Acta August 9, 2021 complexity. We show the benefits and limitations of using major and trace 32 element compositions separately before demonstrating the advantages of a 33 joint inversion. We show that such joint inversion has great sensitivity to 34 mantle temperature, pressure range of melting and composition of the source, 35 even when realistic uncertainties are assigned to both data and predictions. 36 We further test the reliability of the approach on a real dataset from a well37 characterised region: the Rio Grande Rift in western North America. We 38 obtain estimates of mantle potential temperature (∼ 1340 C), lithospheric 39 thickness (∼ 60 km) and source composition that are in excellent agreement 40 with numerous independent geochemical and geophysical estimates. In par41 ticular, this study suggests that the basalts in this region originated from a 42 moderately hot upwelling and include the contribution from a slightly de43 pleted source that experienced a small degree of melt or fluid metasomatism. 44 This component is likely associated with partial melting of the lower portions 45 of the lithosphere. The flexibility of both the melting model and inversion 46 scheme developed here makes the approach widely applicable to assessing 47 the thermo-chemical structure and evolution of the lithosphere-asthenosphere 48 system and paves the way for truly joint geochemical-geophysical inversions. 49


Introduction
Most volcanism observed on Earth is the result of partial melting in the coefficients. Lastly, the role of the mixing function is to describe the spatial distribution of both solid and melt within the melting regime, including how 154 instantaneous melts within the melting zone are aggregated and extracted. 155 Despite often being constructed separately, previous works demonstrate 156 that these three functions are intimately related to one another (e.g. Langisfy mass and energy balance constraints, as they are both controlled by the  Following previous studies such as that of Asimow (2002), we propose an 166 internally-consistent thermodynamic framework for decompression melting 167 in a one-dimensional upwelling mantle column (Fig. 1). We assume an isen-168 tropic fractional fusion model, which in practice is modelled as a sequence 169 of infinitesimal isentropic melting steps along the melting column. At each 170 step, the newly produced melt (Fig. 1, blue lines) is extracted and isolated 171 into an adjacent melt reservoir (Fig. 1, red lines). In order to capture rela-172 tive movement between solid and melt phases during the melting process, we 173 couple the internally-consistent thermodynamic database and formalism for 174 mantle melting of Jennings and Holland (2015) with a two-phase transport 175 model. As such, we integrate the melting, chemistry and mixing functions 176 into a single thermo-chemical-dynamical framework. We use components of 177 the software Perple X (Connolly, 2009) to obtain degree of melting as a func-178 tion of pressure, temperature and bulk composition (i.e. melting function) 179 via Gibbs free-energy minimization. This also allows us to obtain the ma- phases (e.g. Smith and Asimow, 2005). 189 Finally, we use the "residual mantle column" approach  muir, 1992) to evaluate the geochemical effects of mantle flow on the aver-  and solid phases in mantle peridotites (Jennings and Holland, 2015). Given 215 pressure, temperature and solid composition at each incremental depth z i , we 216 minimise Gibbs free-energy using Perple X to obtain solid and melt equilib-217 rium compositions, thermodynamic properties and the equilibrium mineral 218 assemblage. If melting is isentropic, the associated temperature change be-219 tween two consecutive decompression steps is not a free variable and is not 220 known a priori. This problem can be circumvented by either estimating this 221 temperature change using thermodynamic constraints or by minimizing a 222 thermodynamic potential such as enthalpy to guarantee isentropic conditions 223 at each incremental depth (Asimow et al., 1997;Morgan, 2001;Brown and 224 Lesher, 2016). Since we rely on a Gibbs free-energy minimization algorithm, 225 we opt to correct the temperature iteratively at each pressure increment until 226 the difference between the entropy of the system before and after each mini-227 mization falls below a certain tolerance. Once the temperature of the system corresponds to that of an isentropic decompression melting step, the melt 229 fraction and its entropy are removed from the system (which does not affect 230 the specific entropy of the solid residue, nor its composition). We record solid 231 and melt compositions and thermodynamic properties, and proceed with the 232 next isentropic decompression melting step. 233 At every decompression step, we calculate F (z i ) by adding the extracted 234 melt mass fraction to the degree of melting obtained in the previous decom-235 pression melting step F (z i−1 ) as where f refers to the generated melt mass fraction in equilibrium. Similarly, 237 for each mineral phase we define where f j is the mass fraction of mineral phase j. F j describes the mass where φ is volume fraction (subscripts inst, iso, s and j refer to instantaneous 256 liquid, isolated liquid, solid and mineral phases, respectively), ρ is density, z is 257 depth, and w and W are the liquid and solid velocities. Γ corresponds to the 258 mass-transfer from solid to instantaneous melts, and S is the mass-transfer 259 from the instantaneous melts to the isolated melt. Since the total solid mass 260 corresponds to that of all individual mineral phases, we have The total mass-transfer rate is Γ = j Γ j and Γ b = j Γ b j .

262
Here, instantaneous melts are kept in chemical isolation upon formation, 263 which implies that φ inst = 0, S = Γ, and that the total melt fraction is equal 264 to the melt in chemical isolation, φ l = φ iso . With these considerations we 265 simplify Eqs. 3 and 4 and obtain The mass-transfer term of chemical component b, Γ b j , considers both phase 267 changes and diffusion of trace elements from mineral phases to melt as where "0" refers to properties evaluated at the onset of melting, and h is 280 the depth at which melting stops. In steady-state, the mass-transfer rate 281 Γ is related to the degree of melting F obtained from the melting function information we rewrite Eqs. 9-10 and obtain Equations 11 and 12 show that the liquid and solid fluxes (i.e. products 285 of volume fraction, density and velocity) are balanced by melt production.

286
Therefore, once the liquid volume fraction is known, solid and liquid velocities 287 can be determined. Equation 11 is used to obtain the liquid fraction. 288 We relate φ l , w and W with a Darcy-type functional relationship where µ l is the liquid viscosity, k 0 is the permeability constant, n the perme-290 ability exponent, and g is the acceleration of gravity. Equation 13 shows that 291 melt segregation is governed by a balance between the differential buoyancy 292 of the melt and the resistance to flow of melt through the matrix (Darcy 293 resistance). For simplicity, we ignore the resistance of the solid phase to 294 deformation which is of second-order only (no compaction term in Eq. 13).

295
Since we assume a one-dimensional upwelling model, further assumptions 297 are needed to capture the effects of lateral mantle flow on the modelled 298 compositions. Given its applicability to a wide range of tectonic settings, we 299 employ the "residual mantle column" approach to approximate the pooling 300 of melts over the whole melting zone (Plank and Langmuir, 1992).

301
Given our model assumptions, the average major and trace element con-302 centrations of the pooled magma can be calculated as where where d obs is the vector of measured values, n is the dimension of d obs and 430 C D is the covariance matrix. In this expression we implicitly assume that beyond the scope of this paper. Here we simply assume that one can have 439 access to C preparation and C analytical for all data types (e.g. from analytical 440 measurements) and that C sampling can be at least estimated. Errors associated with deterministic numerical models are not random, 443 but rather they are controlled by systematic biases (e.g. underestimation 444 due to simplifying assumptions) and human errors (e.g. 'bugs' in the code).

445
A common example of a systematic error is that related to the coarseness of where e(m) is the sum-of-squares function, which acts as a misfit function 478 between observed and predicted data. It is computed as where d pre represents the prediction from the forward problem for model  boundary layer. We therefore assume that the process of isentropic 501 decompression melting stops at Z top .

502
• The major and trace element composition of a peridotitic mantle source.

503
The mantle trace element content, by which we mean REEs in partic-504 ular, is not always correlated with that of major elements (e.g. an

505
"enriched" signature in traces can co-exist with a "depleted" signature 506 in majors). We therefore treat these two different groups of elements 507 as independent. In contrast, correlations exist within each group. This   Table 2.  Table 1), we created the "incorrect" source by simply     tion coefficients of all elements in clinopyroxene and garnet within a generous, yet realistic, uncertainty of 25% (as 1 STD; thus a total allowed variation 697 of ∼ 100%). As before, we do not vary the partition coefficients in olivine, 698 orthopyroxene and spinel, as their effect is only second-order compared to 699 that of clinopyroxene and garnet.

700
The results of these simulations are summarised in Table 3. In general,  Table 2.

744
For the forward problem we consider partition coefficients from Wood

886
In addition to the improvements discussed above, the present model is    Source compositions are given in Table 1, and partition coefficients are from Oliveira et al.
where φ l is the melt volume fraction, ρ is density, v is velocity, and Γ is the 970 rate of mass exchange between phases (or simply, the melting rate).

971
We consider that the solid phase is comprised of several mineral grains j, velocity. Note also that Γ j refers to the mass lost/gained by each mineral 975 specie j; therefore we have that j Γ j = Γ.

976
For a one-dimensional (depth-dependent only) steady-state case, Eqs.
where W and w are the solid and liquid velocities, respectively, and h is the Similarly, for each mineral grain in the solid residue, we have and mass conservation imposes that j F j = F .

987
Combining Eqs. A.4, A.5 and A.6 with A.7 and A.8, we obtain the final 988 mass-balance equations as functions of F , where d is a symmetric, rheology-dependent, interaction coefficient and mod-997 elled as (Bercovici et al., 2001) where k(φ l ) = k 0 (φ l ) n is the Kozeny-Carman-type permeability law relating 999 permeability and porosity k 0 (φ l ); n is a constant exponent. Here we take and since the solid phase is comprised of several mineral grains j, each is 1006 subjected to a mass conservation equation for its chemical composition, c b j , Because the chemical-mass lost by the solid, Γ b , is the aggregated contribution of each of its constituents, Γ b j , chemical-mass Details on how to model Γ b j are provided both liquid reservoirs, and the generalization of Eq. A.15 leads to, where φ eq and φ iso are the volume liquid fractions allowed to equilibrate and 1025 that in chemical isolation, respectively. Because of mass-balance constraints 1026 φ eq +φ iso = φ l . S is the rate of mass exchange between both liquid reservoirs, where K b j is the usual partition coefficient.

1050
Alternatively, one could fix the solid composition as melting proceeds (i.e. is plenty of scope for further exploration of these laws.

1055
For J b j , we use linear kinetics to approximate mineral-melt finite exchange 1056 that arises from diffusion in minerals and/or dissolution-precipitation, where R b j is the exchange rate constant for the chemical component of interest 1058 between mineral j and the melt. As mentioned above, we consider that solid Chemical-mass conservation in residual solid, s: Chemical-mass conservation in isolated melt, iso: where the instantaneous melt composition c b inst reads Chemical-mass conservation in mineral grains, j: ties. This is done using a standard finite difference approach.

1127
• Last, we compute the averaged major and trace element concentration 1128 of the pooled magma using the mixing function of Eq. 14 for each 1129 melting column, and then interpolating the result for the given set of