# Can one hear supercontinents in the tides of ocean planets?

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

► To cite this version:

Pierre Auclair-Desrotour, Mohammad Farhat, Gwenaël Boué, Mickaël Gastineau, Jacques Laskar. Can one hear supercontinents in the tides of ocean planets?. Astronomy & Astrophysics - A&A, 2023, 680, 10.1051/0004-6361/202347301. insu-04849148

HAL Id: insu-04849148

https://insu.hal.science/insu-04849148v1

Submitted on 19 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Can one hear supercontinents in the tides of ocean planets?

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

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

Received 27 June 2023 / Accepted 9 October 2023

## ABSTRACT

**Context.** Recent observations and theoretical progress made about the history of the Earth-Moon system suggest that tidal dissipation in oceans primarily drives the long-term evolution of orbital systems hosting ocean planets. Particularly, they emphasise the key role played by the geometry of land-ocean distributions in this mechanism. However, the complex way continents affect oceanic tides still remains to be elucidated.

**Aims.** In the present study, we investigate the impact of a single supercontinent on the tidal response of an ocean planet and the induced tidally dissipated energy.

**Methods.** The adopted approach is based on the linear tidal theory. By simplifying the continent to a spherical cap of a given angular radius and position on the globe, we carried out a harmonic analysis of the whole planet's tidal response including the coupling with the solid part due to ocean loading and self-attraction variations. In this framework, tidal flows are formulated analytically in terms of explicitly defined oceanic eigenmodes, as well as the resulting tidal Love numbers, dissipated power, and torque.

**Results.** The analysis highlights the symmetry breaking effect of the continent, which makes the dependence of tidal quantities on the tidal frequency become highly irregular. The metric introduced to quantify this continentality effect reveals abrupt transitions between polar and non-polar configurations, and between small-sized and medium-sized continents. Additionally, it predicts that a continent similar to South America or smaller (~30° angular radius) does not qualitatively alter the tidal response of a global ocean regardless of its position on the planet.

**Key words.** hydrodynamics – planet–star interactions – planets and satellites: oceans – planets and satellites: terrestrial planets – Earth

## 1. Introduction

'Can one hear the shape of a drum?' This question, raised by Kac (1966), emphasises the link connecting the geometry of a vibrating system to the frequencies at which it can vibrate. Mathematically, it defines a broad class of problems referred to as spectral geometry, which aims to establish relationships between the eigenmodes of Riemannian manifolds and their geometric features. In the aforementioned study, Kac investigates whether it is possible or not to infer some information about the shape of a vibrating drumhead from the sound it makes, which amounts to recovering the drumhead's geometry from its acoustic signature in an inverse problem approach. Unfortunately, the answer to this question was shown to be negative in the general case, as two different shapes are able to produce the very same sound (e.g. Gordon et al. 1992). Nevertheless, it remains possible to infer some specific geometric features of the system to a certain extent.

Interestingly, this statement is also applicable to tidally forced ocean planets, which may be regarded as giant celestial vibrating drumheads. In the general case, one calls ‘ocean planets’ rocky bodies that are partly or totally covered by a surface liquid layer<fnref n="1" /> (e.g. Léger et al. 2004). With more than 70% of its surface covered by oceans (e.g. Eakins & Sharman 2010), the Earth thus appears as a typical ocean planet. It is actually the

only planet of this kind that can be found presently in the Solar System, although observational evidences indicate that Mars possibly had large paleo-oceans as well before it dried (Carr & Head 2003; Dohm et al. 2009; Scheller et al. 2021).

This uniqueness of the Earth, however, diminishes if one considers extrasolar systems. The remarkable exoplanetary diversity established thus far suggests that ocean planets are rather ordinary, if not widespread. Particularly, a substantial number of rocky planets were found to orbit in the habitable zone of their host stars, which is the region where temperature conditions can sustain liquid water on the planets' surfaces (Kasting et al. 1993). This is the case of many rocky planets orbiting red and brown dwarfs (e.g. Payne & Lodato 2007; Raymond et al. 2007; Kopparapu et al. 2017) such as the super-Earths GJ 1214 b (e.g. Charbonneau et al. 2009), LHS 1140 b (e.g. Lillo-Box et al. 2020), and TOI-1452 b (Cadieux et al. 2022), or the Earth-sized planets TRAPPIST-1 e, f, and g (e.g. Grimm et al. 2018), with the latter being suspected to harbour ocean-scale volumes of liquid water (e.g. Bolmont et al. 2017; Bourrier et al. 2017; Turbet et al. 2018). Improving our understanding of the surface conditions, climate, and fate of these planets requires to constrain their long-term orbital evolution better (e.g. Pierrehumbert 2010).

Tides are the main mechanism that controls the evolution of planetary systems over timescales ranging from millions to billions of years (e.g. Hut 1981; Correia & Laskar 2001; Levrard et al. 2009). Apart from the atmospheric thermal tides that are induced by insolation variations (e.g. Lindzen & Chapman 1969; Leconte et al. 2015; Auclair-Desrotour et al. 2019b), they result from the mutual gravitational interactions between

<footnote n="1">¹ The extended definition of ocean planets – or ocean worlds – may also refer to bodies containing a substantial amount of water in the form of subsurface oceans, such as Jupiter’s moon Europa (e.g. Kivelson et al. 2000).</footnote>
celestial bodies, which undergo the differential attraction of their neighbours. Combined with dissipative processes, this tidal forcing generates a delayed mass redistribution, thereby leading to exchanges of angular momentum between the orbit of the tidal perturber and the spin of the tidally forced body (e.g. Correia et al. 2014). Additionally, tides are accompanied with energy dissipation where the mechanical energy lost by the orbital system is converted into heat. While it is negligible on Earth, this tidal heating is sometimes able to partly melt the body and to generate surface volcanism, as observed on the Jovian moon Io (e.g. Peale et al. 1979).

Oceanic tides take the form of frequency-resonant surface gravity modes distorted by the planet's rotation (e.g. Auclair-Desrotour et al. 2018), which strongly differs from the smooth visco-elastic elongation of solid bodies (e.g. Correia et al. 2014). As a consequence, tidal dissipation in oceans may be increased by several orders of magnitude while crossing a resonance associated with a predominant mode (e.g. Arbic & Garrett 2010). This resonant behaviour notably explains why, for present day Earth, the oceanic contribution (~2.5 TW) to the semidiurnal component ($M_2$) of dissipated tidal energy is an order of magnitude greater than the solid counterpart (Lambeck 1977; Provost & Lyard 1997; Egbert & Ray 2001; Tyler 2021). Analogously, oceanic tidal flows presumably warm up the icy moons of the Solar System outer planets that are suspected to harbour subsurface oceans, such as Europa, Callisto, or Titan (e.g. Tyler 2008, 2014; Kamata et al. 2015; Beuthe 2016; Matsuyama et al. 2018).

Starting with MacDonald (1964), several authors investigated the role played by oceanic tides in the evolution of the Earth-Moon system, using either semi-analytical approaches with simplified geometries (e.g. Webb 1980, 1982; Bills & Ray 1999; Farhat et al. 2022a) – including global ocean models (e.g. Auclair-Desrotour et al. 2018; Motoyama et al. 2020; Tyler 2021) –, or numerical methods based upon realistic land-ocean distributions (e.g. Le Provost et al. 1994, 1998; Arbic et al. 2010; Kodaira et al. 2016; Green et al. 2017; Blackledge et al. 2020; Daher et al. 2021). For a didactical review of the pioneering studies of the field, the reader is referred to Lambeck (1977) and Bills & Ray (1999). This long series of works particularly emphasises the strong interplay between the resonant oceanic tidal flows and continents, which shapes the four billion year-history of the Earth's length of the day (LOD) and Earth-Moon distance (see e.g. Daher et al. 2021; Farhat et al. 2022a). Moreover, Green et al. (2017) show that the topography significantly increases tidal dissipation, and Blackledge et al. (2020) underline the sensitivity of the latter to coastlines' fractality. However, the way the large-scale geometry of the ocean basin alters the planet's tidal response cannot be easily disentangled from other effects in realistic models due to its complexity.

The present work is an attempt to address this issue by generalising the semi-analytical tidal theory of hemispherical oceans (Longuet-Higgins & Pond 1970; Webb 1980, 1982; Farhat et al. 2022a) to ocean basins of arbitrary sizes. In this approach, the oceanic depth is assumed to be uniform and thin compared to the planet radius. Adopting the so-called 'shallow water' approximation (e.g. Vallis 2006), we ignore both the ocean stratification and the associated baroclinic component of tidal flows – that is the contribution of internal gravity waves (e.g. Gerkema & Zimmerman 2008) – so that the described tide is purely barotropic. Also, the geometry of the continent is simplified to a spherical cap of specified angular radius in order to avoid mathematical complications. Finally, bottom friction is modelled by a standard Rayleigh drag, and the coupling effect between the ocean and the deformable solid surface is taken into account.

Analogously with the examples of the global and hemispherical oceans, this formalism allows the planet's tidal response and the resulting dissipated energy to be formulated in terms of explicit eigenmodes, each of these eigenmodes being resonant for a specific tidal frequency. Additionally, the ocean dynamics is controlled by a small number of dimensionless parameters, which provides an appropriate framework for probing the parameter space and characterising the continentality effect.

In Sect. 2, we detail the theory by successively introducing the Laplace's Tidal Equations (LTEs) that govern the oceanic tidal response, the eigenmodes of the ocean basin, and the expressions of the tidal Love numbers, torque, and power. In Sect. 3, we examine the frequency-behaviour of an Earth-sized planet with an equatorial hemispherical ocean. This reference case is used in Sect. 4 to characterise the sensitivity of the planet's tidal response to the position of the continent on the globe, its size, the oceanic depth, the dissipation timescale, the elasticity of the solid part, and its anelasticity timescale. We elaborate on the continentality effect in Sect. 5 by introducing a metric that quantifies the dependence of the tidal torque on the position and size of the supercontinent. Finally, in Sect. 6, the conclusions of the study are summarised. We stress here that Sect. 2 details technical aspects of the developed model. Thus we invite the reader that might be only interested in the application of the theory to skip this section and to jump directly to Sect. 3. It is also noteworthy that all the notations introduced in the main text can be retrieved in the nomenclature made in Appendix A.

## 2. Oceanic tidal response

We establish here the equations governing the tidal dynamics of the ocean basin harboured by a rocky planet of radius $R_p$. These equations are based on the commonly used shallow water approach, where the ocean is considered as a thin liquid layer of uniform density and depth $H \ll R_p$ (e.g. Vallis 2006). In this approach, the fluid is supposed to be incompressible. The tidal response thus described is said to be barotropic because it does not depend on the vertical structure of the ocean (Vallis 2006). The continental geometry is simplified to a spherical cap of angular radius $\theta_c$, as shown by Fig. 1. Accordingly, the angular radius of the ocean basin is defined as $\theta_0 \equiv \pi - \theta_c$.

### 2.1. Geometry of the ocean basin

To write down the equations that describe the oceanic tidal waves, we shall preliminarily introduce two frames of references and their associated systems of coordinates. First, we denote by $\mathcal{R}: (O, \mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$ the frame of reference rotating with the planet and having its centre of gravity, O, as origin. The Cartesian basis of unit vectors $(\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$ is such that $\mathbf{e}_x$ and $\mathbf{e}_y$ correspond to two orthogonal directions in the planet's equatorial plane, and $\mathbf{e}_z$ to the direction of the spin vector. The unit vector $\mathbf{e}_z$ thus designates the position of the north pole on the unit sphere. The geocentric frame of reference $\mathcal{R}$ is associated with the usual spherical coordinates $(\hat{r}, \hat{\theta}, \hat{\phi})$, where $\hat{r}, \hat{\theta}$, and $\hat{\phi}$ designate the radial, colatitudinal and longitudinal coordinates, respectively. The position of the circular supercontinent on the globe is defined by the coordinates of its centre, $(\hat{\theta}_s, \hat{\phi}_s)$. As the centre of the oceanic basin corresponds to the antipodal point, its coordinates on the globe are given by $(\hat{\theta}_{oc}, \hat{\phi}_{oc}) = (\pi - \hat{\theta}_s, \pi + \hat{\phi}_s)$.

The oceanic frame of reference, $\mathcal{R}_{oc}: (O, \mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$, is such that $\mathbf{e}_z$ points towards the centre of the ocean basin, while $\mathbf{e}_x$ and
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geometry of the studied system and associated parameters. The unit vector $\mathbf{e}'_z$, which indicates the position of the continental centre, is defined as $\mathbf{e}'_z \equiv -\mathbf{e}_z$, with $\mathbf{e}_z$ pointing towards the centre of the ocean basin. In the diagram, the longitude of the continental centre is set to $\hat{\phi}_S = 0$.</figcaption></figure>

$\mathbf{e}_y$ define two orthogonal axes in the plane normal to $\mathbf{e}_z$. This frame of reference is associated with the spherical coordinates $(r, \theta, \varphi)$, where $r = \hat{r}$ is the radial coordinate, $\theta$ the colatitude ($\theta = 0$ at the centre of the ocean basin), and $\varphi$ the longitude. The change of basis vectors $(\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z) \rightarrow (\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$ is expressed as a function of the coordinates of the continental centre as

$$ \begin{bmatrix} \mathbf{e}_x \\ \mathbf{e}_y \\ \mathbf{e}_z \end{bmatrix} = \begin{bmatrix} \cos \hat{\phi}_S \cos \hat{\theta}_S & \sin \hat{\phi}_S \cos \hat{\theta}_S & -\sin \hat{\theta}_S \\ \sin \hat{\phi}_S & -\cos \hat{\phi}_S & 0 \\ -\cos \hat{\phi}_S \sin \hat{\theta}_S & -\sin \hat{\phi}_S \sin \hat{\theta}_S & -\cos \hat{\theta}_S \end{bmatrix} \begin{bmatrix} \mathbf{e}_X \\ \mathbf{e}_Y \\ \mathbf{e}_Z \end{bmatrix}. \quad (1) $$

We note that the longitude of the continental centre does not affect the tidal response since the tidal forcing is periodic in longitude. Therefore, we set this longitude to $\hat{\phi}_S = 0$ in the following, similar to the example shown by Fig. 1, which yields

$$ \begin{bmatrix} \mathbf{e}_x \\ \mathbf{e}_y \\ \mathbf{e}_z \end{bmatrix} = \begin{bmatrix} \cos \hat{\theta}_S & 0 & -\sin \hat{\theta}_S \\ 0 & -1 & 0 \\ -\sin \hat{\theta}_S & 0 & -\cos \hat{\theta}_S \end{bmatrix} \begin{bmatrix} \mathbf{e}_X \\ \mathbf{e}_Y \\ \mathbf{e}_Z \end{bmatrix}. \quad (2) $$

Finally, we introduce the set of unit vectors $(\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_\varphi)$ associated with the spherical coordinates $(r, \theta, \varphi)$, respectively. The change of basis vectors $(\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z) \rightarrow (\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_\varphi)$ is given by the standard relation,

$$ \begin{bmatrix} \mathbf{e}_r \\ \mathbf{e}_\theta \\ \mathbf{e}_\varphi \end{bmatrix} = \begin{bmatrix} \sin \theta \cos \varphi & \sin \theta \sin \varphi & \cos \theta \\ \cos \theta \cos \varphi & \cos \theta \cos \varphi & -\sin \theta \\ -\sin \varphi & \cos \varphi & 0 \end{bmatrix} \begin{bmatrix} \mathbf{e}_x \\ \mathbf{e}_y \\ \mathbf{e}_z \end{bmatrix}. \quad (3) $$

## 2.2. Laplace's tidal equations

The planet is subject to the effects of the tidal gravitational potential generated by the perturber (the Moon or the Sun in the Earth case), which is expressed, in the frame of reference $\mathcal{R}$, as

$$ U_T \equiv \frac{GM_s}{|R_p \mathbf{e}_r - \mathbf{r}_s|} - \frac{GM_s}{r_s^2} R_p \cos \hat{\theta} - \frac{GM_s}{r_s}, \quad (4) $$

where the constant component and the component responsible for the Keplerian dynamics of the two-body system are removed (second and third terms in the right-hand member of Eq. (4)). In this equation, $G$ is the universal gravitational constant, $M_s$ the mass of the perturber, $\mathbf{r}_s$ its position vector, and $r_s \equiv |\mathbf{r}_s|$ the planet-perturber distance. We note that the symbol $\equiv$ is used throughout the text to distinguish between definitions and equalities. The tidal force per unit mass exerted by the perturber on the planet is given by $F \equiv \nabla U_T$. Following the formalism introduced in earlier works (e.g. Tyler 2011; Matsuyama 2014; Auclair-Desrotour et al. 2018, 2019a; Motoyama et al. 2020; Farhat et al. 2022a), we write down the momentum and mass conservation equations, respectively, as

$$ \partial_t V + \sigma_R V + f \times V + g \nabla (\Gamma_D \zeta - \Gamma_T \zeta_{eq}) = 0, \quad (5) $$

$$ \partial_t \zeta + \nabla \cdot (HV) = 0, \quad (6) $$

with $t$ designating the time, $g$ the surface gravity at rest, $\sigma_R$ the Rayleigh drag frequency used to describe the action of dissipative mechanisms, $V$ the horizontal velocity vector – defined from the horizontal displacement vector $\xi$ as $V \equiv \partial_t \xi$, $f$ the Coriolis parameter, $\zeta$ the vertical displacement of the ocean’s surface with respect to the oceanic floor, and $\zeta_{eq} \equiv U_T/g$ the equilibrium displacement corresponding to the equipotential surface induced by the tidal gravitational potential.

In the momentum equation, the notations $\Gamma_D$ and $\Gamma_T$ refer to non-trivial linear operators accounting for the effects of ocean loading, self-attraction, and deformation of the solid regions of the planet (e.g. Hendershott 1972). These operators are called the ‘solid deformation operators’ in the following. They encompass the complex coupling between the oceanic shell and the solid interior described by Poisson’s equation, and the momentum and rheological equations governing the tidal dynamics of the solid part. The two operators simplify to $\Gamma_D = \Gamma_T = 1$ if one neglects both the tidal deformation of the solid part (infinite-rigidity approximation) and the variation of self-attraction induced by the oceanic tidal response (Cowling approximation; see Cowling 1941; Unno et al. 1989). In the general case, they are determined by the rheological behaviour of the solid part in its response to gravitational and surface forcings.

Equations (5) and (6) are known as the Laplace’s tidal equations (hereafter, LTEs) in reference to Laplace’s masterpiece (Laplace 1798). We note that $\zeta_{eq}$ results from the coupled oceanic tidal response and tidal deformation of the solid part, which includes both gravitational and loading interactions, as discussed further. The Coriolis parameter is expressed as a function of the colatitude of the current point in $\mathcal{R}$,

$$ f \equiv 2\Omega \cos \hat{\theta} \mathbf{e}_r, \quad (7) $$

where $\Omega > 0$ is the planet’s spin angular velocity. The horizontal gradient operator and divergence of the horizontal velocity are expressed in $\mathcal{R}_{oc}$, and the associated basis vectors, $(\mathbf{e}_\theta, \mathbf{e}_\varphi)$, as

$$ \nabla \equiv R_p^{-1} [\mathbf{e}_\theta \partial_\theta + \mathbf{e}_\varphi (\sin \theta)^{-1} \partial_\varphi], \quad (8) $$

$$ \nabla \cdot V = (R_p \sin \theta)^{-1} [\partial_\theta (\sin \theta V_\theta) + \partial_\varphi V_\varphi]. \quad (9) $$

The Rayleigh drag term $\sigma_R V$ in Eq. (5) accounts for the cumulated effects of dissipative mechanisms on tidal flows. The frequency $\sigma_R$ is the inverse of the effective dissipation timescale associated with these dissipative mechanisms. It takes values
<table><caption>**Table 1.** Dimensionless control parameters determining the regime of the oceanic tidal response.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Limits</th><th>Asymptotic regimes</th><th>Reference</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">Geometry of the ocean basin</td></tr><tr><td>θ<sub>c</sub></td><td>Size of the supercontinent</td><td>θ<sub>c</sub> = 0°<br>θ<sub>c</sub> = 180°</td><td>Global ocean<br>Dry planet</td><td>Fig. 1</td></tr><tr><td>θ&#770;<sub>S</sub></td><td>Position of the supercontinent on the globe</td><td>θ&#770;<sub>S</sub> = 0°<br>θ&#770;<sub>S</sub> = 90°</td><td>Polar continent<br>Equatorial continent</td><td>Eq. (1)</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Ocean dynamics</td></tr><tr><td>σ&#770; ≡ σ/2Ω</td><td>Distortion of forced tidal waves by Coriolis forces</td><td>σ&#770; ≪ 1<br>σ&#770; ≫ 1</td><td>Sub-inertial regime<br>Super-inertial regime</td><td>Eq. (15)</td></tr><tr><td>σ&#770;<sub>G</sub> ≡ √gH/2ΩR<sub>p</sub></td><td>Deviation of free surface waves by Coriolis forces</td><td>σ&#770;<sub>G</sub> ≪ 1<br>σ&#770;<sub>G</sub> ≫ 1</td><td>Fast rotator regime<br>Slow rotator regime</td><td>Eq. (14)</td></tr><tr><td>σ&#770;<sub>R</sub> ≡ σ<sub>R</sub>/2Ω</td><td>Damping of Coriolis effects by friction</td><td>σ&#770;<sub>R</sub> ≪ 1<br>σ&#770;<sub>R</sub> ≫ 1</td><td>Quasi-adiabatic regime<br>Frictional regime</td><td>Eq. (14)</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Solid deformation</td></tr><tr><td>γ<sup>σ</sup><sub>D;l</sub>, γ<sup>σ</sup><sub>T;l</sub></td><td>Visco-elastic response of the solid part</td><td>γ<sup>σ</sup><sub>D;l</sub>, γ<sup>σ</sup><sub>T;l</sub> = 1<br>else</td><td>Rigid body (with CA)<br>Solid-ocean coupling</td><td>Eq. (47)</td></tr></tbody></table>

**Notes.** The acronym 'CA' refers to the Cowling approximation (Cowling 1941), mentioned in Sect. 2.2.

around ~$10^{-5}$ s⁻¹ for the present day Earth (e.g. Webb 1980; Farhat et al. 2022a)<fnref n="2" />.

## 2.3. Nondimensional tidal equations

The nondimensional momentum and continuity equations are obtained by choosing as reference time and velocity scales the inertial time, $t_0$, and the typical velocity of long-wavelength surface gravity waves, $V_0$, which are defined, respectively, as

$$t_0 \equiv (2\Omega)^{-1}, \qquad V_0 \equiv \sqrt{gH}. \qquad (10)$$

Introducing the normalised time $\tilde{t}$, Coriolis parameter $\tilde{f}$, horizontal gradient $\tilde{\nabla}$, the complex horizontal displacement vector $\tilde{\xi}$, and vertical displacements $\tilde{\zeta}$ and $\tilde{\zeta}_{eq}$, such that

$$\begin{gathered} t = t_0 \tilde{t}, \qquad f = t_0^{-1} \tilde{f} \qquad \nabla = R_p^{-1} \tilde{\nabla}, \\ \xi = \mathfrak{R}(R_p \tilde{\xi}), \qquad \zeta = \mathfrak{R}(H \tilde{\zeta}), \qquad \zeta_{eq} = \mathfrak{R}(H \tilde{\zeta}_{eq}), \end{gathered} \qquad (11)$$

with $\mathfrak{R}$ referring to the real part of a complex number, we end up with the nondimensional complex LTEs given by

$$[\partial_{\tilde{\pi}} + (\tilde{\sigma}_R + \tilde{f} \times) \partial_{\tilde{t}}] \tilde{\xi} + \tilde{\sigma}_G^2 \tilde{\nabla} \tilde{\Xi} = 0, \qquad (12)$$

$$\tilde{\zeta} + \tilde{\nabla} \cdot \tilde{\xi} = 0, \qquad (13)$$

with $\tilde{\Xi} = \Gamma_D \tilde{\zeta} - \Gamma_T \tilde{\zeta}_{eq}$. In the above equations, tidal dynamics are controlled by two dimensionless parameters,

$$\tilde{\sigma}_G \equiv \frac{V_0}{2\Omega R_p} = \frac{\sqrt{gH}}{2\Omega R_p}, \qquad \tilde{\sigma}_R \equiv \frac{\sigma_R}{2\Omega}. \qquad (14)$$

The first parameter, $\tilde{\sigma}_G$, can be considered as a normalised Rossby deformation length since it compares the typical propagation velocity of surface gravity waves ($V_0$) with the Earth's rotation velocity. If $\tilde{\sigma}_G \ll 1$ (fast rotator regime), the inertial forces resulting from Coriolis acceleration predominate with respect to the restoring forces of surface gravity waves (pressure forces and gravity). Conversely, if $\tilde{\sigma}_G \gg 1$ (slow rotator regime), Coriolis terms are not strong enough to significantly deviate free surface gravity waves. As it describes the ratio of drag forces to Coriolis forces, the second dimensionless parameter, $\tilde{\sigma}_R$, may be regarded as an Ekman number. If $\tilde{\sigma}_R \ll 1$, the drag does not alter the tidal response much. Conversely, if $\tilde{\sigma}_R \gg 1$ characterises a frictional (or viscous) regime where inertial effects are annihilated by the strong damping associated with the drag.

In addition to $\tilde{\sigma}_G$ and $\tilde{\sigma}_R$, the time-derivative operators in Eq. (12) introduce a third dimensionless parameter describing the ratio of tidal forces to Coriolis forces,

$$\tilde{\sigma} \equiv \frac{\sigma}{2\Omega}. \qquad (15)$$

The notation $\sigma$ in the above equation designates the typical frequency of tidal flows, which will serve as the tidal frequency of the considered tidal force in the Fourier expansion of tidal quantities detailed further. If $\tilde{\sigma} \ll 1$ (sub-inertial regime), the acceleration term of the momentum equation can be neglected with respect to Coriolis terms. Conversely, if $\tilde{\sigma} \gg 1$ (super-inertial regime), the flow is strongly driven by the tidal forcing and it is thus hardly distorted by the planet's rotation. We note that $\tilde{\sigma}$ is the inverse of the so-called spin parameter that is commonly used to characterise the oscillations of rotating fluids in planetary and stellar hydrodynamics (e.g. Lee & Saio 1997). The dimensionless parameters that control the planet's tidal response are summarised in Table 1.

## 2.4. Helmholtz decomposition

Proudman (1920) demonstrated that Helmholtz's theorem (e.g. Arfken & Weber 2005) can be used to decompose the horizontal

<footnote n="2">² Other values of $\sigma_R$ may be found in the literature. For example, Wunsch et al. (1997) uses a much smaller value to study the dynamics of the long-period tides, $\sigma_R = 2.5 \times 10^{-7}$ s⁻¹, while Motoyama et al. (2020) set this value to $\sigma_R \sim 4 \times 10^{-6}$ s⁻¹ following the prescription given by Schwiderski (1980).</footnote>
displacement vector field into curl-free and divergence-free vector fields,

$$ \tilde{\xi} = \tilde{\nabla}\Phi + \tilde{\nabla}\Psi \times \mathbf{e}_r. \qquad (16) $$

The curl-free $(\tilde{\nabla} \times (\tilde{\nabla}\Phi) = 0)$ and divergence-free $(\tilde{\nabla} \cdot (\tilde{\nabla}\Psi \times \mathbf{e}_r) = 0)$ components of $\tilde{\xi}$ are defined from the divergent displacement potential $\Phi$ and the rotational displacement stream-function $\Psi$, respectively (Webb 1980; Tyler 2011). We emphasise that the Helmholtz decomposition is not unique for bounded domains such as the considered ocean basin. This results from the fact that additional physical constraints on the boundary condition are necessary to define $\Phi$ and $\Psi$ (e.g. Fox-Kemper et al. 2003). However choosing boundary conditions for the two functions is not straightforward given that the two components of Eq. (16) cannot be disentangled in the total flux. For instance, the impermeability condition at the coastline is formulated as $\tilde{\xi} \cdot \mathbf{n} = 0$, where $\mathbf{n}$ designates the outward pointing unit vector defining the normal to the coast. This theoretically requires to find another boundary condition and to solve it for $\tilde{\nabla}\Phi$ and $\tilde{\nabla}\Psi$ together with the first condition, which may lead to significant mathematical complications.

To circumvent this difficulty, it is convenient to adopt a commonly used trick (e.g. Webb 1980, 1982; Gent & McWilliams 1983; Watterson 2001; Han & Huang 2020; Farhat et al. 2022a), which consists in applying the impermeability condition to both components of Eq. (16),

$$ \mathbf{n} \cdot \tilde{\nabla}\Phi = 0, \qquad \mathbf{n} \cdot (\tilde{\nabla}\Psi \times \mathbf{e}_r) = 0. \qquad (17) $$

The first condition of Eq. (17) means that the gradient of $\Phi$ is zero in the direction normal to the coastline, which is equivalent to Neumann condition (specified value of the derivative of the solution applied at the boundary of the domain; e.g. Morse & Feshbach 1953, Sect. 6.1). The second condition may be reformulated as $(\mathbf{e}_r \times \mathbf{n}) \cdot \tilde{\nabla}\Psi = 0$, implying that the streamfunction is a constant along the coastline. This corresponds to a Dirichlet condition (specified value of the solution itself; Morse & Feshbach 1953, Sect. 6.1). Since both $\Phi$ and $\Psi$ are defined to a constant, the value of $\Psi$ at the coastline is set to zero, which simplifies the second condition of Eq. (17) to $\Psi = 0$. Interestingly, this condition induces the orthogonality of the curl-free and divergence-free components of the horizontal displacement (e.g. Farhat et al. 2022a, Appendix E), formulated as

$$ \int_O \overline{\tilde{\nabla}\Phi} \cdot (\tilde{\nabla}\Psi \times \mathbf{e}_r) dS = 0, \qquad (18) $$

with $dS = \sin \theta d\theta d\varphi$ being an infinitesimal surface element of the unit sphere, $O$ the domain occupied by the ocean basin, and $\tilde{\nabla}\Phi$ the complex conjugate of $\tilde{\nabla}\Phi$. This property leads to appreciable simplifications in the LTEs, as discussed further.

The divergent potential function and the streamfunction are defined on a compact connected domain, the ocean basin ($O$), and they satisfy either Neumann or Dirichlet conditions at its boundary, $\partial O$. As a consequence, the two functions can be expanded in terms of the complete sets of orthogonal eigenfunctions $\{\Phi_k\}_{1 \le k \le \infty}$ or $\{\Psi_k\}_{1 \le k \le \infty}$ that are the solutions of the wave equations given by

$$ (\tilde{\nabla}^2 + \lambda) \Phi = 0 \text{ on } O, \qquad \mathbf{n} \cdot \tilde{\nabla}\Phi = 0 \text{ at } \partial O, \qquad (19) $$

$$ (\tilde{\nabla}^2 + \nu) \Psi = 0 \text{ on } O, \qquad \Psi = 0 \text{ at } \partial O, \qquad (20) $$

where $\tilde{\nabla}^2$ designates the normalised Laplacian operator defined, for any function $f$, as

$$ \tilde{\nabla}^2 f = (\sin \theta)^{-2} [\sin \theta \partial_{\theta} (\sin \theta \partial_{\theta} f) + \partial_{\varphi\varphi} f]. \qquad (21) $$

The solutions of the colatitudinal component of Eqs. (19) and (20) are the associated Legendre functions of the first kind (ALFs) of real degrees, which are detailed in Appendix B. The corresponding eigenfunctions $\Phi_k$ are associated with the eigenvalues $\lambda_k \in \mathbb{R}^+$, and the eigenfunctions $\Psi_k$ with the eigenvalues $\nu_k \in \mathbb{R}^+$, the relationship between the real degrees of the ALFs and the eigenvalues being detailed in Appendix C (see Table C.2). These eigenfunctions are obtained by multiplying the ALFs to the solutions of the longitudinal component of Eqs. (19) and (20), namely $\exp(\pm im\varphi)$, with $i$ being the imaginary number. The outcome is a set of functions known as spherical cap harmonics (Haines 1985; Hwang & Chen 1997; Thébault et al. 2004, 2006), and denoted by 'SCHs' in the following (see Eq. (C.19)). Thus, the set $\{\Phi_k\}_{1 \le k \le \infty}$ is referred to as the 'SCHNs' (Neumann condition), and the set $\{\Psi_k\}_{1 \le k \le \infty}$ as the 'SCHDs' (Dirichlet condition). For example, Fig. 2 shows the first basis functions of the two sets for a supercontinent of angular radius $\theta_c = 10^\circ$ (i.e. $\theta_0 = 170^\circ$).

One shall notice that the SCHs and their associated eigenvalues both depend on the boundary condition applied at the coastline, meaning that the sets $\{\Phi_k\}_{1 \le k \le \infty}$ and $\{\Psi_k\}_{1 \le k \le \infty}$ are necessarily different by construction. In each set, the basis functions are orthogonal to each other and normalised so that

$$ \int_O \overline{\Phi_k} \Phi_j dS = \int_O \overline{\Psi_k} \Psi_j dS = \delta_{k,j}, \qquad (22) $$

where $\delta_{k,j}$ designates the Kronecker delta function, such that $\delta_{k,j} = 1$ for $k = j$ and $\delta_{k,j} = 0$ otherwise. However, two functions belonging to different sets are not orthogonal in the general case,

$$ \int_O \overline{\Phi_k} \Psi_j dS \neq 0, \qquad (23) $$

since the sets of basis functions $\{\Phi_k\}_{1 \le k \le \infty}$ and $\{\Psi_k\}_{1 \le k \le \infty}$ are not the same.

Introducing the time-dependent coefficients $p_k$ and $p_{-k}$, we write down the divergent potential function and the streamfunction as

$$ \Phi(\theta, \varphi, \tilde{\tau}) = \sum_{k=1}^{\infty} p_k(\tilde{\tau}) \Phi_k(\theta, \varphi), \qquad (24) $$

$$ \Psi(\theta, \varphi, \tilde{\tau}) = \sum_{k=1}^{\infty} p_{-k}(\tilde{\tau}) \Psi_k(\theta, \varphi). \qquad (25) $$

By substituting Eq. (24) into Eq. (13), the continuity equation becomes

$$ \tilde{\zeta} + \tilde{\nabla}^2 \Phi = 0, \qquad (26) $$

which implies that

$$ \tilde{\zeta} = \sum_{k=1}^{\infty} \lambda_k p_k \Phi_k. \qquad (27) $$

However, such a simple expression cannot be obtained for the equilibrium displacement ($\tilde{\zeta}_{eq}$) in the general case since it is
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Eigenfunctions describing the oceanic tidal response for a supercontinent of angular radius $\theta_c = 10^\circ$ (i.e. $\theta_0 = 170^\circ$), and associated tidal flows. Top: set $\{\Phi_k\}$ (SCHNs). Bottom: set $\{\Psi_k\}$ (SCHDs). The orange disk designates the continent. The eigenfunctions are computed from the expression given by Eq. (C.19), and their real parts are plotted for $l_n$ such that $n = 0, \dots, 3$ (vertical axis) and $-n \le m \le n$ (horizontal axis). Bright or dark colours designate positive or negative values of the eigenfunctions, respectively. Streamlines indicate the tidal flows corresponding to $\tilde{\nabla}\Phi_k$ for the set $\{\Phi_k\}$ and to $\tilde{\nabla}\Psi_k \times e_r$ for the set $\{\Psi_k\}$.</figcaption></figure>

naturally expanded in series of spherical harmonics (hereafter, SPHs) in the coordinate system associated with the planet's spin, which do not satisfy Dirichlet or Neumann conditions at coastlines. The form of the expression given by Eq. (27) can be obtained for $\tilde{\zeta}_{eq}$ only in the particular case of the hemispheric ocean ($\theta_0 = 90^\circ$), where the set of eigenfunctions $\{\Phi_k\}_{1\le k \le \infty}$ actually corresponds to a subset of the SPHs for the hemispherical domain, as shown by Farhat et al. (2022a).

## 2.5. Temporal equations

The method used to establish the temporal differential equations for the coefficients $p_k$ and $p_{-k}$ introduced in Eqs. (24) and (25) closely follows that used for the hemispheric ocean configuration (e.g. Webb 1980, 1982; Farhat et al. 2022a). As a first step, we

substitute $\Phi$ and $\Psi$ by their series expansions in Eq. (16) and in the momentum equation given by Eq. (12). As a second step, we do the dot product of the equation thus obtained by $\tilde{\nabla}\Phi_k$, and we integrate it over the domain of the ocean basin. It follows

$$ \sum_{j=1}^{\infty} \left\{ \begin{aligned} & [\partial_{\tilde{\pi}} + \tilde{\sigma}_R \partial_{\tilde{\tau}}] p_j \int_O \overline{\tilde{\nabla} \Phi_k} \cdot \tilde{\nabla} \Phi_j \, dS \\ & + (\partial_{\tilde{\pi}} + \tilde{\sigma}_R \partial_{\tilde{\tau}}) p_{-j} \int_O \overline{\tilde{\nabla} \Phi_k} \cdot (\tilde{\nabla} \Psi_j \times e_r) \, dS \\ & + \partial_{\tilde{\tau}} p_{-j} \int_O \overline{\tilde{\nabla} \Phi_k} \cdot [\tilde{f} \times (\tilde{\nabla} \Psi_j \times e_r)] \, dS \\ & + \partial_{\tilde{\tau}} p_j \int_O \overline{\tilde{\nabla} \Phi_k} \cdot (\tilde{f} \times \tilde{\nabla} \Phi_j) \, dS \end{aligned} \right\} + \tilde{\sigma}_G^2 \int_O \overline{\tilde{\nabla} \Phi_k} \cdot \tilde{\nabla} \Xi \, dS = 0. \quad (28) $$
given by Eq. (4) can be expanded in Fourier series of time and series of SPHs. In the general case, $U_T$ is expressed as

$$U_T = V_0^2 \mathcal{R} \left\{ \sum_{\sigma} \sum_{l=2}^{+\infty} \sum_{m=-l}^{l} \tilde{U}_{T;l}^{m,\sigma} \hat{Y}_l^m (\hat{\theta}, \hat{\phi}) e^{i\sigma\tau} \right\}, \quad (44)$$

where $\sigma$ is the tidal frequency, $V_0$ the reference velocity introduced in Eq. (10), $\hat{Y}_l^m$ the complex SPH of degree $l$ and order $m$ associated with the coordinate system $(\hat{\theta}, \hat{\phi})$, defined in Eq. (C.9), and $\tilde{U}_{T;l}^{m,\sigma} \in \mathbb{C}$ the associated frequency-dependent normalised component of the tidal gravitational potential. In the following, these notations are shortened for convenience. Similarly as the eigenfunctions of the ocean basin, $\Phi_k$ and $\Psi_k$, the SPHs are simply denoted by $\hat{Y}_k = \hat{Y}_{l_k}^{m_k}$, and the corresponding coefficients by $\tilde{U}_{T;k}^{\sigma} = U_{T;l_k}^{m_k,\sigma}$, the index $k$ referring to an element of the set of SPHs, $\{\hat{Y}_k\}_{1\le k \le \infty}$. Moreover, all the indices used from now on $(k, j, n, q)$ are supposed to run from one to infinity when bounds are not specified. It is noteworthy that the Fourier series expansion in Eq. (44) can always be written in term of positive tidal frequencies as long as $m$ runs from $-l$ to $l$. Therefore, we assume that $\sigma \ge 0$, and that the frequencies are all different from each other (no resonance).

Since the components associated with two different tidal frequencies are not correlated in the linear tidal theory, the resulting tidal responses can be treated separately. We thus consider the contribution of the component associated with a given tidal frequency, $\tilde{U}_T^\sigma$, and the associated forcing term, $\tilde{\Xi}^\sigma$, expressed as

$$\tilde{U}_T^\sigma = \sum_j \tilde{U}_{T;j}^\sigma \hat{Y}_j (\hat{\theta}, \hat{\phi}), \quad \tilde{\Xi}^\sigma = \sum_j \tilde{\Xi}_j^\sigma \hat{Y}_j (\hat{\theta}, \hat{\phi}). \quad (45)$$

As highlighted in earlier studies (e.g. Matsuyama 2014; Matsuyama et al. 2018; Auclair-Desrotour et al. 2019a), the combined contribution of ocean loading, self-attraction variation, and deformation of solid regions may be formulated as

$$\tilde{\Xi}_j^\sigma = \gamma_{D;l_j}^\sigma \left( \sum_n \langle \hat{Y}_j, \Phi_n \rangle \lambda_n p_n \right) - \gamma_{T;l_j}^\sigma \tilde{U}_{T;j}^\sigma, \quad (46)$$

where we have made use of the scalar product defined by Eq. (C.12) and introduced the solid deformation factors, $\gamma_{D;l}^\sigma$ and $\gamma_{T;l}^\sigma$ (see e.g. Hendershott 1972),

$$\gamma_{T;l}^\sigma \equiv 1 + k_l^\sigma - h_l^\sigma, \quad (47)$$

$$\gamma_{D;l}^{\sigma} \equiv 1 - (1 + k_{L;l}^{\sigma} - h_{L;l}^{\sigma}) \frac{3}{(2l+1)} \frac{\rho_{w}}{\bar{\rho}}. \quad (48)$$

The factors $\gamma_{D;l}^\sigma$ and $\gamma_{T;l}^\sigma$ are independent of the order $m$ since they are expressed as functions of the Love numbers of the solid part, which is commonly assumed to be spherically symmetric in the tidal theory. The parameters $k_l^\sigma$ and $h_l^\sigma$ in the first expression are known as the tidal gravitational and displacement Love numbers, respectively. They describe the response of the solid body to a gravitational tidal forcing, which takes the form of a self-attraction variation and a surface displacement. By analogy, the loading Love numbers $k_{L;l}^\sigma$ and $h_{L;l}^\sigma$ account for the gravitational and mechanical responses of solid regions to the cumulated gravitational and pressure forces generated by the oceanic mass redistribution, respectively. The second solid deformation factor, $\gamma_{T;l}^\sigma$, also depends on the ratio of seawater

density to the mean density of the solid regions, $\bar{\rho}$. We remark that $\gamma_{D;l}^\sigma$ and $\gamma_{T;l}^\sigma$ are sometimes considered as real factors (e.g. Matsuyama 2014; Motoyama et al. 2020), which corresponds to an adiabatic elastic response of the solid part, where dissipative mechanisms are ignored. However, the correspondence principle established by Biot (1954) makes it possible to treat dissipative anelastic cases similarly as long as the anelasticity is linear. As a consequence, the expressions given by Eqs. (47) and (48) can be extended to any rheological model including dissipative processes. In that case, the solid Love numbers are complex transfer functions describing the visco-elastic response of solid regions when subjected to a harmonic tidal force (Remus et al. 2012; Auclair-Desrotour et al. 2019a; Farhat et al. 2022b).

Using the expression given by Eq. (46) and proceeding to a change of basis functions between the oceanic eigenfunctions and the SPHs associated with the coordinate system $(\hat{\theta}, \hat{\phi})$, we expand the integrals depending on the forcing in Eqs. (38) and (39) in series of $p_j$ and components of the tidal gravitational potential,

$$\int_O \overline{\Phi_k \tilde{\Xi} dS} = \sum_n \left[ \sum_j \langle \Phi_k, \hat{Y}_j \rangle \gamma_{D;l_j}^\sigma \langle \hat{Y}_j, \Phi_n \rangle \right] \lambda_n p_n \quad (49) \\ - \sum_j \langle \Phi_k, \hat{Y}_j \rangle \gamma_{T;l_j}^\sigma \tilde{U}_{T;j}^\sigma,$$

$$\int_O \mathbf{e}_r \cdot (\overline{\tilde{\nabla} \Psi_k} \times \tilde{\nabla} \tilde{\Xi}) dS = \sum_n \left[ \sum_{j,q} \nu_{k,q} \langle Y_q, \hat{Y}_j \rangle \gamma_{D;l_j}^\sigma \langle \hat{Y}_j, \Phi_n \rangle \right] \lambda_n p_n \\ - \sum_{j,q} \nu_{k,q} \langle Y_q, \hat{Y}_j \rangle \gamma_{T;l_j}^\sigma \tilde{U}_{T;j}^\sigma, \quad (50)$$

where $\nu_{k,q}$ designates the coefficients defined as

$$\nu_{k,q} \equiv \int_O \mathbf{e}_r \cdot (\overline{\tilde{\nabla} \Psi_k} \times \tilde{\nabla} Y_q) dS. \quad (51)$$

These coefficients are all zero as shown in Appendix G, which implies that

$$\int_O \mathbf{e}_r \cdot (\overline{\tilde{\nabla} \Psi_k} \times \tilde{\nabla} \tilde{\Xi}) dS = 0. \quad (52)$$

We remark that the quadrupolar component of the tidal gravitational potential (the component associated with the $\hat{Y}_2^2$ SPH) is far greater than the components of higher degrees if the size of the planet is small compared with the planet-perturber distance. These components can thus be ignored in standard two-body systems. The series in the formulation of the tidal gravitational potential given by Eq. (45) then reduces to one term only. This removes the summation over $j$ in the second term of the right-hand member of Eq. (49).

## 2.7. The tidal solution

As a last step, the tidal equations given by Eqs. (38) and (39) are written down in the Fourier domain. To compute the solution numerically, all the sets of basis functions are truncated, meaning that the infinite spaces of functions they describe are approximated by finite spaces of functions. As the tidal gravitational potential varies over planetary length scales, the tidal response is essentially described by the eigenmodes of largest wavelengths. As a consequence, the truncation does not alter much the solution as long as the number of eigenmodes of the set is sufficiently
large. The numbers of basis functions for the sets {$\Phi_j$}, {$\Psi_j$}, and {$\hat{\Psi}_j$} are denoted by $M$, $N$, and $K$, respectively, while the vectors describing the tidal response and the tidal gravitational potential in these sets are expressed as

$$ \mathbf{\Phi} = [p_1, \dots, p_j, \dots, p_M]^T, \quad (53) $$

$$ \mathbf{\Psi} = [p_{-1}, \dots, p_{-j}, \dots, p_{-N}]^T, \quad (54) $$

$$ \tilde{\mathbf{U}}_{\mathrm{T}} = [\tilde{U}_{\mathrm{T},1}^{\sigma}, \dots, \tilde{U}_{\mathrm{T},j}^{\sigma}, \dots, \tilde{U}_{\mathrm{T},K}^{\sigma}]^{\mathrm{T}}, \quad (55) $$

with $^T$ designating the transpose of a matrix. Converted into an algebraic form, the forcing term given by Eq. (49) is written as

$$ \int_{\mathcal{O}} \overline{\Phi_k \Xi} dS = \mathbf{E}_{M,k}^{\mathrm{T}} (\mathbf{A}_{\Phi} \mathbf{\Phi} + \mathbf{R}_{\Phi} \tilde{\mathbf{U}}_{\mathrm{T}}), \quad (56) $$

where we have introduced the unit vectors $E_{M,k}$ defined as

$$ \mathbf{E}_{M,k} \equiv [\delta_{1,k}, \dots, \delta_{j,k}, \dots, \delta_{M,k}]^T, \quad (57) $$

and the matrices $\mathbf{A}_{\Phi}$, and $\mathbf{R}_{\Phi}$,

$$ \mathbf{A}_{\Phi} = \overline{\mathbf{P}_{\hat{\Phi},\Phi}}^{\mathrm{T}} \mathbf{\Gamma}_{\mathrm{D}} \mathbf{P}_{\hat{\Phi},\Phi} \mathbf{\Lambda}, \quad \mathbf{R}_{\Phi} = - \overline{\mathbf{P}_{\hat{\Phi},\Phi}}^{\mathrm{T}} \mathbf{\Gamma}_{\mathrm{T}}. \quad (58) $$

In the above expressions, $\mathbf{P}_{\hat{\Phi},\Phi}$ designates the transition matrix from the SPHs $\hat{\mathbf{Y}}_k$ to the eigenfunctions $\Phi_k$. The other matrices are defined as

$$ \mathbf{\Gamma}_{\mathrm{D}} \equiv \begin{bmatrix} \gamma_{\mathrm{D};l_1}^{\sigma} & 0 & \dots & \dots & 0 \\ 0 & \ddots & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & \dots & 0 & \gamma_{\mathrm{D};l_K}^{\sigma} \end{bmatrix}, \quad (59) $$

$$ \mathbf{\Gamma}_{\mathrm{T}} \equiv \begin{bmatrix} \gamma_{\mathrm{T};l_1}^{\sigma} & 0 & \dots & \dots & 0 \\ 0 & \ddots & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & \dots & 0 & \gamma_{\mathrm{T};l_K}^{\sigma} \end{bmatrix}, \quad (60) $$

$$ \mathbf{\Lambda} \equiv \begin{bmatrix} \lambda_1 & 0 & \dots & \dots & 0 \\ 0 & \ddots & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & \dots & 0 & \lambda_M \end{bmatrix}, \quad (61) $$

$$ \mathbf{N} \equiv \begin{bmatrix} \nu_1 & 0 & \dots & \dots & 0 \\ 0 & \ddots & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & \dots & 0 & \nu_N \end{bmatrix}. \quad (62) $$

Also, we introduce the matrix of gyroscopic coefficients

$$ \mathbf{B} \equiv \begin{bmatrix} \Lambda^{-1} \mathbf{B}_{\Phi,\Phi} & \Lambda^{-1} \mathbf{B}_{\Phi,\Psi} \\ N^{-1} \mathbf{B}_{\Psi,\Phi} & N^{-1} \mathbf{B}_{\Psi,\Psi} \end{bmatrix}, \quad (63) $$

where the block matrices $\mathbf{B}_{\Phi,\Phi}$, $\mathbf{B}_{\Phi,\Psi}$, $\mathbf{B}_{\Psi,\Phi}$, and $\mathbf{B}_{\Psi,\Psi}$ are defined in terms of the $\beta_{k,j}$ defined in Eqs. (40)–(43) as

$$ \mathbf{B}_{\Phi,\Phi} \equiv \begin{bmatrix} \beta_{1,1} & \dots & \beta_{1,j} & \dots & \beta_{1,M} \\ \vdots & & \vdots & & \vdots \\ \beta_{k,1} & \dots & \beta_{k,j} & \dots & \beta_{k,M} \\ \vdots & & \vdots & & \vdots \\ \beta_{M,1} & \dots & \beta_{M,j} & \dots & \beta_{M,M} \end{bmatrix}, \quad (64) $$

$$ \mathbf{B}_{\Phi,\Psi} \equiv \begin{bmatrix} \beta_{1,-1} & \dots & \beta_{1,-j} & \dots & \beta_{1,-N} \\ \vdots & & \vdots & & \vdots \\ \beta_{k,-1} & \dots & \beta_{k,-j} & \dots & \beta_{k,-N} \\ \vdots & & \vdots & & \vdots \\ \beta_{M,-1} & \dots & \beta_{M,-j} & \dots & \beta_{M,-N} \end{bmatrix}, \quad (65) $$

$$ \mathbf{B}_{\Psi,\Phi} \equiv \begin{bmatrix} \beta_{-1,1} & \dots & \beta_{-1,j} & \dots & \beta_{-1,M} \\ \vdots & & \vdots & & \vdots \\ \beta_{-k,1} & \dots & \beta_{-k,j} & \dots & \beta_{-k,M} \\ \vdots & & \vdots & & \vdots \\ \beta_{-N,1} & \dots & \beta_{-N,j} & \dots & \beta_{-N,M} \end{bmatrix} = - \overline{\mathbf{B}_{\Phi,\Psi}}^{\mathrm{T}}, \quad (66) $$

$$ \mathbf{B}_{\Psi,\Psi} \equiv \begin{bmatrix} \beta_{-1,-1} & \dots & \beta_{-1,-j} & \dots & \beta_{-1,-N} \\ \vdots & & \vdots & & \vdots \\ \beta_{-k,-1} & \dots & \beta_{-k,-j} & \dots & \beta_{-k,-N} \\ \vdots & & \vdots & & \vdots \\ \beta_{-N,-1} & \dots & \beta_{-N,-j} & \dots & \beta_{-N,-N} \end{bmatrix}. \quad (67) $$

Finally, the matrix accounting for the gravitational and pressure effects of the oceanic surface displacement, $\mathbf{A}$, the matrix describing the coupling between the forcing tidal gravitational potential and the oceanic eigenmodes, $\mathbf{R}$, and the solution vector, $\mathbf{X}$, are respectively expressed as

$$ \mathbf{A} \equiv \begin{bmatrix} \mathbf{A}_{\Phi} & 0 \\ 0 & 0 \end{bmatrix}, \quad \mathbf{R} \equiv \begin{bmatrix} \mathbf{R}_{\Phi} \\ 0 \end{bmatrix}, \quad \mathbf{X} \equiv \begin{bmatrix} \mathbf{\Phi} \\ \mathbf{\Psi} \end{bmatrix}. \quad (68) $$

The temporal tidal equations given by Eqs. (38) and (39) thus lead to the solution

$$ \mathbf{X} = \mathbf{H}(\tilde{\sigma}, \tilde{\sigma}_R, \tilde{\sigma}_G, \theta_c, \hat{\theta}_S; \gamma_{D;l}^{\sigma}, \gamma_{T;l}^{\sigma}) \tilde{\mathbf{U}}_{\mathrm{T}}, \quad (69) $$

$$ \mathbf{H} \equiv \tilde{\sigma}_G^2 [\tilde{\sigma}(\tilde{\sigma} - i\tilde{\sigma}_R)\mathbf{I} - i\tilde{\sigma}\mathbf{B} - \tilde{\sigma}_G^2\mathbf{A}]^{-1} \mathbf{R}, \quad (70) $$

with $\mathbf{I}$ being the identity matrix. In the above expression, the oceanic tidal response to the tidal perturbation in the Fourier domain is expressed as a complex frequency-dependent matrix ($\mathbf{H}$), which links the tidal gravitational potential ($\tilde{\mathbf{U}}_{\mathrm{T}}$) to the potential and stream functions that describe the horizontal tidal displacement ($\mathbf{X}$).
## 2.8. Tidal Love numbers, torque, and power

To close this theoretical section, it remains to introduce the parameters that quantify the tidally dissipated energy and its effect on the evolution of planetary systems. The tidal mass redistribution itself is quantified by self-attraction variations, which is derived from the gravitational potential of the tidally distorted body, $U_D$, defined at its surface as

$$U_D = V_0^2 \Re \left\{ \sum_{\sigma} \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \tilde{U}_{D;l}^{m,\sigma} \hat{Y}_l^m (\hat{\theta}, \hat{\phi}) e^{i\sigma t} \right\}. \quad (71)$$

The components $\tilde{U}_{D;l}^{m,\sigma}$ of this gravitational potential are expressed as

$$\tilde{U}_{D;l}^{m,\sigma} = k_l^{\sigma} \tilde{U}_{T;l}^{m,\sigma} + \left(1 + k_{L;l}^{\sigma}\right) \frac{3}{2l+1} \frac{\rho_w}{\bar{\rho}} \tilde{\zeta}_l^{m,\sigma}, \quad (72)$$

where $\tilde{U}_{T;l}^{m,\sigma}$ is the corresponding component of the forcing tidal potential introduced in Eq. (44), and $\tilde{\zeta}_l^{m,\sigma}$ the corresponding component of the normalised surface elevation,

$$\tilde{\zeta}_l^{m,\sigma} = \sum_{k=1}^{\infty} \langle \hat{Y}_l^m, \Phi_k \rangle \lambda_k p_k. \quad (73)$$

The first term of Eq. (72) is the self-attraction variation due to the distortion of the solid part generated by the gravitational tidal force, while the second term is the contribution of the ocean loading, which includes both the potential generated by the oceanic mass redistribution and the potential resulting from the distortion of the solid part induced by ocean loading.

The complex Love numbers are the transfer functions accounting for the intrinsic response of an extended body undergoing a tidal gravitational forcing (e.g. Ogilvie 2014). They are defined, for each SPH, as the ratio of the gravitational potential of the tidal response to the forcing tidal potential evaluated at the body's surface,

$$k_{D;l}^{m,\sigma} (\tilde{\sigma}, \tilde{\sigma}_R, \tilde{\sigma}_G, \theta_c, \hat{\theta}_S; \gamma_{D;l}^{\sigma}, \gamma_{T;l}^{\sigma}) = \frac{\tilde{U}_{D;l}^{m,\sigma}}{\tilde{U}_{T;l}^{m,\sigma}}, \quad (74)$$

which, using the expression of the tidal potential given by Eq. (72), yields

$$k_{D;l}^{m,\sigma} = k_l^{\sigma} + \left(1 + k_{L;l}^{\sigma}\right) \frac{3}{2l+1} \frac{\rho_w}{\bar{\rho}} \frac{\tilde{\zeta}_l^{m,\sigma}}{\tilde{U}_{T;l}^{m,\sigma}}. \quad (75)$$

Considering the above equation, we note that the Love number depends both on the order $m$ and the degree $l$ of the associated SPH in the general case contrary to the solid Love numbers, which only depend on the degrees (see Eqs. (47) and (48)). This is due to the fact that the oceanic tidal response is not spherically symmetric owing to the geometry of coastlines and Coriolis effects.

The long-term effect of tides on the evolution of the planet-perturber system is quantified by the generated tidal torques. The main contributor to this evolution is the time-averaged tidal torque exerted about the spin axis of the planet (Zahn 1966),

$$\mathcal{T}_z = \left\langle \int_{V} \partial_{\hat{\phi}} U_T \delta \rho dV \right\rangle, \quad (76)$$

which affects the long-term variation rate of the planet-perturber distance and the planet's spin angular velocity (e.g. Ogilvie 2014). In the above equation, $\mathcal{V}$ designates the volume filled by the planet, $dV$ an infinitesimal volume parcel, and $\delta\rho$ the local density variations associated with the tidal response. These density variations are related through Poisson's equation to the tidal gravitational potential of the distorted planet, $U_D$, namely

$$\nabla^2 U_D = -4\pi G \delta\rho. \quad (77)$$

This allows the tidal torque to be expressed in terms of the forcing and deformation tidal gravitational potentials,

$$\mathcal{T}_z = -\frac{1}{4\pi G} \left\langle \int_V \partial_{\hat{\phi}} U_T \nabla^2 U_D dV \right\rangle. \quad (78)$$

The latter integral is defined over the volume filled by the tidally deformed planet ($\mathcal{V}$), but it can actually be carried out over any region that includes this volume. By making use of this property, Ogilvie (2013) demonstrates that both the tidal torque and the tidally dissipated power are straightforwardly related to the Fourier components of the two gravitational potentials. As a first step, one shall notice that $\nabla^2 U_T = 0$ outside of the perturber, by virtue of Poisson's equation, which implies that

$$\nabla^2 (\partial_{\hat{\phi}} U_T) = \partial_{\hat{\phi}} (\nabla^2 U_T) = 0, \quad (79)$$

since the Laplacian and $\partial_{\hat{\phi}}$ operators can be permuted. Let $\mathcal{V}_*$ be any region that is simply connected and does not include the perturber. Then, the fact that $\nabla^2 (\partial_{\hat{\phi}} U_T) = 0$ in $\mathcal{V}_*$ allows the integral of Eq. (78) to be rewritten as

$$\int_{V_*} \partial_{\hat{\phi}} U_T \nabla^2 U_D dV = \int_{V_*} [\partial_{\hat{\phi}} U_T \nabla^2 U_D - U_D \nabla^2 (\partial_{\hat{\phi}} U_T)] dV. \quad (80)$$

By virtue of Green's second identity (e.g. Arfken & Weber 2005, Sect. 1.11), Eq. (80) yields

$$\int_{V_*} \partial_{\hat{\phi}} U_T \nabla^2 U_D dV = \oint_{\partial V_*} [\partial_{\hat{\phi}} U_T \nabla U_D - U_D \nabla (\partial_{\hat{\phi}} U_T)] \cdot dS, \quad (81)$$

with the notation $\partial \mathcal{V}_*$ referring to the boundary of the domain $\mathcal{V}_*$, and $dS$ to and outward pointing infinitesimal surface element vector.

By applying the time-averaging operator to the integral we obtain

$$\mathcal{T}_z = \sum_{\sigma} \frac{V_0^4}{8\pi G} \Re \left\{ \oint_{\partial V_*} \left[ \overline{\tilde{U}_D^\sigma} \nabla (\partial_{\hat{\phi}} \tilde{U}_T^\sigma) - \overline{\partial_{\hat{\phi}} \tilde{U}_T^\sigma} \nabla \tilde{U}_D^\sigma \right] \cdot dS \right\}, \quad (82)$$

with $\tilde{U}_T^\sigma$ and $\tilde{U}_D^\sigma$ defined as

$$U_T (r, \hat{\theta}, \hat{\phi}, t) = V_0^2 \Re \left\{ \sum_{\sigma} \tilde{U}_T^\sigma (r, \hat{\theta}, \hat{\phi}) e^{i\sigma t} \right\}, \quad (83)$$

$$U_D (r, \hat{\theta}, \hat{\phi}, t) = V_0^2 \Re \left\{ \sum_{\sigma} \tilde{U}_D^\sigma (r, \hat{\theta}, \hat{\phi}) e^{i\sigma t} \right\}. \quad (84)$$

In practice, $\mathcal{V}_*$ can be chosen such that its boundary is a sphere of some radius $r_*$ centred on the planet.

In the vicinity of $\partial \mathcal{V}_*$, the forcing tidal potential and the deformation tidal potential can be represented using interior and exterior multipole expansions in SPHs, respectively,

$$\tilde{U}_T^\sigma (r, \hat{\theta}, \hat{\phi}) = \sum_{l=2}^{+\infty} \sum_{m=-l}^{l} \tilde{U}_{T;l}^{m,\sigma} \left( \frac{r}{R_p} \right)^l \hat{Y}_l^m (\hat{\theta}, \hat{\phi}), \quad (85)$$
$$ \tilde{U}_{\mathrm{D}}^{\sigma}(r, \hat{\theta}, \hat{\varphi}) = \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \tilde{U}_{\mathrm{D};l}^{m,\sigma} \left(\frac{r}{R_{\mathrm{p}}}\right)^{-(l+1)} \hat{Y}_{l}^{m}(\hat{\theta}, \hat{\varphi}). \quad (86) $$

By substituting Eqs. (85) and (86) into Eq. (82), and for arbitrary radius $r_*$, we end up with the formulation of the torque as a function of the imaginary part of Love numbers,

$$ \mathcal{T}_z = \sum_{\sigma} \sum_{l=2}^{+\infty} \sum_{m=-l}^{l} \frac{m(2l+1)R_p V_0^4}{8\pi G} |\tilde{U}_{T;l}^{m,\sigma}|^2 \mathfrak{I} \{k_{D;l}^{m,\sigma}\}. \quad (87) $$

For the quadrupolar tidal gravitational potential (e.g. Ogilvie 2014),

$$ \tilde{U}_{T;2}^{2,\sigma} = \sqrt{\frac{6\pi}{5}} \left( \frac{GM_s}{aV_0^2} \right) \left( \frac{R_p}{a} \right)^2, \quad (88) $$

with $a$ being the semi-major axis of the perturber, we recover the well-known relationship between the tidal torque and the quadrupolar Love number (e.g. Makarov & Efroimsky 2013; Auclair-Desrotour et al. 2019a),

$$ \mathcal{T}_{z;2} = \frac{3}{2}GM_s^2 \frac{R_p^5}{a^6} \mathfrak{I} \{k_{D;2}^{2,\sigma}\}. \quad (89) $$

The time-averaged power input by the tidal force can be determined in a similar way. This power is defined as the work per unit time done by the tidal force on the tidal motions (e.g. Ogilvie 2013),

$$ \mathcal{P}_T = \left\langle \int_{V_*} \rho V \cdot \nabla U_T dV \right\rangle, \quad (90) $$

where $\rho$ is the local density of the planet. We note that $V$ designates here the velocity field of tidal motions in the whole planet and not only in the ocean basin. By proceeding to an integration by parts, we rewrite the above expression as

$$ \mathcal{P}_T = \oint_{\partial V_*} U_T \rho V \cdot dS - \int_{V_*} U_T \nabla \cdot (\rho V) dV. \quad (91) $$

The first term of the right-hand member in Eq. (91) is zero since the boundary $\partial V_*$ is outside of the planet, where $\rho = 0$. By combining together the continuity equation,

$$ \partial_t \delta\rho + \nabla \cdot (\rho V) = 0, \quad (92) $$

and Poisson's equation, given by Eq. (77), we express $\nabla \cdot (\rho V)$ as a function of the deformation tidal gravitational potential,

$$ \nabla \cdot (\rho V) = \frac{1}{4\pi G} \nabla^2 (\partial_t U_D), \quad (93) $$

which yields

$$ \mathcal{P}_T = -\frac{1}{4\pi G} \left\langle \int_{V_*} U_T \nabla^2 (\partial_t U_D) dV \right\rangle. \quad (94) $$

Using Green's second identity, the above expression becomes

$$ \mathcal{P}_T = \sum_{\sigma} \frac{V_0^4 \sigma}{8\pi G} \mathfrak{I} \left\{ \oint_{\partial V_*} \left[ \overline{\tilde{U}_T^{\sigma}} \nabla \tilde{U}_D^{\sigma} - \tilde{U}_D^{\sigma} \nabla \tilde{U}_T^{\sigma} \right] \cdot dS \right\}. \quad (95) $$

Finally we use the expressions of $\tilde{U}_T^\sigma$ and $\tilde{U}_D^\sigma$ given by Eqs. (85) and (86), and the definition of the Love number given by Eq. (74). It follows

$$ \mathcal{P}_T = - \sum_{\sigma} \sum_{l=2}^{+\infty} \sum_{m=-l}^{l} \frac{(2l+1) R_p V_0^4 \sigma}{8\pi G} |\tilde{U}_{T;l}^{m,\sigma}|^2 \mathfrak{I} \{k_{D;l}^{m,\sigma}\}. \quad (96) $$

The tidal input power in the ocean is expressed as a function of the harmonics of the surface tidal elevation only,

$$ \mathcal{P}_{T;\mathrm{oc}} = -\frac{3R_p V_0^4}{8\pi G} \left(\frac{\rho_w}{\bar{\rho}}\right) \sum_{\sigma} \sum_{l=2}^{+\infty} \sum_{m=-l}^{l} \sigma \mathfrak{I} \left\{ \overline{\tilde{U}_{T;l}^{m,\sigma}} \tilde{\zeta}_l^{m,\sigma} \right\}. \quad (97) $$

By noting that $\tilde{U}_T = \tilde{\zeta}_{eq}$ and expanding the component of the forcing tidal potential in series of oceanic eigenmodes,

$$ \tilde{U}_{T;l}^{m,\sigma} = \sum_{k=1}^{\infty} \langle \hat{Y}_l^m, \Phi_k \rangle \tilde{\zeta}_{eq;k}, \quad (98) $$

we can rewrite the sum on degrees and orders in Eq. (97) as

$$ \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \overline{\tilde{U}_{T;l}^{m,\sigma}} \tilde{\zeta}_l^{m,\sigma} = \sum_{k,j} \left[ \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \langle \Phi_k, \hat{Y}_l^m \rangle \langle \hat{Y}_l^m, \Phi_j \rangle \right] \overline{\tilde{\zeta}_{eq;k}} \lambda_j p_j. \quad (99) $$

Therefore, owing to the orthogonality of the $\Phi_k$ eigenfunctions, Eq. (97) simplifies to

$$ \mathcal{P}_{T;\mathrm{oc}} = \frac{1}{2} R_p^2 H^2 g \rho_w \sum_{\sigma} \sum_{k=1}^{+\infty} \sigma \mathfrak{I} \left\{ \lambda_k p_k \overline{\tilde{\zeta}_{eq;k}} \right\}, \quad (100) $$

which corresponds to the usual expression of the work per unit time done by the tidal force on the oceanic tidal motions (e.g. Webb 1980; Farhat et al. 2022a).

Similarly, the tidally dissipated power within the oceanic layer is defined as

$$ \mathcal{P}_{\mathrm{diss};\mathrm{oc}} \equiv \left\langle \int_O \rho_w H \sigma_R V \cdot V dS \right\rangle. \quad (101) $$

By making use of the orthogonality properties of the oceanic eigenmodes described by Eqs. (30) and (36), this power is expressed as (e.g. Webb 1980)

$$ \mathcal{P}_{\mathrm{diss};\mathrm{oc}} = \frac{1}{2} \sigma_R \rho_w R_p^4 H \sum_{\sigma} \sum_{k=1}^{+\infty} \sigma^2 (\lambda_k |p_k|^2 + \nu_k |p_{-k}|^2). \quad (102) $$

The dissipated power $\mathcal{P}_{\mathrm{diss};\mathrm{oc}}$ is equal to the input tidal power given by Eq. (97) if the solid regions are rigid, as assumed in earlier studies (e.g. Webb 1980; Auclair-Desrotour et al. 2018).

When the solid part is deformable, the total input tidal power given by Eq. (96), $\mathcal{P}_T$, and the total tidally dissipated power, $\mathcal{P}_{\mathrm{diss}}$, are the sum of the solid and oceanic contributions,

$$ \mathcal{P}_T = \mathcal{P}_{T;\mathrm{sol}} + \mathcal{P}_{T;\mathrm{oc}}, \quad (103) $$

$$ \mathcal{P}_{\mathrm{diss}} = \mathcal{P}_{\mathrm{diss;sol}} + \mathcal{P}_{\mathrm{diss;oc}}, \quad (104) $$

where $\mathcal{P}_{T;\mathrm{sol}} \neq \mathcal{P}_{\mathrm{diss;sol}}$ and $\mathcal{P}_{T;\mathrm{oc}} \neq \mathcal{P}_{\mathrm{diss;oc}}$ in the general case. Particularly, we note that the two components of the input tidal power can take negative values due to the energy exchanges between the ocean and solid regions resulting from the loading
<table><caption>**Table 2.** Values of parameters used in the reference case.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td colspan="4">Planet's solid part</td></tr><tr><td>Planet mass</td><td>M<sub>p</sub></td><td>1.0 M<sub>E</sub></td><td>1</td></tr><tr><td>Planet radius</td><td>R<sub>p</sub></td><td>1.0 R<sub>E</sub></td><td>1</td></tr><tr><td>Effective shear modulus</td><td>μ</td><td>25.1189 GPa</td><td>2</td></tr><tr><td>Maxwell time</td><td>τ<sub>M</sub></td><td>685 yr</td><td>2</td></tr><tr><td>Andrade time</td><td>τ<sub>A</sub></td><td>12 897.1 yr</td><td>2</td></tr><tr><td>Rheological parameter</td><td>α</td><td>0.25</td><td>2</td></tr><tr><td colspan="4">Planet's ocean basin</td></tr><tr><td>Oceanic depth</td><td>H</td><td>4.0 km</td><td></td></tr><tr><td>Continental angular radius</td><td>θ<sub>c</sub></td><td>90°</td><td></td></tr><tr><td>Colatitude of the continental centre</td><td>θ̂<sub>S</sub></td><td>90°</td><td></td></tr><tr><td>Longitude of the continental centre</td><td>φ̂<sub>S</sub></td><td>0°</td><td></td></tr><tr><td>Seawater density</td><td>ρ<sub>w</sub></td><td>1022 kg m<sup>-3</sup></td><td></td></tr><tr><td>Rayleigh drag frequency</td><td>σ<sub>R</sub></td><td>1.0 × 10<sup>-5</sup> s<sup>-1</sup></td><td></td></tr><tr><td colspan="4">Perturber</td></tr><tr><td>Mass of the perturber</td><td>M<sub>s</sub></td><td>1.0 M<sub>☉</sub></td><td>1</td></tr><tr><td>Orbital frequency of the perturber</td><td>n<sub>s</sub></td><td>1.0 n<sub>☉</sub></td><td>1</td></tr><tr><td colspan="4">Spectral method</td></tr><tr><td>Truncation degree of the SPHs</td><td>K</td><td>100</td><td></td></tr><tr><td>Truncation degree of the oceanic SCHs</td><td>M</td><td>30</td><td></td></tr></tbody></table>

**References.** (1) Mamajek et al. (2015); (2) Bolmont et al. (2020).

and gravitational coupling, while $P_{\text{diss;sol}} \ge 0$ and $P_{\text{diss;oc}} \ge 0$. In the following, the input and tidally dissipated powers of the solid regions, $P_{\text{T;sol}}$ and $P_{\text{diss;sol}}$, are evaluated from the oceanic contributions given by Eqs. (97) and (102), and the total input tidal power, expressed in Eq. (96), by using the fact that $P_{\text{diss}} = P_{\text{T}}$.

At this stage, we have established the expressions of all the tidal quantities that characterise the dissipative tidal response of the planet in the general case. In the following, we apply the theory to specific cases. The above expressions are thus used to evaluate the tidally dissipated power as a function of the forcing frequency.

### 3. Study of a reference case

Before making any attempt to explore the parameter space, it seems appropriate to elucidate the tidal response of the planet for a well-chosen reference configuration. In this section, we detail the physical setup of this configuration and we compute the frequency spectra of the quantities that characterise the semidiurnal tidal response of the planet. The obtained results are used as a starting point in the parametric study achieved in Sect. 4.

#### 3.1. Physical setup

Since the Earth-Moon system is the planet-satellite system that we know best, it appears as a privileged option for the aforementioned reference configuration. Therefore, we consider a simplified Earth-Moon system where the Moon orbits an Earth-sized planet with a circular and coplanar motion. In this framework, the Moon's trajectory is a circle in the planet's equatorial plane. The values used for the model parameters in this reference case are summarised in Table 2.

Following Webb (1980), the planet's surface is assumed to be divided into continental and oceanic hemispheres centred at the equator. Namely, the continental angular radius is set to $\theta_c = 90^\circ$ and the colatitude of the continental centre to $\theta_S = 90^\circ$. In this

configuration, the coastlines correspond exactly to meridians. The ocean basin may be regarded as a geometrically simplified version of the Pacific ocean. Its depth is set to $H = 4$ km and its Rayleigh drag frequency to $\sigma_R = 10^{-5}$ s$^{-1}$, consistently with the values found in earlier studies from constraints on the actual tidal dissipation rate of the Earth-Moon system and its age (Webb 1980, 1982; Farhat et al. 2022a). The mean density of seawater is set to $\rho_w = 1022$ kg m$^{-3}$, which is a typical value of seawater density on Earth (Gerkema & Zimmerman 2008). Finally, the solid part of the planet is assumed to behave as a visco-elastic body described by the Andrade model (see Appendix H). In this model, the solid tidal response is parametrised by an effective shear modulus, $\mu$, the Maxwell time associated with viscous friction $\tau_M$, the Andrade time associated with the inelastic component of the solid elongation, $\tau_A$, and the dimensionless parameter $\alpha$. The values used for these rheological parameters are adopted from Bolmont et al. (2020), Table 2.

In the considered system, the tidal perturber is a dimensionless body of mass $M_s = M_{\odot}$ and mean motion $n_s = n_{\odot}$, with $M_{\odot}$ and $n_{\odot}$ being the Lunar mass and mean motion, respectively. Since its orbit is circular and coplanar, the induced tidal perturbation reduces to the semidiurnal tide, which is associated with the quadrupolar tidal potential given by Eq. (88) and the forcing frequency $\sigma = 2(\Omega - n_s)$ (e.g. Ogilvie 2014). This quadrupole is coupled with the SCHs describing the oceanic tidal response, which are determined by the geometry of the ocean basin. To account for this coupling properly, the truncation degrees of the two sets of basis functions must be such that the spatial resolution reached with SPHs is higher than that reached with SCHs. Besides, the truncation degree of the SCHs shall be chosen sufficiently high to describe all the excited oceanic eigenmodes.

Preliminary tests with various truncation degrees show that convergence is reached, in all cases, for truncation degrees $M$ less than ~30. Therefore, we set the truncation degree of the SCHs to $M = 30$. Similarly, a truncation degree of $K = 100$ is
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Frequency spectra of tidal quantities and their components for the reference case (Table 2) and the semidiurnal tide. **Top:** Tidal input power (left panel) and tidally dissipated power (right panel). **Bottom:** Imaginary part of the quadrupolar Love number (left panel) and tidal torque exerted about the planet's spin axis (right panel). All quantities are plotted in logarithmic scale (vertical axis) as functions of the normalised semidiurnal tidal frequency $\chi = (\Omega - n_s)/\Omega_E$ (horizontal axis). In each plot, the orange, blue and black lines designate the response of the solid part, ocean, and whole planet, respectively. Solid lines indicate decelerating tidal torque ($\mathcal{F}(k_{D,2}^{2,\sigma}) \le 0$ and $T_{z,2} \le 0$) and positive tidal input power ($P_T \ge 0$), while dashed lines indicate accelerating tidal torque ($\mathcal{F}(k_{D,2}^{2,\sigma}) > 0$ and $T_{z,2} > 0$) and negative tidal input power ($P_T < 0$). In the plot of the tidal torque, the magenta vertical dashed line and arrow indicate the position of the present day Earth in the spectrum and the direction of evolution of the semidiurnal tidal frequency, respectively. The red dot designates the corresponding tidal torque.</figcaption></figure>

sufficient for the SPHs. We note that an estimate of $K$ can be obtained by considering the fact that the degrees of the SCHs approximately scale as the inverse of the ocean's angular radius (see Appendix C). Typically, the SCHs exactly correspond to the hemispherical harmonics in the reference case, meaning that their degrees are roughly twice the degrees of the SPHs. Since the number of SPHs needed to represent accurately the sets of oceanic eigenfunctions increases as the size of the basin decreases, configurations dominated by the continental part would require truncation degrees greater than 100 for the SPHs.

## 3.2. Dissipation frequency spectra

By making use of the TRIP computer algebra system (Gastineau & Laskar 2011), we numerically solved the algebraic LTEs established in Sect. 2.7 for 1001 uniformly sampled values of the normalised semidiurnal frequency $\chi = (\Omega - n_s)/\Omega_E$, with $\Omega_E$ being the spin rate of the present day Earth. For simplicity, we let the values of the spin period $P_{\text{rot}} = 2\pi/\Omega$ vary with the tidal frequency while the mean motion of the satellite is fixed. In reality, both $\Omega$ and $n_s$ would vary in the meantime due to exchanges of angular momentum between the planet's spin and the satellite's orbit. However, this approximation is

relevant as long as $n_s \ll \Omega$. The adopted frequency interval of the tidal forcing ranges between 0 (spin-orbit synchronization; $P_{\text{rot}} \approx 27.2$ days) and 4 (fast spin rotation; $P_{\text{rot}} \approx 5.9$ hr). We note that the present Earth-Moon system corresponds to the normalised frequency $\chi \approx 0.963$.

Figure 3 shows the obtained frequency spectra in the reference case for four tidal quantities and their components for each layer (solid part and ocean basin): (i) the input tidal power, (ii) the tidally dissipated power, (iii) the imaginary part of the quadrupolar Love number, and (iv) the tidal torque exerted on the planet about its spin axis. These quantities are evaluated by making use of the expressions given by Eqs. (96), (102), (75), and (87), respectively.

First of all, one observes that the plotted quantities are simply related to each other, the input and dissipated powers differing from the other two quantities by a scaling factor of $\sigma$, as established in the tidal theory (e.g. Ogilvie 2014). Moreover, the tidally dissipated power predicted by the model for the Lunar semidiurnal tide ($M_2$) in the actual Earth-Moon system is close to the ~2.5 TW estimated from altimetric measurements and lunar laser ranging (e.g. Egbert & Ray 2001, 2003). Considering this point, we shall emphasise that we do not expect to recover exactly the estimates obtained from more sophisticated
numerical models owing to the numerous simplifications made in the present approach. Nevertheless, as demonstrated in the next section, the value obtained from measurements could actually be retrieved by slightly tuning the free parameters of the model, though this is not the purpose of this work.

According to the frequency spectra of Fig. 3, the present parameters of the Earth-Moon system (dashed pink line at $\chi \approx 0.963$) places the oceanic tidal response close to a resonance, which corresponds to a maximum of the tidal dissipation rate. This feature of the spectra is in agreement with the increasing dissipation rate inferred from geological data in the past tens of millions of years (see e.g. Farhat et al. 2022a, and references therein) and the predictions of more sophisticated models including realistic bathymetries and land-ocean distributions (Green et al. 2017; Daher et al. 2021). In the low-frequency range, the tidal response of the planet is dominated by the oceanic dynamical tide, which is characterised by resonant peaks. Conversely, in the high-frequency range, the resonances associated with oceanic modes are attenuated by the visco-elastic adjustment of the solid part, and the latter becomes the predominant contributor to the tidally dissipated energy. This change of regime can be simply understood by considering the fact that the solid tidal torque scales as $\propto \sigma^{-\alpha}$ with $\alpha = 0.25$ (see Table 2), while the non-resonant background of the oceanic tidal torque scales as $\propto \sigma^{-3}$ (e.g. Auclair-Desrotour et al. 2019a, Eqs. (47) and (55), respectively).

We note that the tidal input power and tidally dissipated power are not equivalent for a given layer – solid part or oceanic basin – whereas they are strictly equal for the whole planet. This discrepancy is due to the solid-ocean energy transfer allowed by the elasticity of the solid part through ocean loading and self-attraction variations. The energy transfer strongly affects the tidal elongation of the solid part when the oceanic resonances predominate because this elongation is essentially controlled by the ocean loading in this regime. As a consequence, the tidal torque exerted on the solid part can be accelerating instead of being decelerating as observed in the signs of the tidal input power, imaginary part of the Love number, and tidal torque.

Another change of sign can be noticed near the zero-frequency limit ($\sigma \to 0$). This behaviour actually results from the maximum reached by the solid tidal torque for tidal periods close to the Maxwell and Andrade times, which are both much larger than the typical timescales associated with the oceanic waves (e.g. Bolmont et al. 2020). At these timescales, the ocean basin is no longer resonant and the dissipative component of its tidal response thus tends to vanish, which allows the solid tide to be predominant in the very low-frequency range. However, the solid tidal torque also tends to zero for $\sigma \ll \min(\tau_M^{-1}, \tau_A^{-1})$, and $\mathcal{T}_z = 0$ at $\sigma = 0$.

Finally, Fig. 3 highlights the symmetry breaking effect induced by the continent. Whereas the tidal response of a global ocean takes the form of a regular series of harmonics (see e.g. Auclair-Desrotour et al. 2018, Fig. 5), the continent interferes with the forced wave by deviating tidal flows. As a consequence, additional oceanic eigenmodes are excited by the tidal forcing and the energy concentrated into one mode of the global ocean is spread over several modes. As they are each associated with one resonant peak, these new modes alter the general aspect of the spectra by breaking their regularity and by flattening the resonances of the dominating tidal modes. In the next sections, we show quantitatively that the former effect can actually be identified as the typical signature of a supercontinent in the tidal response of an ocean planet.

# 4. Parametric study

As the reference case is now well understood, we can apply the harmonic analysis performed in Sect. 3 to various configurations. The fact that the tidal response is controlled by a relatively small number of physical parameters in our approach (see Table 2) makes it possible to explore the parameter space in an exhaustive way. This is the purpose of the present section.

## 4.1. Preliminary physical analysis

Before all, we shall briefly proceed to a selection of key parameters representing well the parameter space. We recall that the oceanic tidal response is fully determined by the dimensionless control parameters summarised in Table 1. Notwithstanding the deformation of the solid part, which intervenes as a corrective effect through the deformation factors ($\gamma_{D,l}^{\sigma}, \gamma_{T,l}^{\sigma}$), this set of parameters can essentially be reduced to five numbers: three numbers governing the dynamics of oceanic tidal flows ($\tilde{\sigma}, \tilde{\sigma}_G, \tilde{\sigma}_R$) and two numbers defining the geometry of the ocean basin ($\theta_c, \theta_S$). The parameter space of the oceanic tidal response thus simplifies to five dimensions.

Moreover, $\tilde{\sigma} \approx 1$ for the semidiurnal tide as long as $n_s \ll \Omega$, indicating that this number can be roughly considered as invariant over time for planet-satellite systems similar to the Earth-Moon system. Therefore, the parameter space can be limited to the set $(\tilde{\sigma}_G, \tilde{\sigma}_R, \theta_c, \theta_S)$ for the semidiurnal oceanic tide of a rapidly rotating Earth-like planet. For convenience, instead of varying the dimensionless parameters governing the ocean dynamics, we choose to vary the physical parameters they depend on, namely the oceanic depth $H$ for $\tilde{\sigma}_G$, and the Rayleigh drag frequency $\sigma_R$ for $\tilde{\sigma}_R$. The first parameter defines the water volume of the ocean, while the second one characterises the efficiency of the dissipative mechanisms.

Concerning the geometry of the ocean basin, we note that the size of the continent is linked to the oceanic depth through the volume of sea water. Since one expects the water volume to be conserved as the size of the continent evolves, it seems more appropriate to vary $\theta_c$ at constant volume rather than at constant depth. Consequently, when we vary the continental radius, we shall vary the oceanic depth by using the expression of the seawater volume in the thin layer approximation,

$$V = 2\pi H R_p^2 (1 + \cos \theta_c). \qquad (105)$$

Contrary to its size, the position of the continent on the globe can be changed independently of the other parameters.

A set of four additional parameters is brought by the tidal response of the solid part through the Andrade model: the effective shear modulus ($\mu$), the Maxwell time ($\tau_M$), the Andrade time ($\tau_A$), and the dimensionless rheological parameter that accounts for the inelastic component of the solid deformation ($\alpha$). A quick dimensional analysis allows the Maxwell and Andrade times to be eliminated from the list of parameters affecting the regime of the planet's tidal response. These times are indeed much larger than the typical times characterising the ocean dynamics and tidal forcing. Owing to this net separation of scales, varying them over several orders of magnitude would not induce any change of regime in the tidal dynamics. Typically, for $|\sigma| \gg \tau_M^{-1}, \tau_A^{-1}$, the imaginary part of the quadrupolar Love number defined by the Andrade model for the solid part scales as (e.g. Auclair-Desrotour et al. 2019a, Eq. (47))

$$\mathfrak{F}\{k_2^\sigma\} \sim \frac{3}{2} \frac{A_2}{(1+A_2)^2} \Gamma(1+\alpha) \sin\left(\frac{\alpha\pi}{2}\right) (|\sigma|\tau_A)^{-\alpha}, \qquad (106)$$
with $A_2 = 38\pi R_p^2\mu / (3GM_p^2)$ being a constant, and $\Gamma$ the Gamma function (Abramowitz & Stegun 1972, Chapter 6). The above expression entails that, in the limit of $|\sigma| \gg \tau_M^{-1}$, $\tau_A^{-1}$, the solid tidal response is dependent on $\tau_A$, but independent of $\tau_M$. We thus consider only the former in the parametric study.

Similarly as the Maxwell time, the rheological parameter $\alpha$ can almost be considered as a constant parameter owing to the strong constraints provided by laboratory experiments on olivine minerals and ices, which suggest it takes a value between 0.20 and 0.40 (Castillo-Rogez et al. 2011). Due to these constraints, we choose to let it fixed to the value used in the reference case, $\alpha = 0.25$ (see Table 2). Nevertheless, we note from Eq. (106) that varying $\alpha$ would essentially modify the rheological behaviour of the solid part in the high-frequency range by acting on the scaling law $\mathfrak{S}\{k_2^\alpha\} \propto |\sigma|^{-\alpha}$. Finally, the effective shear modulus of the solid part, $\mu$, accounts for the elasticity of the solid part, which is responsible for the solid-ocean coupling. Thus, as discussed for $\tau_A$, we shall consider various values of $\mu$ in the parametric study.

In the above analysis we have examined the role played by the main parameters of the model, and we have selected the six of them that seem to be the most relevant to map the parameter space, ($\hat{\theta}_S$, $\theta_c$, $H$, $\sigma_R$, $\mu$, and $\tau_A$). Other parameters occur in the model, such as the planet mass ($M_p$) or radius ($R_p$), the seawater density ($\rho_w$), or the mass of the perturber ($M_s$). However, they either affect the same dimensionless control numbers as the selected parameters, or they act as obvious scaling factors. For instance, the tidally dissipated power and torque scale as $\propto M_s^2$ since the tidal response is proportional to the forcing gravitational potential – given by Eq. (88) – in the used linear tidal theory. Similarly, we note that the planet mass and radius act on the speed of gravity waves through the planet's surface gravity, which shifts the frequencies of resonant oceanic modes. Such dependences can be characterised analytically from the derivations detailed in Sect. 2.

## 4.2. Calculations

Using the reference case defined in Table 2 (hemispherical equatorial ocean) as a template configuration, we performed calculations for various values of the colatitude of the continental centre ($\hat{\theta}_S$), angular radius of the continent at constant volume ($\theta_c$), oceanic depth ($H$), Rayleigh drag frequency ($\sigma_R$), effective shear modulus of the solid part ($\mu$), and Andrade time ($\tau_A$). In these computations, the parameter values are defined so as to describe the transitions between the asymptotic regimes identified in Table 1. The obtained results are shown by Fig. 4, where the imaginary part of the quadrupolar Love number is plotted as a function of the normalised semidiurnal frequency for variations of each parameter. The values of the reference case (Table 2) are superscripted by*: $\hat{\theta}_S = 90^\circ$, $\theta_c^* = 90^\circ$, $H^* = 4.0$ km, $\sigma_R^* = 10^{-5} \text{ s}^{-1}$, $\mu^* = 25.1189 \text{ GPa}$, and $\tau_A^* = 12897.1 \text{ yr}$.

**Effect of the continental position.** First, we consider the role played by the geometry of the ocean basin. As the colatitude of the continental centre evolves from 0° (polar continent) to 90° (equatorial continent), the frequency spectrum of the Love number rapidly becomes irregular (see Fig. 4, top left panel). For $\hat{\theta}_S = 0^\circ$, the coastline of the hemispherical continent exactly corresponds to the planet's equator. As a consequence, it does not interfere with the latitudinal component of tidal flows, which is already zero at the equator in the absence of a continent. This

explains why we recover approximately the spectrum obtained for the global ocean (Auclair-Desrotour et al. 2018) up to a factor of 2, the noticed slight discrepancy being due to the contribution of the solid part. This spectrum is composed of a regular series of harmonics that results from the coupling – caused by Coriolis forces – between the oceanic normal modes and the quadrupolar tidal potential. For a non-rotating planet, only one resonance would be observed in that configuration. A slight inclination of the continental centre with respect to the pole (Fig. 4, top left panel, $\hat{\theta}_S = 10^\circ$) appears to be sufficient to break the axis symmetry and to alter significantly the regularity of the spectrum. The harmonics of the polar case are split into more modes, which induces intermediate relative maxima and decreases the dissipation rate reached while crossing the main resonances.

**Effect of the continental size.** Such an effect can also be observed as the angular radius of the continent varies from 0° (global ocean) to 90° (hemispherical ocean) at constant water volume (Fig. 4, top right panel). As the size of the continent increases, the peaks of the global ocean get progressively altered by the interferences of coastlines with tidal flows. However, we remark that this effect really becomes significant only for $\theta_c \ge 30^\circ$, which suggests that continents smaller than South America hardly alter the tidal response of the global ocean even if they are located at the equator. Moreover, the spectra highlight the fact that oceanic resonances are shifted to the right as the size of the continent increases, which also means resonances become of smaller amplitudes, so less dissipative. This effect is directly related to the ocean's geometry. Acoustically, the ocean basin is analogous to a vibrating drumhead, with its size and depth standing for the drumhead's diameter and tension, respectively. On one hand, the eigenfrequencies of the oceanic normal modes increase as the size of the basin decreases since the latter corresponds to the wavelength of the lowest tidal modes. On the other hand, the speed of gravity waves scales as $\propto \sqrt{H}$, meaning that it increases as the size of the ocean basin decreases at constant volume.

**Effect of the oceanic depth.** Consistently with the above drumhead analogy, increasing the oceanic depth essentially induces a dilation of the frequency spectrum without changing its global aspect (Fig. 4, middle left panel). This is related to the way the resonances are controlled by the oceanic depth. Both the heights and frequencies of resonant peaks scale as $\propto \sqrt{H}$ (e.g. Auclair-Desrotour et al. 2019a, Eq. (55)), which explains why their relative positions and amplitudes are conserved as $H$ varies.

**Effect of dissipative mechanisms.** Similarly, we analyse the evolution of the frequency spectrum with the Rayleigh drag frequency (Fig. 4, middle right panel) using the closed-form solutions provided by the analytical theory (e.g. Auclair-Desrotour et al. 2018). The plot highlights the transition from the quasi-adiabatic regime ($\sigma_R = 10^{-6} \text{ s}^{-1}$) to the frictional regime ($\sigma_R = 10^{-3} \text{ s}^{-1}$). In the former, the resonances of oceanic modes are weakly damped, while they are completely annihilated in the latter.

The two asymptotic regimes can be considered separately. In the quasi-adiabatic regime, the frequencies of resonances are almost independent of the Rayleigh drag frequency, whereas their heights scale as $\propto \sigma_R^{-1}$, their widths as $\propto \sigma_R$, and the non-resonant background of the oceanic dissipation rate as
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Frequency spectra of the imaginary part of the quadrupolar Love number associated with the semidiurnal tide, $k_{D,2}^{2\sigma}$, for the considered reference case and various values of the model key parameters. **Top:** variation of the colatitude of the continental centre (left panel) and continental angular radius at constant seawater volume (right panel). **Middle:** variation of the oceanic depth (left panel) and Rayleigh drag frequency (right panel). **Bottom:** variation of the effective shear modulus (left panel) and Andrade time (right panel) of the solid part. In each panel, the imaginary part of the Love number is plotted in logarithmic scale (vertical axis) as a function of the normalised semidiurnal frequency $\chi = (\Omega - n_s)/\Omega_E$ (horizontal axis), with $\Omega_E$ being the spin angular velocity of the actual Earth. The solid black line designates the reference case defined by Table 2: $\dot{\theta}_S^* = 90^\circ$, $\theta_c^* = 90^\circ$, $H^* = 4.0$ km, $\sigma_R^* = 10^{-5} \text{ s}^{-1}$, $\mu^* = 25.1189 \text{ GPa}$, and $\tau_A^* = 12897.1 \text{ yr}$.</figcaption></figure>

$\propto \sigma_R$ (Auclair-Desrotour et al. 2015, 2018, 2019a). As a consequence, the number of visible peaks increases as $\sigma_R$ decreases. The new peaks appear in the high-frequency range since they are associated with oceanic harmonics of high degrees. In the zero-drag limit ($\sigma_R \to 0$), the number of resonant peaks exactly corresponds to the number of harmonics with eigenfrequencies less than the upper bound of the considered frequency interval. In the frictional regime, the spectrum is regular and only one maximum is visible, for a frequency

scaling as $\propto \sigma_R^{-1}$. This behaviour is described by a closed-form solution that will be detailed in a forthcoming article. Considering the estimate of ~$10^{-5} \text{ s}^{-1}$ found for the Earth's Rayleigh drag frequency from model adjustments to the current dissipation rate (Webb 1980; Farhat et al. 2022a), we note that this value is just slightly less than the critical value marking the transition between the quasi-adiabatic and frictional regimes in Fig. 4 (middle right panel), namely $\sigma_R \approx 3.2 \times 10^{-5} \text{ s}^{-1}$ (see also Sect. 5).
<table><caption>**Table 3.** Evolution of the frequency spectrum of the quadrupolar tidal Love number as one parameter increases while the values of other parameters are fixed.</caption><thead><tr><th>Altered features of spectra</th><th>θ̂<sub>S</sub></th><th>θ<sub>c</sub></th><th>H</th><th>σ<sub>R</sub></th><th>μ</th><th>τ<sub>A</sub></th></tr></thead><tbody><tr><td>Frequencies of resonant peaks</td><td>-</td><td>↗</td><td>↗</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Heights of resonant peaks</td><td>-</td><td>-</td><td>↗</td><td>↘</td><td>↗</td><td>-</td></tr><tr><td>Widths of resonant peaks</td><td>-</td><td>-</td><td>-</td><td>↗</td><td>-</td><td>-</td></tr><tr><td>Non-resonant background</td><td>-</td><td>-</td><td>-</td><td>↗</td><td>↗</td><td>↘</td></tr><tr><td>Irregularity of spectra</td><td>↗</td><td>↗</td><td>-</td><td>↘</td><td>-</td><td>-</td></tr></tbody></table>

**Notes.** The table is limited to the quasi-adiabatic regime, where oceanic resonances are not completely damped by dissipative mechanisms. In the frictional regime, there is no oceanic resonances, and the quadrupolar Love number varies smoothly with the forcing tidal frequency.

**Effect of the solid elasticity and anelasticity.** We finally examine the role played by the parameters of the solid part (Fig. 4, bottom panels). As discussed above, the elasticity of the solid part acts to attenuate the resonances of oceanic tidal modes by elastic adjustment. The load caused by a local water mass surplus lowers the oceanic floor, which decreases the gravitational potential of the planet's tidal response. If the solid part is rigid, this tidal response corresponds to the oceanic tide, as observed for $\mu = 300$ GPa. The Andrade time only affects the non-resonant background of the spectrum through the scaling factor $\tau_A^{-\alpha}$ of the solid dissipation rate, as shown by Eq. (106). As $\tau_A$ increases, the non-resonant background decreases. This effect is visible in the high-frequency range, in particular. Table 3 summarises the effects of all the studied parameters on specific features of the frequency spectra.

**Ignored effects and model limitations.** As the spectral features detailed above for tidal dissipation have been obtained from the simplified theoretical framework of Sect. 2, we shall end this section by discussing the main ignored effects that may alter them. First, although the single continent configuration appears as a convenient setup to introduce the anisotropic effect induced by continentality, realistic land-ocean distributions are more complex. Consequently, the oceanic tidal response cannot be simply related to the geometry of the ocean basin in the general case, as highlighted by numerical studies performed for Earth and exoplanetary land-ocean distributions (Egbert & Ray 2001, 2003; Green et al. 2017; Blackledge et al. 2020). However, since the oceanic tidal waves are planetary scale, one may reasonably expect that they shall mainly interfere with large scale continental features. The predominant modes that shape the tidal dissipation spectrum of the planet are thus related to the largest length scales of the continental distribution.

That said, as shown by Green et al. (2017), the smoothness of the continental geometry tends to yield overestimates of the tidally dissipated energy. The coastline fractality thus both acts to attenuate and broaden the resonances. Similarly, we remark that the eigenfrequencies of the oceanic tidal modes vary with the oceanic depth, given that the phase speed of gravity waves scales as $\propto H^{1/2}$. As a consequence, planetary scale variations of the oceanic depth presumably act to attenuate and broaden the resonant peaks. Complex bathymetries are associated with spatially dependent dissipation rates, the dissipation coefficients being strongly determined by small-scale topographic features (e.g. Egbert & Ray 2001; Carter et al. 2008; Arbic et al. 2010; Green & Huber 2013; Merdith et al. 2021). It is also noteworthy that predominant dissipative mechanisms such as turbulent bottom drag scale non-linearly with the speed of tidal flows, which precludes resonances to reach very high amplitudes. Analogously, the stabilising effect of an hypothetic overlying ice shell

would damp the resonant peaks, thus reducing tidal dissipation, as observed in the case of icy satellites (e.g. Beuthe 2016; Matsuyama et al. 2018).

## 5. A metric for the signature of supercontinents

Until now, we have been investigating the direct problem that consists in generating frequency spectra of the tidal quantities from sets of values assigned to the model parameters. Particularly, in Sect. 4, we made an attempt to infer the tidal response of planets hosting ocean basins from knowledge of their geometry. The relationships between spectral and geometric features actually correspond to a very broad field of mathematical problems known as spectral geometry. Through popular case studies of this field, such as Kac's drum (Kac 1966), it was shown that the geometry of a vibrating system can be predicted to a certain extent if its spectral features are known. This is the so-called inverse problem that we aim to address in the present section by investigating the extent to which one can infer the geometry of the ocean basin from the spectral features of the planet's tidal response. Particularly, we define further a metric that quantifies the specific signature of a supercontinent in the evolution of the tidal torque with the tidal frequency. This frequency dependence of the tidal torque is straightforwardly connected to the orbital evolution of the satellite or the rotational evolution of the planet, as any crossed resonance induces a rapid variation of the orbital parameters over a short time period (see e.g. Auclair-Desrotour et al. 2014; Farhat et al. 2022a).

The parametric study of Sect. 4 shows evidence of the numerous degeneracies characterising the dependences of the frequency spectra on the system's parameters (see Table 3). Several features – such as the height of resonant peaks or the non-resonant background – appear to be sensitive to three key parameters at the same time. Particularly, these results highlight the predominant role played by the Rayleigh drag frequency, which affects almost all of the spectral features listed in Table 3. Nevertheless, we also remark that the only parameters having an impact on the regularity of the spectra – notwithstanding the effect of the Rayleigh drag frequency – are those defining the geometry of the ocean basin, namely the colatitude of the continental centre and the angular radius of the continent. This argues for the existence of spectral quantities sensitive to these parameters only. Such quantities are metrics of the geometry. Ideally, they shall vary with $\theta_S$, $\theta_c$, and $\sigma_R$, while being insensitive to the other parameters.

In Fig. 4, the irregularity of spectra seems to be mainly characterised by a significant dispersion of the frequency intervals separating consecutive relative maxima of the tidal torque, denoted by $\Delta\sigma_k \equiv \sigma_{k+1} - \sigma_k$ in the following, with $\sigma_k$ being the frequency of the k-th maximum. In order to make these
<figure><img src="image_6.png" /><figcaption>Fig. 5. CDF of the frequency intervals separating two consecutive relative maxima in the spectra of Fig. 4. Top: variation of the colatitude of the continental centre (left panel) and continental angular radius at constant seawater volume (right panel). Middle: variation of the oceanic depth (left panel) and Rayleigh drag frequency (right panel). Bottom: variation of the effective shear modulus (left panel) and Andrade time (right panel) of the solid part. For each spectrum, the frequency intervals are normalised by their average, and the CDF is normalised. The solid black line designates the reference case defined by Table 2: $\hat{\theta}_S^* = 90^\circ$, $\theta_c^* = 90^\circ$, $H^* = 4.0$ km, $\sigma_R^* = 10^{-5} \text{ s}^{-1}$, $\mu^* = 25.1189 \text{ GPa}$, and $\tau_A^* = 12897.1 \text{ yr}$.</figcaption></figure>

frequency intervals insensitive to the dilation-contraction effect induced by oceanic depth variations (see Fig. 4, middle left panel), they are normalised by their average, given by

$$ \overline{\Delta\sigma} = \frac{1}{N} \sum_{k=1}^{N} \Delta\sigma_k, \qquad (107) $$

where N designates the total number of intervals. We thus denote by $\hat{\Delta\sigma}_k \equiv \Delta\sigma_k / \overline{\Delta\sigma}$ the normalised intervals. Finally, sorting the intervals in ascending order and using the subscript j such that $\Delta\sigma_j < \Delta\sigma_{j+1}$ for $j = 1, \dots, N$, we introduce the cumulative distribution function (CDF) of the normalised frequency intervals,

$$ \text{CDF}(\hat{\Delta\sigma}) = \frac{1}{N} \underset{\Delta\hat{\sigma}_j \le \Delta\sigma}{\operatorname{argmax}} \{\hat{\Delta\sigma}_j\}. \qquad (108) $$

Figure 5 shows the CDFs of the frequency intervals separating two consecutive maxima for the spectra plotted in Fig. 4. In these plots, the quasi-uniformly spaced peaks of the global and polar hemispherical oceans are characterised by a sharp step in the CDF, which abruptly switches from 0 to 1. Conversely, the symmetry breaking effect of continental shelves induces a smooth transition. As noted qualitatively with the spectra of Fig. 4, the slope of the CDF essentially varies with the position of the continent on the globe, its size, and the dissipative
time scale. This leads us to consider, as an appropriate metric for the continentality effect, the normalised standard deviation of the frequency intervals separating consecutive maxima of the tidal torque,

$$ \hat{\sigma} = \sqrt{\frac{1}{N} \sum_{k=1}^{N} \left( \frac{\Delta\sigma_k}{\Delta\sigma} - 1 \right)^2} \quad (109) $$

The thus defined normalised standard deviation is such that $0 \le \hat{\sigma} \le 1$, where $\hat{\sigma} = 0$ corresponds to uniformly spaced relative maxima or to $N = 1$ (two maxima), and $\hat{\sigma} = 1$ to an extreme dispersion of the resonant peaks. By convention, we set $\hat{\sigma}$ to zero when one unique maximum is observed.

By reproducing the methodology detailed in Sect. 4, we computed the tidal response of the reference planet defined in Table 2 for various values of the colatitude of the continental centre, angular radius of the continent at constant seawater volume, oceanic depth, Rayleigh drag frequency, effective shear modulus of the solid part, and Andrade time, each parameter being modified independently of the others. In all cases, we calculated the values of $\hat{\sigma}$ from the frequency spectra obtained for the imaginary part of the quadrupolar Love number. These values are plotted in Fig. 6.

The quantity $\hat{\sigma}$ appears to be almost insensitive to variations of $H$, $\mu$, and $\tau_A$, whereas it varies over several orders of magnitude when $\hat{\theta}_S$, $\theta_c$, or $\sigma_R$ are modified. Moreover, we recover quantitatively the tendencies identified in the parametric study (Table 3), namely $\hat{\sigma}$ really accounts for the role played by the geometry of the ocean basin and the dissipative timescale in the planet's tidal response. Thus, the normalised standard deviation defined by Eq. (109) can be considered as a spectral metric of the planet's geometric features, both on the global scale pertaining to the location or spread of land, if any, along with the local topographical scale dictating the dissipative timescale in the ocean.

Interestingly, the dependence of $\hat{\sigma}$ on $\hat{\theta}_S$, $\theta_c$, or $\sigma_R$, is not even. Instead, threshold effects may be observed, similar to those of Fig. 5. First, the evolution of the metric with $\hat{\theta}_S$ (Fig. 6, top left panel) illustrates the fact that a small inclination of the continent on the globe with respect to the pole is sufficient to break the axi-symmetry. By switching from $\hat{\theta}_S = 0^\circ$ (polar continent) to $\hat{\theta}_S = 10^\circ$, the metric skyrockets to reach a plateau around $\hat{\sigma} \sim 0.4$. This plateau suggests that changing the position of the continental centre from $10^\circ$ to $90^\circ$ does not affect the regularity of the frequency spectrum. Basically, one observes here a saturation regime where the signature of the supercontinent through the chosen metric is the same regardless of its position on the globe.

A similar behaviour is observed as the angular radius of the continent grows from $0^\circ$ (global ocean) to $90^\circ$ (hemispherical ocean) at constant water volume (Fig. 6, top right panel). While the metric does not evolve much with $\theta_c$ in general, it increases drastically around $\theta_c \approx 35^\circ$, which is the critical size identified from the plotted spectra in the parametric study. This value is consistent with the highest value of $\theta_c$ for which the frequency spectrum still exhibits a regular series of resonant peaks that resembles the tidal response of the global ocean (Fig. 4, top right panel). Beyond $35^\circ$, the spectrum becomes complex and irregular due to the strong interference of the continent with tidal flows. We note that the metric remains stable for $\theta_c \ge 40^\circ$, which suggests that the regularity of the frequency spectrum is insensitive to the size of the continent beyond this value. This highlights another saturation regime, similar to that observed for the position of the continent, in which the signature of the supercontinent

is the same regardless of its size. It is noteworthy that the way the continent interferes with the oceanic tidal flows is related to its size with respect to the typical horizontal wavelengths of the predominating tidal modes, which are planetary scale. As a consequence, the critical angular radius marking the regime transition for the studied Earth-sized planet would be similar for a super-Earth or a smaller planet.

Finally, the metric decays monotonically as the Rayleigh drag frequency increases, with a sharp variation occurring around $\sigma_R \sim 3.2 \times 10^{-5} \text{ s}^{-1}$ (Fig. 6, middle right panel). This critical value marks the transition between the quasi-adiabatic and frictional regimes (see Table 1). In the latter, the resonances associated with oceanic modes are damped by the strong friction of tidal flows with the oceanic floor, which makes the tidal torque evolve smoothly with the tidal frequency, as shown by Fig. 4 (middle right panel). Therefore, the regularity of the spectrum conveys no information about the geometry of the ocean basin in this regime. Conversely, in the quasi-adiabatic regime, $\hat{\sigma}$ increases as $\sigma_R$ decays due to the growing number of visible resonant peaks.

## 6. Conclusions

In this paper we have examined the linear response of an ocean planet hosting a supercontinent to gravitational tidal forcing. This problem may be regarded as a simplified model of the tides raised on Earth by the Moon and their interactions with the land-ocean distribution, which predominantly drives the evolution of the Earth-Moon system over long time scales (e.g. Bills & Ray 1999; Daher et al. 2021; Tyler 2021; Farhat et al. 2022a). The theory developed in the present study expands on earlier works that investigated the tidal response of hemispherical or global oceans (Longuet-Higgins 1968; Longuet-Higgins & Pond 1970; Webb 1980, 1982; Auclair-Desrotour et al. 2018; Farhat et al. 2022a). In this approach, the continent is represented by a spherical cap of arbitrary size and position, and the LTEs are written in the shallow water approximation for an oceanic layer of uniform depth. These simplifications allow the problem to be formulated analytically in terms of explicitly defined oceanic eigenmodes, which provides a deep insight into the physics governing the planet's tidal response. Additionally, such a formalism appears to be well suited to the harmonic analysis of the tidal dynamics since it relies on a small set of key parameters.

As a first step, we established the equations describing the dynamics of the forced tidal flows. These equations include the gravitational and surface couplings that result from the deformation of the solid part generated by the ocean loading and the tidal forces exerted by the perturber. The visco-elastic response of the solid part to tidal forcings is described in a generic way by gravitational and load Love numbers. We computed it in practice by making use of the closed-form solutions provided by the Andrade model for a homogeneous body, following the prescriptions obtained from 1D models taking the planet's internal structure into account. The LTEs were expressed as an algebraic system describing the evolution of coupled eigenmodes in the frequency domain. This system is controlled by a small set of clearly identified dimensionless numbers that characterise the regime of oceanic tidal waves and the geometry of the ocean basin. We finally used the harmonic expansion of tidal flows to express the corresponding input tidal power, tidally dissipated power, tidal Love numbers, and tidal torque exerted about the spin axis.

As a second step, we applied the theory to a reference case, which is essentially defined as a simplified Earth-Moon system.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Evolution of the normalised standard deviation of the frequency interval separating two relative maxima with the parameters of Fig. 4. The metric, given by Eq. (109), is plotted in logarithmic scale (vertical axis) as a function of the parameters (horizontal axis). **Top:** Dependence on the colatitude of the continental centre (**left panel**) and angular radius of the continent (**right panel**). **Middle:** Dependence on the oceanic depth (**left panel**) and Rayleigh drag frequency (**right panel**). **Bottom:** Dependence on the effective shear modulus (**left panel**) and Andrade time (**right panel**) of the solid part. Blue dots indicate the values obtained for various values of the parameters, while black squares designate the reference case defined in Table 2: $\dot{\theta}_S = 90^\circ$, $\theta_c^* = 90^\circ$, $H^* = 4.0$ km, $\sigma_R^* = 10^{-5}$ s$^{-1}$, $\mu^* = 25.1189$ GPa, and $\tau_A^* = 12$ 897.1 yr.</figcaption></figure>

In this configuration, the continent is hemispherical and located at the equator. By computing the evolution of the aforementioned tidal quantities with the forcing tidal frequency, we recovered the solutions obtained in Farhat et al. (2022a). We showed that the planet's tidal dissipation rate is dominated by the contribution of the ocean, while the solid part determines it mainly in the high-frequency range, that is for a fast rotating planet. Moreover, the obtained results highlight the symmetry breaking effect of the supercontinent, which splits the modes of the global ocean into more modes by interfering with tidal flows. As a result, the

spectra of tidal quantities are flattened and they exhibit additional relative maxima corresponding to the resonances of the basin modes.

As a third step, we carried out a parametric study for a set of six key parameters, by using the previously defined reference case as a template configuration. The obtained results reveal the degeneracies affecting the evolution of the main spectral features with the selected parameters. However, they also indicate that the specific signature of the geometry can be unravelled to a certain extent by characterising the erraticity of the spectra. As the
position or angular radius of the continent increase, the frequency spectrum of the tidal torque becomes irregular, while the other parameters only alter the frequencies of the resonant peaks, their widths, their heights, and the non-resonant background. The only exception to this statement is the Rayleigh drag frequency, which affects the regularity of frequency spectra as well.

The above qualitative results were refined in a quantitative way as a last step. By following an inverse problem approach inspired from spectral geometry, we showed that some geometric features of the ocean basin may be inferred from the spectral irregularity, which is quantified by the normalised standard deviation of the frequency intervals separating consecutive relative maxima in the frequency spectrum of the tidal Love number. The sensitivity of this quantity to the geometrical parameters is characterised by threshold effects. As regards the position of the continent on the globe, the metric only allows the poles to be distinguished from other locations. Similarly, it seems to indicate that the signatures of continents with angular radii larger than 40° are indistinguishable from each other.

Nevertheless, the harmonic analysis also suggests that a continent similar to South America (~30°-angular radius) or smaller would not alter, qualitatively, the oceanic tidal response even if it is located at the equator, where the interference of the coastlines with tidal flows tends to be maximal. In such configurations, the tidally dissipated energy predicted by the harmonic analysis is approximately the same as in the absence of continent. Finally the metric appears to be sensitive to the transition between the quasi-adiabatic and frictional regimes, which occurs for $\sigma_R \approx 3.2 \times 10^{-5} \text{ s}^{-1}$. This critical value is slightly greater than that estimated for Earth (~$10^{-6}-10^{-5} \text{ s}^{-1}$, e.g. Garrett & Munk 1971; Webb 1973; Tyler 2021; Farhat et al. 2022a), which is consistent with the frequency-resonant behaviour of the Earth's tidal dissipation rate inferred from Lunar laser ranging and geological data (e.g. Bills & Ray 1999; Tyler 2021).

Both the spectral method and the metric for continentality effects elaborated in the present work are highly relevant to studies dealing with the long-term tidal evolution of bodies with surface oceans, such as extrasolar rocky planets orbiting in the habitable zone of their host star. Particularly, the studied spectral irregularity of the tidal torque has direct implications on the history of planet-satellite systems. Tidal resonances cause irregular evolutions of the satellite's orbit and planet's rotation, thus leading to abrupt LOD or obliquity variations that might result in significant climatic effects. We will examine specifically the consequences of continentality on the long-term evolution of the Earth-Moon system in a separate dedicated study by using the theory established in the present work.

As this theory forms a basis for more sophisticated descriptions of continentality effects, we will consider refining the geometry in the future in order to account for the role played by the fractality of non-smooth coastlines, the complex bathymetry of the ocean basin, or the spatial dependence of energy dissipation. Also, while we have characterised the signature of one unique supercontinent in the planet's tidal response, the combined effects of several continents still needs to be elucidated. All these questions require further work improving meanwhile the theory and the quality of observational constraints on the past history of the Earth-Moon system. We will address them through several forthcoming studies.

*Acknowledgements.* The authors thank the anonymous referee for his/her insightful comments that significantly enhanced the manuscript. 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 research has made use of NASA's Astrophysics Data System.

## References

Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions, National Bureau of Standards: Applied Mathematics Series, 55

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

Arbic, B. K., Wallcraft, A. J., & Metzger, E. J. 2010, Ocean Model., 32, 175

Arfken, G. B., & Weber, H. J. 2005, Mathematical Methods for Physicists, 6th edn. (Elsevier Academic Press)

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

Auclair-Desrotour, P., Mathis, S., & Le Poncin-Lafitte, C. 2015, A&A, 581, A118

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

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

Auclair-Desrotour, P., Leconte, J., & Mergny, C. 2019b, A&A, 624, A17

Beuthe, M. 2016, Icarus, 280, 278

Bills, B. G., & Ray, R. D. 1999, Geophys. Res. Lett., 26, 3045

Biot, M. A. 1954, J. Appl. Phys., 25, 1385

Blackledge, B. W., Green, J. A. M., Barnes, R., & Way, M. J. 2020, Geophys. Res. Lett., 47, e85746

Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728

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

Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121

Cadieux, C., Doyon, R., Plotnykov, M., et al. 2022, AJ, 164, 96

Carr, M. H., & Head, J. W. 2003, J. Geophys. Res. (Planets), 108, 5042

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

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

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

Charbonneau, D., Berta, Z. K., Irwin, J., et al. 2009, Nature, 462, 891

Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767

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

Cowling, T. G. 1941, MNRAS, 101, 367

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

Dohm, J. M., Baker, V. R., Boynton, W. V., et al. 2009, Planet. Space Sci., 57, 664

Eakins, B. W., & Sharman, G. F. 2010, NOAA National Geophysical Data Center, Boulder, CO

Efroimsky, M. 2012, ApJ, 746, 150

Egbert, G. D., & Ray, R. D. 2001, J. Geophys. Res., 106, 22

Egbert, G. D., & Ray, R. D. 2003, Geophys. Res. Lett., 30, 1907

Farhat, M., Auclair-Desrotour, P., Boué, G., & Laskar, J. 2022a, A&A, 665, A1

Farhat, M., Laskar, J., & Boué, G. 2022b, J. Geophys. Res. (Solid Earth), 127, e2021JB023323

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

Garrett, C. J. R., & Munk, W. H. 1971, Deep Sea Res. A, 18, 493

Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194

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

Gerkema, T., & Zimmerman, J. 2008, Lecture Notes, Royal NIOZ, Texel

Gordon, C., Webb, D., & Wolpert, S. 1992, Invent. Math., 110, 1

Green, J. A. M., & Huber, M. 2013, Geophys. Res. Lett., 40, 2707

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

Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68

Haines, G. V. 1985, J. Geophys. Res., 90, 2583

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

Hendershott, M. C. 1972, Geophys. J., 29, 389

Hut, P. 1981, A&A, 99, 126

Hwang, C., & Chen, S.-K. 1997, Geophys. J. Int., 129, 450

Johansson, F. 2016, ArXiv e-prints [arXiv:1606.06977]

Kac, M. 1966, Am. Math. Monthly, 73, 1

Kamata, S., Matsuyama, I., & Nimmo, F. 2015, J. Geophys. Res. (Planets), 120, 1528

Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108

Kivelson, M. G., Khurana, K. K., Russell, C. T., et al. 2000, Science, 289, 1340

Kodaira, T., Thompson, K. R., & Bernier, N. B. 2016, J. Geophys. Res. (Oceans), 121, 6159
Kopparapu, R. K., Wolf, E. T., Arney, G., et al. 2017, ApJ, 845, 5

Lambeck, K. 1977, Philos. Trans. Roy. Soc. Lond. Ser. A, 287, 545

Laplace, P. S. 1798, Traité de mécanique céleste (Duprat J. B. M.), Tome 5, Livre XIII, 145

Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632

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

Léger, A., Selsis, F., Sotin, C., et al. 2004, Icarus, 169, 499

Le Provost, C., Genco, M. L., Lyard, F., Vincent, P., & Canceil, P. 1994, J. Geophys. Res., 99, 24,777

Le Provost, C., Lyard, F., Molines, J. M., Genco, M. L., & Rabilloud, F. 1998, J. Geophys. Res., 103, 5513

Levrard, B., Winisdoerffer, C., & Chabrier, G. 2009, ApJ, 692, L9

Lillo-Box, J., Figueira, P., Leleu, A., et al. 2020, A&A, 642, A121

Lindzen, R. S., & Chapman, S. 1969, Space Sci. Rev., 10, 3

Longuet-Higgins, M. S. 1968, Philos. Trans. Roy. Soc. Lond. Ser. A, 262, 511

Longuet-Higgins, M. S., & Pond, G. S. 1970, Philos. Trans. Roy. Soc. Lond. Ser. A, 266, 193

MacDonald, G. J. F. 1964, Rev. Geophys. Space Phys., 2, 467

Makarov, V. V., & Efroimsky, M. 2013, ApJ, 764, 27

Mamajek, E. E., Prsa, A., Torres, G., et al. 2015, ArXiv e-prints [arXiv:1510.07674]

Matsuyama, I. 2014, Icarus, 242, 11

Matsuyama, I., Beuthe, M., Hay, H. C. F. C., Nimmo, F., & Kamata, S. 2018, Icarus, 312, 208

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

Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics (McGraw-Hill Book Comp.)

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

Munk, W. H., & MacDonald, G. J. F. 1960, J. Geophys. Res., 65, 2169

Ogilvie, G. I. 2013, MNRAS, 429, 613

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

Payne, M. J., & Lodato, G. 2007, MNRAS, 381, 1597

Peale, S. J., Cassen, P., & Reynolds, R. T. 1979, Science, 203, 892

Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge: Cambridge University Press)

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes, 3rd edn.: The Art of Scientific Computing (Cambridge University Press)

Proudman, J. 1920, Proc. Lond. Math. Soc., 2, 51

Provost, C. L., & Lyard, F. 1997, Progr. Oceanogr., 40, 37

Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, ApJ, 669, 606

Remus, F., Mathis, S., Zahn, J.-P., & Lainey, V. 2012, A&A, 541, A165

Scheller, E. L., Ehlmann, B. L., Hu, R., Adams, D. J., & Yung, Y. L. 2021, Science, 372, 56

Schwiderski, E. W. 1980, Rev. Geophys. Space Phys., 18, 243

Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337

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

Takeuchi, H., & Saito, M. 1972, Methods Computat. Phys., 11, 217

Thébault, E., Schott, J. J., Mandea, M., & Hoffbeck, J. P. 2004, Geophys. J. Int., 159, 83

Thébault, E., Schott, J. J., & Mandea, M. 2006, J. Geophys. Res. (Solid Earth), 111, B01102

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

Torta, J. M. 2019, Surv. Geophys., 41, 201

Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86

Tyler, R. H. 2008, Nature, 456, 770

Tyler, R. 2011, Icarus, 211, 770

Tyler, R. 2014, Icarus, 243, 358

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

Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: University of Tokyo Press)

Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics, 770

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

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

Webb, D. J. 1973, Deep Sea Res. A, 20, 847

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

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

Wunsch, C., Haidvogel, D. B., Iskandarani, M., & Hughes, R. 1997, Progr. Oceanogr., 40, 81

Zahn, J. P. 1966, Ann. Astrophys., 29, 313
Appendix A: Nomenclature

The notations introduced in the main text are listed below in order of appearance.

LOD Length of the day p. 2

LTEs Laplace's Tidal Equations p. 2

$R_p$ Planet radius p. 2

$H$ Oceanic depth p. 2

$\theta_c$ Angular radius of the supercontinent p. 2

$\theta_0$ Angular radius of the ocean basin p. 2

$\mathcal{R}$ Geocentric frame of reference rotating with the planet p. 2

$\mathbf{e}_X$ Cartesian unit vector associated with $\mathcal{R}$ p. 2

$\mathbf{e}_Y$ Cartesian unit vector associated with $\mathcal{R}$ p. 2

$\mathbf{e}_Z$ Cartesian unit vector associated with $\mathcal{R}$ p. 2

$\hat{\mathbf{r}}$ Radial coordinate in $\mathcal{R}$ p. 2

$\hat{\theta}$ Colatitude in $\mathcal{R}$ p. 2

$\hat{\varphi}$ Longitude in $\mathcal{R}$ p. 2

$\hat{\theta}_S$ Colatitude of the continental centre in $\mathcal{R}$ p. 3

$\hat{\varphi}_S$ Longitude of the continental centre in $\mathcal{R}$ p. 3

$\hat{\theta}_{oc}$ Colatitude of the oceanic centre in $\mathcal{R}$ p. 3

$\hat{\varphi}_{oc}$ Longitude of the oceanic centre in $\mathcal{R}$ p. 3

$\mathcal{O}$ Planet's centre of mass p. 3

$\mathcal{R}_{oc}$ Frame of reference rotating with the planet and oriented by the centre of the ocean basin p. 3

$\mathbf{e}_x$ Cartesian unit vector associated with $\mathcal{R}_{oc}$ p. 3

$\mathbf{e}_y$ Cartesian unit vector associated with $\mathcal{R}_{oc}$ p. 3

$\mathbf{e}_z$ Cartesian unit vector associated with $\mathcal{R}_{oc}$ p. 3

$\mathbf{r}$ Radial coordinate in $\mathcal{R}_{oc}$ p. 3

$\theta$ Colatitude in $\mathcal{R}_{oc}$ p. 3

$\varphi$ Longitude in $\mathcal{R}_{oc}$ p. 3

$\mathbf{e}_r$ Radial unit vector associated with $\mathcal{R}_{oc}$ p. 3

$\mathbf{e}_\theta$ Colatitudinal unit vector associated with $\mathcal{R}_{oc}$ p. 3

$\mathbf{e}_\varphi$ Longitudinal unit vector associated with $\mathcal{R}_{oc}$ p. 3

$U_T$ Tidal gravitational potential generated by the perturber p. 3

$G$ Universal gravitational constant p. 3

$M_s$ Mass of the perturber p. 3

$\mathbf{r}_s$ Position vector of the perturber with respect to $\mathcal{O}$ p. 3

$r_s$ Planet-perturber distance p. 3

$\equiv$ Symbol used for definitions in equations p. 3

$\mathbf{F}$ Tidal force per unit mass induced by $U_T$ p. 3

$t$ Time p. 3

$g$ Surface gravity at rest p. 3

$\sigma_R$ Rayleigh drag frequency p. 3

$\mathbf{V}$ Horizontal velocity vector p. 3

$\xi$ Horizontal displacement vector p. 3

$\mathbf{f}$ Coriolis parameter p. 3

$\zeta$ Vertical displacement of the ocean's surface with respect to the oceanic floor p. 3

$\zeta_{eq}$ Equilibrium displacement corresponding to the equi potential surface p. 3

$\Gamma_D, \Gamma_T$ Solid deformation operators p. 3

$\Omega$ Planet's spin angular velocity p. 3

$\mathbf{CA}$ Cowling approximation p. 4

$t_0$ Reference time for normalisation p. 3

$V_0$ Reference velocity for normalisation p. 3

$\tilde{t}$ Normalised time p. 4

$\tilde{f}$ Normalised Coriolis parameter p. 4

$\tilde{\nabla}$ Normalised gradient operator p. 4

$\tilde{\xi}$ Normalised horizontal displacement vector p. 4

$\tilde{\zeta}$ Normalised oceanic surface elevation p. 4

$\tilde{\zeta}_{eq}$ Normalised equilibrium surface elevation p. 4

$\tilde{\Xi}$ Normalised forcing term p. 4

$\tilde{\sigma}_G$ Normalised Rossby deformation length p. 4

$\tilde{\sigma}_R$ Normalised friction parameter p. 4

$\sigma$ Tidal frequency p. 4

$\tilde{\sigma}$ Rossby number of forced waves p. 4

$\Phi$ Divergent displacement potential p. 4

$\Psi$ Rotational displacement streamfunction p. 4

$\mathbf{n}$ Outward pointing unit vector defining the normal to the coast p. 4

$dS$ Infinitesimal surface element p. 5

$\mathcal{O}$ Domain occupied by the ocean basin p. 5

$\partial\mathcal{O}$ Boundary of the ocean basin p. 5

$\bar{z}$ Conjugate of a complex number (z) p. 5

$\tilde{\nabla}^2$ Normalised Laplacian operator p. 5

$\mathbf{ALFs}$ Associated Legendre Functions of the first kind p. 5

$\Phi_k$ k-th eigenfunction of the set of SCHNs p. 5

$\lambda_k$ Eigenvalue associated with $\Phi_k$ p. 5

$\Psi_k$ k-th eigenfunction of the set of SCHDs p. 5

$\nu_k$ Eigenvalue associated with $\Psi_k$ p. 5

$i$ Imaginary number p. 5

$\mathbf{SCHs}$ Spherical cap harmonics p. 5

$\mathbf{SCHNs}$ SCHs with Neumann boundary conditions p. 5

$\mathbf{SCHDs}$ SCHs with Dirichlet boundary conditions p. 5

$\delta_{k,j}$ Kronecker delta function p. 5

$p_k$ Weighting coefficient of $\Phi_k$ p. 5

$p_{-k}$ Weighting coefficient of $\Psi_k$ p. 5

$\mathbf{SPHs}$ Spherical harmonics p. 5

$\beta_{k,j}$ Gyroscopic coefficient p. 6

$\hat{Y}_l^m$ Complex SPH associated with the coordinates of $\mathcal{R}$ p. 7

$l$ Degree (or latitudinal wavenumber) of the SPHs p. 7

$m$ Order (or longitudinal wavenumber) of the SPHs p. 7

$\tilde{U}_{T,l}^{m,\sigma}$ $\hat{Y}_l^m$-component of the normalised tidal gravitational potential p. 7

$\hat{Y}_k$ k-th function of the set of SPHs associated with $\mathcal{R}$ p. 7

$\tilde{U}_{T,k}^{\sigma}$ k-th component of the normalised tidal potential p. 7

$\tilde{U}_T^{\sigma}$ $\sigma$-component of the normalised tidal gravitational potential p. 7

$\tilde{\Xi}^{\sigma}$ $\sigma$-component of the forcing term $\tilde{\Xi}$ p. 7

$\tilde{U}_{T,j}^{\sigma}$ $\hat{Y}_j$-component of $\tilde{U}_T^{\sigma}$ p. 7

$\tilde{\Xi}_j^{\sigma}$ $\hat{Y}_j$-component of $\tilde{\Xi}$ p. 7

$\gamma_{D,l}^{\sigma}, \gamma_{T,l}^{\sigma}$ Solid deformation factors p. 7

$k_j^{\sigma}$ Degree-l solid tidal gravitational Love number p. 7

$h_l^{\sigma}$ Degree-l solid tidal displacement Love number p. 7

$k_{l,l}^{\sigma}$ Degree-l solid loading tidal Love number p. 7

$h_{l,l}^{\sigma}$ Degree-l solid loading displacement Love number p. 7

$\bar{\rho}$ Mean density of the solid regions p. 7
$v_{k,q}$ Coupling coefficients appearing in the forcing terms of the LTEs p. 7

$M$ Number of basis functions or truncation degree of the set {$\Phi_j$} p. 8

$N$ Number of basis functions or truncation degree of the set {$\Psi_j$} p. 8

$K$ Number of basis functions or truncation degree of the set {$\hat{Y}_j$} p. 8

$\Phi$ Vector defining $\Phi$ in terms of the $\Phi_k$ p. 8

$\Psi$ Vector defining $\Psi$ in terms of the $\Psi_k$ p. 8

$\tilde{U}_T$ Vector defining $\tilde{U}_T^\sigma$ in terms of the SPHs p. 8

$T$ Transpose of a matrix p. 8

$E_{M,k}$ Unit vector of size $M$ made of Kronecker delta functions p. 8

$A_\Phi, R_\Phi$ Matrices accounting for the deformation of the solid regions p. 8

$P_{\hat{Y},\Phi}$ Transition matrix from the SPHs $\hat{Y}_k$ to the SCHs $\Phi_k$ p. 8

$\Gamma_D, \Gamma_T$ Diagonal matrices of solid deformation factors p. 8

$\Lambda, N$ Diagonal matrices of eigenvalues p. 8

$B$ Matrix of gyroscopic coefficients p. 8

$B_{\Phi,\Phi}$ First block matrix of $B$ p. 8

$B_{\Phi,\Psi}$ Second block matrix of $B$ p. 8

$B_{\Psi,\Phi}$ Third block matrix of $B$ p. 8

$B_{\Psi,\Psi}$ Fourth block matrix of $B$ p. 8

$A$ Matrix accounting for the gravitational and pressure effects p. 8

$R$ Matrix describing the coupling between the forcing tidal gravitational potential and the oceanic eigenmodes p. 8

$X$ Vector describing the horizontal displacement in terms of the oceanic SCHs p. 8

$I$ Identity matrix p. 9

$H$ Tidal transfer function matrix p. 9

$U_D$ Gravitational potential of the tidally distorted body p. 9

$\tilde{U}_{D,l}^{m,\sigma}$ Complex $(\sigma, l, m)$-component of $U_D$ p. 9

$\zeta_l^{m,\sigma}$ Complex $(\sigma, l, m)$-component of the normalised surface elevation p. 9

$k_{D,l}^{m,\sigma}$ Complex Love numbers p. 9

$T_z$ Tidal torque exerted about the planet's spin axis p. 9

$dV$ Infinitesimal volume parcel p. 9

$\delta\rho$ Local tidal density variations p. 9

$V_*$ Simply connected region that includes the planet but not the perturber p. 9

$\partial V_*$ Boundary of $V_*$ p. 9

$r_*$ Radius of the spherical volume $V_*$ p. 10

$a$ Semi-major axis of the perturber p. 10

$P_T$ Time-averaged power input by the tidal force p. 10

$\rho$ Local density p. 10

$P_{T,oc}$ Time-averaged tidal input power in the ocean p. 10

$P_{diss;oc}$ Time-averaged tidally dissipated power in the ocean p. 10

$P_{diss}$ Total time-averaged tidally dissipated power p. 11

$P_{T;sol}$ Time-averaged tidal input power in the solid regions p. 11

$P_{diss;sol}$ Time-averaged tidally dissipated power in the solid regions p. 11

$\mu$ Effective shear modulus of the solid regions p. 11

$\tau_M$ Maxwell time p. 11

$\tau_A$ Andrade time p. 11

$\alpha$ Parameter determining the duration of the transient response in the primary creep (Andrade model) p. 11

$n_s$ Mean motion of the perturber p. 12

$M_\zeta$ Lunar mass p. 12

$n_\zeta$ Lunar mean motion p. 12

$\chi$ Normalised semidiurnal tidal frequency p. 13

$\Omega_E$ Spin rate of the actual Earth p. 13

$P_{rot}$ Spin period p. 13

$A_2$ Dimensionless constant parametrising the quadrupolar solid tidal Love number p. 15

$\Gamma$ Gamma function p. 15

* Superscript referring to the values of the reference case defined in Table 2 p. 15

$\hat{\theta}_S^*$ Colatitude of the continental centre in the reference case p. 15

$\theta_c^*$ Angular radius of the continent in the reference case p. 15

$H^*$ Oceanic depth in the reference case p. 15

$\sigma_R^*$ Rayleigh drag frequency in the reference case p. 15

$\mu^*$ Effective shear modulus in the reference case p. 15

$\tau_A^*$ Andrade time in the reference case p. 15

$\sigma_k$ Tidal frequency of the k-th maximum p. 17

$\Delta\sigma_k$ k-th frequency interval between two consecutive maxima p. 17

$\Delta\sigma$ Mean frequency interval separating two consecutive relative maxima of the tidal torque p. 18

$\hat{\Delta}\sigma_k$ k-th frequency interval normalised by the mean interval p. 18

$CDF$ Cumulative distribution function p. 18

$\hat{\sigma}$ Normalised standard deviation of the intervals separating consecutive relative maxima of the tidal torque p. 18

## Appendix B: Associated Legendre functions (ALFs)

The ALFs, usually denoted by $P_l^m$, are the solutions of the equation

$$ \frac{1}{\sin \theta} \frac{d}{d\theta} \left( \sin \theta \frac{dh}{d\theta} \right) + \left[ l(l+1) - \frac{m^2}{\sin^2 \theta} \right] h = 0, \quad (B.1) $$

which is found by separating the longitude ($\varphi$) and the colatitude ($\theta$) in the eigenvalues problems associated with the Laplace equation in spherical coordinates (see Eq. (C.1)). The function $h(\theta)$ in Eq. (B.1) is assumed to be a regular function of the colatitude. The parameters $l$ and $m$ are commonly referred to as the degree and the order of the Legendre function, respectively. They may be integral, non-integral, or complex depending on the eigenvalues of the problem, $m^2$ and $l(l+1)$, which are determined both by the boundary conditions and the geometry of the domain of definition of the solutions on the unit sphere.

In the general case, the ALFs are defined, for $z \in \mathbb{C}$ such that $|z| < 1$ and arbitrary complex constants, $l$ and $m$, as (Abramowitz & Stegun 1972, Eq. (8.1.2))

$$ P_l^m(z) \equiv \frac{1}{\Gamma(1-m)} \left( \frac{1+z}{1-z} \right)^{\frac{m}{2}} {}_2F_1 \left( -l, l+1; 1-m; \frac{1-z}{2} \right). \quad (B.2) $$
In the above equation, $\Gamma$ is the Gamma function introduced in Eq. (106), which is defined for $\mathfrak{R}(z) > 0$ as (Abramowitz & Stegun 1972, Eq. (6.1.1))

$$ \Gamma(z) \equiv \int_{0}^{\infty} t^{z-1} e^{-t} dt. \qquad (B.3) $$

The hypergeometric function ${}_2F_1$ is defined as (Abramowitz & Stegun 1972, Eq. (15.1.1))

$$ {}_2F_1(a, b; c; d) = \sum_{j=0}^{\infty} \frac{(a)_j (b)_j}{(c)_j} \frac{x^j}{j!}, \qquad (B.4) $$

where $(a)_j$ designates the Pochhammer symbol (Abramowitz & Stegun 1972, Eq. (6.1.22)),

$$ (a)_0 \equiv 1, \qquad (B.5) $$

$$ (a)_j \equiv a(a+1)\dots(a+j-1) = \frac{\Gamma(a+j)}{\Gamma(a)}. \qquad (B.6) $$

The derivatives of the ALFs with respect to $z$ are given by

$$ \partial_z P_l^m = \frac{1}{z^2 - 1} \left[ z l P_l^m(z) - (m+l) P_{l-1}^m(z) \right]. \qquad (B.7) $$

We note that the functions $\Gamma$ and ${}_2F_1$ in Eq. (B.2) are defined only if their arguments satisfy certain conditions. In practice, the ALFs are computed for any values of these arguments by using the symmetry properties of the regularised hypergeometric function ${}_2\tilde{F}_1$, which is defined as (e.g. Johansson 2016)

$$ {}_2\tilde{F}_1(a, b; c; z) \equiv \frac{{}_2F_1(a, b; c; z)}{\Gamma(c)}. \qquad (B.8) $$

For $m \in \mathbb{Z}$, $l \in \mathbb{R}$ and $z = x \in \mathbb{R}$ such that $-1 \le x \le 1$, which corresponds to the case of the SCHs, the general definition given by Eq. (B.2) simplifies to (e.g. Thébault et al. 2006)

$$ P_l^m(x) = \frac{(-1)^{\frac{m+|m|}{2}}}{2^{|m|} |m|!} \frac{\Gamma(l + |m| + 1)}{\Gamma(l - m + 1)} (1 - x^2)^{|m|/2} \\ \times {}_2\tilde{F}_1 \left( |m| - l, l + |m| + 1; |m| + 1; \frac{1-x}{2} \right). \qquad (B.9) $$

Here, the hypergeometric function is expressed as

$$ {}_2\tilde{F}_1 \left( |m| - l, l + |m| + 1; |m| + 1; \frac{1-x}{2} \right) = \sum_{j=0}^{\infty} a_j \left( \frac{1-x}{2} \right)^j, \qquad (B.10) $$

with the coefficients $a_j$ recursively defined as

$$ a_0 = 1, \qquad (B.11) $$

$$ a_{j+1} = \frac{(j + |m| - l)(j + |m| + l + 1)}{(j + 1)(j + 1 + |m|)} a_j. \qquad (B.12) $$

In the general case, the sum of Eq. (B.9) converges only for $-1 < x \le 1$, which corresponds to the inequality $0 \le \theta < \pi$ in the change of coordinates $x = \cos \theta$. However the sum of Eq. (B.9) is not infinite any more when $l$ is an integer since $a_{j+1} = 0$ as long as $j > l - |m|$. As a consequence it converges also for $\theta = \pi$ in this case.

Since the degrees $l$ tend to integral values ($l = 0, 1, 2, \dots$) as $\theta \to \pi$, the generalised ALFs given by Eqs. (B.2) and (B.9) continuously converge towards the well-known ALFs of the SPHs, which are defined for $-1 \le x \le 1$, and $l$ and $m$ integers such that $l \ge |m|$, as (Abramowitz & Stegun 1972; Arfken & Weber 2005)

$$ P_l^m(x) \equiv (-1)^m (1 - x^2)^{m/2} \frac{d^m}{dx^m} P_l(x), \qquad (B.13) $$

the $P_l$ designating the Legendre polynomials, given by

$$ P_l(x) \equiv \frac{1}{2^l l!} \frac{d^l}{dx^l} \left[ (x^2 - 1)^l \right]. \qquad (B.14) $$

Owing to the large number of terms that is sometimes required to compute the sum of Eq. (B.10) at machine precision, the numerical evaluation of the generalised ALFs given by Eq. (B.9) may be tricky in some cases. Particularly, the series converges very slowly when $x \approx -1$, which occurs if the ocean is almost global. In practice, the evaluation is achieved by making use of the build-in hypergeometric function ${}_2F_1$ of the TRIP algebraic manipulator (Gastineau & Laskar 2011).

## Appendix C: Spherical cap harmonics (SCHs)

In the same way that the SPHs are the harmonics of the entire sphere, the SCHs designate the harmonics of the sphere truncated at a colatitude $\theta_0$ (e.g. Haines 1985; Thébault et al. 2004, 2006; Torta 2019, and references therein). Figure C.1 shows the two geometrical configurations. We note that the natural coordinate system of the SCHs is such that the z-axis is the axis going through the centre of the bounded ocean, which is the symmetry axis of the system, while the coordinate system of the SPHs does not need to be specified due to the spherical symmetry of the global ocean.

The SPHs and the SCHs are found by solving – for specific boundary conditions discussed further – the eigenvalue-eigenfunction problem described by the Helmholtz equation (e.g. Haines 1985),

$$ \nabla^2 f = -\lambda f, \qquad (C.1) $$

where $\nabla^2$ designates the Laplace operator, defined for any function $f(\theta, \varphi)$ as

$$ \nabla^2 f \equiv \frac{1}{\sin \theta} \partial_{\theta} (\sin \theta \partial_{\theta} f) + \frac{1}{\sin^2 \theta} \partial_{\varphi} f, \qquad (C.2) $$

the notation $\partial_x$ or $\partial_{xx}$ referring to the first or second partial derivatives with respect to the $x$ coordinate. In both cases, the variables $\theta$ and $\varphi$ can be separated, which allows the problem to be solved as two decoupled eigenvalue problems.

Since the functions $f(\theta, \varphi)$ in Eq. (C.1) are defined for all $\varphi$, the boundary conditions on $\varphi$ are those of continuity,

$$ f(\theta, \varphi) = f(\theta, \varphi + 2\pi), \qquad (C.3) $$

$$ \partial_{\varphi} f(\theta, \varphi) = \partial_{\varphi} f(\theta, \varphi + 2\pi), \qquad (C.4) $$

which restricts the values of $m$ to be integral ($m \in \mathbb{Z}$). These conditions hold both for the ordinary SPHs and for the SCHs. Analogously, the two sets of basis functions share the same regularity condition on $\theta$ at the pole ($\theta = 0$),

$$ \partial_{\theta} f(0, \varphi) = 0 \qquad m = 0, \qquad (C.5) $$

$$ f(0, \varphi) = 0 \qquad m \neq 0. \qquad (C.6) $$
<figure><img src="image_8.png" /><figcaption>**Fig. C.1.** Geometry of the land-ocean distribution in the model. Left: global ocean distribution where the eigenfunctions are the SPHs. Right: bounded circular ocean distribution where the eigenfunctions are the SCHs. In the bounded case, the frame of reference $\mathcal{R}_{oc}: (O, \mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$ introduced in Sect. 2.1 is such that $\mathbf{e}_z$ corresponds to the symmetry axis going through the centre of the ocean basin, while it is not specified in the global case due to the spherical symmetry.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. C.2.** Real degrees $l_k$ of the SCHs for doublets ($k, |m|$) with $0 \le |m| \le k \le 2$ as functions of the angular radius of the supercontinent. The set of basis functions obtained for the Neumann condition given by Eq. (C.15) (SCHNs) are designated by solid lines, and those obtained for the Dirichlet condition given by Eq. (C.16) (SCHDs) by dashed lines.</figcaption></figure>

However, they differ as regards the second boundary condition. The SPHs are obtained by solving the Helmholtz equation given by Eq. (C.1) on the entire sphere, $S: (\theta, \varphi) \in [0, \pi] \times [0, 2\pi[$, with regularity boundary conditions at $\theta = \pi$ similar to those applied at $\theta = 0$,

$$ \partial_{\theta} f(\pi, \varphi) = 0 \qquad m = 0, \qquad (C.7) $$

$$ f(\pi, \varphi) = 0 \qquad m \neq 0. \qquad (C.8) $$

They are associated with integral degrees $l \in \mathbb{N}$ such that $l \ge |m|$.

Following Varshalovich et al. (1988), we use, for the complex SPHs, the convention given by (Varshalovich et al. 1988, Sect. 5.2, Eq. (1))

$$ Y_l^m(\theta, \varphi) = \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} P_l^m(\cos\theta) e^{im\varphi}, \qquad (C.9) $$

or, equivalently,

$$ Y_l^m(\theta, \varphi) = (-1)^{\frac{|m|-m}{2}} \sqrt{\frac{(2l+1)(l-|m|)!}{4\pi(l+|m|)!}} P_l^{|m|}(\cos\theta) e^{im\varphi}, \qquad (C.10) $$

where the $P_l^m$ are the standard ALFs defined in Eq. (B.13). The eigenvalues associated with the SPHs only depend on the degree $l$, and are expressed as

$$ \lambda_l = l(l+1). \qquad (C.11) $$

In the used convention (Eq. (C.9)), the SPHs form a set of orthonormal basis functions through the scalar product defined, for any complex functions $f$ and $g$, as

$$ \begin{aligned} \langle f, g \rangle &= \int_S \bar{f} g dS, \qquad (C.12) \\ &= \int_{\varphi=0}^{2\pi} \int_{\theta=0}^{\pi} \overline{f(\theta, \varphi)} g(\theta, \varphi) \sin\theta d\theta d\varphi. \end{aligned} $$

the notation $dS$ referring to an infinitesimal surface element of the unit sphere, and $\bar{f}$ to the conjugate of $f$. This scalar product allows the $l^2$-norm of the space of functions of $S$ to be defined,

$$ \|f\|_2 = \sqrt{\langle f, f \rangle}. \qquad (C.13) $$

Hence, for any pair of SPHs $Y_l^m$ and $Y_p^q$,

$$ \langle Y_l^m, Y_p^q \rangle = \delta_{l,p} \delta_{m,q}. \qquad (C.14) $$

The SCHs are obtained by solving the Helmholtz equation given by Eq. (C.1) on the domain filled by the spherical cap, $C: (\theta, \varphi) \in [0, \theta_0] \times [0, 2\pi[$, together with the boundary conditions specified at $\partial C$, namely $\theta = \theta_0$. In the present model, the applied boundary conditions are either Neumann (fixed value for $\partial_{\theta}f$ at $\theta = \theta_0$) or Dirichlet (fixed value for the function $f$ itself) conditions, which are respectively formulated as

$$ \partial_{\theta} f(\theta_0, \varphi) = 0 \qquad (\text{Neumann condition}), \qquad (C.15) $$
<table><caption>**Table C.1.** Values of the real degrees $l_n$ of the SCHs with Neumann (SCHNs) or Dirichlet (SCHDs) boundary conditions for various angular radii of the supercontinent ($\theta_c$).</caption><thead><tr><th>θ<sub>c</sub> (deg)</th><th>(0, 0)</th><th>(1, 0)</th><th>(1, 1)</th><th>(2, 0)</th><th>(2, 1)</th><th>(2, 2)</th></tr></thead><tbody><tr><td colspan="7">Neumann boundary conditions (SCHNs)</td></tr><tr><td>0</td><td>0.000000</td><td>1.000000</td><td>1.000000</td><td>2.000000</td><td>2.000000</td><td>2.000000</td></tr><tr><td>10</td><td>0.000000</td><td>1.014533</td><td>0.984857</td><td>2.040535</td><td>1.957224</td><td>1.999318</td></tr><tr><td>20</td><td>0.000000</td><td>1.054431</td><td>0.945331</td><td>2.139648</td><td>1.889718</td><td>1.989826</td></tr><tr><td>30</td><td>0.000000</td><td>1.115646</td><td>0.900510</td><td>2.277970</td><td>1.885701</td><td>1.956707</td></tr><tr><td>40</td><td>0.000000</td><td>1.196823</td><td>0.867186</td><td>2.450068</td><td>1.948680</td><td>1.901867</td></tr><tr><td>50</td><td>0.000000</td><td>1.298769</td><td>0.851964</td><td>2.657126</td><td>2.062763</td><td>1.851211</td></tr><tr><td>60</td><td>0.000000</td><td>1.424123</td><td>0.856313</td><td>2.904340</td><td>2.220869</td><td>1.827443</td></tr><tr><td>70</td><td>0.000000</td><td>1.577470</td><td>0.880927</td><td>3.200567</td><td>2.423975</td><td>1.840734</td></tr><tr><td>80</td><td>0.000000</td><td>1.765876</td><td>0.927569</td><td>3.559164</td><td>2.679212</td><td>1.896224</td></tr><tr><td>90</td><td>0.000000</td><td>2.000000</td><td>1.000000</td><td>4.000000</td><td>3.000000</td><td>2.000000</td></tr><tr><td>100</td><td>0.000000</td><td>2.296162</td><td>1.105011</td><td>4.553241</td><td>3.408276</td><td>2.162574</td></tr><tr><td>110</td><td>0.000000</td><td>2.680266</td><td>1.254256</td><td>5.266546</td><td>3.939640</td><td>2.402474</td></tr><tr><td>120</td><td>0.000000</td><td>3.195691</td><td>1.467987</td><td>6.219529</td><td>4.654189</td><td>2.752588</td></tr><tr><td>130</td><td>0.000000</td><td>3.920704</td><td>1.783375</td><td>7.555648</td><td>5.660659</td><td>3.273969</td></tr><tr><td>140</td><td>0.000000</td><td>5.012005</td><td>2.275396</td><td>9.561941</td><td>7.176982</td><td>4.090469</td></tr><tr><td>150</td><td>0.000000</td><td>6.835398</td><td>3.119597</td><td>12.908284</td><td>9.712069</td><td>5.492825</td></tr><tr><td>160</td><td>0.000000</td><td>10.488504</td><td>4.843239</td><td>19.604445</td><td>14.793105</td><td>8.355251</td></tr><tr><td>170</td><td>0.000000</td><td>21.459763</td><td>10.083479</td><td>39.699467</td><td>30.056683</td><td>17.051968</td></tr><tr><td colspan="7">Dirichlet boundary conditions (SCHDs)</td></tr><tr><td>0</td><td>0.000000</td><td>1.000000</td><td>1.000000</td><td>2.000000</td><td>2.000000</td><td>2.000000</td></tr><tr><td>10</td><td>0.201220</td><td>1.286874</td><td>1.014533</td><td>2.359807</td><td>2.040535</td><td>2.000661</td></tr><tr><td>20</td><td>0.274502</td><td>1.424749</td><td>1.054431</td><td>2.561533</td><td>2.139648</td><td>2.009153</td></tr><tr><td>30</td><td>0.346184</td><td>1.568297</td><td>1.115646</td><td>2.777734</td><td>2.277970</td><td>2.037781</td></tr><tr><td>40</td><td>0.422281</td><td>1.726592</td><td>1.196823</td><td>3.019705</td><td>2.450068</td><td>2.095353</td></tr><tr><td>50</td><td>0.506293</td><td>1.905895</td><td>1.298769</td><td>3.296147</td><td>2.657126</td><td>2.186686</td></tr><tr><td>60</td><td>0.601509</td><td>2.112863</td><td>1.424123</td><td>3.616943</td><td>2.904340</td><td>2.315699</td></tr><tr><td>70</td><td>0.711801</td><td>2.355851</td><td>1.577470</td><td>3.994876</td><td>3.200567</td><td>2.487779</td></tr><tr><td>80</td><td>0.842253</td><td>2.646183</td><td>1.765876</td><td>4.447503</td><td>3.559164</td><td>2.711416</td></tr><tr><td>90</td><td>1.000000</td><td>3.000000</td><td>2.000000</td><td>5.000000</td><td>4.000000</td><td>3.000000</td></tr><tr><td>100</td><td>1.195606</td><td>3.441373</td><td>2.296162</td><td>5.690014</td><td>4.553241</td><td>3.374653</td></tr><tr><td>110</td><td>1.445566</td><td>4.008030</td><td>2.680266</td><td>6.576628</td><td>5.266546</td><td>3.869485</td></tr><tr><td>120</td><td>1.777288</td><td>4.762779</td><td>3.195691</td><td>7.758259</td><td>6.219529</td><td>4.542151</td></tr><tr><td>130</td><td>2.240037</td><td>5.818619</td><td>3.920704</td><td>9.412016</td><td>7.555648</td><td>5.497155</td></tr><tr><td>140</td><td>2.932276</td><td>7.401491</td><td>5.012005</td><td>11.892079</td><td>9.561941</td><td>6.944195</td></tr><tr><td>150</td><td>4.083687</td><td>10.038551</td><td>6.835398</td><td>16.024836</td><td>12.908284</td><td>9.373283</td></tr><tr><td>160</td><td>6.383235</td><td>15.311204</td><td>10.488504</td><td>24.289410</td><td>19.604445</td><td>14.255320</td></tr><tr><td>170</td><td>13.275607</td><td>31.126398</td><td>21.459763</td><td>49.081367</td><td>39.699467</td><td>28.946234</td></tr></tbody></table>

Notes. The values are given for the doublets $(n, |m|)$ such that $n = 0, 1, 2$ and $|m| \le n$.

$$f(\theta_0, \varphi) = 0 \qquad (\text{Dirichlet condition}). \qquad (\text{C.16})$$

Each condition leads to a set of orthogonal basis functions. These sets are referred to as the 'SCHNs' for Eq. (C.15) and the 'SCHDs' for Eq. (C.16).

As the Helmholtz equation is the same for both the SCHs and the SPHs, solving the eigenfunction-eigenvalue problem defined by Eqs. (C.1) and (C.15) (or Eq. (C.16)) amounts to finding the set of ALFs (see Eq. (B.9)) that satisfy the specified boundary conditions for a given $m \in \mathbb{Z}$. As a consequence, the eigenfunction-eigenvalue problem becomes a root-finding problem where one computes the series of real degrees $l$ such that, alternately,

$$\frac{dP_l^m}{d\theta} (\cos \theta_0) = 0, \qquad (\text{Neumann condition}), \qquad (\text{C.17})$$

$$P_l^m (\cos \theta_0) = 0 \qquad (\text{Dirichlet condition}). \qquad (\text{C.18})$$

This can be achieved numerically using a standard Newton-Raphson algorithm (e.g. Press et al. 2007, Sect. 9.4) since the ALFs defined by Eq. (B.9) and their derivatives are regular functions of $l$.

For each condition (Eq. (C.15) or Eq. (C.16)), we end up with the degrees of the distorted ALFs, $l_n(m)$, which are subscripted with indices $n \in \mathbb{N}$ such that $n \ge |m|$. These indices<fnref n="3" /> correspond to the integral degrees of the SPHs towards which the SCHs converge for the $l^2$-norm defined by Eq. (C.13) in the global ocean

<footnote n="3">³ Other conventions exist for the indices $n$. For instance, Haines (1985) introduces a convention where the parity of $k - |m|$ allows the two sets of basis functions (with Neumann or Dirichlet conditions) to be distinguished. In this convention, the roots of Eq. (C.17) (Neumann condition) are designated by even $n - |m|$, while those of Eq. (C.18) (Dirichlet condition) are designated by odd $n - |m|$.</footnote>
<table><caption>**Table C.2.** Eigenvalues of the SCHs with Neumann (SCHNs) or Dirichlet (SCHDs) boundary conditions for various angular radii of the supercontinent.</caption><thead><tr><th>θ<sub>c</sub> (deg)</th><th>(0, 0)</th><th>(1, 0)</th><th>(1, 1)</th><th>(2, 0)</th><th>(2, 1)</th><th>(2, 2)</th></tr></thead><tbody><tr><td colspan="7">Neumann boundary conditions (SCHNs)</td></tr><tr><td>0</td><td>0.000000</td><td>2.000000</td><td>2.000000</td><td>6.000000</td><td>6.000000</td><td>6.000000</td></tr><tr><td>10</td><td>0.000000</td><td>2.043810</td><td>1.954799</td><td>6.204319</td><td>5.787950</td><td>5.996589</td></tr><tr><td>20</td><td>0.000000</td><td>2.166257</td><td>1.838982</td><td>6.717744</td><td>5.460750</td><td>5.949235</td></tr><tr><td>30</td><td>0.000000</td><td>2.360311</td><td>1.711429</td><td>7.467118</td><td>5.441569</td><td>5.785407</td></tr><tr><td>40</td><td>0.000000</td><td>2.629208</td><td>1.619197</td><td>8.452903</td><td>5.746033</td><td>5.518966</td></tr><tr><td>50</td><td>0.000000</td><td>2.985569</td><td>1.577808</td><td>9.717446</td><td>6.317753</td><td>5.278194</td></tr><tr><td>60</td><td>0.000000</td><td>3.452250</td><td>1.589586</td><td>11.339534</td><td>7.153126</td><td>5.166990</td></tr><tr><td>70</td><td>0.000000</td><td>4.065883</td><td>1.656959</td><td>13.444194</td><td>8.299630</td><td>5.229035</td></tr><tr><td>80</td><td>0.000000</td><td>4.884195</td><td>1.787952</td><td>16.226814</td><td>9.857388</td><td>5.491887</td></tr><tr><td>90</td><td>0.000000</td><td>6.000000</td><td>2.000000</td><td>20.000000</td><td>12.000000</td><td>6.000000</td></tr><tr><td>100</td><td>0.000000</td><td>7.568523</td><td>2.326061</td><td>25.285246</td><td>15.024623</td><td>6.839299</td></tr><tr><td>110</td><td>0.000000</td><td>9.864090</td><td>2.827415</td><td>33.003048</td><td>19.460404</td><td>8.174357</td></tr><tr><td>120</td><td>0.000000</td><td>13.408133</td><td>3.622974</td><td>44.902072</td><td>26.315661</td><td>10.329330</td></tr><tr><td>130</td><td>0.000000</td><td>19.292626</td><td>4.963802</td><td>64.643464</td><td>37.703724</td><td>13.992842</td></tr><tr><td>140</td><td>0.000000</td><td>30.132199</td><td>7.452826</td><td>100.992648</td><td>58.686045</td><td>20.822406</td></tr><tr><td>150</td><td>0.000000</td><td>53.558065</td><td>12.851483</td><td>179.532083</td><td>104.036347</td><td>35.663952</td></tr><tr><td>160</td><td>0.000000</td><td>120.497214</td><td>28.300199</td><td>403.938705</td><td>233.629059</td><td>78.165464</td></tr><tr><td>170</td><td>0.000000</td><td>481.981173</td><td>111.760030</td><td>1615.747119</td><td>933.460872</td><td>307.821577</td></tr><tr><td colspan="7">Dirichlet boundary conditions (SCHDs)</td></tr><tr><td>0</td><td>0.000000</td><td>2.000000</td><td>2.000000</td><td>6.000000</td><td>6.000000</td><td>6.000000</td></tr><tr><td>10</td><td>0.241710</td><td>2.942918</td><td>2.043810</td><td>7.928496</td><td>6.204319</td><td>6.003307</td></tr><tr><td>20</td><td>0.349854</td><td>3.454660</td><td>2.166257</td><td>9.122983</td><td>6.717744</td><td>6.045851</td></tr><tr><td>30</td><td>0.466027</td><td>4.027851</td><td>2.360311</td><td>10.493540</td><td>7.467118</td><td>6.190333</td></tr><tr><td>40</td><td>0.600603</td><td>4.707711</td><td>2.629208</td><td>12.138321</td><td>8.452903</td><td>6.485855</td></tr><tr><td>50</td><td>0.762626</td><td>5.538331</td><td>2.985569</td><td>14.160734</td><td>9.717446</td><td>6.968283</td></tr><tr><td>60</td><td>0.963323</td><td>6.577054</td><td>3.452250</td><td>16.699223</td><td>11.339534</td><td>7.678163</td></tr><tr><td>70</td><td>1.218462</td><td>7.905884</td><td>4.065883</td><td>19.953910</td><td>13.444194</td><td>8.676820</td></tr><tr><td>80</td><td>1.551642</td><td>9.648467</td><td>4.884195</td><td>24.227788</td><td>16.226814</td><td>10.063195</td></tr><tr><td>90</td><td>2.000000</td><td>12.000000</td><td>6.000000</td><td>30.000000</td><td>20.000000</td><td>12.000000</td></tr><tr><td>100</td><td>2.625081</td><td>15.284420</td><td>7.568523</td><td>38.066269</td><td>25.285246</td><td>14.762935</td></tr><tr><td>110</td><td>3.535228</td><td>20.072333</td><td>9.864090</td><td>49.828657</td><td>33.003048</td><td>18.842401</td></tr><tr><td>120</td><td>4.936042</td><td>27.446847</td><td>13.408133</td><td>67.948839</td><td>44.902072</td><td>25.173291</td></tr><tr><td>130</td><td>7.257803</td><td>39.674941</td><td>19.292626</td><td>97.998058</td><td>64.643464</td><td>35.715869</td></tr><tr><td>140</td><td>11.530519</td><td>62.183561</td><td>30.132199</td><td>153.313622</td><td>100.992648</td><td>55.166033</td></tr><tr><td>150</td><td>20.760187</td><td>110.811047</td><td>53.558065</td><td>272.820192</td><td>179.532083</td><td>97.231717</td></tr><tr><td>160</td><td>47.128920</td><td>249.744161</td><td>120.497214</td><td>614.264854</td><td>403.938705</td><td>217.469478</td></tr><tr><td>170</td><td>189.517352</td><td>999.979043</td><td>481.981173</td><td>2458.061916</td><td>1615.747119</td><td>866.830706</td></tr></tbody></table>

**Notes.** The values of $\lambda_l$ are were computed from Eq. (C.20) for the doublets $(n, |m|)$ such that $n = 0, 1, 2$ and $|m| \le n$ using the values of real degrees given by Table C.1.

limit: $l_n \to n$ as $\theta_0 \to 180^\circ$ independently of the applied boundary condition. Similarly as integral degrees, the real degrees $l_n$ are ranked in ascending order. For any $n$ and $j$ such that $n < j$, the inequality $l_n < l_j$ is verified. The values of the real degrees are given by Table C.1 for $0 \le |m| \le 2$ and various values of angular radius of the continent, $\theta_c = 180 - \theta_0$ (in degrees). They are plotted on Fig. C.2.

We note that $l_n \ge |m|$ when Dirichlet conditions are applied, similarly as in the case of SPHs (e.g. Thébault et al. 2006). However, when the Neumann conditions are applied, this is only valid for $\theta_0 \le 90^\circ$ since there exists degrees $l_{|m|} < |m|$ as long as $\theta_0 > 90^\circ$. Asymptotic scaling laws may be used to estimate the values of the degrees $l_n \gg |m|$ as functions of the indices $n$ and the angle $\theta_0$ (see e.g. Thébault et al. 2006, Eq. (40)). Particularly, these laws show that $l_n \propto \theta_0^{-1}$ for $\theta_0 \ll 1$ consistently with what

may be observed in Fig. C.2. This results from the fact that the degree $l$ of an harmonic scales as the inverse of the corresponding wavelength, which is determined by the size of the spherical cap.

The resulting SCHs are defined on the unit sphere, for $|m| = 0, \dots, +\infty$ and $l = l_{|m|}, l_{|m|+1}, \dots,$ as

$$ \Theta_l^m(\theta, \varphi) \equiv \begin{cases} \alpha_{l,m} P_l^{|m|} (\cos \theta) e^{im\varphi} & \text{if } (\theta, \varphi) \in C, \\ 0 & \text{if } (\theta, \varphi) \in S \setminus C \end{cases} \quad (C.19) $$

where the normalisation coefficients $\alpha_{l,m} \in \mathbb{R}$ are to be specified. The corresponding eigenvalues $\lambda_l$ are expressed as

$$ \lambda_l(m) = l(m) [l(m) + 1], \quad (C.20) $$

similarly as in the case of SPHs, though they now also depend on the order $m$ through $l$. The eigenvalues associated with the
SCHNs and SCHDs are given in Table C.2 for $0 \le |m| \le n \le 2$ and various values of the continental angular radius ($\theta_c$). They are computed straightforwardly from the values of the real degrees given in Table C.1 using the expression given by Eq. (C.20).

The SCHs obtained for a given boundary condition are orthogonal basis functions of the space of regular functions on C through the scalar product introduced in Eq. (C.12). Since the harmonics are set to zero for $\theta > \theta_0$, this scalar product reads

$$ \langle f, g \rangle = \int_{\varphi=0}^{2\pi} \int_{\theta=0}^{\theta_0} \overline{f(\theta, \varphi)} g(\theta, \varphi) \sin \theta \, d\theta \, d\varphi. \quad (C.21) $$

For any pair of SCHs belonging to the same set, $\Theta_l^m$ and $\Theta_p^q$, the property given by Eq. (C.14) for the SPHs,

$$ \langle \Theta_l^m, \Theta_p^q \rangle = \delta_{l,p} \delta_{m,q}. \quad (C.22) $$

is enforced by defining the normalisation coefficients $\alpha_{l,m} \in \mathbb{R}$ as

$$ \alpha_{l,m} = (-1)^{\frac{|m|-|n|}{2}} \left[ 2\pi \int_{0}^{\theta_0} \left[ P_l^{|m|} (\cos \theta) \right]^2 \sin \theta \, d\theta \right]^{-1/2} . \quad (C.23) $$

However, the functions in one set are not orthogonal to those in the other since (e.g. Haines 1985)

$$ \langle \Theta_l^m, \Theta_p^q \rangle = \delta_{m,q} \frac{2\pi \sin \theta_0 P_l^m (\cos \theta_0) dP_p^q}{(p-l)(p+l+1) \sin \theta_0} \quad (C.24) $$

for $\Theta_l^m$ and $\Theta_p^q$ taken in the sets of SCHNs and SCHDs, respectively. In the general case, $P_l^m(\cos \theta_0) \neq 0$ and $\frac{dP_p^q}{d\theta}(\cos \theta_0) \neq 0$, which implies $\langle \Theta_l^m, \Theta_p^q \rangle \neq 0$. This highlights the fact that the SCHDs and the SCHNs are two distinct complete sets of basis functions for the space of functions defined on C.

In the case of real eigenfunctions defined as

$$ \Theta_l^m(\theta, \varphi) = \begin{cases} \alpha_{l,m} P_l^{|m|} (\cos \theta) \cos (m\varphi + \varphi_m) & \text{if } (\theta, \varphi) \in C, \\ 0 & \text{if } (\theta, \varphi) \in S \setminus C, \end{cases} \quad (C.25) $$

with $\varphi_m = 0$ if $m \ge 0$ and $\varphi_m = \pi/2$ if $m < 0$, the orthogonality property given by Eq. (C.22) still holds and the normalisation coefficients defined in Eq. (C.23) become

$$ \alpha_{l,m} = (-1)^{\frac{|m|-|n|}{2}} \left[ (1 + \delta_{m,0}) \pi \int_{0}^{\theta_0} \left[ P_l^{|m|} (\cos \theta) \right]^2 \sin \theta \, d\theta \right]^{-1/2} . \quad (C.26) $$

We note that the factors $(-1)^k$ with $k = (|m| - m)/2$ or $k = (|m| + m)/2$ in the expressions of the normalisation coefficients are introduced to make the SCHs consistent with the SPHs in the used convention (Eq. (C.9)).

As $\theta_0$ increases, the ALFs smoothly evolve from the Bessel functions ($\theta_0 \to 0$) to the standard ALFs of integral parameters ($\theta_0 \to \pi$). Similarly, the SCHs evolve from the cylindrical harmonics, where there is no curvature effect (small-sea configuration), to the SPHs given by Eq. (C.9) (global-ocean configuration). Nevertheless, we emphasise that both the eigenvalues and the eigenfunctions depend on the specified boundary condition through the degrees $l$, different conditions leading to different

<figure><img src="image_10.png" /><figcaption>Fig. D.1. Euler angles corresponding to the standard Euler rotation matrix defined by Eq. (D.1).</figcaption></figure>

sets of eigenfunctions. Besides, the limit $\theta_0 \to \pi$ is not strictly equivalent to the global-ocean configuration since boundary conditions are still applied at $\partial C$, instead of the periodic conditions of the SPHs. In other words, the presence of a small island still alters the eigenmodes of the oceanic tide locally although this effect is negligible at a global scale. The impact of the island on the horizontal structure of tidal flows is quantified by the differences of degrees, $|l_n - n|$, and overlap coefficients, $|1 - \langle \Theta_l^m, Y_n^m \rangle|$, both of them tending to zero as the size of the island decreases.

## Appendix D: Euler rotation matrix

We consider the frame of reference $\mathcal{R}: (O, e_X, e_Y, e_Z)$ having the planet's centre of gravity as origin and associated with the Cartesian coordinates $(X, Y, Z)$. As shown by Fig. D.1, any rotation $\mathcal{R}: (O, e_X, e_Y, e_Z) \to \mathcal{R}_{\text{rot}}: (O, e_X, e_Y, e_Z)$ of this frame of reference may be performed by three rotations about the coordinate axes (e.g. Varshalovich et al. 1988, Sect. 1.4): (i) a rotation about the Z-axis through an angle $\alpha$ such that $0 \le \alpha < 2\pi$, (ii) a rotation about the new Y'-axis through an angle $\beta$ such that $0 \le \beta \le \pi$, and (iii) a rotation about the new axis Z through an angle $\gamma$ such that $0 \le \gamma \le 2\pi$. The frame of reference $\mathcal{R}_{\text{rot}}: (O, e_X, e_Y, e_Z)$ is associated with the Cartesian coordinates $(x, y, y)$. The three rotation angles, $\alpha, \beta$, and $\gamma$, are the so-called Euler angles.

Introducing the notations $C_\phi = \cos \phi$ and $S_\phi = \sin \phi$ for $\phi = \alpha, \beta$, or $\gamma$, the standard Euler rotation matrix is defined in terms of the Euler angles as (e.g. Varshalovich et al. 1988, Eq. (54) p. 30)

$$ \mathbf{R}_{\alpha,\beta,\gamma} = \begin{bmatrix} C_{\alpha}C_{\beta}C_{\gamma} - S_{\alpha}S_{\gamma} & -C_{\alpha}C_{\beta}S_{\gamma} - S_{\alpha}C_{\gamma} & C_{\alpha}S_{\beta} \\ S_{\alpha}C_{\beta}C_{\gamma} + C_{\alpha}S_{\gamma} & -S_{\alpha}C_{\beta}S_{\gamma} + C_{\alpha}C_{\gamma} & S_{\alpha}S_{\beta} \\ -S_{\beta}C_{\gamma} & S_{\beta}S_{\gamma} & C_{\beta} \end{bmatrix}, \quad (\text{D.1}) $$
which is formed by the product of three rotation matrices,

$$ \mathbf{R}_{\alpha,\beta,\gamma} = \mathbf{R}_{3,\alpha} \mathbf{R}_{2,\beta} \mathbf{R}_{3,\gamma} \qquad (\text{D.2}) $$

$$ = \begin{bmatrix} C_{\alpha} & -S_{\alpha} & 0 \\ S_{\alpha} & C_{\alpha} & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} C_{\beta} & 0 & S_{\beta} \\ 0 & 1 & 0 \\ -S_{\beta} & 0 & C_{\beta} \end{bmatrix} \begin{bmatrix} C_{\gamma} & -S_{\gamma} & 0 \\ S_{\gamma} & C_{\gamma} & 0 \\ 0 & 0 & 1 \end{bmatrix}. $$

The inverse rotation matrix $\mathbf{R}^{-1}$ is given by

$$ \mathbf{R}_{\alpha,\beta,\gamma}^{-1} = \mathbf{R}_{\alpha,\beta,\gamma}^{\mathrm{T}} = \mathbf{R}_{-\gamma,-\beta,-\alpha}. \qquad (\text{D.3}) $$

The Euler rotation matrix and its inverse correspond to the change of basis matrices between the initial and rotated vectorial bases, $(\mathbf{e}_X, \mathbf{e}_Y, \mathbf{e}_Z) \rightarrow (\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$, which is expressed as

$$ \begin{bmatrix} \mathbf{e}_x \\ \mathbf{e}_y \\ \mathbf{e}_z \end{bmatrix} = \mathbf{R}_{\alpha,\beta,\gamma}^{-1} \begin{bmatrix} \mathbf{e}_X \\ \mathbf{e}_Y \\ \mathbf{e}_Z \end{bmatrix}, \qquad \begin{bmatrix} \mathbf{e}_X \\ \mathbf{e}_Y \\ \mathbf{e}_Z \end{bmatrix} = \mathbf{R}_{\alpha,\beta,\gamma} \begin{bmatrix} \mathbf{e}_x \\ \mathbf{e}_y \\ \mathbf{e}_z \end{bmatrix}. \qquad (\text{D.4}) $$

Besides, the rotation $U \rightarrow U_{\text{rot}}$ of any vector $U$ through the Euler angles $\alpha, \beta$, and $\gamma$ in the Cartesian coordinate system $S \{X, Y, Z\}$ is expressed as $U_{\text{rot}} = \mathbf{R}_{\alpha,\beta,\gamma}U$.

In spherical coordinates, the transformation relations are deduced from Eq. (D.4) by making use of trigonometric identities. We denote by $(\hat{r}, \hat{\theta}, \hat{\phi})$ the spherical coordinates associated with the initial frame of reference, $\mathcal{R}$, where $\hat{r}$ designates the radial coordinate, $\hat{\theta}$ the colatitude, and $\hat{\phi}$ the longitude. The analogous coordinates in the rotated frame of reference, $\mathcal{R}_{\text{rot}}$, are denoted by $(r, \theta, \varphi)$, with $r = \hat{r}$. In the forward transformation, $(\hat{r}, \hat{\theta}, \hat{\phi}) \rightarrow (r, \theta, \varphi)$, the relations between angles $(\hat{\theta}, \hat{\phi})$ and $(\theta, \varphi)$ are given by (Varshalovich et al. 1988, Sect. 1.4.1, Eqs (2-3))

$$ \cos \theta = \cos \hat{\theta} \cos \beta + \sin \hat{\theta} \sin \beta \cos (\hat{\phi} - \alpha), \qquad (\text{D.5}) $$

$$ \cot(\varphi + \gamma) = \cot(\hat{\varphi} - \alpha) \cos \beta - \frac{\cot \hat{\theta} \sin \beta}{\sin(\hat{\varphi} - \alpha)}, \qquad (\text{D.6}) $$

while the inverse relations are

$$ \cos \hat{\theta} = \cos \theta \cos \beta - \sin \theta \sin \beta \cos (\varphi + \gamma), \qquad (\text{D.7}) $$

$$ \cot(\hat{\varphi} - \alpha) = \cot(\varphi + \gamma) \cos \beta + \frac{\cot \theta \sin \beta}{\sin(\varphi + \gamma)}. \qquad (\text{D.8}) $$

As an example, we consider the rotation generating the coordinate system of the ocean basin from the geocentric coordinate system rotating with the planet. From the definition of the ocean centre in the usual geocentric frame of reference, we obtain

$$ \mathbf{e}_z = \cos \hat{\phi}_{\text{oc}} \sin \hat{\theta}_{\text{oc}} \mathbf{e}_X + \sin \hat{\phi}_{\text{oc}} \sin \hat{\theta}_{\text{oc}} \mathbf{e}_Y + \cos \hat{\theta}_{\text{oc}} \mathbf{e}_Z. \qquad (\text{D.9}) $$

In addition, by considering the same vector as the image of $\mathbf{e}_Z$ through a rotation of Euler angles $\alpha, \beta$, and $\gamma$, we can write $\mathbf{e}_z = \mathbf{R}_{\alpha,\beta,\gamma} (0, 0, 1)^{\mathrm{T}}$, where we made use of the standard Euler rotation matrix defined in Eq. (D.1). It follows

$$ \mathbf{e}_z = \cos \alpha \sin \beta \mathbf{e}_X + \sin \alpha \sin \beta \mathbf{e}_Y + \cos \beta \mathbf{e}_Z. \qquad (\text{D.10}) $$

By comparing Eqs. (D.9) and (D.10), we identify $\beta = \hat{\theta}_{\text{oc}}$ and $\alpha = \hat{\phi}_{\text{oc}}$. Considering the inverse transformation, $\mathbf{e}_Z = \mathbf{R}_{\alpha,\beta,\gamma}^{\mathrm{T}} (0, 0, 1)^{\mathrm{T}}$, we get

$$ \mathbf{e}_Z = -\sin \beta \cos \gamma \mathbf{e}_X + \sin \beta \sin \gamma \mathbf{e}_Y + \cos \beta \mathbf{e}_Z, \qquad (\text{D.11}) $$

which shows that $\mathbf{e}_Z$ does not depend on the rotation through the angle $\alpha$ in the oceanic frame of reference.

We note that the angle $\gamma$ describes a rotation of coordinates around the axis defined by $\mathbf{e}_z$. Since the ocean is circular, the geometry remains unchanged through this rotation. The angle $\gamma$ can therefore be set to $\gamma = 0$, which allows the position of the North pole in the ocean's coordinate system to be simplified to $\mathbf{e}_Z = -\sin \beta \mathbf{e}_X + \cos \beta \mathbf{e}_Z$. The coordinates of the North pole in the coordinate system defined from the centre of the ocean basin are thus given by $(\theta_{\text{NP}}, \varphi_{\text{NP}}) = (\beta, \pi)$. At the end, the ocean's coordinate system is simply obtained from the usual geocentric coordinate system through a rotation of Euler angles $(\alpha, \beta, \gamma) = (\hat{\phi}_{\text{oc}}, \hat{\theta}_{\text{oc}}, 0)$. We shall emphasise that the gyroscopic coefficients of the tidal theory only depend on the colatitude of the supercontinent's centre since $\beta$ is the only angle necessary to define the position of the spin axis in the coordinate system of the ocean basin. Moreover, the rotation of angle $\alpha$ can be ignored owing to the periodicity of the tidal response across the longitudinal coordinate, $\hat{\phi}$, which implies that the longitude of the supercontinent's centre can be set to $\hat{\phi}_S = 0$ without any loss of generality.

## Appendix E: Rotation of spherical harmonics

We consider a rotation $\mathcal{R}: (O, \mathbf{e}_X, \mathbf{e}_Y, \mathbf{e}_Z) \rightarrow \mathcal{R}_{\text{rot}}: (O, \mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$ described by the Euler angles introduced in Appendix D, which are denoted by $\alpha, \beta$, and $\gamma$ (see Fig. D.1). The initial frame of reference, $\mathcal{R}$, is associated with the spherical coordinate system $(\hat{r}, \hat{\theta}, \hat{\phi})$, while the rotated frame of reference, $\mathcal{R}_{\text{rot}}$, is associated with the coordinate system $(r, \theta, \varphi)$. In the convention defined by Eq. (C.9), the complex SPHs of the rotated coordinate system are expressed as functions of those defined in the initial coordinate system (Varshalovich et al. 1988, Sect. 4.1),

$$ Y_l^m(\theta, \varphi) = \sum_{q=-l}^{l} D_{q,m}^l(\alpha, \beta, \gamma) Y_l^q(\hat{\theta}, \hat{\phi}), \qquad (\text{E.1}) $$

and, conversely, the SPHs of the initial coordinate system are expressed as

$$ Y_l^m(\hat{\theta}, \hat{\phi}) = \sum_{q=-l}^{l} \overline{D_{m,q}^l(\alpha, \beta, \gamma)} Y_l^q(\theta, \varphi). \qquad (\text{E.2}) $$

The notation $D_{q,m}^l$ refers to the Wigner D-functions (Varshalovich et al. 1988, Sect. 4.3, Eq. (1)),

$$ D_{q,m}^l(\alpha, \beta, \gamma) \equiv e^{-iq\alpha} d_{q,m}^l(\beta) e^{-im\gamma}, \qquad (\text{E.3}) $$

where $d_{q,m}^l$ is a real function expressed as (Varshalovich et al. 1988, Sect. 4.3.1, Eq. (2))

$$ d_{q,m}^l(\beta) = (-1)^{l-m} \sqrt{(l+q)!(l-q)!(l+m)!(l-m)!} \times \sum_k (-1)^k \frac{(\cos \frac{\beta}{2})^{q+m+2k} (\sin \frac{\beta}{2})^{2l-q-m-2k}}{k!(l-q-k)!(l-m-k)!(q+m+k)!} \qquad (\text{E.4}) $$

In the above equation, $k$ runs over all integer values for which the factorial arguments are non-negative, namely $\max(0, -q-m) \le k \le \min(l-q, l-m)$.
<figure><img src="image_11.png" /><figcaption>Fig. E.1. Rotation matrix $D$ of Eq. (E.5) for Euler angles $(\alpha, \beta, \gamma) = (120^{\circ}, 60^{\circ}, 270^{\circ})$. Left: Modulus of the matrix coefficients. Right: Argument of the matrix coefficients (degrees). The truncation degree of the set of SPHs is set to $L = 4$. The SPHs are sorted in ascending order of degrees ($l = 0, 1, 2, 3, 4$), with orders ($-l \le m \le l$) grouped together.</figcaption></figure>

The Wigner D-functions introduced in Eqs. (E.1) and (E.2) are the matrix elements of the rotation operator,

$$D = \begin{bmatrix} D_{0,0}^0 & 0 & \cdots & 0 \\ 0 & D_{q,m}^1 & \cdots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & \cdots & D_{q,m}^L \end{bmatrix}, \quad (E.5)$$

where $L$ is the chosen truncation degree, and $D_{q,m}^l$ is the matrix defined as

$$D_{q,m}^l = \begin{bmatrix} D_{-l,-l}^l & \cdots & D_{-l,m}^l & \cdots & D_{-l,l}^l \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ D_{q,-l}^l & \cdots & D_{q,m}^l & \cdots & D_{q,l}^l \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ D_{l,-l}^l & \cdots & D_{l,m}^l & \cdots & D_{l,l}^l \end{bmatrix}. \quad (E.6)$$

Thus, the set of SPHs associated with the rotated coordinate system, denoted by $Y = [Y_1, \dots, Y_K]^T$, is expressed as a function of the initial set, $\hat{Y} = [\hat{Y}_1, \dots, \hat{Y}_K]^T$, through the algebraic relation

$$Y = D^T \hat{Y}, \quad (E.7)$$

which is analogous to the change of basis equation given by Eq. (D.4). The inverse transformation is given by

$$\hat{Y} = \overline{D} Y. \quad (E.8)$$

We remark that $D^T \overline{D} = I$, the symbol $I$ referring to the identity matrix. As a consequence, the rotation transformation is conservative: going back and forth through the transformation does not induce any loss of information. Figure E.1 shows the rotation matrix obtained for a rotation of Euler angles $(\alpha, \beta, \gamma) = (120^{\circ}, 60^{\circ}, 270^{\circ})$ with a truncation degree set to $L = 4$. As previously observed in Eq. (E.1), the rotated SPHs are expressed in terms of the non-rotated SPHs having the same degree ($l$). The orders ($m$) only are mixed up by the rotation of the coordinate system.

## Appendix F: Transition matrices

In this appendix, we establish the general transition relations between sets of basis functions of various types (SPHs, SCHNs, SCHDs; see Sect. 2.4 and Appendix C), and associated with different coordinate systems. These relations are used to compute the overlap coefficients between the forcing gravitational tidal potential and the tidal response in the calculation of the tidal torque.

Let be two sets of basis functions, $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$, associated with the same coordinate system, $\mathcal{R}$, but not necessarily the same definition domains$^4$, denoted by $C_{\hat{\mathcal{A}}}$ and $C_{\hat{\mathcal{B}}}$, respectively. Both $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ can be either the SPHs, SCHNs, or SCHDs, whose expressions are given by Eqs. (C.9) and (C.19). The sets undergo the rotation transformations $\hat{\mathcal{A}} \to \mathcal{A}$ and $\hat{\mathcal{B}} \to \mathcal{B}$, where $\mathcal{A}$ and $\mathcal{B}$ are assumed to have the same definition domains as $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$, respectively. The expressions of the standard Euler rotation matrices $R_{\mathcal{A}}$ and $R_{\mathcal{B}}$ describing these transformations are given by Eq. (D.1), and the corresponding Wigner rotation matrices for the functional bases, $D_{\mathcal{A}}$ and $D_{\mathcal{B}}$, are given by Eq. (E.5). We note that the second Euler angle can only be set to $\beta = 0$ in the case of the SCHs since the definition domains of the initial and rotated function would not overlap otherwise in that case. However, the rotation matrix of the SPHs still holds for the rotations of angles $\alpha$ and $\gamma$ since these angles just introduce a phase lag in the $e^{im\varphi}$ component of the SCHs.

The initial and rotated functional basis vectors are denoted by $Y_{\hat{\mathcal{A}}}, Y_{\hat{\mathcal{B}}}, Y_{\mathcal{A}},$ and $Y_{\mathcal{B}}$, the vector $Y_{\mathcal{H}}$ of any set $\mathcal{H}$ being expressed as

$$Y_{\mathcal{H}} = [Y_{\mathcal{H},1}, Y_{\mathcal{H},2}, \dots, Y_{\mathcal{H},k}, \dots, Y_{\mathcal{H},K_{\mathcal{H}}}]^T, \quad (F.1)$$

where the $Y_{\mathcal{H},k}$ designate the orthogonal functions of set $\mathcal{H}$, and $K_{\mathcal{H}}$ its total number of elements, which depends on the truncation degree $L_{\mathcal{H}}$ introduced in Eq. (E.5). To lighten expressions, we assume that the sets of basis functions $\hat{\mathcal{A}}$ and $\mathcal{A}$ contain $M$ elements, and that $\hat{\mathcal{B}}$ and $\mathcal{B}$ contain $N$ elements, so that $K_{\hat{\mathcal{A}}} = K_{\mathcal{A}} = M$ and $K_{\hat{\mathcal{B}}} = K_{\mathcal{B}} = N$. Using the above notations, we can write down the relations

$$\begin{align} Y_{\mathcal{A}} &= D_{\mathcal{A}}^{\mathrm{T}} Y_{\hat{\mathcal{A}}}, & Y_{\mathcal{B}} &= D_{\mathcal{B}}^{\mathrm{T}} Y_{\hat{\mathcal{B}}}, \nonumber \\ Y_{\hat{\mathcal{A}}} &= \overline{D_{\mathcal{A}} Y_{\mathcal{A}}}, & Y_{\hat{\mathcal{B}}} &= \overline{D_{\mathcal{B}} Y_{\mathcal{B}}}. \tag{F.2} \end{align}$$

<footnote>$^4$By definition domain, we refer to the domain where the functions are non-zero, the SCHs being defined on the full unit sphere similarly as the SPHs in the used convention (see Eq. (C.19)).</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. F.1. Transition matrices between the SCHNs, SCHDs and SPHs for an ocean basin of angular radius $\theta_0 = 130^\circ$. The transition matrices are defined in Eq. (F.10). The SCHs (SCHNs and SCHDs) are calculated from the expressions given by Eq. (C.19), and the SPHs from the expression given by Eq. (C.9), using the same coordinate system for all the sets of basis functions and the truncation degree $L = 4$. In every set, the basis functions are sorted in ascending order of orders ($m = -4, -3, \dots, 3, 4$), with degrees ($l \ge |m|$) grouped together.</figcaption></figure>

As a next step, we establish the relationship between the bases $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$, which share the same coordinate system. The functions $Y_{\hat{\mathcal{B}},j}$ and $Y_{\hat{\mathcal{A}},k}$ are expressed in terms of each other through the relations

$$Y_{\hat{\mathcal{B}},j} = \sum_{k=1}^{M} \frac{\langle Y_{\hat{\mathcal{A}},k}, Y_{\hat{\mathcal{B}},j} \rangle}{\langle Y_{\hat{\mathcal{A}},k}, Y_{\hat{\mathcal{A}},k} \rangle} Y_{\hat{\mathcal{A}},k} + \varepsilon_{\hat{\mathcal{B}},j}, \quad (F.3)$$

$$Y_{\hat{\mathcal{A}},k} = \sum_{j=1}^{N} \frac{\langle Y_{\hat{\mathcal{B}},j}, Y_{\hat{\mathcal{A}},k} \rangle}{\langle Y_{\hat{\mathcal{B}},j}, Y_{\hat{\mathcal{B}},j} \rangle} Y_{\hat{\mathcal{B}},j} + \varepsilon_{\hat{\mathcal{A}},k}. \quad (F.4)$$

In the above equations, the functions $\varepsilon_{\hat{\mathcal{B}},j}$ and $\varepsilon_{\hat{\mathcal{A}},k}$ are the residuals of the functions $Y_{\hat{\mathcal{B}},j}$ or $Y_{\hat{\mathcal{A}},k}$ when the latter are expanded in terms of the functions of set $\hat{\mathcal{A}}$ or $\hat{\mathcal{B}}$, respectively. These residuals result both from the truncation of the sets $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ – which are mathematically defined as infinite-dimension sets of basis functions for the definition domains $C_{\hat{\mathcal{A}}}$ and $C_{\hat{\mathcal{B}}}$ –, and from the fact that $C_{\hat{\mathcal{A}}}$ and $C_{\hat{\mathcal{B}}}$ are different from each other in the general case.

For instance, if $C_{\hat{\mathcal{A}}} \subset C_{\hat{\mathcal{B}}}$, the functions of set $\hat{\mathcal{B}}$ can be expressed as linear combinations of the elements of set $\hat{\mathcal{A}}$ on $C_{\hat{\mathcal{A}}}$ only. This is the case for SPHs, which cannot be approximated using SCHs while the opposite is true. Therefore, the

function $Y_{\hat{\mathcal{B}},j}$ tends to be fully described by the functions of set $\hat{\mathcal{A}}$ as $M \to \infty$ if $C_{\hat{\mathcal{B}}} \subset C_{\hat{\mathcal{A}}}$, and consequently $\|\varepsilon_{\hat{\mathcal{B}},j}\|_2 \to 0$ in that case. Conversely, $\|\varepsilon_{\hat{\mathcal{B}},j}\|_2 > 0$ for $M = \infty$ in the general case if $C_{\hat{\mathcal{A}}} \subset C_{\hat{\mathcal{B}}}$. The residual function $\varepsilon_{\hat{\mathcal{A}},k}$ behaves in a similar way if the domains $C_{\hat{\mathcal{A}}}$ and $C_{\hat{\mathcal{B}}}$ are interchanged: $\|\varepsilon_{\hat{\mathcal{A}},k}\|_2 \to 0$ as $N \to \infty$ if $C_{\hat{\mathcal{A}}} \subset C_{\hat{\mathcal{B}}}$, while $\|\varepsilon_{\hat{\mathcal{A}},k}\|_2 > 0$ if $C_{\hat{\mathcal{B}}} \subset C_{\hat{\mathcal{A}}}$.

Since the functions of the two sets are assumed to be orthonormal by construction (see Appendix C), $\langle Y_{\hat{\mathcal{A}},k}, Y_{\hat{\mathcal{A}},k} \rangle = \langle Y_{\hat{\mathcal{B}},j}, Y_{\hat{\mathcal{B}},j} \rangle = 1$ for all $k \in \{1, \dots, M\}$ and $j \in \{1, \dots, N\}$. As a consequence, the expressions given by Eqs. (F.3) and (F.4) simplify to

$$Y_{\hat{\mathcal{B}},j} = \sum_{k=1}^{M} \langle Y_{\hat{\mathcal{A}},k}, Y_{\hat{\mathcal{B}},j} \rangle Y_{\hat{\mathcal{A}},k} + \varepsilon_{\hat{\mathcal{B}},j} \quad \text{on } C_{\hat{\mathcal{B}}}, \quad (F.5)$$

$$Y_{\hat{\mathcal{A}},k} = \sum_{j=1}^{N} \langle Y_{\hat{\mathcal{B}},j}, Y_{\hat{\mathcal{A}},k} \rangle Y_{\hat{\mathcal{B}},j} + \varepsilon_{\hat{\mathcal{A}},k} \quad \text{on } C_{\hat{\mathcal{A}}}, \quad (F.6)$$

which, introducing the vectors

$$\varepsilon_{\hat{\mathcal{B}}} = [\varepsilon_{\hat{\mathcal{B}},1}, \dots, \varepsilon_{\hat{\mathcal{B}},j}, \dots, \varepsilon_{\hat{\mathcal{B}},N}]^T, \quad (F.7)$$

$$\varepsilon_{\hat{\mathcal{A}}} = [\varepsilon_{\hat{\mathcal{A}},1}, \dots, \varepsilon_{\hat{\mathcal{A}},k}, \dots, \varepsilon_{\hat{\mathcal{A}},M}]^T, \quad (F.8)$$
<figure><img src="image_13.png" /><figcaption>Fig. F.2. Real part of complex SCH of order m = 1 and degree l₃ = 6.806057 corresponding to an ocean basin of angular radius θ₀ = 80° (i.e. θ_c = 100°) with Dirichlet boundary conditions (Eq. (C.16)). Left: Function plotted in its coordinate system using the expression given by Eq. (C.19). Middle: Function expanded in SPHs up to the truncation degree L = 6, 8, 16 (from top to bottom) in the same coordinate system using the transition relation of Eq. (F.9). Right: Function expanded in rotated SPHs (α = 120°, β = 60°, γ = 270°) following the relation given by Eq. (F.11). All functions are plotted in polar coordinates, the radial coordinate corresponding to the colatitude of the considered coordinate system, and the angular coordinate to its longitude (in degrees).</figcaption></figure>

may be rewritten as

$$Y_{\hat{B}} = P_{\hat{A}, \hat{B}}^{T} Y_{\hat{A}} + \varepsilon_{\hat{B}}, \qquad Y_{\hat{A}} = \overline{P_{\hat{A}, \hat{B}}} Y_{\hat{B}} + \varepsilon_{\hat{A}}. \qquad (F.9)$$

In the above equations, the transition matrix is expressed as

$$P_{\hat{A}, \hat{B}} \equiv \begin{bmatrix} \langle Y_{\hat{A},1}, Y_{\hat{B},1} \rangle & \dots & \langle Y_{\hat{A},1}, Y_{\hat{B},j} \rangle & \dots & \langle Y_{\hat{A},1}, Y_{\hat{B},N} \rangle \\ \vdots & & \vdots & & \vdots \\ \langle Y_{\hat{A},k}, Y_{\hat{B},1} \rangle & \dots & \langle Y_{\hat{A},k}, Y_{\hat{B},j} \rangle & \dots & \langle Y_{\hat{A},k}, Y_{\hat{B},N} \rangle \\ \vdots & & \vdots & & \vdots \\ \langle Y_{\hat{A},M}, Y_{\hat{B},1} \rangle & \dots & \langle Y_{\hat{A},M}, Y_{\hat{B},j} \rangle & \dots & \langle Y_{\hat{A},M}, Y_{\hat{B},N} \rangle \end{bmatrix}, \qquad (F.10)$$

where $\langle Y_{\hat{A},k}, Y_{\hat{B},j} \rangle$ is the scalar product of the k-th element of $\hat{A}$ with the j-th element of $\hat{B}$ following the formulation given

by Eq. (C.12). We note that the integral of the scalar product is actually performed on the intersection of the domains associated with $\hat{A}$ and $\hat{B}$, namely $C = C_{\hat{A}} \cap C_{\hat{B}}$.

With the used normalisations, the transition matrix satisfies the properties $P_{\hat{B},\hat{A}} = \overline{P_{\hat{A},\hat{B}}^T}$, $P_{\hat{A},\hat{A}} = I_M$, and $P_{\hat{B},\hat{B}} = I_N$, the notation $I_M$ referring to the identity matrix of size M. However, one should pay attention to the fact that $P_{\hat{A},\hat{B}}P_{\hat{B},\hat{A}} \neq P_{\hat{A},\hat{A}}$ and $P_{\hat{B},\hat{A}}P_{\hat{A},\hat{B}} \neq P_{\hat{B},\hat{B}}$ due to the inability of the sets of basis functions to fully describe the space of functions defined on the unit sphere, as discussed above. In other words, going back and forth between the two sets necessarily induces a loss of information contrary to the rotation operator of the SPHs detailed in Appendix E, which is conservative.

Figure F.1 shows the transition matrices $P_{\hat{A},\hat{B}}$ between the SCHNs, SCHDs, and SPHs for an ocean of angular radius $\theta_0 =$
130° (i.e. a continent of angular radius $\theta_c = 50°$). The basis functions are sorted in ascending order of degrees and orders in every set, with degrees grouped together, which differs from the SPHs of Fig. E.1, where orders are grouped together. Only basis functions with the same order overlap, due to the fact that the SCHs and SPHs have the same longitudinal component. However, some functions belonging to one set are poorly described by the basis functions of the other sets.

By combining together the relations given by Eqs. (F.2) and (F.9), we finally obtain

$$Y_B = P_{A,B}^T Y_A + D_B^T \epsilon_B, \quad Y_A = \overline{P_{A,B}} Y_B + D_A^T \epsilon_A, \quad (F.11)$$

where the transition matrix $P_{A,B}$ is given by

$$P_{A,B} = \overline{D_A}^T P_{A,B} D_B. \quad (F.12)$$

A given function $f$ defined on the sphere can thus be written both in terms of the basis functions of $\mathcal{A}$ and those of $\mathcal{B}$ as

$$f = y_{\mathcal{A}}^{T} Y_{\mathcal{A}} + \epsilon_{\mathcal{A},f}, \quad f = y_{\mathcal{B}}^{T} Y_{\mathcal{B}} + \epsilon_{\mathcal{B},f}, \quad (F.13)$$

where we have introduced the residual functions $\epsilon_{\mathcal{A},f}$ and $\epsilon_{\mathcal{B},f}$, and the coordinate vectors of $f$ in the two sets of basis functions, $y_{\mathcal{A}}$ and $y_{\mathcal{B}}$, defined as

$$y_{\mathcal{A}} = [y_{\mathcal{A},1}, \dots, y_{\mathcal{A},k}, \dots, y_{\mathcal{A},M}]^T, \quad (F.14)$$

$$y_{\mathcal{B}} = [y_{\mathcal{B},1}, \dots, y_{\mathcal{B},j}, \dots, y_{\mathcal{B},N}]^T. \quad (F.15)$$

The vectors $y_{\mathcal{A}}$ and $y_{\mathcal{B}}$ contain the coefficients of the expansion of $f$ on the set of basis functions $\mathcal{A}$ or $\mathcal{B}$, respectively. They are related to each other by the transition matrix $P_{\mathcal{A},\mathcal{B}}$ introduced in Eq. (F.12) through the relation equations

$$y_{\mathcal{A}} = P_{\mathcal{A},\mathcal{B}} y_{\mathcal{B}}, \quad y_{\mathcal{B}} = \overline{P_{\mathcal{A},\mathcal{B}}}^T y_{\mathcal{A}}. \quad (F.16)$$

Figure F.2 shows the three steps of the transition from SCHs to rotated SPHs through the example of the eigenfunction $\Theta_{l_3}^1$ obtained with Dirichlet conditions (SCHD). First the function is evaluated in its associated coordinate system using the expression given by Eq. (C.19) (left panels). We note that the region where the function is set to zero is not represented in the plot. Second, the function is expanded in series of SPHs in the same coordinate system (middle panels). Finally, we apply an Euler rotation to the SPHs, as described by Eq. (E.1), and $\Theta_{l_3}^1$ is evaluated in terms of the rotated SPHs in the rotated coordinate system (right panel). These three steps are shown for various truncation degrees of the SPHs, $L = 6, 8, 16$ (from top to bottom). We observe that the discontinuity of the gradient of the function at $\theta = \theta_0$ gives rise to 'ringing' when $\Theta_{l_3}^1$ is expanded in series of SPHs. This effect also alters the expansions of the SCHNs due to their discontinuous transition at $\theta = \theta_0$. However it vanishes as the truncation degrees $L$ increases.

## Appendix G: Gyroscopic coefficients

In this section, we detail the calculation of the gyroscopic coefficients that couple the LTEs. Following earlier studies (Webb 1980, 1982; Farhat et al. 2022a), these coefficients are denoted by $\beta$. They are defined as

$$\beta_{k,j} = - \int_C \cos \hat{\theta} \mathbf{e}_r \cdot (\overline{\nabla \Phi_k} \times \nabla \Phi_j) dS, \quad (G.1)$$

$$\beta_{k,-j} = \int_C \cos \hat{\theta} (\overline{\nabla \Phi_k} \cdot \nabla \Psi_j) dS, \quad (G.2)$$

$$\beta_{-k,j} = - \int_C \cos \hat{\theta} (\overline{\nabla \Psi_k} \cdot \nabla \Phi_j) dS, \quad (G.3)$$

$$\beta_{-k,-j} = - \int_C \cos \hat{\theta} \mathbf{e}_r \cdot (\overline{\nabla \Psi_k} \times \nabla \Psi_j) dS, \quad (G.4)$$

where the $\Phi_k$ and $\Psi_k$ are the basis functions generating the potential and stream functions, respectively (see Sect. 2.4 and Appendix C).

We emphasise that $k$ and $j$ are the indices of the functions in their respective sets. One shall be cautious not to mix up these indices with the orders and degrees of the functions. The third coefficient, given by Eq. (G.3), is straightforwardly deduced from the second one, given by Eq. (G.2), since $\beta_{-k,j} = -\beta_{j,-k}$. In some particular cases, such as the hemispherical ocean configuration considered in previous studies (e.g. Webb 1980; Farhat et al. 2022a), the recurrence relations of the ALFs can be used to obtain the analytical expressions of the gyroscopic coefficients (see e.g. Farhat et al. 2022a, Appendix G). Unfortunately, these recurrence relations do not hold for non-integral degrees, which requires to evaluate the integrals of Eqs. (G.1–G.4) numerically. However, it is still possible to write down these integrals as products of separated coordinates.

Introducing the coordinates of the north pole in the ocean's coordinate system, $(\theta_{NP}, \varphi_{NP}) = (\hat{\theta}_{oc}, \pi)$, and using trigonometric identities, we express $\cos \hat{\theta}$ as a function of the colatitude of the ocean's centre $(\hat{\theta}_{oc})$ and the ocean's coordinate system,

$$\cos \hat{\theta} = \cos \theta_{NP} \cos \theta + \sin \theta_{NP} \sin \theta \cos(\varphi - \varphi_{NP}) \quad (G.5)$$

$$= \cos \hat{\theta}_{oc} \cos \theta - \sin \hat{\theta}_{oc} \sin \theta \cos \varphi. \quad (G.6)$$

Considering the definitions of the SPHs and SCHs given by Eqs. (C.9) and (C.19), the functions $\Phi_k$ and $\Psi_k$ are expressed in the general case as

$$\Phi_k(\theta, \varphi) = \alpha_k F_k(\theta) e^{im_k \varphi}, \quad (G.7)$$

$$\Psi_k(\theta, \varphi) = \gamma_k G_k(\theta) e^{im_k \varphi}, \quad (G.8)$$

where $F_k(\theta) = P_{l_k}^{m_k}(\cos \theta)$ designates the ALF of the $k$-th function of the set of SCHNs (Neumann conditions) used to expand the potential function, and $G_k(\theta) = P_{l_k}^{m_k}(\cos \theta)$ the ALF of the $k$-th function of the set of SCHDs (Dirichlet conditions) used to expand the stream function, as detailed in Appendix C. The gyroscopic coefficients given by Eqs. (G.1–G.4) are therefore expressed as

$$\beta_{k,j} = -i2\pi\alpha_k\alpha_j \left\{ \begin{aligned} & \cos \hat{\theta}_{oc} \delta_{m_k,m_j} m_k \int_0^{\theta_0} \cos \theta \left[ \frac{dF_k}{d\theta} F_j + F_k \frac{dF_j}{d\theta} \right] d\theta \\ & - \frac{\sin \hat{\theta}_{oc}}{2} \delta_{|m_j-m_k|,1} \int_0^{\theta_0} \sin \theta \left[ m_j \frac{dF_k}{d\theta} F_j + m_k F_k \frac{dF_j}{d\theta} \right] d\theta \end{aligned} \right\}, \quad (G.9)$$

$$\beta_{k,-j} = \pi \alpha_j \gamma_k \left\{ \begin{aligned} & \cos \hat{\theta}_{oc} \delta_{m_k,m_j} \int_0^{\theta_0} \sin 2\theta \left[ \frac{dF_j}{d\theta} \frac{dG_k}{d\theta} + m_k^2 \frac{F_j G_k}{\sin^2 \theta} \right] d\theta \\ & - \sin \hat{\theta}_{oc} \delta_{|m_j-m_k|,1} \int_0^{\theta_0} \sin^2 \theta \left[ \frac{dF_j}{d\theta} \frac{dG_k}{d\theta} + m_k m_j \frac{F_j G_k}{\sin^2 \theta} \right] d\theta \end{aligned} \right\}, \quad (G.10)$$
<figure><img src="image_14.png" /><figcaption>**Fig. G.1.** Gyroscopic coefficients of an ocean basin of angular radius $\theta_0 = 130^\circ$ (i.e. $\theta_c = 50^\circ$) for various values of the colatitude of the ocean's centre. Top: $\hat{\theta}_{oc} = 0^\circ$ (polar ocean). Middle: $\hat{\theta}_{oc} = 45^\circ$ (mid-latitude ocean). Bottom: $\hat{\theta}_{oc} = 90^\circ$ (equatorial ocean). The gyroscopic coefficients are computed from the expressions given by Eqs. (G.9–G.12) with the truncation degree set to $L = 4$. The basis functions of indices $k$ and $j$ are sorted in ascending order of orders ($m = -4, -3, \dots, 3, 4$), with degrees ($l \ge |m|$) grouped together. For convenience, the minimum and maximum of the colour bar are set to 6 in absolute value although the maximum values reached by gyroscopic coefficients are around ~15 in this case.</figcaption></figure>

$$ \beta_{-k,j} = -\beta_{j,-k}, \qquad (G.11) $$

$$ \beta_{-k,-j} = -i2\pi\gamma_k\gamma_j \left\{ \begin{aligned} & \cos \hat{\theta}_{oc} \delta_{m_k,m_j} m_k \int_0^{\theta_0} \cos \theta \left[ \frac{dG_k}{d\theta} G_j + G_k \frac{dG_j}{d\theta} \right] d\theta \\ & - \frac{\sin \hat{\theta}_{oc}}{2} \delta_{|m_j - m_k|, 1} \int_0^{\theta_0} \sin \theta \left[ m_j \frac{dG_k}{d\theta} G_j + m_k G_k \frac{dG_j}{d\theta} \right] d\theta \end{aligned} \right\}. \qquad (G.12) $$

Interestingly, the coefficients $\beta_{\pm k, \pm j}$ are zero if $|m_j - m_k| > 1$, meaning that Coriolis effects couple the eigenmodes having either the same order or neighbouring orders in the general case. If the ocean is polar ($\hat{\theta}_{oc} = 0, \pi$), the gyroscopic coefficients only mix up eigenmodes with the same order ($m_j = m_k$). In this configuration, $\beta_{\pm k, \pm j}$ ($\hat{\theta}_{oc} = \pi$) = $-\beta_{\pm k, \pm j}$ ($\hat{\theta}_{oc} = 0$). If the ocean is equatorial ($\hat{\theta}_{oc} = \pi/2$), the gyroscopic coefficients only mix up eigenmodes of neighbouring orders ($|m_j - m_k| = 1$). Between these two extremal configurations, the gyroscopic coefficients

are the sum of one component in $\cos \hat{\theta}_{oc}$ associated with eigenmodes having the same order, and one component in $\sin \hat{\theta}_{oc}$ associated with eigenmodes having neighbouring orders. The gyroscopic coefficients that do not mix up potential and stream functions ($\beta_{k,-j}$ and $\beta_{-k,-j}$) are pure imaginary numbers. Conversely, the coefficients that mix up the functions ($\beta_{k,-j}$ and $\beta_{-k,j}$) are pure real numbers.

By comparing the expressions given by Eqs. (G.9) and (G.12), we remark that the $\beta_{k,j}$ and $\beta_{-k,-j}$ converge towards the same asymptotic values as $\theta_0 \to 180^\circ$ since both the SCHNs and SCHDs converge towards the standard SPHs for the $l^2$-norm introduced in Eq. (C.13), as illustrated by Fig. C.2. To our knowledge, no recurrence formula exist between the ALFs of real degrees of the SCHNs and SCHDs except in the specific case of the hemispherical ocean, where the degrees are integers (see Webb 1980, 1982; Farhat et al. 2022a, and Table C.1). As a consequence, the integrals in Eqs. (G.9–G.12) are evaluated numerically in practice.

Figure G.1 shows the gyroscopic coefficients characterising an ocean basin of angular radius $\theta_0 = 130^\circ$ for three typical configuration: $\hat{\theta}_{oc} = 0^\circ$ (polar ocean), $\hat{\theta}_{oc} = 45^\circ$ (mid-latency
ocean), and $\hat{\theta}_{oc} = 90^{\circ}$ (equatorial ocean). The eigenmodes of indices $k$ and $j$ are sorted in ascending order of degrees and orders in the displayed plots, which emphasises the fact that gyroscopic coefficients only couple eigenmodes with the same or neighbouring orders.

The coefficients $v_{k,q}$ introduced in Eq. (51) are evaluated similarly as the gyroscopic coefficients. These coefficients are defined as

$$v_{k,q} = \int_{O} \mathbf{e}_r \cdot (\overline{\tilde{\nabla} \Psi_k} \times \tilde{\nabla} Y_q) dS, \quad (G.13)$$

where the SPHs $Y_q$ are written as products of functions of separated variables,

$$Y_q(\theta, \varphi) = \eta_q H_q(\theta) e^{im_q \varphi}. \quad (G.14)$$

After an integration by part, the coefficients are simply expressed as

$$v_{k,q} = i2\pi\gamma_k\eta_q m_k \delta_{m_k,m_q} [G_k(\theta_0) H_q(\theta_0) - G_k(0) H_q(0)]. \quad (G.15)$$

We note that the $v_{k,q}$ are zero if $m_k = 0$. Besides, if $|m_k| > 0$, the Dirichlet condition satisfied by the $\Psi_k$ at $\theta_0$ implies that $G_k(\theta_0) = 0$, and the polar condition of $Y_q$ implies that $H_q(0) = 0$. As a consequence, $v_{k,q} = 0$ for all $k$ and $q$.

## Appendix H: Solid tidal Love numbers

The gravitational and load Love numbers introduced in Eqs. (47) and (48) describe the tidal response of the planet's solid part in the general case as far as the latter is spherically symmetric. For the present study, these parameters were formulated analytically from the prescriptions given by Bolmont et al. (2020) for solid Earth (see Table 2 of the article). This formulation is based on the tidal model described in Tobie et al. (2005), which integrates the equations of the elasto-gravitational theory (Takeuchi & Saito 1972) for realistic radial profiles of background quantities. As a first step, Bolmont et al. (2020) used this tidal model with the background profiles computed from the internal structure model of Sotin et al. (2007). As a second step, they fitted the parameters of simplified models to the obtained solutions. These simplified models describe the rheology of equivalent bodies with uniform interiors, for which the Love numbers can be explicitly formulated in terms of the degrees of the associated SPHs (Munk & MacDonald 1960).

The degree-$l$ Love numbers of a homogeneous solid body of mass $M_{core}$, radius $R_{core}$, and uniform shear modulus $\mu$, are expressed as (e.g. Munk & MacDonald 1960)

$$\{k_l^\sigma, h_l^\sigma, k_{L,l}^\sigma, h_{L,l}^\sigma\} = \frac{1}{1 + \tilde{\mu}_l^\sigma} \left\{ \frac{3}{2(l-1)}, \frac{2l+1}{2(l-1)}, -1, -\frac{2l+1}{3} \right\}. \quad (H.1)$$

In the above equation, the degree-$l$ dimensionless shear modulus $\tilde{\mu}_l^\sigma$ is defined as

$$\tilde{\mu}_l^\sigma \equiv A_l \tilde{\mu}, \quad (H.2)$$

where $\tilde{\mu}$ designates the dimensionless visco-elastic rigidity accounting for the rheology of the material, and $A_l$ the dimensionless parameter given by (Castillo-Rogez et al. 2011; Efraimsky 2012)

$$A_l \equiv \frac{4(2l^2 + 4l + 3)\pi R_{core}^4 \mu}{3lGM_{core}^2}. \quad (H.3)$$

Since the oceanic layer is thin in the present study, the mass and radius of the solid part can be approximated by the planet mass and radius, respectively, which yields $M_{core} \approx M_p$ and $R_{core} \approx R_p$.

The normalised rigidity $\tilde{\mu}$ is a complex parameter that describes the rheological behaviour of the solid part in the frequency domain. This parameter is determined by the adopted rheological model. For standard models such as the Maxwell visco-elastic rheology (e.g. Correia et al. 2014), it is expressed as a simple function of the tidal frequency, $\sigma$. Following Bolmont et al. (2020), we opt for an Andrade rheology. Similarly as the Maxwell rheology, the Andrade model describes the visco-elastic deformation of the material. However, it additionally accounts for the unrecoverable creep forming the anelastic component of the response, which affects the frequency-dependence of tidal dissipation in the high-frequency range (e.g. Efraimsky 2012). In the Andrade rheology, $\tilde{\mu}$ is expressed as (e.g. Castillo-Rogez et al. 2011)

$$\tilde{\mu} = \frac{1}{1 + (i\sigma\tau_A)^{-\alpha} \Gamma(1 + \alpha) + (i\sigma\tau_M)^{-1}}, \quad (H.4)$$

where $\tau_M$ designate the Maxwell time, $\tau_A$ the Andrade time, $\Gamma$ the Gamma function introduced in Eq. (B.3), and $\alpha$ the parameter determining the duration of the transient response in the primary creep (Castelnau et al. 2008; Castillo-Rogez et al. 2011).