Rapid solidification of Earth's magma ocean limits early lunar recession

The early evolution of the Earth-Moon system prescribes the tidal environment of the Hadean Earth and holds the key to the formation mechanism of the Moon and its thermal evolution. Estimating its early state by backtracking from the present, however, suffers from substantial uncertainties associated with ocean tides. Tidal evolution during the solidification of Earth's magma ocean, on the other hand, has the potential to provide robust constraints on the Earth-Moon system before the appearance of a water ocean. Here we show that energy dissipation in a solidifying magma ocean results in considerably more limited lunar recession than previously thought, and that the Moon was probably still at the distance of $\sim$7-9 Earth radii at the end of solidification. This limited early recession aggravates the often overlooked difficulty of modeling tidal dissipation in Earth's first billion years, but it also offers a new possibility of resolving the lunar inclination problem by allowing the operation of multiple excitation mechanisms.


Introduction
Since the pioneering work of George Darwin (Darwin, 1879(Darwin, , 1880, the tidal evolution of the Earth-Moon system has been studied by a number of scientists (MacDonald, 1964;Kaula, 1964;Goldreich, 1966;Lambeck, 1975;Webb, 1982;Touma and Wisdom, 1994;Bills and Ray, 1999;Abe and Ooe, 2001;Zahnle et al., 2015;Tian and Wisdom, 2020;Farhat et al., 2022). Tidal friction caused by lunar tides slows down Earth's rotation and results in the gradual recession of the Moon. It follows that Earth was rotating faster in the past and the lunar semimajor axis was shorter than its present-day value, ∼60 Earth radii (R E ). Estimates on Earth's past rotation rate from tidal deposits are consistent with such orbital evolution (Williams, 2000;Coughenour et al., 2009), with the oldest inference being 53.5±0.4 R E at 1.4 Ga (Meyers and Malinverno, 2018). Going deeper in time with theory is subject to considerable uncertainties. Ocean tides, which are overwhelmingly more important than earth tides (Lambeck, 1975;Webb, 1982), are sensitive to rotation rate, ocean depth, seafloor topography, and continental configuration (Webb, 1982;Abe and Ooe, 2001;Motoyama et al., 2020), and the Hadean and the Archean are the period when continents may have grown rapidly (Harrison, 2009;Rosas and Korenaga, 2018;Guo and Korenaga, 2020), and the ocean likely experienced non-monotonic volume changes (Korenaga, 2021a;Miyazaki and Korenaga, 2022). Whereas the treatment of ocean tides in recent studies on the Earth-Moon evolution has become sophisticated (Tyler, 2021;Daher et al., 2021;Farhat et al., 2022), the effects of changing ocean volume as well as changing continental mass are yet to be considered, so the reliability of such recent results is unclear for the Hadean and the Archean.
The early evolution of the Earth-Moon distance plays a central role in various mechanisms proposed for the origin of the lunar inclination, which in turn is critical to the formation mechanism of the Moon itself (Goldreich, 1966;Touma and Wisdom, 1998;Ward and Canup, 2000;Pahlevan and Morbidelli, 2015;Chen and Nimmo, 2016;Ćuk et al., 2016). Furthermore, the lunar distance controls the magnitude of tidal dissipation within the Moon, thereby affecting its early thermal evolution, such as the duration of the lunar magma ocean (Elkins-Tanton et al., 2011). It also determines the amplitude of ocean tides on Earth, which has important ramifications for the origin of life. Some of major abiogenesis hypotheses require the presence of exposed landmasses, on which wet-dry cycles can promote the polymerization of organic compounds (Campblel et al., 2019;Damer and Deamer, 2020), but the stable operation of wet-dry cycles would be jeopardized if exposed landmasses are frequently swamped by huge tides (Russell and Hall, 2006).
Complications associated with ocean tides can be avoided if we track the early tidal evolution forward in time, starting from the aftermath of the Moon-forming giant impact. Though the de-tails of the Moon-forming giant impact are still actively debated (Canup, 2012;Ćuk and Stewart, 2012;Asphaug, 2014;Rufu et al., 2017), the Moon is likely to have formed just outside the Roche limit, at the distance of ∼3-5 R E (Salmon and Canup, 2012;Ćuk and Stewart, 2012). Tidal dissipation is limited in a completely molten Earth, but it can be substantial in a solidifying magma ocean, allowing rapid lunar recession. An order of magnitude calculation based on an isothermal magma ocean suggests that, while a magma ocean solidifies within one hundred million years after the giant impact, the lunar distance may expand to ∼40 R E (Zahnle et al., 2015). The isothermal approximation is, however, unlikely to represent some important aspects of the physics of a solidifying magma ocean (Abe, 1993;Solomatov and Stevenson, 1993;Abe, 1997;Solomatov, 2015).
To constrain the early history of the Earth-Moon system, it becomes essential to use a more realistic model of Earth's magma ocean. In this paper, we will show that, in contrast to the earlier suggestion by Zahnle et al. (2015), lunar recession during magma ocean solidification is likely to have been severely limited, reaching only ∼7-9 Earth radii, which has important implications for the origin of the lunar inclination. The structure of this paper is as follows. In the Methods section, we first describe how we parameterize the mantle phase diagram and how we calculate mantle adiabats. Mantle rheology is then described, and the timescale of Rayleigh-Taylor instability is estimated, with which the use of adiabats can be justified for tidal dissipation calculation. With these preparations, we describe our procedure to calculate tidal dissipation, thermal evolution, and orbital evolution. The results of our orbital calculation are then presented, followed by the discussion of their implications. The main focus of this paper is to lay out how to take into account the physics of magma ocean solidification into the tidal dissipation of Earth, and as such, the other aspects of tidal evolution are kept simple.

Phase diagram
Based on the experimental and theoretical studies of high-pressure mantle melting (Zhang and Herzberg, 1994;Fiquet et al., 2010;Miyazaki and Korenaga, 2019a), the mantle phase diagram is parameterized using cubic splines as follows: where T s and T l are solidus and liquidus temperatures, respectively, in Kelvin, z is depth in km, H() is the Heaviside step function, and S() is the cubic spline function. The spline function is defined as where h i = z i+1 − z i , z i is the ith knot, and T i and b i are the values of the function and its second derivative at the ith knot. We use four knots, located at the depths of 0, 410, 670, and We also define T 40 as the temperature at which the melt fraction is 40%, and it is parameterized as where f 40 (z) is linearly interpolated from 0.4 at 0 km, 0.4 at 410 km, 0.6 at 670 km, and 0.845 at 2900 km. That is, the melt fraction is more nonlinearly distributed between the solidus and the liquidus in the lower mantle (Miyazaki and Korenaga, 2019b). This is our standard choice of the mantle phase diagram and is referred to as 'MK19'. For comparison, we also use the version called 'const', in which f 40 is 0.4 throughout the mantle. At a given depth, the melt fraction between T s and T 40 is assumed to be distributed linearly from 0 to 0.4, and that between T 40 and T l is also linearly from 0.4 to 1.0.

Mantle adiabat
When temperature is above the liquidus or below the solidus, a mantle adiabat for a given surface temperature is calculated simply by integrating the isentropic thermal gradient where T is temperature, P is pressure, α is thermal expansivity (K −1 ), ρ is density (kg m −3 ), and C p is isobaric specific heat (J kg −1 K −1 ). Based on the thermodynamic calculations conducted by Miyazaki and Korenaga (2019b), these material properties are parameterized as α(T, P) = 3.622 × 10 −5 exp(−2.377 × 10 −5 T − 0.0106P), ρ(T, P) = 2870 − 0.082T + 162P 0.58 , where T is in Kelvin and P is in GPa.
When temperature is between the solidus and the liquidus, the effect of the latent heat of fusion must be taken into account as where H f is the latent heat of fusion and φ is the melt fraction. The term (∂ φ /∂ P) S can be obtained by requiring self-consistency with φ (T, P) given by the chosen mantle phase diagram (Langmuir et al., 1992). To approximate the adiabat calculations of Miyazaki and Korenaga (2019b), which are based on Gibbs energy minimization, the latent heat of fusion is linearly interpolated from 6 × 10 5 J kg −1 at 0 GPa and 9 × 10 6 J kg −1 at 136 GPa.
As seen in Figure 1a, upon cooling, the mantle adiabat crosses the liquidus at the base of the mantle. Crossing the liquidus first at some mid-mantle depth, i.e., the formation of a basal magma ocean, has been suggested on the basis of the thermodynamics of simple silicate melt (Labrosse et al., 2007;Mosenfelder et al., 2009;Stixrude et al., 2009), but its possibility appears to be small for mantle-like melt compositions (Miyazaki and Korenaga, 2019b).

Rheology
The rheology of a magma ocean is a strong function of the melt fraction, and when φ > φ c , where φ c is the critical melt fraction, we simply assume that viscosity takes the value of melt viscosity, η L .
The viscosity of a solid-melt mixture increases by a few orders of magnitude when φ approaches φ c from the above (Abe, 1993;Solomatov, 2015), but such an increase is considerably smaller than a viscosity jump associated with rheological transition at φ = φ c . Our simplified treatment does not affect our calculation of tidal dissipation, in which a layer with φ > φ c is treated as a fluid. We set φ c to the traditional value of 0.4 (Abe, 1993;Solomatov, 2015) for our main results.
An experimental study with olivine aggregates suggests φ c ∼ 0.3 (Scott and Kohlstedt, 2006), and experiments with other silicate systems suggest that φ c may vary from 0.2 to 0.6 (Costa et al., 2009). As shown later, however, the effect of varying φ c on lunar recession is minor. Melt viscosity is 0.1 Pa s for peridotitic melt composition (Dingwell et al., 2004); this is most appropriate for Earth's magma ocean, so it is used in the reference case. At the late stage of magma ocean solidification, higher melt viscosity corresponding to basaltic melt composition (5-10 Pa s (McBirney, 1993)) may become more appropriate.
When φ <= φ c , the viscosity of a solid-melt mixture (or just a solid when φ = 0) is calculated where η s (z) is the background solid viscosity at the solidus, E is the activation energy, and R is the universal gas constant. Experimental studies suggest α η ∼ 26 for diffusion creep and α η ∼ 31 −32 for dislocation creep (Mei et al., 2002;Scott and Kohlstedt, 2006). For the reference case, α η is set to 26. The activation energy is set to 300 kJ mol −1 (Korenaga, 2006). For the reference case, log 10 η s is linearly interpolated from 19 at the top and 20 at the bottom of the mantle; this background viscosity is based on an extrapolation from the present-day viscosity structure (Miyazaki and Korenaga, 2019b). For softer background viscosity, we interpolate log 10 η s from 18 and 19, and for harder background viscosity, from 20 to 21. The viscosity of the present-day convecting mantle is lowered by its water content (Hirth and Kohlstedt, 1996), and the solid phase coexisting with the melt phase is dry because water is mostly partitioned into the melt phase. This consideration is absent in the extrapolation done by Miyazaki and Korenaga (2019b), so the harder background viscosity may be more appropriate (Karato, 1986).
For viscoelastic rheology used in the calculation of tidal dissipation, we use the temperaturedependent model of Takei (2017), which offers a more complete treatment of viscoelasticity compared to earlier attempts (Faul and Jackson, 2005;Jackson and Faul, 2010). In this model, the complex compliance is parameterized as with and where ω is the tide-raising frequency, J U is the unrelaxed compliance (= 1/µ E , where µ E is the elastic shear modulus), and p ′ is the period normalized by the Maxwell relaxation time (τ M = η/µ E ), i.e., p ′ = (ωτ M ) −1 , and other parameters are defined as in Yamauchi and Takei (2016) and Takei (2017). The elastic shear modulus is calculated as the Voigt-Reuss-Hill average of the shear moduli of the solid and melt phases. The shear modulus of the solid phase is calculated as a function of depth from the Preliminary Reference Earth Model (PREM (Dziewonski and Anderson, 1981)), and that of the melt phase (at the infinite frequency limit) is assigned a value of 10 GPa (Webb and Dingwell, 1995). The bulk viscosity is practically zero (Schubert et al., 2001), so we assume that the bulk modulus remains purely elastic. The complex shear modulus,μ, is simply the reciprocal of the complex compliance.
The temperature-dependent model of Takei (2017) may be considered as a modified Andrade model, characterized by a broad relaxation spectrum. For comparison, we also consider a Maxwell body, for which the complex shear modulus is calculated as (Tobie et al., 2005): Tidal dissipation, and thus lunar recession, is more limited with Maxwell body.

Rayleigh-Taylor instability
When a magma ocean solidifies, it solidifies upwards from the bottom (Solomatov, 2015), and because viscosity increases by many orders of magnitude at the rheological transition at T = T 40 , the thermal structure of a solidifying magma ocean initially follows T 40 (z) (e.g., Figure 2a of Miyazaki and Korenaga (2019b)). However, dT 40 /dP is greater than the adiabatic gradient, meaning that the solid-melt mixture lying on T 40 (z) is gravitationally unstable, susceptible to the Rayleigh-Taylor instability (Solomatov, 2015;Miyazaki and Korenaga, 2019b). The time scale of this instability may be estimated as (Turcotte and Schubert, 1982) where η eff is the effective viscosity, ∆ρ is the density difference across the superadiabatic layer, g is gravitational acceleration, L is the thickness of the superadiabatic layer. After the Rayleigh-Taylor instability, the thermal structure becomes adiabatic, so the effective viscosity relevant to this transition may be calculated from the logarithmic average of two viscosity profiles, one along T 40 (z) and the other along the final adiabat. The density difference can be calculated from the difference between these two thermal profiles, using density and thermal expansivity parameterized in the subsection "Mantle adiabat". Some representative results are shown in Figure 2. For the reference case, the time scale becomes less than one year as soon as the surface temperature drops to ∼2800 K and remains below 100 years until the end of the magma ocean phase. Even with the harder background viscosity, the time scale is less than 10 3 years, which is still two orders of magnitude shorter than the time scale of magma ocean solidification, justifying the assumption of an adiabatic thermal structure throughout solidification.
Using adiabats calculated with a chosen mantle phase diagram is equivalent to assuming equilibrium crystallization, i.e., the composition of a solid-melt mixture does not evolve during magma ocean solidification. Fractional crystallization is equally likely (Solomatov, 2015), and it has been suggested that fractional crystallization is necessary for the rapid sequestration of atmospheric carbon in the early Earth (Miyazaki and Korenaga, 2022). The assumption of equilibrium crystallization in this study is merely for the sake of simplicity, as the mode of crystallization has only limited influence on the time scale of magma ocean solidification (Miyazaki and Korenaga, 2019b). We also do not consider the temporal evolution of melt fraction caused by melt percolation. Melt percolation in the solid phase reduces the melt fraction in one part, which increases viscosity thus reduces tidal dissipation, and creates a melt layer in another part, which does not contribute to tidal dissipation. Thus, the assumption of no melt percolation in a solidifying magma ocean maximizes its tidal dissipation.

Tidal dissipation
The tidal deformation of a compressible self-gravitating Earth caused by the Moon is calculated by solving the following coupled ordinary differential equations (Takeuchi and Saito, 1972;Sabadini et al., 2016): where r is the radius, ω is the tide-raising frequency, y 1 and y 2 are the radial and tangential displacements, respectively, y 3 and y 4 are the radial and tangential stresses, respectively, y 5 is the gravitational potential, and y 6 is the potential stress. Our formulation is identical to that of Sabadini et al.
(2016) except that we do not neglect inertia, which is important when considering a fluid layer.
The non-zero elements of the matrix A are: where l is the order of spherical harmonics, G is the universal gravitational constant, ρ o is the reference density, g is the corresponding gravitational acceleration, K is the bulk modulus, andλ andβ are defined asλ We use PREM to calculate ρ 0 , g, and K.
A region with φ > φ c is regarded as a fluid layer, in which the elastic shear modulus vanishes and the tangential stress y 4 is uniformly zero. Because A 24 diverges, y 2 is obtained by making use of the equation for dy 4 /dr (Takeuchi and Saito, 1972).
In case of the semidiurnal lunar tide (l = 2), these differential equations are subject to the following boundary conditions at the surface: y 3 = y 4 = 0 and where g s is the surface gravity, M M and M E are the masses of the Moon and Earth, respectively, R E is the Earth radius, and a is the lunar distance, and at the core-mantle boundary (CMB; r = r C ) (Sabadini et al., 2016): where g C is the gravity at CMB, ρ C is the core density (we assume a uniform density for the core), and C i are the constants to be determined to satisfy the boundary conditions at the surface. We use the shooting method with the 4th-order Runge-Kutta integration scheme to solve the above differential equations (Tobie et al., 2005).
To calculate tidal dissipation, the strain rate tensor is first obtained from y 1 and y 2 as where θ is the colatitude, φ is the longitude, and Y lm is the unnormalized spherical harmonics of degree l and order m. In case of a semidiurnal tide, l = m = 2. Note that the symbol φ is also used for melt fraction in this paper, but which one is used should be clear from the context. Then, the stress tensor is calculated asσ Finally, tidal dissipation is obtained by evaluating the following volume integral (Tobie et al., 2005): We validate our code using the case of "liquid Fe-FeS core + silicate mantle" in Tobie et al. (2005).
The dissipation of earth tide, as a function of the lunar distance, is shown in Figure 1c for the 'MK19' phase diagram and in Figure 3c for the 'const' phase diagram. As T 40 is lower for the latter, significant tidal dissipation takes place only after the potential temperature becomes lower than ∼2550 K.
Following Zahnle et al. (2015), tidal Q is calculated as Q tidal =Ė pot /Ė earth , with (Kaula, 1968) where k 2 is the Love number for Earth (= 0.3). During the magma ocean phase, there is no surface water ocean, so there is no need to consider the ocean tide. At the early phase of magma ocean solidification when the adiabat is still above T 40 ,Ė earth is calculated asĖ pot /Q max tidal where Q max tidal is set to 10 4 in the reference case.Tidal dissipation in an entirely fluid sphere is estimated to have Q tidal of 10 4 to 10 6 (Zahnle et al., 2015). Even after part of the magma ocean starts to experience the rheological transition, dissipation due to earth tide during the magma ocean phase is set to the the value given by this Q max tidal if the dissipation based on equation (53) is smaller.

Thermal evolution
The thermal evolution of Earth is tracked by solving the following energy balance: where t is time, T p is the potential temperature of the mantle, C m (T p ) is the heat capacity of the mantle as a function of its potential temperature, H m is radiogenic heating within the mantle, Q c is the core heat flux, and Q s is the surface heat flux. Time-dependent radiogenic heating is scaled to reproduce the present-day Urey ratio of 0.3 (Korenaga, 2006), and the core heat flux is linearly interpolated from 15 TW at the beginning and 10 TW at the present (O'Rourke et al., 2017). Being much smaller than other terms, however, the effects of radiogenic heating and core heat flux are negligible during the magma ocean phase. To take into account the release of the latent heat of crystallization upon cooling, the heat capacity function is constructed from the material properties used to calculate mantle adiabats (subsection "Mantle adiabat", Figure 4). In the petrology literature, the mantle potential temperature is defined as the temperature of the mantle adiabatically brought up to the surface without melting (McKenzie and Bickle, 1988), but in the magma ocean literature, the effect of melting is not suppressed (Solomatov, 2015) as doing so is obviously awkward. We follow the latter convention in this study. In all cases, we start with the initial potential temperature of 4500 K; model results are not very sensitive to the initial temperature as long as it is above the liquidus.
Surface heat flux from a magma ocean is calculated as (Solomatov, 2015): where A E is the surface area of Earth, k m is the melt thermal conductivity (2 W K −1 m −1 ), α s is the melt thermal expansivity (5 × 10 −5 K −1 ), ρ s is the surface melt density (3000 kg m −3 ), C p is the specific heat (10 3 J kg −1 K −1 ), g s is the surface gravity, and T s is the surface temperature. Earth's magma ocean is considered to be covered by a dense (∼100 bar) CO 2 -rich atmosphere (Zahnle et al., 2007), and based on the 1-D radiative-convective atmosphere modeling of Miyazaki and Korenaga (2022), the surface temperature is parameterized as where c 0 = 546 and c 1 = 0.63. Equation (56) is for soft turbulence, and heat flow scaling for hard turbulence, which may be more appropriate for magma oceans, predicts even higher heat flow (Solomatov, 2015). Here we choose the scaling for soft turbulence to provide a conservative estimate on the time scale of magma ocean solidification. The magma ocean stage ends when the mantle potential temperature reaches T 40 (z = 0) (∼1613 K), and the surface temperature is assumed to collapse to 500 K. Before subsolidus mantle convection eventually takes over, residual melt in the shallow upper mantle migrates upward and degas, forming a water ocean (Miyazaki and Korenaga, 2022).
Owing to contrasting solubilities in magma (Elkins-Tanton, 2008), most of carbon dioxide is in the atmosphere whereas most of water is in the magma ocean, and this situation is maintained during magma ocean solidification (Hier-Majumder and Hirschmann, 2017;Miyazaki and Korenaga, 2022). Given the likely water budget of the bulk Earth (∼1-3 ocean) (Hirschmann and Dasgupta, 2009;Korenaga et al., 2017), therefore, the amount of water in the atmosphere is too small to form a water ocean during the magma ocean phase, and the Komabayashi-Ingersoll limit (or the run-away greenhouse threshold) is not relevant because the limit requires the coexistence of two phases (water vapor and liquid water) (Pierrehumbert, 2011).

Orbital evolution
To track the lunar distance, we adopt the simplest assumption that the Moon evolves away from Earth maintaining a circular orbit (Zahnle et al., 2015). We also assume zero obliquity and zero inclination, so we are concerned only with the semidiurnal tide with the angular frequency of 2(Ω − n 1 ), where Ω is the angular velocity of Earth's rotation and n 1 is the angular velocity of the Moon on its orbit. This facilitates the comparison of our results with those of Zahnle et al. (2015), who made the same assumptions. We have tested with the effect of non-zero lunar inclination on tidal dissipation, with negligible differences.
With these assumptions, the orbital evolution is described by and where I E is Earth's moment of inertia and L is the total angular momentum of the Earth-Moon system: The lunar orbital angular velocity is updated as Time-stepping is controlled so that a relative change in the lunar distance change is less than 0.1 % and that in mantle potential temperature is less than 2 K at every time step.

Results
One of the unequivocal features of magma ocean solidification is that it solidifies from the lower mantle ( Figure 1). This is because the pressure dependence of the mantle liquidus is greater than that of the mantle adiabat. Consequently, tidal dissipation in a solidifying magma ocean takes place while the surface of the magma ocean is well above the liquidus, releasing a large amount of heat. In addition, the thermal structure of a solidifying magma ocean is always close to adiabatic, as opposed to isothermal, because of the Rayleigh-Taylor instability (Solomatov, 2015;Miyazaki and Korenaga, 2019b). This limits the depth extent of optimal viscosity that is important for tidal resonance (Figure 1b), and thus the magnitude of tidal dissipation ( Figure 1c). As such, tidal dissipation in a solidifying magma ocean is always substantially lower than its surface heat loss, being ineffective to slow down the solidification process. The completion of the magma ocean stage, defined by the time when the surface melt fraction reaches to ∼0.4, can take place within only ∼10 5 years (Figure 5a), and the lunar distance is still at ∼8 R E in the reference case ( Figure 5c).
As noted earlier, carbon dioxide and water, two dominant volatile species in Earth's magma ocean (Sossi et al., 2020), have contrasting solubilities in magma, the latter being more soluble than the former (Elkins-Tanton, 2008). When its melt fraction is reduced to ∼0.4, a solid-melt mixture experiences the rheological transition, jumping from melt viscosity to solid-like viscosity (Abe, 1993;Solomatov, 2015) (Figure 1b). The combination of the rheological transition and the Rayleigh-Taylor instability through a solidifying magma ocean results in the efficient entrapment of water in the mantle (Hier-Majumder and Hirschmann, 2017;Miyazaki and Korenaga, 2022), and after the completion of the magma ocean stage, a water ocean will form gradually first by the upward percolation of remaining melt in the shallow upper mantle, and then by degassing associated with subsolidus mantle convection (Miyazaki and Korenaga, 2022).
As soon as a water ocean emerges, the effect of ocean tides needs be considered, but given our limited understanding of the landscape of the early Earth, it is still impractical to model the details of ocean tides in the early Earth. To guide our discussion, therefore, hypothetical orbital evolution with a range of tidal Q/k 2 , where k 2 is the Love number, is shown in Figure 5c. To match the orbital reconstruction of Tyler (2021) or Daher et al. (2021) at the time of 10 7 years, for example, Q/k 2 has to be consistently as low as ∼2, which seems unrealistic. These studies do not consider continental growth nor variable ocean volume in the early Earth, so if we only aim to converge to their model predictions at the end of the Hadean (5 × 10 8 years) or the end of the Archean (2 × 10 9 years), Q/k 2 can take values that are usually assumed in previous studies (∼50-100) (Ward and Canup, 2000;Pahlevan and Morbidelli, 2015;Ćuk et al., 2021). It is important to note that the details of lunar recession after the solidification of the magma ocean (∼10 5 years) until the end of the Archean remain highly uncertain, because the period corresponds to the part of Earth history during which both continental mass and ocean volume can change dramatically (Korenaga, 2021b) and because tidal dissipation in the ocean can be very sensitive to such changes (Motoyama et al., 2020). Nevertheless, the Hadean tidal environment is unlikely to have been intense enough to challenge the possibility of wet-dry cycles on stable landmasses. The tidal amplitude is inversely proportional to the cube of the lunar distance; with the present-day value of ∼0.5 m, the amplitude can be up to only several meters after ∼4.3 Ga (i.e., after 200 Myr in Figure 6).
It is difficult to increase the extent of lunar recession during the magma ocean stage (Figures 5c   and 7). The maximum tidal Q during this stage is set to 10 4 in the reference case (Methods), which is already at the lower end of the acceptable range (Zahnle et al., 2015). Starting closer to Earth, i.e., with greater tidal dissipation, has only a self-limiting effect (Figure 5c). Varying mantle rheology within a reasonable range of its uncertainty can amount to a difference of only ∼±0.5 R E (Figure 7). The lunar distance at the end of the magma ocean stage can be pushed up to ∼8.5 R E with a softer mantle background, but this choice of background viscosity may not be appropriate given the partitioning behavior of water (Karato, 1986) (see the subsection "Rheology"). The lunar distance at the end of the magma ocean stage is most limited in the case of the harder mantle with the initial distance of 3.5 R E (dashed blue curve; ∼7.5 R E ). The effect of phase diagram choice, whether 'MK19' or 'const', is trivial. The most effective way is to increase melt viscosity from a peridotitic value (0.1 Pa s (Dingwell et al., 2004)) to a basaltic one (10 Pa s (McBirney, 1993)), which lowers surface heat flux by a factor of ∼5 and thus lengthens the duration of the magma ocean stage. Fractional crystallization, if it takes place, does modify the composition of surface magma, but the composition remains largely peridotitic during the bulk of the solidification (Miyazaki and Korenaga, 2019b); thus, the case with basaltic melt viscosity (Figures 5c and 7) should be regarded as an unlikely end-member.
The effect of varying the critical melt fraction φ c is also minor on lunar recession ( Figure 8).
Also, the use of Maxwell viscoelasticity leads to more reduced lunar recession as expected (Figure 8). Our calculation of orbital evolution does not take into account the effect of tidal dissipation within the Moon. However, doing so would only result in reducing the lunar semimajor axis (Murray and Dermott, 2000), thereby strengthening the case of limited lunar recession during magma ocean solidification.

Discussion
Differences from the model of Zahnle et al. (2015), which is also based on the time-forward approach with magma ocean solidification, are striking (Figure 5c), and most of them can be traced back to their isothermal approximation. This approximation was necessary for their analysis because their calculation of tidal dissipation assumed homogeneous material properties for the entire Earth. In our model, a magma ocean starts to experience the rheological transition and a significant amount of tidal dissipation as soon as the surface temperature drops below ∼3000 K (Figure 1), but with the isothermal approximation, in which the solidus and liquidus temperatures do not change with depth (Figure 9), such a transition has to wait until the surface melt fraction reaches ∼0.4, which happens at the surface temperature of ∼1560 K in their model. Such low magma temperature limits surface heat flow, and at the same time, tidal dissipation is maximized as the entire mantle has low viscosities around the rheological transition. Tidal heating is thus effectively trapped in Earth's magma ocean, extending its lifetime by tens of million years, during which substantial lunar recession can take place (Figure 5c). The isothermal approximation introduces about three orders of magnitude difference in the solidification time scale, which results in a factor of ∼5 difference in the Earth-Moon distance. This may not be surprising given the drastic nature of the isothermal approximation ( Figure 9).
As discussed above, the isothermal approximation made by Zahnle et al. (2015) maximizes the total amount of tidal dissipation (1) (Miyazaki and Korenaga, 2019b), solidification from a mid-mantle level could still occur using a thermodynamic database different from that of Miyazaki and Korenaga (2019a), and if it happens, a basal magma ocean would form (Labrosse et al., 2007;Mosenfelder et al., 2009;Stixrude et al., 2009). Even with bottom-up solidification, a basal magma ocean could still form if matrix compaction is efficient (Miyazaki and Korenaga, 2019b). A basal magma ocean would solidify from the top, with a limited extent of a partially molten boundary layer (Labrosse et al., 2007), so tidal dissipation associated with a basal magma ocean is unlikely to be significant. In general, matrix compaction (or equivalently, melt percolation) increases the vertical extent of a fluid layer (φ > φ c ) and reduces the melt fraction of other solid parts of the mantle (φ ≤ φ c ). Tidal dissipation within a fluid layer is considerably lower than that within a solid layer, even using the maximum fluid tidal Q of 10 4 (Figures 5b), and reduction in the melt fraction drives the solid layer away from the rheological transition, which acts to reduce tidal dissipation. In other words, the solidification scenario of our model probably corresponds to the case of the maximum tidal dissipation, and the most important finding of this study is that, even with this optimal scenario, lunar recession during magma ocean solidification is severely limited.
Limited lunar recession during the early Hadean, seen in our model results, may resurrect the early origin of the lunar inclination, either by the evection resonance (Touma and Wisdom, 1998) or by interaction with remanent proto-lunar disk materials (Ward and Canup, 2000). The presentday lunar orbit has an inclination of ∼5 • with respect to the ecliptic plane, and the origin of this inclination has been the most serious dynamical problem for the formation mechanism of the Moon (Goldreich, 1966;Touma and Wisdom, 1998;Ward and Canup, 2000;Pahlevan and Morbidelli, 2015;Ćuk et al., 2016;Tian and Wisdom, 2020). To preserve the primordial lunar inclination, it is suggested that the lunar magma ocean has to solidify before the onset of the Cassini state transition at ∼30 R E (Chen and Nimmo, 2016). As the lunar magma ocean is generally considered to have solidified with a time scale of 10-100 Myr (Elkins-Tanton et al., 2011;Chen and Nimmo, 2016), such preservation appears impossible with the lunar recession model of Zahnle et al. (2015) ( Figure 5c). A proper treatment of magma ocean solidification, however, suggests that lunar recession is limited to only ∼8 R E at the end of the magma ocean phase, and the preservation of primordial lunar inclination is certainly within the bounds of possibility (Figure 5c).
The late origin of the lunar inclination, by collisionless encounters with left-over planetesimals (Pahlevan and Morbidelli, 2015), remains viable although the efficiency of collisionless excitation is likely to be reduced. This mechanism, which operates on the time scale of 10 8 years, is more efficient for a larger cross-section of the Earth-Moon system, and with the recession model of Zahnle et al. (2015), the lunar distance is already at ∼40 R E at the time of 10 7 years. In light of our study, such recession requires an unrealistically low Q/k 2 for an ocean-covered Earth. With more likely values of Q/k 2 (∼ 100), the efficiency of collisionless excitation would be reduced by ∼50% than originally suggested (Figure 5c). Perhaps most important is that, with limited early lunar recession, the early and late origins of the lunar inclination are not mutually exclusive; all of them could contribute, facilitating the resolution of the lunar inclination problem.
It is more difficult to assess the impact of our results on the lunar orbital evolution with high initial obliquity (Ćuk et al., 2016, 2021). The duration of such a scenario extends over ∼100 Myr, so most of the suggested evolution must take place in the presence of a water ocean (Figure 5c).
It is common to use constant Q/k 2 for tidal dissipation in previous studies on lunar orbital evolution (Ward and Canup, 2000;Pahlevan and Morbidelli, 2015;Ćuk et al., 2016, but ocean tides are expected to have been highly variable on the early Earth (Motoyama et al., 2020), owing to the competing effects of the growing ocean and continents (Korenaga, 2021b). For the high initial obliquity scenario to work, the Earth-Moon system needs to remove a considerable fraction of the initial angular momentum through a series of resonances (Ćuk et al., 2021), but the removal of angular momentum by resonance is efficient only for a narrow range of tidal parameters (Rufu and Canup, 2020). A more incisive evaluation of the lunar orbital evolution and the formation mechanism of the Moon, therefore, hinges on an improved understanding of how early surface environments were shaped by the dynamics of Earth's interior.   of the day (red) for Earth's first one billion years, with a range of tidal Q/k 2 : 1 (solid), 10 (dashed), 10 2 (dot-dashed), and 10 3 (dotted).   Figure 7: Early history of lunar recession with a range of parameter combinations: four choices of mantle rheology (red -reference, purple -softer background with α η = 31 (i.e., softest solid viscosity), orangesofter background with α η = 31 and melt viscosity of 10 Pa s, and blue -harder background), the initial lunar distance of 5 R E (solid) and 3.5 R E (dashed), and the 'MK19' (solid) and 'const' (dotted) phase diagrams.   Figure 9: The phase diagram of the mantle and 'adiabats' assumed in Zahnle et al. (2015).