Formation of the Lunar Primary Crust From a Long-Lived Slushy Magma Ocean.

Classical fractional crystallization scenarios of early lunar evolution suggest crustal formation by the flotation of light anorthite minerals from a liquid magma ocean. However, this model is challenged by the > 200 Myr age range of primitive ferroan anorthosites, their concordance with Mg-suite magmatism and by the compositional diversity observed in lunar anorthosites. Here, we propose a new model of slushy magma ocean crystallization in which crystals remain suspended in the lunar interior and crust formation only begins once a critical crystal content is reached. Thereafter crustal formation occurs by buoyant melt extraction and magmatism. The mixture viscosity strongly depends on temperature and solid fraction driving the development of a surface stagnant lid where enhanced solidification and buoyant ascent of melt lead to an anorthite-enriched crust. This model explains lunar anorthosites heterogeneity and suggests a crustal formation timescale of 100s Ma, reconciling anorthosite ages with an early age of the Moon.

Introduction

The Moon likely formed following a giant impact with the proto‐Earth, dated at ∼4.48 ± 0.02 Ga from Rb‐Sr isotopic data (Halliday, 2008). This theory reconciles the Moon's relatively small core and compositional similarity with the Earth (Canup, 2012; Cuk & Stewart, 2012). Such a highly energetic impact would have produced a hot magma ocean on both the Earth and Moon (Pritchard & Stevenson, 2000; Thompson & Stevenson, 1988). While, on Earth, the processes of early differentiation and crustal formation have been heavily overprinted by ongoing tectonic activity, the Moon retains a unique record of its early evolution. Analyses of Apollo‐era samples have demonstrated that the lunar crust is composed of low‐density ferroan anorthosites (FAN) and is indeed primitive with ages up to 4.55 ± 0.12 Ga (Shearer et al., 2006). The GRAIL gravity survey has produced a high resolution map of the lunar gravitational field which, in conjunction with detailed topography and seismic data, has revealed the lunar crustal structure (Wieczorek et al., 2013). This data indicates that the lunar crust has a low average density of ∼2,550 kg m−3, consistent with an anorthite‐rich composition and a large porosity, and an average thickness of 34–43 km (Wieczorek et al., 2013). The Highlands, composed of more than 90% anorthite, constitute the bulk of the lunar far‐side. Classical models of magma ocean fractional crystallization propose that the crust formed by flotation of buoyant Ca‐rich plagioclase crystals on top of the liquid ocean (Wood et al., 1970). On the nearside, the oldest basaltic lava plain, the Oceanus Procellarum, concentrates incompatible elements and is interpreted as the eruptive product of the residual liquid from the crystallizing magma ocean ‐ denoted urKREEP for enrichments in potassium (K), rare‐earth elements (REEs), and phosphorus (P) ‐ that accumulated below the anorthite‐rich crust (Wieczorek & Phillips, 2000). Although this model reconciles many previous observations, it is currently being challenged by more recent observations and analysis showing a greater duration and temporal variability of lunar crust formation.

Models of magma ocean fractionation with a floating crust predict a relatively short timescale for crust formation of Ma (Elkins‐Tanton et al., 2011) while the dating of the FAN from the lunar Highlands indicates that they form over more than 250 Ma with ages ranging from 4.29 to 4.55 Ga (Alibert et al., 1994; Borg et al., 2020; Marks et al., 2019; Nyquist et al., 2006). The variability and young ages of primitive anorthosites has previously been explained by radiometric clock resetting by metamorphism or impacts (Carlson, 2019). Yet the few samples that show concordant ages from different radiometric systems are the youngest, 4.360 ± 0.003 Ga for FAN 60 025 to 4.304 ± 0.012 Ga for FAN 60 016 (Borg et al., 2020), which argues against thermal or metamorphic resetting (Borg et al., 2011; Carlson et al., 2014; Marks et al., 2019). It has thus been proposed that the Moon was formed late, up to 200 Myr after solar system formation (Borg et al., 2011).

Recently, Maurice et al. (2020) have shown that, by accounting for the low thermal conductivity of anorthosites (Branlund & Hofmeister, 2012) and the heat brought by convection of the deep underlying mafic cumulates, the magma ocean lifetime is prolonged and may last up to 150–200 Myr. Using a fractionation model, the authors show that the flotation scenario however still requires a young age for the Moon of 4.425 ± 0.025 Gyr in order to explain the required isotopic Sm/Nd and Lu/Hf signature of the supposed residual urKREEP liquid (Maurice et al., 2020). However, the small difference in 182W abundances between the Earth and Moon argues for an advanced lunar age of 45–60 Ma after Solar System formation (Thiemens et al., 2019), which is dynamically plausible if the Moon was formed from a large impactor still present in the early Solar System, and is in agreement with the Rb‐Sr age of the Moon (Halliday, 2008).

In comparison with the sum of these constraints, the ∼200 Myr anorthosite age range appears robust. Since only a limited number of lunar anorthosite samples are available to be dated, it is likely that the full spectrum of anorthosite ages has not yet been sampled and that the total duration of anorthosite crustal construction lasted even longer. Thus the minimum ∼200 Myr age range for early lunar anorthosites still poses a challenge to traditional lunar magma ocean crystallization models.

Another challenge for the classical flotation scenario is the concordant ages of primitive anorthosites and Mg‐suite rocks. Classical flotation models suggest that Mg‐suite magmatism would post‐date the anorthosite crust and urKREEP formation, while observations suggest that instead the ages for Mg‐suite rocks range from 4.18 to 4.51 Ga, overlapping those of FANs and continuing for another  Myr (Borg et al., 2020).

Finally, different types of recent observations suggest that FAN may not be globally distributed and that anorthosite rocks forming the lunar crust may be the product of more complex crustal processes than originally envisioned. Lunar feldspathic meteorites appear chemically heterogeneous and in general more magnesian‐rich than FAN. In particular, they do not appear to have crystallized from a common magma source (Gross et al., 2014; Russell et al., 2014). From remote‐sensing studies, it appears that primitive highland rocks are more magnesium rich relative to iron on the farside than on the nearside (Arai et al., 2008; Lucey & Cahill, 2006; Ohtake et al., 2012).

Although the lunar magma ocean concept itself appears necessary if the Moon originated from a giant impact (Pritchard & Stevenson, 2000), a reevaluation of its crystallization history seems essential. Here we propose a new, slushy model of magma ocean solidification, capable of resolving the large age range observed in lunar crustal rocks and explaining the diversity in anorthosites. In contrast to the usual crystal fractionation scenario, in which crystals are assumed to be extracted from the liquid, instantaneously sedimenting up or down according to their buoyancy, we argue that a crystal suspension persists up to and throughout crustal formation since convection is sufficiently vigorous to suppress crystal sedimentation (Solomatov & Stevenson, 1993b). We verify that this slushy scenario is indeed a possibility for the Moon for relevant mixture rheology and crystal sizes (Text S2 in the Supporting Information S1). This slushy scenario is particularly relevant for the lunar mantle because of its low gravity and associated small adiabat and pressure range (Elkins‐Tanton, 2012). Differentiation of the lunar interior instead occurs by upward extraction of a buoyant melt in the upper thermal boundary layer that becomes stagnant once the crystal packing limit is reached. For terrestrial compositions, the melt is buoyant relative to olivine for pressures up to 5–10 GPa (Agee & Walker, 1988), suggesting that the lunar magma ocean could indeed undergo efficient differentiation by melt extraction (Solomatov & Stevenson, 1993a). This scenario of slushy convection and crustal magmatism has not yet been considered for the Moon. Here we model the evolution of a slushy lunar magma ocean and characterize four different dynamical regimes during its crystallization, the results of which constrain the range of time scales of solidification and crust formation.

Model of Slushy Magma Ocean Solidification

The initial, vigorous convection within a magma ocean may act to suppress the formation of a stable crust, with cooling instead driving crystallization throughout the mainly liquid magma ocean. We assume that crystals remain uniformly suspended in the magma ocean because of vigorous convection, which requires crystal smaller than 1 cm in radius according to our analyses (Text S2 in the Supporting Information S1). This agrees well with Suckale et al. (2012) who found that sedimentation in the nonturbulent boundary layers of the magma ocean requires crystals larger than 1 mm in radius, although sedimentation of smaller crystals may be possible by clustering. The transition from the initial stage of rapid resurfacing and a nearly uniform crystal suspension to the subsequent stage of crustal production by buoyant liquid extraction occurs only after the crystal fraction reaches the critical concentration, ϕ  ≃ 0.5–0.6, throughout the depth of the magma ocean. On reaching this critical crystal fraction, the mixture viscosity dramatically increases which may result in a prolongated mushy magma ocean stage as argued by Lebrun et al. (2013) and Monteux et al. (2020) for Mars and the Earth. The resultant mixture viscosity is highly temperature dependent, and we presume that a stagnant‐lid regime of convection develops (Davaille & Jaupart, 1993; Solomatov, 1995), contrary to Lebrun et al. (2013) and Monteux et al. (2020) who instead consider a mobile‐lid regime of convection. Thus, contrary to the classical scenario of crystallization which assumes instantaneous sedimentation and floating of crystals and where the lid thickness is that of accumulating plagioclase, here we assume that crystal suspension persists in the mantle while buoyant liquid extraction takes place in a thermal lid whose thickness is determined by the mixture rheology and cooling history. Since the initial phase of resurfacing is rapid, ∼1,000 years (Elkins‐Tanton et al., 2011; Lebrun et al., 2013), we start our simulations at the rheological transition with an initial well‐mixed crystal content equal to the critical solid fraction . The viscosity of the slush μ is controlled by the temperature‐dependent solid viscosity (Karato & Wu, 1993) and the solid fraction ϕ (Solomatov, 2007), which itself is temperature and concentration dependent;

Here μ is the mixture viscosity at the reference eutectic temperature T , E and E are the thermal and effective solid fraction activation energies, and R is the ideal gas constant.

Previous estimates suggest that at lunar pressures and initial terrestrial bulk silicate compositions, olivine and pyroxene crystals comprise up to 70%–80% volume fraction before plagioclase formation (Elkins‐Tanton et al., 2011; Lin et al., 2016; Longhi, 2003), which indeed suggests a slushy interior existed at the time of crust formation provided that crystals are sufficiently small, as assumed here.

Given the low lunar gravity, the magma ocean pressure range is small. Furthermore, we anticipate that the low pressure phase plagioclase crystallizes within the upper thermal boundary layer, and hence we neglect the pressure‐dependence of the phase diagram. We therefore take the slushy interior to lie along the liquidus of a simplified binary‐eutectic, anorthite‐olivine phase diagram (O’Driscoll et al., 2010) (Figure 1). Hence, we assume that, within the well‐mixed interior magma ocean, the solid fraction , average temperature and bulk concentration in anorthite component C (t) are all constrained by the phase diagram, where is the anorthite composition of the liquid part, T the melting temperature of olivine and m is the liquidus slope (Figure 1 and Table 1 for parameters).

Figure 1

Simplified phase diagram and schematic of a slushy magma ocean. (a) Simplified olivine‐anorthite phase diagram with liquidus T (C) = T  − mC as a function of anorthite concentration C, where T  = 2150 K is the melting point of olivine (forsterite), T  = 1600 K is the eutectic temperature, the eutectic composition C  = 0.7 and the slope of the liquidus m = 785 K (adapted from O’Driscoll et al., 2010). We start with an initial bulk anorthite content between C ≃ 0.05–0.15 and an initial crystal fraction ϕ 0 = 0.5. (b) Schematic of the model with eutectic crust, stagnant lid with melt segregation and a convecting slushy mantle.

Table 1

Model Parameters

Description	Symbol	Value
Moon radius	R L	1,737 km
Basal radius of the slushy magma ocean	c	350–1,237 km
Gravity	g	1.62 m s−2
Mantle density	ρ m	3,000 kg m−3
Thermal conductivity	k	3 Wm−1 K−1
Thermal diffusivity	κ	1 × 10−6 m2 s−1
Heat capacity	c p	1000 J kg−1 K−1
Latent heat of crystallization	L	4 × 105 J kg−1
Thermal expansivity	α	10−5 K−1
Thermal activation energy	E T	300 kJ mol−1
Gas constant	R g	8.314 J mol−1 K−1
Effective solid fraction activation energy	E ϕ	26
Surface temperature	T s	250 K
Olivine melting point	T m0	2150 K
Eutectic temperature	T E	1,600 K
Liquidus slope	m	785 K
Eutectic composition	C E	0.7
Initial solid fraction	ϕ 0	0.5
Heat flux constant	λ 1	0.5
Rheological temperature scale constant	λ 2	2.24
Reference permeability	k 0	10−10−10−12 m2
Liquid viscosity	μ l	1 Pa s
Solid‐Liquid density difference	Δρ	500 kg m−3
Solidus viscosity	μ s	1017−1019 Pa s
Initial bulk anorthite content	C B0	0.05–0.15
magma ocean bulk anorthite content	C B
Stagnant lid thickness	b
Eutectic crust thickness	a
Magma ocean mean temperature	T¯
Magma ocean mean solid fraction	ϕ¯
Liquid anorthite content	C¯

We extend previous temperature‐dependent stagnant lid theories (Davaille & Jaupart, 1993; Solomatov, 1995) to include the evolution of the mean solid fraction. In particular, we write the heat flow F and the temperature difference between the well‐mixed mantle and the base of the stagnant lid as a function of a rheological temperature scale ΔT defined by analogy with the experimental study of Davaille & Jaupart (1993) (Text S1 in the Supporting Information S1);

Here ρ is the magma ocean density, α the thermal expansion expansion coefficient, κ the diffusivity, g the gravity, all assumed constant, , ΔT is the full temperature difference across the lid, 1/γ is the fraction of ΔT controlling convection and Equation 5 is obtained from Equations 1 and 3, λ 1 ≈ 0.5 and λ 2 ≈ 2.24 are constants.

We solve for heat conservation in the convective slushy magma ocean. The properties, and hence the viscosity, of the well‐mixed interior evolve through a balance between secular cooling, a convective heat flux, and the latent heat associated with changes in ,

Here c is the specific heat, is the volume of the convecting magma ocean between depth z and the cumulate or core radius c = 350–1,237 km, is the surface area at depth z, L is latent heat of crystallization, and t is time. Radiogenic heating is examined in Text S1 in the Supporting Information S1.

At the surface, cooling drives enhanced solidification and the development of a very viscous and stagnant lid, where heat is transported by conduction;

Here T is the temperature in the lid, r the radial coordinate. The growth of a lid of thickness b limits crustal growth and is itself controlled by a balance between the heat flux through the lid, which depends on the secular cooling of the lid and is driven by conductive heat losses to the surface, and the heat flux F from slushy convection below (Breuer & Spohn, 2003; Thiriet et al., 2019), where c is the heat capacity and R the lunar radius. If the lid is sufficiently thin and its growth rate sufficiently slow compared to the rate of thermal diffusion in the lid, the temperature profile is approximately linear,

Differentiation of the buoyant anorthite‐rich melt occurs within this stable upper boundary layer through compaction and magmatism, driving the formation of an enriched, near‐eutectic composition crust (Figure 1). The extraction of heat and anorthite‐rich melt through this stagnant lid drives the evolution of the bulk composition (Equation 11) and solid fraction (Equation 2) of the interior, resulting in an ever increasing interior viscosity. For simplicity, we assume that melt is extracted from the slushy magma ocean at the Darcy velocity u , but only if the lid thickness is larger than the crust thickness;

Here we assume that the liquid viscosity μ and the density difference between liquid and solid Δρ are constant and take k 0 as a reference permeability value.

Within the stagnant lid, an enriched, buoyant and eutectic concentration crust is formed through compaction and Darcy flow of the extracted melt. Conservation of eutectic composition melt therefore determines the evolution of crustal thickness during magma ocean solidification, with crustal magmatism driving enhanced fractionation of crustal material;

To account for potential sedimentation of dense cumulates during the initial, turbulent, magma ocean phase, we vary the slushy magma ocean depth, while we assume no incoming heat flux from its inner radius c, for simplicity.

Results

We find four distinct regimes in the evolution of the stagnant lid and stable buoyant crust, determined by the thermal and rheological properties of the slushy mixture (Figure 2) that are most clearly and simply demonstrated in absence of radiogenic heating and assuming a linear temperature profile in the lid (Equation 9).

Figure 2

An example of slushy magma ocean evolution. Evolution of stagnant lid and crustal thickness, crystal fraction and magma ocean temperature for μ  = 1018 Pa s, C  = 0.1, k 0 = 10−11 m2, neglecting radiogenic heating and assuming a linear temperature profile in the lid (bold lines) or solving for the full transient spherical temperature profile in the lid (thin lines). The four dynamical regimes are identified along with their characteristic asymptotic behaviors (dash‐dotted colored lines). Evolution of (a) stagnant lid thickness, (b) stagnant lid and crust thicknesses, (c) bulk crystal fraction and (d) magma ocean bulk temperature.

In the first regime, the interior properties remain approximately constant and a stagnant lid grows by conductive cooling. Using the linear thermal profile (Equation 9), to describe the evolution of the lid, Equation 8 gives (Figure 2, first regime). This regime lasts until the convective heat flux from the interior balances conduction across the lid and a quasi steady‐state lid thickness is reached (Figure 2, second regime) at time and where F (0) is the convective heat flux from the slushy interior with initial solid fraction and temperature, and .

During the second regime, which lasts from Ma, a significant fraction of the eutectic crust quickly forms by melt extraction at a large and approximately constant melt fraction () and interior temperature . The crust grows at a constant rate, which may be predicted by integrating the extracted eutectic component at the Darcy velocity,

In the third crystallization buffering regime, which lasts from Ma and which follows extraction of the eutectic component from the slushy mixture, the convective heat flux drives crystallization of the interior which increases the interior magma ocean viscosity (Equation 1) by increasing the solid fraction. Equilibration of the terms on the right of Equation (6) as the interior temperature remains approximately constant gives

As a result, the mean solid fraction evolves approximately as with a characteristic timescale τ cb = ρLV(b)3/[A(b)F (ϕ(0))E ] ≃ 1.5–8 Ma, where we use an initial solid fraction (see Figure 2). During this third regime, a balance between the conductive heat flux through the stagnant lid and the convective heat flux, , determines the stagnant lid thickness

The crustal thickness in this regime is determined by the melt extracted, and given by where we approximate a(t → 0) = 0 and assume (Figure 2b, see Text S1 in the Supporting Information S1 for more details). The bulk of the crust is formed during this regime. A timescale to form the bulk crust is given by the time to completely crystallize the interior, , which from Equation (17) is and is typically between Myr (400 Myr for the example in Figure 2). This result mainly depends on the slushy magma ocean depth, reference viscosity and bulk composition.

Finally, during the fourth shutdown regime, which lasts Myr, the interior temperature approaches the eutectic and the solid fraction , further increasing the viscosity (Equation 1). In this phase, the last eutectic liquid is extracted to form the crust and the crustal thickness is set, since any subsequent solid‐state convection of the lunar mantle does not lead to differentiation. Magma ocean cooling to the eutectic is described by which we may rewrite approximately as and where and is a Rayleigh number at the initiation of the stagnant lid, written here in terms of the equivalent magma ocean depth D accounting for sphericity, the rheological temperature difference (ΔT(0) ≃ 10–20 K) driving convection in the slushy magma ocean, (Figure S2 in the Supporting Information S1), and the solid reference viscosity μ which is decreased by a factor accounting for the presence of partial melt, (see Figure 3). In this final shutdown phase the lunar mantle temperature is described by

Figure 3

Crystallization time and final crustal thickness for different slush depths, reference viscosity, initial bulk anorthite content and reference permeability. Top and middle rows: the slushy magma ocean makes up the whole mantle, R  − c = 1,387 km, Top: k 0 = 10−11 m2, middle: μ  = 3 × 1018 Pa s. Bottom row: R  − c = 500 km, that is, the magma ocean is fully crystallized up to a radius c = 1,237 km, k 0 = 10−11 m−2. Time to reach the eutectic temperature in Ma (left) and final crust thickness in km (right) as a function of initial bulk anorthite content C and solidus viscosity μ or reference permeability k 0. We assume no radiogenic heating (Q = 0) and no incoming heat flux at the bottom of the magma ocean.

Since , this expression suggests a timescale for complete lunar highland formation of

The heat flux through the stagnant lid in this final phase balances the mantle heat flow, and provides a scaling for the evolution of the stagnant lid thickness with time in this last regime (see Figure 2 and Supporting Information S1).

These four cooling regimes are always present for a wide range of model parameters (Text S1 in the Supporting Information S1). Accounting for sphericity, and in particular secular cooling, in the lid thermal profile (solving thermal diffusion 7 rather than assuming a linear thermal profile 9) decreases the lid thickness (Figure 2a, thin solid lines) because the lid growth rate is larger than the rate at which heat is transported by conduction and heat remains stored in the deepest part of the lid. Consequently, the steady‐state lid thickness in regime two is not reached, and the lid is hotter and hence thinner which can further limit crustal growth during the crystal buffering regime (Figure 2b). Although the lid thickness is significantly reduced by solving for thermal diffusion (Equation 7), the evolution of the crustal thickness, magma ocean temperature and crystal fraction are not significantly changed, and the different cooling regimes are clearly visible on Figures 2b–2d. Since we neglect the heat released by crystallization in the lid as well as radiogenic heating, these timescales provide a lower bound on slushy magma ocean solidification. Including radiogenic heating in the model delays crust formation further, by 100s Ma, although radioelement extraction with the melt likely limits this effect (Figure S1 in the Supporting Information S1).

The solidification timescale of the slushy magma ocean of several hundreds of Myr is much longer than that of a liquid ocean overlaid by a floating anorthosite crust in particular because of the large, solid‐like viscosity of the slush. As expected, this timescale increases with the reference solidus viscosity and the initial depth of the slush (Figure 3). The solidification timescale also increases as the bulk composition becomes more depleted in anorthite‐component, since the magma ocean starts farther from the eutectic temperature (see phase diagram Figure 1). However, although crustal formation can be prolonged over a timescale that can even reach Gyr for large solidus viscosities, the bulk of the crust is formed during the crystallization buffering regime, that is in the first hundreds of Myr of evolution (Figure 2).

Mineralogical models of magma ocean crystallization from a bulk silicate composition suggest crystallization of ∼ 7%–10% plagioclase (Elkins‐Tanton et al., 2011; Snyder et al., 1992) from which we set our initial range of values for anorthite bulk composition, C  = [0.06–0.12], for a slush encompassing the whole lunar mantle (c = 350 km). If cumulates form during the first turbulent phase, the shallower slushy magma ocean becomes enriched in anorthite, thus to a much smaller slushy ocean depth of 500 km, we associate a range C  = [0.10–015] for anorthite conservation (Figure 3). In contrast to classical models of lunar magma ocean crystallization, the final crustal thickness not only depends on the initial composition and magma ocean depth but also depends on the rate of melt extraction and magma ocean cooling history, and hence on the characteristic permeability and viscosity (Figure 3). Our model can thus reconcile a relatively thin crustal thickness of ∼40 km with an initially fully melted lunar mantle hence a fully liquid Moon (Pritchard & Stevenson, 2000) since a non‐negligible fraction of anorthite‐component can remain solid in the mantle while classical models produce an overabundance of plagioclase from a fully melted mantle (Charlier et al., 2018; Elkins‐Tanton et al., 2011; Lin et al., 2016).

Discussion

The predicted timescale of slushy magma ocean solidification is in good agreement with the estimated Myr age range of lunar FAN, suggesting that a slushy magma ocean is not incompatible with an ancient age for the Moon.

Solidification and differentiation of a slushy magma ocean results in the formation of a eutectic‐concentration crust that contains of anorthite, while the anorthite concentration of the Highlands may be closer to 90%–95% as suggested by Apollo anorthosite samples (Tompkins & Pieters, 1999). Such an extreme enrichment in anorthite remains challenging to explain. Even in scenarios envisioning a flotation crust, the extraction of sufficient residual melt from between anorthite crystals to such low values is challenging (Piskorz & Stevenson, 2014). In our present model, melt extraction is simplified, and in practice crustal differentiation should lead to the formation of magma reservoirs in the deep crust by accumulation of eutectic liquids enriched in radioactive incompatible elements. Differentiation of these liquids would produce buoyant anorthite‐rich diapirs which could lose their mafic residue by deformation and squeezing while ascending through the crust (Jolliff & Haskin, 1995; Longhi, 2003), leaving behind different types of residual mafic liquids potentially enriched in Magnesium, as the Mg‐suite rocks, and in incompatible elements.

Since, in the slushy model, anorthosite would be formed by the extraction of evolving mantle melts, this process would likely result in various types of anorthosites. Anorthosite samples are indeed heterogeneous (Nyquist et al., 2006; Shearer & Floss, 2000) which is compatible with a magmatic origin (Russell et al., 2014). Their chemistry and petrology seems to require a two‐stage process formation (Jolliff & Haskin, 1995; Longhi & Ashwal, 1985). Furthermore, recent analysis of lunar highland meteorites, that may represent a more random sampling of the lunar surface than Apollo samples, shows that crustal anorthosites may be more magnesian and similar in composition to terrestrial massif anorthosites formed by magmatism (Gross et al., 2014). Serial magmatism has been proposed by many authors as a viable way for lunar crust formation for several decades. The slushy lunar model reconciles this view with the magma ocean scenario explaining the complexities observed in lunar anorthosite samples (Gross et al., 2014; Haskin & Warren, 1991; Jolliff & Haskin, 1995; Longhi & Ashwal, 1985; Nyquist et al., 2006; Pritchard & Stevenson, 2000; Russell et al., 2014; Shearer & Floss, 2000).

Supporting Information S1

Click here for additional data file.