Effective buoyancy ratio : a new parameter for characterizing thermochemical mixing in the Earth ’ s mantle

Numerical modeling has been carried out in a 2-D cylindrical shell domain to quantify the evolution of a primordial dense layer around the core–mantle boundary. Effective buoyancy ratio, Beff was introduced to characterize the evolution of the two-layer thermo-chemical convection in the Earth’s mantle. Beff decreases with time due to (1) warming of the compositionally dense layer, (2) cooling of the overlying mantle, (3) eroding of the dense layer through thermal convection in the overlying mantle and (4) diluting of the dense layer through inner convection. When Beff reaches the instability point, Beff = 1, effective thermo-chemical convection starts, and the mantle will be mixed (Beff = 0) over a short time period. A parabolic relationship was revealed between the initial density difference of the layers and the mixing time. Morphology of large low-shear-velocity provinces and results from seismic tomography and normal mode data suggest a value of Beff ≥ 1 for the mantle.


Introduction
The most prominent feature of the lowermost part of the Earth's mantle is the two seismically slow domains beneath the Pacific and Africa (e.g., Dziewonski et al., 1993;Garnero et al., 2007a).The nearly antipodal large low-shear-velocity provinces (LLSVPs) are characterized by −2 to −4 % shear wave and −1 to −2 % primary wave anomaly, several thousand kilometers lateral extent and 800-1000 km elevation from the core-mantle boundary (CMB) (Mégnin and Romanowicz, 2000;Masters et al., 2000;Lay, 2005;Zhao, 2009).The margins of the anomalies, where the lateral shear wave velocity gradients are the most pronounced, have sharp sides (Ni et al., 2002;Wang and Wen, 2004;Ford et al., 2006;Garnero and McNamara, 2008) and correlate with hot spot volcanism (Thorne et al., 2004;Torsvik et al., 2010).The existence and the morphology of LLSVPs cannot be satisfactorily explained by the variation in temperature, mineralogical phases or melts.Compositionally dense and thus stable material accumulated above the CMB is necessary in a consistent mantle model (Trampert et al., 2004;Ishi and Tromp, 2004;Garnero et al., 2007b;Bull et al., 2009).
A compositionally dense layer around the core is expected to hinder the mantle convection through reducing the heat transport from the Earth's core (Nakagawa and Tackley, 2004).Thus, a chemically dense layer at the base of the mantle has a stabilizing role (Sleep, 1988;Deschamps and Tackley, 2009).On the other hand, the heat coming from the core is trapped in the dense layer leading to a hot and unstable bottom thermal boundary layer.The dominant process of the two opposite effects can be predicted by the buoyancy ratio (Davaille et al., 2002), which is the ratio of the stabilizing chemical density difference and the destabilizing thermal density difference.β denotes the relative chemical density difference between the layers, α is the thermal expansion coefficient and T m is the temperature difference across the mantle.When B is larger than 1, the dense layer is thought to be stable, but in the case of B < 1, the density decrease through thermal expansion is strong enough to break up and mix it with the overlying mantle through thermo-chemical convection (TCC).
Published by Copernicus Publications on behalf of the European Geosciences Union.
As early as the 1980s, pioneer numerical simulations were made to investigate the effect of the compositionally dense lower layer on the mantle dynamics (Christensen and Yuen, 1984;Hansen and Yuen, 1988).Laboratory experiments and numerical models of mantle convection have shown that a chemically dense primordial layer can survive for the age of the Earth if B is large enough (e.g., Davaille et al., 2002;Jellinek and Manga, 2002;Lin and Van Keken, 2006).Depending on the density contrast and the initial thickness of the dense layer thermo-chemical domes/piles are formed in these models which morphologically resemble the seismological LLSVPs (Trampert et al., 2004;Bull et al., 2009).Deschamps andTackley (2008, 2009) investigated systematically the influence of some important parameters (depth-, temperature-and concentration-dependent viscosity, internal heating, chemical density contrast, mineralogical phase change at 660 km) on the evolution of the initial dense layer and compared the power spectra of density and thermal anomalies obtained from seismic tomography and numerical models.They mapped the parameter space of the thermochemical convection and suggested the essential ingredients for a successful mantle convection model.
In these thermo-chemical models, B is time-independent during the simulations.However, the primordial dense layer might change greatly due to the heat from the core and possibly from the decay of enriched radioactive elements, the surface erosion of dense material through convection occurring in the overlying mantle, internal convection within the dense layer and termination of subducted slabs at CMB (Nakagawa and Tackley, 2004;Lay, 2005;McNamara and Zhong, 2005;Lay et al., 2006;Garnero et al., 2007a).In this paper we present the results of numerical model calculations made with different values of B including values larger than 1.We studied the evolution of the convection, and we suggest the introduction of the time-dependent effective buoyancy ratio which characterizes better the dynamics of the TCC.

Model description
Boussinesq approximation of the equation system governing the thermo-chemical convection was applied (Chandrasekhar, 1961;Hansen and Yuen, 1988;Čížková and Matyska, 2004).The dimensional equations expressing the conservation of mass, momentum as well as the heat and the mass transport are where the unknown variables are the density, the pressure, the flow velocity, the temperature of the fluid and the concentration of the dense material, ρ, p, u i , T and c, respectively.In a two-dimensional model domain there are five equations to determine six variables.Therefore, a simple linear relationship is given among the density, the temperature and the concentration by the equation of state, where ρ R and T S denote the reference density and the surface temperature, β is the initial relative density difference between the dense layer and the overlying mantle.Q and σ ij are the internal heat production and the deviatoric stress tensor for incompressible Newtonian fluid, respectively.The space coordinates and the time are denoted by x i and t, respectively; e i shows the direction of the gravitational acceleration, downwards.According to the Boussinesq approximation, other parameters in Eqs.(2-6) are supposed to be constant (Table 1) (Van Keken, 2001).Focusing on the processes of the disintegration and homogenization of the primordial dense layer above the core-mantle boundary, dynamic viscosity was chosen to characterize the deep mantle instead, resulting in less intense convection with a thermal Rayleigh number of about 6 × 10 6 .Finite element method was applied to solve the partial differential equation system of Eqs.(2-5) using COMSOL Multiphysics software package (Zimmerman, 2006).A field method was applied to calculate the concentration distribution of dense material.The applied numerical scheme was tested, and the comparison with the benchmark study of Van Keken et al. (1997) is presented in the Supplement.Twodimensional cylindrical shell geometry was used to approximate the shape of the Earth's mantle.Geometrical scaling was adopted from Van Keken (2001) to maintain the ratio of the CMB and the Earth's surface ( ∼ = 0.3) and not to overstate the role of the deep mantle; thus, the outer and inner radius of the mantle were 4123 km and 1238 km, respectively.The boundaries were isothermal as well as symmetrical and impermeable with respect to the velocity and the concentration.
Simulation was started from a quasi-stationary state of the temperature field obtained from a chemically homogeneous, purely thermal convection model.Concentration of dense material was set to 1 for the dense layer and 0 above; the transition was adjusted using a smoothed Heaviside function with continuous first derivative and interval thickness of 50 km.The initial thickness of the dense layer was 300 km around the core.Maximum element size was 50 km within the model domain, 30 km along the surface and 15 km along the CMB and the surface of the initial dense layer (300 km above the CMB) to ensure sharp variation in the thermal and/or chemical boundary layer.During the systematical model calculations, the mantle was taken as isoviscous without internal heating.The only parameter modified during the simulation was the initial relative density difference between the dense layer and the light overlying mantle, β; it ranged between 0-8 % (B = 0−1.33).We investigated the effect of β on the monitoring parameters; heat flux, velocity, temperature and concentration time series were calculated in the upper and the lower layer.Hence we use the lower and upper layer expression in geometrical meaning as the deepest 300 km thick part of the mantle and the overlying zone, respectively.For example, denote the concentration of the dense material in the upper and the lower layer, respectively, and A 0 and A 1 the area of the mantle above and beneath 300 km above the CMB.
Other volumetric parameters such as the temperature and the velocity are calculated in the same manner (Table 2).The concentration and temperature difference between the layers is calculated as Indices S, D and CMB denote the values at the surface, at the depth of 300 km above the CMB and at the CMB, respectively.Time is defined as non-dimensional diffusion time.
In addition, we compiled a model with complex rheology (depth-, temperature-and composition-dependent viscosity) and composition-dependent internal heating to test their influence on the effective buoyancy ratio.

Results
Figure 1 illustrates the influence of a basal dense layer on the heat flux, velocity, temperature and concentration time series (left) as well as on the evolution of the concentration and temperature field (right).The initial density difference was β = 6 % between the layers that results in B = 1 for the buoyancy ratio.The initial state (stage a) is given by a temperature field obtained from a purely thermal convection calculation and a compositionally dense basal layer Effective buoyancy ratio placed instantaneously above the CMB.After approximately 1 Gyr (stage b) two-layer convection is being evolved separately in the upper and the lower layers.Inner convection within the dense layer and cold downwellings in the overlying mantle deform the surface of the dense layer.At this stage, the temperature of the dense layer reaches its maximum (T 1 ), and the heat flux (q S , q CMB , q D ) decreases to a low quasi-stationary level.The erosion of the dense layer through thermal convection in the overlying mantle reduces the concentration of the dense material in the lower layer (c 1 ) and increases it in the upper one (c 0 ).The concentration variation shows a linear trend.A similar linear reduction in the volume of the dense layer was found by Zhong and Hager (2003), who studied the entrainment of the dense material by examining one stationary thermal plume.4.5 Gyr later (stage c) the dense layer disintegrates, it becomes unstable and effective thermo-chemical convection (TCC) starts.
The TCC mixes the layers quickly, the flow accelerates (v 0 , v 1 ), the heat flux increases (q S , q CMB , q D ), the dense layer cools (T 1 ), while the upper layer warms (T 0 ).The mass flux of the dense material (q DC ) starts up and the heterogeneity of the concentration (c het , normalized standard deviation of the concentration) decreases suddenly.In other words, the thermal energy of the dense layer transforms to kinetic energy in a short period of time.After 5.1 Gyr (stage d) the dense layer ceased, having been mixed in the mantle, the system reached the stable state.Time series converge to the values characterizing the pure thermal convection, concentration time series tend to the average value of 0.0538.The heat flux (q S , q CMB , q D ) and velocity (v 0 , v 1 ) time series have higher values and larger fluctuations than in the two-layer convection regime (from stage a to d) that underlines the retaining role of the  chemically dense bottom layer.Of course, the homogenization continues protractedly, and after 7.8 Gyr (stage e) the heterogeneity (c het ) decreases below 1 %.The heating of the mantle (T ) requires billions of years.
Figure 1 illustrates that although the buoyancy ratio is B = 1 -that is, the stabilizing chemical density difference and the destabilizing thermal density difference is balanced -the dense layer evolves considerably, moreover disappears during disappears after about 5 Gyr.Additional model calculations revealed that mixing of the layers occurred for both B < 1 (β < 6 %) and B > 1 (β > 6 %).Therefore, we suggest introducing the effective buoyancy ratio in order to characterize the evolution of the dense layer and the dynamics of the thermo-chemical convection.The effective buoyancy ratio, is time-dependent and includes c concentration and T temperature differences between the bottom layer (i.e., the lower 300 km of the mantle) and the overlying mantle.
Figure 2 shows the concentration and temperature differences between the layers as well as the calculated effective buoyancy ratio at different values of β.As the dense layer warms up from heat coming from the core and the overlying mantle cools down from retained heat transport due to twolayer convection, the temperature difference increases.It results in the initial rapid decrease of B eff .The concentration difference is decreased monotonically through the erosion of the dense material that later becomes the dominant process in reduction of B eff .Gonnermann et al. (2002) observed a similar decrease of the time-dependent buoyancy ratio in their laboratory experiments.However, they used the temperature drop across the interface of the two layers and the temperature difference between the bottom of the tank and the lower layer for the estimation of the chemical and the thermal density difference, respectively, to define the buoyancy ratio.When the effective buoyancy ratio reaches the value of B eff = 1, that is, the instability point of the system (stage c in Fig. 1), one-layer thermo-chemical convection (mixing) starts.Mixing results in the quick reduction of the temperature and concentration differences.When the effective buoyancy ratio reaches the value of B eff = 0 (stage d in Fig. 1), the dense layer ceases, the mantle becomes mixed.Overturns of dense material cause temporarily negative values in B eff , especially in the cases of lower initial density contrast (β).It is obvious that larger initial density contrast entails more stable layering, however the mixing occurs in each model even for B > 1.
We attribute the occurrence of the effective thermochemical convection in each model to four main physical processes: 1. Heat coming from the core warms up the dense layer reducing its density through thermal expansion.
2. The overlying mantle cools down through retained heat transport due to two-layer convection.
3. Thermal convection forming in the upper layer erodes the surface of the dense layer through viscous drag.
4. Inner convection within the dense layer intermixes light material from the overlying mantle.
Processes (1) and (2) result in the increase of the temperature difference between the layers, the processes (3) and ( 4) cause a decrease in the concentration difference.While the first two phenomena are constrained by the total temperature drop across the mantle (practically T m /2, see Fig. 2b), the latter two are not.Erosion (3) and dilution (4) gradually reduce the chemical density difference between the layers until the system reaches the instability point (B eff = 1) when mixing begins.Mixing occurs in every case, even if the time might exceed the Earth's age (B ≥ 1). Figure 3 illustrates the phenomena of the erosion and dilution of the dense layer in the concentration and the temperature fields.Black arrows denote (a) the mass flux of the dense material and (b) the velocity of the flow.We investigated how the occurrence time of the two most characteristic events (the onset and the end of the effective TCC) depends on the initial chemical density difference, β (Fig. 4a).Obviously, larger β results in a more stable, longlived dense layer and larger occurrence time.A parabolic relationship was found between the occurrence time of B eff = 1 (onset of mixing) and β.Davaille (1999a) observed a similar relationship in their laboratory experiments studying the effect of the buoyancy ratio (and other parameters) on the entrainment rate.Parabolic function fits well on data of B eff = 0 (end of mixing) too.
As shown in Fig. 2, both the erosion/dilution phase (leading up to stage c) and the effective TCC phase (between stage c and d) can be characterized by a linear decrease in c.The effective buoyancy ratio displays a similar feature apart from its initial phase, which is due to the transient heating of the dense layer and the cooling of the overlying mantle (from stage a to b). Figure 4b illustrates the slope of the linear curves fitted on c and B eff time series during the erosion/dilution phase.It is established that larger initial density contrast (β or B) entails more stable layering owing to the less effective erosion/dilution process.Figure 4b presents a power function relationship between the slopes of time series and β ( c ∼ β −1.91 or B eff ∼ β −2.35 ).Both the parabolic relationship in Fig. 4a and the power function relationship in Fig. 4b support the idea that mixing of the layers occurs for arbitrary density contrast.It is worth noting that the effective TCC phase demonstrates also a linear decrease in c and B eff , but with steeper slope (Fig. 2).The slope of the linear curves fitted on the time series shows a slight decrease as β increases (not shown).
Focusing on the erosion phase (from stage b to c) the linear decrease of the concentration difference between the layers, or in other words, the constant entrainment rate has been shown earlier (e.g., Davaille, 1999b;Gonnermann et al., 2002;Zhong and Hager, 2003).The entrainment of the dense material through a thin tendril can be calculated using the balance between the viscous force dragging the dense material upward and the buoyancy force drawing it back (Jellinek and Manga, 2002), where v p , l and ρ denote the plume velocity, the thickness of the tendril and the density difference between the layers, respectively.The tendril thickness from Eq. ( 8) is The concentration of the dense material in the lower layer decreases through the dense material entrained by hot plumes, where N p is the number of plumes and t is the elapsed time.Thus, the entrainment rate which is the concentration variation in the lower layer is At the relative density difference between the layers of β = 0.06 the plume velocity varies from 3×10 −11 m s −1 to 10 −10 m s −1 in the numerical model during the erosion phase just above the lower layer resulting in a tendril thickness of 10-20 km (Table 1).Applying Eq. ( 13) the concentration decrease caused by hot plumes of about N p = 7 (Fig. 1b and c) gives a value of 0.85-5.2×10−18 s −1 .On the other hand, in Fig. 1 the concentration decrease of the lower layer during the erosion phase has a slope of 3.3×10 −18 s −1 , or in nondimensional form, 27.2.These conclusions were drawn from a simple isoviscous model.However, the TCC leading to the dissolution of the dense layer strongly depends on the viscosity.Therefore, a more complex model including depth-temperatureand composition-dependent viscosity and compositiondependent internal heating was calculated in order to investigate the dynamics of the TCC and the variation of the effective buoyancy ratio.Parameters controlling the viscosity and the internal heating were assigned based on the results of Deschamps andTackley (2008, 2009).An Arrhenius-type law determined the depth-, and temperature-dependence of the viscosity, which increased 1 order of magnitude from the surface to the CMB and decreased 6 orders of magnitude with the temperature.A viscosity jump with a factor of 30 was superimposed at the depth of 660 km, reflecting the effect of mineralogical phase change on the viscosity.The viscosity of the dense material (c = 1) is half of that of the light material (c = 0) with a linear transition.Internal heating was adjusted to produce 65 mW m −2 average heat flux on the surface, but the heat production of the dense material was increased by a factor of 10 due to the higher abundance of radioactive elements.The initial compositional density contrast between the layers was β = 6 %, which corresponds with the model presented in Fig. 1.The simulation started from a quasi-stationary state of the temperature field obtained from a chemically homogeneous, purely thermal convection model with depth-and temperature-dependent viscosity and homogeneous internal heating.
Figure 5 illustrates the pattern of the TCC for the complex model at 3.5 Gyr after the inset of the dense layer when the effective buoyancy ratio is approximately 1.13.Over a period of 3.5 Gyr the dense layer disintegrated and two hot, compositionally dense, nearly antipodal piles formed with sharp sides.Due to the concentration-dependent internal heating, the temperature within piles exceeds the CMB temperature; thus, the viscosity decreases considerably.The concentration and velocity field attest to a sluggish internal convection forming within the piles.A stagnant lid regime evolved owing to the strongly temperature-dependent viscosity (Solomatov, 1995), which does not participate in the convection.Beneath the stagnant lid, vivid small-scale convection occurs in the upper mantle (Kuslits et al., 2014).Due to the lack of the endothermic phase transition, advective mass and heat transport exists between the upper and lower mantle.
Figure 2 displays the variation of the concentration and temperature differences between the layers and the effective buoyancy ratio for the "mantle-like" model (mm_6 %).As a consequence of the stagnant lid regime, T decreased compared to the isoviscous case, but the character of the curve remained similar.The rate of the decrease in c through erosion and dilution processes became steeper owing to the reduced viscosity of the hot, dense thermo-chemical layer.As a result, the effective buoyancy ratio has a similar nature with a steeper erosion/dilution phase and less steep mixing phase.In summary, the stability of the dense layer in the complex model with varying viscosity and internal heating was reduced compared to isoviscous model by about 20 %, but the physical processes acting in the two models were the same.

Discussion and conclusions
A new parameter, the effective buoyancy ratio, B eff was defined to characterize the dynamics of thermo-chemical convection occurring in the Earth's mantle.Buoyancy ratio, B, in its classical meaning (Davaille et al., 2002) is a constant and predicts the resistance of the lower compositionally dense but hot layer against mixing, and so it is insensitive to the behavior of the system.Additionally, our calculations show that mixing also occurs in the case of B > 1 suggesting the instability of two-layer convection for arbitrary value of B (Davaille, 1999a).On the other hand, the effective buoyancy ratio, B eff , is a time-dependent parameter which represents the instantaneous stability of the thermo-chemical system.During the numerical modeling, B eff decreases monotonically and when its value reaches 1 (the instability point, when stabilizing chemical and the destabilizing thermal buoyancy is balanced), the pattern of the flow system changes considerably; the two-layer convection is replaced with a one-layer thermo-chemical mixing.Thus, the effective buoyancy ratio is a good diagnostic tool to determine the actual state of two miscible fluid-like layers.Additionally, the time of B eff = 0 might define the time when the two layers are mixed.B eff illustrates well the evolution of the initial dense layer above the CMB consisting of four phases: (i) transition phase of warming dense layer; (ii) erosion and dilution of the dense layer; (iii) effective thermo-chemical convection (mixing of layers); (iv) homogenization.
There is no exact knowledge of the early stage of the mantle evolution, that is, the initial condition of the numerical and laboratory models is pending.Still, there are some hypotheses which make it plausible that the compositionally dense layer might have formed during the first 100-200 Myr of the Earth's history.Deschamps and Tackley (2008) and Tackley (2012) detail this problem and mention that (1) mixing between the liquid iron of the outer core and silicate deep mantle (Mao et al., 2006); (2) early differentiation of the magmatic mantle (Agee and Walker, 1988;Kellogg et al., 1999); (3) crystallization of the basal magma ocean (Labrosse et al., 2007;Lee et al., 2010) forming iron-rich and therefore dense material would result in a compositionally dense layer above the core-mantle boundary in the early phase of the Earth's evolution.On the contrary, there are arguments for mechanisms that the dense material generates over time, e.g., through the recycling and segregation of oceanic crust (e.g., Christensen and Hofmann, 1994).Undoubtedly, there are arguments for and against the formation of a dense layer in the deepest mantle at the beginning of the Earth's evolution, though the dense material accumulated in the Archean deepest mantle is a geologically realistic initial condition.Additionally, it is popular because it can be implemented easily both in numerical (e.g., Van Summeren et al., 2009;Li et al, 2014) and laboratory models (e.g., Jellinek and Manga, 2002;Gonnermann et al., 2002).
In order to make a comparison among different numerical models, Tackley (2012) rescaled the results for the heat expansion of α = 10 −5 K −1 as a more realistic value in the deep, compressible mantle (Mosenfelder et al., 2009).Applying smaller heat expansion requires less initial compositional density contrast to obtain the same B eff .Re-scaling our model (Fig. 1) for reduced heat expansion minimum β = 3 %, initial compositional density contrast is needed to maintain the dense layer over the age of the Earth.It is in accordance with the results of Tackley (2012), who arrived at a density difference of 2-3 % based on different model calculations.
Using tomographic likelihoods, Trampert et al. (2004) separated the total density variation in the mantle into temperature and chemical density variation.They established that the present compositional density variation is dominant in the lower 1000 km of mantle and it is likely to exceed 2 %.It corresponds to our models with initial density contrast of β = 3 % assuming reduced heat expansion because the density difference decreases gradually due to erosion and dilution processes (Fig. 2).
Several normal modes of the Earth show a significant sensitivity to the density/shear velocity ratio in the deep mantle (Koelemeijer et al., 2012).Ishi and Tromp (2004) revealed a total density increment of approximately 0.5 % beneath Africa and the Pacific in which the opposite effect of temperature and compositional variation is superimposed.Taking into account the compositional density increase of more than 2 % and the total density increase of only 0.5 %, a rough esti-mate of the effective buoyancy ratio gives a value of slightly above 1.Based on our model results at this stage, the TCC system in the Earth's mantle might be just before the instability point.It agrees well with the present strongly deformed, disintegrated morphology of LLSVPs (e.g., Garnero et al., 2007a).
The Supplement related to this article is available online at doi:10.5194/se-15-93-2015-supplement.

Figure 1 .
Figure 1.Five stages characterizing the evolution of the thermo-chemical convection as a funktion of non-dimensional time.Left: time series of monitoring parameters (heat flux, velocity, temperature, concentration, see in Table 2), vertical lines denote the stages shown in the right side.Right: the evolution of the concentration of the dense material and the temperature field.

Figure 2 .
Figure 2. (a) The concentration, (b) the temperature differences between the lower and upper layers and (c) the effective buoyancy ratio as a function of time at different values of the initial compositional density contrast, β.Dashed blue line denotes the complex model (see in text).

Figure 3 .
Figure 3. (a) Concentration of the dense material and (b) temperature field demonstrating the processes of (3) erosion and (4) dilution of the dense layer.Black arrows denote the logarithmically scaled (a) mass flux of the dense material and (b) flow velocity.

Figure 4 .
Figure 4. (a) Occurrence time of the two most characteristic events: B eff = 1 (onset of mixing) and B eff = 0 (end of mixing) as well as (b) slope of the decrease of the concentration difference and the effective buoyancy ratio during the erosion/dilution phase as a function of β.

Figure 5 .
Figure 5.A quasi-stationary state of the (a) temperature, (b) viscosity, (c) concentration of the dense material and (d) velocity for the complex model (see text) at 3.5 Gyr (t = 0.01325).Viscosity is non-dimensionalized as log10(η/η S ) where η S denotes the surface viscosity.