Issue 
A&A
Volume 665, September 2022



Article Number  L1  
Number of page(s)  20  
Section  Letters to the Editor  
DOI  https://doi.org/10.1051/00046361/202243445  
Published online  05 September 2022 
Letter to the Editor
The resonant tidal evolution of the EarthMoon distance
IMCCE, CNRS, Observatoire de Paris, PSL University, Sorbonne Université, 77 Avenue DenfertRochereau, 75014 Paris, France
email: mohammad.farhat@obspm.fr
Received:
1
March
2022
Accepted:
13
August
2022
Due to tidal interactions in the EarthMoon system, the spin of the Earth slows down and the Moon drifts away. This recession of the Moon can now be measured with great precision, but it was noticed more than fifty years ago that simple tidal models extrapolated back in time lead to an age of the Moon that is largely incompatible with the geochronological and geochemical evidence. In order to evade this problem, more elaborate models have been proposed, taking into account the oceanic tidal dissipation. However, these models have not been able to fit both the estimated lunar age and the present rate of lunar recession simultaneously. In the present work, we present a physical model that reconciles these two constraints and yields a unique solution for the tidal history. This solution fits the available geological proxies for the history of the EarthMoon system well and it consolidates the cyclostratigraphic method. Our work extends the lineage of earlier works on the analytical treatment of fluid tides on varying bounded surfaces that is further coupled with solid tidal deformations. This allows us to take into account the timevarying continental configuration on Earth by considering hemispherical and global ocean models. The resulting evolution of the EarthMoon system involves multiple crossings of resonances in the oceanic dissipation that are associated with significant and rapid variations in the lunar orbital distance, the length of an Earth day and the Earth’s obliquity.
Key words: Earth / Moon / planets and satellites: dynamical evolution and stability / planets and satellites: oceans
© M. Farhat et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Due to the tidal interplay in the EarthMoon system, the spin of the Earth slows down over time and the EarthMoon distance increases (Darwin 1879) at a current rate of 3.830 ± 0.008 cm yr^{−1}, measured using Lunar Laser Ranging (LLR; Williams & Boggs 2016). A rich narrative has been cultivated, exploring the longterm evolution of the system (Goldreich 1966; Mignard 1979; Touma & Wisdom 1994; Neron de Surgy & Laskar 1997) and the dynamical constraints on the origin of the Moon (Touma & Wisdom 1998; Ćuk et al. 2019). Overall, it has been established that simple tidal models starting with the present recession rate and integrated backward in time predict a close encounter in the EarthMoon system within less than 1.6 billion years (Ga; Gerstenkorn 1967; MacDonald 1967). This assumption is clearly incompatible with the estimated age of the Moon of 4.425 ± 0.025 Ga (Maurice et al. 2020), which suggests that the present rate of rotational energy dissipation is much greater than it has typically been over the Earth’s history. To bypass this difficulty, empirical models have been fitted to the available geological evidences of the past rotational state of the Earth, (Walker & Zahnle 1986; Waltham 2015), acquired through the analysis of paleontological data (e.g., Williams 2000), sedimentary records of tidal rhythmites, (Williams 1997, 2000; Sonett & Chan 1998; Eriksson & Simpson 2000; de Azarevich & Azarevich 2017), or Milankovitch cyclostratigraphic sequences (Meyers & Malinverno 2018; Huang et al. 2020; Sørensen et al. 2020; Lantink et al. 2021). However, empirically fitting these geological data brings little insight into how different physical components have actually contributed to the evolution of the EarthMoon system.
Major progress has been achieved with the elaboration of oceanic tidal models. These models offer a tidal frequencydependent dissipation behavior (LonguetHiggins 1968; Platzman 1984; Müller 2008a) that allows for the encounter of highdissipation resonant states over the course of Earth history (Webb 1980; AuclairDesrotour et al. 2018; Tyler 2021). However, the literature lacks an effective model stemming from controlled analytical formulations that fits both the currently measured rate of lunar recession and the estimated age of the Moon.
Alongside its dependence on the Earth’s rotation rate, the varying continental configuration has also played a role in enhancing the oceanic resonances or even exciting additional ones (Platzman 1983; Ooe 1989; Tyler 2021). Paleodissipation might have also varied significantly during ice ages, as areas of continental shelves vary with sea level (Griffiths & Peltier 2009; Arbic & Garrett 2010). However, both ice ages and basin geometry cycles have much smaller periodicities compared to the Earth’s age (Boulila et al. 2018; Farhat et al. 2022). Moreover, accurately accounting for such level of realism is hindered by the accumulating uncertainty in deeptime modeling. It is thus necessary to compromise between the practicality of effective models with simplified geometries (Webb 1980; Hansen 1982; Tyler 2021) and the realism of costly numerical models that depend on paleogeographic reconstructions (Green et al. 2017; Daher et al. 2021).
Here, we undertake a systematic exploration of the timevarying tidal dissipation in the oceans. We propose a physical model that reconciles the two aforementioned limits, described in Sect. 2. With a minimum number of free parameters, we constrain our model to only fit the two most certain points of lunar evolution history: the present rate of lunar recession and the lunar age, presented in Sect. 3. This provides a unique solution to the EarthMoon separation history, described in Sect. 4. We focus on the computation of the tidal response of the Earth, considering a reconstruction of the continental drift up to one billion years ago, followed by a smooth transition toward a global ocean planet. We supplement this computation with a reduced dynamical model of the system that captures the skeletal structure of the longterm evolution based on robust features of the tidal response. However, we anticipate that this model may serve as the backbone of a fully spatial dynamical evolution in the system (more on that in Appendix A). The orbital solution that we produce demonstrates the robustness of the cyclostratigraphic machinery and further suggests interesting intervals for future investigations, as shown in Sect. 5.
2. Oceanic model
In our model, we compute the tidal response of the oceans and the solidEarth to lunisolar semidiurnal forcing, both combined with mimetic continental drift driven by plate tectonics. We focus on the dependence of dissipation on the Earth’s spin rate. We combine two analytical approaches that describe longwavelength barotropic tidal flows over shallow spherical and hemispherical shells. The spherical shell describes a global ocean that we assume had existed in the earliest eons of the lifetime of the Earth (Motoyama et al. 2020). The existence of an early ocean is supported by evidence from the analysis of detrital zircon around 4.4 Ga (Wilde et al. 2001), from the interaction between the ocean and continental crust 4 billion years ago (Mojzsis et al. 1996), and from records of the oxygen isotope composition of seawater (Peck et al. 2001; Johnson & Wing 2020). The “globality” of this ocean is justified by the analysis of continental crust growth curves based on geochemical evidence in zircon crystallization ages (Dhuime et al. 2012; Hawkesworth et al. 2020). In compliance with these curves, we consider that a hemispherical oceanic shell has taken over in the most recent times. In our model, the center of this hemispheric continental cap follows the evolution of the paleogeographic center. In doing so, we emphasize on the role of “continentality” in the tidal response, while avoiding the undersampling of geometric scenarios due to theoretical limitations (Hansen 1982; Tyler 2021) or due to uncertainties in plate tectonic models (Matthews et al. 2016; Daher et al. 2021). To compute this evolution, we adopt the recently developed paleogeographic reconstructions that cover the past billion years (Merdith et al. 2021). A postprocessing of these reconstructions allows us to produce the latitudinal evolution of the center of the continental cap captured in Fig. 1. The tidal frequencies at which oceanic resonances are excited and the amplitudes of these resonances vary with the surface position of the hemispherical ocean (Fig. B.1). Supercontinental formations and breakups thus have their mark on the predicted lunar recession rate.
Fig. 1. Temporal evolution of the latitude of the surface “paleobarycenter” over the last one billion years. The plate tectonics reconstruction is adopted from Merdith et al. (2021), which establishes the first kinematically continuous tectonic motion model across multiple supercontinental cycles. The evolution is smoothed in red using a moving polynomial regression filter with a window of 200 Myr. In our effective model, this curve maps the evolution of the center of the hemispheric continental cap that transitions from being symmetric about the equator during the Mesozoic, to being almost polar during the Paleozoic. 
The dynamical evolution of the EarthMoon system is coupled to the computation of the tidal flows, which is also dependent on the chosen oceanic geometry. For the global ocean, the tidal torque is computed by solving the modified Laplace tidal equation using the Hough functions as eigenfunctions (see Appendix K and AuclairDesrotour et al. 2019). The oceanic dissipation is parameterized by an effective frequency, σ_{R}, which globally models the bottom friction and the conversion of barotropic flows into internal gravity waves, both mechanisms amounting to ∼91% of the total dissipation (Carter et al. 2008). This frequency, σ_{R}, can also be interpreted as the inverse of a dissipation timescale, τ, that quantifies the time needed to deplete the kinetic energy budget of tidal oscillations after switching off the forcing. Although σ_{R} is probably a function of local topography, its spatial variation can be averaged out longitudinally over the Earth’s rotation and latitudinally over precession and plate tectonics. The second free parameter in our model is the uniform effective oceanic thickness, H. The imprints of these two parameters on the tidal response spectrum are distinguishable: variations in H smoothly shifts the positions of the resonant peaks while slightly varying their amplitudes. In contrast, variations in σ_{R} can completely reshape the tidal spectrum, amplifying the resonant peaks by several orders of magnitude when σ_{R} decreases or, otherwise, completely absorbing the resonant peaks into the background spectrum (Fig. K.1). For the hemispherical geometry, we adopt the analytical approach of Webb (1980) (see Appendix E), in which the tidal solution is expanded in spherical harmonics (Fig. E.1). In both geometries, we take into account the effect of the deformation of the solid part of the Earth adopting an Andrade rheology (Anderson & Minster 1979; CastilloRogez et al. 2011; Lau & Faul 2019; Appendix F).
3. Constraining the effective parameters
Assuming a reduced planar orbital model (described in Appendix A), we compute the evolution of the EarthMoon system that results from the lunisolar semidiurnal tidal torque for ranges of values of our effective parameters (H, σ_{R}). We do so for three models that ascend in realism: a global ocean model across the full geological history (similar to Tyler 2021); an “average” hemispherical ocean model across the full geological history (similar to Webb 1982) for which the response at any tidal frequency is averaged over all possible oceanic positions on the sphere; and our combined model that starts at the present with the hemispherical ocean evolving with the mimetic continental drift, then switching to the global ocean. For every constructed history of the EarthMoon separation, we compute the chisquared χ^{2}, taking only two data points into account: the wellconstrained lunar age of 4.425 ± 0.025 Ga (Maurice et al. 2020) and the currently measured rate of lunar recession of 3.830 ± 0.008 cm yr^{−1} (Williams & Boggs 2016). Misfit surfaces of χ^{2} for the three models are shown in Fig. 2. Two χ^{2} local minima exist for the global oceanic response; however, one of them corresponds to an unreasonably large average oceanic depth H ≈ 5500 m, leaving us with a global minimum of (H, log_{10}σ_{R})=(2273 m, −4.89), where σ_{R} is in s^{−1}. The global minimum in the “average” hemispherical ocean model corresponds to (H, log_{10}σ_{R})=(3816 m, −4.54), which is close to the average depth of the Pacific Ocean (Amante & Eakins 2009). For the combined model, the global minimum corresponds to (H, log_{10}σ_{R})=(4674 m, −5.19), where H is the thickness for the hemispherical phase of the model – which is twice that of the global ocean phase during earlier eons (Appendix B). The switch between the two geometries occurs at t_{switch}, which is implicitly determined by the dynamical integrator (Appendix B). For the bestfit solution, we have t_{switch} = 3.25 Ga, which is in agreement with assumptions from the literature of the existence of a global ocean until ∼2.5 Ga (Dhuime et al. 2012; Hawkesworth et al. 2020; Motoyama et al. 2020). If we assume that the oceanic volume is conserved over time, the bestfit value of H for the combined model corresponds to a volume of 1.19 × 10^{18} m^{3}, which is only 10% off from the currently estimated value of 1.33 × 10^{18} m^{3} based on global relief models (Amante & Eakins 2009). The fitted dissipation frequency, σ_{R}, corresponds to a decay time τ = 43.1 h, which is consistent with real oceanic studies (Garrett & Munk 1971; Webb 1973) that offer a range between 24 and 60 h (or log_{10}σ_{R} ∈ [ − 4.93, −5.33]). The bestfit values for the combined model correspond to a lunar trajectory characterized by a current rate of recession ̇a_{0} = 3.829 cm yr^{−1} and an impact time at 4.431 Ga (Table C.1).
Fig. 2. Misfit surfaces of χ^{2} for the three studied geometric models. The past dynamical evolution of the EarthMoon system is reconstructed for the shown ranges of our two free model parameters H and σ_{R}. The misfit is established using the currently measured lunar recession rate via LLR, and the lunar age (Appendix C). The three models differ in the imposed geometry of the oceanic shell over the geological history, with the combined model featuring more physical realism that the other two. The numerical results of this analysis are summarized in (Table C.1). The dynamical evolution associated with each of the misfit minima is plotted in terms of: the lunar semimajor axis in Fig. 3, length of the day in Fig. 5, and obliquity and precession frequency in Fig. 6. 
4. EarthMoon separation: A history of surfing resonances
For each of the global minima of the misfit parametric studies, we plot the evolution of the EarthMoon distance in Fig. 3. At the top of the evolution, we spread a compilation of geological proxies from tidal rhythmites (Walker & Zahnle 1986; Sonett & Chan 1998; Williams 2000; Eriksson & Simpson 2000; de Azarevich & Azarevich 2017) and cyclostratigraphy (Meyers & Malinverno 2018; Huang et al. 2020; Sørensen et al. 2020; Lantink et al. 2021; Tables D.1 and D.2). The three models are constrained at the end points, thus differences arise mostly in between. To better elaborate on the models’ discrepancies, we plot in Fig. 4 the temporal evolution of the tidal torque (normalized by its present value) associated with the combined model. The corresponding evolution of the length of an Earth day (LOD), precession frequency, and obliquity are plotted in Figs. 5 and 6. As it is directly proportional to tidal dissipation, the longterm evolution of the torque is characterized by a nonmonotonic increase, characteristic of the shrinking EarthMoon separation, that is interrupted by multiple crossings of resonances. The distribution of resonances in the hemispherical configuration, t < t_{switch}, is less regular than that in the global configuration, t > t_{switch}, (see also Figs. K.1 and K.2 for a global description of the tidal response spectrum). Each resonance crossing in the torque generates an inflection point in the evolution of a_{M}, which depends on the width and to a lesser degree on the amplitude of the resonance peak (AuclairDesrotour et al. 2014). Figure 4 depicts a critical feature of the combined model: starting with the hemispherical geometry at present, the torque is located around a resonance peak, which provides a higher dissipation rate than for the global ocean configuration. This models the anomalous present rate of dissipation attributed to the blocking of westward tidal propagation by the current continental distribution and the effect of enhanced dissipation by continental shelves (Arbic et al. 2009). The first phase of the model involves two major resonances between the present and 700 Ma, resulting in cascade falls of a_{M} of 2.8R_{E} within 330 Myr. These resonances are associated with rapid variations of the Earth’s obliquity (Fig. 6) that could have triggered major climatic events. We observe that the first resonance overlaps with the Paleozoic oxygenation event (∼350 Ma), while the second overlaps with the Neoproterozoic major oxygenation event (∼600 Ma) and the Cambrian Explosion (Wood et al. 2019). Possible correlation between the Earth’s LOD and the benthic ecosystem should thus be considered (Klatt et al. 2021). The second resonance peak is almost half an order of magnitude lower than in the global configuration. This is an essential feature of the combined model for preserving the lunar angular momentum budget at this stage to better match the cyclostratigraphic proxy estimates at 1.4 and 2.5 Ga, which clearly cannot be explained by the other, more dissipative models considered in Fig. 2.
Fig. 3. Evolution of the lunar semimajor axis over time. The EarthMoon separation, a_{M}, is plotted for the three studied models, taking the bestfit values of the free parameters (H, σ_{R}) as described in Fig. 2 and in the main text. Plotted on top of the evolution curves: Geological inferences of a_{M} from cyclostratigraphy and tidal laminae data (Tables D.1 and D.2). The shaded envelope corresponds to 2σuncertainty in the fitted parameters of the combined model (Appendix C). In the narrow window, we zoom over the most recent 250 Myr of the evolution and make a comparison with the evolution corresponding to explicit numerical tidal modeling using paleogeographic reconstructions (Green et al. 2017) and the prediction of the numerical solution La2004 (Laskar et al. 2004). We note that the integration of a_{E} extends to 3R_{E}, but the yaxis is trimmed to start at 15R_{E} for a better visualization of the geological data. 
Fig. 4. History of the tidal torque. The logarithm of the semidiurnal tidal torque of the Earth (normalized by its present value: ) is plotted as a function of time. The solid curve corresponds to the torque of the combined model that involves three phases: in the first phase, a hemispherical ocean migrates on the surface of the Earth following the evolution of the continental barycenter of Fig. 1. Given we lack a continuous plate tectonics model beyond 1 Ga, in Phase 2, we fix the hemispherical ocean to its configuration at 1 Ga to avoid discontinuities in the modeling. It is noteworthy that the attenuated tidal torque over this phase is not due to the fixed oceanic position but due to the tidal response occupying the nonresonant background of the spectrum for the tidal frequencies associated with this interval. Beyond t_{switch}, we enter Phase 3 of the model with the global ocean configuration. The dashed and dasheddotted curves correspond, respectively, to the global and hemispherical oceanic torques that are ignored over the specified intervals by the selective combined model. 
Fig. 5. Evolution of the Earth’s length of the day with time. Similar to Fig. 3, but here for the LOD evolution associated with the three studied oceanic models. Geological data on the LOD are summarized in Tables D.1 and D.2. The minimal value reached for the LOD when the integration is terminated at 3 Earth radii is 5.25 h. 
Fig. 6. Evolution of the Earth’s obliquity, precession frequency, and precession period with time. The evolution of a_{M} (Fig. 3) and LOD (Fig. 5) are used to compute the evolution of obliquity and precession by Eqs. (A.6) and (A.7). The geological data of the precession frequency from tidal rhythmites and cyclostratigraphy are also plotted on top of the curve (Tables D.1 and D.2). We note that the precession frequency is the directly measured observable in cyclostratigraphy. 
Following these resonances, the torque enters a long nonresonant interval associated with the intrinsic tidal response occupying the background of the spectrum (Fig. K.1). This “dormant” torque phase covers the interval of the socalled “boring billion years” associated with stabilized rates of atmospheric oxygenation (Alcott et al. 2019). Entering the oceanic global geometry phase of the combined model occurs at 3.25 Ga, namely after covering all significant supercontinental cycles, although t_{switch} is implicitly determined by the dynamical integrator (Appendix B). Samples of continental growth curves predict a fast decay in continental crust volume beyond t_{switch} (Sun et al. 2019; Hawkesworth et al. 2020). After switching to the global ocean response spectrum, the torque passes through a major resonance around 3.35 Ga, resulting in a significant and abrupt drop in a_{M} of 6.5R_{E} within 250 Myr. Beyond this age, the evolution again follows the tidal dissipation background spectrum before terminating with the impact.
5. New target for geological studies
In this work, we present the first semianalytical physical model that fits the most accurate constraints in the EarthMoon evolution: the present tidal dissipation rate and the age of the Moon. We deliberately avoided fitting our model to any of the available geological data. The unique solution of our combined model is a nearly perfect match to a large set of those geological data (Figs. 3, 5 and 6). This solution will provide a new target for geological studies, as it clearly validates the cyclostratigraphic approach, which estimates the Earth’s precession frequency from stratigraphic sequences (Meyers & Malinverno 2018; Huang et al. 2020; Sørensen et al. 2020; Lantink et al. 2021; Table D.1). In particular, the cyclostratigraphic evaluation of the EarthMoon distance at 2459 ± 1.3 Ma in the Joffre banded iron formations (BIF, Lantink et al. 2021) is in remarkable agreement with our model, compared to the equivalent estimates deciphering tidal rhythmites in the (∼2450 Ma) Weeli Wooli BIF in Australia (Walker & Zahnle 1986; Williams 2000). Our target curve can probably now be used to elaborate robust procedures for the analysis of these tidal rhythmites that led sometimes to divergent interpretations (Walker & Zahnle 1986; Sonett & Chan 1998; Williams 2000; Table D.2).
We obtained a striking fit with the estimate of a_{M} at 3.2 Ga obtained through the analysis of the Moodies group rhythmites (Eriksson & Simpson 2000; de Azarevich & Azarevich 2017), but we do not deny that this agreement could be coincidental and a new analysis of these sections in association with cyclostratigraphic estimates is certainly welcome. We expect that substantial progress will be made in the near future with the analysis of many cyclostratigraphic records, which could then be used to further constrain our physical model. The sequences that occur during the resonant states (or in their vicinity), corresponding to the steep slopes in Fig. 3, are of particular interest. Finally, as this model provides a coherent history of the EarthMoon distance, it can also be used to constrain the timescale of the lunar formation scenario (Ćuk et al. 2016). This coherence between the geological data and the present scenario for the EarthMoon evolution will also promote the use of these geological data and, in particular, the cyclostratigraphic geological data as a standard observational window for recovering the past history of the solar system.
Acknowledgments
We thank Maëlis Arnould for her help with the plate tectonics model and the GPlates software and Matthias Sinnesael for discussions on the geological data. We are grateful to Margriet Lantink and coworkers for the communication of their results on the Joffre sequence before publication and allowing us to include their data point in the present work. This project has been supported by the French Agence Nationale de la Recherche (AstroMeso ANR19CE31000201) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo885250). This work was granted access to the HPC resources of MesoPSL financed by the Region ÎledeFrance and the project Equip@Meso (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. Contributions: MF conducted the simulations and drafted the paper. PAD brought expertise in oceanic tides and GB in solid tides. JL initiated the study and supervised it. All contributed to the study at all stages. All contributed to the writing of the paper.
References
 Abramowitz, M., Stegun, I. A., & Romer, R. H. 1988, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables [Google Scholar]
 Alcott, L. J., Mills, B. J., & Poulton, S. W. 2019, Science, 366, 1333 [NASA ADS] [CrossRef] [Google Scholar]
 Amante, C., & Eakins, B. W. 2009, NOAA Technical Memorandum NESDIS NGDC24 [Google Scholar]
 Anderson, D. L., & Minster, J. B. 1979, Geophys. J. Int., 58, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Andrade, E. N. D. C. 1910, Proc. R. Soc. London. Ser. A, 84, 1 [CrossRef] [Google Scholar]
 Arbic, B. K., & Garrett, C. 2010, Cont. Shelf Res., 30, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Arbic, B. K., Karsten, R. H., & Garrett, C. 2009, Atmos. Ocean, 47, 239 [CrossRef] [Google Scholar]
 Arfken, G. B., & Weber, H. J. 1999, Mathematical Methods for Physicists [Google Scholar]
 AuclairDesrotour, P., Le PoncinLafitte, C., & Mathis, S. 2014, A&A, 561, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Mathis, S., Laskar, J., & Leconte, J. 2018, A&A, 615, A23 [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Leconte, J., Bolmont, E., & Mathis, S. 2019, A&A, 629, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bagheri, A., Khan, A., AlAttar, D., et al. 2019, J. Geophys. Res. Planets, 124, 2703 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, T., Jr 1975, J. Geophys. Res., 80, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Bolmont, E., Breton, S. N., Tobie, G., et al. 2020, A&A, 644, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boué, G., & Laskar, J. 2006, Icarus, 185, 312 [Google Scholar]
 Boulila, S., Laskar, J., Haq, B. U., Galbrun, B., & Hara, N. 2018, Global Planet. Change, 165, 128 [CrossRef] [Google Scholar]
 Boyden, J. A., Müller, R. D., & Gurnis, M. 2011, Nextgeneration Platetectonic Reconstructions using GPlates (Cambridge University Press) [Google Scholar]
 Carter, G. S., Merrifield, M., Becker, J. M., et al. 2008, J. Phys. Oceanogr., 38, 2205 [NASA ADS] [CrossRef] [Google Scholar]
 Castelnau, O., Duval, P., Montagnat, M., & Brenner, R. 2008, J. Geophys. Res. (Solid Earth), 113, B11203 [NASA ADS] [CrossRef] [Google Scholar]
 CastilloRogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. (Planets), 116, E09008 [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010, in Exoplanets (Tucson, AZ: University of Arizona Press), 239 [Google Scholar]
 Correia, A. C., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crease, J. 1966, Tables of the Integral (National Institute of Oceanography) [Google Scholar]
 Ćuk, M., Hamilton, D. P., Lock, S. J., & Stewart, S. T. 2016, Nature, 539, 402 [CrossRef] [Google Scholar]
 Ćuk, M., Hamilton, D. P., & Stewart, S. T. 2019, J. Geophys. Res. Planets, 124, 2917 [Google Scholar]
 Daher, H., Arbic, B. K., Williams, J. G., et al. 2021, J. Geophys. Res. Planets, 126, e2021JE006875 [CrossRef] [Google Scholar]
 Darwin, G. H. 1879, Philos. Trans. R. Soc. London Ser. I, 170, 447 [NASA ADS] [Google Scholar]
 de Azarevich, V. L. L., & Azarevich, M. B. 2017, GeoMarine Lett., 37, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Dhuime, B., Hawkesworth, C. J., Cawood, P. A., & Storey, C. D. 2012, Science, 335, 1334 [NASA ADS] [CrossRef] [Google Scholar]
 Dong, S.H., & Lemus, R. 2002, Appl. Math. Lett., 15, 541 [CrossRef] [Google Scholar]
 Dziewonski, A. M., & Anderson, D. L. 1981, Phys. Earth Planet. Inter., 25, 297 [CrossRef] [Google Scholar]
 Efroimsky, M. 2012, ApJ, 746, 150 [Google Scholar]
 Efroimsky, M., & Williams, J. G. 2009, Celestial Mech. Dyn. Astron., 104, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksson, K. A., & Simpson, E. L. 2000, Geology, 28, 831 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, J., Wu, H., Fang, Q., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 540, 109530 [NASA ADS] [CrossRef] [Google Scholar]
 Farhat, M. A., & Touma, J. R. 2021, MNRAS, 507, 6078 [NASA ADS] [CrossRef] [Google Scholar]
 Farhat, M., Laskar, J., & Boué, G. 2022, J. Geophys. Res. Solid Earth, 127, e2021JB023 [CrossRef] [Google Scholar]
 Farrell, W. 1972, Rev. Geophys., 10, 761 [NASA ADS] [CrossRef] [Google Scholar]
 Fienga, A., Deram, P., Di Ruscio, A., et al. 2021, Notes Scientifiques et Techniques de l’Institut de Mécanique Céleste, 110 [Google Scholar]
 Findley, W. N., Lai, J. S., Onaran, K., & Christensen, R. M. 1977, J. Appl. Mech., 44, 364 [NASA ADS] [CrossRef] [Google Scholar]
 FoxKemper, B., Ferrari, R., & Pedlosky, J. 2003, J. Phys. Oceanogr., 33, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Garrett, C., & Munk, W. 1971, Deep Sea Res. Oceanogr. Abstracts, 18, 493 [CrossRef] [Google Scholar]
 Gent, P. R., & McWilliams, J. C. 1983, Dyn. Atmos. Oceans, 7, 67 [CrossRef] [Google Scholar]
 Gerkema, T., & Zimmerman, J. 2008, in Lecture Notes (Texel: Royal NIOZ), 207 [Google Scholar]
 Gerstenkorn, H. 1967, Icarus, 7, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P. 1966, Rev. Geophys., 4, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Green, J., Huber, M., Waltham, D., Buzan, J., & Wells, M. 2017, Earth Planet. Sci. Lett., 461, 46 [CrossRef] [Google Scholar]
 Griffiths, S. D., & Peltier, W. R. 2009, J. Clim., 22, 2905 [NASA ADS] [CrossRef] [Google Scholar]
 Gurnis, M., Turner, M., Zahirovic, S., et al. 2012, Comput. Geosci., 38, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Han, L., & Huang, R. X. 2020, J. Phys. Oceanogr., 50, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, K. S. 1982, Rev. Geophys., 20, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Hawkesworth, C., Cawood, P. A., & Dhuime, B. 2020, Front. Earth Sci., 8 [CrossRef] [Google Scholar]
 Hough, S. S. 1898, Philos.Trans. R. Soc. London Ser. A, 191, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, H., Gao, Y., Jones, M. M., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 550, 109735 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, B. W., & Wing, B. A. 2020, Nat. Geosci., 13, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Kaula, W. M. 2013, Theory of Satellite Geodesy: Applications of Satellites to Geodesy (Courier Corporation) [Google Scholar]
 Klatt, J. M., Chennu, A., Arbic, B. K., Biddanda, B., & Dick, G. J. 2021, Nat. Geosci., 14, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Lantink, M., Davies, J., & Hilgen, F. 2021, Nat. Rev., 12, 369 [Google Scholar]
 Laskar, J. 2005, Celestial Mech. Dyn. Astron., 91, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lau, H. C., & Faul, U. H. 2019, Earth Planet. Sci. Lett., 508, 18 [CrossRef] [Google Scholar]
 Lau, H. C., Yang, H.Y., Tromp, J., et al. 2015, Geophys. J. Int., 202, 1392 [NASA ADS] [CrossRef] [Google Scholar]
 Lau, H. C., Faul, U., Mitrovica, J. X., et al. 2016a, Geophys. J. Int., 208, 368 [Google Scholar]
 Lau, H. C., Mitrovica, J. X., Austermann, J., et al. 2016b, J. Geophys. Res. Solid Earth, 121, 6991 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [Google Scholar]
 Levrard, B., & Laskar, J. 2003, Geophys. J. Int., 154, 970 [NASA ADS] [CrossRef] [Google Scholar]
 LonguetHiggins, M. S. 1968, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci., 262, 511 [NASA ADS] [CrossRef] [Google Scholar]
 LonguetHiggins, M. S., & Pond, G. S. 1970, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci., 266, 193 [NASA ADS] [Google Scholar]
 MacDonald, G. 1967, Proc. R. Soc. London Ser. A. Math. Phys. Sci., 296, 298 [Google Scholar]
 Matsuyama, I. 2014, Icarus, 242, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Matthews, K. J., Maloney, K. T., Zahirovic, S., et al. 2016, Global Planet. Change, 146, 226 [NASA ADS] [CrossRef] [Google Scholar]
 Maurice, M., Tosi, N., Schwinger, S., Breuer, D., & Kleine, T. 2020, Sci. Adv., 6, 28 [CrossRef] [Google Scholar]
 Mavromatis, H., & Alassar, R. 1999, Appl. Math. Lett., 12, 101 [CrossRef] [Google Scholar]
 Merdith, A. S., Williams, S. E., Collins, A. S., et al. 2021, EarthSci. Rev., 214, 103477 [NASA ADS] [CrossRef] [Google Scholar]
 Meyers, S. R., & Malinverno, A. 2018, Proc. Nat. Acad. Sci., 115, 6363 [NASA ADS] [CrossRef] [Google Scholar]
 Mignard, F. 1979, Moon Planets, 20, 301 [Google Scholar]
 Mojzsis, S. J., Arrhenius, G., McKeegan, K., et al. 1996, Nature, 384, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Motoyama, M., Tsunakawa, H., & Takahashi, F. 2020, Icarus, 335, 113382 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, M. 2008a, A Large Spectrum of Free Oscillations of the World Ocean Including the Full Ocean Loading and Selfattraction Effects (Springer Science& Business Media), 14 [Google Scholar]
 Müller, M. 2008b, Ocean Model., 20, 207 [CrossRef] [Google Scholar]
 Munk, W. H., & MacDonald, G. J. 1960, The Rotation of the Earth: A Geophysical Discussion (Cambridge University Press) [Google Scholar]
 Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
 Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
 Ooe, M. 1989, J. Phys. Earth, 37, 345 [CrossRef] [Google Scholar]
 Palmer, T., Shutts, G., & Swinbank, R. 1986, Q. J. R. Meteorol. Soc., 112, 1001 [NASA ADS] [CrossRef] [Google Scholar]
 Peck, W. H., Valley, J. W., Wilde, S. A., & Graham, C. M. 2001, Geochim. Cosmochim. Acta, 65, 4215 [NASA ADS] [CrossRef] [Google Scholar]
 Petit, G., & Luzum, B. 2010, IERS conventions (2010), Tech. rep. (France: Bureau International des Poids et mesures sevres) [Google Scholar]
 Platzman, G. W. 1983, Science, 220, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Platzman, G. W. 1984, J. Phys. Oceanogr., 14, 1532 [NASA ADS] [CrossRef] [Google Scholar]
 Proudman, J. 1920a, Proc. London Math. Soc., 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Proudman, J. 1920b, Proc. London Math. Soc., 2, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Regge, T. 1958, Il Nuovo Cimento (1955–1965), 10, 544 [CrossRef] [Google Scholar]
 Renaud, J. P., & Henning, W. G. 2018, ApJ, 857, 98 [Google Scholar]
 Riley, K. F., Hobson, M. P., & Bence, S. J. 1999, Mathematical Methods for Physics and Engineering [Google Scholar]
 Ross, M., & Schubert, G. 1989, J. Geophys. Res. Solid Earth, 94, 9533 [NASA ADS] [CrossRef] [Google Scholar]
 Rufu, R., & Canup, R. M. 2020, J. Geophys. Res. Planets, 125, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Sonett, C., & Chan, M. A. 1998, Geophys. Res. Lett., 25, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Sørensen, A. L., Nielsen, A. T., Thibault, N., et al. 2020, Earth Planet. Sci. Lett., 548, 116475 [CrossRef] [Google Scholar]
 Strauss, W. A. 2007, in Partial Differential Equations: An introduction (John Wiley& Sons) [Google Scholar]
 Sun, C., Xu, W., Cawood, P. A., et al. 2019, Sci. Rep., 9, 1 [Google Scholar]
 Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534 [Google Scholar]
 Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A&A, 630, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Touma, J., & Wisdom, J. 1994, AJ, 108, 1943 [NASA ADS] [CrossRef] [Google Scholar]
 Touma, J., & Wisdom, J. 1998, AJ, 115, 1653 [NASA ADS] [CrossRef] [Google Scholar]
 Touma, J., & Wisdom, J. 2001, AJ, 122, 1030 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706 [Google Scholar]
 Tyler, R. 2011, Icarus, 211, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Tyler, R. H. 2021, Planet. Sci. J., 2, 70 [CrossRef] [Google Scholar]
 Vallis, G. K. 2017, Atmospheric and Oceanic Fluid Dynamics (Cambridge University Press) [CrossRef] [Google Scholar]
 Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1988, Quantum Theory of Angular Momentum (World Scientific) [CrossRef] [Google Scholar]
 Walker, J. C. G., & Zahnle, K. J. 1986, Nature, 320, 600 [NASA ADS] [CrossRef] [Google Scholar]
 Waltham, D. 2015, J. Sediment. Res., 85, 990 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, H., Boyd, J. P., & Akmaev, R. A. 2016, Geosci. Model Dev., 9, 1477 [Google Scholar]
 Watterson, I. G. 2001, J. Atmos. Oceanic Technol., 18, 691 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. 1973, in Deep Sea Research and Oceanographic Abstracts (Elsevier), 20, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. 1980, Geophys. J. Int., 61, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. 1982, Geophys. J. Int., 70, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Wilde, S. A., Valley, J. W., Peck, W. H., & Graham, C. M. 2001, Nature, 409, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, G. E. 1990, J. Phys. Earth, 38, 475 [CrossRef] [Google Scholar]
 Williams, G. E. 1997, Geophys. Res. Lett., 24, 421 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, G. E. 2000, Rev. Geophys., 38, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, J. G., & Boggs, D. H. 2016, Celestial Mech. Dyn. Astron., 126, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Wood, R., Liu, A. G., Bowyer, F., et al. 2019, Nat. Ecol. Evol., 3, 528 [CrossRef] [Google Scholar]
 Zahel, W. 1980, Phys. Earth Planet. Inter., 21, 202 [CrossRef] [Google Scholar]
 Zahnle, K. J., Lupu, R., Dobrovolskis, A., & Sleep, N. H. 2015, Earth Planetary Sci. Lett., 427, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Zhong, Y., Wu, H., Fan, J., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 540, 109520 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Orbital dynamics
For the reconstruction of the EarthMoon distance, we used a reduced secular dynamical model describing the exchange of angular momentum between the Earth’s rotation and the lunar orbital motion, ignoring the Earth’s obliquity, lunar eccentricity, and lunar inclination, and mainly focusing on terrestrial tides. This simplification allows for a systematic understanding of the hierarchically complex contributions of multiple intervening players. In particular, the contribution of the eccentricity tides becomes significant when the orbit of the Moon was highly eccentric. Accounting for lunar tides and lunar coremantle boundary dissipation would counteract the effect of terrestrial tides and increase the lunar eccentricity when going backwards in time, but only to moderate values (e ≤ 0.1) (Daher et al. 2021). By attaining a highly eccentric lunar orbit, is possible through evection resonance^{1}. Touma & Wisdom (1998) studied this regime and showed that capture into such a resonance could have been encountered for a_{M} ≈ 4.6R_{E}, exciting the lunar eccentricity to e ≈ 0.5. However, the timescale of capture and escape from this evection resonance is 10^{4} ∼ 10^{5} yr, after which the lunar orbit tends to circularization again (Rufu & Canup 2020).
Dissipation within the moon is also rendered significant at the early stage of the system with the Earth fully molten and the moon having little to no atmosphere, forcing it to quickly cool into a highly dissipative body. However, this also occurs over a relatively short time interval – specifically, when a_{M} < 20R_{E} (Zahnle et al. 2015). The timescale of these mechanisms is much smaller than that associated with the longterm tidal evolution that we model in this work. Furthermore, a key feature of the lunar distance evolution is the runaway effect, with a_{M} dropping rapidly from 30R_{E} to the formation site, slightly beyond the Roche limit, within few million years (see Figure 3). The same figure shows that the Moon spends 97% of its lifetime with a_{M} > 30R_{E}. Thus, as much as early stage mechanisms are essential to constrain the formation scenarios of the system, the robustness of the runaway evolution beyond 30R_{E} renders the reduced dynamical model a safe and sufficient approach for the longterm study. This model should provide the skeleton of the secular evolution in the system around which full spatial dynamics can flesh; the latter would require us to extend the oceanic tidal model also to the obliquity component, which could be the task of a next stage of this work. Other effects such as climate friction (Levrard & Laskar 2003) and coremantle coupling (Neron de Surgy & Laskar 1997; Touma & Wisdom 2001) are ignored, as well as halts of tidal interaction due to Laplace plane transitions (Ćuk et al. 2016).
Under these assumptions, the governing dynamical system of equations is expressed as:
where 𝒯_{M} is the lunar semidiurnal tidal torque coupling between the oceanic and the solid response of the Earth, and 𝒯_{S} is its solar counterpart. The orbital angular momentum of the Moon , where β = M_{E}M_{M}/(M_{E} + M_{M}) is the EarthMoon system’s reduced mass. The rotational angular momentum of the Earth is defined as L_{Ω} = C(Ω)Ω, with the timevarying principal moment of inertia given by (Goldreich 1966)
Here, is the seconddegree fluid Love number of centrifugal/tidal deformation and G is the gravitational constant. The differential equation is integrated backwards in time using the RungeKutta 9(8) method, starting from the present and stopping at a_{M} = 3R_{E}. The tidal torque computation is coupled to the orbital integrator and is computed simultaneously at each step. It takes the model parameters (H, σ_{R}) as input, and the system’s variables, a_{M} and Ω, to compute the tidal frequency and, consequently, the coupled tidal response.
Once the lunar semimajor axis (a_{M}) and the rotation speed of the Earth (Ω) are determined, we compute the obliquity of the Earth (ϵ) and the precession frequency (p) as derived quantities (Laskar et al. 2004). Starting with Equations (40) and (46) from Correia & Laskar (2010) in the case of zero eccentricity, we obtain:
and
that is
We note that the tidal response parameter, K, disappears from the equations. This would also be the case if K depended on Ω. The obliquity evolution equation (A.6) is integrated using the values of a_{M} and Ω that result from the tidal flows and orbital dynamics coupled system. The precession frequency, p, is then derived using Equations (6) and (8) from Laskar et al. (2004) with zero eccentricity and inclination, that is:
In (A.6) and (A.7), the constant values taken for the Earth’s radius R_{E}, the gravitational constant of the Moon GM_{M}, and the Sun GM_{S}, the mass ratio M_{E}/M_{M}, the rotational velocity Ω_{0}, the Earth’s semimajor axis a_{E}, and the inertia parameter are adopted from INPOP21 (Fienga et al. 2021). The dynamical ellipticity at the origin of date, E_{d}(Ω_{0})=0.003243, is determined from the initial conditions for the obliquity (ϵ_{0}) and precession (p_{0}) adopted from the La2004 solution (Laskar et al. 2004). All the values of applied parameters are summarized in Table D.3.
This derivation of the obliquity and precession frequency evolutions is only valid in the limit of a distant Moon, that is, when the Moon is beyond its Laplace radius (Boué & Laskar 2006; Farhat & Touma 2021) and its Laplace plane is the ecliptic rather than the Earth’s equatorial plane (Tremaine et al. 2009). In our a_{M} evolution of Figure 3, the Laplace regime transition occurs very early in the evolution (t > 4Ga); thus, in Figure 6, we plot the evolution of the precession frequency and obliquity between the present and 3.5 Ga. We also add the geological inferences of the precession frequency as scatter on the curve; in the case of cyclostratigraphy, this is the direct observable (Table D.1).
Appendix B: Continental drift and oceanic geometry shifting
Our dynamical integrator allows for variations in the oceanic geometry, be it a variation in the position of the oceanic hemisphere or a shift between the hemispheric and global oceanic configurations. In the combined model, the first phase (Figure 4) starts with the center of the continental cap following the evolution of the geographic center over the recent billion years. For this purpose, we adopted a recent model that reconstructs a kinetically continuous history of plate tectonics (Merdith et al. 2021). The geographic center is traced by computing the surface projection of the “barycenter” of the continental distribution. This allows for a higher level of realism in oceanic modeling. Beyond 1 Ga (in part due to the lack of plate tectonic data), the model continues with the position of the ocean at 1 Ga. To postprocess the continental drift evolution and to produce the timesliced sketches of Figure 1, we used the GPlates opensource reconstruction software (Boyden et al. 2011; Gurnis et al. 2012).
Fig. B.1. Drifting effect of the continental cap on the oceanic response: the tidal torque of a hemispheric ocean is plotted as a function of the forcing semidiurnal frequency for different positions of the center of the ocean. With longitudinal symmetry, the latter is defined by the latitude of the oceanic center, which evolves according to Figure 1. The drifting effect on the resonances ranges from position shifting and attenuation for small forcing frequencies to major distortion in the spectrum at larger frequencies. Extreme distortion occurs in the polar oceanic scenario: the major resonance around 11 rad/day reaches a maximum relative to other configurations and the rest of the resonances are absorbed into the background leaving a unimodal spectrum. This behavior makes it important to take into account the position of the hemispherical cap into the model (Figure 1). 
At 1.5 Ga, the integrator starts computing simultaneously the tidal response of a global oceanic geometry, with uniform thickness H_{global} = H/2, in order to guarantee oceanic volume conservation when we switch between the geometries. However, the hemispherical response remains the one that is accounted for in the dynamical evolution. While simultaneously computing both, the code detects when they equate, and switches to the global configuration identifying this time as t_{switch}. The physical outcome of this process is guaranteeing a better compliance with continental crust growth curves (Dhuime et al. 2012; Hawkesworth et al. 2020), thus avoiding effects arising from blocking of westward tidal propagation or enhanced continental dissipation at continental shelves (Arbic et al. 2009). The mathematical outcome of this process is evident in Figure 4 in terms of guaranteeing a smooth dynamical evolution of a_{M} (Figure 3) without any discontinuities and modeling artifacts. For the misfit minimum of our combined model, we have t_{switch} = 3.25 Ga.
Appendix C: Parameter fits
To construct the misfit surfaces of Figure 2, we computed the evolution for each pair of (H, σ_{R}) on the twodimensional grid. The present rate of lunar recession, ̇a_{0}, and the impact time, t_{f}, are then extracted for each evolution sample and the mean square weighted deviation, χ^{2}, (Table 1) is then computed as:
where we use Lunar Laser Ranging (LLR) estimates of lunar orbital recession (Williams & Boggs 2016): mm/year; as well as geochemical estimates of lunar formation time (Maurice et al. 2020): Ga. The maximum likelihood detection problem is further optimized by fitting the surface around the minimum by a parabola to avoid a limitation of the grid resolution.
Misfit analysis summary.
We evaluated the uncertainties on the fitted parameters from those on the observables following the standard propagation of uncertainty method. Because of the absence of correlation between the two data ȧ_{0} and t_{f}, the entries of the variance matrix,
are given by
The partial derivatives entering in these formulae are computed numerically from the fit of (ȧ_{0} ± σ^{LLR},t_{f}) and (ȧ_{0},t_{f} ± σ^{geo}). The marginal uncertainties on parameters H and σ_{R} are σ_{H} = var(H)^{1/2} = 32.75 m and σ_{σR} = var(σ_{R})^{1/2} = 0.2631 × 10^{−6} s^{−1}, respectively. We used the variance matrix Σ to evaluate the 2σconfidence ellipsoid around the bestfit parameters (H, σ_{R}). The EarthMoon distance, a_{M}, and the length of the day LOD
have been integrated for 25 pairs of parameters (H, σ_{R}) chosen at the boundary of this 2σconfidence region. Their envelope represents the 2σuncertainty area plotted in shaded blue in (Figure 3) and (Figs. 5 and 6).
Appendix D: Geological data
In tables D.1 and D.2, we compile geological data sets that provide historical snapshots of the past rotational state of the Earth and the lunar orbital distance.
Cyclostratigraphic data.
Tidal rhythmites data.
Values of constant parameters used in the numerical implementation of the theory.
Appendix E: Tidal response of a hemispherical ocean
This appendix is aimed at obtaining a computation of the tidal response of a hemispherical ocean on the surface of the Earth. The formalism is heavily based on earlier works (LonguetHiggins & Pond 1970; Webb 1980) describing the free oscillations and the tidal response of a hemispherical ocean symmetric about the equator; we expand upon it here by adopting the true polar wander scenario (Webb 1982) to solve for a general oceanic position. We note that the mathematical formulation of the referenced works (LonguetHiggins & Pond 1970; Webb 1980, 1982) contains several misprints that we correct here.
In the coplanar problem under study (ignoring the Earth’s obliquity and the lunar orbital inclination), we define a frame of reference corotating with the Earth, with a spin vector of , Ω being the Earth’s spin rate and ẑ as the unit vector along its figure axis. In this frame, we used the spherical coordinates (r, θ, λ) denoting the radius, the colatitude, and the longitude respectively, and their corresponding unit vectors . We start with the linearized system of equations that describe the conservation of momentum and mass in a tidally forced shallow oceanic layer (Matsuyama 2014):
where is the horizontal velocity field, g is the gravitational acceleration at the surface, ζ is the oceanic depth variation, ζ_{eq} is the equilibrium depth variation, H is the uniform oceanic thickness (the first of only two free parameters in our model), and σ_{R} is the Rayleigh (or linear) drag frequency (Matsuyama 2014; AuclairDesrotour et al. 2018); the latter is an effective dissipation parameter characterizing the damping of the oceanic tidal response by dissipative mechanisms (the second free parameter in our model). On Earth, σ_{R} mainly accounts for the conversion of barotropic tidal flows into internal gravity waves, which represents nearly 85% of the total dissipation for the actual lunar semidiurnal oceanic tide (see, e.g., Carter et al. 2008). For this mechanism, the Rayleigh drag frequency can actually be related to physical parameters such as the BruntVäisälä frequency, which quantifies the stability of the ocean’s stratification against convection (see, e.g., Gerkema & Zimmerman 2008), or the lengthscale of topographical patterns at the oceanic floor (Bell 1975; Palmer et al. 1986). In Eq. (E.1), the Coriolis parameter f is given by:
the horizontal gradient operator ∇ is defined as:
and the horizontal divergence of the velocity field ∇ ⋅ u as:
with R being the Earth’s radius. Finally, we remark that the interaction of tidal flows with the mean flows of the oceanic circulation are ignored in the momentum equation, namely, Eq. (E.1a).
For u = ∂_{t}x, where x is the horizontal tidal displacement field, we have:
Following Proudman (1920b), we use Helmholtz’s theorem at this step (e.g., Arfken & Weber 1999, Chapter 1) to decompose the horizontal displacement vector field into:
where ∇Φ is a curlfree vector field (∇×(∇Φ) = 0) and is a divergencefree vector field (). In the above equation, we introduce the divergent displacement potential Φ and the rotational displacement streamfunction Ψ (e.g., Gent & McWilliams 1983; Webb 1980; Tyler 2011), with the latter accounting for the vortical component of the tidal displacement field (e.g., Vallis 2017). As discussed by FoxKemper et al. (2003), while the Helmholtz decomposition is unique for infinite domains, this is not true for bounded domains such as hemispherical oceanic shells due to lack of additional physical constraints on the boundary conditions for either of the components of the sum. There are boundary conditions only on the total flux at coastlines. Impermeability is a typical boundary condition: the net flux normal to the coast is zero, which is expressed as x · n̂ = 0, where n̂ designates the outward pointing unit vector defining the normal to the coast. Following Webb (1980, 1982), we assume that both components of the flux satisfy this condition, namely:
We note that the second condition of the above equation can be rewritten as , which implies that Ψ is constant along the coastline (in the following, we set Ψ = 0 at the oceanic boundary). This condition thus means that the coastline corresponds to a closed streamfunction contour, which depicts a distinct gyre of the tidal flow.
Although arbitrary, the assumption that both components of the flux satisfy the impermeability condition has been widely used to study the dynamics of ocean basins because of its convenience relative to other possible conditions (e.g., Gent & McWilliams 1983; Watterson 2001; Han & Huang 2020). In particular, this assumption provides a unique decomposition apart from an arbitrary additive constant to each function, Φ and Ψ. Moreover, the second condition given by Eq. (E.7) enforces the orthogonality of the curlfree and divergencefree components of the tidal flow. By combining together the identity and Gauss’ theorem (e.g., Arfken & Weber 1999):
with dA and dℓ being the infinitesimal area element of the hemispherical oceanic domain, 𝒪, and length element of the coastline, ∂𝒪, respectively, we obtain
As the second condition of Eq. (E.7) enforces along the coastline, it follows that
meaning that the components ∇Φ and each belong to one of the two orthogonal subspaces that form the space of horizontal displacements satisfying the assumed boundary conditions. We remark that the orthogonality of the Helmholtz decomposition is not necessarily verified in the general case since it is itself a consequence of the specific boundary condition chosen for the divergencefree component of the tidal flow.
The functions Φ and Ψ are expanded in terms of complete sets of eigenfunctions over the domain, 𝒪, such that
The eigenfunctions (ϕ_{r}, ψ_{r}) satisfy, over the oceanic domain (𝒪), the Helmholtz equations (e.g., Riley et al. (1999), Chapter 21):
and, along the coastline (∂𝒪), the boundary conditions given by Eq. (E.7):
where we introduced the horizontal Laplacian:
and the real eigenvalues, μ_{r} and ν_{r}, associated with the eigenfunctions ϕ_{r} and ψ_{r}, respectively. We note that the eigenfunctions are normalized such that
where the notation δ_{rs} referring to the Kronecker δsymbol is δ_{rs} = 1 for r = s and 0 otherwise. Using these conditions, the eigenfunctions are defined as:
with eigenvalues of μ_{r} = ν_{r} = n(n + 1)/R^{2} and the normalization coefficient of:
In Eqs. (E.18) and (E.19), each harmonic index, r, of the eigenfunctions is associated with a degree, n, and order, m, and the expansion functions are the associated Legendre functions (Abramowitz et al. 1988). In the definition of ϕ_{r} Eq. (E.18), n ∈ ℕ and m = 0, 1, ..., n while in the expression of ψ_{r} Eq. (E.19), n ∈ ℕ^{*} and m = 1, 2, ..., n. By convention, we set ψ_{0} = 0 hereafter. Figure E.1 shows the eigenfunctions ϕ_{r} and ψ_{r} for 1 ≤ m ≤ n ≤ 4 and the streamlines of the associated tidal flows.
Fig. E.1. Eigenfunctions ϕ_{r} (left) and ψ_{r} (right), and the associated tidal flows. The eigenfunctions defined by Eqs. (E.18) and (E.19) are plotted over the hemispherical oceanic domain for 0 ≤ n ≤ 4 (from top to bottom) and 0 ≤ m ≤ n (from left to right). Bright or dark colors designate positive or negative values of the eigenfunctions, respectively. Streamlines indicate the tidal flows corresponding to ∇ϕ_{r} for the set {ϕ_{r}} and to for the set {ψ_{r}}. 
The eigenfunctions (ϕ_{r}, ψ_{r}) can be split into two sets describing tidal solutions that are symmetric or antisymmetric about the equator, and thus one can decide, based on the symmetry of the tidal forcing, on the associated set of eigenfunctions that need to be considered using a classification scheme (LonguetHiggins 1968; LonguetHiggins & Pond 1970) for the pairs (n, m). However, in our model, where the ocean is no longer symmetric about the equator, both symmetric and antisymmetric eigenfunctions are required. In substituting the definitions of Eqs. (E.11), (E.12), and (E.13) into the continuity equation (E.1b), we find that
What is left to complete the solution is finding the coefficients p_{r} and p_{−r} by substituting the series expansions in the momentum equation (E.1a) and multiplying by ∇ϕ_{r} and , then integrating over the oceanic area. Starting with the former, we get
The product of the gradients of two eigenfunctions is computed using Green’s first identity (e.g., Strauss (2007), Chapter 7):
The first term on the righthand side vanishes as it includes the boundary condition at the coast (Eq. E.15). The second term is computed using the eigenvalue equation (Eq. E.13) and the normalization condition, thus:
Rearranging the other terms using vector identities, we can write Eq.(E.22) as:
The third term vanishes due to orthogonality (see Eq. (E.10)) and upon replacing the Coriolis term by its definition in Eq. (E.2), we are left with
To close the system, we multiply the momentum equation with to get
Integrating Eq. (E.27) over the area of the ocean and using basic vector product identities, we get
where upon replacing the Coriolis term as before we have
We identify in Eqs. (E.26) and (E.29) the socalled “gyroscopic coefficients" (e.g., Proudman (1920a,b)) that are defined as
with β_{−s, r} = −β_{r, −s}. These coefficients carry the effect of rotational distortion to the tidal waves and the boundary conditions imposed by the coastlines. Using these definitions, Eqs. (E.26) and (E.29) form an infinite linear system in the coefficients p_{r}(t) and p_{−r}(t) that is expressed as:
This system will be transformed to the frequency domain (F), and then truncated and solved spectrally as a function of the tidal forcing frequency. However, we shift here to extend the theory to the effects of selfattraction and loading between the ocean and the deforming mantle in order to have a complete selfconsistent tidal response of the Earth.
Appendix F: Coupling the hemispheric oceanic response with solid deformation
In the tidal theory under study, the solid component of the Earth is subject to viscoelastic deformation as a result of three contributions: the direct tidal effect of the tidal perturber, the loading effect of the perturbed oceanic shell, and the effect of gravitational selfattraction between the oceanic shell and the solid part (Farrell 1972; Zahel 1980). If we were to take these into account when studying oceanic tides, ζ becomes a function of two moving surfaces: the free oceanic surface, ζ_{os}, and the vertically deforming solid surfacem ζ_{ss}.
In the frame corotating with the Earth as defined in Appendix E, the gravitational potential is expressed as (e.g., AuclairDesrotour et al. 2019)
where G is the gravitational constant, M is the mass of the tidal perturber (the Sun or the Moon), and r′ is the distance between the Earth and the perturber. In the shallow ocean approximation, the tidal potential at the Earth’s surface (r = R) is:
where a constant offset was removed as it does not contribute to the tidal force. In the frequency domain, the tidal potential U^{T} is expanded spectrally in Fourier series and spatially in spherical harmonics, with complex coefficients , as (Kaula 2013; AuclairDesrotour et al. 2019):
where s is an integer and the tidal forcing frequency , the frequency n_{orb} being the orbital mean motion of the tidal perturber. In the absence of obliquity, the n^{th} harmonic of the tidal potential is given by (Ogilvie 2014):
where a is the semimajor axis of this perturber and A_{n, m, s}(e) are dimensionless functions of the orbital eccentricity of the perturber e computed via the Hansen coefficients (Laskar 2005; Correia et al. 2014). In our study, we restricted the tidal potential to the dominant contribution of the semidiurnal component identified by n = m = s = 2 and corresponding to the tidal frequency . For this component, while neglecting the small orbital eccentricity of the Sun and the Moon, . Hereafter, we use to represent a single harmonic (n, m, s) of the tidal potential. This harmonic of degree n is defined as
Moreover, in the following, we write σ instead of to simplify the notation. Subject to only, the equilibrium oceanic depth would be . However, the loading effect of the deforming oceanic shell adds to the tidal potential and they both affect the ocean surface and the ocean floor corresponding to the solid surface . The former is expressed as (Matsuyama 2014):
where ρ_{oc} and ρ_{se} stand for the uniform oceanic and solid Earth densities, respectively. In this equation, the oceanic depth variation, ζ, is decomposed into spherical harmonics defined over the full sphere as:
Although ζ given in (E.22) is only defined over the oceanic hemisphere, this decomposition over the whole sphere is required when applying the Love numbers. Using the orthogonality of spherical harmonics, and the fact that ζ(θ, λ, t)=0 over the continental hemisphere, Eq. (F.7) can also be expressed as:
where the integral is computed over the solid angle 2π spanned by the ocean. The second contribution to the equilibrium tide, which is due to the solid redistribution of mass, is
In equations (F.6) and (F.9), we used the tidal Love numbers and , and the surfaceloading Love numbers and , where the first of each set is the transfer function corresponding to the gravitational response, and the second codes for the vertical displacement. We emphasize that here the Love numbers are defined in the Fourier domain, therefore, they correspond to the intrinsic mechanical impedances of the solid part that relate its viscoelastic tidal response to tidal forcings in the permanent regime and they characterise both the elastic deformation of the body and its anelastic deformation resulting from energy dissipation due to viscous friction in the Earth’s interior. In the general case, the four Love numbers (, and ) can be computed from internal structure models (e.g., Tobie et al. (2005, 2019), Bolmont et al. (2020)). In the present study, for the sake of simplicity, we use the closedform solutions derived for a uniform solid interior (Munk & MacDonald 1960):
where is a complex dimensionless effective shear modulus, with a form that is dependent on the chosen solid rheology (Efroimsky 2012; Renaud & Henning 2018). To specify , we consider an Andrade rheology (Andrade 1910; Bagheri et al. 2019), which has an advantage over the commonly used Maxwell rheology in attenuating the rapid decay of the anelastic component of the deforming solid Earth for high tidal frequencies (CastilloRogez et al. 2011; AuclairDesrotour et al. 2019). This is particularly useful with regard to avoiding an overestimation of the tidally dissipated energy of the solid part during early eons. For this rheology, takes the following form (Findley et al. 1977; Efroimsky 2012):
where M_{E} is the mass of the Earth, μ_{E} its average rigidity, and Γ is the gamma function (Abramowitz et al. 1988); α_{A} is a dimensionless rheological exponent determined experimentally (Castelnau et al. 2008; Petit & Luzum 2010); τ_{A} is the anelastic Andrade timescale; and τ_{M} the Maxwell relaxation time defined as the ratio of viscosity to rigidity. For a volumetric average of the mantle’s shear modulus μ_{E} = 17.3 × 10^{10} Pa, and volumetric average of viscosity deduced from inversions of Lau et al. (2016b), we have τ_{M} = 685 yrs. The values α_{A} = 0.25 and τ_{A} = 2.19 × 10^{4} yrs that we use in our model are adopted from AuclairDesrotour et al. (2019). All the applied parameter values are summarized in Table D.3.
Taking the effect of solid Earth deformation into account, we replace the equilibrium tide ζ_{eq} in the momentum equation (E.1a) by the difference of Eqs.(F.6) and (F.9) and we resolve it in the Fourier domain using the forcing tidal frequency, σ. The modified momentum conservation equations is now expressed as:
with , and where the loading and tidal tilt factors are defined as (Matsuyama 2014):
Just like the Love numbers, and are complex in the general case and tend toward unity as the deformability of the solid and oceanic layers decreases. Now we get to the added contribution of the oceansolid coupling to the linear system of p_{r} and p_{−r}. Multiplying the added contribution of loading and selfattraction effects by ∇ϕ_{r} and , then resolving the added terms in the frequency domain, after a number of manipulations, we can finally rewrite the expression of the system of Eqs. (E.31) and (E.32) as:
where we define
with corresponding to the Gram matrix of the ALFs,
for which the method of computation is detailed in Appendix H. Coupled to the orbital dynamical evolution, the tidal frequency σ is determined at each point in time in the hemispherical phase of the model, then the system is truncated at r_{max} and solved numerically (see Appendix J). We rewrite the linear system as:
where the first equation is for r > 0 and the second is for r < 0 – and where we have introduced the coefficients:
Appendix G: Gyroscopic coefficients
The gyroscopic coefficients introduced in Eq. (E.30) characterize the rotational distortion of tidal waves via the Coriolis force term and the effect of boundary conditions imposed by the oceanic geometry. This coupling is dependent on the position of the ocean on the sphere and the relative position of the tidal perturber with respect to the tidally forced ocean. Since we are after a generic configuration describing the response of the oceanic hemisphere at any position, the expressions of Eq. (E.30) should be written for any frame rotating with the ocean. We start with the definition of the ALFs (Chapter 8 of Abramowitz et al. (1988)):
which are solutions to the Legendre equation,
Following differentiation, we obtain
Substituting Eq. (G.3) in Eq. (G.2) we get the recurrence relation
which gives the useful relation
From Eqs. (G.1G.5), it is straightforward to obtain the ALFs recurrence relations that are necessary to compute the integral equations of the gyroscopic coefficients,
Fig. G.1. Adopted transformation scheme that allows recovering the tidal response of a hemispheric ocean with an arbitrary center on the sphere. We use an Eulerian transformation of the form ℛ_{3}(α)ℛ_{2}(β)ℛ_{3}(γ) with γ = 0, allowing us to shift the latitude of the oceanic center O by shifting the spin axis from ŝ to ŝ^{″} in a true polar wander scenario (Webb 1982). 
The theory of the hemispherical tidal response is based on an ocean bound by two meridians. Thus for an oceanic center moving on the sphere, we instead rotate the spin axis relative to the center of the ocean, and accordingly the frame of the tidal perturber to maintain the coplanar configuration of the dynamical system. These rotations will enter the system through the Coriolis term, specifically through the gyroscopic coefficients, along with the tidal forcing term. We define an arbitrary rotation {θ, λ}→{θ′,λ′} using an Eulerian rotation matrix of the form ℛ_{3}(α)ℛ_{2}(β)ℛ_{3}(γ), with (0 ≤ α ≤ 2π) and (0 ≤ β ≤ π), and we fix γ = 0 (see Fig. G.1). For a vector J defined as:
the transformed gyroscopic coefficients are:
where, for r associated with the harmonic pair of integers (n, m) and s associated with (u, v), we introduce the coefficients:
with and
Under the terms of this transformation, the latitude of the center of the ocean in the rotating frame is given by:
To compute the integrals involved in the gyroscopic coefficients, we make use of the essential condition^{1} (e.g., LonguetHiggins & Pond (1970))
and we use the overlap integral of two ALFs (Eq. F.17), which we compute using the closed form relations provided in the following section. Now we have at hand all the elements to compute the gyroscopic coefficients harmonically. The final form of the three coefficients with superscript ^{(1)} are identical to those in Webb (1980) and similar to those in LonguetHiggins & Pond (1970) up to certain misprints. For the rest of the terms, the expressions given in Webb (1982) involve numerous typographical errors and inconsistencies, so we provide here the full set of the gyroscopic coefficients. The coefficients and are expressed as:
where we used Eq. (G.8) for each integrand in Eq. (G.31a) to obtain Eq. (G.31b). The coefficients , , and read as:
where we used the recurrence relations of Eq. (G.8) and Eq. (G.9 to compute the first integral of Eq. (G.34b), and the relation of Eq. (G.7) to compute the second integral. Finally, the remaining terms , , , and are expressed as:
Appendix H: Overlap integral
Here, we provide a closed form solution for the computation of the overlap integral of Eq. (F.17). The procedure is assimilated from tools of angular momentum quantization (Varshalovich et al. 1988). Following Dong & Lemus (2002) and introducing the notation q = v − m, we have:
where the factors are given by:
and the coefficients 𝒟(q,l) by
We note here that the phase κ introduced in Dong & Lemus (2002) as
corrects the phase given in Mavromatis & Alassar (1999) and Crease (1966), where the latter was used for the computation of the gyroscopic coefficients in LonguetHiggins & Pond (1970) and Webb (1980, 1982).
In Eq. (H.1), the summation over l runs for n − u≤l ≤ (n + u);l ≥ q; and l + n + u is even. Finally, the Wigner 3jm symbols are determined from Varshalovich et al. (1988) by
where J = a + b + c, and R_{ij} are the elements of the socalled Regge ℜsymbol (Regge 1958) defined as
The summation in Eq. (H.5) runs over all integer values of z for which all the factorial arguments are nonnegative. Finally, we note that using Eq. (H.1), when v = m. In that case, we alternatively use
This method for the computation of the overlap integral was verified numerically using MATLAB’s ALFs package.
Appendix I: The tidal forcing term
As in Webb (1980), considering the equilibrium tide to have a unit root mean square amplitude, and to be driven by a spherical harmonics term:
with angular frequency σ, we have:
Under the rotation of the coordinate system described by the Euler angles (α, β, γ) (see Appendix G and Fig. G.1), the spherical harmonics transform as (Varshalovich et al. 1988):
or
where designate the Wigner Dfunctions. These functions are themselves the product of three functions (Varshalovich et al. 1988), each depending on one argument: α, β, or γ,
In this expression, the functions are given by
with j running over all integer values for which the factorial arguments are positive. This sum involves N + 1 terms, where N is the minimum of (p + q),(p − q),(p + s), and (p − s). Since we are studying the semidiurnal tide (p = q = 2), we are left with one term only. Expanding the harmonic factor of Eq. (I.4) in terms of the basis eigenfunctions, we get the expression of the equilibrium oceanic depth variation in the rotated frame of reference,
Then, invoking the definition of the component ,
we get its expression in the rotated frame of reference,
where the dot products of the eigenfunctions are simplified as:
and
We note that as the index s takes negative values, we use
in the overlap integral of Eq. (F.17).
Appendix J: Tidal torque of a hemispherical ocean
Once the gyroscopic coefficients are computed, the linear system of the coefficients (p_{r}, p_{−r}) in Eq. (F.18) is truncated and solved numerically. What we are after is the tidal torque that enters in the dynamical equations of the EarthMoon system. Two torques are involved as explained in the main text, and they depend on the rotational angular velocity of the Earth and the orbital frequency of the tidal perturber. Defining the tidal torque transferring power from the Earth’s rotational momentum to the perturber’s orbital angular momentum by 𝒯, the power lost by the Earth would be 𝒯Ω, and the power gained by the perturber is 𝒯n_{orb}. The difference between them is the dissipative work of the total tidal mass redistribution 𝒲_{diss}, thus
The total dissipative work is the sum of two contributions: the dissipative work of oceanic tidal currents, , and dissipation in the deforming viscoelastic mantle. In the formalism established thus far, we calculated the selfconsistently coupled tidal responses of the ocean and solid part for the Earth’s half hosting the hemispherical ocean, which corresponds to the effective tidal response of the planet for this hemisphere. The tidal response of the continental hemisphere is simply described by the solid Love numbers introduced in Eq. (F.10), since there is no oceanic tide in that case. The coupled solidoceanic tidal response accounts for both the direct gravitational tidal forcing generated by the perturber on the solid part and ocean and the mutual forcings of the two layers through the variations of the loading exerted by the latter on the former, and the variations of the Earth’s selfgravitational potential due to mass redistribution. For simplicity, we ignore the energy dissipated in the solid part in the calculation of the tidal torque and we only consider that which is occurring within the oceanic shell, namely, . This is justified by the predominance of the oceanic response over the solid part over the time interval covered by the hemispherical ocean configuration in our model. This hierarchy of contributions is only jeopardized by the emerging significance of the solid dissipation when moving backwards in time and increasing the Earth’s rotational velocity Ω. Solid Earth dissipation would also be amplified with an early less viscous mantle due to higher HadeanArchean temperatures (Ross & Schubert 1989). Eventually, a regime transition may lead to the predominance of the mantle’s elastic response (Lau et al. 2015, 2016a). In our nominal model of the main text, the switch from the hemispherical ocean configuration to the global ocean configuration occurs midArchean, beyond which we selfconsistently account for the dissipative contribution of the solid part (Appendix K). Thus, we have only ignored the dissipative contribution of the mantle when it is insignificant.
Fig. J.1. Numerical analysis on the dependence of the tidal response computation on the truncation order r_{max}. The response is quantified by the root mean square tidal amplitude, ζ_{rms}, (Eq. J.5) and the dissipative work, , (Eq. J.2), and plotted for three tidal frequencies: 7.3, 11.4, and 22 rad/day that correspond to the vicinity of a tidal resonance, the peak of a resonance, and the background spectrum respectively. 
The oceanic dissipative work is given by (Webb 1980):
where ⟨⟩ denotes time averaging over the tidal period. This work should be equal to the work done by the tidal force on the ocean:
Hence, the tidal torque associated with the lunar semidiurnal frequency σ = 2(Ω − n_{M}), n_{M} being the lunar mean motion, is
and we obtain a similar expression for the solar tides 𝒯_{S} when solving the system with the solar tidal frequency component σ = 2(Ω − n_{S}), n_{S} being the solar mean motion that generates the solar tidal work.
Besides the tidal torque, the tidal response can also be quantified by the root mean square tidal height variation, ζ_{rms}, given as
As these quantities are computed numerically, a truncation order, r_{max}, is required. In Fig. J.1, we show the numerical dependence of the tidal response on r_{max}. Since the response is dominated by gravity modes, the tidal solution converges fast enough with r_{max}. To avoid any truncation effect in our computation, and to properly account for the resonances, we adopted r_{max} = 50.
Appendix K: Modeling the tidal response of a global ocean
When the global oceanic geometry is encountered in our model, the tidal response is computed based on the analytical formalism described in AuclairDesrotour et al. (2018, 2019). We refer to these references for a complete development of the theory and we only offer a brief here on the essential steps that lead to a computation of the tidal response described in the main text^{3}. In this approach, solving the governing system in Eq. (E.1) is done by expanding the velocity field, the tidal elevation, and the forcing gravitational tidal potential in Fourier series of time and longitude, with the tidal frequency serving as the expansion frequency. Thus, we have:
Defining the complex tidal frequency and the complex spin parameter as in AuclairDesrotour et al. (2018) by:
and replacing the tidal quantities by their expansions, the governing system reduces to an eigenvalueeigenfunction problem, known classically (when ignoring friction) as the Laplace tidal equation (Lee & Saio 1997). We assume that the Fourier components can be expanded spatially using a set of the latitudinal complex Hough functions (Hough 1898) {}, associated with a set of eigenvalues {}. To compute these functions and their associated eigenvalues, we adopt the method developed in Wang et al. (2016), where Hough functions are expanded in terms of Associated Legendre Functions
with and being complex change of basis coefficients. Using the change of basis coefficients , the tidal displacement solution is expressed as:
where the components are solutions to the linear algebraic system
Fig. K.1. Tidal torque between the Earth and the Moon corresponding to the coupled oceanicsolid response in two configurations: a global oceanic shell of thickness H = 4000 m (shown in red), and a hemispherical ocean with the same thickness, symmetric around the equator and bounded by longitudes λ = 0 and λ = π (in blue). Energy dissipation is quantified by the linear Rayleigh drag frequency σ_{R}. The logarithm of the torque is plotted as a function of the normalized frequency ω = (Ω − n_{orb})/Ω_{0}, where the Earth’s spin rate varies with the tidal forcing frequency Ω = n_{orb} + σ/2 at fixed n_{orb}, and Ω_{0} being the present spin rate of the Earth. 
Fig. K.2. Similar to Fig.K.1, but comparing the torque of the hemispherical ocean model between a pure oceanic response, and the response of the ocean when accounting for loading and selfattraction effects arising from solid Earth deformation assuming an Andrade rheology. The procedure of this coupling for the hemispheric configuration is detailed in Appendix F. We recall that the energy tidally dissipated in the solid part is ignored in the hemispherical configuration. 
In this linear system, I_{N} denotes the identity matrix of size N × N, the forcing terms of the studied tidal potential (Eq. F.3) are expressed as
and the complex characteristic frequencies σ_{n, k} as
where the horizontal wavenumber of the degreen mode , and the coupling coefficients and are defined in Eq. (F.13). Once the solution of this algebraic system is obtained, the selfconsistent tidal response of the Earth is quantified by the total frequency dependent complex Love number defined, for each order, m, and degree, l, as
The first term of the above expression accounts for the direct tidal gravitational forcing of the solid part by the perturber. The second term is related to the oceanic tidal response, which is coupled to that of the solid part through gravitational and surface loading interactions. We remark that the effective Love number characterizing the full tidal response of the planet (Eq. K.8) depends on both the latitudinal and longitudinal harmonic degrees, l and m, in contrast with the solid Love number, . This results from the fact that Coriolis forces alter the oceanic tidal response, which is not the case for the solid tidal response. The contribution of the component of the tidal potential to the total tidal torque exerted on the Earth scales as the imaginary part of the associated Love number and is expressed as (Efroimsky & Williams 2009; Correia et al. 2014)
Since we restricted our analysis to the study of the dominant semidiurnal tide, we only consider the quadrupolar potential with l = m = 2.
In Fig. K.1, we compute the tidal torque for both the hemispheric and global oceanic geometries for a fixed value of H and different orders of magnitude of σ_{R}. We consider the semidiurnal lunar gravitational forcing exerted on the Earth. The spectrum of the torque is plotted against the normalized frequency ω = (Ω − n_{orb})/Ω_{0}, where Ω_{0} is the present spin rate of the Earth. The distribution of resonances associated with surfacegravity modes distorted by rotation is clearly visible for both geometries. In the global ocean case, these resonances are each characterised by the pair of complex frequencies given by (AuclairDesrotour et al. 2018)
which depicts explicitly the predominance of friction over the rotational distortion of tidal waves in a strong friction regime. This can be verified by visual inspection of Fig. K.1. The spectral coverage of the nonresonant background of the torque increases with increasing σ_{R}. In the opposite limit, resonant peaks are spelled out intensifying in amplitude as friction is weakened. Besides, when σ_{R} → 0, the frequencies become real and positive, and we recover the eigenfrequencies of highwavelength surface gravity modes travelling around the sphere. It can be also clearly seen from the high friction regime that the torque of the global ocean is twice that of the hemispherical one, consistent with the simple argument of dissipation increasing proportionally with oceanic area. The same can be deduced if we consider the nonresonant background of the weak friction regime. Comparing the two spectra in this limit reveals the highly irregular nature of the waveforms in the hemispheric response against the fairly regular resonance periodicity in the global configuration. Several resonances can be encountered in the hemispherical configuration spectrum in between two resonant peaks of the global configuration. Applied to the EarthMoon system evolution studied in the main text, we start at the present with a hemispheric ocean, then we switch to a global one. The fitted parameters of H and σ_{R} would place the present torque around a resonant peak, then multiple resonances are crossed before settling into the nonresonant background of the hemispherical ocean response. The switch between the configurations occurs right before surfing the next major resonance.
In Fig. K.2, we plot the torque of the hemispheric configuration for two scenarios: accounting for the oceanic response only and accounting for both the oceanic and solid responses selfconsistently. As explained in Appendix F, the effects of selfattraction and loading interactions between the solid mantle and the oceanic shell are evident in attenuating the amplitude of the response and slightly shifting the position of resonances. This delay effect is due to the influence of this coupling on the phase of resonance depths of nearresonant free oscillations (Müller 2008b).
All Tables
Values of constant parameters used in the numerical implementation of the theory.
All Figures
Fig. 1. Temporal evolution of the latitude of the surface “paleobarycenter” over the last one billion years. The plate tectonics reconstruction is adopted from Merdith et al. (2021), which establishes the first kinematically continuous tectonic motion model across multiple supercontinental cycles. The evolution is smoothed in red using a moving polynomial regression filter with a window of 200 Myr. In our effective model, this curve maps the evolution of the center of the hemispheric continental cap that transitions from being symmetric about the equator during the Mesozoic, to being almost polar during the Paleozoic. 

In the text 
Fig. 2. Misfit surfaces of χ^{2} for the three studied geometric models. The past dynamical evolution of the EarthMoon system is reconstructed for the shown ranges of our two free model parameters H and σ_{R}. The misfit is established using the currently measured lunar recession rate via LLR, and the lunar age (Appendix C). The three models differ in the imposed geometry of the oceanic shell over the geological history, with the combined model featuring more physical realism that the other two. The numerical results of this analysis are summarized in (Table C.1). The dynamical evolution associated with each of the misfit minima is plotted in terms of: the lunar semimajor axis in Fig. 3, length of the day in Fig. 5, and obliquity and precession frequency in Fig. 6. 

In the text 
Fig. 3. Evolution of the lunar semimajor axis over time. The EarthMoon separation, a_{M}, is plotted for the three studied models, taking the bestfit values of the free parameters (H, σ_{R}) as described in Fig. 2 and in the main text. Plotted on top of the evolution curves: Geological inferences of a_{M} from cyclostratigraphy and tidal laminae data (Tables D.1 and D.2). The shaded envelope corresponds to 2σuncertainty in the fitted parameters of the combined model (Appendix C). In the narrow window, we zoom over the most recent 250 Myr of the evolution and make a comparison with the evolution corresponding to explicit numerical tidal modeling using paleogeographic reconstructions (Green et al. 2017) and the prediction of the numerical solution La2004 (Laskar et al. 2004). We note that the integration of a_{E} extends to 3R_{E}, but the yaxis is trimmed to start at 15R_{E} for a better visualization of the geological data. 

In the text 
Fig. 4. History of the tidal torque. The logarithm of the semidiurnal tidal torque of the Earth (normalized by its present value: ) is plotted as a function of time. The solid curve corresponds to the torque of the combined model that involves three phases: in the first phase, a hemispherical ocean migrates on the surface of the Earth following the evolution of the continental barycenter of Fig. 1. Given we lack a continuous plate tectonics model beyond 1 Ga, in Phase 2, we fix the hemispherical ocean to its configuration at 1 Ga to avoid discontinuities in the modeling. It is noteworthy that the attenuated tidal torque over this phase is not due to the fixed oceanic position but due to the tidal response occupying the nonresonant background of the spectrum for the tidal frequencies associated with this interval. Beyond t_{switch}, we enter Phase 3 of the model with the global ocean configuration. The dashed and dasheddotted curves correspond, respectively, to the global and hemispherical oceanic torques that are ignored over the specified intervals by the selective combined model. 

In the text 
Fig. 5. Evolution of the Earth’s length of the day with time. Similar to Fig. 3, but here for the LOD evolution associated with the three studied oceanic models. Geological data on the LOD are summarized in Tables D.1 and D.2. The minimal value reached for the LOD when the integration is terminated at 3 Earth radii is 5.25 h. 

In the text 
Fig. 6. Evolution of the Earth’s obliquity, precession frequency, and precession period with time. The evolution of a_{M} (Fig. 3) and LOD (Fig. 5) are used to compute the evolution of obliquity and precession by Eqs. (A.6) and (A.7). The geological data of the precession frequency from tidal rhythmites and cyclostratigraphy are also plotted on top of the curve (Tables D.1 and D.2). We note that the precession frequency is the directly measured observable in cyclostratigraphy. 

In the text 
Fig. B.1. Drifting effect of the continental cap on the oceanic response: the tidal torque of a hemispheric ocean is plotted as a function of the forcing semidiurnal frequency for different positions of the center of the ocean. With longitudinal symmetry, the latter is defined by the latitude of the oceanic center, which evolves according to Figure 1. The drifting effect on the resonances ranges from position shifting and attenuation for small forcing frequencies to major distortion in the spectrum at larger frequencies. Extreme distortion occurs in the polar oceanic scenario: the major resonance around 11 rad/day reaches a maximum relative to other configurations and the rest of the resonances are absorbed into the background leaving a unimodal spectrum. This behavior makes it important to take into account the position of the hemispherical cap into the model (Figure 1). 

In the text 
Fig. E.1. Eigenfunctions ϕ_{r} (left) and ψ_{r} (right), and the associated tidal flows. The eigenfunctions defined by Eqs. (E.18) and (E.19) are plotted over the hemispherical oceanic domain for 0 ≤ n ≤ 4 (from top to bottom) and 0 ≤ m ≤ n (from left to right). Bright or dark colors designate positive or negative values of the eigenfunctions, respectively. Streamlines indicate the tidal flows corresponding to ∇ϕ_{r} for the set {ϕ_{r}} and to for the set {ψ_{r}}. 

In the text 
Fig. G.1. Adopted transformation scheme that allows recovering the tidal response of a hemispheric ocean with an arbitrary center on the sphere. We use an Eulerian transformation of the form ℛ_{3}(α)ℛ_{2}(β)ℛ_{3}(γ) with γ = 0, allowing us to shift the latitude of the oceanic center O by shifting the spin axis from ŝ to ŝ^{″} in a true polar wander scenario (Webb 1982). 

In the text 
Fig. J.1. Numerical analysis on the dependence of the tidal response computation on the truncation order r_{max}. The response is quantified by the root mean square tidal amplitude, ζ_{rms}, (Eq. J.5) and the dissipative work, , (Eq. J.2), and plotted for three tidal frequencies: 7.3, 11.4, and 22 rad/day that correspond to the vicinity of a tidal resonance, the peak of a resonance, and the background spectrum respectively. 

In the text 
Fig. K.1. Tidal torque between the Earth and the Moon corresponding to the coupled oceanicsolid response in two configurations: a global oceanic shell of thickness H = 4000 m (shown in red), and a hemispherical ocean with the same thickness, symmetric around the equator and bounded by longitudes λ = 0 and λ = π (in blue). Energy dissipation is quantified by the linear Rayleigh drag frequency σ_{R}. The logarithm of the torque is plotted as a function of the normalized frequency ω = (Ω − n_{orb})/Ω_{0}, where the Earth’s spin rate varies with the tidal forcing frequency Ω = n_{orb} + σ/2 at fixed n_{orb}, and Ω_{0} being the present spin rate of the Earth. 

In the text 
Fig. K.2. Similar to Fig.K.1, but comparing the torque of the hemispherical ocean model between a pure oceanic response, and the response of the ocean when accounting for loading and selfattraction effects arising from solid Earth deformation assuming an Andrade rheology. The procedure of this coupling for the hemispheric configuration is detailed in Appendix F. We recall that the energy tidally dissipated in the solid part is ignored in the hemispherical configuration. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.