A Long-lived Lunar Magnetic Field Powered by Convection in the Core and a Basal Magma Ocean

Almost any scenario for the formation of the Moon involves enough energy to melt much of the newly formed Moon (e.g., Hartmann & Davis 1975; Warren 1985; Elkins-Tanton et al. 2011; Canup 2012; Ćuk & Stewart 2012; Nakajima & Stevenson 2014). The resulting magma ocean is often modeled as solidifying in three primary stages (e.g., Elardo et al. 2011; Wieczorek et al. 2006; Hess & Parmentier 1995; Hamid & O'Rourke 2022). As the lunar magma ocean cooled, dense mafic cumulates (e.g., olivine and pyroxene) formed and sank toward the bottom. Once most of the lunar magma ocean solidified, anorthositic plagioclase with lower density began to crystallize, rising to form the lunar crust. The final, highly evolved liquids, "ur-KREEP" (enriched in uranium, thorium, potassium, rare earth elements, and phosphorus), alongside ilmenite-rich cumulates, would be gravitationally unstable because of their high densities. Some fraction of this ur-KREEP–ilmenite mixture eventually sank to the base of the mantle, ponding as a layer above the core–mantle boundary (CMB). Radiogenic heat from elements present in this fallen ur-KREEP layer, such as uranium, thorium, and potassium (with concentrations up to ∼12 times higher than the bulk mantle), could fully melt this layer (e.g., Scheinberg et al. 2018). The result is a basal magma ocean (BMO) that persists until convective heat loss into the overlying mantle causes solidification. The nominal model of Scheinberg et al. (2018) had a 301 km peak thickness BMO; less conservative models had BMO thicknesses up to 450 km.

Models are equivocal about the lifetime of a BMO. For example, a small compositional density contrast between the BMO and the overlying mantle could make the BMO short-lived (Stegman et al. 2003). In this scenario, thermal expansion of the BMO can overcome the compositional density contrast between the BMO and the overlying mantle, causing the BMO to buoyantly rise and remix with the mantle. Conversely, the persistence of interstitial fluid trapped within the solidified cumulates could leave the BMO sequestered at the CMB (Elkins-Tanton et al. 2011; Scheinberg et al. 2018). Indeed, interpretations of geophysical (Khan et al. 2014), seismic (Weber et al. 2011), and gravity (Williams et al. 2014) data have indicated that a deep-seated zone of partial melt at the CMB may exist today. This partial melt could be the last remnant of a once-thicker BMO.

The lunar BMO could have sustained a dynamo if it was vigorously convecting and had an electrical conductivity, σ, of several thousand siemens per meter (S m⁻¹) (Scheinberg et al. 2018). Such a BMO dynamo would have an advantage over the core in terms of generating strong crustal fields because it is closer to the surface (e.g., Ziegler & Stegman 2013). Magnetic fields attenuate rapidly with distance, so a magnetic field generated in the BMO would appear stronger at the surface than a magnetic field generated with the same strength in the core (e.g., Scheinberg et al. 2018; Stevenson 1983; Christensen 2010). While sufficiently high conductivity is a challenge for this hypothesis, thermal coupling between the BMO and core can fortunately be explored regardless of this uncertainty.

Our study is built on the whole-Moon models presented in Scheinberg et al. (2018). That study focused on the thermal evolution of the solid mantle and BMO to explain the early, strong (i.e., >10 μT) lunar dynamo. Both the BMO and the core were assumed to be well mixed on the timescales of the overlying solid mantle convection and to have an adiabatic temperature gradient, except during the phase in which the magma ocean increases in temperature. That study further tested the sensitivity of their model to the reference viscosity in the solid mantle, the fraction of the KREEP layer that remained near the surface, and the fraction of radioactive material concentrated in the BMO. At the start of their simulations, the BMO exhibited a rapid increase in heat flow from radiogenic heating, followed by a steady decline to its solidus temperature. A detailed model of the core was not included because the core is relatively small and does not strongly affect the thermal evolution of the BMO and solid mantle. In this study, we do not directly model the BMO-hosted dynamo but rather focus on the core to test if models of lunar evolution that feature a BMO as a boundary condition can explain both the strong, early dynamo and the later dynamo that produced much weaker fields (Figure 1).

Zoom In Zoom Out Reset image size

We assume that the lunar core is an iron alloy that starts fully liquid with no chemical or thermal stratification. To build our models, we assume that sulfur is the major light element in the core, given its siderophile behavior and cosmochemical abundance (e.g., Pommier 2018; Cameron 1973). Our models also include trace amounts of potassium as a source of radiogenic heating. Other studies have speculated about the possible roles of other light elements in the lunar core, including carbon (e.g., Dasgupta et al. 2009), silicon (e.g., Berrada et al. 2020), and phosphorous (e.g., Yin et al. 2019). However, the complexities of a core with multiple light elements are beyond the scope of this study.

A 1D parameterized description of the structure of the core is the foundation of our models. As described in Appendix A, we used the hydrostatic equilibrium and equations of state detailed in Khan et al. (2018) to calculate the radial profiles of density, pressure, temperature, and gravitational acceleration within the core. Our fiducial structural model assumes that the core contains 6 wt% sulfur and has a central pressure and temperature of 5.15 GPa and 1800 K, respectively, to match the core parameters described in Scheinberg et al. (2018). The radius of the core is then 350 km, which is also the same as in Scheinberg et al. (2018) and in agreement with available observational constraints. However, Scheinberg et al. (2018) used an average density for the core appropriate to a composition of pure iron, which would increase the total mass of the core by ∼20%. Fortunately, most of the structural parameters that are key to our thermodynamic calculations (e.g., K₀, K₁, L_ρ , and A_ρ in Table D1) are not sensitive to the bulk composition of the core. Sulfur is most important to the thermal evolution of the lunar core via its effect on the bulk liquidus. Using a fixed sulfur content to calculate other parameters (e.g., ρ₀, P₀, and M_C ) should only introduce inaccuracies that are smaller than the observational uncertainties.

The overlying BMO controls the evolution of the core. In the models of Scheinberg et al. (2018), the BMO is initially set to 1700 K at 4.2 Ga, heats up for ∼200 Myr due to radiogenic heating, and subsequently cools until it reaches the initial temperature when the models are stopped. We start tracking the evolution of the core at the time when the BMO starts cooling again. At that time, we assume the core is fully molten and has an adiabatic temperature gradient throughout. We set the "initial" temperature at the top of the core equal to that at the bottom of the BMO. From the results of Scheinberg et al. (2018), we know the total heat flow across the CMB (Q_CMB) over time:

$$\begin{eqnarray}&&{Q}_{{\rm{CMB}}}={Q}_{B}-{Q}_{{S}_{{\rm{BMO}}}}-{H}_{{\rm{BMO}}}.\end{eqnarray} \tag{ 1 }$$

Here Q_B is the heat flow outward from the BMO into the solid mantle, ${{\rm{Q}}}_{{{\rm{S}}}_{\mathrm{BMO}}}$ is the heat associated with secular cooling, and H_BMO is the radiogenic heating in the BMO. In order to model the magnetic history of the Moon until the present day (i.e., after the BMO model has stopped), we further assume that Q_CMB changes linearly to a specified present-day value, which could be the same or (much) less than the value of Q_CMB when the BMO solidifies. With the boundary condition provided by the BMO model, we then use a well-established method, developed to study Earth's core (e.g., Labrosse 2015), to model the thermodynamic evolution of the lunar core once it starts cooling again. First, we can calculate the global heat budget of the core:

$$\begin{eqnarray}&&{Q}_{{\rm{CMB}}}={Q}_{L}+{Q}_{G}+{Q}_{R}+{Q}_{S}.\end{eqnarray} \tag{ 2 }$$

Here Q_S represents the secular cooling of the core and is proportional to the core's specific heat. We assume that trace amounts of potassium produce radiogenic heating (Q_R ). The remaining two terms are only relevant once the inner core nucleates: energy from latent heat (Q_L ) and gravitational energy from the exclusion of light elements into the outer core (Q_G ) that are released as the inner core grows.

Given the total heat flow, we solve for the rate of change in the CMB temperature. As shown in Appendix B, most of the terms on the right side of Equation (2) are products of d T_CMB/d t and a term ($\tilde{Q}$) that depends only on the thermodynamic properties of the core and its structural parameters. Each of those terms is calculable using polynomial functions. We can thus rearrange Equation (2):

$$\begin{eqnarray}&&\displaystyle \frac{{{dT}}_{{\rm{CMB}}}}{{dt}}=\displaystyle \frac{{Q}_{{\rm{CMB}}}-{Q}_{R}}{\widetilde{{Q}_{S}}+\widetilde{{Q}_{G}}+\widetilde{{Q}_{L}}}.\end{eqnarray} \tag{ 3 }$$

The growth rate of the inner core is directly proportional to d T_CMB/d t also (see Appendix B). Because Equation (2) does not include any secular cooling of the inner core, we are implicitly assuming that the inner core is perfectly insulating (i.e., with zero thermal conductivity). We could also model a conductive inner core with infinite thermal conductivity, but the associated heat flow is a minor contribution to the global heat budget if the inner core extends to only <75% of the core radius, as expected at the present day. Technically, Equations (2) and (3) are only valid if the liquid portion of the core is convective and thus maintaining a nearly adiabatic thermal profile. This assumption is not valid at the present day since thermal stratification probably exists, because the core heat flux was likely lower than the heat flux that can be conducted along the adiabat for most of the Moon's evolution (e.g., Laneuville et al. 2014).

Our models use a liquidus for the core that depends on the bulk composition. We adapted Equation (29) from Buono & Walker (2011), in which the Fe–FeS liquidus is fit to a polynomial that is fourth order in both pressure and sulfur content. Our model uses an approximation of the liquidus that is first order in both pressure and sulfur content. Specifically, we estimated the approximate pressure derivative (d T_L /d P) based on the difference in the liquidus temperatures at 5.15 GPa at the center of the core versus 4.43 GPa at the CMB for 6 wt% sulfur. We found the approximate compositional derivative (d T_L /d c) based on the difference in liquidus temperatures for zero versus 25 wt% sulfur at 5 GPa (Table D1).

Vigorous convection in the core can produce a dynamo through the conversion of kinetic to magnetic energy. In general, there are two types of power sources for convection in the core. First, the buoyancy of light elements released from inner core solidification can drive compositional convection. Second, thermal buoyancy from secular cooling of the core, freezing of the inner core, and/or radiogenic heating can power thermal convection. For thermal convection to occur from secular cooling alone, Q_CMB must exceed the adiabatic heat flow (Q_AD), which equals the product of the thermal conductivity of the core and the adiabatic temperature gradient (see Appendix B). Once the inner core nucleates, the critical heat flow above which convection occurs is lowered.

We combined the energy and entropy budgets for the core to calculate the total dissipation available to power a dynamo (e.g., Labrosse 2015):

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{CMB}}}={{\rm{\Phi }}}_{L}+{{\rm{\Phi }}}_{G}+{{\rm{\Phi }}}_{R}+{{\rm{\Phi }}}_{S}-\ {{\rm{\Phi }}}_{K}.\end{eqnarray} \tag{ 4 }$$

Here Φ_L , Φ_G , Φ_R , and Φ_S are the dissipation terms associated with Q_L , Q_G , Q_R , and Q_S , respectively. The last term (Φ_K ) corresponds to the entropy sink associated with thermal conduction in the core. Appendix B contains the polynomial expressions for each dissipation term, which, like the energy terms, depend on the thermophysical properties of the core and its overall cooling rate. Critically, we assume a dynamo exists if the dissipation is positive (i.e., if Φ_CMB > 0 W). This criterion yields similar predictions as another often-used criterion, which is that the magnetic Reynolds number (defined below) exceeds a critical value of 50–100 (e.g., Roberts 2007).

Several scaling laws are available to convert the dissipation (in Watts) into the strength of the magnetic field at the equatorial surface of the Moon (in Teslas). First, we use a scaling law based on core energetics (see Appendix B) to calculate the total dipole moment (D_M) of the Moon (units of A m²). In this case, assuming the lunar magnetic field is dipolar, the surface field strength at the magnetic equator is

$$\begin{eqnarray}&&B=\displaystyle \frac{{\mu }_{0}{D}_{{\rm{M}}}}{4\pi {R}_{{\rm{M}}}^{3}},\end{eqnarray} \tag{ 5 }$$

where R_M is the radius of the Moon, and μ₀ is the permeability of free space. Additionally, we estimate the magnetic field intensity using three scaling laws that relate the associated convective power to the anticipated convective velocities (e.g., Christensen 2010). These scaling laws use different force balances to calculate the strength of the magnetic field in the core (B_C ). First, mixing-length (ML) theory assumes a balance between inertial and buoyancy forces,

$$\begin{eqnarray}&&{B}_{{\rm{ML}}}={\left[2c{\mu }_{0}{\left({\rho }_{0}{R}_{C}^{2}{{\rm{\Phi }}}_{{\rm{CMB}}}^{2}\right)}^{\tfrac{1}{3}}\right]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 6 }$$

where c ∼ 0.63 is a constant of proportionality, ρ₀ is the central density in the core, and R_C is the radius of the core. Second, assuming a balance of Coriolis, inertial, and gravitational (Archimedes; CIA) forces yields

$$\begin{eqnarray}&&{B}_{{\rm{CIA}}}={\left[2c{\mu }_{0}{\left({\rho }_{0}^{2}{R}_{C}^{4}{\rm{\Omega }}{{\rm{\Phi }}}_{{\rm{CMB}}}^{3}\right)}^{\tfrac{1}{5}}\right]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 7 }$$

where Ω is the present-day angular velocity of the Moon, which may underestimate the field strength, since the Moon likely rotated faster in the past. Third, the magneto-Archimedes-Coriolis (MAC) scaling assumes a balance between Lorentz, gravitational, and Coriolis forces:

$$\begin{eqnarray}&&{B}_{{\rm{MAC}}}={\left[2c{\mu }_{0}{\left({\rho }_{0}{R}_{C}^{2}{\rm{\Omega }}{{\rm{\Phi }}}_{{\rm{CMB}}}\right)}^{\tfrac{1}{2}}\right]}^{\tfrac{1}{2}}.\end{eqnarray} \tag{ 8 }$$

With these three scaling laws, we calculate the surface field strength of the dipole component as

$$\begin{eqnarray}&&{B}_{S}=\displaystyle \frac{1}{7}{B}_{C}{\left(\displaystyle \frac{{R}_{C}}{{R}_{M}}\right)}^{3}.\end{eqnarray} \tag{ 9 }$$

The ratio of the Moon's core radius to the Moon's radius (R_M) accounts for the fact that the dipole field at the surface is smaller than the dipole field at the core (Scheinberg et al. 2018). The prefactor of 1/7 assumes an Earth-like power spectrum for the magnetic field and accounts for the fact that not all of the energy in the magnetic field is partitioned into the poloidal components that can reach the surface (e.g., Christensen et al. 2009; Scheinberg et al. 2018). Note that the core field is assumed to diffuse across an electrically insulating mantle in this approach, thus neglecting the contribution of the BMO. Because the BMO is argued to have a relatively large electrical conductivity, our surface field strength calculations may be considered as lower-bound estimates (discussed further in Section 4.3).

We further assess the dipolarity of the Moon's magnetic field, particularly whether a dipole-dominated or multipolar dynamo may be preferred. Although there are numerous hypotheses for what controls the breakdown of the dipole (e.g., Soderlund & Stanley 2020), we consider here the local Rossby number,

$$\begin{eqnarray}&&\begin{array}{c}{{\rm{Ro}}}_{l}=\displaystyle \frac{U}{2\ {\rm{\Omega }}\ l},\end{array}\end{eqnarray} \tag{ 10 }$$

where Ω = 2π/T is the angular velocity of the Moon, T is the rotation period in seconds, l is the characteristic length scale of the flow, and U is the characteristic fluid velocity. This dimensionless parameter measures the relative importance of inertial to Coriolis forces at convective length scales. Numerical models of planetary dynamos indicate that dipole-dominated solutions tend to be found approximately when Ro_l < 0.1 (i.e., when inertial effects are relatively weak compared to rotation), with multipolar solutions occurring for larger Ro_l values (e.g., Christensen & Aubert 2006).

In order to estimate this parameter, we assume a characteristic fluid velocity and length scale following scaling law predictions as done for the magnetic field strengths (e.g., Christensen 2010). The ML scaling yields

$$\begin{eqnarray}&&\begin{array}{l}{U}_{{\rm{ML}}}={\left(\displaystyle \frac{{\phi }_{{\rm{CMB}}}{R}_{C}}{{\rho }_{0}}\right)}^{\tfrac{1}{3}},\ \ \ \ \ \ {l}_{{\rm{ML}}}={R}_{C},\\ {{\rm{Ro}}}_{l,{\rm{ML}}}={\left(\displaystyle \frac{{\phi }_{{\rm{CMB}}}{R}_{C}}{{\rho }_{0}}\right)}^{\tfrac{1}{3}}{(2\ {\rm{\Omega }}\ {R}_{C})}^{-1};\end{array}\end{eqnarray} \tag{ 11 }$$

the CIA scaling yields

$$\begin{eqnarray}&&\begin{array}{c}{U}_{{\rm{CIA}}}={\left(\displaystyle \frac{{\phi }_{{\rm{CMB}}}}{{\rho }_{0}}\right)}^{\tfrac{2}{5}}{\left(\displaystyle \frac{{R}_{C}}{{\rm{\Omega }}}\right)}^{\tfrac{1}{5}},\ \ \ \ \ \ {l}_{{\rm{CIA}}}={\left(\displaystyle \frac{{{U}_{{\rm{CIA}}}\ R}_{C}}{{\rm{\Omega }}}\right)}^{\tfrac{1}{2}},\\ {{\rm{Ro}}}_{l,{\rm{CIA}}}={\left(\displaystyle \frac{{\phi }_{{\rm{CMB}}}}{{\rho }_{0}}\right)}^{\tfrac{2}{5}}{\left(\displaystyle \frac{{R}_{C}}{{\rm{\Omega }}}\right)}^{\tfrac{1}{5}}{\left(4\ {\rm{\Omega }}\ {U}_{{\rm{CIA}}}\ {R}_{C}\right)}^{-\tfrac{1}{2}};\end{array}\end{eqnarray} \tag{ 12 }$$

and the MAC scaling yields

$$\begin{eqnarray}&&\begin{array}{l}{U}_{{\rm{MAC}}}={\left(\displaystyle \frac{{\phi }_{{\rm{CMB}}}}{{\rho }_{0}{\rm{\Omega }}}\right)}^{\tfrac{1}{2}},\ \ \ \ \ {l}_{{\rm{MAC}}}={R}_{C}\\ {{\rm{Ro}}}_{l,{\rm{MAC}}}={\left(\displaystyle \frac{{\phi }_{{\rm{CMB}}}}{{\rho }_{0}{\rm{\Omega }}}\right)}^{\tfrac{1}{2}}\ {(2\ {\rm{\Omega }}\ {R}_{C})}^{-1}.\end{array}\end{eqnarray} \tag{ 13 }$$

Here ϕ = Φ_CMB/V_OC is the volumetric thermodynamically available power over the fluid core. We could also use these velocity scalings to confirm that the magnetic Reynolds number, which relates the Ohmic diffusion timescale to the convective timescale, exceeds the critical value of ∼50 for magnetic field generation to occur (e.g., Roberts 2007). With the definition

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{Re}}}_{m}\,=\,{\mu }_{0}\sigma {Ul},\end{array}\end{eqnarray} \tag{ 14 }$$

a flow velocity faster than ∼0.1–1 mm s⁻¹ produces Re_m > 50 if we assume that the length scale is equal to the core radius and the electrical conductivity is σ ∼ 10⁵–10⁶ S m⁻¹ (e.g., Berrada et al. 2020; Pommier et al. 2020).

Our model ingests the BMO model outputs from Scheinberg et al. (2018) and calculates the energy and dissipation budgets for the core to determine when the core may host a dynamo (see Table D2). Following the nomenclature of Scheinberg et al. (2018), naming of the BMO models corresponds to the parameters chosen to describe the mantle and the initial solidification of its magma ocean. For example, V19 indicates a reference mantle viscosity of 10¹⁹ Pa s, K50 indicates that 50% of the KREEP layer remained trapped near the surface, and p54 indicates that 54% of the internal radiogenic heating occurs in the sunken KREEP material. We focus on the BMO models that generate magnetic fields with lifetimes of <2.9 Gyr for consistency with the paleomagnetic record (e.g., Mighani et al. 2020). We adopt the nominal BMO case, V19K50p54, as the basis for our nominal model of the core, as it assumes moderate yet reasonable values for the mantle parameters. To test the sensitivity of our models to the properties of the core, we scan across four different parameters: the abundances of sulfur and potassium in the core, the thermal conductivity of the core, and the present-day heat flux at the CMB.

As with other planets, the Moon's core is expected to be an alloy of iron and light elements, such as sulfur (e.g., Steenstra et al. 2016). Properties of the FeS system are relatively well known (e.g., Fei et al. 1997, 2000; Chudinovskikh & Boehler 2007; Morard et al. 2007, 2008; Buono & Walker 2011; Chen et al. 2008; Pommier 2018), and concentrations of sulfur in the lunar core are likely <6–8 wt% based on interpretations of seismic data (e.g., Weber et al. 2011) and models of the lunar core (e.g., Scheinberg et al. 2015; Laneuville et al. 2014). We vary the sulfur abundance, [S], in the bulk core from 1 to 9 wt% in increments of 0.5 wt%.

Potassium is a potential heat source in planetary cores and soluble in iron alloys at planetary conditions (Murthy et al. 2003; Lee et al. 2004). However, the potassium content of the lunar core remains uncertain. Based on previous studies (e.g., Laneuville et al. 2014; Scheinberg et al. 2015), we test a lower limit of 0 ppm, which assumes a complete lack of radiogenic heating in the lunar core. Although the lower pressures and temperatures in the lunar interior might lead to lower amounts of potassium in the lunar core (e.g., Steenstra et al. 2018), we use plausible concentrations of potassium in Earth's core as an upper limit (e.g., Hirose et al. 2013). In our models, we assume that potassium is incompatible in the inner core, meaning that the outer core becomes enriched in potassium as the inner core grows. We vary the bulk potassium abundance, [K], from zero to 50 ppm in increments of 25 ppm.

The thermal conductivity, k_C , of iron alloys defines the adiabatic heat flux of the core. We assume that the maximum plausible value of k_C is ∼50 W m⁻¹ K⁻¹, cited from thermal conductivity experiments on Fe–FeS alloys in the lunar pressure and temperature range (e.g., Pommier 2018). Small amounts of impurities, such as sulfur, can cause a large reduction in the thermal conductivity. We investigate k_C and [S] independently in our models to isolate the effects of each parameter, but they are coupled in reality. A minimum value of 10 W m⁻¹ K⁻¹ is selected to represent relatively large impurities of sulfur (e.g., Pommier 2018). Other proposed compositions for the lunar core, such as Fe–Si alloys, have thermal conductivities that are intermediate between these upper and lower bounds (Berrada et al. 2020). Overall, we vary k_C from 10 to 50 W m⁻¹ K⁻¹ in increments of 10 W m⁻¹ K⁻¹.

The present-day heat flux at the CMB is highly uncertain and may have been susceptible to higher heat fluxes out of the lower mantle from the enrichment of water and other incompatible elements during solidification of the lunar magma ocean (e.g., Elkins-Tanton & Grove 2011; Khan et al. 2014; Evans et al. 2014; Weiss & Tikoo 2014; Dygert et al. 2017; Greenwood et al. 2018). To monitor how the core's temperature evolves given a certain heat flow, we test a range of values using thermal evolution models as a guide (e.g., Laneuville et al. 2014). After the BMO solidifies, we assume that Q_CMB decreases linearly from the final BMO simulation output (∼0.90–3.70 GW) to a heat flux value specified at present. We therefore vary the present-day heat flow, Q_C , from zero to 2 GW in increments of 1 GW. While the lower limit of 0 GW may represent an extreme scenario, we want to explore a full range of modeling possibilities to account for multiple scenarios for the lunar solid mantle. Furthermore, 1D models for small planetary bodies typically indicate that the heat flux varies slightly during most of the core's evolution (e.g., Laneuville et al. 2014). We find that model outputs from simulations with a Q_CMB equal to the final BMO simulation output are similar to those from models where the Q_CMB slightly decreases.

Astute readers will realize that our modeling approach makes the cooling rate of the core seem artificially smooth over time after the BMO solidifies. While the BMO exists, we use Q_CMB from the 3D solid mantle models of Scheinberg et al. (2018), which contain realistic time variability and fluctuations. Once the BMO has presumably solidified, our parameterized model is effectively 1D and uses a simplified approach for Q_CMB to capture the average field strength and lifetime of the core dynamo. In reality, some smaller-scale temporal variations in Q_CMB should be expected, and the very last time step is not necessarily representative of the end of the time series.

We ran a total of ∼800 simulations to test the sensitivity of the core model to [S], [K], k_C , and Q_C using BMO model outputs from Scheinberg et al. (2018) as boundary conditions.

Our nominal values for the core parameters are [S] = 7.5 wt%, [K] = 25 ppm, k_C = 40 W m⁻¹ K⁻¹, and Q_C = 0 GW for the V19K50p54 BMO boundary condition (Table 1). Figure 2 details the outputs of our nominal model for the core coupled to the nominal BMO model (i.e., V19K50p54). The temperature at the CMB begins at ∼1760 K and quickly spikes to ∼1940 K due to radiogenic heating in the BMO (Figure 2(a)). The BMO then begins solidifying as radiogenic heating declines over time, followed by the core cooling in tandem with the BMO. Once the BMO solidifies, an inner core forms at ∼2.2 Gyr as relatively pure iron crystallizes from the inside out (Figure 2(b)), expelling sulfur into the outer core (Figure 2(c)). Our models assume that the lunar core always contains sub-eutectic amounts of sulfur. We verified that this assumption is consistent with our results, which track the sulfur content of the outer core over time (e.g., Figure 2(c)). The liquidus temperature of the outer core is lowered as it is progressively enriched in sulfur. The result is a molten outer core and a growing inner core. The heat flow is always less than that transported by thermal conduction along the core adiabat, Q_AD. After the inner core nucleates, most extracted heat from the core arises from the release of Q_G and Q_L (Figure 2(d)). The release of Q_G is nonzero but small compared to Q_L. Following the release of Q_G and Q_L , there is a reduction in the core cooling rate due to these heat sources acting as a buffer to secular cooling. We note that the Q_CMB is much lower than the heat flow across the upper boundary of the BMO (Q_B = ∼100 GW at 2.6 Ga) in Scheinberg et al. (2018) because Q_B includes radiogenic heating in the BMO and the heat associated with cooling the BMO itself.

Zoom In Zoom Out Reset image size

Abundant sulfur influences the core's ability to drive a magnetic field by lowering its solidus temperature and controlling the onset of inner core crystallization (discussed further in Section 3.2.1). The nominal model produces an inner core radius of 250 km at the present day (Figure 2(b)) and is consistent with core radii derived from calculated models of lunar gravity data (Williams et al. 2014) and reanalyzed Apollo seismic data (Weber et al. 2011).

The lunar BMO suppresses convection in the core by lowering its cooling rate. The core produces a dynamo that begins near the cessation of the nominal BMO-hosted dynamo and ends at ∼1 Ga, consistent with the lower estimate on the cessation of the lunar dynamo derived from radiometric dating of Apollo 15 samples (e.g., Mighani et al. 2020; Figure 3(a)). The relatively weak surface magnetic field strength of ≲2.55 μT is also consistent with the paleomagnetic data and intensities from previous models of the lunar core dynamo (e.g., Laneuville et al. 2014; Tikoo et al. 2014, 2017; Evans et al. 2018; Mighani et al. 2020).

Zoom In Zoom Out Reset image size

We next consider different BMO conditions for our core model. Table 1 presents the nominal core input parameters for each BMO boundary condition used in this study. The BMOs with a smaller fraction of KREEP that remained near the surface (i.e., V19K25p54 and V18K00p100 in Table 1) have greater initial thicknesses and tend to require lower sulfur abundances (6.5–7 wt%) in the bulk core to initiate dynamo action during the observed timing of the low-intensity epoch. Because a BMO with a greater thickness will have a longer lifetime (e.g., Scheinberg et al. 2018), the core will begin crystallizing at a later time, when the BMO eventually solidifies. Conversely, models with shallower BMOs (i.e., 301 km) mostly require higher sulfur abundances in the core (7–8.5 wt%) to achieve a core dynamo during the same period. The BMO boundary conditions with greater lifetimes additionally suppress inner core growth for longer periods, resulting in smaller inner core radii at the present day. Furthermore, models that contain shallower BMOs match the estimated timing of the lunar dynamo if balanced by less radiogenic heating in the core (i.e., ≤25 ppm of potassium). In general, BMO boundary conditions typically require the core to have a higher thermal conductivity (i.e., ≥30 W m⁻¹ K⁻¹) to match the estimated timing of the lunar dynamo.

3.2.1. Influence of Sulfur in the Core

An inverse relationship exists between the sulfur content and the solidus temperature of the core. As the sulfur content increases, the solidus temperature of the Fe–S system decreases, delaying core solidification until lower temperatures are reached. Therefore, the timing of inner core growth, and thus the start time of compositional convection in our models, depends on the sulfur content of the bulk core (Figure 4(a)). The sulfur concentration is viable when the end of the core-hosted dynamo matches the lower estimate on the cessation of the lunar dynamo at ∼1 Gyr (e.g., Mighani et al. 2020). Initial sulfur abundances of 1–6.5 wt% result in inner core nucleation at higher temperatures, causing the core to solidify rapidly early in its history (Figure 4(b)). Sulfur abundances from 7 to 8.5 wt% result in the inner core nucleating near the cessation of the BMO-hosted dynamo. Increasing the bulk sulfur content to >8.5 wt% further delays inner core growth and generally results in temporal gaps between the BMO- and core-hosted dynamos, a complete lack of core dynamo action, or contradictions with timing estimates derived from paleomagnetic data (Figure 4(c)). However, if the BMO model assumes a lower solid mantle viscosity (i.e., V18.5K50p54), then convective heat transfer is more efficient and results in shorter BMO lifetimes (Scheinberg et al. 2018). As a result, the inner core begins crystallizing earlier, and a bulk sulfur content of up to 12 wt% can produce results consistent with the lower estimate on the cessation of the lunar dynamo (e.g., Mighani et al. 2020). The trends outlined in Figure 4 that arise from variations in k_C , [K], [S], and Q_C continue under all other BMO boundary conditions.

Zoom In Zoom Out Reset image size

The duration and intensity of the core dynamo are also sensitive to k_C , Q_C , and [K] (Figure 4). A potassium abundance of 50 ppm in the core contributes thermal energy to the dynamo but suppresses growth of the inner core, which can decrease the predicted intensity of the magnetic field overall. Decreasing [K] has a minimal effect on the field intensity because radiogenic heating is nearly equivalent to secular cooling in the dissipation budget. In contrast, increasing the total heat flow to 1–2 GW increases the duration and strength of the core-hosted dynamo, unless a low sulfur abundance leads to rapid core solidification. Furthermore, the duration and intensity of the field generally increase with decreasing thermal conductivity values. We find that purely thermal convection typically occurs before the onset of inner core crystallization if the thermal conductivity is low (i.e., 10–30 W m⁻¹ K⁻¹ as in Figure 3(c)). As thermal conductivity decreases, the superadiabatic heat flow increases, leading to a stronger, more long-lived dynamo. Thermal convection–driven dynamos typically occur simultaneously with BMO-hosted dynamos, as the core is still hot and fully molten. Compared to the abundance of sulfur in the bulk core, our simulations reveal that small variations in parameters such as k_C , Q_C , and [K] play an overall negligible role in the onset of a compositionally driven dynamo, whereas a thermal convection–driven dynamo is largely dictated by k_C .

Depending on the BMO boundary condition, inner core crystallization can produce fields of ∼0.7–2.3 Gyr in duration, with peak magnetic fields of 0.16, 0.45, 0.77, and 4 μT for the ML, CIA, BF, and MAC scaling laws, respectively (Table 2). A general issue arises in the case of the CIA, ML, and BF scaling laws in which the strengths are not sufficient to reproduce the first period of decline to ∼4–7 μT by 3.19 Ga (e.g., Strauss et al. 2021) or the second period of decline to ∼5 ±2 μT by ∼1–2 Ga (e.g., Tikoo et al. 2017; Mighani et al. 2020). However, intensities ranging from ∼1.7 to 4 μT can be achieved under all BMO boundary conditions if the MAC scaling law is assumed. In particular, an intensity of 4 μT is achieved if the BMO boundary condition contains a lower fraction of radioactive material concentrated in the BMO (i.e., V19K50p36). However, this magnetic field weakens to a maximum of 4 μT ∼2 Gyr after accretion, which is ∼0.7 Gyr later than what is observed in the lunar paleomagnetic record (Strauss et al. 2021).

Surface magnetic fields are weaker if they are driven by thermal convection rather than inner core crystallization. The peak surface magnetic field driven by thermal convection in the nominal core model is 3 × 10⁻⁶, 3 × 10⁻⁷, 0.001, and 0.003 μT for the ML, CIA, BF, and MAC scaling laws, respectively (Table 2). For all BMO boundary conditions, thermal convection in the core is initiated ∼3.7 Gyr ago (albeit briefly in some models; e.g., Figure 3(a)). Furthermore, depending on the BMO boundary condition, thermal convection can persist intermittently for up to ∼1.7 Gyr, resulting in an overlap with the BMO-hosted field (e.g., Figure 3(a)). Thermal convection produces intensities that are consistent with previous modeled estimates of the core (e.g., Laneuville et al. 2014; Evans et al. 2018; Scheinberg et al. 2015) but inconsistent with paleomagnetic analyses constraining the initial and final decline of the lunar dynamo (e.g., Tikoo et al. 2017; Mighani et al. 2020; Strauss et al. 2021). Furthermore, these results are consistent with a low-intensity epoch that persisted from ∼1.9 to 0.8 Ga (e.g., Mighani et al. 2020; Tikoo et al. 2017, Tikoo et al. 2014; Strauss et al. 2021).

An uneven heat flow across the CMB may make the magnetic field intermittent because dynamos can be sensitive to slight variations in heat flow (Scheinberg et al. 2015). As an artifact of our modeling approach, early magnetic fields produced via thermal convection are discontinuous due to fluctuations in the Q_CMB from mantle dynamics. In some cases, thermal convection generates fields that are predicted to drop to zero multiple times before rising again from inner core crystallization. The duration of these gaps in the magnetic field is much longer than the magnetic diffusion time (Appendix C). Using the nominal models but with core conductivity lowered to k_C = 30 W m⁻¹ K⁻¹ as an example case (i.e., Figure 3(c)), gaps in thermal convection last, on average, ∼140 Myr, whereas the magnetic diffusion time is only a few hundred years. Alternatively, dynamos induced by thermal convection can transition directly into those induced by inner core crystallization, compounding the resultant fields.

In order to make initial predictions for the magnetic field morphologies in our models, we estimate the local Rossby number as a proxy for whether the core dynamos would be dipole-dominated or multipolar, as has been done for Ganymede's dynamo (Ruckriemen et al. 2015), for example. The CIA scaling law predicts higher values of the local Rossby number (Ro_l ∼ 10⁻²) relative to the ML (Ro_l ∼ 10⁻³) and MAC (Ro_l ∼ 10⁻⁴) scaling laws, since inertia plays a larger role in the force balance (Christensen & Aubert 2006). However, for the nominal core model, all scaling laws predict that the local Rossby number is below the threshold value of ∼0.1 throughout the lifetime of the core dynamo, suggesting a prevailing dipole-dominated magnetic field (Table 2 and Figure 5).

Zoom In Zoom Out Reset image size

Future studies could model more complex modes of crystallization in the lunar core. To recap, we made two relevant assumptions. First, we assumed that the core always contains sub-eutectic amounts of sulfur, which most of our models indeed predict (Section 3.2.1). Second, we assumed that the core solidifies from the center outward. We set the liquidus temperature to increase faster than the adiabatic temperature with pressure (e.g., with gradients of 30 versus ∼23–25 K GPa^–1, respectively).

Future studies could relax these two assumptions, which would produce more complicated behavior in models (e.g., Hauck et al. 2006). First, FeS rather than Fe could crystallize from the outer core as it cooled if the sulfur content were supereutectic. Being sulfur-rich compared to the residual liquid, solid FeS would float to the top of the liquid rather than sink to form an inner core like solid Fe. Second, solidification could occur at the top or middle of the outer core, rather than at its bottom. For example, "iron snow" could occur in metallic cores if the liquidus crosses the adiabat above the base of the outer core. This process could help drive a dynamo as the solidified iron sinks and remelts in the warmer fluid below, leading to compositional convection (e.g., Williams 2009; Breuer et al. 2015). Whether the Moon's core entered an FeS crystallization or Fe snow regime at any time remains an ongoing question.

Scientists might make more realistic models of the thermal evolution of sulfur-rich cores if they include these processes. Such models require detailed phase diagrams for the Fe–FeS system. The neglect of Fe snow and FeS crystallization in our models does not change our takeaway message, however, that the presence of a BMO overlying the core may influence the timing and intensity of the core dynamo. Our models may interface with future, more detailed descriptions of Fe snow and FeS in the core.

The geometry and paleo-orientation of the Moon's magnetic field remain largely uncertain. Estimates of paleoinclinations from five Apollo samples suggest the existence of a dipolar field and a paleopole located at ∼75°N between 3.8 and 3.3 billion yr ago (e.g., Cournède et al. 2012). These findings are possibly best explained with a paleofield geometry close to the present-day rotation axis of the Moon. Assumptions of the paleopole were made based on the location of the Apollo samples; samples collected from the northern hemisphere were given a positive declination, while samples collected from the southern hemisphere were given a negative declination. However, the sign of the inclination remains largely unknown, and more data are required to confirm interpretations made from lunar samples. Studies of Apollo 17 mare basalts estimated an inclination of ∼34° based on the layering of its parent boulder (Nichols et al. 2021). This inclination is consistent with, but does not require, a dipole in the center of the Moon aligned along its rotation axis.

Conversely, Olson & Christensen (2006) hypothesized that the Moon's magnetic field may have been multipolar rather than dipole-dominated. The critical difference between our studies is the amplitude of buoyancy flux (BF) in the core. Their study assumed that the average BF associated with convection in the lunar core was 0.3 times the terrestrial value. That is, F_Moon = 0.3 F_Earth, where F = αgQ/(ρ C_p ) with thermal expansivity α, gravitational acceleration g, convective heat flux Q, density ρ, and specific heat capacity C_p . This assumption was based on the idea that tidal dissipation could add several TW of power to the ancient lunar core (e.g., Williams et al. 2001). This larger heat flow leads to larger estimates of the local Rossby number (e.g., Ro_l ∼ 2), which would shift the lunar dynamo into a multipolar regime. In contrast, our models do not include tidal heating in the lunar core. So the total power available for convection is only several GW in our models, as shown in Figure 2(d).

It is also possible that the directional magnetization of lunar rocks does not record a long-term orientation of the lunar magnetic field, since differential rotations between the mantle and core would cause a core dipole field to drift across the lunar surface (e.g., Ćuk et al. 2019). Relative motions of the core and mantle or misalignment between the lunar dynamo and spin axis may further explain the great variability in the inferred orientation of the lunar dynamo from proposed paleopole locations (e.g., Oliveira & Wieczorek 2017; Nayak et al. 2017). Future missions sampling the lunar bedrock along varying latitudes will allow for more precise geometric determinations of the Moon's magnetic field.

The effects of an electrically conducting lower mantle on the core dynamo is not considered in our study. This limitation is significant for several reasons. First, as noted in Section 2.3, the relatively large conductivity of the BMO, especially when it is fully liquid, will likely cause our estimates of surface magnetic field strengths to be artificially small compared to if this conductivity were taken into account. Our estimates for the core field strength assume that the entire mantle, including the BMO, is electrically insulating such that the core-generated magnetic field becomes a potential field that diffuses upward through the mantle. Given the anticipated higher conductivity of metalliferous silicate melts compared to solid mantle rocks (e.g., Scheinberg et al. 2018), the top of the dynamo region may effectively be the top of the BMO, rather than the top of the core, even if the BMO is subcritical for dynamo action.

Second, fluid flows within the BMO may also modulate the core field itself (e.g., Gómez-Pérez et al. 2010). Conversely, if the BMO fluid is stably stratified, its presence may still filter out small-scale components of the core field that rapidly vary via the magnetic skin effect (e.g., Christensen 2006). Third, the BMO may have resulted in a larger magnetic coupling between the core and mantle in the past, relevant to studies of the Moon's rotational dynamics over time (e.g., Dumberry & Wieczorek 2016). Further work, such as numerical dynamo modeling, is needed to better understand the full degree of coupling between the BMO and core of the Moon.

The effects of thermal stratification in the lunar core are not considered in this study. The inclusion of thermal stratification can have several effects on the heat flux at the CMB. Studies of Mercury's core (e.g., Knibbe & van Westrenen 2018; Knibbe & Van Hoolst 2021) found that thermal stratification can lead to an increased inner core size, higher temperatures, and a larger heat flux at the CMB, which together result in an early start to the magnetic field. Subsequent heat released upon core solidification would enable slow core growth and an active magnetic field until the present day. Future work could apply these models of Mercury to the Moon.

Thermal stratification is probably inevitable at the present day but could also exist early in the Moon's history. In this study, we assume that the core was initially fully molten and had an adiabatic temperature gradient. If radiogenic heating in the BMO ever made the bottom of the BMO hotter than the top of the core, then heat would move from the BMO into the core, which would cause thermal stratification at the top of the core that might delay the start of a core-hosted dynamo. However, the Moon could have formed with "superheat" (such that the core was initially hotter than the BMO; e.g., Evans et al. 2018), in which case the core could deliver heat to the BMO even while the BMO heats up radiologically. Neither our study nor Scheinberg et al. (2018) modeled these two countervailing possibilities in detail. Further work is thus needed to better understand the formation and early evolution of the Moon.