## Abstract

The lack of contraction features on the Moon has been used to argue that the Moon underwent limited secular cooling, and thus had a relatively cool initial state. A cool early state in turn limits the depth of the lunar magma ocean. Recent GRAIL gravity measurements, however, suggest that dikes were emplaced in the lower crust, requiring global lunar expansion. Starting from the magma ocean state, we show that solidification of the lunar magma ocean would most likely result in expansion of the young lunar crust, and that viscous relaxation of the crust would prevent early tectonic features of contraction or expansion from being recorded permanently. The most likely process for creating the expansion recorded by the dikes is melting during cumulate overturn of the newly solidified lunar mantle.

## 1. Introduction

Since the first Apollo samples were returned in 1969, the Moon has been posited to have begun in at least a partially molten state [1,2]. The Moon's surface, however, does not show large-scale global features created by contraction [3,4]. Earlier geodynamic models of lunar contraction suggest that the lunar magma ocean was constrained in depth, because freezing of a deep magma ocean would accumulate more crustal strain than is observed on the Moon [5,6]. More recent work by Zhang *et al.* [7] indicated that magma ocean depth was not a constraint on contraction later in lunar history during thermal evolution following magma ocean solidification.

Gravity data from the GRAIL mission have recently been interpreted as showing dike intrusions into the lunar lower crust, which require 0.6–4.9 km expansion of the Moon during crustal formation [8]. Here, we test whether the Moon would be expected to expand or contract during its earliest history, and whether these dikes may have been created during that time.

The lack of contraction features on the Moon has been used to argue that the Moon underwent limited secular cooling, and thus had a relatively cool initial state, indicating in turn that the lunar magma ocean had limited depth. Extensive contraction during cooling from a hot initial state had been expected to result in thrust faults and scarps, as on Mercury.

The net change in planet volume owing to freezing the magma ocean is
1.1
where *V*_{l} is the initial melt volume. The associated surface strain would then simply be
1.2

The solidifying Moon is not, however, a simple bulk system that transfers its entire volumetric change directly to surface strain. The critical difference is that while crystallization of mafic phases (whose solids are usually denser than their coexisting liquid) results in contraction, crystallization of plagioclase (which is less dense than the coexisting liquid) is associated with expansion. Thus, freezing of an assemblage with sufficient plagioclase will result in expansion during solidification (see also Kirk & Stevenson [9], who invoke expansion during a late gabbroic melting episode). Thus, two significant modifications have to be made to the simple assessments of lunar contraction or expansion during freezing:

The relative densities of solid and liquid are sensitive to liquid composition, particularly TiO

_{2}, FeO and SiO_{2}contents. Therefore, tracking the composition of the evolving magma ocean liquid and the corresponding equilibrium composition of the solidifying minerals over a number of solidification steps is critical to calculate contraction or expansion of the body.The plagioclase flotation crust, in which the strains might be recorded, is not formed until after about 80% of the magma ocean is solidified. Before that point, the solids are overlain by a liquid magma layer and will not record strain.

Here, we present compositional models of the freezing magma ocean that indicate either contraction or expansion is possible in the flotation crust. Further, we derive wavelengths and timescales for the folding that would predate faulting in a contracting lunar flotation crust, and we show that relaxation would occur so rapidly that any sign of contraction would likely be lost.

## 2. Models of solidification

### (a) Density evolution during magma ocean freezing

We consider three lunar magma ocean scenarios. In two cases (labelled LMO1 and LMO3), the magma ocean has the same bulk composition as was used in Elkins-Tanton *et al.* [10], but with slightly different and equally plausible mineral assemblages during solidification. In the third case, LMO2, the lunar primitive upper mantle composition from Longhi [11] is used.

The solidifying mineral assemblage and the point at which the flotation crust forms are critical to analysing the possible crustal strain that resulted from solidification, and they depend upon both the bulk composition of the lunar magma ocean and whether the solidification process is fractional, bulk or a mixture. Here, we assume pure fractional solidification produced by ongoing separation owing to settling or ascent of crystals. Unfortunately, bulk magma ocean compositions are unknown for any planet, and a fairly wide range of compositions have been proposed for the Moon (tables 1 and 2).

The pressure boundaries and mineral assemblages predicted for a fractionally solidifying lunar magma ocean are broadly consistent among researchers. There are currently no complete experimental fractional crystallization studies of relevant lunar compositions on which to base these mineral assemblages. Previous studies, including those on which we base our analysis, used theoretical calculations of phase stability [12,13,15–18,20–22].

Previous workers [10] investigated the constraints a given bulk composition exerts on the solidifying phase assemblages. They found that in the absence of tightly constraining phase equilibrium data, there is significant leeway—without violating bounds on the bulk magma ocean composition—in both the ratios of phases within a solidifying assemblage (e.g. 40–60% plagioclase is allowed when the magma ocean is between 0.80 and 0.92 solidified) and with the stage at which the solidifying assemblage changes (e.g. plagioclase begins to crystallize when 70% to more than 80% of the magma ocean has solidified).

In our model, the mineral assemblages allow solidification to proceed to 92% before expending any of the major oxides. At high degrees of solidification, the standard mineral assemblages shown in table 2 would no longer pertain, and the liquid would approach a highly evolved KREEP composition enriched in incompatible elements. No satisfactory model for the solidifying assemblage of a very late-stage liquid exists, thus the final 8% is assumed to solidify in bulk.

In each case, the magma ocean is assumed to be 1000 km deep, and to retain 1% interstitial liquid at each step of solidification which proceeds in equal volume increments *δV* ; as shown by Elkins-Tanton *et al.* [10], values of *δV* =*V*_{l}/1000 (where recall *V*_{l} is the total initial magma ocean volume) provide sufficient model convergence. Other details of the solidification model are described fully in Elkins-Tanton *et al.* [10]. At each step of solidification, the fractional change in volume produced by solidification is calculated by the ratio of liquid density to solid density. The liquid density is calculated first at as a function of composition, temperature and pressure. Composition and temperature dependence are from Lange & Carmichael [23]. Pressure dependence, which includes only the first derivative of *V* with respect to temperature *T*, is from [24]. Further updates are from Liu & Lange [25] and from Ochs & Lange [26,27]. These calculations are expected to be accurate only to about 2 GPa, so here we use them for 1 atm calculations only, and continue the density to depth using a third-order Birch–Murnaghan equation of state, after the technique of Delano [28].

To calculate the density of solids, first their composition is calculated in equilibrium with the coexisting magma ocean composition using experimental distribution coefficients [29,30], and then densities are calculated based on their temperature, pressure and composition (databases are given in [30]). The molar volume of the mineral phase is at one atmosphere, and the model temperature is assumed to be the solidus temperature at the depth of solidification. The pressure effect on density is then calculated separately for each phase using appropriate moduli and a third-order Burch–Murnaghan equation of state, which is solved iteratively for the molar volume at a given pressure and temperature by adjusting the pressure *P* from the Burch–Murnaghan equation until it matches the model pressure. For details of the calculation, see Elkins-Tanton *et al.* [29] and for updated thermophysical parameters, see the electronic supplementary material for Elkins-Tanton [30].

In each of the models, the evolving liquid becomes richer in iron and titanium as solidification proceeds, and the densities of the resulting mineral assemblages also increase (figure 1). The slight differences in the solidifying mineral assemblages used for models LMO1 and LMO3 cause the solid densities to diverge towards the end of solidification, even though the models have the same bulk composition. This divergence shows the sensitivity of the system to parameters whose true values we do not know.

The final liquid densities (figure 1) are not simply a product of iron enrichment (table 2). The final liquid for LMO3 is far denser than the other model liquids, in part because its high iron is not balanced with the high quantities of low-density alumina in the LMO1 liquid, nor with the higher MgO in LMO2 liquid. Note also that the smaller fraction of oxides being removed from LMO2 has produced high chromium in that liquid, also affecting its density. Thus, the density (volume) of the fractionating solids is not always constant with respect to the density (volume) of the liquids, because their relationship depends upon the partitioning of each element between the liquid and solid. Thus, a simple Fe–Mg or Fe–Si model will not make an acceptable prediction for the volume change of a solidifying Moon.

The late-stage liquids shown in table 3 are also geochemically pertinent. These liquids represent the candidate predecessors to KREEP, the last, incompatible-element-enriched melts from the magma ocean that appear in trace amounts in some lunar samples [31]. The liquids in table 3 are the final 8% of the magma ocean, and KREEP itself is likely to be a later-stage fractionate, but its original composition can only be modelled from samples. The liquids shown here, however, are sufficiently different from each other that constraints from natural samples may be applied to models using these liquids.

Figure 2 shows the cumulative radius change as the body solidifies. The models show planetary contraction of about 13 km during the olivine-only lowest layer of solidification. During the second layer of solidification, consisting of olivine and orthopyroxene, models LMO1 and LMO2 indicate continued contraction of an additional approximately 7 km, whereas model LMO3 shows expansion of about a kilometre. None of these radius changes would be recorded in the lunar crust, however, because no solid lid has formed on top of the magma ocean, so the liquid uppermost layer simply conforms to whatever contraction or expansion is happening beneath it.

A solid lid prior to plagioclase flotation is unlikely on the Moon. Even a thin nebular atmosphere may hold the surface above its solidus [32] and prevent formation of a quench crust (thus, no quench crust would be expected on larger planets such as the Moon or the Earth). Should surface temperatures exist that allow a quench crust, this crust would be denser than its underlying magma and would rapidly founder [33,34]. A possible quench crust on the Moon is critically different from a lava lake on the Earth: the lava crust clings to the solid sides of the lake, but on the Moon, its density would cause it to sink. Thus, on the Moon, the first solid lid would be the plagioclase flotation crust.

All these models assume that plagioclase begins to form at 80% solidification. We further assume that 5 km of plagioclase flotation is required to make a competent lid (as was assumed in [10]), and this requires about 2% additional solidification beyond the 80% point (this depends, of course, upon the fraction of plagioclase in the solidifying assemblage, and assumes perfect and instantaneous flotation; the choice of 5 km as a competent lid is arbitrary). This is the critical stage, after the plagioclase lid forms and can record strain. Models LMO2 and LMO3 show approximately 4 km expansion over the remaining solidification, whereas LMO1 shows just a hundred metres of further contraction. Models LMO1 and LMO3 produce approximately 46 km of flotation crust, and LMO2 produces approximately 48 km. These estimates are somewhat higher than the current best average crustal thickness estimates, between 34 and 43 km [35].

The differences in volume change over the last 20% of solidification among the models are produced in part by slight differences in mineralogy, but in larger part by the different relative densities of the late-stage liquids. In some cases, even pyroxenes become less dense than the coexisting liquid, which is of possible significance, because there is an enigmatic pyroxene fraction in much of the anorthosite flotation lid on the Moon [36]. The lunar magma ocean models shown here indicate that pyroxene may in some cases naturally float into the plagioclase lid as it forms, and not require any more complex explanation.

We conclude that expansion or contraction of the solidifying Moon was sensitive to both bulk composition and solidifying assemblages, to a degree that we cannot address because these compositions are not sufficiently constrained for any planet. Further, the expansion or contraction that happens after a solid crust is formed is the only expansion or contraction that will be recorded, and in the case of magma ocean processes, this encompasses only the last 20–30% of solidification, and may be too small to be recorded.

Gravity data from the GRAIL mission have recently been interpreted as showing dike intrusions into the lunar lower crust, which would require 0.6–4.9 km expansion of the Moon during crustal formation [8], consistent with the models presented here. The amount of contraction or expansion is not a constraint on the original depth of the magma ocean; as demonstrated here, a wide range of expansion or contraction values can be obtained for a single magma ocean depth.

These models demonstrate that either contraction or expansion is possible within the range of suggested bulk compositions and solidifying assemblages. If, in the limiting case that the Moon experienced contraction, would thrust faults be expected? In §3, we derive the growth rate and wavelength of the folding that would result from contraction and precede thrust faulting, and discuss whether any vestige of this process would be preserved for detection in the present day.

## 3. Crustal folding and relaxation

### (a) Folding during contraction

As a magma ocean freezes, its floating crust could conceivably shorten and undergo folding during planetary contraction, and these folding instabilities would lead to features such as thrusts and scarps. Because the crust is still hot and accumulating, we treat it as a viscous layer, bounded above and below by relatively inviscid medium (vacuum or tenuous atmosphere above, magma below). The crustal layer is of thickness *h*, and we allow the origin plane *z*=0 to be at the mid-depth, so the top and bottom are at *z*=±*h*/2. If the layer is deflected vertically during folding at some horizontal position *x* by a vertical displacement *w*, then the force balance on a thin vertical segment of the layer implies (removing the hydrostatic balance)
3.1
where *σ*_{xz} is the shear-stress on the vertical sides of the segment, *ρ*_{m} is the magma-ocean density, *ρ*_{a} is the ‘atmospheric’ density. We have assumed there is no load to the crustal layer other than its own deflection (i.e. no volcanic masses). The net torque on this segment is
3.2
where
3.3
where the term in [..] is the hydrostatic pressure at some depth *z* inside the crust, *P*_{a} is the atmospheric pressure at an undeflected surface *z*=*h*/2, *ρ*_{c} is crustal density, *τ*_{0} is a uniform imposed normal stress (in the *x*-direction owing to contraction) and the last term on the right-hand side of (3.3) is the non-uniform deviatoric stress, owing to curvature of the layer [37], in which *μ* is crustal viscosity. After carrying through the integration in (3.2) on *σ*_{xx} as given by (3.3), the torque equation becomes (to first order in *w*)
3.4
where
3.5
Substituting (3.4) into (3.1), we arrive at the governing equation for the deflection *w*
3.6
where Δ*ρ*=*ρ*_{m}−*ρ*_{a}.

#### (i) Non-dimensional equations and folding growth rate

Non-dimensionalizing space and time according to
3.7
(3.6) becomes (after dropping the primes)
3.8
Using normal mode analysis, we assume *w*∼*e*^{ikx+st} in that case (3.8) yields a simple relation for the growth rate *s*
3.9
The growth rate is of course negative for *k*<1 and goes to 0 for . But *s* is maximum at wavenumber at which the growth rate is .

Thus, the dimensional wavelength for the fastest growing folding mode is 3.10 and the minimum growth time for folding (i.e. inverse of the growth rate for the fastest growing folding mode) is 3.11

#### (ii) Contraction stresses

The wavelength and time-scale for folding instability (equations (3.10) and (3.11)) depend on the values of *Π* and hence the choice for *τ*_{0}. There are two end-member values for *τ*_{0}, which would be (i) the compressive stresses are due to the magma ocean contracting faster than the crust can shorten and thicken, and (ii) the stresses occur as the shortening/thickening keeps pace with magma ocean contraction. In the first case, we can estimate the lateral compressive stress by considering a cylindrical axisymmetric geometry and how the hoop stress relates to radial or vertical load; in this case, the momentum equation in the radial direction is
3.12
where *r* is radial distance from the centre of the planet, and we have ignored shear stresses, because motion is only radial for net contraction. (Because the folding theory is in two dimensions, this analysis is best compared with an axisymmetric geometry, but the result would be the same if we used a spherically axisymmetric geometry.) Because the crustal displacement *w* is independent of *z* and hence *r*, then *τ*_{rr}=0, and hence the lateral contraction stress is
3.13
If the system is hydrostatic then *τ*_{0}=0.

Condition 1 is if there is no support from hydrostatic pressure from below because contraction is too fast, and then we obtain *τ*_{0}=−*ρ*_{c}*gr*, which is the largest possible contraction stress. In this case, *Π*=*P*_{a}+*ρ*_{c}*g*(*h*/2+*r*)≈*ρ*_{c}*gr*, because *P*_{a}≈0 given the Moon's lack of atmosphere, and *h*≪*r*. Then, the wavelength and timescale for the least stable folding mode are
3.14
respectively. Using *ρ*_{c}≈2700 kg/m^{3}, Δ*ρ*/*ρ*_{c}≈1, *g*=1.6 m/s^{2}, *r*=1740 km, *h*=20−100 km, then *λ*_{m}≈800−2000 km and years.

Condition 2 is if contraction of the crust keeps pace with contraction of the magma ocean. In this case, where is given by the contraction rate. The contraction rate is determined by the fact that the net planetary mass is conserved, where is the planet's mean density, which changes during freezing and contraction (see §2a); we thus obtain , and therefore . In this case, . Using *ρ*_{c}=2700 kg m^{−3}, *g*=1.6 m s^{−2} and *h*=100 km, *μ*≈10^{18} Pa s, then *τ*_{0} is only comparable to *ρ*_{c}*gh*/2 if , i.e. if the contraction timescale is shorter than about 100 years, which is unlikely; thus, we assume and thus, *Π*≈*ρ*_{c}*gh*/2. In this case, the wavelength and timescale for the least-stable folding mode are
3.15
respectively. Using the same parameters as above then *λ*_{m}≈60−300 km, and .

Although the folding times for either case are short, it is important to compare them with the relaxation time over which the warped crust flattens viscously.

#### (iii) Crustal viscous relaxation

The viscous folding theory assumes that the crustal layer is folded without changing its net thickness *h*. However, gravitational sliding would cause this thickness to adjust spatially in order to erase the deflection *w*. A complete folding theory would allow for variations in *h*, however, that is beyond the scope of this simple analysis. Nevertheless, we can estimate the relaxation time for *w* given pressure gradients and gravitational sliding owing to the folding instability.

The balance of forces parallel to the crustal surface on a segment of crust that is d*x* wide and *h* thick is
3.16
where is the stress tensor (and assuming stress-free top and bottom), and is the unit vector parallel to the surface given by *z*=*w*(*x*), i.e. (assuming ). Using then to first order in *w*, we arrive at
3.17
where *v*_{x} is the horizontal velocity of a segment of crust (which we assume is independent of *z* given the free-slip top and bottom). The hydrostatic pressure at position (*x*,*z*) is approximately
3.18
where *P*_{0} is determined by the boundary condition that *P*(*x*,*h*/2+*w*)=*P*_{a}−*ρ*_{a}*gw*, where *P*_{a} is the atmospheric pressure at *z*=*h*/2. (Equation (3.18) comes about by saying that at the surface *P*(*x*,*h*/2+*w*)=*P*_{a}−*ρ*_{a}*gw*, and at the base of the crust *P*(*x*,−*h*/2+*w*)=*P*_{a}+*ρ*_{c}*gh*−*ρ*_{m}*gw*, which assumes that the depth *z*=−*h*/2 is a constant pressure level in that the magma ocean always equilibrates pressure below this depth). In this case, (3.17) becomes
3.19
where . Integrating in *x* implies
3.20
where *τ*_{0} is the background in-plane stress. As the crustal deflection relaxes viscously, conservation of mass requires that
3.21
where the first term on the right-hand side represents the decay of the deflection, whereas the second term represents the uniform thickening owing to contraction stresses *τ*_{0}<0, as discussed above. The relaxation timescale is thus
3.22
Given the above parameter values this timescale is about 300–1500 years, which is only a bit larger than the longest growth time for the folding instability, and is also very short. This implies the folding instability is fairly rapid and is established quickly, but then relaxes away almost as quickly. Thus, given the short timescales for folding and relaxation, it is very unlikely that compressional features would persist in a floating crust on a contracting lunar magma ocean.

## 4. Discussion and conclusions

Based on both extensive petrologic evidence and modelling, the Moon accreted with a magma ocean of some depth, and the magma solidified via fractional crystallization, differentiating its silicates into a heterogeneous cumulate mantle, late-stage incompatible element-enriched materials, and an anorthositic flotation crust [1,2,38–43].

No crust exists to record strain owing to volumetric change within the Moon until the anorthositic crust is developed to some thickness, which is not likely to be obtained until after at least 80% solidification by volume. In some cases, during the last stages of magma ocean solidification, both pyroxenes and plagioclase become less dense than the coexisting liquid. The simultaneous flotation of both phases may explain the enigmatic pyroxene fraction in much of the anorthosite flotation lid on the Moon [36].

Expansion or contraction during solidification is controlled sensitively by the thermochemical parameters of the solidifying phases compared with the liquid; the underlying solid mantle will remain at its solidus temperature and not thermally contract until solidification is complete and the conductive cooling front reaches the mantle. Unfortunately, bulk silicate composition and solidifying assemblages are not sufficiently well known on any planet to make accurate forward predictions of expansion or contraction. The analysis of magma ocean solidification shown here, however, favours expansion over contraction. If the lunar crust did contract, then the folding timescale and relaxation timescale are both very short, on the order of hundreds to thousands of years, and thus any folds caused by contraction would relax away and not appear in the present day as either faults or folds.

Andrews-Hanna *et al.* [8] suggest that intrusions in the lunar lower crust imaged with GRAIL gravity data require 0.6–4.9 km expansion of the Moon during crustal formation, consistent with the models presented here. The amount of contraction or expansion is not a constraint on the original depth of the magma ocean: as demonstrated here, a wide range of expansion or contraction values can be obtained for a single magma ocean depth. These dikes were not likely formed at the end of magma ocean solidification, however. Although the expected expansion of the Moon at this stage matches the calculations from Andrews-Hanna *et al.* [8], the lower crust was almost certainly too soft to allow the immediately underlying liquids to form brittle dike structures. The temperature of solidification of the late magma ocean liquids has been estimated at above 1100°C [44], and this would also be the temperature of the bottom of the anorthosite lid at that time.

Following the bulk of cumulate formation, the lunar mantle would be unstable to gravitational overturn and is expected to have reordered in the solid state over a period of a few million years [44]. During this reordering, the flotation lid would have been cooling and becoming more rigid as dense late-stage cumulates sank away beneath it, but hot buoyant orthopyroxene + olivine cumulates from the lower mantle are expected to have risen and melted adiabatically. Melting parts of the orthopyroxene + olivine layer is expected to produce less than 5 km of global expansion (figure 2). This overturn and melting event would thus be expected to produce global expansion and a dense mafic melt commensurate with the dikes inferred with Andrews-Hanna *et al.* [8], and to occur at a time when the flotation crust may be rigid enough to respond to these events by forming linear dikes. Further, the lunar upper mantle would then consist of materials with a strong orthopyroxene (not olivine) mode, consistent with the observations of Melosh *et al.* [45].

Although cumulate overturn offers a process consistent with expansion, dike formation and the compositions of the dikes and the upper mantle, it is also possible that the expansion required for dike formation may have occurred at the end of magma ocean solidification, or it could have occurred through later heating within cumulates enriched in radiogenic elements [7]. Further evidence for crustal tectonism during magma ocean freezing and overturn, as well as the early stages of radiogenic heating, was likely eliminated by the late heavy bombardment. Thus, the evidence gleaned from compositions within and without impact basins and gravity must form the tests for theory. Dike emplacement in the lower crust likely coincided with near post-magma ocean processes of cumulate overturn or radiogenic heating, and not with later processes where the requisite expansion is not predicted.

## Funding statement

Support was provided by the DTM Merle A. Tuve Senior Fellowship, and by NSF Astronomy CAREER grant no. 0747154.

## Footnotes

One contribution of 19 to a Discussion Meeting Issue ‘Origin of the Moon’.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.