# The resonant tidal evolution of the Earth-Moon distance

Mohammad Farhat, Pierre Auclair-Desrotour, Gwenaël Boué, Jacques Laskar

► To cite this version:

Mohammad Farhat, Pierre Auclair-Desrotour, Gwenaël Boué, Jacques Laskar. The resonant tidal evolution of the Earth-Moon distance. Astronomy and Astrophysics - A&A, 2022, 665, 10.1051/0004-6361/202243445. insu-03850267

HAL Id: insu-03850267

https://insu.hal.science/insu-03850267

Submitted on 13 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER TO THE EDITOR

# The resonant tidal evolution of the Earth-Moon distance

Mohammad Farhat, Pierre Auclair-Desrotour, Gwenaël Boué, and Jacques Laskar

IMCCE, CNRS, Observatoire de Paris, PSL University, Sorbonne Université, 77 Avenue Denfert-Rochereau, 75014 Paris, France
e-mail: mohammad.farhat@obspm.fr

Received 1 March 2022 / Accepted 13 August 2022

ABSTRACT

Due to tidal interactions in the Earth-Moon 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 Earth-Moon 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 time-varying continental configuration on Earth by considering hemispherical and global ocean models. The resulting evolution of the Earth-Moon 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

1. Introduction

Due to the tidal interplay in the Earth-Moon system, the spin of the Earth slows down over time and the Earth-Moon distance increases (Darwin 1879) at a current rate of 3.830 ± 0.008 cm yr⁻¹, measured using Lunar Laser Ranging (LLR; Williams & Boggs 2016). A rich narrative has been cultivated, exploring the long-term 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 Earth-Moon 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 Earth-Moon system.

Major progress has been achieved with the elaboration of oceanic tidal models. These models offer a tidal frequency-dependent dissipation behavior (Longuet-Higgins 1968; Platzman 1984; Müller 2008a) that allows for the encounter of high-dissipation resonant states over the course of Earth history (Webb 1980; Auclair-Desrotour 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). Paleo-dissipation 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 deep-time 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 time-varying tidal dissipation in the oceans. We propose a physical model that reconciles the two aforementioned limits, described
<figure><img src="image_2.png" /></figure>

**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 super-continental 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 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 Earth-Moon 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 long-term 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 solid-Earth to luni-solar semi-diurnal 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 long-wavelength 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 under-sampling 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). Super-continental formations and breakups thus have their mark on the predicted lunar recession rate.

The dynamical evolution of the Earth-Moon 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 Auclair-Desrotour et al. 2019). The oceanic dissipation is parameterized by an effective frequency, $\sigma_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, $\sigma_R$, can also be interpreted as the inverse of a dissipation timescale, $\tau$, that quantifies the time needed to deplete the kinetic energy budget of tidal oscillations after switching off the forcing. Although $\sigma_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 $\sigma_R$ can completely reshape the tidal spectrum, amplifying the resonant peaks by several orders of magnitude when $\sigma_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; Castillo-Rogez 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 Earth-Moon system that results from the luni-solar semi-diurnal tidal torque for ranges of values of our effective parameters ($H, \sigma_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 Earth-Moon separation, we compute the chi-squared $\chi^2$, taking only two data points into account: the well-constrained lunar age of 4.425 ± 0.025 Ga
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Misfit surfaces of $\chi^2$ for the three studied geometric models. The past dynamical evolution of the Earth-Moon system is reconstructed for the shown ranges of our two free model parameters $H$ and $\sigma_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 than 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 semi-major axis in Fig. 3, length of the day in Fig. 5, and obliquity and precession frequency in Fig. 6.</figcaption></figure>

(Maurice et al. 2020) and the currently measured rate of lunar recession of $3.830 \pm 0.008 \text{ cm yr}^{-1}$ (Williams & Boggs 2016). Misfit surfaces of $\chi^2$ for the three models are shown in Fig. 2. Two $\chi^2$ local minima exist for the global oceanic response; however, one of them corresponds to an unreasonably large average oceanic depth $H \approx 5500 \text{ m}$, leaving us with a global minimum of ($H, \log_{10} \sigma_R$) = (2273 m, -4.89), where $\sigma_R$ is in $\text{s}^{-1}$. The global minimum in the “average” hemispherical ocean model corresponds to ($H, \log_{10} \sigma_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} \sigma_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_{\text{switch}}$, which is implicitly determined by the dynamical integrator (Appendix B). For the best-fit solution, we have $t_{\text{switch}} = 3.25 \text{ 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 best-fit value of $H$ for the combined model corresponds to a volume of $1.19 \times 10^{18} \text{ m}^3$, which is only 10% off from the currently estimated value of $1.33 \times 10^{18} \text{ m}^3$ based on global relief models (Amante & Eakins 2009). The fitted dissipation frequency, $\sigma_R$, corresponds to a decay time $\tau = 43.1 \text{ 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} \sigma_R \in [-4.93, -5.33]$). The best-fit values for the combined model correspond to a lunar trajectory characterized by a current rate of recession $\dot{a}_0 = 3.829 \text{ cm yr}^{-1}$ and an impact time at 4.431 Ga (Table C.1).

# 4. Earth-Moon separation: A history of surfing resonances

For each of the global minima of the misfit parametric studies, we plot the evolution of the Earth-Moon 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 long-term evolution of the torque is characterized by a non-monotonic increase, characteristic of the shrinking Earth-Moon separation, that is interrupted by multiple crossings of resonances. The distribution of resonances in the hemispherical configuration, $t < t_{\text{switch}}$, is less regular than that in the global configuration, $t > t_{\text{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 (Auclair-Desrotour 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
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Evolution of the lunar semi-major axis over time. The Earth-Moon separation, $a_M$, is plotted for the three studied models, taking the best-fit values of the free parameters ($H, \sigma_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 y-axis is trimmed to start at $15R_E$ for a better visualization of the geological data.

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.

Following these resonances, the torque enters a long non-resonant 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 so-called “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 super-continental cycles, although $t_{switch}$ is implicitly determined by the dynamical inte-

<figure><img src="image_5.png" /></figure>

**Fig. 4.** History of the tidal torque. The logarithm of the semi-diurnal tidal torque of the Earth (normalized by its present value: $\bar{T} = T/T(t=0)$) 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 non-resonant 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 dashed-dotted curves correspond, respectively, to the global and hemispherical oceanic torques that are ignored over the specified intervals by the selective combined model.

grator (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 semi-analytical physical model that fits the most accurate constraints in the Earth-Moon 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 Earth-Moon 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
<figure><img src="image_6.png" /></figure>

**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.

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 Earth-Moon 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 Earth-Moon 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.

**Acknowledgements.** We thank Maëlis Arnould for her help with the plate tec-
tonics 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
ANR-19-CE31-0002-01) and by the European Research Council (ERC) under
the European Union's Horizon 2020 research and innovation program (Advanced
Grant AstroGeo-885250). This work was granted access to the HPC resources of
MesoPSL financed by the Region Île-de-France and the project Equip@Meso
(reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir
supervised by the Agence Nationale pour la Recherche. Contributions: MF con-
ducted 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 con-
tributed to the study at all stages. All contributed to the writing of the paper.

<figure><img src="image_7.png" /></figure>

**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.

References

Abramowitz, M., Stegun, I. A., & Romer, R. H. 1988, Handbook of
Mathematical Functions with Formulas, Graphs, and Mathematical Tables

Alcott, L. J., Mills, B. J., & Poulton, S. W. 2019, Science, 366, 1333

Amante, C., & Eakins, B. W. 2009, NOAA Technical Memorandum NESDIS NGDC-24

Anderson, D. L., & Minster, J. B. 1979, Geophys. J. Int., 58, 431

Andrade, E. N. D. C. 1910, Proc. R. Soc. London. Ser. A, 84, 1

Arbic, B. K., & Garrett, C. 2010, Cont. Shelf Res., 30, 564

Arbic, B. K., Karsten, R. H., & Garrett, C. 2009, Atmos. Ocean, 47, 239

Arfken, G. B., & Weber, H. J. 1999, Mathematical Methods for Physicists

Auclair-Desrotour, P., Le Poncin-Lafitte, C., & Mathis, S. 2014, A&A, 561, L7

Auclair-Desrotour, P., Mathis, S., Laskar, J., & Leconte, J. 2018, A&A, 615, A23

Auclair-Desrotour, P., Leconte, J., Bolmont, E., & Mathis, S. 2019, A&A, 629, A132

Bagheri, A., Khan, A., Al-Attar, D., et al. 2019, J. Geophys. Res. Planets, 124, 2703

Bell, T., Jr 1975, J. Geophys. Res., 80, 320

Bolmont, E., Breton, S. N., Tobie, G., et al. 2020, A&A, 644, A165

Boué, G., & Laskar, J. 2006, Icarus, 185, 312

Boulila, S., Laskar, J., Haq, B. U., Galbrun, B., & Hara, N. 2018, Global Planet. Change, 165, 128

Boyden, J. A., Müller, R. D., & Gurnis, M. 2011, Next-generation Plate-tectonic Reconstructions using GPlates (Cambridge University Press)

Carter, G. S., Merrifield, M., Becker, J. M., et al. 2008, J. Phys. Oceanogr., 38, 2205

Castelnau, O., Duval, P., Montagnat, M., & Brenner, R. 2008, J. Geophys. Res. (Solid Earth), 113, B11203

Castillo-Rogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. (Planets), 116, E09008

Correia, A. C. M., & Laskar, J. 2010, in Exoplanets (Tucson, AZ: University of Arizona Press), 239

Correia, A. C., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50

Crease, J. 1966, Tables of the Integral (National Institute of Oceanography)

Čuk, M., Hamilton, D. P., Lock, S. J., & Stewart, S. T. 2016, Nature, 539, 402
Ćuk, M., Hamilton, D. P., & Stewart, S. T. 2019, J. Geophys. Res. Planets, 124, 2917

Daher, H., Arbic, B. K., Williams, J. G., et al. 2021, J. Geophys. Res. Planets, 126, e2021JE006875

Darwin, G. H. 1879, Philos. Trans. R. Soc. London Ser. I, 170, 447

de Azarevich, V. L. L., & Azarevich, M. B. 2017, Geo-Marine Lett., 37, 333

Dhuime, B., Hawkesworth, C. J., Cawood, P. A., & Storey, C. D. 2012, Science, 335, 1334

Dong, S.-H., & Lemus, R. 2002, Appl. Math. Lett., 15, 541

Dziewonski, A. M., & Anderson, D. L. 1981, Phys. Earth Planet. Inter., 25, 297

Efroimsky, M. 2012, ApJ, 746, 150

Efroimsky, M., & Williams, J. G. 2009, Celestial Mech. Dyn. Astron., 104, 257

Eriksson, K. A., & Simpson, E. L. 2000, Geology, 28, 831

Fang, J., Wu, H., Fang, Q., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 540, 109530

Farhat, M. A., & Touma, J. R. 2021, MNRAS, 507, 6078

Farhat, M., Laskar, J., & Boué, G. 2022, J. Geophys. Res. Solid Earth, 127, e2021JB023

Farrell, W. 1972, Rev. Geophys., 10, 761

Fienga, A., Deram, P., Di Ruscio, A., et al. 2021, Notes Scientifiques et Techniques de l'Institut de Mécanique Céleste, 110

Findley, W. N., Lai, J. S., Onaran, K., & Christensen, R. M. 1977, J. Appl. Mech., 44, 364

Fox-Kemper, B., Ferrari, R., & Pedlosky, J. 2003, J. Phys. Oceanogr., 33, 478

Garrett, C., & Munk, W. 1971, Deep Sea Res. Oceanogr. Abstracts, 18, 493

Gent, P. R., & McWilliams, J. C. 1983, Dyn. Atmos. Oceans, 7, 67

Gerkema, T., & Zimmerman, J. 2008, in Lecture Notes (Texel: Royal NIOZ), 207

Gerstenkorn, H. 1967, Icarus, 7, 160

Goldreich, P. 1966, Rev. Geophys., 4, 411

Green, J., Huber, M., Waltham, D., Buzan, J., & Wells, M. 2017, Earth Planet. Sci. Lett., 461, 46

Griffiths, S. D., & Peltier, W. R. 2009, J. Clim., 22, 2905

Gurnis, M., Turner, M., Zahirovic, S., et al. 2012, Comput. Geosci., 38, 35

Han, L., & Huang, R. X. 2020, J. Phys. Oceanogr., 50, 679

Hansen, K. S. 1982, Rev. Geophys., 20, 457

Hawkesworth, C., Cawood, P. A., & Dhuime, B. 2020, Front. Earth Sci., 8

Hough, S. S. 1898, Philos. Trans. R. Soc. London Ser. A, 191, 139

Huang, H., Gao, Y., Jones, M. M., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 550, 109735

Johnson, B. W., & Wing, B. A. 2020, Nat. Geosci., 13, 243

Kaula, W. M. 2013, Theory of Satellite Geodesy: Applications of Satellites to Geodesy (Courier Corporation)

Klatt, J. M., Chennu, A., Arbic, B. K., Biddanda, B., & Dick, G. J. 2021, Nat. Geosci., 14, 564

Lantink, M., Davies, J., & Hilgen, F. 2021, Nat. Rev., 12, 369

Laskar, J. 2005, Celestial Mech. Dyn. Astron., 91, 351

Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261

Lau, H. C., & Faul, U. H. 2019, Earth Planet. Sci. Lett., 508, 18

Lau, H. C., Yang, H.-Y., Tromp, J., et al. 2015, Geophys. J. Int., 202, 1392

Lau, H. C., Faul, U., Mitrovica, J. X., et al. 2016a, Geophys. J. Int., 208, 368

Lau, H. C., Mitrovica, J. X., Austermann, J., et al. 2016b, J. Geophys. Res. Solid Earth, 121, 6991

Lee, U., & Saio, H. 1997, ApJ, 491, 839

Levrard, B., & Laskar, J. 2003, Geophys. J. Int., 154, 970

Longuet-Higgins, M. S. 1968, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci., 262, 511

Longuet-Higgins, M. S., & Pond, G. S. 1970, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci., 266, 193

MacDonald, G. 1967, Proc. R. Soc. London Ser. A. Math. Phys. Sci., 296, 298

Matsuyama, I. 2014, Icarus, 242, 11

Matthews, K. J., Maloney, K. T., Zahirovic, S., et al. 2016, Global Planet. Change, 146, 226

Maurice, M., Tosi, N., Schwinger, S., Breuer, D., & Kleine, T. 2020, Sci. Adv., 6, 28

Mavromatis, H., & Alassar, R. 1999, Appl. Math. Lett., 12, 101

Merdith, A. S., Williams, S. E., Collins, A. S., et al. 2021, Earth-Sci. Rev., 214, 103477

Meyers, S. R., & Malinverno, A. 2018, Proc. Nat. Acad. Sci., 115, 6363

Mignard, F. 1979, Moon Planets, 20, 301

Mojzsis, S. J., Arrhenius, G., McKeegan, K., et al. 1996, Nature, 384, 55

Motoyama, M., Tsunakawa, H., & Takahashi, F. 2020, Icarus, 335, 113382

Müller, M. 2008a, A Large Spectrum of Free Oscillations of the World Ocean Including the Full Ocean Loading and Self-attraction Effects (Springer Science & Business Media), 14

Müller, M. 2008b, Ocean Model., 20, 207

Munk, W. H., & MacDonald, G. J. 1960, The Rotation of the Earth: A Geophysical Discussion (Cambridge University Press)

Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975

Ogilvie, G. I. 2014, ARA&A, 52, 171

Ooe, M. 1989, J. Phys. Earth, 37, 345

Palmer, T., Shutts, G., & Swinbank, R. 1986, Q. J. R. Meteorol. Soc., 112, 1001

Peck, W. H., Valley, J. W., Wilde, S. A., & Graham, C. M. 2001, Geochim. Cosmochim. Acta, 65, 4215

Petit, G., & Luzum, B. 2010, IERS conventions (2010), Tech. rep. (France: Bureau International des Poids et mesures sevres)

Platzman, G. W. 1983, Science, 220, 602

Platzman, G. W. 1984, J. Phys. Oceanogr., 14, 1532

Proudman, J. 1920a, Proc. London Math. Soc., 2, 1

Proudman, J. 1920b, Proc. London Math. Soc., 2, 51

Regge, T. 1958, Il Nuovo Cimento (1955-1965), 10, 544

Renaud, J. P., & Henning, W. G. 2018, ApJ, 857, 98

Riley, K. F., Hobson, M. P., & Bence, S. J. 1999, Mathematical Methods for Physics and Engineering

Ross, M., & Schubert, G. 1989, J. Geophys. Res. Solid Earth, 94, 9533

Rufu, R., & Canup, R. M. 2020, J. Geophys. Res. Planets, 125, 8

Sonett, C., & Chan, M. A. 1998, Geophys. Res. Lett., 25, 539

Sørensen, A. L., Nielsen, A. T., Thibault, N., et al. 2020, Earth Planet. Sci. Lett., 548, 116475

Strauss, W. A. 2007, in Partial Differential Equations: An introduction (John Wiley & Sons)

Sun, C., Xu, W., Cawood, P. A., et al. 2019, Sci. Rep., 9, 1

Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534

Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A&A, 630, A70

Touma, J., & Wisdom, J. 1994, AJ, 108, 1943

Touma, J. 1998, AJ, 115, 1653

Touma, J., & Wisdom, J. 2001, AJ, 122, 1030

Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706

Tyler, R. 2011, Icarus, 211, 770

Tyler, R. H. 2021, Planet. Sci. J., 2, 70

Vallis, G. K. 2017, Atmospheric and Oceanic Fluid Dynamics (Cambridge University Press)

Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1988, Quantum Theory of Angular Momentum (World Scientific)

Walker, J. C. G., & Zahnle, K. J. 1986, Nature, 320, 600

Waltham, D. 2015, J. Sediment. Res., 85, 990

Wang, H., Boyd, J. P., & Akmaev, R. A. 2016, Geosci. Model Dev., 9, 1477

Watterson, I. G. 2001, J. Atmos. Oceanic Technol., 18, 691

Webb, D. 1973, in Deep Sea Research and Oceanographic Abstracts (Elsevier), 20, 847

Webb, D. 1980, Geophys. J. Int., 61, 573

Webb, D. 1982, Geophys. J. Int., 70, 261

Wilde, S. A., Valley, J. W., Peck, W. H., & Graham, C. M. 2001, Nature, 409, 175

Williams, G. E. 1990, J. Phys. Earth, 38, 475

Williams, G. E. 1997, Geophys. Res. Lett., 24, 421

Williams, G. E. 2000, Rev. Geophys., 38, 37

Williams, J. G., & Boggs, D. H. 2016, Celestial Mech. Dyn. Astron., 126, 89

Wood, R., Liu, A. G., Bowyer, F., et al. 2019, Nat. Ecol. Evol., 3, 528

Zahel, W. 1980, Phys. Earth Planet. Inter., 21, 202

Zahnle, K. J., Lupu, R., Dobrovolskis, A., & Sleep, N. H. 2015, Earth Planetary Sci. Lett., 427, 74

Zhong, Y., Wu, H., Fan, J., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 540, 109520
## Appendix A: Orbital dynamics

For the reconstruction of the Earth-Moon 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 inter-vening 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 core-mantle 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 \le 0.1$) (Daher et al. 2021). By attaining a highly eccentric lunar orbit, it is possible through ejection resonance<fnref n="1" />. Touma & Wisdom (1998) studied this regime and showed that capture into such a resonance could have been encountered for $a_M \approx 4.6R_E$, exciting the lunar eccentricity to $e \approx 0.5$. However, the timescale of capture and escape from this ejection resonance is $10^4 \sim 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 long-term 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 long-term 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 core-mantle coupling (Neron de Surgy & Laskar 1997; Touma & Wisdom 2001) are ignored, as well as halts of tidal interaction due to Laplace plane transitions (Cuk et al. 2016).

Under these assumptions, the governing dynamical system of equations is expressed as:

$$ \frac{dL_{\Omega}}{dt} = -(\mathcal{T}_{M} + \mathcal{T}_{S}), \qquad (A.1) $$

$$ \frac{dL_{M}}{dt} = \mathcal{T}_{M}, \qquad (A.2) $$

where $\mathcal{T}_M$ is the lunar semi-diurnal tidal torque coupling between the oceanic and the solid response of the Earth, and $\mathcal{T}_S$ is its solar counterpart. The orbital angular momentum of the Moon $L_M = \beta \sqrt{G(M_E + M_M)a_M}$, where $\beta = M_E M_M / (M_E + M_M)$ is the Earth-Moon system's reduced mass. The rotational angular momentum of the Earth is defined as $L_{\Omega} = C(\Omega)\Omega$, with the time-varying principal moment of inertia given by

(Goldreich 1966)

$$ C(\Omega) = C(\Omega_0) + \frac{2k_2^f R_E^5}{9G}(\Omega^2 - \Omega_0^2). \qquad (A.3) $$

Here, $k_2^f$ is the second-degree fluid Love number of centrifugal/tidal deformation and $G$ is the gravitational constant. The differential equation is integrated backwards in time using the Runge-Kutta 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, \sigma_R$) as input, and the system's variables, $a_M$ and $\Omega$, to compute the tidal frequency and, consequently, the coupled tidal response.

Once the lunar semi-major axis ($a_M$) and the rotation speed of the Earth ($\Omega$) are determined, we compute the obliquity of the Earth ($\epsilon$) 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:

$$ \frac{d\epsilon}{dt} = \frac{Kn}{C(\Omega)\Omega} \sin\epsilon \left(\frac{\Omega}{2n_M} \cos\epsilon - 1\right), \qquad (A.4) $$

and

$$ \frac{da_M}{dt} = \frac{2K}{\beta a_M} \left( \frac{\Omega}{n_M} \cos \epsilon - 1 \right), \qquad (A.5) $$

that is

$$ \frac{d\epsilon}{da_M} = \frac{\beta n_M a_M}{4C(\Omega)\Omega} \sin\epsilon \frac{\Omega \cos\epsilon - 2n_M}{\Omega \cos\epsilon - n_M}. \qquad (A.6) $$

We note that the tidal response parameter, $K$, disappears from the equations. This would also be the case if $K$ depended on $\Omega$. The obliquity evolution equation (A.6) is integrated using the values of $a_M$ and $\Omega$ 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:

$$ p = \frac{3}{2} \left( \frac{GM_S}{a_E^3} + \frac{GM_M}{a_M^3} \right) E_d(\Omega_0) \frac{\Omega}{\Omega_0^2} \cos \epsilon. \qquad (A.7) $$

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 $\Omega_0$, the Earth's semi-major axis $a_E$, and the inertia parameter $C(\Omega_0)/M_E R_E^2$ are adopted from INPOP21 (Fienga et al. 2021). The dynamical ellipticity at the origin of date, $E_d(\Omega_0) = 0.003243$, is determined from the initial conditions for the obliquity ($\epsilon_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 > 4$Ga); 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).

<footnote n="1">¹ occurs when the precession period of the perigee of the Moon equals 1 year, the orbital period of the Earth.</footnote>
## Appendix B: Continental drift and oceanic geometry shifting

<figure><img src="image_8.png" /><figcaption>**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 semi-diurnal 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).</figcaption></figure>

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 time-sliced sketches of Figure 1, we used the GPlates open-source reconstruction software (Boyden et al. 2011; Gurnis et al. 2012).

At 1.5 Ga, the integrator starts computing simultaneously the tidal response of a global oceanic geometry, with uniform thickness $H_{\text{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_{\text{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_{\text{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, \sigma_R$) on the two-dimensional grid. The present rate of lunar recession, $\dot{a}_0$, and the impact time, $t_f$, are then extracted for each evolution sample and the mean square weighted deviation, $\chi^2$, (Table 1) is then computed as:

$$ \chi^2 = \frac{1}{2} \left[ \left( \frac{\dot{a}_0 - \dot{a}_0^{\text{LLR}}}{\sigma^{\text{LLR}}} \right)^2 + \left( \frac{t_f - t_f^{\text{geo}}}{\sigma^{\text{geo}}} \right)^2 \right], \quad (C.1) $$

where we use Lunar Laser Ranging (LLR) estimates of lunar orbital recession (Williams & Boggs 2016): $\dot{a}_0^{\text{LLR}} \pm \sigma^{\text{LLR}} = 38.30 \pm 0.08$ mm/year; as well as geochemical estimates of lunar formation time (Maurice et al. 2020): $t_f^{\text{geo}} \pm \sigma^{\text{geo}} = 4.425 \pm 0.025$ 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.

<table><caption>**Table C.1.** Misfit analysis summary.</caption><thead><tr><th>Model</th><th>Global Ocean</th><th>Hemispherical Ocean</th><th>Combined Model</th></tr></thead><tbody><tr><td>σ<sub>R</sub> [s<sup>-1</sup>]</td><td>1.2770 × 10<sup>-5</sup></td><td>2.8860 × 10<sup>-5</sup></td><td>6.4417 × 10<sup>-6</sup></td></tr><tr><td>H [m]</td><td>2273</td><td>3816</td><td>4674</td></tr><tr><td>ā<sub>0</sub>[cm yr<sup>-1</sup>]</td><td>3.833</td><td>3.828</td><td>3.829</td></tr><tr><td>t<sub>f</sub> [Ga]</td><td>4.422</td><td>4.432</td><td>4.431</td></tr><tr><td>χ<sup>2</sup></td><td>0.0775</td><td>0.0705</td><td>0.0345</td></tr></tbody></table>

**Notes.** Best-fit values of the two free parameters, $\sigma_R$ and $H$, for each of the three studied models are shown, along with the corresponding value of $\chi^2$, as well as the resulting lunar recession rate at the present $\dot{a}_0$ and the impact time, $t_f$.

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 $\dot{a}_0$ and $t_f$, the entries of the variance matrix,

$$ \Sigma = \begin{bmatrix} \text{var}(H) & \text{cov}(H, \sigma_R) \\ \text{cov}(H, \sigma_R) & \text{var}(\sigma_R) \end{bmatrix}, \quad (C.2) $$

are given by

$$ \text{var}(H) = \left(\frac{\partial H}{\partial \dot{a}_0}\right)^2 (\sigma^{\text{LLR}})^2 + \left(\frac{\partial H}{\partial t_f}\right)^2 (\sigma^{\text{geo}})^2, \quad (C.3a) $$

$$ \text{var}(\sigma_R) = \left(\frac{\partial \sigma_R}{\partial \dot{a}_0}\right)^2 (\sigma^{\text{LLR}})^2 + \left(\frac{\partial \sigma_R}{\partial t_f}\right)^2 (\sigma^{\text{geo}})^2, \quad (C.3b) $$

$$ \text{cov}(H, \sigma_R) = \frac{\partial H}{\partial \dot{a}_0} \frac{\partial \sigma_R}{\partial \dot{a}_0} (\sigma^{\text{LLR}})^2 + \frac{\partial H}{\partial t_f} \frac{\partial \sigma_R}{\partial t_f} (\sigma^{\text{geo}})^2 . \quad (C.3c) $$

The partial derivatives entering in these formulae are computed numerically from the fit of ($\dot{a}_0 \pm \sigma^{\text{LLR}}$, $t_f$) and ($\dot{a}_0$, $t_f \pm \sigma^{\text{geo}}$). The marginal uncertainties on parameters $H$ and $\sigma_R$ are $\sigma_H = \text{var}(H)^{1/2} = 32.75$ m and $\sigma_R = \text{var}(\sigma_R)^{1/2} = 0.2631 \times 10^{-6}$ s$^{-1}$, respectively. We used the variance matrix $\Sigma$ to evaluate the $2\sigma$-confidence ellipsoid around the best-fit parameters ($H, \sigma_R$). The Earth-Moon distance, $a_M$, and the length of the day LOD
have been integrated for 25 pairs of parameters ($H, \sigma_R$) chosen at the boundary of this $2\sigma$-confidence region. Their envelope represents the $2\sigma$-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.

<table><caption>**Table D.1.** Cyclostratigraphic data.</caption><thead><tr><th>dataset</th><th>Reference</th><th>T [Ga]</th><th>p [arcsec/yr]</th><th>a<sub>M</sub> [R<sub>E</sub>]</th><th>a<sub>M</sub> [km]</th><th>LOD [hr]</th></tr></thead><tbody><tr><td>IC</td><td>Laskar et al. (2004)</td><td>0.000</td><td><strong>50.467718</strong></td><td>60.142611</td><td>383598</td><td>24.00</td></tr><tr><td>Walvis Ridge</td><td>Meyers & Malinverno (2018)</td><td>0.055</td><td><strong>51.28 ± 1.02</strong></td><td>59.94 ± 0.26</td><td>382284 ± 1650</td><td>23.80 ± 0.24</td></tr><tr><td>Lucaogou(a)</td><td>Huang et al. (2020)</td><td>0.268</td><td><strong>57.01 ± 1.37</strong></td><td>58.57 ± 0.31</td><td>373538 ± 2000</td><td>22.56 ± 0.27</td></tr><tr><td>Lucaogou(b)</td><td>Huang et al. (2020)</td><td>0.270</td><td><strong>55.36 ± 0.51</strong></td><td>58.95 ± 0.12</td><td>375967 ± 750</td><td>22.90 ± 0.11</td></tr><tr><td>Yangtze Block</td><td>Zhong et al. (2020)</td><td>0.445</td><td><strong>57.19 ± 0.53</strong></td><td>58.52 ± 0.12</td><td>373277 ± 750</td><td>22.53 ± 0.10</td></tr><tr><td>Alum shale</td><td>Sørensen et al. (2020)</td><td>0.493</td><td><strong>60.11 ± 1.59</strong></td><td>57.88 ± 0.34</td><td>369153 ± 2200</td><td>21.99 ± 0.28</td></tr><tr><td>Luoyixi</td><td>Fang et al. (2020)</td><td>0.500</td><td><strong>61.06 ± 0.94</strong></td><td>57.67 ± 0.20</td><td>367854 ± 1300</td><td>21.82 ± 0.16</td></tr><tr><td>Xiamaling</td><td>Meyers & Malinverno (2018)</td><td>1.400</td><td><strong>85.79 ± 2.72</strong></td><td>53.27 ± 0.41</td><td>339777 ± 2600</td><td>18.70 ± 0.25</td></tr><tr><td>Joffre</td><td>Lantink et al. (2021)</td><td>2.460</td><td><strong>108.91 ± 8.28</strong></td><td>50.24 ± 0.96</td><td>320452 ± 6100</td><td>16.98 ± 0.50</td></tr></tbody></table>

**Notes.** In boldface, we display the direct observables we used: the precession frequency *p* in arcsec/yr. The time of observation, *T*, is in the second column. The semi-major axis of the Moon (*a*<sub>M</sub>) given in Earth radius (*R*<sub>E</sub>) or in km, and the length of the solar day (LOD), in hours, are derived from the observed quantities using the model that is presented in the text. These values may thus differ from the corresponding values published in the referenced publications; IC denotes the initial conditions (Laskar et al. 2004). The two values (a,b) for the Lucaogou data set correspond to different analyses (a): TimeOptMCMC (Meyers & Malinverno 2018); (b) obliquity and precession cycle counting (Huang et al. 2020). Whenever it is specified in the original publication, the uncertainty in *p* is set to $2\sigma$. The uncertainty of the other variables is propagated through the nominal solution of the present study.

<table><caption>**Table D.2.** Tidal rhythmites data.</caption><thead><tr><th>dataset</th><th>Reference</th><th>T [Ga]</th><th>smo/yr</th><th>p [arcsec/yr]</th><th>a<sub>M</sub> [R<sub>E</sub>]</th><th>LOD [hr]</th></tr></thead><tbody><tr><td>IC</td><td>Laskar et al. (2004)</td><td>0.000</td><td>13.4289</td><td>50.467718</td><td>60.142611</td><td>24.00</td></tr><tr><td>Mansfield</td><td>Sonett & Chan (1998)</td><td>0.310</td><td><strong>13.86 ± 0.21</strong></td><td>55.60 ± 2.55</td><td>58.89 ± 0.59</td><td>22.85 ± 0.52</td></tr><tr><td>Elatina</td><td>Sonett & Chan (1998)</td><td>0.620</td><td><strong>14.93 ± 0.01</strong></td><td>69.24 ± 0.13</td><td>56.04 ± 0.03</td><td>20.56 ± 0.02</td></tr><tr><td>Elatina</td><td>Williams (1997, 2000)</td><td>0.620</td><td><strong>14.10 ± 0.10</strong></td><td>58.55 ± 1.24</td><td>58.22 ± 0.28</td><td>22.27 ± 0.23</td></tr><tr><td>Cottonwood</td><td>Sonett & Chan (1998)</td><td>0.900</td><td><strong>15.33 ± 0.60</strong></td><td>74.68 ± 8.32</td><td>55.06 ± 1.44</td><td>19.87 ± 0.99</td></tr><tr><td>Weeli Wolli</td><td>Williams (1990, 2000)</td><td>2.450</td><td><strong>16.70 ± 1.10</strong></td><td>94.71 ± 17.19</td><td>52.01 ± 2.29</td><td>17.95 ± 1.32</td></tr><tr><td>Weeli Wolli</td><td>Walker & Zahnle (1986)</td><td>2.450</td><td><strong>15.50 ± 0.50</strong></td><td>77.04 ± 7.03</td><td>54.66 ± 1.18</td><td>19.59 ± 0.79</td></tr><tr><td>Moodies Group</td><td>Eriksson & Simpson (2000)<br>de Azarevich & Azarevich (2017)</td><td>3.200</td><td></td><td>148.36 ± 18.61</td><td><strong>46.45 ± 1.50</strong></td><td>15.17 ± 0.65</td></tr></tbody></table>

**Notes.** In boldface, we display the observables. In general, the observable is the number of synodic lunar months per year or in an equivalent way, as quoted here, the number of sidereal lunar months per year (col. 3). The values are issued from the referenced publications (col. 1). For the Moodies Group, we could not infer this quantity from the original publication, and the corresponding estimate of the lunar semi-major axis was taken from de Azarevich & Azarevich (2017). The semi-major axis, $a_M$, is obtained through Kepler's law ($n_M^2 a_M^3 = G(M_M + M_E)$). As for the cyclostratigraphic data (Table.D.1), all other quantities (*p*, LOD) are derived from the observed quantities using the model that is presented in the text. These values may thus differ from the corresponding values published in the referenced publications. The uncertainty of the observables are propagated to the derived variables through the nominal solution of the present study. The values at the origin ($T = 0$) are from Laskar et al. (2004). It should be noted that the present value of sidereal lunar months per year and lunar semi-major axis provided here for $T = 0$ differs from some published value because here we consider averaged values, which should be the case for such long-term studies (see Figure 18 from Laskar et al. (2004)).
<table><caption>**Table D.3.** Values of constant parameters used in the numerical implementation of the theory.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Surface gravitational acceleration (g)</td><td>9.81 m s<sup>-2</sup></td></tr><tr><td>Earth radius (R)</td><td>6378.1366 km</td></tr><tr><td>Solar gravitational constant (GM<sub>S</sub>)</td><td>2.959122082853813556248×10<sup>-4</sup> AU<sup>3</sup> day<sup>-2</sup></td></tr><tr><td>Earth-Moon gravitational constant (G[M<sub>E</sub> + M<sub>M</sub>])</td><td>8.997011395221144381906×10<sup>-10</sup> AU<sup>3</sup> day<sup>-2</sup></td></tr><tr><td>Earth to Moon mass ratio (M<sub>E</sub>/M<sub>M</sub>)</td><td>81.30056789872074318737</td></tr><tr><td>Uniform oceanic density (ρ<sub>oc</sub>)</td><td>1022 kg m<sup>-3</sup></td></tr><tr><td>Andrade characteristic time (τ<sub>A</sub>)</td><td>2.19 × 10<sup>4</sup> yr</td></tr><tr><td>Andrade rheological exponent (α<sub>A</sub>)</td><td>0.25</td></tr><tr><td>Average rigidity of the deformable mantle (μ<sub>E</sub>)</td><td>17.3 × 10<sup>10</sup> Pa</td></tr><tr><td>Average viscosity of the deformable mantle (η<sub>E</sub>)</td><td>3.73 × 10<sup>21</sup> Pa s</td></tr><tr><td>Present day mean lunar semi-major axis (a<sub>0</sub>)</td><td>60.142611 R<sub>E</sub></td></tr><tr><td>Present day mean sidereal length of the day (LOD<sub>s</sub>)</td><td>23.934468 hr</td></tr><tr><td>Present day mean obliquity (ε<sub>0</sub>)</td><td>23.2545°</td></tr><tr><td>Present day mean precession frequency (p<sub>0</sub>)</td><td>50.467718 arcsec yr<sup>-1</sup></td></tr><tr><td>Earth's semi-major axis (a<sub>E</sub>)</td><td>1.495978707 × 10<sup>8</sup> km</td></tr><tr><td>Earth's inertia parameter (C<sub>0</sub>/(M<sub>E</sub>R<sup>2</sup>))</td><td>0.3306947357075918999972</td></tr><tr><td>Earth's fluid Love number (k<sub>2</sub><sup>f</sup>)</td><td>0.93</td></tr></tbody></table>

**Notes.** Astronomical values are adopted from INPOP21 (Fienga et al. 2021). Oceanic and rheological parameters are adopted from Auclair-Desrotour et al. (2019) and Gerkema & Zimmerman (2008). The average rigidity is computed from the PREM model (Dziewonski & Anderson 1981), while the average viscosity is computed from mantle viscosity inversions in Lau et al. (2016b). The initial conditions of the orbital integration are the mean elements from the La2004 astronomical solution (Laskar et al. 2004).

## 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 (Longuet-Higgins & 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 (Longuet-Higgins & Pond 1970; Webb 1980, 1982) contains several misprints that we correct here.

In the co-planar problem under study (ignoring the Earth's obliquity and the lunar orbital inclination), we define a frame of reference co-rotating with the Earth, with a spin vector of $\Omega = \hat{\Omega}$, $\Omega$ being the Earth's spin rate and $\hat{\Omega}$ as the unit vector along its figure axis. In this frame, we used the spherical coordinates $(r, \theta, \lambda)$ denoting the radius, the co-latitude, and the longitude respectively, and their corresponding unit vectors $(\hat{r}, \hat{\theta}, \hat{\lambda})$. 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):

$$ \partial_t u + \sigma_R u + f \times u + g \nabla \zeta = g \nabla \zeta_{eq}, \quad (E.1a) $$

$$ \partial_t \zeta + \nabla \cdot (Hu) = 0, \quad (E.1b) $$

where $u = u_{\theta}\hat{\theta} + u_{\lambda}\hat{\lambda}$ is the horizontal velocity field, $g$ is the gravitational acceleration at the surface, $\zeta$ is the oceanic depth variation, $\zeta_{eq}$ is the equilibrium depth variation, $H$ is the uniform oceanic thickness (the first of only two free parameters

in our model), and $\sigma_R$ is the Rayleigh (or linear) drag frequency (Matsuyama 2014; Auclair-Desrotour 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, $\sigma_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 semi-diurnal 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 Brunt-Väisälä frequency, which quantifies the stability of the ocean's stratification against convection (see, e.g., Gerkema & Zimmerman 2008), or the length-scale of topographical patterns at the oceanic floor (Bell 1975; Palmer et al. 1986). In Eq. (E.1), the Coriolis parameter $f$ is given by:

$$ f = 2\Omega \cos \theta \hat{r}, \quad (E.2) $$

the horizontal gradient operator $\nabla$ is defined as:

$$ \nabla = R^{-1} [\hat{\theta}\partial_{\theta} + \hat{\lambda}(\sin \theta)^{-1} \partial_{\lambda}], \quad (E.3) $$

and the horizontal divergence of the velocity field $\nabla \cdot u$ as:

$$ \nabla \cdot u = (R \sin \theta)^{-1} [\partial_{\theta} (\sin \theta u_{\theta}) + \partial_{\lambda} u_{\lambda}], \quad (E.4) $$

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 = \partial_t x$, where $x$ is the horizontal tidal displacement field, we have:
$$[\partial_t^2 + (\sigma_R + f\mathbf{x})\partial_t]\mathbf{x} + g(\nabla\zeta - \nabla_{\zeta_{eq}}) = 0. \qquad (E.5)$$

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:

$$\mathbf{x} = \nabla\Phi + \nabla\Psi \times \hat{\mathbf{r}}, \qquad (E.6)$$

where $\nabla\Phi$ is a curl-free vector field ($\nabla \times (\nabla\Phi) = 0$) and $\nabla\Psi \times \hat{\mathbf{r}}$ is a divergence-free vector field ($\nabla \cdot (\nabla\Psi \times \hat{\mathbf{r}}) = 0$). In the above equation, we introduce the divergent displacement potential $\Phi$ and the rotational displacement streamfunction $\Psi$ (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 Fox-Kemper 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 $\mathbf{x} \cdot \hat{\mathbf{n}} = 0$, where $\hat{\mathbf{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:

$$\hat{\mathbf{n}} \cdot \nabla\Phi = 0, \qquad \hat{\mathbf{n}} \cdot (\nabla\Psi \times \hat{\mathbf{r}}) = 0. \qquad (E.7)$$

We note that the second condition of the above equation can be rewritten as $(\hat{\mathbf{r}} \times \hat{\mathbf{n}}) \cdot \nabla\Psi = 0$, which implies that $\Psi$ is constant along the coastline (in the following, we set $\Psi = 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, $\Phi$ and $\Psi$. Moreover, the second condition given by Eq. (E.7) enforces the orthogonality of the curl-free and divergence-free components of the tidal flow. By combining together the identity $\nabla \cdot (\Phi\nabla\Psi \times \hat{\mathbf{r}}) = (\nabla\Psi \times \hat{\mathbf{r}}) \cdot \nabla\Phi$ and Gauss' theorem (e.g., Arfken & Weber 1999):

$$\int_O \nabla \cdot (\Phi \nabla\Psi \times \hat{\mathbf{r}}) dA = \oint_{\partial O} \Phi (\nabla\Psi \times \hat{\mathbf{r}}) \cdot \hat{\mathbf{n}} d\ell, \qquad (E.8)$$

with $dA$ and $d\ell$ being the infinitesimal area element of the hemispherical oceanic domain, $O$, and length element of the coastline, $\partial O$, respectively, we obtain

$$\int_O (\nabla\Phi) \cdot (\nabla\Psi \times \hat{\mathbf{r}}) dA = \oint_{\partial O} \Phi (\nabla\Psi \times \hat{\mathbf{r}}) \cdot \hat{\mathbf{n}} d\ell. \qquad (E.9)$$

As the second condition of Eq. (E.7) enforces $(\nabla\Psi \times \hat{\mathbf{r}}) \cdot \hat{\mathbf{n}} = 0$ along the coastline, it follows that

$$\int_O (\nabla\Phi) \cdot (\nabla\Psi \times \hat{\mathbf{r}}) dA = 0, \qquad (E.10)$$

meaning that the components $\nabla\Phi$ and $\nabla\Psi \times \hat{\mathbf{r}}$ 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 divergence-free component of the tidal flow.

The functions $\Phi$ and $\Psi$ are expanded in terms of complete sets of eigenfunctions over the domain, $O$, such that

$$\Phi(\theta, \lambda, t) = \sum_{r=1}^{\infty} p_r(t)\phi_r(\theta, \lambda), \qquad (E.11)$$

$$\Psi(\theta, \lambda, t) = \sum_{r=1}^{\infty} p_{-r}(t)\psi_r(\theta, \lambda). \qquad (E.12)$$

The eigenfunctions $(\phi_r, \psi_r)$ satisfy, over the oceanic domain ($O$), the Helmholtz equations (e.g., Riley et al. (1999), Chapter 21):

$$\nabla^2 \phi_r + \mu_r \phi_r = 0, \qquad (E.13)$$

$$\nabla^2 \psi_r + v_r \psi_r = 0, \qquad (E.14)$$

and, along the coastline ($\partial O$), the boundary conditions given by Eq. (E.7):

$$\hat{\mathbf{n}} \cdot \nabla \phi_r = 0, \qquad \psi_r = 0, \qquad (E.15)$$

where we introduced the horizontal Laplacian:

$$\nabla^2 = (R \sin \theta)^{-2} [\sin \theta \partial_{\theta} (\sin \theta \partial_{\theta}) + \partial_{\lambda\lambda}], \qquad (E.16)$$

and the real eigenvalues, $\mu_r$ and $v_r$, associated with the eigenfunctions $\phi_r$ and $\psi_r$, respectively. We note that the eigenfunctions are normalized such that

$$\int_O \phi_r \phi_s dA = \int_O \psi_r \psi_s dA = \delta_{rs}, \qquad (E.17)$$

where the notation $\delta_{rs}$ referring to the Kronecker $\delta$-symbol is $\delta_{rs} = 1$ for $r = s$ and 0 otherwise. Using these conditions, the eigenfunctions are defined as:

$$\phi_r = \frac{\alpha_{n,m}}{R} P_n^m(\cos \theta) \cos m\lambda, \qquad (E.18)$$

$$\psi_r = \frac{\alpha_{n,m}}{R} P_n^m(\cos \theta) \sin m\lambda, \qquad (E.19)$$

with eigenvalues of $\mu_r = v_r = n(n+1)/R^2$ and the normalization coefficient of:

$$\alpha_{n,m} = \sqrt{\frac{2n+1}{\pi} \frac{(n-m)!}{(n+m)!} \frac{1}{1+\delta_{m0}}}. \qquad (E.20)$$

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 $\phi_r$ Eq. (E.18), $n \in \mathbb{N}$ and $m = 0, 1, ..., n$ while in the expression of $\psi_r$ Eq. (E.19), $n \in \mathbb{N}^*$ and $m = 1, 2, ..., n$. By convention, we set $\psi_0 = 0$ hereafter. Figure E.1 shows the eigenfunctions $\phi_r$ and $\psi_r$ for $1 \le m \le n \le 4$ and the streamlines of the associated tidal flows.

The eigenfunctions $(\phi_r, \psi_r)$ can be split into two sets describing tidal solutions that are symmetric or anti-symmetric 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 (Longuet-Higgins
<figure><img src="image_9.png" /><figcaption>Fig. E.1. Eigenfunctions $\phi_r$ (left) and $\psi_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 \le n \le 4$ (from top to bottom) and $0 \le m \le 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 $\nabla\phi_r$ for the set $\{\phi_r\}$ and to $\nabla\psi_r \times \hat{r}$ for the set $\{\psi_r\}$.</figcaption></figure>

1968; Longuet-Higgins & Pond 1970) for the pairs $(n, m)$. However, in our model, where the ocean is no longer symmetric about the equator, both symmetric and anti-symmetric 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

$$ \zeta = H \sum_{r=1}^{\infty} \mu_r p_r \phi_r. \qquad (E.21) $$

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 $\nabla\phi_r$ and $\nabla\psi_r \times \hat{r}$, then integrating over the oceanic area. Starting with the former, we get

$$ \sum_{s=0}^{\infty} \left( \begin{aligned} & (\partial_t^2 p_s + \sigma_R \partial_t p_s + gH p_s \mu_s - g\zeta_{eq,s}) \nabla\phi_s \cdot \nabla\phi_r \\ & + (\partial_t^2 p_{-s} + \sigma_R \partial_t p_{-s}) (\nabla\psi_s \times \hat{r}) \cdot \nabla\phi_r \\ & + \partial_t p_s (f \times \nabla\phi_s) \cdot \nabla\phi_r + \partial_t p_{-s} [f \times (\nabla\psi_s \times \hat{r})] \cdot \nabla\phi_r \end{aligned} \right) = 0. \qquad (E.22) $$

The product of the gradients of two eigenfunctions is computed using Green's first identity (e.g., Strauss (2007), Chapter 7):

$$ \int_O \nabla\phi_s \cdot \nabla\phi_r dA = \int_{\partial O} \phi_s (\nabla\phi_r \cdot \hat{n}) d\ell - \int_O \phi_s \nabla^2 \phi_r dA. \qquad (E.23) $$

The first term on the right-hand 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:

$$ \int_O \nabla\phi_s \cdot \nabla\phi_r dA = \mu_r \delta_{rs}. \qquad (E.24) $$

Rearranging the other terms using vector identities, we can write Eq. (E.22) as:

$$ \begin{aligned} & \sum_{s=0}^{\infty} \left( \begin{aligned} & (\partial_t^2 p_s + \sigma_R \partial_t p_s + gH p_s \mu_s - g\zeta_{eq,s}) \mu_r \delta_{rs} \\ & + (\partial_t^2 p_{-s} + \sigma_R \partial_t p_{-s}) \int_O (\nabla\psi_s \times \hat{r}) \cdot \nabla\phi_r dA \\ & + \partial_t p_s \int_O f \cdot (\nabla\phi_s \times \nabla\phi_r) dA \\ & + \partial_t p_{-s} \int_O (f \cdot \hat{r}) (\nabla\phi_r \cdot \nabla\psi_s) dA \end{aligned} \right) = 0. \end{aligned} \qquad (E.25) $$

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

$$ \begin{aligned} & (\partial_t^2 p_r + \sigma_R \partial_t p_r + gH p_r \mu_r - g\zeta_{eq,r}) \mu_r \\ & - 2\Omega \sum_{s=1}^{\infty} \partial_t p_s \int \cos \theta \hat{r} \cdot (\nabla\phi_r \times \nabla\phi_s) dA \\ & + 2\Omega \sum_{s=1}^{\infty} \partial_t p_{-s} \int \cos \theta (\nabla\phi_r \cdot \nabla\psi_s) dA = 0. \end{aligned} \qquad (E.26) $$
To close the system, we multiply the momentum equation with $\nabla\psi_r \times \hat{r}$ to get

$$ \sum_{s=0}^{\infty} \left( \partial_t^2 p_s + \sigma_R \partial_t p_s + gH p_s \mu_s - g\zeta_{eq,s} \right) \nabla\phi_s \cdot (\nabla\psi_r \times \hat{r}) \\ + \partial_t p_s (\mathbf{f} \times \nabla\phi_s) \cdot (\nabla\psi_r \times \hat{r}) + \partial_t p_{-s} [\mathbf{f} \times (\nabla\psi_s \times \hat{r})] \cdot (\nabla\psi_r \times \hat{r}) \\ + (\partial_t^2 p_{-s} + \sigma_R \partial_t p_{-s}) (\nabla\psi_s \times \hat{r}) \cdot (\nabla\psi_r \times \hat{r}) = 0. \quad (E.27) $$

Integrating Eq. (E.27) over the area of the ocean and using basic vector product identities, we get

$$ \sum_{s=0}^{\infty} \left( \partial_t^2 p_s + \sigma_R \partial_t p_s \right) \int \nabla \psi_s \cdot \nabla \psi_r dA \\ - \partial_t p_s \int (\mathbf{f} \cdot \hat{r}) (\nabla \phi_s \cdot \nabla \psi_r) dA \\ - \partial_t p_{-s} \int (\mathbf{f} \cdot \hat{r}) \hat{r} \cdot (\nabla \psi_r \times \nabla \psi_s) dA = 0, \quad (E.28) $$

where upon replacing the Coriolis term as before we have

$$ \partial_t^2 p_{-r} + \sigma_R \partial_t p_{-r} - \frac{2\Omega}{\nu_r} \sum_{s=1}^{\infty} \partial_t p_s \int \cos \theta (\nabla \phi_s \cdot \nabla \psi_r) dA \\ - \frac{2\Omega}{\nu_r} \sum_{s=0}^{\infty} \partial_t p_{-s} \int \cos \theta \hat{r} \cdot (\nabla \psi_r \times \nabla \psi_s) dA = 0. \quad (E.29) $$

We identify in Eqs. (E.26) and (E.29) the so-called “gyroscopic coefficients” (e.g., Proudman (1920a,b)) that are defined as

$$ \begin{aligned} \beta_{r,s} &= -\int_O \cos \theta \hat{r} \cdot (\nabla \phi_r \times \nabla \phi_s) dA, \\ \beta_{r,-s} &= \int_O \cos \theta \nabla \phi_r \cdot \nabla \psi_s dA, \\ \beta_{-r,s} &= -\int_O \cos \theta \nabla \psi_r \cdot \nabla \phi_s dA, \\ \beta_{-r,-s} &= -\int_O \cos \theta \hat{r} \cdot (\nabla \psi_r \times \nabla \psi_s) dA, \end{aligned} \quad (E.30) $$

with $\beta_{-s,r} = -\beta_{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:

$$ \partial_t^2 p_r + \sigma_R \partial_t p_r + gH \mu_r p_r - g\zeta_{eq,r} + \frac{2\Omega}{\mu_r} \sum_{s=-\infty}^{s=\infty} \beta_{r,s} \partial_t p_s = 0, \quad (E.31) $$

$$ \partial_t^2 p_{-r} + \sigma_R \partial_t p_{-r} + \frac{2\Omega}{\nu_r} \sum_{s=-\infty}^{s=\infty} \beta_{-r,s} \partial_t p_s = 0. \quad (E.32) $$

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 self-attraction and loading between the ocean and the deforming mantle in order to have a complete self-consistent 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 self-attraction between the oceanic shell and the solid part (Farrell 1972; Zahel 1980). If we were to take these into account when studying oceanic tides, $\zeta$ becomes a function of two moving surfaces: the free oceanic surface, $\zeta_{os}$, and the vertically deforming solid surfacem $\zeta_{ss}$.

In the frame co-rotating with the Earth as defined in Appendix E, the gravitational potential is expressed as (e.g., Auclair-Desrotour et al. 2019)

$$ U(\mathbf{r}, \mathbf{r}') = \frac{GM}{|\mathbf{r} - \mathbf{r}'|} - \frac{GM}{r'^2} r \cos \theta, \quad (F.1) $$

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:

$$ U^T(\theta, \lambda, \mathbf{r}') = U(R, \theta, \lambda, \mathbf{r}') - \frac{GM}{r'}, \quad (F.2) $$

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 $U_n^{m,s}$, as (Kaula 2013; Auclair-Desrotour et al. 2019):

$$ U^T = \sum_{n=2}^{\infty} \sum_{m=-n}^{n} \sum_{s=-\infty}^{\infty} U_n^{m,s} P_n^m(\cos \theta) \exp\{i(\sigma_m^s t + m\lambda)\}, \quad (F.3) $$

where $s$ is an integer and the tidal forcing frequency $\sigma_m^s = m\Omega - sn_{orb}$, 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 $U_n^{m,s}$ is given by (Ogilvie 2014):

$$ U_n^{m,s} = \frac{GM}{a} \left(\frac{R}{a}\right)^n A_{n,m,s}(e), \quad (F.4) $$

where $a$ is the semi-major 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 semi-diurnal component identified by $n = m = s = 2$ and corresponding to the tidal frequency $\sigma_2^2 = 2(\Omega - n_{orb})$. For this component, while neglecting the small orbital eccentricity of the Sun and the Moon, $A_{2,2,2}(0) = \sqrt{3}/5$. Hereafter, we use $U_n^T$ to represent a single harmonic ($n, m, s$) of the tidal potential. This harmonic of degree $n$ is defined as

$$ U_n^T = U_n^{m,s} P_n^m(\cos \theta) \exp\{i(\sigma_m^s t + m\lambda)\}. \quad (F.5) $$

Moreover, in the following, we write $\sigma$ instead of $\sigma_m^s$ to simplify the notation. Subject to $U_n^T$ only, the equilibrium oceanic depth would be $\bar{\zeta} = U_n^T/g$. However, the loading effect of the deforming oceanic shell adds to the tidal potential and they both affect the ocean surface $\bar{\zeta}_{os}$ and the ocean floor corresponding to the solid surface $\bar{\zeta}_{ss}$. The former is expressed as (Matsuyama 2014):

$$ \bar{\zeta}_{os} = \frac{h_n^T U_n^T}{g} + \sum_l \frac{3\rho_{oc}}{(2l+1)\rho_{se}} h_l^L \zeta_l, \quad (F.6) $$
where $\rho_{oc}$ and $\rho_{se}$ stand for the uniform oceanic and solid Earth densities, respectively. In this equation, the oceanic depth variation, $\zeta$, is decomposed into spherical harmonics defined over the full sphere as:

$$ \zeta_l(\theta, \lambda, t) = \sum_{m=-l}^{l} \zeta_l^m(t) P_l^m(\cos \theta) \exp(im\lambda). \qquad (F.7) $$

Although $\zeta$ given in (E.21) 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 $\zeta(\theta, \lambda, t) = 0$ over the continental hemisphere, Eq. (F.7) can also be expressed as:

$$ \zeta_l(\theta, \lambda, t) = \frac{1}{2} \sum_{m=0}^{l} \alpha_{lm}^2 \int_O \zeta(\theta', \lambda', t) P_l^m(\cos \theta) P_l^m(\cos \theta') \\ \times \cos(m(\lambda - \lambda')) d\Omega, \qquad (F.8) $$

where the integral is computed over the solid angle $2\pi$ spanned by the ocean. The second contribution to the equilibrium tide, which is due to the solid redistribution of mass, is

$$ \bar{\zeta}_{ss} = (1 + k_n^T) \frac{U_n^T}{g} + \sum_l \frac{3\rho_{oc}}{(2l+1)\rho_{se}} (1 + k_l^L) \zeta_l. \qquad (F.9) $$

In equations (F.6) and (F.9), we used the tidal Love numbers $k_n^T$ and $h_n^T$, and the surface-loading Love numbers $k_n^L$ and $h_n^L$, 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 visco-elastic 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 ($k_n^T, h_n^T, k_n^L, h_n^L$) 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 closed-form solutions derived for a uniform solid interior (Munk & MacDonald 1960):

$$ \{k_n^T, h_n^T, k_n^L, h_n^L\} = \frac{1}{(1 + \tilde{\mu}_n)} \left\{ \frac{3}{2(2n-1)}, \frac{2n+1}{2(2n-1)}, -1, -\frac{2n+1}{3} \right\}, \qquad (F.10) $$

where $\tilde{\mu}_n$ 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 $\tilde{\mu}_n$, 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 (Castillo-Rogez et al. 2011; Auclair-Desrotour 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, $\tilde{\mu}_n$ takes the following form (Findley et al. 1977; Efroimsky 2012):

$$ \tilde{\mu}_n = \frac{4(2n^2 + 4n + 3)\pi R^4}{3nGM_E^2} \frac{\mu_E}{1 + (i\sigma\tau_A)^{-\alpha_A}\Gamma(1 + \alpha_A) + (i\sigma\tau_M)^{-1}}, \qquad (F.11) $$

where $M_E$ is the mass of the Earth, $\mu_E$ its average rigidity, and $\Gamma$ is the gamma function (Abramowitz et al. 1988); $\alpha_A$ is a dimensionless rheological exponent determined experimentally (Castelnau et al. 2008; Petit & Luzum 2010); $\tau_A$ is the anelastic Andrade timescale; and $\tau_M$ the Maxwell relaxation time defined as the ratio of viscosity to rigidity. For a volumetric average of the mantle's shear modulus $\mu_E = 17.3 \times 10^{10}$ Pa, and volumetric average of viscosity deduced from inversions of Lau et al. (2016b), we have $\tau_M = 685$ yrs. The values $\alpha_A = 0.25$ and $\tau_A = 2.19 \times 10^4$ yrs that we use in our model are adopted from Auclair-Desrotour 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 $\zeta_{eq}$ in the momentum equation (E.1a) by the difference $\bar{\zeta}_{os} - \bar{\zeta}_{ss}$ of Eqs. (F.6) and (F.9) and we resolve it in the Fourier domain using the forcing tidal frequency, $\sigma$. The modified momentum conservation equations is now expressed as:

$$ i\sigma u + \sigma_R u + f \times u = -g \nabla \left( -\gamma_n^T \bar{\zeta} + \sum_l \gamma_l^L \zeta_l \right), \qquad (F.12) $$

with $\bar{\zeta} = U_n^T/g$, and where the loading and tidal tilt factors are defined as (Matsuyama 2014):

$$ \gamma_n^T = 1 + k_n^T - h_n^T ; \quad \gamma_l^T = 1 - \frac{3\rho_{oc}}{(2l+1)\rho_{se}}(1 + k_l^L - h_l^L). \qquad (F.13) $$

Just like the Love numbers, $\gamma_n^T$ and $\gamma_l^L$ 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 ocean-solid coupling to the linear system of $p_r$ and $p_{-r}$. Multiplying the added contribution of loading and self-attraction effects by $\nabla\phi_r$ and $\nabla\psi_r \times \hat{r}$, then resolving the added terms in the frequency domain, after a number of manipulations, we can finally re-write the expression of the system of Eqs. (E.31) and (E.32) as:

$$ -\sigma^2 p_r - i\sigma\sigma_R p_r + gH\mu_r \left(1 - \frac{\gamma_r^L}{2}\right) p_r - g\gamma_n^T \bar{\zeta}_r - \frac{2i\sigma\Omega}{\mu_r} \sum_{s=-\infty}^{s=\infty} \beta_{r,s} p_s - \frac{1}{2} gH \sum_{\substack{s'=1 \\ s' \neq r}}^{\infty} \mu_{s'} F_r^{s'} p_{s'} = 0, \qquad (F.14) $$

$$ -\sigma^2 p_{-r} - i\sigma\sigma_R p_{-r} - \frac{2i\sigma\Omega}{\nu_r} \sum_{s=-\infty}^{s=\infty} \beta_{-r,s} \partial_t p_s = 0, \qquad (F.15) $$

where we define

$$ F_r^{s'} = 4\alpha_{n,m}\alpha_{n',m'} \sum_p \sum_q \gamma_p^L q^2 \alpha_{p,q}^{s'} \frac{\mathcal{O}_{p,q}^{n,m} \mathcal{O}_{p,q}^{n',m'}}{(q^2 - m^2)(q^2 - m'^2)}, \qquad (F.16) $$

with $\mathcal{O}_{p,q}^{n,m}$ corresponding to the Gram matrix of the ALFs,

$$ \mathcal{O}_{n,m}^{u,v} = \int_{-1}^{1} P_{n}^{m}(\mu) P_{u}^{w}(\mu) d\mu, \qquad (F.17) $$

for which the method of computation is detailed in Appendix H. Coupled to the orbital dynamical evolution, the tidal frequency $\sigma$ 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 re-write the linear system as:

$$ (a^{(1)} + a_r^{(2)})p_r + a_r^{(3)} \sum_{s=-\infty}^{s=\infty} \beta_{r,s} p_s + a^{(5)} \sum_{\substack{s'=1 \\ s' \neq r}}^{\infty} \mu_{s'} F_r^{s'} p_{s'} = c_r, \quad (\text{F.18}) $$

$$ a^{(1)} p_r + a_r^{(4)} \sum_{s=-\infty}^{s=\infty} \beta_{r,s} p_s = 0, \quad (\text{F.19}) $$

where the first equation is for $r > 0$ and the second is for $r < 0$ – and where we have introduced the coefficients:

$$ \begin{aligned} a^{(1)} &= -\sigma^2 - i\sigma\sigma_R, & a_r^{(2)} &= gH\mu_r(1 - \gamma_r^L/2), \\ a_r^{(3)} &= -2i\sigma\Omega\mu_r^{-1}, & a_r^{(4)} &= -2i\sigma\Omega\nu_r^{-1}, \\ a^{(5)} &= -\frac{1}{2}gH, & c_r &= g\gamma_n^T \bar{\zeta}_r. \end{aligned} \quad (\text{F.20}) $$

## 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)):

$$ P_n^m(\mu) = \frac{(-1)^m}{2^n n!} (1 - \mu^2)^{m/2} \partial_\mu^{n+m} (\mu^2 - 1)^n, \quad (\text{G.1}) $$

which are solutions to the Legendre equation,

$$ \partial_\mu [(1 - \mu^2)\partial_\mu P_n^m] + \left[ n(n + 1) - \frac{m^2}{1 - \mu^2} \right] P_n^m = 0. \quad (\text{G.2}) $$

Following differentiation, we obtain

$$ \partial_\mu P_n^m = - \frac{m\mu}{1 - \mu^2} P_n^m - \frac{P_n^{m+1}}{\sqrt{1 - \mu^2}}. \quad (\text{G.3}) $$

Substituting Eq. (G.3) in Eq. (G.2) we get the recurrence relation

$$ P_n^{m+2} - \frac{2m\mu^2(m+1)}{1-\mu^2} P_n^m - 2\mu(m+1)\partial_\mu P_n^m + (n(n+1)-m(m+1))P_n^m = 0, \quad (\text{G.4}) $$

which gives the useful relation

$$ \mu \partial_\mu P_n^m = \frac{P_n^{m+2}}{2(m+1)} + \left[ \frac{n(n+1)+m(m+1)}{2(m+1)} - \frac{m}{1-\mu^2} \right] P_n^m. \quad (\text{G.5}) $$

From Eqs. (G.1-G.5), it is straightforward to obtain the ALFs recurrence relations that are necessary to compute the integral equations of the gyroscopic coefficients,

$$ \frac{\mu P_n^m}{\sqrt{1 - \mu^2}} = -\frac{1}{2m} \left( P_n^{m+1} + (n - m + 1)(n + m)P_n^{m-1} \right), \quad (\text{G.6}) $$

$$ \frac{P_n^m}{\sqrt{1 - \mu^2}} = -\frac{1}{2m} \left( P_{n-1}^{m+1} + (n + m - 1)(n + m)P_{n-1}^{m-1} \right), \quad (\text{G.7}) $$

$$ \sqrt{1 - \mu^2} \partial_\mu P_n^m = -\frac{1}{2} P_n^{m+1} + \frac{1}{2} (n + m)(n - m + 1) P_n^{m-1}, \quad (\text{G.8}) $$

$$ (1 - \mu^2) \partial_\mu P_n^m = \frac{1}{2n + 1} \left( (n + 1)(n + m)P_{n-1}^m - n(n - m + 1)P_{n+1}^m \right). \quad (\text{G.9}) $$

<figure><img src="image_10.png" /><figcaption>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 $R_3(\alpha)R_2(\beta)R_3(\gamma)$ with $\gamma = 0$, allowing us to shift the latitude of the oceanic center $O$ by shifting the spin axis from $\hat{s}$ to $\hat{s}''$ in a true polar wander scenario (Webb 1982).</figcaption></figure>

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 $\{\theta, \lambda\} \rightarrow \{\theta', \lambda'\}$ using an Eulerian rotation matrix of the form $R_3(\alpha)R_2(\beta)R_3(\gamma)$, with $(0 \le \alpha \le 2\pi)$ and $(0 \le \beta \le \pi)$, and we fix $\gamma = 0$ (see Fig. G.1). For a vector $J$ defined as:

$$ \begin{aligned} J &= R_3(\alpha)R_2(\beta) \begin{pmatrix} 0 & 0 & 1 \end{pmatrix}^T \\ &= (\sin\beta \cos\alpha \quad \sin\beta \sin\alpha \quad \cos\beta)^T, \end{aligned} \quad (\text{G.10}) $$

the transformed gyroscopic coefficients are:

$$ \frac{R^2}{\alpha_r \alpha_s} \beta_{r,s} = J_z \beta_{r,s}^{(1)} + J_x \beta_{r,s}^{(2)} + J_y \beta_{r,s}^{(3)}, \quad (\text{G.11}) $$

$$ \frac{R^2}{\alpha_r \alpha_s} \beta_{r,-s} = J_z \beta_{r,-s}^{(1)} + J_x \beta_{r,-s}^{(2)} + J_y \beta_{r,-s}^{(3)}, \quad (\text{G.12}) $$

$$ \frac{R^2}{\alpha_r \alpha_s} \beta_{-r,-s} = J_z \beta_{-r,-s}^{(1)} + J_x \beta_{-r,-s}^{(2)} + J_y \beta_{-r,-s}^{(3)}, \quad (\text{G.13}) $$

$$ \beta_{-r,s} = -\beta_{s,-r}, \quad (\text{G.14}) $$

where, for $r$ associated with the harmonic pair of integers $(n, m)$ and $s$ associated with $(u, v)$, we introduce the
coefficients:

$$ \beta_{r,s}^{(1)} = \frac{2}{m^2 - v^2} \int [v^2 P_u^v \partial_\mu P_n^m + m^2 P_n^m \partial_\mu P_u^v] \mu d\mu, \quad (G.15) $$

$$ \beta_{r,s}^{(2)} = \frac{\pi}{4} \int [m K_{m,v}^{(1)} P_n^m \partial_\mu P_u^v \bar{\mu}^{1/2} - v K_{m,v}^{(2)} \partial_\mu P_n^m P_u^v \bar{\mu}^{1/2}] d\mu, \quad (G.16) $$

$$ \beta_{r,s}^{(3)} = \int [m K_{m,v}^{(3)} P_n^m \partial_\mu P_u^v \bar{\mu}^{1/2} - v K_{m,v}^{(4)} \partial_\mu P_n^m P_u^v \bar{\mu}^{1/2}] d\mu, \quad (G.17) $$

$$ \beta_{r,-s}^{(1)} = \frac{-2v}{m^2 - v^2} \int [\partial_\mu P_u^v \partial_\mu P_n^m \mu \bar{\mu} + m^2 P_n^m P_u^v \bar{\mu}] d\mu, \quad (G.18) $$

$$ \beta_{r,-s}^{(2)} = \frac{\pi}{4} \int [K_{m,v}^{(2)} \partial_\mu P_n^m \partial_\mu P_u^v \bar{\mu}^{3/2} - m v K_{m,v}^{(1)} P_n^m P_u^v \bar{\mu}^{-1/2}] d\mu, \quad (G.19) $$

$$ \beta_{r,-s}^{(3)} = \int [K_{m,v}^{(4)} \partial_\mu P_n^m \partial_\mu P_u^v \bar{\mu}^{3/2} - m v K_{m,v}^{(3)} P_n^m P_u^v \bar{\mu}^{-1/2}] d\mu, \quad (G.20) $$

$$ \beta_{-r,-s}^{(1)} = \frac{2mv}{m^2 - v^2} \int [P_u^v \partial_\mu P_n^m + P_n^m \partial_\mu P_u^v] \mu d\mu, \quad (G.21) $$

$$ \beta_{-r,-s}^{(2)} = \frac{\pi}{4} \int [v K_{m,v}^{(1)} \partial_\mu P_n^m P_u^v \bar{\mu}^{1/2} - m K_{m,v}^{(2)} P_n^m \partial_\mu P_u^v \bar{\mu}^{1/2}] d\mu, \quad (G.22) $$

$$ \beta_{-r,-s}^{(3)} = \int [v K_{m,v}^{(3)} \partial_\mu P_n^m P_u^v \bar{\mu}^{1/2} - m K_{m,v}^{(4)} P_n^m \partial_\mu P_u^v \bar{\mu}^{1/2}] d\mu, \quad (G.23) $$

with $\bar{\mu} = 1 - \mu^2$ and

$$ K_{m,v}^{(1)} = (1 + \delta_{v,0}) \delta_{m-v,1} - \delta_{m-v,-1}, \quad (G.24) $$

$$ K_{m,v}^{(2)} = K_{v,m}^{(1)}, \quad (G.25) $$

$$ K_{m,v}^{(3)} = m \left( \frac{1}{m^2 - (v^2 + 1)^2} + \frac{1}{m^2 - (v^2 - 1)^2} \right), \quad (G.26) $$

$$ K_{m,v}^{(4)} = K_{v,m}^{(3)}, \quad (G.27) $$

Under the terms of this transformation, the latitude of the center of the ocean in the rotating frame is given by:

$$ \cos \theta' = \cos \theta \cos \beta + \sin \theta \sin \beta \cos(\lambda - \alpha). \quad (G.28) $$

To compute the integrals involved in the gyroscopic coefficients, we make use of the essential condition<fnref n="1" /> (e.g., Longuet-Higgins & Pond (1970))

$$ P_n^m P_u^v \Big|_{\mu=\pm 1} = 0, \quad (G.29) $$

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 Longuet-Higgins & 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 $\beta_{r,s}^{(1)}$ and $\beta_{r,s}^{(2)}$ are expressed as:

$$ \beta_{r,s}^{(1)} = \left[ \frac{u(u+1)+v(v+1)}{v+1} - 2v \frac{n(n+1)-u(u+1)+v}{m^2-v^2} \right] O_{n,m}^{u,v} + \frac{1}{v+1} O_{n,m}^{u,v+2}, \quad (G.30) $$

$$ \beta_{r,s}^{(2)} = \frac{\pi}{4} \left[ m K_{m,v}^{(1)} \int P_n^m \partial_\mu P_u^v \bar{\mu}^{1/2} d\mu - v K_{m,v}^{(2)} \int \partial_\mu P_n^m P_u^v \bar{\mu}^{1/2} d\mu \right], \quad (G.31a) $$

$$ = \frac{\pi}{8} \left\{ m K_{m,v}^{(1)} \left[ (u+v)(u-v+1) O_{n,m}^{u,v-1} - O_{n,m}^{u,v+1} \right] - v K_{m,v}^{(2)} \left[ (n+m)(n-m+1) O_{n,m-1}^{u,v} - O_{n,m+1}^{u,v} \right] \right\}, \quad (G.31b) $$

where we used Eq. (G.8) for each integrand in Eq. (G.31a) to obtain Eq. (G.31b). The coefficients $\beta_{r,s}^{(3)}$, $\beta_{r,-s}^{(1)}$, and $\beta_{r,-s}^{(2)}$ read as:

$$ \beta_{r,s}^{(3)} = \frac{1}{2} \left\{ m K_{m,v}^{(3)} \left[ (u+v)(u-v+1) O_{n,m}^{u,v-1} - O_{n,m}^{u,v+1} \right] - v K_{m,v}^{(4)} \left[ (n+m)(n-m+1) O_{n,m-1}^{u,v} - O_{n,m+1}^{u,v} \right] \right\}, \quad (G.32) $$

$$ \beta_{r,-s}^{(1)} = \frac{-2v}{(m^2 - v^2)(2n + 1)} \left\{ \begin{aligned} & (n + 1)(n - 1)(n + m) O_{n-1,m}^{u,v} \\ & + n(n + 2)(n - m + 1) O_{n+1,m}^{u,v} \end{aligned} \right\}, \quad (G.33) $$

$$ \beta_{r,-s}^{(2)} = \frac{\pi}{4} \left[ K_{m,v}^{(2)} \int \partial_\mu P_n^m \partial_\mu P_u^v \bar{\mu}^{3/2} d\mu - m v K_{m,v}^{(1)} \int P_n^m P_u^v \bar{\mu}^{-1/2} d\mu \right], \quad (G.34a) $$

$$ = \frac{\pi}{4} \left[ K_{m,v}^{(2)} \int \partial_\mu P_n^m \bar{\mu} \partial_\mu P_u^v \bar{\mu}^{1/2} d\mu - m v K_{m,v}^{(1)} \int P_n^m P_u^v \bar{\mu}^{-1/2} d\mu \right], \quad (G.34b) $$

$$ = \frac{\pi K_{m,v}^{(2)}}{8(2n+1)} \left\{ \begin{aligned} & (n+1)(n+m) \left[ (u+v)(u-v+1) O_{n-1,m}^{u,v-1} - O_{n-1,m}^{u,v+1} \right] \\ & + n(n-m+1) \left[ O_{n+1,m}^{u,v+1} - (u+v)(u-v+1) O_{n+1,m}^{u,v-1} \right] \end{aligned} \right\} + \frac{\pi v K_{m,v}^{(1)}}{8} \left\{ O_{n-1,m+1}^{u,v} + (n+m-1)(n+m) O_{n-1,m-1}^{u,v} \right\}, \quad (G.34c) $$

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 $\beta_{r,-s}^{(3)}$, $\beta_{-r,-s}^{(1)}$, $\beta_{-r,-s}^{(2)}$, and $\beta_{-r,-s}^{(3)}$ are expressed as:

$$ \beta_{r,-s}^{(3)} = \frac{K_{m,v}^{(4)}}{2(2n+1)} \left\{ \begin{aligned} & (n+1)(n+m) \left[ (u+v)(u-v+1) O_{n-1,m}^{u,v-1} - O_{n-1,m}^{u,v+1} \right] \\ & + n(n-m+1) \left[ O_{n+1,m}^{u,v+1} - (u+v)(u-v+1) O_{n+1,m}^{u,v-1} \right] \end{aligned} \right\} + \frac{v K_{m,v}^{(3)}}{2} \left\{ O_{n-1,m+1}^{u,v} + (n+m-1)(n+m) O_{n-1,m-1}^{u,v} \right\}, \quad (G.35) $$

$$ \beta_{-r,-s}^{(1)} = \frac{-2mv}{m^2 - v^2} O_{n,m}^{u,v}, \quad (G.36) $$

<footnote n="1">¹ we note that this general condition is invalid in the case where $n = m = u = v = 0$. However, this case is excluded here by the definition of the eigenfunctions in Eqs. (E.18) and (E.19).</footnote>
$$
\beta_{-r,-s}^{(2)} = \frac{\pi}{8} \left\{ -m K_{m,v}^{(2)} \left[ (u+v)(u-v+1)O_{n,m}^{u,v-1} - O_{n,m}^{u,v+1} \right] \right. \\
\left. \qquad + v K_{m,v}^{(1)} \left[ (n+m)(n-m+1)O_{n,m-1}^{u,v} - O_{n,m+1}^{u,v} \right] \right\}, \qquad (G.37)
$$

$$
\beta_{-r,-s}^{(3)} = \frac{1}{2} \left\{ -m K_{m,v}^{(4)} \left[ (u+v)(u-v+1)O_{n,m}^{u,v-1} - O_{n,m}^{u,v+1} \right] \right. \\
\left. \qquad + v K_{m,v}^{(3)} \left[ (n+m)(n-m+1)O_{n,m-1}^{u,v} - O_{n,m+1}^{u,v} \right] \right\}. \qquad (G.38)
$$

**Appendix H: Overlap integral** $O_{n,m}^{u,v}$

Here, we provide a closed form solution for the computa-
tion 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:

$$
O_{n,m}^{u,v} = C_{n,m}^{u,v} \sum_l (2l+1) \mathcal{D}(|q|, l) \cdot \begin{pmatrix} n & u & l \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} n & u & l \\ -m & v & m-v \end{pmatrix}, \quad (H.1)
$$

where the factors $C_{n,m}^{u,v}$ are given by:

$$
C_{n,m}^{u,v} = (-1)^\kappa 2^{|q|-2} |q| \sqrt{\frac{(n+m)!(u+v)!}{(n-m)!(u-v)!}}, \quad (H.2)
$$

and the coefficients $\mathcal{D}(|q|, l)$ by

$$
\mathcal{D}(|q|, l) = \left[1 + (-1)^{l+|q|}\right] \sqrt{\frac{(l - |q|)!}{(l + |q|)!} \frac{\Gamma(l/2)\Gamma((l + |q| + 1)/2)}{(l - |q|)/2!\Gamma((l + 3)/2)}. \quad (H.3)
$$

We note here that the phase $\kappa$ introduced in Dong & Lemus (2002) as

$$
\kappa = \begin{cases} m & \text{if } v \ge m, \\ v & \text{otherwise,} \end{cases} \tag{H.4}
$$

corrects the phase given in Mavromatis & Alassar (1999) and
Crease (1966), where the latter was used for the computation of
the gyroscopic coefficients in Longuet-Higgins & 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 3-jm symbols are determined from Varshalovich et al. (1988) by

$$
\begin{pmatrix} a & b & c \\ d & e & f \end{pmatrix} = (-1)^{R_{21}+R_{31}+R_{32}} \left[ \frac{R_{31}!R_{32}!R_{33}!R_{13}!R_{23}!}{(J+1)!R_{11}!R_{12}!R_{21}!R_{22}!} \right]^{1/2} \\
\times \sum_{z} (-1)^z \frac{(R_{21}+z)!(R_{11}+R_{31}-z)!}{z!(R_{31}-z)!(R_{23}-z)!(R_{13}-R_{31}+z)!} \tag{H.5}
$$

where $J = a + b + c$, and $R_{ij}$ are the elements of the so-called
Regge $\Re$-symbol (Regge 1958) defined as

$$
\begin{align*}
R_{11} &= -a + b + c, & R_{12} &= a - b + c, & R_{13} &= a + b - c, \\
R_{21} &= a + d, & R_{22} &= b + e, & R_{23} &= c + f, \\
R_{31} &= a - d, & R_{32} &= b - e, & R_{33} &= c - f. \tag{H.6}
\end{align*}
$$

The summation in Eq. (H.5) runs over all integer values of z for
which all the factorial arguments are non-negative. Finally, we
note that using Eq. (H.1), $O_{n,m}^{u,v} = 0$ when $v = m$. In that case, we
alternatively use

$$
O_{n,m}^{u,m} = \frac{2}{2n+1} \frac{(n+m)!}{(n-m)!} \delta_{n,u}. \qquad (H.7)
$$

This method for the computation of the overlap integral was ver-
ified numerically using MATLAB's ALFs package.

**Appendix I: The tidal forcing term** $\tilde{\zeta}_r$

As in Webb (1980), considering the equilibrium tide $\tilde{\zeta}$ to have a
unit root mean square amplitude, and to be driven by a spherical
harmonics term:

$$
Y_p^q(\theta, \lambda) = \sqrt{\frac{2p+1}{4\pi} \frac{(p-q)!}{(p+q)!}} P_p^q(\cos\theta) \exp(iq\lambda), \quad (I.1)
$$

with angular frequency $\sigma$, we have:

$$
\bar{\zeta} = \sqrt{2\pi} Y_p^q(\theta, \lambda) \exp(i\sigma t). \quad (I.2)
$$

Under the rotation of the coordinate system described by the
Euler angles (α, β, γ) (see Appendix G and Fig. G.1), the spheri-
cal harmonics transform as (Varshalovich et al. 1988):

$$
Y_p^s(\theta', \lambda') = \sum_{q=-p}^{p} Y_p^q(\theta, \lambda) D_{s,q}^p(\alpha, \beta, \gamma), \quad (I.3)
$$

or

$$
Y_p^q(\theta, \lambda) = \sum_{s=-p}^{p} Y_p^s(\theta', \lambda') D_{s,q}^{p*}(\alpha, \beta, \gamma), \quad (I.4)
$$

where $D_{q,s}^p$ designate the Wigner D-functions. These functions
are themselves the product of three functions (Varshalovich et al.
1988), each depending on one argument: $\alpha, \beta$, or $\gamma$,

$$
D_{s,q}^{p}(\alpha, \beta, \gamma) = e^{-iq\alpha} d_{sq}^{p}(\beta)e^{-is\gamma}. \quad (I.5)
$$

In this expression, the functions $d_{sq}^p(\beta)$ are given by

$$
d_{sq}^{p}(\beta) = (-1)^{p-s} \left[ (p+q)!(p-q)!(p+s)!(p-s)! \right]^{1/2} \\
\times \sum_{j} (-1)^{j} \frac{(\cos \beta/2)^{q+s+2j} (\sin \beta/2)^{2p-q-s-2j}}{j!(p-q-j)!(p-s-j)!(q+s+j)!}, \quad (1.6)
$$

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 semi-diurnal tide (*p* = *q* = 2), we are left with one term only. Expanding the harmonic factor $Y_p^s(\theta', \lambda')$ 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,

$$
\bar{\zeta} = \sqrt{\frac{\pi R^2}{2}} \exp(i\sigma t) \sum_{s=-p}^{p} D_{s,q}^{p*}(\alpha, \beta, \gamma) (1 + \delta_{s,0})^{1/2} [\phi_p^s + i\psi_p^s]. \quad (1.7)
$$

Then, invoking the definition of the component $\tilde{\zeta}_r$,

$$
\bar{\zeta}_r = \int_O \phi_r \bar{\zeta} dA, \qquad (1.8)
$$
we get its expression in the rotated frame of reference,

$$ \tilde{\zeta}_r = \sqrt{\frac{\pi R^2}{2}} \exp(i\sigma t) \sum_{s=-p}^{p} D_{s,q}^{p*}(\alpha, \beta, \gamma) (1 + \delta_{s,0})^{1/2} \times \left[ \int_O \phi_p^s \phi_n^m dA + i \int_O \psi_p^s \phi_n^m dA \right], \quad (I.9) $$

where the dot products of the eigenfunctions are simplified as:

$$ \int_O \phi_p^s \phi_n^m dA = \begin{cases} \delta_{n,p}, & \text{if } s = m, \\ (-1)^m \delta_{n,p}, & \text{if } s = -m, \\ 0, & \text{otherwise,} \end{cases} \quad (I.10) $$

and

$$ \int_O \psi_p^s \phi_n^m dA = \begin{cases} 0, & \text{if } m + s = \text{even}, \\ \alpha_{p,s} \alpha_{n,m} \frac{2s}{s^2 - m^2} O_{p,s}^{n,m} & \text{otherwise.} \end{cases} \quad (I.11) $$

We note that as the index s takes negative values, we use

$$ P_p^{-s} = (-1)^s \frac{(p-s)!}{(p+s)!} P_p^s \quad (I.12) $$

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 Earth-Moon 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 $\mathcal{T}$, the power lost by the Earth would be $\mathcal{T}\Omega$, and the power gained by the perturber is $\mathcal{T}n_{\text{orb}}$. The difference between them is the dissipative work of the total tidal mass redistribution $\mathcal{W}_{\text{diss}}$, thus

$$ \mathcal{T} = \frac{\mathcal{W}_{\text{diss}}}{\Omega - n_{\text{orb}}}. \quad (J.1) $$

The total dissipative work is the sum of two contributions: the dissipative work of oceanic tidal currents, $\mathcal{W}_{\text{diss}}^{\text{oc}}$, and dissipation in the deforming viscoelastic mantle. In the formalism established thus far, we calculated the self-consistently 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 solid-oceanic 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 self-gravitational 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, $\mathcal{W}_{\text{diss}}^{\text{oc}}$. 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

<figure><img src="image_11.png" /><figcaption>Fig. J.1. Numerical analysis on the dependence of the tidal response computation on the truncation order $r_{\text{max}}$. The response is quantified by the root mean square tidal amplitude, $\zeta_{\text{rms}}$, (Eq. J.5) and the dissipative work, $\mathcal{W}_{\text{diss}}^{\text{oc}}$, (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.</figcaption></figure>

contributions is only jeopardized by the emerging significance of the solid dissipation when moving backwards in time and increasing the Earth's rotational velocity $\Omega$. Solid Earth dissipation would also be amplified with an early less viscous mantle due to higher Hadean-Archean 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 mid-Archean, beyond which we self-consistently 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.

The oceanic dissipative work is given by (Webb 1980):

$$ \begin{aligned} \mathcal{W}_{\text{diss}}^{\text{oc}} &= \left\langle \int_O \mathbf{u}(t) \cdot \sigma_{\text{R}} \mathbf{u}(t) dA \right\rangle \\ &= \frac{1}{2} \sigma_{\text{R}} \sigma^2 \sum_{r=1}^{\infty} (\mu_r p_r p_r^* + v_r p_{-r} p_{-r}^*), \end{aligned} \quad (J.2) $$

where $\langle \rangle$ denotes time averaging over the tidal period. This work should be equal to the work done by the tidal force on the ocean:

$$ \begin{aligned} \mathcal{W}_{\text{tide}}^{\text{oc}} &= \left\langle \rho_{\text{oc}} g H \int_O \nabla \bar{\zeta}(t) \cdot \mathbf{u}(t) dA \right\rangle \\ &= \frac{1}{2} \rho_{\text{oc}} g H \sigma \operatorname{Im} \left\{ \sum_{r=1}^{\infty} p_r \mu_r \bar{\zeta}_r^* \right\}. \end{aligned} \quad (J.3) $$

Hence, the tidal torque associated with the lunar semi-diurnal frequency $\sigma = 2(\Omega - n_M)$, $n_M$ being the lunar mean motion, is

$$ \mathcal{T}_M = \rho_{\text{oc}} g H \operatorname{Im} \left\{ \sum_{r=1}^{\infty} p_r \mu_r \bar{\zeta}_r^* \right\}, \quad (J.4) $$

and we obtain a similar expression for the solar tides $\mathcal{T}_S$ when solving the system with the solar tidal frequency component $\sigma = 2(\Omega - 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, $\zeta_{rms}$, given as

$$ \zeta_{rms} = \sqrt{\frac{H}{\pi R^2} \sum_{r=1}^{\infty} \mu_r^2 p_r^* p_r}. \qquad (J.5) $$

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 Auclair-Desrotour 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<fnref n="2" />. 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:

$$ \begin{aligned} \boldsymbol{u} &= \sum_{m, \sigma} \boldsymbol{u}^{m, \sigma}(\theta) \exp\{i(\sigma t + m\lambda)\}, \\ \zeta &= \sum_{m, \sigma} \zeta^{m, \sigma}(\theta) \exp\{i(\sigma t + m\lambda)\}, \\ \zeta_{\text{eq}} &= \sum_{m, \sigma} \zeta_{\text{eq}}^{m, \sigma}(\theta) \exp\{i(\sigma t + m\lambda)\}. \end{aligned} \qquad (K.1) $$

Defining the complex tidal frequency $\tilde{\sigma}$ and the complex spin parameter $\tilde{\nu}$ as in Auclair-Desrotour et al. (2018) by:

$$ \tilde{\sigma} = \sigma - i\sigma_R \quad \text{and} \quad \tilde{\nu} = \frac{2\Omega}{\tilde{\sigma}}, \qquad (K.2) $$

and replacing the tidal quantities by their expansions, the governing system reduces to an eigenvalue-eigenfunction 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) {$\Theta_n^{m,\tilde{\nu}}(\theta)$}, associated with a set of eigenvalues {$A_n^{m,\tilde{\nu}}$}. 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

$$ \begin{aligned} \Theta_n^{m,\tilde{\nu}}(\theta) &= \sum_{m \le l} A_{n,l}^{m,\tilde{\nu}} P_l^m(\cos \theta), \\ P_l^m(\cos \theta) &= \sum_n B_{l,n}^{m,\tilde{\nu}} \Theta_n^{m,\tilde{\nu}}(\theta), \end{aligned} \qquad (K.3) $$

with $A_{n,l}^{m,\tilde{\nu}}$ and $B_{l,n}^{m,\tilde{\nu}}$ being complex change of basis coefficients. Using the change of basis coefficients $A_{n,l}^{m,\tilde{\nu}}$, the tidal displacement solution is expressed as:

$$ \zeta_l^{m,\sigma} = \sum_n A_{n,l}^{m,\tilde{\nu}} \zeta_n^{m,\sigma}, \qquad (K.4) $$

where the components $\zeta_n^{m,\sigma}$ are solutions to the linear algebraic system

$$ \left\{ \begin{array}{c} \begin{bmatrix} \sigma_{1,1}^2 & \dots & \sigma_{1,n}^2 & \dots & \sigma_{1,N}^2 \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \sigma_{n,1}^2 & \dots & \sigma_{n,n}^2 & \dots & \sigma_{n,N}^2 \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \sigma_{N,1}^2 & \dots & \sigma_{N,n}^2 & \dots & \sigma_{N,N}^2 \end{bmatrix} \begin{bmatrix} \zeta_1^{m,\sigma} \\ \vdots \\ \zeta_n^{m,\sigma} \\ \vdots \\ \zeta_N^{m,\sigma} \end{bmatrix} = \begin{bmatrix} \mathcal{F}_1^{m,\sigma} \\ \vdots \\ \mathcal{F}_n^{m,\sigma} \\ \vdots \\ \mathcal{F}_N^{m,\sigma} \end{bmatrix} \\ \sigma \tilde{\sigma} \mathbf{I}_N - \sum_{l=1}^{N} \mathcal{F}_l^{m,\sigma} = -\frac{H \Lambda_n^{m,\tilde{\nu}}}{R^2} \sum_{m \le l} B_{l,n}^{m,\tilde{\nu}} \gamma_l^T U_l^{m,\sigma}, \end{array} \right. \qquad (K.5) $$

In this linear system, $\mathbf{I}_N$ denotes the identity matrix of size $N \times N$, the forcing terms of the studied tidal potential $U_l^{m,\sigma}$ (Eq. F.3) are expressed as

$$ \mathcal{F}_n^{m,\sigma} = -\frac{H \Lambda_n^{m,\tilde{\nu}}}{R^2} \sum_{m \le l} B_{l,n}^{m,\tilde{\nu}} \gamma_l^T U_l^{m,\sigma}, \qquad (K.6) $$

and the complex characteristic frequencies $\sigma_{n,k}$ as

$$ \sigma_{n,k} = \sqrt{gH\hat{k}_n^2 \sum_{l \ge m} \gamma_l^T A_{k,l}^{m,\tilde{\nu}} B_{l,n}^{m,\tilde{\nu}}}, \qquad (K.7) $$

where the horizontal wave-number of the degree-$n$ mode $\hat{k}_n = \sqrt{\Lambda_n^{m,\tilde{\nu}}/R}$, and the coupling coefficients $\gamma_l^T$ and $\gamma_l^T$ are defined in Eq. (F.13). Once the solution of this algebraic system is obtained, the self-consistent tidal response of the Earth is quantified by the total frequency dependent complex Love number defined, for each order, $m$, and degree, $l$, as

$$ \hat{k}_l^{m,\sigma} = \hat{k}_l^T + (1 + \hat{k}_l^T) \frac{3g}{2l+1} \frac{\rho_{\text{oc}}}{\rho_{\text{se}}} \frac{\zeta_l^{m,\sigma}}{U_l^{m,\sigma}}. \qquad (K.8) $$

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, $\hat{k}_l^T$. 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 $U_l^{m,\sigma}$ 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)

$$ \mathcal{T}_l^m = \frac{3}{2} GM^2 \frac{R^5}{a^6} \operatorname{Im}\{\hat{k}_l^{m,\sigma}\}. \qquad (K.9) $$

Since we restricted our analysis to the study of the dominant semi-diurnal 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$

<footnote n="2">² We remind the reader that we can proceed with this theory as such only because we are studying dynamics in the coplanar setting. The theory requires further development if one were to account for the Earth's obliquity and lunar inclination.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. K.1.** Tidal torque between the Earth and the Moon corresponding to the coupled oceanic-solid 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 $\lambda = 0$ and $\lambda = \pi$ (in blue). Energy dissipation is quantified by the linear Rayleigh drag frequency $\sigma_R$. The logarithm of the torque is plotted as a function of the normalized frequency $\omega = (\Omega - n_{orb})/\Omega_0$, where the Earth's spin rate varies with the tidal forcing frequency $\Omega = n_{orb} + \sigma/2$ at fixed $n_{orb}$, and $\Omega_0$ being the present spin rate of the Earth.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**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 self-attraction 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.</figcaption></figure>

and different orders of magnitude of $\sigma_R$. We consider the semi-diurnal lunar gravitational forcing exerted on the Earth. The spectrum of the torque is plotted against the normalized frequency $\omega = (\Omega - n_{orb})/\Omega_0$, where $\Omega_0$ is the present spin rate of the Earth. The distribution of resonances associated with surface-gravity 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 (Auclair-Desrotour et al. 2018)

$$ \sigma_n^{\pm} = i \frac{\sigma_R}{2} \pm \sqrt{gHk_n^2 - \left(\frac{\sigma_R}{2}\right)^2}, \qquad (K.10) $$

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 non-resonant background of the torque increases with increasing $\sigma_R$. In the opposite limit, resonant peaks are spelled out intensifying in amplitude as friction is weakened. Besides, when $\sigma_R \to 0$, the frequencies $\sigma_n^{\pm}$ become real and positive, and we recover the eigenfrequencies of high-wavelength 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 proportion-

ally with oceanic area. The same can be deduced if we consider the non-resonant background of the weak friction regime. Comparing the two spectra in this limit reveals the highly irregular nature of the waveforms in the hemispherical 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 Earth-Moon 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 $\sigma_R$ would place the present torque around a resonant peak, then multiple resonances are crossed before settling into the non-resonant 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 self-consistently. As explained in Appendix F, the effects of self-attraction 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 near-resonant free oscillations (Müller 2008b).