Abstract
The early lunar mantle overturn, associated with the sinking of the dense ilmenite-bearing cumulate (IBC) crystallized at the shallow lunar mantle, provides satisfactory explanations for the origination of high-Ti basalt, the abnormally strong magnetic field between ∼ 3.9 and ∼ 3.6 Ga and the low-viscosity zone in the deep lunar mantle, but still poses a debate regarding the initial state of IBC in the early lunar mantle. If the sinking of IBC initiated before the end of lunar magma ocean crystallization, the solidified IBC can acquire a greater thickness and a higher initial velocity at the IBC-mantle boundary. The variation of initial velocity can affect the strain rate of IBC and, correspondingly, the dislocation creep components at the shallow lunar mantle. In this work, we analyze the effects of initial velocity on the dynamics of early lunar mantle by using the theory of Rayleigh-Taylor instability. To couple the effects of diffusion creep and dislocation creep for all major minerals in the lunar mantle, we exploit an improved Minimized Power Geometric (IMPG) model and isostress mixing model to characterize the upper limit and lower limit for the viscosity of the lunar mantle comprising four major minerals, i.e. olivine, orthopyroxene, clinopyroxene and ilmenite. The modeling results suggest that a high initial velocity, in any case, can shorten the onset time, tending to promote the early lunar mantle overturn even in a rheologically-strong lunar mantle. The effect of initial velocity on the overturn wavelength shows a strong dependence on the rheological mixing model. For the isostress mixing model, the increase of initial velocity tends to elongate the overturn wavelength. For the IMPG mixing model, the overturn wavelength is insensitive to the variation of initial velocity. As the actual lunar mantle rheology sandwiches between the rheologies predicted by isostress mixing model and IMPG model, it can be anticipated that the increase of initial velocity tends to elongate the overturn wavelength. In consideration of the importance of the initial velocity on the dynamics of early lunar mantle, future investigations should focus on the dynamics of the solid IBC in the solidifying lunar magma ocean.
Export citation and abstract BibTeX RIS
1. Introduction
The crystallization of lunar magma ocean determined the structure in the early lunar mantle, but could generate a dense ilmenite-bearing cumulate (IBC) at the shallow lunar mantle that was prone to sink downwards, which in turn caused a chemical-composition overturn in the early lunar mantle. Such an early lunar mantle overturn provides satisfactory explanations for the origination of high-Ti basalt (Ringwood & Kesson 1976; Zhong et al. 2000; Zhang et al. 2013), the time evolution of lunar paleomagnetism (Stegman et al. 2003; Li et al. 2019) and the low-viscosity zone in the deep lunar mantle (Harada et al. 2014, 2016), but still poses a problem regarding the dynamic condition needed for the overturn itself. Stagnant lid, i.e. the cold and stiff part on the top of lunar mantle, can entrap the IBC and prevent this dense layer from sinking downwards. In order to allow the overturn to take place, the rheology of early lunar mantle, or the IBC at least, must be weak enough (Yu et al. 2019; Zhao et al. 2019; Li et al. 2019).
Nevertheless, the previous works also posed a debate on the initial state of IBC in the early lunar mantle. The dynamic modeling for the early lunar mantle overturn was initiated at the time for which the lunar magma ocean just solidified, i.e. ∼ 4.38 and ∼ 4.42 Ga (e.g. Boyet et al. 2015). The initial conditions, in most cases, were prescribed according to the thermodynamic modeling for the lunar magma ocean crystallization (Snyder et al. 1992; Elkins-Tanton et al. 2011; Lin et al. 2017). Such kind of thermodynamic modelings just accounted for the chemistry in the magma ocean and prescribed the dense cumulate piles to be stationary all the time. Correspondingly, those dynamic modelings also prescribed implicitly a stationary IBC at the initial time (de Vries et al. 2010; Yu et al. 2019; Zhao et al. 2019). In contrast, some investigators suggested that the IBC would have mobilized before the end of lunar magma ocean crystallization (Parmentier & Hess 1999; Boukaré et al. 2018; Li et al. 2019). When ilmenite crystallizes, the coexisting anorthosite can contribute to the thickening of lunar crust. Such an effect can significantly decelerate the crystallization of lunar magma ocean (Solomatov 2007), thus allowing a sufficiently long time for the crystallized IBC to develop some small-scale diapirs that are prone to sink downwards at a slow rate (Li et al. 2019). Until the whole lunar magma ocean solidifies, such a presolidification sinking tends to thicken the solidified IBC and possibly assigns a high initial velocity for the IBCmantle boundary.
The initial velocity at the IBC-mantle boundary would significantly affect the dynamics of early lunar mantle in particular for a mixed mantle rheology coupling the diffusion creep and dislocation creep of lunar mantle minerals. In the actual lunar mantle, the deformation of minerals can take place in two different ways, diffusion creep and dislocation creep. Diffusion creep is dominated by the migration of vacancies on the crystal lattice, which results in a linear relation between stress and strain rate. In contrast, dislocation creep, associated with the migration of dislocations on the crystal lattice, results in a nonlinear relation between stress and strain rate, which in turn causes a feedback between the mobility of mantle materials and the change of lunar mantle rheology (Schubert et al. 2001; Karato 2008). The initial velocity at the IBC-mantle boundary can significantly affect the initial strain rate of IBC and, correspondingly, the dislocation creep components at the shallow lunar mantle. Hence, the variation of the initial velocity at the IBC-mantle boundary can also modify the rheology of IBC and, correspondingly, perturbs the onset of early lunar mantle overturn. Nevertheless, such an effect was never constrained in the past. This work aims at examining the effects of the initial velocity at the IBC-mantle boundary on the onset of early lunar mantle overturn by using the theory of Rayleigh-Taylor instability.
2. Methods
2.1. Theory of Rayleigh-Taylor Instability
We initiate the modeling from the lunar mantle before the early lunar mantle overturn and, then, continue the modeling by tracking the development of the sinking of IBC. The sinking of IBC in the early lunar mantle can be described by utilizing the theory of Rayleigh-Taylor instability. According to the theory (Whitehead 1988; Hess & Parmentier 1995), the downwelling velocity of the IBC-mantle boundary at the predominant wavelength can be written as a function of time, i.e.
where w is the downwelling velocity at the IBC-mantle boundary, w0 is the initial velocity, t is the time and tonset is the onset time. The onset time is given by
where ηm is the viscosity of the underlying lunar mantle, ηibc is the viscosity of IBC, Δ ρ is the density difference between IBC and the underlying lunar mantle, i.e. Δ ρ = ρibc – ρm where ρibc is the density of IBC and ρm is the density of the underlying lunar mantle, gibc is the gravitational acceleration at the IBC-mantle boundary and Dibc is the thickness of IBC. The predominant wavelength over which the overturn tends to develop is given by
We assign ηm with the viscosity underneath the IBC and ηibc with the volumetrically averaged viscosity of IBC. Before the onset of lunar mantle overturn, the heat transfer in the lunar mantle is dominated by heat conduction. The temperature profile, needed for constraining the rheology of lunar mantle, is thus determined from the heat conduction equation, i.e.
where ρ is the density, cpm is the heat capacity of lunar silicate portion, T is the temperature, t is the time, Km is the heat capacity of lunar silicate portion and Q is the heat production rate. Here we set the lunar surface temperature as a constant 250 K and the core-mantle boundary (CMB) temperature to be identical to the local solidus temperature. The model above has been applied successfully in predicting the rheological conditions under which the early lunar mantle overturn can take place, which shows a good agreement with the results obtained from geodynamic modeling (Yu et al. 2019).
Strain rate is needed to determine the viscosities of the dislocation creep components (see Sect. 2.3). The lunar crust remains buoyant on top of the lunar silicate portion all the time and, hence, results in a zero velocity at the IBC-crust boundary. Therefore, the strain rate of IBC can be given by
Because of the continuity at the IBC-mantle boundary, we assume the strain rate underneath the IBC, i.e. , to be identical to .
2.2. Parameters
All parameter values are summarized in Table 2.2. Here we consider a Moon with a core radius of 390 km, in agreement with the latest estimation from lunar laser ranging experiment, i.e. 384 ± 93 km (Viswanathan et al. 2019). The thermal conductivity and heat capacity of the lunar silicate portion are set to their typical values of 4.0 W m−1 K−1 and 1600 J kg−1 K−1 respectively. In a previous companion work (Yu et al. 2019), the software alphaMELTS was used to model the crystallization of a lunar magma ocean. The initial conditions of density, temperature and heat production rate are plotted in Figure 1. By prescribing a core radius of 390 km following the latest estimation of Viswanathan et al. (2019), a bulk composition of O'Neill (1991) and an elevated TiO2 content of 0.4 wt. %, i.e. the maximum TiO2 content for the bulk Moon (Buck & Toksoz 1980; Elkins-Tanton et al. 2011; Morgan et al. 1978; Snyder et al. 1992), the modeling yielded a crust thickness of 43 km and an IBC thickness of 36 km. By averaging volumetrically the modeling results by introducing a three-layer configuration, the densities of lunar crust, IBC and the underlying lunar mantle were estimated to be 2715, 4102 and 3204 kg m3 respectively. The heat production rate of the three layers is given by
where QQ is the initial heat production rate and τ is the half-decaying period. The values of QQ are 3.11 × 10−10, 4.35 × 10−11 and 2.77 × 10−12 W kg−1 for the lunar crust, IBC and the underlying lunar mantle respectively. These data are utilized to characterize the chemical-composition stratification after a chemical crystallization of lunar magma ocean. The modeling also yielded a temperature profile by the end of lunar magma ocean crystallization (see Fig. 1(b)). We initiate the temperature modeling by this temperature profile. The time evolution of the temperature profile can be referred to in figure 8 of Yu et al. (2019).
Table 1. Values of Parameters Applied in the Modeling
Parameter | Meaning | Value | Unit |
---|---|---|---|
rm | lunar radius | 1740 | km |
rc | core radius | 390 | km |
Dcr | crust thickness | 43 | km |
Dibc | IBC thickness | 36–150 | km |
cpm | heat capacity of lunar silicate portion | 1600 | J kg−1 K−1 |
Km | thermal conductivity of lunar silicate portion | 4.0 | W m−1 K−1 |
gibc | gravitational acceleration | 1.6 | m s−2 |
ρcr | lunar crust density | 2715 | kg m−3 |
ρibc | IBC density | 3479–4102 | kg m−3 |
ρm | mantle density | 3204 | kg m−3 |
Qcr | heat production rate of lunar crust | 3.11 × 10−10 | W kg−1 |
Qibc | heat production rate of IBC | 1.75–4.35 × 10−11 | W kg−1 |
Qm | heat production rate of lunar mantle | 2.77 × 10−12 | W kg−1 |
τ | half-decaying period | 1.64 | Ga |
w0 | initial downwelling velocity | 10−20–10−5 | m s−1 |
To achieve an accurate estimation on the viscosity of IBC and the viscosity of underlying lunar mantle, a mineralogy-dependent rheology is considered in this work. Hence, the exact mineralogy of IBC and of the underlying lunar mantle need to be known. We repeat the modeling as done by Yu et al. (2019) and then estimate the mineralogy of these two parts. The composition of the chemically-crystallized IBC is determined as 41.97 vol.% olivine + 45.16 vol.% clinopyroxene + 10.39 vol.% ilmenite + 2.48 vol.% minor minerals (whitelockite, quartz, etc.). The composition of the underlying lunar mantle is determined as 63.41 vol.% olivine + 27.20 vol.% orthopyroxene + 9.2 vol.% clinopyroxene + 0.19 vol.% minor minerals.
We then consider the thickened IBC associated with the pre-solidification sinking. Exactly speaking, constraining the dynamics of such a pre-solidification sinking requires tracking the behaviors of a two-phase fluid, which is still technically challenging at present. Alternatively, Li et al. (2019) evaluated the thickening effect based on Stokes law for the velocity of the tiny diapirs. Given a two-order viscosity contrast between IBC and the underlying lunar mantle, the estimation yielded a maximum thickness of 150 km for the final solidified IBC. The viscosity of IBC applied here coincides with the experiments on the rheology of olivine-ilmenite mixture, i.e. the weak rheology of ilmenite can reduce the viscosity of IBC by up to ∼ 2 orders of magnitude in the case of isostress mixing, the weakest case for the rheology of multiphase rock (Dygert et al. 2016). We note that this thickness value would be overestimated because the authors used the timescale of lunar crust formation, i.e. ∼ 200 Myr, relating to the anorthosite crystallization (Elkins-Tanton et al. 2011), as the time span for the pre-solidification sinking. Based on the crystallization modeling, the initiation of anorthosite crystallization is earlier than the initiation of ilmenite crystallization (Snyder et al. 1992; Elkins-Tanton et al. 2011; Lin et al. 2017). Hence, the actual time allowing the pre-solidification sinking would be shorter. Nevertheless, the time span of ilmenite crystallization was never constrained through isotopic chronology. For this reason, we still treat 150 km as a plausible upper limit for the thickness of IBC. Based on volume conservation, the density of the thickened IBC is determined as 3479 kg m−3. The heat production rate of the thickened IBC is determined as 1.75 × 10−11 W kg−1. The composition of the thick IBC is determined as 56.86 vol.% olivine + 18.89 vol.% orthopyroxene + 20.19 vol.% clinopyroxene + 3.17 vol.% ilmenite + 0.89 vol.% minor minerals.
The initial velocity w0 depends on the states of solidified IBC in the solidifying lunar magma ocean. If the sinking of IBC has initiated before the end of lunar magma ocean crystallization, the solidified IBC underneath the residual melt would firstly develop some small-scale diapirs that are prone to sink downwards (Hess & Parmentier 1995). Until the whole lunar magma ocean solidifies, the IBC-mantle boundary can acquire a higher initial velocity from those early-formed diapirs. By prescribing the weakest case for the rheology of IBC, Li et al. (2019) estimated the velocity of those tiny diapirs to be 0.65 km Myr−1, or ∼ 10−4 m s−1 based on Stokes law. This velocity likely characterizes the maximum case for the velocity of those small-scale diapirs. Correspondingly, we use a much smaller value of 10−20 m s−1 to characterize the velocity of IBC-mantle boundary if the final solidified IBC is nearly stationary.
2.3. Rheology: the Weakest End Member
Based on the discussion in Section 2.2, IBC and the underlying lunar mantle are actually composed of multiple minerals. In each mineral, diffusion creep and dislocation creep can act simultaneously. To estimate the effective viscosity of a multiphase rock, the so-called rheological mixing model, which determines how stress and strain rate are partitioned among different deformation components, must be exploited. Nevertheless, there is still no perfect rheological mixing model to describe the rheology of a multiphase rock. Instead, we follow the bounding method and provide the upper limit and lower limit for the viscosity of lunar mantle. In this section, we will discuss the isostress mixing model used to characterize the lower limit.
For a single deformation component, the relationship between stress and strain rate can be determined by the Arrhenius power-law (e.g. Hirth & Kohlstedt 2003), i.e.
where the subscript k specifies the deformation component, e.g. the dislocation creep of clinopyroxene, is the strain rate, Ak is a pre-factor, σk is the stress, d is the grain size, nk and mk are two empirical coefficients, Ek is the activation energy, p is the pressure, Vk is the activation energy, R is the gas constant and T is the temperature. Based on the definition of viscosity, i.e. , the effective viscosity of a single deformation component can be determined as
Because of the low pressure in the lunar mantle, i.e. ∼ 4 GPa at maximum, the pressure has an insignificant effect on the mantle rheology and the dynamics of early lunar mantle overturn (Yu et al. 2019). Hence, we assume here pV = 0.
For the isostress mixing model, every deformation component bears the same strain, but partitions the total stress based on the arithmetic average (Reuss 1929), i.e.
Based on the definition of viscosity, i.e. , the bulk effective viscosity of isostress mixing model can then be determined as
In reality, isostress mixing model corresponds to a transverse loading style and yields the weakest end member for the rheology of a multiphase rock (Tullis et al. 1991; Tullis & Wenk 1994). However, transverse loading requires the orientation of all minerals to be perpendicular to the direction of stress, which is too idealized for the actual condition in the lunar mantle. Hence, we rely on the isostress mixing model to estimate the lower limit for the viscosity of lunar mantle.
2.4. Rheology: A Mathematically-defined Upper Limit
Here we consider the other rheological mixing model improved from the Minimized Power Geometric (MPG) mixing model to characterize the upper limit for the viscosity of a multiphase rock. The MPG model was firstly proposed to provide an approximation for inter-mineralogical mixing (Huet et al. 2014). This model has two versions, MPGe and MPGs. For the MPGe model, the total strain rate is partitioned based on the law of geometric mean, i.e.
where is the strain rate of the i-th mineral and ϕi is the volumetric fraction of the i-th mineral. The viscosity of the MPGe model can then be determined by minimizing the deformation power, i.e.
where ηi is the viscosity of the i-th mineral. For the MPGs model, the total stress is partitioned based on the law of geometric mean, i.e.
where σi is the stress partitioned by the i-th mineral. Similarly, the bulk viscosity of MPGs model can also be obtained by minimizing the deformation power. Based on the discussion of Huet et al. (2014), MPGe model and MPGs model yield nearly identical results for the viscosity of a multiphase rock and, thus, can be considered to be equivalent. Correspondingly, the geometric mean partitioning for the total strain rate results in the geometric mean partitioning for the total stress under the minimization of deformation power and vice versa. Based on this equivalency, the MPG model yields a bulk viscosity in the form of
By comparing the bulk viscosity estimated by the MPG model with the experimental data for the rheology of multiphase rocks, Huet et al. (2014) suggested that the MPG model can predict well the rheology of a multiphase rock on the level of inter-mineralogical mixing.
For a single mineral, its deformation can take place in diffusion creep and dislocation creep at the same time. As these two deformation components relate to different imperfect structures on the crystal lattice, their strain rates are essentially additive (Karato 2008). For this reason, if the MPG mixing model is reliable in predicting the inter-mineralogical mixing, the total strain rate and total stress need to be partitioned based on
where and σi,dif are the strain rate and stress partitioned by the diffusion creep of the i-th mineral respectively, and and σi,dis are the strain rate and stress partitioned by the dislocation creep of the i-th mineral respectively. As the MPG model can provide a good approximation for the rheological mixing at the inter-mineralogical level, the accurate bulk viscosity of a multiphase rock obeys a form like
where the subscripts dif and dis still specify diffusion creep and dislocation creep respectively.
Equation (18) corresponds to a complex strain rate partitioning law that is hard to be solved analytically. Hence, a numerical minimization algorithm is needed to search for the best combination of strain rates for every deformation component, yielding the minimum deformation power. Nevertheless, the minimization algorithm with multiple strain rates as variables is usually very time consuming, which limits its applications in the numerical modeling. Alternatively, we extend the MPG model into the level of deformation component mixing and yield the improved MPG (IMPG) model to estimate the upper limit for the bulk viscosity of a multiphase rock. We firstly write the total strain rate as
In order to exploit the mathematical schema for deriving the MPGe model, we define here a parameter and a target function
Note that and P' do not have any physical meanings, but are just defined to fulfill the mathematical form allowing the usage of the MPGe mixing model. By minimizing the target function P' through the Lagrangian multiplier method, we can obtain a bulk viscosity
where and is the viscosity of the k-th deformation component with the total strain rate.
Here we illustrate the mathematical meaning of Equation (21). Due to the equivalency between MPGe model and MPGs model, Equation (19) results in under the minimization of P'. Hence, the bulk viscosity estimated by the IMPG model essentially obeys the form
In comparison to Equation (18), Equation (22) still preserves the form of MPG mixing model on the level of inter-mineralogical mixing, but estimates the viscosity of the i-th mineral by taking the geometric mean for the viscosity of diffusion creep component and the viscosity of dislocation creep component. Hence, the IMPG model tends to overestimate the viscosity of the i-th mineral. Correspondingly, Equation (21) also overestimates the viscosity of a multiphase rock. For this reason, we can use ηIMPG as a plausible upper limit. Here we note that the IMPG model does not have any physical meaning, but just yields a mathematically-defined upper limit for the viscosity of a multiphase rock.
2.5. Rheological Parameters
In this work, we just consider a dry lunar mantle because the rheology of wet ilmenite was never measured through laboratory experiments. As we initiate the modeling from the time at which the lunar magma ocean solidifies, we also assume the lunar mantle to be melt-free. To simplify the setting of rheological parameters, we only consider the diffusion creep and dislocation creep of four major minerals in the lunar mantle, i.e. olivine, orthopyroxene, clinopyroxene and ilmenite, in the water-free case. Other minor minerals are considered as olivine and are assigned with the rheological parameters of olivine. Table 2 lists the rheological parameters for the diffusion creep and dislocation creep of these four major minerals. The diffusion creep of ilmenite was never measured via a laboratory experiment. Instead, we use the rheological parameter of magnetite diffusion-creep because of its similar crystal structure with respect to ilmenite (Dygert et al. 2016).
Table 2. Rheological Parameters of Olivine, Orthopyroxene, Clinopyroxene, Magnetite and Ilmenite
Component (k) | Ak (s−1 μmmk MPa−nk ) | mk | nk | Ek (kJ mol−1) | References |
---|---|---|---|---|---|
Olivine diffusion creep | 1.59 × 109 | 3.0 | 1.0 | 375 | Hirth & Kohlstedt (2003) |
Olivine dislocation creep | 1.10 × 105 | 0 | 3.5 | 530 | Hirth & Kohlstedt (2003) |
Orthopyroxene diffusion creep | 1.92 × 104 | 2.0 | 1.0 | 320 | Tasaka et al. (2013) |
Orthopyroxene dislocation creep | 6.6 × 1013 | 0 | 3.8 | 820 | Mackwell (1991) |
Clinopyroxene diffusion creep | 1.26 × 1015 | 3.0 | 1.0 | 560 | Bystricky & Mackwell (2001) |
Clinopyroxene dislocation creep | 6.31 × 109 | 0 | 4.7 | 760 | Bystricky & Mackwell (2001) |
Magnetite diffusion creep | 19.95 | 3.0 | 1.0 | 188 | Till & Moskowitz (2013) |
Ilmenite dislocation creep | 6.36 | 0 | 3.0 | 307 | Dygert et al. (2016) |
The grain size in the early lunar mantle was poorly constrained in the past. However, by prescribing a melt-free and olivine-dominated lunar mantle, the interpretations on lunar free libration suggest a grain size of ∼ 10 mm at present (Nimmo et al. 2012). Here we consider the grain size to be three values, 2.6 mm, 5.6 mm and 12.0 mm. We also define a reference viscosity at 1600 K as an indicator for the grain size based on the diffusion creep of olivine, i.e.
where Tr is the reference temperature 1600 K. The three grain sizes above yield the reference viscosities 1019, 1020 and 1021 Pa s respectively, in line with the condition of reference viscosity, i.e. 5 × 1019–1021 Pa s, that can reproduce the present lunar mantle temperature in the geodynamic modeling (Zhang et al. 2013).
As two examples, Figures 2 and 3 plot the variation of bulk viscosity over temperature for the chemically-crystallized thin IBC and thickened IBC with different combinations of reference viscosity, strain rate and rheological mixing model. For both of the two mixing models, the viscosity of IBC always tends to decrease with the increase of temperature. Increasing the strain rate can significantly decrease their viscosities for all temperatures, implying a positive feedback between the mobility and the weakening of IBC's rheology. The difference between the upper limit viscosity and lower limit viscosity is 2–4 orders for a temperature of 1000 K, but decreases to 1–2 orders for a temperature of 1500 K.
Download figure:
Standard imageDownload figure:
Standard image3. Results
The modeling begins from the time at which the entire lunar magma ocean just solidifies. We run the modeling by tracking the time evolution of onset time and overturn wavelength. The reference viscosity is varied between 1019 and 1021 Pa s. The initial velocity is varied between 10−20 and 10−5 m s−1. To ensure the stability of numerical modeling, we apply a time step of 10−4 Myr. However, for the isostress mixing model, a high initial velocity can cause a very high strain rate and correspondingly, unrealistically low viscosities for the IBC and the underlying lunar mantle, which in turn results in a numerical breakdown. For this reason, we just assign two initial velocity values, 10−20 and 10−15, for the cases utilizing isostress mixing model.
Figure 4 displays the modeling results for the sinking of the thin IBC with the IMPG mixing model. For all three reference viscosities, the increase of initial velocity tends to shorten the onset time and the sinking of IBC turns out to be possible when the initial velocity is as high as 10−5 m s−1. The overturn wavelength is insensitive to the variation of initial velocity and, furthermore, the variation of reference viscosity. In all cases, the overturn wavelength falls at the level of ∼ 4 × 10−3 at the initial time, and then tends to decrease over time from ∼ 10−1 Myr. Such a trend of normalized wavelength variation suggested that the IMPG mixing model tends to facilitate an overturn developing in small-scale instabilities.
Download figure:
Standard imageFigure 5 shows the modeling results for the sinking of thick IBC with the IMPG mixing model. For all three reference viscosities, the increase of initial velocity still tends to shorten the onset time. The sinking of thick IBC turns out to be possible when the initial velocity increases to 10−10 m s−1 or higher values. The overturn wavelength is still insensitive to the variation of initial velocity and the variation of reference viscosity. In all cases, the normalized wavelength falls to the level of ∼ 8 × 10−3 initially and then tends to decrease over time from ∼ 0.1 Myr. Such a trend of wavelength variation still favors an overturn developing in small-scale instabilities.
Download figure:
Standard imageFigure 6 depicts the time evolution of onset time and wavelength for the sinking of thin IBC with the isostress mixing model. For all three reference viscosities, the increase of the initial velocity tends to shorten the onset time of overturn. The time evolution of onset time always shows a stable value at the initial phase, then tends to increase over time and finally tends to decrease over time in a very rapid manner. For a reference viscosity of 1021 Pa s and an initial velocity of 10−20 m s−1, the onset time is close to the regime allowing the sinking of IBC at ∼ 10 Myr. Afterwards, the onset time tends to increase again, but finally tends to decrease rapidly at ∼ 102 Myr. Hence, such a combination of reference viscosity and initial velocity allows a late lunar mantle overturn. In other modeling cases, the onset time becomes lower than the onset time at the very early phase and then does not exceed the instant time any more during the consequential evolution, which implies an early lunar mantle overturn. The overturn wavelength manifests an increase with the increase of initial velocity. For an initial velocity of 10−20 m s−1, the initial normalized overturn wavelengths are ∼ 8 × 10−3, ∼ 10−2 and ∼ 2 × 10−2 respectively for the reference viscosities 1019, 1020 and 1021 Pa s. For a higher initial velocity 10−15 m s−1, the initial normalized wavelengths are ∼ 4 × 10−2, ∼ 7 × 10−2 and ∼ 10−1 respectively for the three reference viscosities. In all modeling cases, we observe an increase of overturn wavelength from ∼ 10 Myr.
Download figure:
Standard imageFigure 7 features the time evolution of onset time and wavelength for the sinking of thickened IBC with the isostress mixing model. The increase of the initial velocity tends to shorten the onset time of overturn. For all three reference viscosities, isostress mixing model always allows the overturn to take place even for a low initial velocity between 10−20 and 10−15 m s−1. The overturn wavelength tends to increase with the increase of initial velocity. For a reference viscosity of 1019–1021 Pa s, an initial velocity of 10−15 m s−1 results in a normalized overturn wavelength between ∼ 6 × 10−3 and ∼ 10−2, whereas an initial velocity of 10−15 m s−1 results in a normalized overturn wavelength between ∼ 2 × 10−2 and ∼ 6 × 10−2. For the same initial velocity, the overturn wavelength still tends to increase with the increase of reference viscosity.
Download figure:
Standard image4. Discussions
Our modeling results suggest that the initial velocity at the IBC-mantle boundary can significantly affect the onset of early lunar mantle overturn, but its effects depend strongly on the rheological mixing model. For both isostress mixing model and IMPG mixing model, increasing the initial velocity can significantly decrease the onset time. As the actual mantle rheology sandwiches between the upper limit predicted by the IMPG model and the lower limit predicted by the isostress mixing model, it can be anticipated that the increase of initial velocity always tends to shorten the onset time, indicating the trend to promote the early lunar mantle overturn. For the weakest isostress mixing model, an initial velocity of 10−20 m s−1 allows the overturn to take place unless (1) the reference viscosity is as high as 1021 Pa s, and (2) the thickness of IBC is as small as 36 km, while an initial velocity > 10−15 m s−1 always allows the overturn to take place for all three reference viscosities. For the IMPG mixing model, the overturn is only able to take place if the initial velocity is as high as 10−5 m s−1 for the chemically-crystallized thin IBC, or is higher than 10−10 m s−1 for the thickened IBC. Such a result implies the necessity of a high initial velocity in activating the early lunar mantle overturn in the case of IMPG mixing model.
In the previous works, a rheologically-weak lunar mantle with a reference viscosity ≤ 1020 Pa s was suggested to be a prerequisite for activating the early lunar mantle overturn (Yu et al. 2019; Li et al. 2019; Zhao et al. 2019). Based on the Rayleigh-Taylor instability analysis as done by Yu et al. (2019), the weakening of lunar mantle rheology can shorten the onset time needed for the development of overturn, which allows the chance for the onset time to be lower than the instant time. In this case, the downwelling velocity can increase continuously over time, and then the early lunar mantle overturn tends to be possible. Our modelings above suggest that a high initial velocity at the IBC-mantle boundary, owing to the pre-solidification sinking of IBC, can also generate a similar effect in the context of a mixed rheology coupling the diffusion creep and dislocation creep of all major minerals in the lunar mantle. Moreover, even for the IMPG mixing model and a reference viscosity of 1021 Pa s, which yields the upper limit of mantle viscosity, a high initial velocity can still facilitate the early lunar mantle overturn.
The high initial velocity at the IBC-mantle boundary, as we referred to in Section 2.2, is associated with the sinking of IBC before the end of lunar magma ocean crystallization. Based on the estimation of Li et al. (2019), such a mechanism can assign an initial velocity of ∼ 10−4 m s−1 at maximum. Nevertheless, Li et al. (2019) used the weakest rheology of olivine-ilmenite mixture associated with the isostress mixing model. As isostress mixing model relates to the transverse loading style that is too idealized for the lunar mantle, the actual initial velocity would be lower than 10−4 m s−1. Nevertheless, the IBC-mantle boundary can also acquire a high initial velocity through the early onset of lunar mantle convection, which was not considered well in this work. Maurice et al. (2017) reported the onset of mantle convection during the magma ocean phase. Similarly, if the lunar mantle convection had begun before the end of lunar magma ocean crystallization, the convection would entrain the solidified IBC and, in turn, assign the IBC-mantle boundary with a high initial velocity. At the upper boundary layer of the convective lunar mantle, the downwelling velocity correlates with the Nusselt number of the lunar mantle, i.e. Nu, via the scaling relationship (Melosh 2011)
where κm is the thermal diffusivity of the lunar mantle and Dm is the thickness of lunar mantle. The Nusselt number of the lunar mantle can be estimated by a scaling relationship
where Ra is the Rayleigh number of the lunar mantle and Racr is the critical Rayleigh number, i.e. Racr ∼ 10–2 based on Schubert et al. (2001). The Rayleigh number of lunar mantle can be given by
where αm is the thermal expansivity of lunar mantle. Given αm ∼ 10−5 K−1, κm ∼ 10−6 m2 s−1, g ∼ 1 m s−2, Δ T ∼ 103 K, Dm ∼ 106 m and ηr ∼ 1019–1021 Pa s, the Rayleigh number of the lunar mantle varies between ∼ 104 and ∼ 106. Hence, the initial velocity at the IBC-mantle boundary would be elevated to be 10−11–10−9 m s−1 by the early onset of lunar mantle convection. In this case, we can anticipate that the sinking of thick IBC can always take place, whereas the sinking of thin IBC may still need an isostress mixing model.
The effect of initial velocity on the overturn wavelength also depends on the rheological mixing model. For the IMPG model, the overturn wavelength is insensitive to the variation of the initial velocity at the IBC-mantle boundary. We also observe a normalized wavelength of ∼ 4 × 10−3 m s−1 or lower values for the thin IBC, whereas a normalized wavelength of ∼ 8 × 10−3 m s−1 or lower values for the thick IBC. This result implies that regardless of reference viscosity or initial velocity, the IMPG model always facilitates the overturn developing in small-scale instabilities. For the isostress mixing model, the overturn wavelength tends to increase with the increase of initial velocity. As the actual mantle rheology sandwiches between the rheologies predicted by the two rheological mixing models, we can anticipate that the overturn wavelength tends to increase with the increase of initial velocity. On the other hand, we also observe that for the same reference viscosity and initial velocity, a thicker IBC even tends to shorten the overturn wavelength. A possible reason is that for the same initial velocity, a thicker IBC tends to reduce the strain rate of IBC and the strain rate underneath the IBC. In this case, the dislocation creep components can be strengthened, which in turn affects the viscosity contrast over the IBC-mantle boundary.
The overturn wavelength determines the extent over which the sinking of IBC takes place in the lunar interior. Some investigators suggested that the sinking of IBC would take place in a hemisphere of the Moon (Parmentier et al. 2002). This so-called degree-one overturn, i.e. λ* ≈ 0.5, provides a possible explanation for the focusing of KREEP materials on the lunar nearside and the formation of lunar global asymmetry (Haskin et al. 1998; Korotev 2000; Wieczorek & Phillips 2000; Jolliff et al. 2000). In our modeling results, we do not observe the long-wavelength overturn for the IMPG mixing model. The overturn wavelength is short at the initial time, but even tends to decrease during the consequential evolution. Nevertheless, isostress mixing model allows a short overturn wavelength at the initial time, but the overturn wavelength tends to increase over time during the consequential evolution, which then provides the possibility to generate the degree-one overturn. Because of the uncertainty of lunar mantle rheology, it is still unknown whether degree-one overturn is a robust feature or not. This point still needs to be verified by using geodynamic modeling and a more accurate model for the rheology of lunar mantle.
A much more robust reason for facilitating the long-wavelength overturn, which is not considered well in this work, is the effect of water. As the rheology of wet ilmenite was never measured through experiments, it is still not possible to carry out the Rayleigh-Taylor instability analysis for the sinking of IBC in a hydrous lunar mantle. Nevertheless, recent studies suggest a water content of ∼100 ppm for the bulk Moon (e.g. Saal et al. 2008; Karato 2013). The detection of water in the lunar exosphere, which hints at an active water cycle, also favors the idea that the Moon is releasing its primordial water (Benna et al. 2019). The water retained in the lunar interior is believed to be indigenous, i.e., to have been obtained before solidification of the lunar magma ocean. Otherwise, the rapid formation of the lunar lithosphere would effectively prevent any external sources from delivering water to the deep lunar interior (Hauri et al. 2015). Hence, a wet lunar magma ocean can be expected. Being highly incompatible, water would concentrate in the upper part of the lunar mantle by the end of magma ocean crystallization. In that case, the IBC would be very enriched in water, whereas the underlying lunar mantle would be depleted in water. Such a difference of water content can form a large viscosity contrast over the IBC-mantle boundary, which tends to elongate the overturn wavelength based on Equation (3). Additionally, the viscosity of clinopyroxene strongly decreases with the increasing water content in both diffusion and dislocation creep regimes (Hier-Majumder et al. 2005; Chen et al. 2006). As the solidified IBC contains ∼ 20–45 vol.% clinopyroxene depending on its thickness (see Sect. 2.2), the rheology of wet IBC could be further weakened. In future works, the sinking of a hydrous IBC will have to be studied utilizing experimental constraints on the rheologies of wet olivine, clinopyroxene, orthopyroxene and ilmenite.
5. Conclusions
In this work, we discuss the effects of the initial velocity at the IBC-mantle boundary on the dynamics of early lunar mantle overturn in the context of a mixed rheology coupling the diffusion creep and dislocation creep of all major minerals in the lunar mantle. In order to bound the viscosity of lunar mantle, we design a new model, i.e. the IMPG model, to provide a mathematically-defined upper limit for the bulk viscosity of lunar mantle. The isostress mixing model is applied to characterize the lower limit. We then examine the influences of the initial velocity at the IBC-mantle boundary on the dynamics of early lunar mantle overturn by considering the theory of Rayleigh-Taylor instability. Our modelings suggest that in any cases, increasing the initial velocity at the IBC-mantle boundary tends to shorten the onset time, which implies the trend to promote the early lunar mantle overturn. A high initial velocity, e.g. 10−5 m s−1, allows the overturn to take place all the time, even in a rheologically strong lunar mantle with a reference viscosity > 1020 Pa s. On the other hand, our modelings also imply that increasing the initial velocity also tends to elongate the overturn wavelength. In consideration of the importance of the initial velocity at the IBC-mantle boundary in the early lunar mantle overturn, future investigations should focus on the dynamics of IBC in the solidifying lunar magma ocean.
Acknowledgements
Authors are thankful for the constructive comments from Dr. Dan Zhu. Dr. Sabrina Schwinger at the German Aerospace Center (DLR) provided assistance on repeating the magma ocean crystallization modeling. This work is sponsored by the Fund for Development of Science and Technology of Macau SAR (121/2017/A3).