EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 665 (September 2022) A&A, 665 (2022) L1 Full HTML
Open Access
Issue
A&A
Volume 665, September 2022
Article Number L1
Number of page(s) 20
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202243445
Published online 05 September 2022

© M. Farhat et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

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−1, 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 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.

thumbnail Fig. 1.

Temporal evolution of the latitude of the surface “paleo-barycenter” 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.

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, σR, which globally models the bottom friction and the conversion of barotropic flows into internal gravity waves, both mechanisms amounting to ∼91% of the total dissipation (Carter et al. 2008). This frequency, σR, can also be interpreted as the inverse of a dissipation timescale, τ, that quantifies the time needed to deplete the kinetic energy budget of tidal oscillations after switching off the forcing. Although σR is probably a function of local topography, its spatial variation can be averaged out longitudinally over the Earth’s rotation and latitudinally over precession and plate tectonics. The second free parameter in our model is the uniform effective oceanic thickness, H. The imprints of these two parameters on the tidal response spectrum are distinguishable: variations in H smoothly shifts the positions of the resonant peaks while slightly varying their amplitudes. In contrast, variations in σR can completely reshape the tidal spectrum, amplifying the resonant peaks by several orders of magnitude when σR decreases or, otherwise, completely absorbing the resonant peaks into the background spectrum (Fig. K.1). For the hemispherical geometry, we adopt the analytical approach of Webb (1980) (see Appendix E), in which the tidal solution is expanded in spherical harmonics (Fig. E.1). In both geometries, we take into account the effect of the deformation of the solid part of the Earth adopting an Andrade rheology (Anderson & Minster 1979; 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, σ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 χ2, taking only two data points into account: the well-constrained lunar age of 4.425 ± 0.025 Ga (Maurice et al. 2020) and the currently measured rate of lunar recession of 3.830 ± 0.008 cm yr−1 (Williams & Boggs 2016). Misfit surfaces of χ2 for the three models are shown in Fig. 2. Two χ2 local minima exist for the global oceanic response; however, one of them corresponds to an unreasonably large average oceanic depth H ≈ 5500 m, leaving us with a global minimum of (H, log10σR)=(2273 m, −4.89), where σR is in s−1. The global minimum in the “average” hemispherical ocean model corresponds to (H, log10σ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, log10σ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 tswitch, which is implicitly determined by the dynamical integrator (Appendix B). For the best-fit solution, we have tswitch = 3.25 Ga, which is in agreement with assumptions from the literature of the existence of a global ocean until ∼2.5 Ga (Dhuime et al. 2012; Hawkesworth et al. 2020; Motoyama et al. 2020). If we assume that the oceanic volume is conserved over time, the best-fit value of H for the combined model corresponds to a volume of 1.19 × 1018 m3, which is only 10% off from the currently estimated value of 1.33 × 1018 m3 based on global relief models (Amante & Eakins 2009). The fitted dissipation frequency, σR, corresponds to a decay time τ = 43.1 h, which is consistent with real oceanic studies (Garrett & Munk 1971; Webb 1973) that offer a range between 24 and 60 h (or log10σR ∈ [ − 4.93, −5.33]). The best-fit values for the combined model correspond to a lunar trajectory characterized by a current rate of recession ̇a0 = 3.829 cm yr−1 and an impact time at 4.431 Ga (Table C.1).

thumbnail Fig. 2.

Misfit surfaces of χ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 σR. The misfit is established using the currently measured lunar recession rate via LLR, and the lunar age (Appendix C). The three models differ in the imposed geometry of the oceanic shell over the geological history, with the combined model featuring more physical realism that the other two. The numerical results of this analysis are summarized in (Table C.1). The dynamical evolution associated with each of the misfit minima is plotted in terms of: the lunar semi-major axis in Fig. 3, length of the day in Fig. 5, and obliquity and precession frequency in Fig. 6.

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 <  tswitch, is less regular than that in the global configuration, t >  tswitch, (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 aM, 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 aM of 2.8RE within 330 Myr. These resonances are associated with rapid variations of the Earth’s obliquity (Fig. 6) that could have triggered major climatic events. We observe that the first resonance overlaps with the Paleozoic oxygenation event (∼350 Ma), while the second overlaps with the Neoproterozoic major oxygenation event (∼600 Ma) and the Cambrian Explosion (Wood et al. 2019). Possible correlation between the Earth’s LOD and the benthic ecosystem should thus be considered (Klatt et al. 2021). The second resonance peak is almost half an order of magnitude lower than in the global configuration. This is an essential feature of the combined model for preserving the lunar angular momentum budget at this stage to better match the cyclostratigraphic proxy estimates at 1.4 and 2.5 Ga, which clearly cannot be explained by the other, more dissipative models considered in Fig. 2.

thumbnail Fig. 3.

Evolution of the lunar semi-major axis over time. The Earth-Moon separation, aM, is plotted for the three studied models, taking the best-fit values of the free parameters (H, σR) as described in Fig. 2 and in the main text. Plotted on top of the evolution curves: Geological inferences of aM 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 aE extends to 3RE, but the y-axis is trimmed to start at 15RE for a better visualization of the geological data.

thumbnail Fig. 4.

History of the tidal torque. The logarithm of the semi-diurnal tidal torque of the Earth (normalized by its present value: T ∼ = T / T ( t = 0 ) $ \tilde{\mathcal{T}}= \mathcal{T}/\mathcal{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 tswitch, 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.

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

thumbnail Fig. 6.

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

Following these resonances, the torque enters a long 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 tswitch is implicitly determined by the dynamical integrator (Appendix B). Samples of continental growth curves predict a fast decay in continental crust volume beyond tswitch (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 aM of 6.5RE 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 sometimes to divergent interpretations (Walker & Zahnle 1986; Sonett & Chan 1998; Williams 2000; Table D.2).

We obtained a striking fit with the estimate of aM 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.

1

occurs when the precession period of the perigee of the Moon equals 1 year, the orbital period of the Earth.

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

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.

Acknowledgments

We thank Maëlis Arnould for her help with the plate tectonics model and the GPlates software and Matthias Sinnesael for discussions on the geological data. We are grateful to Margriet Lantink and coworkers for the communication of their results on the Joffre sequence before publication and allowing us to include their data point in the present work. This project has been supported by the French Agence Nationale de la Recherche (AstroMeso 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 conducted the simulations and drafted the paper. PAD brought expertise in oceanic tides and GB in solid tides. JL initiated the study and supervised it. All contributed to the study at all stages. All contributed to the writing of the paper.

References

  1. Abramowitz, M., Stegun, I. A., & Romer, R. H. 1988, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables [Google Scholar]
  2. Alcott, L. J., Mills, B. J., & Poulton, S. W. 2019, Science, 366, 1333 [NASA ADS] [CrossRef] [Google Scholar]
  3. Amante, C., & Eakins, B. W. 2009, NOAA Technical Memorandum NESDIS NGDC-24 [Google Scholar]
  4. Anderson, D. L., & Minster, J. B. 1979, Geophys. J. Int., 58, 431 [NASA ADS] [CrossRef] [Google Scholar]
  5. Andrade, E. N. D. C. 1910, Proc. R. Soc. London. Ser. A, 84, 1 [CrossRef] [Google Scholar]
  6. Arbic, B. K., & Garrett, C. 2010, Cont. Shelf Res., 30, 564 [NASA ADS] [CrossRef] [Google Scholar]
  7. Arbic, B. K., Karsten, R. H., & Garrett, C. 2009, Atmos. Ocean, 47, 239 [CrossRef] [Google Scholar]
  8. Arfken, G. B., & Weber, H. J. 1999, Mathematical Methods for Physicists [Google Scholar]
  9. Auclair-Desrotour, P., Le Poncin-Lafitte, C., & Mathis, S. 2014, A&A, 561, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Auclair-Desrotour, P., Mathis, S., Laskar, J., & Leconte, J. 2018, A&A, 615, A23 [EDP Sciences] [Google Scholar]
  11. Auclair-Desrotour, P., Leconte, J., Bolmont, E., & Mathis, S. 2019, A&A, 629, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bagheri, A., Khan, A., Al-Attar, D., et al. 2019, J. Geophys. Res. Planets, 124, 2703 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bell, T., Jr 1975, J. Geophys. Res., 80, 320 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bolmont, E., Breton, S. N., Tobie, G., et al. 2020, A&A, 644, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Boué, G., & Laskar, J. 2006, Icarus, 185, 312 [Google Scholar]
  16. Boulila, S., Laskar, J., Haq, B. U., Galbrun, B., & Hara, N. 2018, Global Planet. Change, 165, 128 [CrossRef] [Google Scholar]
  17. Boyden, J. A., Müller, R. D., & Gurnis, M. 2011, Next-generation Plate-tectonic Reconstructions using GPlates (Cambridge University Press) [Google Scholar]
  18. Carter, G. S., Merrifield, M., Becker, J. M., et al. 2008, J. Phys. Oceanogr., 38, 2205 [NASA ADS] [CrossRef] [Google Scholar]
  19. Castelnau, O., Duval, P., Montagnat, M., & Brenner, R. 2008, J. Geophys. Res. (Solid Earth), 113, B11203 [NASA ADS] [CrossRef] [Google Scholar]
  20. Castillo-Rogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. (Planets), 116, E09008 [CrossRef] [Google Scholar]
  21. Correia, A. C. M., & Laskar, J. 2010, in Exoplanets (Tucson, AZ: University of Arizona Press), 239 [Google Scholar]
  22. Correia, A. C., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Crease, J. 1966, Tables of the Integral (National Institute of Oceanography) [Google Scholar]
  24. Ćuk, M., Hamilton, D. P., Lock, S. J., & Stewart, S. T. 2016, Nature, 539, 402 [CrossRef] [Google Scholar]
  25. Ćuk, M., Hamilton, D. P., & Stewart, S. T. 2019, J. Geophys. Res. Planets, 124, 2917 [Google Scholar]
  26. Daher, H., Arbic, B. K., Williams, J. G., et al. 2021, J. Geophys. Res. Planets, 126, e2021JE006875 [CrossRef] [Google Scholar]
  27. Darwin, G. H. 1879, Philos. Trans. R. Soc. London Ser. I, 170, 447 [NASA ADS] [Google Scholar]
  28. de Azarevich, V. L. L., & Azarevich, M. B. 2017, Geo-Marine Lett., 37, 333 [NASA ADS] [CrossRef] [Google Scholar]
  29. Dhuime, B., Hawkesworth, C. J., Cawood, P. A., & Storey, C. D. 2012, Science, 335, 1334 [NASA ADS] [CrossRef] [Google Scholar]
  30. Dong, S.-H., & Lemus, R. 2002, Appl. Math. Lett., 15, 541 [CrossRef] [Google Scholar]
  31. Dziewonski, A. M., & Anderson, D. L. 1981, Phys. Earth Planet. Inter., 25, 297 [CrossRef] [Google Scholar]
  32. Efroimsky, M. 2012, ApJ, 746, 150 [Google Scholar]
  33. Efroimsky, M., & Williams, J. G. 2009, Celestial Mech. Dyn. Astron., 104, 257 [NASA ADS] [CrossRef] [Google Scholar]
  34. Eriksson, K. A., & Simpson, E. L. 2000, Geology, 28, 831 [NASA ADS] [CrossRef] [Google Scholar]
  35. Fang, J., Wu, H., Fang, Q., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 540, 109530 [NASA ADS] [CrossRef] [Google Scholar]
  36. Farhat, M. A., & Touma, J. R. 2021, MNRAS, 507, 6078 [NASA ADS] [CrossRef] [Google Scholar]
  37. Farhat, M., Laskar, J., & Boué, G. 2022, J. Geophys. Res. Solid Earth, 127, e2021JB023 [CrossRef] [Google Scholar]
  38. Farrell, W. 1972, Rev. Geophys., 10, 761 [NASA ADS] [CrossRef] [Google Scholar]
  39. Fienga, A., Deram, P., Di Ruscio, A., et al. 2021, Notes Scientifiques et Techniques de l’Institut de Mécanique Céleste, 110 [Google Scholar]
  40. Findley, W. N., Lai, J. S., Onaran, K., & Christensen, R. M. 1977, J. Appl. Mech., 44, 364 [NASA ADS] [CrossRef] [Google Scholar]
  41. Fox-Kemper, B., Ferrari, R., & Pedlosky, J. 2003, J. Phys. Oceanogr., 33, 478 [NASA ADS] [CrossRef] [Google Scholar]
  42. Garrett, C., & Munk, W. 1971, Deep Sea Res. Oceanogr. Abstracts, 18, 493 [CrossRef] [Google Scholar]
  43. Gent, P. R., & McWilliams, J. C. 1983, Dyn. Atmos. Oceans, 7, 67 [CrossRef] [Google Scholar]
  44. Gerkema, T., & Zimmerman, J. 2008, in Lecture Notes (Texel: Royal NIOZ), 207 [Google Scholar]
  45. Gerstenkorn, H. 1967, Icarus, 7, 160 [NASA ADS] [CrossRef] [Google Scholar]
  46. Goldreich, P. 1966, Rev. Geophys., 4, 411 [NASA ADS] [CrossRef] [Google Scholar]
  47. Green, J., Huber, M., Waltham, D., Buzan, J., & Wells, M. 2017, Earth Planet. Sci. Lett., 461, 46 [CrossRef] [Google Scholar]
  48. Griffiths, S. D., & Peltier, W. R. 2009, J. Clim., 22, 2905 [NASA ADS] [CrossRef] [Google Scholar]
  49. Gurnis, M., Turner, M., Zahirovic, S., et al. 2012, Comput. Geosci., 38, 35 [NASA ADS] [CrossRef] [Google Scholar]
  50. Han, L., & Huang, R. X. 2020, J. Phys. Oceanogr., 50, 679 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hansen, K. S. 1982, Rev. Geophys., 20, 457 [NASA ADS] [CrossRef] [Google Scholar]
  52. Hawkesworth, C., Cawood, P. A., & Dhuime, B. 2020, Front. Earth Sci., 8 [CrossRef] [Google Scholar]
  53. Hough, S. S. 1898, Philos.Trans. R. Soc. London Ser. A, 191, 139 [NASA ADS] [CrossRef] [Google Scholar]
  54. Huang, H., Gao, Y., Jones, M. M., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 550, 109735 [NASA ADS] [CrossRef] [Google Scholar]
  55. Johnson, B. W., & Wing, B. A. 2020, Nat. Geosci., 13, 243 [NASA ADS] [CrossRef] [Google Scholar]
  56. Kaula, W. M. 2013, Theory of Satellite Geodesy: Applications of Satellites to Geodesy (Courier Corporation) [Google Scholar]
  57. Klatt, J. M., Chennu, A., Arbic, B. K., Biddanda, B., & Dick, G. J. 2021, Nat. Geosci., 14, 564 [NASA ADS] [CrossRef] [Google Scholar]
  58. Lantink, M., Davies, J., & Hilgen, F. 2021, Nat. Rev., 12, 369 [Google Scholar]
  59. Laskar, J. 2005, Celestial Mech. Dyn. Astron., 91, 351 [NASA ADS] [CrossRef] [Google Scholar]
  60. Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Lau, H. C., & Faul, U. H. 2019, Earth Planet. Sci. Lett., 508, 18 [CrossRef] [Google Scholar]
  62. Lau, H. C., Yang, H.-Y., Tromp, J., et al. 2015, Geophys. J. Int., 202, 1392 [NASA ADS] [CrossRef] [Google Scholar]
  63. Lau, H. C., Faul, U., Mitrovica, J. X., et al. 2016a, Geophys. J. Int., 208, 368 [Google Scholar]
  64. Lau, H. C., Mitrovica, J. X., Austermann, J., et al. 2016b, J. Geophys. Res. Solid Earth, 121, 6991 [NASA ADS] [CrossRef] [Google Scholar]
  65. Lee, U., & Saio, H. 1997, ApJ, 491, 839 [Google Scholar]
  66. Levrard, B., & Laskar, J. 2003, Geophys. J. Int., 154, 970 [NASA ADS] [CrossRef] [Google Scholar]
  67. Longuet-Higgins, M. S. 1968, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci., 262, 511 [NASA ADS] [CrossRef] [Google Scholar]
  68. Longuet-Higgins, M. S., & Pond, G. S. 1970, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci., 266, 193 [NASA ADS] [Google Scholar]
  69. MacDonald, G. 1967, Proc. R. Soc. London Ser. A. Math. Phys. Sci., 296, 298 [Google Scholar]
  70. Matsuyama, I. 2014, Icarus, 242, 11 [NASA ADS] [CrossRef] [Google Scholar]
  71. Matthews, K. J., Maloney, K. T., Zahirovic, S., et al. 2016, Global Planet. Change, 146, 226 [NASA ADS] [CrossRef] [Google Scholar]
  72. Maurice, M., Tosi, N., Schwinger, S., Breuer, D., & Kleine, T. 2020, Sci. Adv., 6, 28 [CrossRef] [Google Scholar]
  73. Mavromatis, H., & Alassar, R. 1999, Appl. Math. Lett., 12, 101 [CrossRef] [Google Scholar]
  74. Merdith, A. S., Williams, S. E., Collins, A. S., et al. 2021, Earth-Sci. Rev., 214, 103477 [NASA ADS] [CrossRef] [Google Scholar]
  75. Meyers, S. R., & Malinverno, A. 2018, Proc. Nat. Acad. Sci., 115, 6363 [NASA ADS] [CrossRef] [Google Scholar]
  76. Mignard, F. 1979, Moon Planets, 20, 301 [Google Scholar]
  77. Mojzsis, S. J., Arrhenius, G., McKeegan, K., et al. 1996, Nature, 384, 55 [NASA ADS] [CrossRef] [Google Scholar]
  78. Motoyama, M., Tsunakawa, H., & Takahashi, F. 2020, Icarus, 335, 113382 [NASA ADS] [CrossRef] [Google Scholar]
  79. 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 [Google Scholar]
  80. Müller, M. 2008b, Ocean Model., 20, 207 [CrossRef] [Google Scholar]
  81. Munk, W. H., & MacDonald, G. J. 1960, The Rotation of the Earth: A Geophysical Discussion (Cambridge University Press) [Google Scholar]
  82. Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
  83. Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
  84. Ooe, M. 1989, J. Phys. Earth, 37, 345 [CrossRef] [Google Scholar]
  85. Palmer, T., Shutts, G., & Swinbank, R. 1986, Q. J. R. Meteorol. Soc., 112, 1001 [NASA ADS] [CrossRef] [Google Scholar]
  86. Peck, W. H., Valley, J. W., Wilde, S. A., & Graham, C. M. 2001, Geochim. Cosmochim. Acta, 65, 4215 [NASA ADS] [CrossRef] [Google Scholar]
  87. Petit, G., & Luzum, B. 2010, IERS conventions (2010), Tech. rep. (France: Bureau International des Poids et mesures sevres) [Google Scholar]
  88. Platzman, G. W. 1983, Science, 220, 602 [NASA ADS] [CrossRef] [Google Scholar]
  89. Platzman, G. W. 1984, J. Phys. Oceanogr., 14, 1532 [NASA ADS] [CrossRef] [Google Scholar]
  90. Proudman, J. 1920a, Proc. London Math. Soc., 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
  91. Proudman, J. 1920b, Proc. London Math. Soc., 2, 51 [NASA ADS] [CrossRef] [Google Scholar]
  92. Regge, T. 1958, Il Nuovo Cimento (1955–1965), 10, 544 [CrossRef] [Google Scholar]
  93. Renaud, J. P., & Henning, W. G. 2018, ApJ, 857, 98 [Google Scholar]
  94. Riley, K. F., Hobson, M. P., & Bence, S. J. 1999, Mathematical Methods for Physics and Engineering [Google Scholar]
  95. Ross, M., & Schubert, G. 1989, J. Geophys. Res. Solid Earth, 94, 9533 [NASA ADS] [CrossRef] [Google Scholar]
  96. Rufu, R., & Canup, R. M. 2020, J. Geophys. Res. Planets, 125, 8 [NASA ADS] [CrossRef] [Google Scholar]
  97. Sonett, C., & Chan, M. A. 1998, Geophys. Res. Lett., 25, 539 [NASA ADS] [CrossRef] [Google Scholar]
  98. Sørensen, A. L., Nielsen, A. T., Thibault, N., et al. 2020, Earth Planet. Sci. Lett., 548, 116475 [CrossRef] [Google Scholar]
  99. Strauss, W. A. 2007, in Partial Differential Equations: An introduction (John Wiley& Sons) [Google Scholar]
  100. Sun, C., Xu, W., Cawood, P. A., et al. 2019, Sci. Rep., 9, 1 [Google Scholar]
  101. Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534 [Google Scholar]
  102. Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A&A, 630, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Touma, J., & Wisdom, J. 1994, AJ, 108, 1943 [NASA ADS] [CrossRef] [Google Scholar]
  104. Touma, J., & Wisdom, J. 1998, AJ, 115, 1653 [NASA ADS] [CrossRef] [Google Scholar]
  105. Touma, J., & Wisdom, J. 2001, AJ, 122, 1030 [NASA ADS] [CrossRef] [Google Scholar]
  106. Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706 [Google Scholar]
  107. Tyler, R. 2011, Icarus, 211, 770 [NASA ADS] [CrossRef] [Google Scholar]
  108. Tyler, R. H. 2021, Planet. Sci. J., 2, 70 [CrossRef] [Google Scholar]
  109. Vallis, G. K. 2017, Atmospheric and Oceanic Fluid Dynamics (Cambridge University Press) [CrossRef] [Google Scholar]
  110. Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1988, Quantum Theory of Angular Momentum (World Scientific) [CrossRef] [Google Scholar]
  111. Walker, J. C. G., & Zahnle, K. J. 1986, Nature, 320, 600 [NASA ADS] [CrossRef] [Google Scholar]
  112. Waltham, D. 2015, J. Sediment. Res., 85, 990 [NASA ADS] [CrossRef] [Google Scholar]
  113. Wang, H., Boyd, J. P., & Akmaev, R. A. 2016, Geosci. Model Dev., 9, 1477 [Google Scholar]
  114. Watterson, I. G. 2001, J. Atmos. Oceanic Technol., 18, 691 [NASA ADS] [CrossRef] [Google Scholar]
  115. Webb, D. 1973, in Deep Sea Research and Oceanographic Abstracts (Elsevier), 20, 847 [NASA ADS] [CrossRef] [Google Scholar]
  116. Webb, D. 1980, Geophys. J. Int., 61, 573 [NASA ADS] [CrossRef] [Google Scholar]
  117. Webb, D. 1982, Geophys. J. Int., 70, 261 [NASA ADS] [CrossRef] [Google Scholar]
  118. Wilde, S. A., Valley, J. W., Peck, W. H., & Graham, C. M. 2001, Nature, 409, 175 [NASA ADS] [CrossRef] [Google Scholar]
  119. Williams, G. E. 1990, J. Phys. Earth, 38, 475 [CrossRef] [Google Scholar]
  120. Williams, G. E. 1997, Geophys. Res. Lett., 24, 421 [NASA ADS] [CrossRef] [Google Scholar]
  121. Williams, G. E. 2000, Rev. Geophys., 38, 37 [NASA ADS] [CrossRef] [Google Scholar]
  122. Williams, J. G., & Boggs, D. H. 2016, Celestial Mech. Dyn. Astron., 126, 89 [NASA ADS] [CrossRef] [Google Scholar]
  123. Wood, R., Liu, A. G., Bowyer, F., et al. 2019, Nat. Ecol. Evol., 3, 528 [CrossRef] [Google Scholar]
  124. Zahel, W. 1980, Phys. Earth Planet. Inter., 21, 202 [CrossRef] [Google Scholar]
  125. Zahnle, K. J., Lupu, R., Dobrovolskis, A., & Sleep, N. H. 2015, Earth Planetary Sci. Lett., 427, 74 [NASA ADS] [CrossRef] [Google Scholar]
  126. Zhong, Y., Wu, H., Fan, J., et al. 2020, Palaeogeogr. Palaeoclimatol. Palaeoecol., 540, 109520 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Orbital dynamics

For the reconstruction of the 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 intervening players. In particular, the contribution of the eccentricity tides becomes significant when the orbit of the Moon was highly eccentric. Accounting for lunar tides and lunar 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 ≤ 0.1) (Daher et al. 2021). By attaining a highly eccentric lunar orbit, is possible through evection resonance1. Touma & Wisdom (1998) studied this regime and showed that capture into such a resonance could have been encountered for aM ≈ 4.6RE, exciting the lunar eccentricity to e ≈ 0.5. However, the timescale of capture and escape from this evection resonance is 104 ∼ 105 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 aM <  20RE (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 aM dropping rapidly from 30RE 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 aM >  30RE. Thus, as much as early stage mechanisms are essential to constrain the formation scenarios of the system, the robustness of the runaway evolution beyond 30RE 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 (Ćuk et al. 2016).

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

d L Ω dt = − ( T M + T S ) , $$ \begin{aligned} \frac{dL_\Omega }{dt}&= -\left(\mathcal{T} _{\rm M} + \mathcal{T} _{\rm S}\right) \ , \end{aligned} $$(A.1)

d L M dt = T M , $$ \begin{aligned} \frac{dL_{\rm M}}{dt}&= \mathcal{T} _{\rm M} \ , \end{aligned} $$(A.2)

where 𝒯M is the lunar semi-diurnal tidal torque coupling between the oceanic and the solid response of the Earth, and 𝒯S is its solar counterpart. The orbital angular momentum of the Moon L M = β G ( M E + M M ) a M $ L_M=\beta\sqrt{G(M_{\mathrm{E}}+M_{\mathrm{M}}) a_{\mathrm{M}}} $, where β = MEMM/(ME + MM) is the Earth-Moon system’s reduced mass. The rotational angular momentum of the Earth is defined as LΩ = C(Ω)Ω, with the time-varying principal moment of inertia given by (Goldreich 1966)

C ( Ω ) = C ( Ω 0 ) + 2 k 2 f R E 5 9 G ( Ω 2 − Ω 0 2 ) . $$ \begin{aligned} C(\Omega )= C(\Omega _0) + \frac{2 k_2^f R_{\rm E}^5 }{9G}(\Omega ^2-\Omega _0^2) . \end{aligned} $$(A.3)

Here, k 2 f $ 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 aM = 3RE. The tidal torque computation is coupled to the orbital integrator and is computed simultaneously at each step. It takes the model parameters (H, σR) as input, and the system’s variables, aM and Ω, to compute the tidal frequency and, consequently, the coupled tidal response.

Once the lunar semi-major axis (aM) and the rotation speed of the Earth (Ω) are determined, we compute the obliquity of the Earth (ϵ) and the precession frequency (p) as derived quantities (Laskar et al. 2004). Starting with Equations (40) and (46) from Correia & Laskar (2010) in the case of zero eccentricity, we obtain:

d ϵ dt = Kn C ( Ω ) Ω sin ϵ ( Ω 2 n M cos ϵ − 1 ) , $$ \begin{aligned} \frac{d\epsilon }{dt} = \frac{K n}{C(\Omega )\Omega } \sin \epsilon \left( \frac{\Omega }{2n_{\rm M}} \cos \epsilon -1 \right) , \end{aligned} $$(A.4)

and

d a M dt = 2 K β a M ( Ω n M cos ϵ − 1 ) , $$ \begin{aligned} \frac{da_{\rm M}}{dt} = \frac{2K }{\beta a_{\rm M}} \left( \frac{\Omega }{n_{\rm M}}\cos \epsilon -1 \right) , \end{aligned} $$(A.5)

that is

d ϵ d a M = β n M a M 4 C ( Ω ) Ω sin ϵ Ω cos ϵ − 2 n M Ω cos ϵ − n M . $$ \begin{aligned} \frac{d \epsilon }{d a_{\rm M}} = \frac{\beta n_{\rm M} a_{\rm M}}{4C(\Omega )\Omega } \sin \epsilon \frac{ \Omega \cos \epsilon -2n_{\rm M} }{ \Omega \cos \epsilon -n_{\rm M}} . \end{aligned} $$(A.6)

We note that the tidal response parameter, K, disappears from the equations. This would also be the case if K depended on Ω. The obliquity evolution equation (A.6) is integrated using the values of aM and Ω that result from the tidal flows and orbital dynamics coupled system. The precession frequency, p, is then derived using Equations (6) and (8) from Laskar et al. (2004) with zero eccentricity and inclination, that is:

p = 3 2 ( G M S a E 3 + G M M a M 3 ) E d ( Ω 0 ) Ω Ω 0 2 cos ϵ . $$ \begin{aligned} p=\frac{3}{2}\left(\frac{GM_{\rm S}}{a_{\rm E}^3} + \frac{GM_{\rm M}}{a_{\rm M}^3}\right)E_d(\Omega _0)\frac{\Omega }{\Omega _0^2}\cos \epsilon . \end{aligned} $$(A.7)

In (A.6) and (A.7), the constant values taken for the Earth’s radius RE, the gravitational constant of the Moon GMM, and the Sun GMS, the mass ratio ME/MM, the rotational velocity Ω0, the Earth’s semi-major axis aE, and the inertia parameter C ( Ω 0 ) / M E R E 2 $ C(\Omega_0)/M_{\rm E} R_{\rm E}^2 $ are adopted from INPOP21 (Fienga et al. 2021). The dynamical ellipticity at the origin of date, Ed(Ω0)=0.003243, is determined from the initial conditions for the obliquity (ϵ0) and precession (p0) 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 aM evolution of Figure 3, the Laplace regime transition occurs very early in the evolution (t >  4Ga); thus, in Figure 6, we plot the evolution of the precession frequency and obliquity between the present and 3.5 Ga. We also add the geological inferences of the precession frequency as scatter on the curve; in the case of cyclostratigraphy, this is the direct observable (Table D.1).

Appendix B: Continental drift and oceanic geometry shifting

Our dynamical integrator allows for variations in the oceanic geometry, be it a variation in the position of the oceanic hemisphere or a shift between the hemispheric and global oceanic configurations. In the combined model, the first phase (Figure 4) starts with the center of the continental cap following the evolution of the geographic center over the recent billion years. For this purpose, we adopted a recent model that reconstructs a kinetically continuous history of plate tectonics (Merdith et al. 2021). The geographic center is traced by computing the surface projection of the “barycenter” of the continental distribution. This allows for a higher level of realism in oceanic modeling. Beyond 1 Ga (in part due to the lack of plate tectonic data), the model continues with the position of the ocean at 1 Ga. To postprocess the continental drift evolution and to produce the time-sliced sketches of Figure 1, we used the GPlates open-source reconstruction software (Boyden et al. 2011; Gurnis et al. 2012).

thumbnail 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).

At 1.5 Ga, the integrator starts computing simultaneously the tidal response of a global oceanic geometry, with uniform thickness Hglobal = 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 tswitch. 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 aM (Figure 3) without any discontinuities and modeling artifacts. For the misfit minimum of our combined model, we have tswitch = 3.25 Ga.

Appendix C: Parameter fits

To construct the misfit surfaces of Figure 2, we computed the evolution for each pair of (H, σR) on the two-dimensional grid. The present rate of lunar recession, ̇a0, and the impact time, tf, are then extracted for each evolution sample and the mean square weighted deviation, χ2, (Table 1) is then computed as:

χ 2 = 1 2 [ ( a 0 − a 0 LLR σ LLR ) 2 + ( t f − t f geo σ geo ) 2 ] , $$ \begin{aligned} \chi ^2 = \frac{1}{2} \left[\left( \frac{{a}_0-{a}_0^\mathrm{LLR}}{\sigma ^\mathrm{LLR}} \right)^2 + \left( \frac{t_{\rm f}-t_{\rm f}^\mathrm{geo}}{\sigma ^\mathrm{geo}} \right)^2\right]\ , \end{aligned} $$(C.1)

where we use Lunar Laser Ranging (LLR) estimates of lunar orbital recession (Williams & Boggs 2016): a ˙ 0 L L R ± σ L L R = 38.30 ± 0.08 $ \dot{a}_0^{\rm LLR}\pm\sigma^{\rm LLR}= 38.30 \pm 0.08 $ mm/year; as well as geochemical estimates of lunar formation time (Maurice et al. 2020): t f geo ± σ geo = 4.425 ± 0.025 $ t_{\mathrm{f}}^{\mathrm{geo}}\pm\sigma^{\mathrm{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 C.1.

Misfit analysis summary.

We evaluated the uncertainties on the fitted parameters from those on the observables following the standard propagation of uncertainty method. Because of the absence of correlation between the two data ȧ0 and tf, the entries of the variance matrix,

Σ = [ var ( H ) cov ( H , σ R ) [ 0.5 e m ] cov ( H , σ R ) var ( σ R ) ] , $$ \begin{aligned} \Sigma = \begin{bmatrix} {{\,\mathrm{var}\,}}(H)&{{\,\mathrm{cov}\,}}(H,\sigma _{\rm R}) \\[0.5em] {{\,\mathrm{cov}\,}}(H,\sigma _{\rm R})&{{\,\mathrm{var}\,}}(\sigma _{\rm R}) \end{bmatrix}\,, \end{aligned} $$(C.2)

are given by

var ( H ) = ( ∂ H ∂ a ˙ 0 ) 2 ( σ LLR ) 2 + ( ∂ H ∂ t f ) 2 ( σ geo ) 2 , $$ \begin{aligned} {{\,\mathrm{var}\,}}(H)&= \left(\frac{\partial H}{\partial \dot{a}_0}\right)^2\left(\sigma ^\mathrm{LLR}\right)^2 + \left(\frac{\partial H}{\partial t_\mathrm{f} }\right)^2\left(\sigma ^\mathrm{geo}\right)^2\,, \end{aligned} $$(C.3a)

var ( σ R ) = ( ∂ σ R ∂ a ˙ 0 ) 2 ( σ LLR ) 2 + ( ∂ σ R ∂ t f ) 2 ( σ geo ) 2 , $$ \begin{aligned} {{\,\mathrm{var}\,}}(\sigma _{\rm R})&= \left(\frac{\partial \sigma _{\rm R}}{\partial \dot{a}_0}\right)^2\left(\sigma ^\mathrm{LLR}\right)^2 + \left(\frac{\partial \sigma _{\rm R}}{\partial t_\mathrm{f} }\right)^2\left(\sigma ^\mathrm{geo}\right)^2\,, \end{aligned} $$(C.3b)

cov ( H , σ R ) = ∂ H ∂ a ˙ 0 ∂ σ R ∂ a ˙ 0 ( σ LLR ) 2 + ∂ H ∂ t f ∂ σ R ∂ t f ( σ geo ) 2 . $$ \begin{aligned} {{\,\mathrm{cov}\,}}(H,\sigma _{\rm R})&= \frac{\partial H}{\partial \dot{a}_0}\frac{\partial \sigma _{\rm R}}{\partial \dot{a}_0} \left(\sigma ^\mathrm{LLR}\right)^2 + \frac{\partial H}{\partial t_\mathrm{f} }\frac{\partial \sigma _{\rm R}}{\partial t_\mathrm{f} }\left(\sigma ^\mathrm{geo}\right)^2\,. \end{aligned} $$(C.3c)

The partial derivatives entering in these formulae are computed numerically from the fit of (ȧ0 ± σLLR,tf) and (ȧ0,tf ± σgeo). The marginal uncertainties on parameters H and σR are σH = var(H)1/2 = 32.75 m and σσR = var(σR)1/2 = 0.2631 × 10−6 s−1, respectively. We used the variance matrix Σ to evaluate the 2σ-confidence ellipsoid around the best-fit parameters (H, σR). The Earth-Moon distance, aM, and the length of the day LOD

have been integrated for 25 pairs of parameters (H, σR) chosen at the boundary of this 2σ-confidence region. Their envelope represents the 2σ-uncertainty area plotted in shaded blue in (Figure 3) and (Figs. 5 and 6).

Appendix D: Geological data

In tables D.1 and D.2, we compile geological data sets that provide historical snapshots of the past rotational state of the Earth and the lunar orbital distance.

Table D.1.

Cyclostratigraphic data.

Table D.2.

Tidal rhythmites data.

Table D.3.

Values of constant parameters used in the numerical implementation of the theory.

Appendix E: Tidal response of a hemispherical ocean

This appendix is aimed at obtaining a computation of the tidal response of a hemispherical ocean on the surface of the Earth. The formalism is heavily based on earlier works (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 Ω = Ω z ̂ $ \boldsymbol{\Omega}=\Omega \hat{z} $, Ω being the Earth’s spin rate and ẑ as the unit vector along its figure axis. In this frame, we used the spherical coordinates (r, θ, λ) denoting the radius, the co-latitude, and the longitude respectively, and their corresponding unit vectors ( r ̂ , θ ̂ , λ ̂ ) $ (\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):

∂ t u + σ R u + f × u + g ∇ ζ = g ∇ ζ eq , $$ \begin{aligned}&\partial _t\boldsymbol{u}+\sigma _{\rm R}\boldsymbol{u}+\boldsymbol{f}\boldsymbol{\times }\boldsymbol{u}+g\boldsymbol{\nabla }\zeta = g\boldsymbol{\nabla }\zeta _\mathrm{eq} ,\end{aligned} $$(E.1a)

∂ t ζ + ∇ · ( H u ) = 0 , $$ \begin{aligned}&\partial _t\zeta +\boldsymbol{\nabla }\cdot {\left( H\boldsymbol{u} \right)}=0, \end{aligned} $$(E.1b)

where u = u θ θ ̂ + u λ λ ̂ $ \boldsymbol{u}= u_\theta\hat{\theta}+u_\lambda \hat{\lambda} $ is the horizontal velocity field, g is the gravitational acceleration at the surface, ζ is the oceanic depth variation, ζeq is the equilibrium depth variation, H is the uniform oceanic thickness (the first of only two free parameters in our model), and σR is the Rayleigh (or linear) drag frequency (Matsuyama 2014; 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, σ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 Ω cos θ r ̂ , $$ \begin{aligned} \boldsymbol{f}=2\Omega \cos \theta \hat{r}, \end{aligned} $$(E.2)

the horizontal gradient operator ∇ is defined as:

∇ = R − 1 [ θ ̂ ∂ θ + λ ̂ ( sin θ ) − 1 ∂ λ ] , $$ \begin{aligned} \boldsymbol{\nabla } = R^{-1} \left[ \hat{\theta } \partial _{\theta } + \hat{\lambda } \left(\sin \theta \right)^{-1} \partial _{\lambda } \right], \end{aligned} $$(E.3)

and the horizontal divergence of the velocity field ∇ ⋅ u as:

∇ · u = ( R sin θ ) − 1 [ ∂ θ ( sin θ u θ ) + ∂ λ u λ ] , $$ \begin{aligned} \boldsymbol{\nabla } \cdot \boldsymbol{u} = \left( R \sin \theta \right)^{-1} \left[ \partial _\theta \left( \sin \theta u_\theta \right) + \partial _\lambda u_\lambda \right] , \end{aligned} $$(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 = ∂tx, where x is the horizontal tidal displacement field, we have:

[ ∂ t 2 + ( σ R + f × ) ∂ t ] x + g ( ∇ ζ − ∇ ζ eq ) = 0 . $$ \begin{aligned} \left[\partial _t^2 +(\sigma _{\rm R}+\boldsymbol{f}\boldsymbol{\times })\partial _t \right]\boldsymbol{x}+g\left(\boldsymbol{\nabla }\zeta -\boldsymbol{\nabla }{\zeta _\mathrm{eq} }\right) =0. \end{aligned} $$(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:

x = ∇ Φ + ∇ Ψ × r ̂ , $$ \begin{aligned} \boldsymbol{x}=\boldsymbol{\nabla }\Phi +\boldsymbol{\nabla }\Psi \boldsymbol{\times }\hat{r}, \end{aligned} $$(E.6)

where ∇Φ is a curl-free vector field (∇×(∇Φ) = 0) and ∇ Ψ × r ̂ $ \boldsymbol{\nabla}\Psi \boldsymbol{\times} \hat{r} $ is a divergence-free vector field ( ∇ · ( ∇ Ψ × r ̂ ) = 0 $ \boldsymbol{\nabla} \cdot \left( \boldsymbol{\nabla}\Psi\boldsymbol{\times}\hat{r} \right) = 0 $). In the above equation, we introduce the divergent displacement potential Φ and the rotational displacement streamfunction Ψ (e.g., Gent & McWilliams 1983; Webb 1980; Tyler 2011), with the latter accounting for the vortical component of the tidal displacement field (e.g., Vallis 2017). As discussed by 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 x · n̂ = 0, where n̂ designates the outward pointing unit vector defining the normal to the coast. Following Webb (1980, 1982), we assume that both components of the flux satisfy this condition, namely:

n ̂ · ∇ Φ = 0 , n ̂ · ( ∇ Ψ × r ̂ ) = 0 . $$ \begin{aligned}&\hat{n} \cdot \boldsymbol{\nabla } \Phi = 0,&\hat{n} \cdot \left( \boldsymbol{\nabla } \Psi \boldsymbol{\times } \hat{r} \right) = 0. \end{aligned} $$(E.7)

We note that the second condition of the above equation can be rewritten as ( r ̂ × n ̂ ) · ∇ Ψ = 0 $ \left( \hat{r} \boldsymbol{\times} \hat{n} \right) \cdot \boldsymbol{\nabla} \Psi = 0 $, which implies that Ψ is constant along the coastline (in the following, we set Ψ = 0 at the oceanic boundary). This condition thus means that the coastline corresponds to a closed streamfunction contour, which depicts a distinct gyre of the tidal flow.

Although arbitrary, the assumption that both components of the flux satisfy the impermeability condition has been widely used to study the dynamics of ocean basins because of its convenience relative to other possible conditions (e.g., Gent & McWilliams 1983; Watterson 2001; Han & Huang 2020). In particular, this assumption provides a unique decomposition apart from an arbitrary additive constant to each function, Φ and Ψ. Moreover, the second condition given by Eq. (E.7) enforces the orthogonality of the curl-free and divergence-free components of the tidal flow. By combining together the identity ∇ · ( Φ ∇ Ψ × r ̂ ) = ( ∇ Ψ × r ̂ ) · ∇ Φ $ \boldsymbol{\nabla} \cdot \left( \Phi \boldsymbol{\nabla} \Psi \boldsymbol{\times} \hat{r} \right) = \left( \boldsymbol{\nabla} \Psi \boldsymbol{\times} \hat{r} \right) \cdot \boldsymbol{\nabla} \Phi $ and Gauss’ theorem (e.g., Arfken & Weber 1999):

∫ O ∇ · ( Φ ∇ Ψ × r ̂ ) d A = ∮ ∂ O Φ ( ∇ Ψ × r ̂ ) · n ̂ d ℓ , $$ \begin{aligned} \int _{\mathcal{O} } \boldsymbol{\nabla } \cdot \left( \Phi \boldsymbol{\nabla } \Psi \boldsymbol{\times } \hat{r} \right) d A = \oint _{\partial \mathcal{O} } \Phi \left( \boldsymbol{\nabla } \Psi \boldsymbol{\times } \hat{r} \right) \cdot \hat{n} \, d \ell , \end{aligned} $$(E.8)

with dA and dℓ being the infinitesimal area element of the hemispherical oceanic domain, 𝒪, and length element of the coastline, ∂𝒪, respectively, we obtain

∫ O ( ∇ Φ ) · ( ∇ Ψ × r ̂ ) d A = ∮ ∂ O Φ ( ∇ Ψ × r ̂ ) · n ̂ d ℓ . $$ \begin{aligned} \int _{\mathcal{O} } \left( \boldsymbol{\nabla } \Phi \right) \cdot \left( \boldsymbol{\nabla } \Psi \boldsymbol{\times } \hat{r} \right) d A = \oint _{\partial \mathcal{O} } \Phi \left( \boldsymbol{\nabla } \Psi \boldsymbol{\times } \hat{r} \right) \cdot \hat{n} \, d \ell . \end{aligned} $$(E.9)

As the second condition of Eq. (E.7) enforces ( ∇ Ψ × r ̂ ) · n ̂ = 0 $ \left( \boldsymbol{\nabla} \Psi \boldsymbol{\times} \hat{r} \right) \cdot \hat{n} = 0 $ along the coastline, it follows that

∫ O ( ∇ Φ ) · ( ∇ Ψ × r ̂ ) d A = 0 , $$ \begin{aligned} \int _{\mathcal{O} } \left( \boldsymbol{\nabla } \Phi \right) \cdot \left( \boldsymbol{\nabla } \Psi \boldsymbol{\times } \hat{r} \right) d A = 0, \end{aligned} $$(E.10)

meaning that the components ∇Φ and ∇ Ψ × r ̂ $ \boldsymbol{\nabla} \Psi \boldsymbol{\times} \hat{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 Φ and Ψ are expanded in terms of complete sets of eigenfunctions over the domain, 𝒪, such that

Φ ( θ , λ , t ) = ∑ r = 1 ∞ p r ( t ) ϕ r ( θ , λ ) , $$ \begin{aligned} \Phi (\theta ,\lambda ,t)= \sum _{r=1}^{\infty }p_r(t)\phi _r(\theta ,\lambda ), \end{aligned} $$(E.11)

Ψ ( θ , λ , t ) = ∑ r = 1 ∞ p − r ( t ) ψ r ( θ , λ ) . $$ \begin{aligned} \Psi (\theta ,\lambda ,t)= \sum _{r=1}^{\infty }p_{-r}(t)\psi _r(\theta ,\lambda ). \end{aligned} $$(E.12)

The eigenfunctions (ϕr, ψr) satisfy, over the oceanic domain (𝒪), the Helmholtz equations (e.g., Riley et al. (1999), Chapter 21):

∇ 2 ϕ r + μ r ϕ r = 0 , $$ \begin{aligned} \boldsymbol{\nabla }^2\phi _r+\mu _r\phi _r&=0, \end{aligned} $$(E.13)

∇ 2 ψ r + ν r ψ r = 0 , $$ \begin{aligned} \boldsymbol{\nabla }^2\psi _r+\nu _r\psi _r&=0, \end{aligned} $$(E.14)

and, along the coastline (∂𝒪), the boundary conditions given by Eq. (E.7):

n ̂ · ∇ ϕ r = 0 , ψ r = 0 , $$ \begin{aligned} \hat{n} \cdot \boldsymbol{\nabla }\phi _r=0, \psi _r=0, \end{aligned} $$(E.15)

where we introduced the horizontal Laplacian:

∇ 2 = ( R sin θ ) − 2 [ sin θ ∂ θ ( sin θ ∂ θ ) + ∂ λ λ ] , $$ \begin{aligned} \boldsymbol{\nabla }^2 =\left( R \sin \theta \right)^{-2} \left[ \sin \theta \, \partial _\theta \left( \sin \theta \, \partial _\theta \right) + \partial _{\lambda \lambda } \right] , \end{aligned} $$(E.16)

and the real eigenvalues, μr and νr, associated with the eigenfunctions ϕr and ψr, respectively. We note that the eigenfunctions are normalized such that

∫ O ϕ r ϕ s d A = ∫ O ψ r ψ s d A = δ rs , $$ \begin{aligned} \int _\mathcal{O} \phi _r\phi _sdA=\int _\mathcal{O} \psi _r\psi _sdA= \delta _{rs}, \end{aligned} $$(E.17)

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

ϕ r = α n , m R P n m ( cos θ ) cos m λ , $$ \begin{aligned} \phi _r&=\frac{\alpha _{n,m}}{R}P_n^m(\cos \theta )\cos m\lambda , \end{aligned} $$(E.18)

ψ r = α n , m R P n m ( cos θ ) sin m λ , $$ \begin{aligned} \psi _r&=\frac{\alpha _{n,m}}{R}P_n^m(\cos \theta )\sin m\lambda , \end{aligned} $$(E.19)

with eigenvalues of μr = νr = n(n + 1)/R2 and the normalization coefficient of:

α n , m = 2 n + 1 π ( n − m ) ! ( n + m ) ! 1 1 + δ m 0 . $$ \begin{aligned} \alpha _{n,m} = \sqrt{\frac{2n+1}{\pi }\frac{(n-m)!}{(n+m)!}\frac{1}{1+\delta _{m0}}}{.} \end{aligned} $$(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 ϕr Eq. (E.18), n ∈ ℕ and m = 0, 1, ..., n while in the expression of ψr Eq. (E.19), n ∈ ℕ* and m = 1, 2, ..., n. By convention, we set ψ0 = 0 hereafter. Figure E.1 shows the eigenfunctions ϕr and ψr for 1 ≤ m ≤ n ≤ 4 and the streamlines of the associated tidal flows.

thumbnail Fig. E.1.

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

The eigenfunctions (ϕr, ψ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 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

ζ = H ∑ r = 1 ∞ μ r p r ϕ r . $$ \begin{aligned} \zeta =H\sum _{r=1}^\infty \mu _r p_r\phi _r. \end{aligned} $$(E.21)

What is left to complete the solution is finding the coefficients pr and p−r by substituting the series expansions in the momentum equation (E.1a) and multiplying by ∇ϕr and ∇ ψ r × r ̂ $ \boldsymbol{\nabla}\psi_r\boldsymbol{\times}\hat{r} $, then integrating over the oceanic area. Starting with the former, we get

∑ s = 0 ∞ ( ∂ t 2 p s + σ R ∂ t p s + g H p s μ s − g ζ eq , s ) ∇ ϕ s · ∇ ϕ r + ( ∂ t 2 p − s + σ R ∂ t p − s ) ( ∇ ψ s × r ̂ ) · ∇ ϕ r + ∂ t p s ( f × ∇ ϕ s ) · ∇ ϕ r + ∂ t p − s [ f × ( ∇ ψ s × r ̂ ) ] · ∇ ϕ r = 0 . $$ \begin{aligned} \nonumber \sum _{s=0}^\infty&\left(\partial _t^2p_s+\sigma _{\rm R}\partial _t p_s +gHp_s\mu _s -g{\zeta }_{\mathrm{eq} ,s}\right)\boldsymbol{\nabla }\phi _s\cdot \boldsymbol{\nabla }\phi _r \\\nonumber&+\left(\partial _t^2p_{-s}+\sigma _{\rm R}\partial _tp_{-s}\right)\left(\boldsymbol{\nabla }\psi _s\boldsymbol{\times }\hat{r}\right)\cdot \boldsymbol{\nabla }\phi _r \\ &+ \partial _t p_s\left(\boldsymbol{f}\boldsymbol{\times }\boldsymbol{\nabla }\phi _s\right)\cdot \boldsymbol{\nabla }\phi _r+\partial _t p_{-s}\left[\boldsymbol{f}\boldsymbol{\times }\left(\boldsymbol{\nabla }\psi _s\boldsymbol{\times }\hat{r}\right)\right]\cdot \boldsymbol{\nabla }\phi _r=0. \end{aligned} $$(E.22)

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

∫ O ∇ ϕ s · ∇ ϕ r d A = ∫ ∂ O ϕ s ( ∇ ϕ r · n ̂ ) d ℓ − ∫ O ϕ s ∇ 2 ϕ r d A . $$ \begin{aligned} \int _\mathcal{O} \boldsymbol{\nabla }\phi _s\cdot \boldsymbol{\nabla }\phi _r dA =\int _{\partial \mathcal{O} } \phi _s\left(\boldsymbol{\nabla }\phi _r\cdot \hat{n}\right) d\ell - \int _\mathcal{O} \phi _s\boldsymbol{\nabla }^2\phi _rdA. \end{aligned} $$(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:

∫ O ∇ ϕ s · ∇ ϕ r d A = μ r δ rs . $$ \begin{aligned} \int _\mathcal{O} \boldsymbol{\nabla }\phi _s\cdot \boldsymbol{\nabla }\phi _r dA = \mu _r \delta _{rs}. \end{aligned} $$(E.24)

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

∑ s = 0 ∞ ( ∂ t 2 p s + σ R ∂ t p s + g H p s μ s − g ζ eq , s ) μ r δ r , s + ( ∂ t 2 p − s + σ R ∂ t p − s ) ∫ O ( ∇ ψ s × r ̂ ) · ∇ ϕ r d A + ∂ t p s ∫ O f · ( ∇ ϕ s × ∇ ϕ r ) d A + ∂ t p − s ∫ O ( f · r ̂ ) ( ∇ ϕ r · ∇ ψ s ) d A = 0 . $$ \begin{aligned} \nonumber \sum _{s=0}^\infty&\left(\partial _t^2p_s+\sigma _{\rm R}\partial _t p_s +gHp_s\mu _s -g{\zeta }_{\mathrm{eq} ,s}\right)\mu _r \delta _{r,s} \\\nonumber&+\left(\partial _t^2p_{-s}+\sigma _{\rm R}\partial _tp_{-s}\right) \int _\mathcal{O} \left(\boldsymbol{\nabla }\psi _s\boldsymbol{\times }\hat{r}\right)\cdot \boldsymbol{\nabla }\phi _r \,dA \\\nonumber&+ \partial _t p_s \int _\mathcal{O} \boldsymbol{f}\cdot \left(\boldsymbol{\nabla }\phi _s\boldsymbol{\times }\boldsymbol{\nabla }\phi _r\right)dA\\&+\partial _t p_{-s}\int _\mathcal{O} \left(\boldsymbol{f}\cdot \hat{r}\right) \left(\boldsymbol{\nabla }\phi _r\cdot \boldsymbol{\nabla }\psi _s\right)dA=0. \end{aligned} $$(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

( ∂ t 2 p r + σ R ∂ t p r + g H p r μ r − g ζ eq , r ) μ r − 2 Ω ∑ s = 1 ∞ ∂ t p s ∫ cos θ r ̂ · ( ∇ ϕ r × ∇ ϕ s ) d A + 2 Ω ∑ s = 1 ∞ ∂ t p − s ∫ cos θ ( ∇ ϕ r · ∇ ψ s ) d A = 0 . $$ \begin{aligned} \nonumber&\left(\partial _t^2p_r+\sigma _{\rm R}\partial _t p_r +gHp_r\mu _r -g{\zeta }_{\mathrm{eq} ,r}\right)\mu _r \\\nonumber&- 2\Omega \sum _{s=1}^\infty \partial _t p_s \int \cos \theta \hat{r}\cdot \left(\boldsymbol{\nabla }\phi _r\boldsymbol{\times }\boldsymbol{\nabla }\phi _s\right)dA\\&+2\Omega \sum _{s=1}^\infty \partial _t p_{-s}\int \cos \theta \left(\boldsymbol{\nabla }\phi _r\cdot \boldsymbol{\nabla }\psi _s\right)dA=0. \end{aligned} $$(E.26)

To close the system, we multiply the momentum equation with ∇ ψ r × r ̂ $ \boldsymbol{\nabla}\psi_r\boldsymbol{\times}\hat{r} $ to get

∑ s = 0 ∞ ( ∂ t 2 p s + σ R ∂ t p s + g H p s μ s − g ζ eq , s ) ∇ ϕ s · ( ∇ ψ r × r ̂ ) + ∂ t p s ( f × ∇ ϕ s ) · ( ∇ ψ r × r ̂ ) + ∂ t p − s [ f × ( ∇ ψ s × r ̂ ) ] · ( ∇ ψ r × r ̂ ) + ( ∂ t 2 p − s + σ R ∂ t p − s ) ( ∇ ψ s × r ̂ ) · ( ∇ ψ r × r ̂ ) = 0 . $$ \begin{aligned} \nonumber&\sum _{s=0}^\infty \left(\partial _t^2p_s+\sigma _{\rm R}\partial _t p_s +gHp_s\mu _s -g{\zeta }_{\mathrm{eq} ,s}\right)\boldsymbol{\nabla }\phi _s\cdot \left(\boldsymbol{\nabla }\psi _r\boldsymbol{\times }\hat{r}\right)\\\nonumber&+ \partial _t p_s\left(\boldsymbol{f}\boldsymbol{\times }\boldsymbol{\nabla }\phi _s\right)\cdot \left(\boldsymbol{\nabla }\psi _r\boldsymbol{\times }\hat{r}\right)+\partial _t p_{-s}\left[\boldsymbol{f}\boldsymbol{\times }\left(\boldsymbol{\nabla }\psi _s\boldsymbol{\times }\hat{r}\right)\right]\cdot \left(\boldsymbol{\nabla }\psi _r\boldsymbol{\times }\hat{r}\right) \\&+\left(\partial _t^2p_{-s}+\sigma _{\rm R}\partial _tp_{-s}\right)\left(\boldsymbol{\nabla }\psi _s\boldsymbol{\times }\hat{r}\right)\cdot \left(\boldsymbol{\nabla }\psi _r\boldsymbol{\times }\hat{r}\right)=0. \end{aligned} $$(E.27)

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

∑ s = 0 ∞ ( ∂ t 2 p s + σ R ∂ t p s ) ∫ ∇ ψ s · ∇ ψ r d A − ∂ t p s ∫ ( f · r ̂ ) ( ∇ ϕ s · ∇ ψ r ) d A − ∂ t p − s ∫ ( f · r ̂ ) r ̂ · ( ∇ ψ r × ∇ ψ s ) d A = 0 , $$ \begin{aligned} \nonumber \sum _{s=0}^\infty&\left(\partial _t^2p_s+\sigma _{\rm R}\partial _t p_s \right)\int \boldsymbol{\nabla }\psi _s\cdot \boldsymbol{\nabla }\psi _r dA \\\nonumber&-\partial _t p_s\int \left(\boldsymbol{f}\cdot \hat{r}\right)\left(\boldsymbol{\nabla }\phi _s\cdot \boldsymbol{\nabla }\psi _r\right)dA \\&- \partial _t p_{-s}\int \left(\boldsymbol{f}\cdot \hat{r}\right)\hat{r}\cdot \left(\boldsymbol{\nabla }\psi _r\boldsymbol{\times }\boldsymbol{\nabla }\psi _s\right)dA =0 , \end{aligned} $$(E.28)

where upon replacing the Coriolis term as before we have

∂ t 2 p − r + σ R ∂ t p − r − 2 Ω ν r ∑ s = 1 ∞ ∂ t p s ∫ cos θ ( ∇ ϕ s · ∇ ψ r ) d A − 2 Ω ν r ∑ s = 0 ∞ ∂ t p − s ∫ cos θ r ̂ · ( ∇ ψ r × ∇ ψ s ) d A = 0 . $$ \begin{aligned} \nonumber&\partial _t^2p_{-r}+\sigma _{\rm R}\partial _t p_{-r}- \frac{2\Omega }{\nu _r}\sum _{s=1}^\infty \partial _t p_s \int \cos \theta \left(\boldsymbol{\nabla }\phi _s\cdot \boldsymbol{\nabla }\psi _{r}\right)dA\\&-\frac{2\Omega }{\nu _r}\sum _{s=0}^\infty \partial _t p_{-s}\int \cos \theta \hat{r}\cdot \left(\boldsymbol{\nabla }\psi _r\boldsymbol{\times }\boldsymbol{\nabla }\psi _s\right)dA=0. \end{aligned} $$(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

β r , s = − ∫ O cos θ r ̂ · ( ∇ ϕ r × ∇ ϕ s ) d A , β r , − s = ∫ O cos θ ∇ ϕ r · ∇ ψ s d A , β − r , s = − ∫ O cos θ ∇ ψ r · ∇ ϕ s d A , β − r , − s = − ∫ O cos θ r ̂ · ( ∇ ψ r × ∇ ψ s ) d A , $$ \begin{aligned} \nonumber \beta _{r,s}&=-\int _\mathcal{O} \cos \theta \, \hat{r} \cdot \left( \boldsymbol{\nabla }\phi _r\boldsymbol{\times }\boldsymbol{\nabla }\phi _s \right) dA,\\ \nonumber \beta _{r,-s}&= \ \ \, \int _\mathcal{O} \cos \theta \, \boldsymbol{\nabla }\phi _r\cdot \boldsymbol{\nabla }\psi _s dA, \\ \nonumber \beta _{-r,s}&=-\int _\mathcal{O} \cos \theta \, \boldsymbol{\nabla }\psi _r\cdot \boldsymbol{\nabla }\phi _s dA,\\ \beta _{-r,-s}&=-\int _\mathcal{O} \cos \theta \, \hat{r} \cdot \left( \boldsymbol{\nabla }\psi _r\boldsymbol{\times }\boldsymbol{\nabla }\psi _s \right) dA, \end{aligned} $$(E.30)

with β−s, r = −βr, −s. These coefficients carry the effect of rotational distortion to the tidal waves and the boundary conditions imposed by the coastlines. Using these definitions, Eqs. (E.26) and (E.29) form an infinite linear system in the coefficients pr(t) and p−r(t) that is expressed as:

∂ t 2 p r + σ R ∂ t p r + g H μ r p r − g ζ eq , r + 2 Ω μ r ∑ s = − ∞ s = ∞ β r , s ∂ t p s = 0 , $$ \begin{aligned} \partial _t^2p_r+\sigma _{\rm R}\partial _t p_r +gH\mu _rp_r -g{\zeta }_{\mathrm{eq} ,r} + \frac{2\Omega }{\mu _r}\sum _{s=-\infty }^{s=\infty }\beta _{r,s}\partial _tp_s =0, \end{aligned} $$(E.31)

∂ t 2 p − r + σ R ∂ t p − r + 2 Ω ν r ∑ s = − ∞ s = ∞ β − r , s ∂ t p s = 0 . $$ \begin{aligned} \partial _t^2p_{-r} +\sigma _{\rm R}\partial _tp_{-r} +\frac{2\Omega }{\nu _r}\sum _{s=-\infty }^{s=\infty }\beta _{-r,s}\partial _t p_s =0. \end{aligned} $$(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, ζ becomes a function of two moving surfaces: the free oceanic surface, ζos, and the vertically deforming solid surfacem ζ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 ( r , r ′ ) = GM | r − r ′ | − GM r ′ 2 r cos θ , $$ \begin{aligned} {U}(\boldsymbol{r},\boldsymbol{r}^{\prime })= \frac{GM}{|\boldsymbol{r}-\boldsymbol{r}^{\prime }|} - \frac{GM}{r^{\prime \, 2}}r\cos \theta \, , \end{aligned} $$(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 ( θ , λ , r ′ ) = U ( R , θ , λ , r ′ ) − GM r ′ , $$ \begin{aligned} {U}^\mathrm{T}(\theta ,\lambda , \boldsymbol{r}^{\prime }) = {U}(R,\theta ,\lambda ,\boldsymbol{r}^{\prime }) -\frac{GM}{r^\prime }\, , \end{aligned} $$(F.2)

where a constant offset was removed as it does not contribute to the tidal force. In the frequency domain, the tidal potential UT is expanded spectrally in Fourier series and spatially in spherical harmonics, with complex coefficients U n m ; s $ U_n^{m;s} $, as (Kaula 2013; Auclair-Desrotour et al. 2019):

U T = ∑ n = 2 ∞ ∑ m = − n n ∑ s = − ∞ ∞ U n m ; s P n m ( cos θ ) exp i ( σ m s t + m λ ) , $$ \begin{aligned} U^\mathrm{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}^st+m\lambda )}, \end{aligned} $$(F.3)

where s is an integer and the tidal forcing frequency σ m s = m Ω − s n o r b $ \sigma_{m}^s=m\Omega -s n_{\rm orb} $, the frequency norb being the orbital mean motion of the tidal perturber. In the absence of obliquity, the nth harmonic of the tidal potential U n m ; s $ U_n^{m;s} $ is given by (Ogilvie 2014):

U n m ; s = GM a ( R a ) n A n , m , s ( e ) , $$ \begin{aligned} U_n^{m;s} = \frac{GM}{a}\left(\frac{R}{a}\right)^n A_{n,m,s}(e), \end{aligned} $$(F.4)

where a is the semi-major axis of this perturber and An, 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 σ 2 2 = 2 ( Ω − n o r b ) $ \sigma_{2}^2=2(\Omega - n_{\rm orb}) $. For this component, while neglecting the small orbital eccentricity of the Sun and the Moon, A 2 , 2 , 2 ( 0 ) = 3 / 5 $ A_{2,2,2}(0)= \sqrt{3/5} $. Hereafter, we use U n T $ U^{\rm T}_n $ 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 θ ) exp { i ( σ m s t + m λ ) } . $$ \begin{aligned} U^\mathrm{T}_n = U_n^{m;s}P_n^m(\cos \theta )\exp \left\{ i(\sigma ^s_m t+m\lambda )\right\} \,. \end{aligned} $$(F.5)

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

ζ ¯ os = h n T U n T g + ∑ l 3 ρ oc ( 2 l + 1 ) ρ se h l L ζ l , $$ \begin{aligned} \bar{\zeta }_\mathrm{os} =\frac{ h^\mathrm{T}_nU^\mathrm{T}_n}{g}+ \sum _l \frac{3\rho _{\rm oc}}{(2l+1)\rho _{\rm se}} h^\mathrm{L}_l\zeta _l, \end{aligned} $$(F.6)

where ρoc and ρse stand for the uniform oceanic and solid Earth densities, respectively. In this equation, the oceanic depth variation, ζ, is decomposed into spherical harmonics defined over the full sphere as:

ζ l ( θ , λ , t ) = ∑ m = − l l ζ l m ( t ) P l m ( cos θ ) exp ( i m λ ) . $$ \begin{aligned} \zeta _l(\theta ,\lambda ,t) = \sum _{m=-l}^l \zeta _l^m(t)P_l^m(\cos \theta )\exp (im\lambda )\,. \end{aligned} $$(F.7)

Although ζ given in (E.22) is only defined over the oceanic hemisphere, this decomposition over the whole sphere is required when applying the Love numbers. Using the orthogonality of spherical harmonics, and the fact that ζ(θ, λ, t)=0 over the continental hemisphere, Eq. (F.7) can also be expressed as:

ζ l ( θ , λ , t ) = 1 2 ∑ m = 0 l α lm 2 ∫ O ζ ( θ ′ , λ ′ , t ) P l m ( cos θ ) P l m ( cos θ ′ ) × cos ( m ( λ − λ ′ ) ) d Ω , $$ \begin{aligned} \nonumber \zeta _l(\theta ,\lambda ,t) =&\frac{1}{2}\sum _{m=0}^l \alpha _{lm}^2 \int _{\mathcal{O} } \zeta (\theta ^{\prime },\lambda ^{\prime },t) P_l^m(\cos \theta )P_l^m(\cos \theta ^{\prime }) \\ &\times \cos (m(\lambda -\lambda ^{\prime }))\,d\Omega \,, \end{aligned} $$(F.8)

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

ζ ¯ ss = ( 1 + k n T ) U n T g + ∑ l 3 ρ oc ( 2 l + 1 ) ρ se ( 1 + k l L ) ζ l . $$ \begin{aligned} \bar{\zeta }_\mathrm{ss} = (1+k_n^\mathrm{T})\frac{ U_n^\mathrm{T}}{g} + \sum _l \frac{3\rho _{\rm oc}}{(2l+1)\rho _{\rm se}} (1+k_l^\mathrm{L}) \zeta _l. \end{aligned} $$(F.9)

In equations (F.6) and (F.9), we used the tidal Love numbers k n T $ k_n^{\rm T} $ and h n T $ h_n^{\rm T} $, and the surface-loading Love numbers k n L $ k_n^{\rm L} $ and h n L $ h_n^{\rm 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 $ k_n^{\rm T},h_n^{\rm T},k_n^{\rm L} $, and h n L $ h_n^{\rm 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 } = 1 ( 1 + μ ∼ n ) { 3 2 ( n − 1 ) , 2 n + 1 2 ( n − 1 ) , − 1 , − 2 n + 1 3 } , $$ \begin{aligned} \left\{ k_n^\mathrm{T},h_n^\mathrm{T},k_n^\mathrm{L},h_n^\mathrm{L}\right\} = \frac{1}{(1+\tilde{\mu }_n)}\left\{ \frac{3}{2(n-1)},\frac{2n+1}{2(n-1)},-1, -\frac{2n+1}{3} \right\} , \end{aligned} $$(F.10)

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

μ ∼ n = 4 ( 2 n 2 + 4 n + 3 ) π R 4 3 n G M E 2 μ E 1 + ( i σ τ A ) − α A Γ ( 1 + α A ) + ( i σ τ M ) − 1 , $$ \begin{aligned} \tilde{\mu }_n = \frac{4(2n^2+4n+3)\pi R^4}{3nGM_{\rm E}^2} \frac{\mu _{\rm E}}{1+(i\sigma \tau _{\rm A})^{-\alpha _{\rm A}} \Gamma (1+\alpha _{\rm A})+(i\sigma \tau _{\rm M})^{-1}}, \end{aligned} $$(F.11)

where ME is the mass of the Earth, μE its average rigidity, and Γ is the gamma function (Abramowitz et al. 1988); αA is a dimensionless rheological exponent determined experimentally (Castelnau et al. 2008; Petit & Luzum 2010); τA is the anelastic Andrade timescale; and τM the Maxwell relaxation time defined as the ratio of viscosity to rigidity. For a volumetric average of the mantle’s shear modulus μE = 17.3 × 1010 Pa, and volumetric average of viscosity deduced from inversions of Lau et al. (2016b), we have τM = 685 yrs. The values αA = 0.25 and τA = 2.19 × 104 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 ζeq in the momentum equation (E.1a) by the difference ζ ¯ os − ζ ¯ ss $ \bar{\zeta}_{\mathrm{os}}-\bar{\zeta}_{\mathrm{ss}} $ of Eqs.(F.6) and (F.9) and we resolve it in the Fourier domain using the forcing tidal frequency, σ. The modified momentum conservation equations is now expressed as:

i σ u + σ R u + f × u = − g ∇ ( − γ n T ζ ¯ + ∑ l γ l L ζ l ) , $$ \begin{aligned} i\sigma \boldsymbol{u}+\sigma _{\rm R}\boldsymbol{u}+\boldsymbol{f}\boldsymbol{\times }\boldsymbol{u} = -g\boldsymbol{\nabla }\left(-\gamma _n^\mathrm{T}\bar{\zeta } + \sum _l \gamma _l^\mathrm{L}\zeta _l\right), \end{aligned} $$(F.12)

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

γ n T = 1 + k n T − h n T ; γ l L = 1 − 3 ρ oc ( 2 l + 1 ) ρ se ( 1 + k l L − h l L ) . $$ \begin{aligned} \gamma _n^\mathrm{T}= 1+k_n^\mathrm{T}-h^\mathrm{T}_n ; \gamma _l^\mathrm{L}=1-\frac{3\rho _{\rm oc}}{(2l+1)\rho _{\rm se}}(1+k_l^\mathrm{L} -h_l^\mathrm{L}). \end{aligned} $$(F.13)

Just like the Love numbers, γ n T $ \gamma_n^{\rm T} $ and γ l L $ \gamma_l^{\rm 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 pr and p−r. Multiplying the added contribution of loading and self-attraction effects by ∇ϕr and ∇ ψ r × r ̂ $ \boldsymbol{\nabla}\psi_r\boldsymbol{\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:

− σ 2 p r − i σ σ R p r + g H μ r ( 1 − γ r L 2 ) p r − g γ n T ζ ¯ r − 2 i σ Ω μ r ∑ s = − ∞ s = ∞ β r , s p s − 1 2 g H ∑ s ′ = 1 s ′ ≠ r ∞ μ s ′ F r s ′ p s ′ = 0 , $$ \begin{aligned} \nonumber -\sigma ^2p_r-i\sigma \sigma _{\rm R} p_r +gH\mu _r\left(1-\frac{\gamma _r^\mathrm{L}}{2} \right)p_r -g\gamma ^\mathrm{T}_n\bar{\zeta }_r \\ -\frac{2i\sigma \Omega }{\mu _r}\sum _{s=-\infty }^{s=\infty }\beta _{r,s}p_s -\frac{1}{2}gH \sum _{\begin{matrix} s^\prime =1\\ s^\prime \ne r \end{matrix}}^\infty \mu _{s^\prime }F_r^{s^\prime }p_{s^\prime }&=0,\end{aligned} $$(F.14)

− σ 2 p − r − i σ σ R p − r − 2 i σ Ω ν r ∑ s = − ∞ s = ∞ β − r , s ∂ t p s = 0 , $$ \begin{aligned} -\sigma ^2p_{-r} -i\sigma \sigma _{\rm R}p_{-r} -\frac{2i\sigma \Omega }{\nu _r}\sum _{s=-\infty }^{s=\infty }\beta _{-r,s}\partial _t p_s&=0, \end{aligned} $$(F.15)

where we define

F r s ′ = 4 α n , m α n ′ , m ′ ∑ p ∑ q γ p L q 2 α p , q 2 O p , q n , m O p , q n ′ , m ′ ( q 2 − m 2 ) ( q 2 − m ′ 2 ) , $$ \begin{aligned} F_r^{s^\prime } = 4 \alpha _{n,m}\alpha _{n^\prime ,m^\prime }\sum _p\sum _q\gamma _p^\mathrm{L} q^2\alpha _{p,q}^2\frac{\mathcal{O} _{p,q}^{n,m}\mathcal{O} _{p,q}^{n^\prime ,m^\prime }}{(q^2-m^2)(q^2-m^{\prime 2})} \,, \end{aligned} $$(F.16)

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

O n , m u , v = ∫ − 1 1 P n m ( μ ) P u v ( μ ) d μ , $$ \begin{aligned} \mathcal{O} _{n,m}^{u,v} = \int _{-1}^1 P_n^m(\mu )P_u^v(\mu )d\mu , \end{aligned} $$(F.17)

for which the method of computation is detailed in Appendix H. Coupled to the orbital dynamical evolution, the tidal frequency σ is determined at each point in time in the hemispherical phase of the model, then the system is truncated at rmax and solved numerically (see Appendix J). We re-write the linear system as:

( a ( 1 ) + a r ( 2 ) ) p r + a r ( 3 ) ∑ s = − ∞ s = ∞ β r , s p s + a ( 5 ) ∑ s ′ = 1 s ′ ≠ r ∞ μ s ′ F r s ′ p s ′ = c r , $$ \begin{aligned}&(a^{(1)} + a_r^{(2)})p_r + a_r^{(3)} \sum _{s=-\infty }^{s=\infty }\beta _{r,s}p_s+a^{(5)} \sum _{\begin{matrix} s^\prime =1\\ s^\prime \ne r \end{matrix}}^\infty \mu _{s^\prime }F_r^{s^\prime }p_{s^\prime } = c_r ,\end{aligned} $$(F.18)

a ( 1 ) p r + a r ( 4 ) ∑ s = − ∞ s = ∞ β r , s p s = 0 , $$ \begin{aligned}&a^{(1)} p_r + a_r^{(4)} \sum _{s=-\infty }^{s=\infty }\beta _{r,s}p_s = 0, \end{aligned} $$(F.19)

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

a ( 1 ) = − σ 2 − i σ σ R , a r ( 2 ) = g H μ r ( 1 − γ r L / 2 ) , a r ( 3 ) = − 2 i σ Ω μ r − 1 , a r ( 4 ) = − 2 i σ Ω ν − r − 1 , a ( 5 ) = − 1 2 g H , c r = g γ n T ζ r ¯ . $$ \begin{aligned} \nonumber a^{(1)}&= -\sigma ^2 -i\sigma \sigma _{\rm R}, a^{(2)}_r = gH\mu _r(1-\gamma _r^\mathrm{L}/2), \\ \nonumber a^{(3)}_r&= -2i\sigma \Omega \mu _r^{-1}, a^{(4)}_r = -2i\sigma \Omega \nu _{-r}^{-1},\\ a^{(5)}&= -\frac{1}{2}gH, c_r = g\gamma _n^\mathrm{T}\bar{\zeta _r}. \end{aligned} $$(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 ( μ ) = ( − 1 ) m 2 n n ! ( 1 − μ 2 ) m / 2 ∂ μ n + m ( μ 2 − 1 ) n , $$ \begin{aligned} P_n^m(\mu )= \frac{(-1)^m}{2^n n!}(1-\mu ^2)^{m/2}\partial _\mu ^{n+m}(\mu ^2-1)^n, \end{aligned} $$(G.1)

which are solutions to the Legendre equation,

∂ μ [ ( 1 − μ 2 ) ∂ μ P n m ] + [ n ( n + 1 ) − m 2 1 − μ 2 ] P n m = 0 . $$ \begin{aligned} \partial _\mu \left[(1-\mu ^2)\partial _\mu P_n^m\right]+\left[n(n+1)-\frac{m^2}{1-\mu ^2}\right]P_n^m=0. \end{aligned} $$(G.2)

Following differentiation, we obtain

∂ μ P n m = − m μ 1 − μ 2 P n m − P n m + 1 1 − μ 2 . $$ \begin{aligned} \partial _\mu P_n^m= -\frac{m\mu }{1-\mu ^2}P_n^m - \frac{P_n^{m+1}}{\sqrt{1-\mu ^2}}. \end{aligned} $$(G.3)

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

P n m + 2 − 2 m μ 2 ( m + 1 ) 1 − μ 2 P n m − 2 μ ( m + 1 ) ∂ μ P n m + ( n ( n + 1 ) − m ( m + 1 ) ) P n m = 0 , $$ \begin{aligned} \nonumber&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, \end{aligned} $$(G.4)

which gives the useful relation

μ ∂ μ P n m = P n m + 2 2 ( m + 1 ) + [ n ( n + 1 ) + m ( m + 1 ) 2 ( m + 1 ) − m 1 − μ 2 ] P n m . $$ \begin{aligned} \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 . \end{aligned} $$(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,

μ P n m 1 − μ 2 = − 1 2 m ( P n m + 1 + ( n − m + 1 ) ( n + m ) P n m − 1 ) , $$ \begin{aligned} &\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), \end{aligned} $$(G.6)

P n m 1 − μ 2 = − 1 2 m ( P n − 1 m + 1 + ( n + m − 1 ) ( n + m ) P n − 1 m − 1 ) , $$ \begin{aligned}&\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), \end{aligned} $$(G.7)

1 − μ 2 ∂ μ P n m = − 1 2 P n m + 1 + 1 2 ( n + m ) ( n − m + 1 ) P n m − 1 , $$ \begin{aligned}&\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},\end{aligned} $$(G.8)

( 1 − μ 2 ) ∂ μ P n m = 1 2 n + 1 ( ( n + 1 ) ( n + m ) P n − 1 m − n ( n − m + 1 ) P n + 1 m ) . $$ \begin{aligned}&(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)\!. \end{aligned} $$(G.9)

thumbnail Fig. G.1.

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

The theory of the hemispherical tidal response is based on an ocean bound by two meridians. Thus for an oceanic center moving on the sphere, we instead rotate the spin axis relative to the center of the ocean, and accordingly the frame of the tidal perturber to maintain the coplanar configuration of the dynamical system. These rotations will enter the system through the Coriolis term, specifically through the gyroscopic coefficients, along with the tidal forcing term. We define an arbitrary rotation {θ, λ}→{θ′,λ′} using an Eulerian rotation matrix of the form ℛ3(α)ℛ2(β)ℛ3(γ), with (0 ≤ α ≤ 2π) and (0 ≤ β ≤ π), and we fix γ = 0 (see Fig. G.1). For a vector J defined as:

J = R 3 ( α ) R 2 ( β ) ( 0 0 1 ) T = ( sin β cos α sin β sin α cos β ) T , $$ \begin{aligned} \nonumber J&= \mathcal{R} _3(\alpha )\mathcal{R} _2(\beta )\begin{pmatrix} 0&0&1 \end{pmatrix}^T \\&=\begin{pmatrix} \sin \beta \cos \alpha&\sin \beta \sin \alpha&\cos \beta \end{pmatrix}^T, \end{aligned} $$(G.10)

the transformed gyroscopic coefficients are:

R 2 α r α s β r , s = J z β r , s ( 1 ) + J x β r , s ( 2 ) + J y β r , s ( 3 ) , $$ \begin{aligned}&\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)}, \end{aligned} $$(G.11)

R 2 α r α s β r , − s = J z β r , − s ( 1 ) + J x β r , − s ( 2 ) + J y β r , − s ( 3 ) , $$ \begin{aligned}&\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)},\end{aligned} $$(G.12)

R 2 α r α s β − r , − s = J z β − r , − s ( 1 ) + J x β − r , − s ( 2 ) + J y β − r , − s ( 3 ) , $$ \begin{aligned}&\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)},\end{aligned} $$(G.13)

β − r , s = − β s , − r , $$ \begin{aligned}&\beta _{-r,s}= -\beta _{s,-r}, \end{aligned} $$(G.14)

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

β r , s ( 1 ) = 2 m 2 − v 2 ∫ [ v 2 P u v ∂ μ P n m + m 2 P n m ∂ μ P u v ] μ d μ , $$ \begin{aligned} \beta _{r,s}^{(1)}&= \frac{2}{m^2-v^2} \int \left[ v^2 P_u^v\partial _\mu P_n^m+ m^2P_n^m\partial _\mu P_u^v \right]\mu d\mu , \end{aligned} $$(G.15)

β r , s ( 2 ) = π 4 ∫ [ m K m , v ( 1 ) P n m ∂ μ P u v μ ¯ 1 / 2 − v K m , v ( 2 ) ∂ μ P n m P u v μ ¯ 1 / 2 ] d μ , $$ \begin{aligned} \beta _{r,s}^{(2)}&= \frac{\pi }{4} \int \left[ 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}\right]d\mu ,\end{aligned} $$(G.16)

β r , s ( 3 ) = ∫ [ m K m , v ( 3 ) P n m ∂ μ P u v μ ¯ 1 / 2 − v K m , v ( 4 ) ∂ μ P n m P u v μ ¯ 1 / 2 ] d μ , $$ \begin{aligned} \beta _{r,s}^{(3)}&= \int \left[ 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}\right]d\mu ,\end{aligned} $$(G.17)

β r , − s ( 1 ) = − 2 v m 2 − v 2 ∫ [ ∂ μ P u v ∂ μ P n m μ μ ¯ + m 2 P n m P u v μ μ ¯ ] d μ , $$ \begin{aligned} \beta _{r,-s}^{(1)}&= \frac{-2v}{m^2-v^2} \int \left[ \partial _\mu P_u^v\partial _\mu P_n^m \mu \bar{\mu }+ m^2 P_n^m P_u^v \frac{\mu }{\bar{\mu }} \right]d\mu ,\end{aligned} $$(G.18)

β r , − s ( 2 ) = π 4 ∫ [ K m , v ( 2 ) ∂ μ P n m ∂ μ P u v μ ¯ 3 / 2 − m v K m , v ( 1 ) P n m P u v μ ¯ − 1 / 2 ] d μ , $$ \begin{aligned} \beta _{r,-s}^{(2)}&= \frac{\pi }{4} \int \left[ K_{m,v}^{(2)} \partial _\mu P_n^m \partial _\mu P_u^v \bar{\mu }^{3/2} -mv K_{m,v}^{(1)} P_n^m P_u^v \bar{\mu }^{-1/2}\right]d\mu , \end{aligned} $$(G.19)

β r , − s ( 3 ) = ∫ [ K m , v ( 4 ) ∂ μ P n m ∂ μ P u v μ ¯ 3 / 2 − m v K m , v ( 3 ) P n m P u v μ ¯ − 1 / 2 ] d μ , $$ \begin{aligned} \beta _{r,-s}^{(3)}&= \int \left[ K_{m,v}^{(4)} \partial _\mu P_n^m \partial _\mu P_u^v \bar{\mu }^{3/2} -mv K_{m,v}^{(3)} P_n^m P_u^v \bar{\mu }^{-1/2}\right]d\mu , \end{aligned} $$(G.20)

β − r , − s ( 1 ) = 2 m v m 2 − v 2 ∫ [ P u v ∂ μ P n m + P n m ∂ μ P u v ] μ d μ , $$ \begin{aligned} \beta _{-r,-s}^{(1)}&= \frac{2mv}{m^2-v^2} \int \left[ P_u^v\partial _\mu P_n^m+ P_n^m\partial _\mu P_u^v \right]\mu d\mu , \end{aligned} $$(G.21)

β − r , − s ( 2 ) = π 4 ∫ [ v K m , v ( 1 ) ∂ μ P n m P u v μ ¯ 1 / 2 − m K m , v ( 2 ) P n m ∂ μ P u v μ ¯ 1 / 2 ] d μ , $$ \begin{aligned} \beta _{-r,-s}^{(2)}&= \frac{\pi }{4} \int \left[ 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}\right]d\mu , \end{aligned} $$(G.22)

β − r , − s ( 3 ) = ∫ [ v K m , v ( 3 ) ∂ μ P n m P u v μ ¯ 1 / 2 − m K m , v ( 4 ) P n m ∂ μ P u v μ ¯ 1 / 2 ] d μ , $$ \begin{aligned} \beta _{-r,-s}^{(3)}&= \int \left[ 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}\right]d\mu , \end{aligned} $$(G.23)

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

K m , v ( 1 ) = ( 1 + δ v , 0 ) δ m − v , 1 − δ m − v , − 1 , $$ \begin{aligned} K_{m,v}^{(1)}&= (1+\delta _{v,0})\delta _{m-v,1} -\delta _{m-v,-1}, \end{aligned} $$(G.24)

K m , v ( 2 ) = K v , m ( 1 ) , $$ \begin{aligned} K_{m,v}^{(2)}&= K_{v,m}^{(1)}, \end{aligned} $$(G.25)

K m , v ( 3 ) = m ( 1 m 2 − ( v 2 + 1 ) 2 + 1 m 2 − ( v 2 − 1 ) 2 ) , $$ \begin{aligned} K_{m,v}^{(3)}&= m \left( \frac{1}{m^2-(v^2+1)^2} + \frac{1}{m^2-(v^2-1)^2}\right),\end{aligned} $$(G.26)

K m , v ( 4 ) = K v , m ( 3 ) . $$ \begin{aligned} K_{m,v}^{(4)}&= K_{v,m}^{(3)}. \end{aligned} $$(G.27)

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

cos θ ′ = cos θ cos β + sin θ sin β cos ( λ − α ) . $$ \begin{aligned} \cos \theta ^\prime = \cos \theta \cos \beta + \sin \theta \sin \beta \cos (\lambda -\alpha ). \end{aligned} $$(G.28)

To compute the integrals involved in the gyroscopic coefficients, we make use of the essential condition1 (e.g., Longuet-Higgins & Pond (1970))

P n m P u v | μ = ± 1 = 0 , $$ \begin{aligned} \left. P_n^m P_u^v \right|_{\mu =\pm 1} =0, \end{aligned} $$(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 β r , s ( 1 ) $ \beta_{r,s}^{(1)} $ and β r , s ( 2 ) $ \beta_{r,s}^{(2)} $ are expressed as:

β r , s ( 1 ) = [ u ( u + 1 ) + v ( v + 1 ) v + 1 − 2 v n ( n + 1 ) − u ( u + 1 ) + v m 2 − v 2 ] O n , m u , v + 1 v + 1 O n , m u , v + 2 , $$ \begin{aligned} \nonumber \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]\mathcal{O} _{n,m}^{u,v} \\&+ \frac{1}{v+1}\mathcal{O} _{n,m}^{u,v+2}, \end{aligned} $$(G.30)

β r , s ( 2 ) = π 4 [ m K m , v ( 1 ) ∫ P n m ∂ μ P u v μ ¯ 1 / 2 d μ − v K m , v ( 2 ) ∫ ∂ μ P n m P u v μ ¯ 1 / 2 d μ ] , $$ \begin{aligned} \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] ,\end{aligned} $$(G.31a)

= π 8 { m K m , v ( 1 ) [ ( u + v ) ( u − v + 1 ) O n , m u , v − 1 − O n , m u , v + 1 ] − v K m , v ( 2 ) [ ( n + m ) ( n − m + 1 ) O n , m − 1 u , v − O n , m + 1 u , v ] } , $$ \begin{aligned} =&\frac{\pi }{8}\left\{ m K_{m,v}^{(1)}\left[(u+v)(u-v+1)\mathcal{O} _{n,m}^{u,v-1} - \mathcal{O} _{n,m}^{u,v+1}\right]\right.\nonumber \\&\qquad \left. -vK_{m,v}^{(2)}\left[(n+m)(n-m+1)\mathcal{O} _{n,m-1}^{u,v} - \mathcal{O} _{n,m+1}^{u,v}\right] \right\} , \end{aligned} $$(G.31b)

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

β r , s ( 3 ) = 1 2 { m K m , v ( 3 ) [ ( u + v ) ( u − v + 1 ) O n , m u , v − 1 − O n , m u , v + 1 ] − v K m , v ( 4 ) [ ( n + m ) ( n − m + 1 ) O n , m − 1 u , v − O n , m + 1 u , v ] } , $$ \begin{aligned} \beta _{r,s}^{(3)}= \frac{1}{2}&\left\{ m K_{m,v}^{(3)}\left[(u+v)(u-v+1)\mathcal{O} _{n,m}^{u,v-1} - \mathcal{O} _{n,m}^{u,v+1}\right]\right.\nonumber \\&\qquad \left. -vK_{m,v}^{(4)}\left[(n+m)(n-m+1)\mathcal{O} _{n,m-1}^{u,v} - \mathcal{O} _{n,m+1}^{u,v}\right] \right\} , \end{aligned} $$(G.32)

β r , − s ( 1 ) = − 2 v ( m 2 − v 2 ) ( 2 n + 1 ) { ( n + 1 ) ( n − 1 ) ( n + m ) O n − 1 , m u , v + n ( n + 2 ) ( n − m + 1 ) O n + 1 , m u , v } , $$ \begin{aligned} \beta _{r,-s}^{(1)}=&\frac{-2v}{(m^2-v^2)(2n+1)} \left\{ (n+1)(n-1)(n+m) \mathcal{O} _{n-1,m}^{u,v} \right.\nonumber \\&\qquad \left. + n(n+2)(n-m+1) \mathcal{O} _{n+1,m}^{u,v} \right\} , \end{aligned} $$(G.33)

β r , − s ( 2 ) = π 4 [ K m , v ( 2 ) ∫ ∂ μ P n m ∂ μ P u v μ ¯ 3 / 2 d μ − m v K m , v ( 1 ) ∫ P n m P u v μ ¯ − 1 / 2 d μ ] , $$ \begin{aligned} \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 -mv K_{m,v}^{(1)} \!\!\int \!\! P_n^m P_u^v \bar{\mu }^{-1/2}d\mu \right] ,\end{aligned} $$(G.34a)

= π 4 [ K m , v ( 2 ) ∫ ∂ μ P n m μ ¯ ∂ μ P u v μ ¯ 1 / 2 d μ − m v K m , v ( 1 ) ∫ P n m P u v μ ¯ − 1 / 2 d μ ] , $$ \begin{aligned}&= \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 -mv K_{m,v}^{(1)} \!\! \int \!\! P_n^m P_u^v \bar{\mu }^{-1/2}d\mu \right] ,\end{aligned} $$(G.34b)

= π K m , v ( 2 ) 8 ( 2 n + 1 ) { ( n + 1 ) ( n + m ) [ ( u + v ) ( u − v + 1 ) O n − 1 , m u , v − 1 − O n − 1 , m u , v + 1 ] + n ( n − m + 1 ) [ O n + 1 , m u , v + 1 − ( u + v ) ( u − v + 1 ) O n + 1 , m u , v − 1 ] } + π v K m , v ( 1 ) 8 { O n − 1 , m + 1 u , v + ( n + m − 1 ) ( n + m ) O n − 1 , m − 1 u , v } , $$ \begin{aligned} \nonumber&= \frac{\pi K_{m,v}^{(2)} }{8(2n+1)} \Bigg \{ \\\nonumber&\,\,\,\,\,\,(n+1)(n+m) \left[ (u+v)(u-v+1)\mathcal{O} _{n-1,m}^{u,v-1} - \mathcal{O} _{n-1,m}^{u,v+1}\right] \\\nonumber&\,\,\,+n(n-m+1) \left[\mathcal{O} _{n+1,m}^{u,v+1} - (u+v)(u-v+1)\mathcal{O} _{n+1,m}^{u,v-1} \right]\Bigg \} \\&\,\,\,+\frac{\pi v K_{m,v}^{(1)}}{8}\left\{ \mathcal{O} _{n-1,m+1}^{u,v} + (n+m-1)(n+m)\mathcal{O} _{n-1,m-1}^{u,v}\right\} , \end{aligned} $$(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 β r , − s ( 3 ) $ \beta_{r,-s}^{(3)} $, β − r , − s ( 1 ) $ \beta_{-r,-s}^{(1)} $, β − r , − s ( 2 ) $ \beta_{-r,-s}^{(2)} $, and β − r , − s ( 3 ) $ \beta_{-r,-s}^{(3)} $ are expressed as:

β r , − s ( 3 ) = K m , v ( 4 ) 2 ( 2 n + 1 ) { ( n + 1 ) ( n + m ) [ ( u + v ) ( u − v + 1 ) O n − 1 , m u , v − 1 − O n − 1 , m u , v + 1 ] + n ( n − m + 1 ) [ O n + 1 , m u , v + 1 − ( u + v ) ( u − v + 1 ) O n + 1 , m u , v − 1 ] } + v K m , v ( 3 ) 2 { O n − 1 , m + 1 u , v + ( n + m − 1 ) ( n + m ) O n − 1 , m − 1 u , v } , $$ \begin{aligned} \nonumber \beta _{r,-s}^{(3)}=&\frac{ K_{m,v}^{(4)} }{2(2n+1)} \Bigg \{ \\\nonumber&\,\,\,(n+1)(n+m) \left[ (u+v)(u-v+1)\mathcal{O} _{n-1,m}^{u,v-1} - \mathcal{O} _{n-1,m}^{u,v+1}\right] \\\nonumber&+n(n-m+1) \left[\mathcal{O} _{n+1,m}^{u,v+1} - (u+v)(u-v+1)\mathcal{O} _{n+1,m}^{u,v-1} \right]\Bigg \} \\&+\frac{ v K_{m,v}^{(3)}}{2}\left\{ \mathcal{O} _{n-1,m+1}^{u,v} + (n+m-1)(n+m)\mathcal{O} _{n-1,m-1}^{u,v}\right\} , \end{aligned} $$(G.35)

β − r , − s ( 1 ) = − 2 m v m 2 − v 2 O n , m u , v , $$ \begin{aligned} \beta _{-r,-s}^{(1)}= \frac{-2mv}{m^2-v^2} \mathcal{O} _{n,m}^{u,v}, \end{aligned} $$(G.36)

β − r , − s ( 2 ) = π 8 { − m K m , v ( 2 ) [ ( u + v ) ( u − v + 1 ) O n , m u , v − 1 − O n , m u , v + 1 ] + v K m , v ( 1 ) [ ( n + m ) ( n − m + 1 ) O n , m − 1 u , v − O n , m + 1 u , v ] } , $$ \begin{aligned} \beta _{-r,-s}^{(2)}=&\frac{\pi }{8}\left\{ -m K_{m,v}^{(2)}\left[(u+v)(u-v+1)\mathcal{O} _{n,m}^{u,v-1} - \mathcal{O} _{n,m}^{u,v+1}\right]\right.\nonumber \\&\qquad \left. +vK_{m,v}^{(1)}\left[(n+m)(n-m+1)\mathcal{O} _{n,m-1}^{u,v} - \mathcal{O} _{n,m+1}^{u,v}\right] \right\} , \end{aligned} $$(G.37)

β − r , − s ( 3 ) = 1 2 { − m K m , v ( 4 ) [ ( u + v ) ( u − v + 1 ) O n , m u , v − 1 − O n , m u , v + 1 ] + v K m , v ( 3 ) [ ( n + m ) ( n − m + 1 ) O n , m − 1 u , v − O n , m + 1 u , v ] } . $$ \begin{aligned} \beta _{-r,-s}^{(3)}= \frac{1}{2}&\left\{ -m K_{m,v}^{(4)}\left[(u+v)(u-v+1)\mathcal{O} _{n,m}^{u,v-1} - \mathcal{O} _{n,m}^{u,v+1}\right]\right.\nonumber \\&\qquad \left. +vK_{m,v}^{(3)}\left[(n+m)(n-m+1)\mathcal{O} _{n,m-1}^{u,v} - \mathcal{O} _{n,m+1}^{u,v}\right] \right\} . \end{aligned} $$(G.38)

Appendix H: Overlap integral O n , m u , v $ \mathcal{O}_{n,m}^{u,v} $

Here, we provide a closed form solution for the computation of the overlap integral of Eq. (F.17). The procedure is assimilated from tools of angular momentum quantization (Varshalovich et al. 1988). Following Dong & Lemus (2002) and introducing the notation q = v − m, we have:

O n , m u , v = C n , m u , v ∑ l ( 2 l + 1 ) D ( | q | , l ) · ( n u l 0 0 0 ) ( n u l − m v m − v ) , $$ \begin{aligned} \mathcal{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}, \end{aligned} $$(H.1)

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

C n , m u , v = ( − 1 ) κ 2 | q | − 2 | q | ( n + m ) ! ( u + v ) ! ( n − m ) ! ( u − v ) ! , $$ \begin{aligned} C_{n,m}^{u,v} = (-1)^\kappa 2^{|q|-2}|q| \sqrt{\frac{(n+m)!(u+v)!}{(n-m)!(u-v)!}}, \end{aligned} $$(H.2)

and the coefficients 𝒟(|q|,l) by

D ( | q | , l ) = [ 1 + ( − 1 ) l + | q | ] ( l − | q | ) ! ( l + | q | ) ! Γ ( l / 2 ) Γ ( ( l + | q | + 1 ) / 2 ) ( ( l − | q | ) / 2 ) ! Γ ( ( l + 3 ) / 2 ) . $$ \begin{aligned} \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)}. \end{aligned} $$(H.3)

We note here that the phase κ introduced in Dong & Lemus (2002) as

κ = { m if v ≥ m , v otherwise , $$ \begin{aligned} \kappa = {\left\{ \begin{array}{ll} m&\text{ if}\ \ v \ge m,\\ v&\text{ otherwise}, \end{array}\right.} \end{aligned} $$(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

( a b c d e f ) = ( − 1 ) R 21 + R 31 + R 32 [ R 31 ! R 32 ! R 33 ! R 13 ! R 23 ! ( J + 1 ) ! R 11 ! R 12 ! R 21 ! R 22 ! ] 1 / 2 × ∑ z ( − 1 ) z ( R 21 + z ) ! ( R 11 + R 31 − z ) ! z ! ( R 31 − z ) ! ( R 23 − z ) ! ( R 13 − R 31 + z ) ! , $$ \begin{aligned} \nonumber \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)!}, \end{aligned} $$(H.5)

where J = a + b + c, and Rij are the elements of the so-called Regge ℜ-symbol (Regge 1958) defined as

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 . $$ \begin{aligned} \nonumber&R_{11}=-a+b+c, & R_{12}= a-b+c, &R_{13}= a+b-c, \\ \nonumber&R_{21}= a+d , &R_{22}= b+e,&R_{23}= c +f,\\&R_{31}= a-d , &R_{32}= b-e,&R_{33}=c-f. \end{aligned} $$(H.6)

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 $ \mathcal{O}_{n,m}^{u,v}=0 $ when v = m. In that case, we alternatively use

O n , m u , m = 2 2 n + 1 ( n + m ) ! ( n − m ) ! δ n , u . $$ \begin{aligned} \mathcal{O} _{n,m}^{u,m} = \frac{2}{2n+1}\frac{(n+m)!}{(n-m)!}\delta _{n,u}\,. \end{aligned} $$(H.7)

This method for the computation of the overlap integral was verified numerically using MATLAB’s ALFs package.

Appendix I: The tidal forcing term ζ r ¯ $ \bar{\zeta_r} $

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

Y p q ( θ , λ ) = 2 p + 1 4 π ( p − q ) ! ( p + q ) ! P p q ( cos θ ) exp ( i q λ ) , $$ \begin{aligned} Y_p^q(\theta ,\lambda ) = \sqrt{\frac{2p+1}{4\pi }\frac{(p-q)!}{(p+q)!}}\,P_p^q(\cos \theta )\,\exp (iq\lambda )\,, \end{aligned} $$(I.1)

with angular frequency σ, we have:

ζ ¯ = 2 π Y p q ( θ , λ ) exp ( i σ t ) . $$ \begin{aligned} \bar{\zeta }= \sqrt{2\pi }Y_{p}^{q}(\theta ,\lambda )\exp (i\sigma t). \end{aligned} $$(I.2)

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

Y p s ( θ ′ , λ ′ ) = ∑ q = − p p Y p q ( θ , λ ) D s , q p ( α , β , γ ) , $$ \begin{aligned} Y_p^{s} (\theta ^\prime ,\lambda ^\prime )= \sum _{q=-p}^p Y_p^q(\theta ,\lambda )D_{s,q}^p(\alpha ,\beta ,\gamma ), \end{aligned} $$(I.3)

or

Y p q ( θ , λ ) = ∑ s = − p p Y p s ( θ ′ , λ ′ ) D s , q p ∗ ( α , β , γ ) , $$ \begin{aligned} Y_p^{q} (\theta ,\lambda )= \sum _{s=-p}^p Y_p^{s}(\theta ^\prime ,\lambda ^\prime )D_{s,q}^{p*}(\alpha ,\beta ,\gamma ), \end{aligned} $$(I.4)

where D q , s p $ 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: α, β, or γ,

D s , q p ( α , β , γ ) = e − i q α d sq p ( β ) e − i s γ . $$ \begin{aligned} D_{s,q}^{p}(\alpha ,\beta ,\gamma )= e^{-iq\alpha }d_{sq}^{p}(\beta )e^{-is\gamma }. \end{aligned} $$(I.5)

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

d sq p ( β ) = ( − 1 ) p − s [ ( p + q ) ! ( p − q ) ! ( p + s ) ! ( p − s ) ! ] 1 / 2 × ∑ j ( − 1 ) j ( cos β / 2 ) q + s + 2 j ( sin β / 2 ) 2 p − q − s − 2 j j ! ( p − q − j ) ! ( p − s − j ) ! ( q + s + j ) ! , $$ \begin{aligned} \nonumber 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)!}{,} \end{aligned} $$(I.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 ( θ ′ , λ ′ ) $ Y_p^{s}(\theta^\prime,\lambda^\prime) $ 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,

ζ ¯ = π R 2 2 exp ( i σ t ) ∑ s = − p p D s , q p ∗ ( α , β , γ ) ( 1 + δ s , 0 ) 1 / 2 [ ϕ p s + i ψ p s ] . $$ \begin{aligned} \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} \left[\phi _p^{s} + i\psi _p^{s}\right]{.} \end{aligned} $$(I.7)

Then, invoking the definition of the component ζ r ¯ $ \bar{\zeta_r} $,

ζ r ¯ = ∫ O ϕ r ζ ¯ d A , $$ \begin{aligned} \bar{\zeta _r} = \int _\mathcal{O} \phi _r\bar{\zeta } dA {,} \end{aligned} $$(I.8)

we get its expression in the rotated frame of reference,

ζ r ¯ = π R 2 2 exp ( i σ t ) ∑ s = − p p D s , q p ∗ ( α , β , γ ) ( 1 + δ s , 0 ) 1 / 2 × [ ∫ O ϕ p s ϕ n m d A + i ∫ O ψ p s ϕ n m d A ] , $$ \begin{aligned} \nonumber \bar{\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 \Bigg [\int _\mathcal{O} \phi _p^s\phi _n^mdA + i\int _\mathcal{O} \psi _p^{s}\phi _n^mdA\Bigg ] {,} \end{aligned} $$(I.9)

where the dot products of the eigenfunctions are simplified as:

∫ O ϕ p s ϕ n m d A = { δ n , p , if s = m , ( − 1 ) m δ n , p , if s = − m , 0 , otherwise , $$ \begin{aligned} \int _\mathcal{O} \phi _p^s\phi _n^mdA = {\left\{ \begin{array}{ll} \delta _{n,p} ,&\text{ if}\ \ s=m, \\ (-1)^m\delta _{n,p} ,&\text{ if}\ \ s=-m,\\ 0 ,&\text{otherwise}, \end{array}\right.} \end{aligned} $$(I.10)

and

∫ O ψ p s ϕ n m d A = { 0 , if m + s = even , α p , s α n , m 2 s s 2 − m 2 O p , s n , m otherwise. $$ \begin{aligned} \int _\mathcal{O} \psi _p^s\phi _n^mdA ={\left\{ \begin{array}{ll} 0 ,&\text{ if}\ \ m+s ={\text{ even}},\\ \displaystyle \alpha _{p,s}\,\alpha _{n,m} \frac{2s}{s^2-m^2}\mathcal{O} _{p,s}^{n,m}&\text{otherwise.} \end{array}\right.} \end{aligned} $$(I.11)

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

P p − s = ( − 1 ) s ( p − s ) ! ( p + s ) ! P p s $$ \begin{aligned} P_p^{-s}= (-1)^s \frac{(p-s)!}{(p+s)!}P_p^s \end{aligned} $$(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 (pr, 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 𝒯, the power lost by the Earth would be 𝒯Ω, and the power gained by the perturber is 𝒯norb. The difference between them is the dissipative work of the total tidal mass redistribution 𝒲diss, thus

T = W diss Ω − n orb . $$ \begin{aligned} \mathcal{T} = \frac{ \mathcal{W} _{\rm diss}}{\Omega -n_{\rm orb}}. \end{aligned} $$(J.1)

The total dissipative work is the sum of two contributions: the dissipative work of oceanic tidal currents, W diss oc $ \mathcal{W}_{\mathrm{diss}}^{\mathrm{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, W diss oc $ \mathcal{W}_{\mathrm{diss}}^{\mathrm{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 contributions is only jeopardized by the emerging significance of the solid dissipation when moving backwards in time and increasing the Earth’s rotational velocity Ω. Solid Earth dissipation would also be amplified with an early less viscous mantle due to higher 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.

thumbnail Fig. J.1.

Numerical analysis on the dependence of the tidal response computation on the truncation order rmax. The response is quantified by the root mean square tidal amplitude, ζrms, (Eq. J.5) and the dissipative work, W diss oc $ \mathcal{W}_{\mathrm{diss}}^{\mathrm{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.

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

W diss oc = 〈 ∫ O u ( t ) · σ R u ( t ) d A 〉 = 1 2 σ R σ 2 ∑ r = 1 ∞ ( μ r p r p r ∗ + ν r p − r p − r ∗ ) , $$ \begin{aligned} \nonumber \mathcal{W} _{\rm diss}^\mathrm{oc}&=\left\langle \int _\mathcal{O} \boldsymbol{u}(t)\cdot \sigma _{\rm R}\boldsymbol{u}(t) dA\right\rangle \\&=\frac{1}{2}\sigma _{\rm R}\sigma ^2 \sum _{r=1}^{\infty }\left(\mu _rp_rp_r^*+\nu _rp_{-r}p_{-r}^*\right), \end{aligned} $$(J.2)

where ⟨⟩ denotes time averaging over the tidal period. This work should be equal to the work done by the tidal force on the ocean:

W tide oc = 〈 ρ oc g H ∫ O ∇ ζ ¯ ( t ) · u ( t ) d A 〉 = 1 2 ρ oc g H σ I { ∑ r = 1 ∞ p r μ r ζ ¯ r ∗ } . $$ \begin{aligned} \nonumber \mathcal{W} _{\rm tide}^\mathrm{oc}&=\left\langle \rho _{\rm oc} gH\int _\mathcal{O} \boldsymbol{\nabla }\bar{\zeta }(t)\cdot \boldsymbol{u}(t)dA\right\rangle \\&=\frac{1}{2}\rho _{\rm oc} gH\sigma \mathfrak{I} \left\{ \sum _{r=1}^\infty p_r\mu _r \bar{\zeta }_r^* \right\} . \end{aligned} $$(J.3)

Hence, the tidal torque associated with the lunar semi-diurnal frequency σ = 2(Ω − nM), nM being the lunar mean motion, is

T M = ρ oc g H I { ∑ r = 1 ∞ p r μ r ζ ¯ r ∗ } , $$ \begin{aligned} \mathcal{T} _{\rm M}= \rho _{\rm oc} gH \mathfrak{I} \left\{ \sum _{r=1}^\infty p_r\mu _r \bar{\zeta }_r^* \right\} , \end{aligned} $$(J.4)

and we obtain a similar expression for the solar tides 𝒯S when solving the system with the solar tidal frequency component σ = 2(Ω − nS), nS being the solar mean motion that generates the solar tidal work.

Besides the tidal torque, the tidal response can also be quantified by the root mean square tidal height variation, ζrms, given as

ζ rms = H π R 2 ∑ r = 1 ∞ μ r 2 p r ∗ p r . $$ \begin{aligned} \zeta _{\rm rms} = \sqrt{\frac{H}{\pi R^2} \sum _{r=1}^\infty \mu _r^2p_r^*p_r}. \end{aligned} $$(J.5)

As these quantities are computed numerically, a truncation order, rmax, is required. In Fig. J.1, we show the numerical dependence of the tidal response on rmax. Since the response is dominated by gravity modes, the tidal solution converges fast enough with rmax. To avoid any truncation effect in our computation, and to properly account for the resonances, we adopted rmax = 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 text3. 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:

u = ∑ m , σ u m , σ ( θ ) exp i ( σ t + m λ , ζ = ∑ m , σ ζ m , σ ( θ ) exp i ( σ t + m λ , ζ eq = ∑ m , σ ζ eq m , σ ( θ ) exp i ( σ t + m λ . $$ \begin{aligned} \nonumber&\boldsymbol{u}= \sum _{m,\sigma }\boldsymbol{u}^{m,\sigma }(\theta ) \exp {i(\sigma t+m\lambda }, \\\nonumber&\zeta = \sum _{m,\sigma }\zeta ^{m,\sigma }(\theta ) \exp {i(\sigma t+m\lambda },\\&\zeta _{\rm eq}= \sum _{m,\sigma }\zeta _{\rm eq}^{m,\sigma }(\theta ) \exp {i(\sigma t+m\lambda }. \end{aligned} $$(K.1)

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

σ ∼ = σ − i σ R and ν ∼ = 2 Ω σ ∼ , $$ \begin{aligned} \tilde{\sigma }=\sigma -i\sigma _{\rm R} \text{ and} \tilde{\nu }=\frac{2\Omega }{\tilde{\sigma }}, \end{aligned} $$(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) { Θ n m , ν ∼ ( θ ) $ \Theta_n^{m,\tilde{\nu}}(\theta) $}, associated with a set of eigenvalues { Λ n m , ν ∼ $ \Lambda_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

Θ n m , ν ∼ ( θ ) = ∑ m ≤ l A n , l m , ν ∼ P l m ( cos θ ) , P l m ( cos θ ) = ∑ n B l , n m , ν ∼ Θ n m , ν ∼ ( θ ) , $$ \begin{aligned} \nonumber \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} $$(K.3)

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

ζ l m , σ = ∑ n A n , l m , ν ∼ ζ n m , σ , $$ \begin{aligned} \zeta _l^{m,\sigma } = \sum _n A_{n,l}^{m,\tilde{\nu }} \zeta _n^{m,\sigma }, \end{aligned} $$(K.4)

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

( σ σ ∼ I N − [ σ 1 , 1 2 ⋯ σ 1 , n 2 ⋯ σ 1 , N 2 ⋮ ⋱ ⋮ ⋮ σ n , 1 2 ⋯ σ n , n 2 ⋯ σ n , N 2 ⋮ ⋮ ⋱ ⋮ σ N , 1 2 ⋯ σ N , n 2 ⋯ σ N , N 2 ] ) [ ζ 1 m , σ ⋮ ζ n m , σ ⋮ ζ N m , σ ] = [ F 1 m , σ ⋮ F n m , σ ⋮ F N m , σ ] . $$ \begin{aligned} \left( \sigma \tilde{\sigma } \boldsymbol{I}_N - \begin{bmatrix} \sigma ^2_{1,1}&\dots&\sigma ^2_{1,n}&\dots&\sigma ^2_{1,N} \\ \vdots&\ddots&\vdots&\vdots \\ \sigma ^2_{n,1}&\dots&\sigma ^2_{n,n}&\dots&\sigma ^2_{n,N} \\ \vdots&\vdots&\ddots&\vdots \\ \sigma ^2_{N,1}&\dots&\sigma ^2_{N,n}&\dots&\sigma ^2_{N,N} \end{bmatrix} \right) \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}{.} \end{aligned} $$(K.5)

thumbnail 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 λ = 0 and λ = π (in blue). Energy dissipation is quantified by the linear Rayleigh drag frequency σR. The logarithm of the torque is plotted as a function of the normalized frequency ω = (Ω − norb)/Ω0, where the Earth’s spin rate varies with the tidal forcing frequency Ω = norb + σ/2 at fixed norb, and Ω0 being the present spin rate of the Earth.

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

In this linear system, IN denotes the identity matrix of size N × N, the forcing terms of the studied tidal potential U l m , σ $ U_l^{m,\sigma} $ (Eq. F.3) are expressed as

F n m , σ = − H Λ n m , v ∼ R 2 ∑ m ≤ l B l , n m , ν ∼ γ l T U l m , σ , $$ \begin{aligned} \mathcal{F} _n^{m,\sigma } = -\frac{H \Lambda _n^{m,\tilde{v}}}{R^2} \sum _{m\le l}B_{l,n}^{m,\tilde{\nu }}\gamma _l^T U_l^{m,\sigma }, \end{aligned} $$(K.6)

and the complex characteristic frequencies σn, k as

σ n , k = g H k ̂ n 2 ∑ l ≥ m γ l L A k ; l m , ν ∼ B l , n m , ν ∼ , $$ \begin{aligned} \sigma _{n,k}= \sqrt{gH \hat{k}_n^2\sum _{l\ge m}\gamma _{l}^L A_{k;l}^{m,\tilde{\nu }} B_{l,n}^{m,\tilde{\nu }}}, \end{aligned} $$(K.7)

where the horizontal wave-number of the degree-n mode k ̂ n = Λ n m , ν ∼ / R $ \hat{k}_n= \sqrt{\Lambda_n^{m,\tilde{\nu}}}/R $, and the coupling coefficients γ l T $ \gamma_l^{\rm T} $ and γ l L $ \gamma_l^{\rm L} $ 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

k l m , σ = k l T + ( 1 + k l L ) 3 g 2 l + 1 ρ oc ρ se ζ l m , σ U l m , σ . $$ \begin{aligned} \mathfrak{k} _l^{m,\sigma } = k_l^\mathrm{T} + \left(1+k_l^\mathrm{L}\right)\frac{3g}{2l+1}\frac{\rho _{\rm oc}}{\rho _{\rm se}}\frac{\zeta _l^{m,\sigma }}{U_l^{m,\sigma }}. \end{aligned} $$(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, k l T $ k_l^{\rm 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 , σ $ 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)

T l m = 3 2 G M 2 R 5 a 6 I k l m , σ . $$ \begin{aligned} \mathcal{T} _l^m = \frac{3}{2}GM^2 \frac{R^5}{a^6}\mathfrak{I} {\mathfrak{k} _l^{m,\sigma }}. \end{aligned} $$(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 and different orders of magnitude of σR. We consider the semi-diurnal lunar gravitational forcing exerted on the Earth. The spectrum of the torque is plotted against the normalized frequency ω = (Ω − norb)/Ω0, where Ω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)

σ n ± = i σ R 2 ± g H k ̂ n 2 − ( σ R 2 ) 2 , $$ \begin{aligned} \sigma _n^\pm = i\frac{\sigma _{\rm R}}{2} \pm \sqrt{gH \hat{k}_n^2 - \left(\frac{\sigma _{\rm R}}{2}\right)^2} {,} \end{aligned} $$(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 σR. In the opposite limit, resonant peaks are spelled out intensifying in amplitude as friction is weakened. Besides, when σR → 0, the frequencies σ n ± $ \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 proportionally 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 hemispheric response against the fairly regular resonance periodicity in the global configuration. Several resonances can be encountered in the hemispherical configuration spectrum in between two resonant peaks of the global configuration. Applied to the 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 σ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).

All Tables

Table C.1.

Misfit analysis summary.

In the text
Table D.1.

Cyclostratigraphic data.

In the text
Table D.2.

Tidal rhythmites data.

In the text
Table D.3.

Values of constant parameters used in the numerical implementation of the theory.

In the text

All Figures

thumbnail Fig. 1.

Temporal evolution of the latitude of the surface “paleo-barycenter” 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 the text
thumbnail Fig. 2.

Misfit surfaces of χ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 σR. The misfit is established using the currently measured lunar recession rate via LLR, and the lunar age (Appendix C). The three models differ in the imposed geometry of the oceanic shell over the geological history, with the combined model featuring more physical realism that the other two. The numerical results of this analysis are summarized in (Table C.1). The dynamical evolution associated with each of the misfit minima is plotted in terms of: the lunar semi-major axis in Fig. 3, length of the day in Fig. 5, and obliquity and precession frequency in Fig. 6.
In the text
thumbnail Fig. 3.

Evolution of the lunar semi-major axis over time. The Earth-Moon separation, aM, is plotted for the three studied models, taking the best-fit values of the free parameters (H, σR) as described in Fig. 2 and in the main text. Plotted on top of the evolution curves: Geological inferences of aM 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 aE extends to 3RE, but the y-axis is trimmed to start at 15RE for a better visualization of the geological data.
In the text
thumbnail Fig. 4.

History of the tidal torque. The logarithm of the semi-diurnal tidal torque of the Earth (normalized by its present value: T ∼ = T / T ( t = 0 ) $ \tilde{\mathcal{T}}= \mathcal{T}/\mathcal{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 tswitch, 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.
In the text
thumbnail Fig. 5.

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

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

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

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

Numerical analysis on the dependence of the tidal response computation on the truncation order rmax. The response is quantified by the root mean square tidal amplitude, ζrms, (Eq. J.5) and the dissipative work, W diss oc $ \mathcal{W}_{\mathrm{diss}}^{\mathrm{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.
In the text
thumbnail 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 λ = 0 and λ = π (in blue). Energy dissipation is quantified by the linear Rayleigh drag frequency σR. The logarithm of the torque is plotted as a function of the normalized frequency ω = (Ω − norb)/Ω0, where the Earth’s spin rate varies with the tidal forcing frequency Ω = norb + σ/2 at fixed norb, and Ω0 being the present spin rate of the Earth.
In the text
thumbnail 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.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.