Issue 
A&A
Volume 680, December 2023



Article Number  A13  
Number of page(s)  36  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202347301  
Published online  05 December 2023 
Can one hear supercontinents in the tides of ocean planets?
IMCCE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université,
77 Avenue DenfertRochereau,
75014
Paris, France
email: pierre.auclairdesrotour@obspm.fr
Received:
27
June
2023
Accepted:
9
October
2023
Context. Recent observations and theoretical progress made about the history of the EarthMoon system suggest that tidal dissipation in oceans primarily drives the longterm evolution of orbital systems hosting ocean planets. Particularly, they emphasise the key role played by the geometry of landocean 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 selfattraction 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 nonpolar configurations, and between smallsized and mediumsized 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 / planetstar interactions / planets and satellites: oceans / planets and satellites: terrestrial planets / Earth
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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^{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 paleooceans 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 superEarths GJ 1214 b (e.g. Charbonneau et al. 2009), LHS 1140 b (e.g. LilloBox et al. 2020), and TOI1452 b (Cadieux et al. 2022), or the Earthsized planets TRAPPIST1 e, f, and g (e.g. Grimm et al. 2018), with the latter being suspected to harbour oceanscale 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 longterm 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; AuclairDesrotour et al. 2019b), they result from the mutual gravitational interactions between 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 frequencyresonant surface gravity modes distorted by the planet’s rotation (e.g. AuclairDesrotour et al. 2018), which strongly differs from the smooth viscoelastic 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 EarthMoon system, using either semianalytical approaches with simplified geometries (e.g. Webb 1980, 1982; Bills & Ray 1999; Farhat et al. 2022a)  including global ocean models (e.g. AuclairDesrotour et al. 2018; Motoyama et al. 2020; Tyler 2021) –, or numerical methods based upon realistic landocean 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 yearhistory of the Earth’s length of the day (LOD) and EarthMoon 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 largescale 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 semianalytical tidal theory of hemispherical oceans (LonguetHiggins & 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 socalled ‘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 frequencybehaviour of an Earthsized 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 ≪ 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 θ_{c}, as shown by Fig. 1. Accordingly, the angular radius of the ocean basin is defined as θ_{0} ≡ π − θ_{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 ℛ: (O, e_{X}, e_{Y}, 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 (e_{X}, e_{Y}, e_{Z}) is such that e_{X} and e_{Y} correspond to two orthogonal directions in the planet’s equatorial plane, and e_{Z} to the direction of the spin vector. The unit vector e_{Z} thus designates the position of the north pole on the unit sphere. The geocentric frame of reference ℛ is associated with the usual spherical coordinates (, , ), where , , and 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, (, ). As the centre of the oceanic basin corresponds to the antipodal point, its coordinates on the globe are given by .
The oceanic frame of reference, ℛ_{oc}: (O, e_{x}, e_{y}, e_{z}), is such that e_{z} points towards the centre of the ocean basin, while e_{x} and e_{y} define two orthogonal axes in the plane normal to e_{z}. This frame of reference is associated with the spherical coordinates (r, θ, φ), where is the radial coordinate, θ the colatitude (θ = 0 at the centre of the ocean basin), and φ the longitude. The change of basis vectors (e_{X}, e_{Y}, e_{Z}) → (e_{x}, e_{y}, e_{z}) is expressed as a function of the coordinates of the continental centre as (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 in the following, similar to the example shown by Fig. 1, which yields (2)
Finally, we introduce the set of unit vectors (e_{r}, e_{θ}, e_{ϕ}) associated with the spherical coordinates (r, θ, ϕ), respectively. The change of basis vectors (e_{x}, e_{y}, e_{z}) → (e_{r}, e_{θ}, e_{ϕ}) is given by the standard relation, (3)
Fig. 1 Geometry of the studied system and associated parameters. The unit vector , which indicates the position of the continental centre, is defined as , with e_{z} pointing towards the centre of the ocean basin. In the diagram, the longitude of the continental centre is set to . 
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 ℛ, as (4)
where the constant component and the component responsible for the Keplerian dynamics of the twobody system are removed (second and third terms in the righthand member of Eq. (4)). In this equation, G is the universal gravitational constant, M_{s} the mass of the perturber, r_{s} its position vector, and r_{s} ≡ r_{s} the planetperturber distance. We note that the symbol ≡ 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 ≡ ∇U_{T}. Following the formalism introduced in earlier works (e.g. Tyler 2011; Matsuyama 2014; AuclairDesrotour et al. 2018, 2019a; Motoyama et al. 2020; Farhat et al. 2022a), we write down the momentum and mass conservation equations, respectively, as (5) (6)
with t designating the time, g the surface gravity at rest, σ_{R} the Rayleigh drag frequency used to describe the action of dissipative mechanisms, V the horizontal velocity vector – defined from the horizontal displacement vector ξ as V ≡ ∂_{t}ξ –, f the Coriolis parameter, ζ the vertical displacement of the ocean’s surface with respect to the oceanic floor, and ζ_{eq} ≡ U_{T}/g the equilibrium displacement corresponding to the equipotential surface induced by the tidal gravitational potential.
In the momentum equation, the notations Γ_{D} and Γ_{T} refer to nontrivial linear operators accounting for the effects of ocean loading, selfattraction, 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 Γ_{D} = Γ_{T} = 1 if one neglects both the tidal deformation of the solid part (infiniterigidity approximation) and the variation of selfattraction 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 ζ_{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 ℛ, (7)
where Ω > 0 is the planet’s spin angular velocity. The horizontal gradient operator and divergence of the horizontal velocity are expressed in ℛ_{oc}, and the associated basis vectors, (e_{θ}, e_{ϕ}), as (8) (9)
The Rayleigh drag term σ_{R} V in Eq. (5) accounts for the cumulated effects of dissipative mechanisms on tidal flows. The frequency σ_{R} is the inverse of the effective dissipation timescale associated with these dissipative mechanisms. It takes values around ~ 10^{−5} s^{−1} for the present day Earth (e.g. Webb 1980; Farhat et al. 2022a)^{2}.
Dimensionless control parameters determining the regime of the oceanic tidal response.
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 longwavelength surface gravity waves, V_{0}, which are defined, respectively, as (10)
Introducing the normalised time , Coriolis parameter , horizontal gradient , the complex horizontal displacement vector , and vertical displacements and , such that (11)
with ℜ referring to the real part of a complex number, we end up with the nondimensional complex LTEs given by (12) (13)
with . In the above equations, tidal dynamics are controlled by two dimensionless parameters, (14)
The first parameter, , 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 (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 (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, , may be regarded as an Ekman number. If , the drag does not alter the tidal response much. Conversely, characterises a frictional (or viscous) regime where inertial effects are annihilated by the strong damping associated with the drag.
In addition to and , the timederivative operators in Eq. (12) introduce a third dimensionless parameter describing the ratio of tidal forces to Coriolis forces, (15)
The notation σ 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 (subinertial regime), the acceleration term of the momentum equation can be neglected with respect to Coriolis terms. Conversely, if (superinertial regime), the flow is strongly driven by the tidal forcing and it is thus hardly distorted by the planet’s rotation. We note that is the inverse of the socalled 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 displacement vector field into curlfree and divergencefree vector fields, (16)
The curlfree and divergencefree components of are defined from the divergent displacement potential Φ and the rotational displacement streamfunction Ψ, 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 Φ and Ψ (e.g. FoxKemper 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 , where 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 and 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), (17)
The first condition of Eq. (17) means that the gradient of Φ 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 , 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 Φ and Ψ are defined to a constant, the value of ψ at the coastline is set to zero, which simplifies the second condition of Eq. (17) to Ψ = 0. Interestingly, this condition induces the orthogonality of the curlfree and divergencefree components of the horizontal displacement (e.g. Farhat et al. 2022a, Appendix E), formulated as (18)
with dS = sin θdθdϕ being an infinitesimal surface element of the unit sphere, 𝒪 the domain occupied by the ocean basin, and the complex conjugate of . 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 (𝒪), and they satisfy either Neumann or Dirichlet conditions at its boundary, ∂𝒪. As a consequence, the two functions can be expanded in terms of the complete sets of orthogonal eigenfunctions {Φ_{k}}_{1≤k≤∞} or {Φ_{k}}_{1≤k≤∞} that are the solutions of the wave equations given by (19) (20)
where designates the normalised Laplacian operator defined, for any function f, as (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 Φ_{k} are associated with the eigenvalues λ_{k} ∈ ℝ^{+}, and the eigenfunctions ψ_{k} with the eigenvalues v_{k} ∈ ℝ^{+}, 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 (±_{imϕ}), 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 {Φ_{k}}_{1≤k≤∞} is refered to as the ‘SCHNs’ (Neumann condition), and the set {Φ_{k}}_{1≤k≤∞} as the ‘SCHDs’ (Dirichlet condition). For example, Fig. 2 shows the first basis functions of the two sets for a supercontinent of angular radius θ_{c} = 10° (i.e. θ_{0} = 170°).
One shall notice that the SCHs and their associated eigenvalues both depend on the boundary condition applied at the coastline, meaning that the sets {Φ_{k}}_{1≤k≤∞} and {Φ_{k}}_{1≤k≤∞} are necessarily different by construction. In each set, the basis functions are orthogonal to each other and normalised so that (22)
where δ_{k,j} designates the Kronecker delta function, such that δ_{k,j} = 1 for k = j and δ_{k,j} = 0 otherwise. However, two functions belonging to different sets are not orthogonal in the general case, (23)
since the sets of basis functions {Φ_{k}}_{1≤k≤∞} and {Φ_{k}}_{1≤k≤∞} are not the same.
Introducing the timedependent coefficients p_{k} and p_{−k}, we write down the divergent potential function and the streamfunction as (24) (25)
By substituting Eq. (24) into Eq. (13), the continuity equation becomes (26)
However, such a simple expression cannot be obtained for the equilibrium displacement () in the general case since it is 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 only in the particular case of the hemispheric ocean (θ_{0} = 90°), where the set of eigenfunctions {Φ_{k}}_{1≤k≤∞} actually corresponds to a subset of the SPHs for the hemispherical domain, as shown by Farhat et al. (2022a).
Fig. 2 Eigenfunctions describing the oceanic tidal response for a supercontinent of angular radius θ_{c} = 10° (i.e. θ_{0} = 170°), and associated tidal flows. Top: set {Φ_{k}} (SCHNs). Bottom: set {Ψ_{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,…, 3 (vertical axis) and −n ≤ m ≤ n (horizontal axis). Bright or dark colours designate positive or negative values of the eigenfunctions, respectively. Streamlines indicate the tidal flows corresponding to for the set {Φ_{k}} and to for the set {Ψ_{k}}. 
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 Φ and Ψ 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 , and we integrate it over the domain of the ocean basin. It follows (28)
The dot product of two eigenfunctions that appears in the first term of Eq. (28) is computed by invoking, successively, Green’s first identity (e.g. Strauss 2007, Chapter 7), (29)
the impermeability boundary condition on φ_{k} given by Eq. (17), the fact that the φ_{k} are eigenfunctions of the Laplacian operator as described by Eq. (19), and the orthogonality property given by Eq. (22). This yields (30)
and for the forcing term, (31)
Since the tidal gravitational potential is expressed in terms of the SPHs associated with the coordinates (, ), the calculation of the above integral requires computing the transition matrix between these SPHs and the oceanic eigenfunctions. This calculation is achieved in two steps. First, one computes the rotated SPHs associated with the change of coordinates , as detailed in Appendices D and E. Second, the transition matrices between the nonrotated SPHs and the SCHs are evaluated following the method described in Appendix F.
The integral of the second term in the lefthand member of Eq. (28) actually vanishes owing to the orthogonality property of the curlfree and divergencefree components of the horizontal displacement (Eq. (18)), (32)
Finally, the integrals of the third and fourth terms of Eq. (28) are simplified with the help of vectorial identities, yielding (33)
The above simplifications lead to the first set of temporal equations for the coefficients p_{j} and p−j, (34)
The second set of equations is obtained by repeating the same steps with instead of in the dot product operation, (35)
which, using the orthogonality properties of the eigenfunctions and the relations (36)
Thus, following the formalism used in earlier studies (Webb 1980, 1982; Farhat et al. 2022a), Eqs. (34) and (37) form an infinite linear system in the coefficients p_{k} and p_{−k}, written as (38) (39)
In the above equations, the symbols β_{k,j} designate the socalled ‘gyroscopic coefficients’ (e.g. Proudman 1920; LonguetHiggins & Pond 1970; Webb 1980), (40) (41) (42) (43)
These coefficients account for the coupling effect of the Coriolis terms and the geometry of the ocean basin, which affects how the gradients of the eigenfunctions overlap. Their behaviour and properties are examined in Appendix G. The expression of the forcing term in Eq. (39) is given here for generality. This term is actually zero, as discussed in the next section, meaning that Eq. (39) does not depend on .
2.6 Ocean loading and selfattraction variation
At this stage, we still have to express as a function of p_{j} and p_{−j}. To do so, we consider the fact that the gravitational tidal force acts on the planet as a periodic perturbation oscillating in time and longitude. Therefore, the tidal gravitational potential 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 (44)
where σ is the tidal frequency, V_{0} the reference velocity introduced in Eq. (10), the complex SPH of degree l and order m associated with the coordinate system (, ), defined in Eq. (C.9), and the associated frequencydependent normalised component of the tidal gravitational potential. In the following, these notations are shortened for convenience. Similarly as the eigenfunctions of the ocean basin, Φ_{k} and Ψ_{k}, the SPHs are simply denoted by , and the corresponding coefficients by , the index k referring to an element of the set of SPHs, . 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 σ ≥ 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, , and the associated forcing term, , expressed as (45)
As highlighted in earlier studies (e.g. Matsuyama 2014; Matsuyama et al. 2018; AuclairDesrotour et al. 2019a), the combined contribution of ocean loading, selfattraction variation, and deformation of solid regions may be formulated as (46)
where we have made use of the scalar product defined by Eq. (C.12) and introduced the solid deformation factors, and (see e.g. Hendershott 1972), (47) (48)
The factors and 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 and 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 selfattraction variation and a surface displacement. By analogy, the loading Love numbers and 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, , also depends on the ratio of seawater density to the mean density of the solid regions, . We remark that and 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 viscoelastic response of solid regions when subjected to a harmonic tidal force (Remus et al. 2012; AuclairDesrotour 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 eigen functions and the SPHs associated with the coordinate system (, ), 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, (49) (50)
where υ_{k,q} designates the coefficients defined as (51)
These coefficients are all zero as shown in Appendix G, which implies that (52)
We remark that the quadrupolar component of the tidal gravitational potential (the component associated with the SPH) is far greater than the components of higher degrees if the size of the planet is small compared with the planetperturber distance. These components can thus be ignored in standard twobody systems. The series in the formulation of the tidal gravitational given by Eq. (45) then reduces to one term only. This removes the summation over j in the second term of the righthand 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 {Φ_{j}}, {Ψ_{j}}, and 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 (53) (54) (55)
with ^{T} designating the transpose of a matrix. Converted into an algebraic form, the forcing term given by Eq. (49) is written as (56)
where we have introduced the unit vectors E_{M,k} defined as (57)
and the matrices Α_{Φ}, and R_{Φ}, (58)
In the above expressions, designates the transition matrix from the SPHs to the eigenfunctions Φ_{k}. The other matrices are defined as (59) (60) (61) (62)
Also, we introduce the matrix of gyroscopic coefficients (63)
where the block matrices Β_{φ,φ}, Β_{φ,ψ}, Β_{ψ,φ}, and Β_{ψ,ψ} are defined in terms of the β_{k,j} defined in Eqs. (40)–(43) as (64) (65) (66) (67)
Finally, the matrix accounting for the gravitational and pressure effects of the oceanic surface displacement, A, the matrix describing the coupling between the forcing tidal gravitational potential and the oceanic eigenmodes, R, and the solution vector, X, are respectively expressed as (68)
The temporal tidal equations given by Eqs. (38) and (39) thus lead to the solution (69) (70)
with 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 frequencydependent matrix (H), which links the tidal gravitational potential () to the potential and stream functions that describe the horizontal tidal displacement (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 selfattraction variations, which is derived from the gravitational potential of the tidally distorted body, U_{D}, defined at its surface as (71)
The components of this gravitational potential are expressed as (72)
where is the corresponding component of the forcing tidal potential introduced in Eq. (44), and the corresponding component of the normalised surface elevation, (73)
The first term of Eq. (72) is the selfattraction 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, (74)
which, using the expression of the tidal potential given by Eq. (72), yields (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 longterm effect of tides on the evolution of the planetperturber system is quantified by the generated tidal torques. The main contributor to this evolution is the timeaveraged tidal torque exerted about the spin axis of the planet (Zahn 1966), (76)
which affects the longterm variation rate of the planetperturber distance and the planet’s spin angular velocity (e.g. Ogilvie 2014). In the above equation, 𝒱 designates the volume filled by the planet, dV an infinitesimal volume parcel, and δρ 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 (77)
This allows the tidal torque to be expressed in terms of the forcing and deformation tidal gravitational potentials, (78)
The latter integral is defined over the volume filled by the tidally deformed planet (𝒱), 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 ∇^{2} U_{T} = 0 outside of the perturber, by virtue of Poisson’s equation, which implies that (79)
since the Laplacian and operators can be permuted. Let 𝒱_{*} be any region that is simply conne_{t}cted a_{)}nd does not include the perturber. Then, the fact that in 𝒱_{*} allows the integral of Eq. (78) to be rewritten as (80)
By virtue of Green’s second identity (e.g. Arfken & Weber 2005, Sect. 1.11), Eq. (80) yields (81)
with the notation ∂𝒱_{*} referring to the boundary of the domain 𝒱_{*}, and dS to and outward pointing infinitesimal surface element vector.
By applying the timeaveraging operator to the integral we obtain (82)
In practice, 𝒱_{*} can be chosen such that its boundary is a sphere of some radius r_{*} centred on the planet.
In the vicinity of ∂𝒱_{*}, the forcing tidal potential and the deformation tidal potential can be represented using interior and exterior multipole expansions in SPHs, respectively, (85) (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, (87)
For the quadrupolar tidal gravitational potential (e.g. Ogilvie 2014), (88)
with a being the semimajor axis of the perturber, we recover the wellknown relationship between the tidal torque and the quadrupolar Love number (e.g. Makarov & Efroimsky 2013; AuclairDesrotour et al. 2019a), (89)
The timeaveraged 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), (90)
where ρ 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 (91)
The first term of the righthand member in Eq. (91) is zero since the boundary ∂𝒱_{*} is outside of the planet, where ρ = 0. By combining together the continuity equation, (92)
and Poisson’s equation, given by Eq. (77), we express ∇ · (ρV) as a function of the deformation tidal gravitational potential, (93)
Using Green’s second identity, the above expression becomes (95)
Finally we use the expressions of and given by Eqs. (85) and (86), and the definition of the Love number given by Eq. (74). It follows (96)
The tidal input power in the ocean is expressed as a function of the harmonics of the surface tidal elevation only, (97)
By noting that and expanding the component of the forcing tidal potential in series of oceanic eigenmodes, (98)
we can rewrite the sum on degrees and orders in Eq. (97) as (99)
Therefore, owing to the orthogonality of the Φ_{k} eigenfunctions, Eq. (97) simplifies to (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 (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) (102)
The dissipated power 𝒫_{diss;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; AuclairDesrotour et al. 2018).
When the solid part is deformable, the total input tidal power given by Eq. (96), 𝒫_{T}, and the total tidally dissipated power, 𝒫_{diss}, are the sum of the solid and oceanic contributions, (103) (104)
where 𝒫_{T;sol} ≠ 𝒫_{diss;sol} and 𝒫_{T;oc} ≠ 𝒫_{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 and gravitational coupling, while 𝒫_{diss;sol} ≥ 0 and 𝒫_{diss;oc} ≥ 0. In the following, the input and tidally dissipated powers of the solid regions, 𝒫_{T;sol} and 𝒫_{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 𝒫_{diss} = 𝒫_{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.
Values of parameters used in the reference case.
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 wellchosen 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 EarthMoon system is the planetsatellite system that we know best, it appears as a privileged option for the aforementioned reference configuration. Therefore, we consider a simplified EarthMoon system where the Moon orbits an Earthsized 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 θ_{c:} = 90° and the colatitude of the continental centre to . 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 Η = 4 km and its Rayleigh drag frequency to σ_{R} = 10^{−5} s^{−1}, consistently with the values found in earlier studies from constraints on the actual tidal dissipation rate of the EarthMoon system and its age (Webb 1980, 1982; Farhat et al. 2022a). The mean density of seawater is set to ρ_{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 viscoelastic body described by the Andrade model (see Appendix H). In this model, the solid tidal response is parametrised by an effective shear modulus, μ, the Maxwell time associated with viscous friction τ_{M}, the Andrade time associated with the inelastic component of the solid elongation, τ_{A}, and the dimensionless parameter α. 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_{⦅} and mean motion n_{s} = n_{⦅}, with M_{⦅} and n_{⦅} 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 σ = 2 (Ω − 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 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.
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 χ = (Ω − n_{s}) /Ω_{Ε} (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 ( and 𝒯_{z;2} ≤ 0) and positive tidal input power (𝒫_{T} ≥ 0), while dashed lines indicate accelerating tidal torque and 𝒯_{z;2} > 0) and negative tidal input power (𝒫_{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. 
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 χ = (Ω − n_{s}) /Ω_{Ε}, with Ω_{Ε} being the spin rate of the present day Earth. For simplicity, we let the values of the spin period P_{rot} = 2π/Ω vary with the tidal frequency while the mean motion of the satellite is fixed. In reality, both Ω 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} ≪ Ω. The adopted frequency interval of the tidal forcing ranges between 0 (spinorbit synchronisation; P_{rot} ≈ 27.2 days) and 4 (fast spin rotation; P_{rot} ≈ 5.9 hr). We note that the present EarthMoon system corresponds to the normalised frequency χ ≈ 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 σ, 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 EarthMoon 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 EarthMoon system (dashed pink line at χ ≈ 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 landocean distributions (Green et al. 2017; Daher et al. 2021). In the lowfrequency range, the tidal response of the planet is dominated by the oceanic dynamical tide, which is characterised by resonant peaks. Conversely, in the highfrequency range, the resonances associated with oceanic modes are attenuated by the viscoelastic 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 ∝ σ^{−α} with α = 0.25 (see Table 2), while the nonresonant background of the oceanic tidal torque scales as ∝ σ^{−3} (e.g. AuclairDesrotour 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 solidocean energy transfer allowed by the elasticity of the solid part through ocean loading and selfattraction 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 zerofrequency limit (σ → 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 lowfrequency range. However, the solid tidal torque also tends to zero for , and 𝒯_{z} = 0 at σ = 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. AuclairDesrotour 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 (, ), this set of parameters can essentially be reduced to five numbers: three numbers governing the dynamics of oceanic tidal flows (, , ) and two numbers defining the geometry of the ocean basin (θ_{c}, ). The parameter space of the oceanic tidal response thus simplifies to five dimensions.
Moreover, for the semidiurnal tide as long as n_{s} ≪ Ω, indicating that this number can be roughly considered as invariant over time for planetsatellite systems similar to the EarthMoon system. Therefore, _{)}the parameter space can be limited to the set (, , θ_{c},) for the semidiurnal oceanic tide of a rapidly rotating Earthlike 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 , and the Rayleigh drag frequency σ_{R} for . 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 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, (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 (μ), the Maxwell time (τ_{Μ}), the Andrade time (τ_{Α}), and the dimensionless rheological parameter that accounts for the inelastic component of the solid deformation (α). 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 σ ≫ , , the imaginary part of the quadrupolar Love number defined by the Andrade model for the solid part scales as (e.g. AuclairDesrotour et al. 2019a, Eq. (47)) (106)
with being a constant, and Γ the Gamma function (Abramowitz & Stegun 1972, Chapter 6). The above expression entails that, in the limit of , , the solid tidal response is dependent on τ_{A}, but independent of τ_{M}. We thus consider only the former in the parametric study.
Similarly as the Maxwell time, the rheological parameter α 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 (CastilloRogez et al. 2011). Due to these constraints, we choose to let it fixed to the value used in the reference case, α = 0.25 (see Table 2). Nevertheless, we note from Eq. (106) that varying α would essentially modify the rheological behaviour of the solid part in the highfrequency range by acting on the scaling law . Finally, the effective shear modulus of the solid part, μ, accounts for the elasticity of the solid part, which is responsible for the solidocean coupling. Thus, as discussed for τ_{A}, we shall consider various values of μ 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, (, θ_{c}, H, σ_{R}, μ, and τ_{A}). Other parameters occur in the model, such as the planet mass (M_{p}) or radius (R_{p}), the seawater density (ρ_{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 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 (), angular radius of the continent at constant volume (θ_{c}), oceanic depth (H), Rayleigh drag frequency (σ_{R}), effective shear modulus of the solid part (μ), and Andrade time (τ_{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 , , H* = 4.0 km, , μ* = 25.1189 GPa, and .
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 , 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 (AuclairDesrotour 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 nonrotating 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, ) 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 θ_{c} ≥ 30°, 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 , 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 (e.g. AuclairDesrotour 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 closedform solutions provided by the analytical theory (e.g. AuclairDesrotour et al. 2018). The plot highlights the transition from the quasiadiabatic regime (σ_{R} = 10^{−6} s^{−1}) to the frictional regime (σ_{R} = 10^{−3} 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 quasiadiabatic regime, the frequencies of resonances are almost independent of the Rayleigh drag frequency, whereas their heights scale as , their widths as ∝ σ_{R}, and the nonresonant background of the oceanic dissipation rate as ∝ σ_{R} (AuclairDesrotour et al. 2015, AuclairDesrotour et al. 2018, 2019a). As a consequence, the number of visible peaks increases as σ_{R} decreases. The new peaks appear in the highfrequency range since they are associated with oceanic harmonics of high degrees. In the zerodrag limit (σ_{R} → 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 . This behaviour is described by a closedform solution that will be detailed in a forthcoming article. Considering the estimate of ~10^{−5} 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 quasiadiabatic and frictional regimes in Fig. 4 (middle right panel), namely σ_{R} ≈ 3.2 × 10^{−5} s^{−1} (see also Sect. 5).
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 μ = 300 GPa. The Andrade time only affects the nonresonant background of the spectrum through the scaling factor of the solid dissipation rate, as shown by Eq. (106). As τ_{A} increases, the nonresonant background decreases. This effect is visible in the highfrequency 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 landocean 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 landocean 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 ∝ 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 smallscale 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 nonlinearly 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).
Fig. 4 Frequency spectra of the imaginary part of the quadrupolar Love number associated with the semidiurnal tide, , 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 χ = (Ω − n_{s}) /Ω_{Ε} (horizontal axis), with Ω_{Ε} being the spin angular velocity of the actual Earth. The solid black line designates the reference case defined by Table 2: , , H* = 4.0 km, , μ* = 25.1189 GPa, and . 
Evolution of the frequency spectrum of the quadrupolar tidal Love number as one parameter increases while the values of other parameters are fixed.
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 socalled 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. AuclairDesrotour 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 nonresonant 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 , θ_{c}, and σ_{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 Δσ_{k} = σ_{k+1} − σ_{k} in the following, with σ_{k} being the frequency of the kth maximum. In order to make these frequency intervals insensitive to the dilationcontraction effect induced by oceanic depth variations (see Fig. 4, middle left panel), they are normalised by their average, given by (107)
where N designates the total number of intervals. We thus denote by the normalised intervals. Finally, sorting the intervals in ascending order and using the subscript jsuch that Δσ_{j} < Δσ_{j+1} for j = 1,…, N, we introduce the cumulative distribution function (CDF) of the normalised frequency intervals, (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 quasiuniformly 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, (109)
The thus defined normalised standard deviation is such that , where corresponds to uniformly spaced relative maxima or to N = 1 (two maxima), and to an extreme dispersion of the resonant peaks. By convention, we set 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 from the frequency spectra obtained for the imaginary part of the quadrupolar Love number. These values are plotted in Fig. 6.
The quantity appears to be almost insensitive to variations of H, μ, and τ_{Α}, whereas it varies over several orders of magnitude when , θ_{c} or σ_{R} are modified. Moreover, we recover quantitatively the tendencies identified in the parametric study (Table 3), namely 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 on θ_{c} or σ_{R}, is not even. Instead, threshold effects may be observed, similar to those of Fig. 5. First, the evolution of the metric with (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 axisymmetry. By switching from (polar continent) to , the metric skyrockets to reach a plateau around . This plateau suggests that changing the position of the continental centre from 10^{°} to 90^{°} 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^{°} (global ocean) to 90^{°} (hemispherical ocean) at constant water volume (Fig. 6, top right panel). While the metric does not evolve much with in general, it increases drastically around θ_{c} ≈ 35°, which is the critical size identified from the plotted spectra in the parametric study. This value is consistent with the highest value of θ_{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^{°}, 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 θ_{c} ≳ 40°, 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 Earthsized planet would be similar for a superEarth or a smaller planet.
Finally, the metric decays monotonically as the Rayleigh drag frequency increases, with a sharp variation occurring around σ_{R} ~ 3.2 × 10^{−5} s^{−1} (Fig. 6, middle right panel). This critical value marks the transition between the quasiadiabatic 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 quasiadiabatic regime, increases as σ_{R} decays due to the growing number of visible resonant peaks.
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: , , H* = 4.0 km, , μ* = 25.1189 GPa, and . 
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 landocean distribution, which predominantly drives the evolution of the EarthMoon 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 (LonguetHiggins 1968; LonguetHiggins & Pond 1970; Webb 1980, 1982; AuclairDesrotour 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 viscoelastic 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 closedform 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 EarthMoon system. 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 highfrequency 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 nonresonant 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 quasiadiabatic and frictional regimes, which occurs for σ_{R} ≈ 3.2 × 10^{−5} s^{−1}. This critical value is slightly greater than that estimated for Earth (~10^{−6}−10^{−5} s^{−1}, e.g. Garrett & Munk 1971; Webb 1973; Tyler 2021; Farhat et al. 2022a), which is consistent with the frequencyresonant 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 longterm 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 planetsatellite 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 longterm evolution of the EarthMoon 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 nonsmooth 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 EarthMoon system. We will address them through several forthcoming studies.
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: , , H^{*} = 4.0 km, , μ* = 25.1189 GPa, and . 
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 ANR19CE31000201) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo885250). This research has made use of NASA’s Astrophysics Data System.
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 
θ_{c}  Angular radius of the supercontinent  p. 2 
θ_{0}  Angular radius of the ocean basin  p. 2 
ℛ  Geocentric frame of reference rotating with the planet  p. 2 
e_{X}  Cartesian unit vector associated with ℛ  p. 2 
e_{Y}  Cartesian unit vector associated with ℛ  p. 2 
e_{Z}  Cartesian unit vector associated with ℛ  p. 2 
Radial coordinate in ℛ  p. 2  
Colatitude in ℛ  p. 2  
Longitude in ℛ  p. 2  
Colatitude of the continental centre in ℛ  p. 3  
Longitude of the continental centre in ℛ  p. 3  
Colatitude of the oceanic centre in ℛ  p. 3  
Longitude of the oceanic centre in ℛ  p. 3  
O  Planet’s centre of mass  p. 3 
ℛ_{oc}  Frame of reference rotating with the planet and oriented by the centre of the ocean basin  p. 3 
e_{x}  Cartesian unit vector associated with ℛ_{oc}  p. 3 
e_{y}  Cartesian unit vector associated with ℛ_{oc}  p. 3 
e_{z}  Cartesian unit vector associated with ℛ_{oc}  p. 3 
r  Radial coordinate in ℛ_{oc}  p. 3 
θ  Colatitude in ℛ_{oc}  p. 3 
^{φ}  Longitude in ℛ_{oc}  p. 3 
e_{r}  Radial unit vector associated with ℛ_{oc}  p. 3 
e_{θ}  Colatitudinal unit vector associated with ℛ_{oc}  p. 3 
e_{φ}  Longitudinal unit vector associated with ℛ_{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 
r_{s}  Position vector of the perturber with respect to O  p. 3 
r_{s}  Planetperturber distance  p. 3 
≡  Symbol used for definitions in equations  p. 3 
F  Tidal force per unit mass induced by U_{T}  p. 3 
t  Time  p. 3 
^{g}  Surface gravity at rest  p. 3 
σ_{R}  Rayleigh drag frequency  p. 3 
V  Horizontal velocity vector  p. 3 
ξ  Horizontal displacement vector  p. 3 
f  Coriolis parameter  p. 3 
^{ζ}  Vertical displacement of the ocean’s surface with respect to the oceanic floor  p. 3 
^{ζ}_{eq}  Equilibrium displacement corresponding to the equipotential surface  p. 3 
Γ_{D}, Γ_{T}  Solid deformation operators  p. 3 
Ω  Planet’s spin angular velocity  p. 3 
CA  Cowling approximation  p. 4 
t_{0}  Reference time for normalisation  p. 3 
V_{0}  Reference velocity for normalisation  p. 3 
Normalised time  p. 4  
Normalised Coriolis parameter  p. 4  
Normalised gradient operator  p. 4  
Normalised horizontal displacement vector  p. 4  
Normalised oceanic surface elevation  p. 4  
Normalised equilibrium surface elevation  p. 4  
Normalised forcing term  p. 4  
Normalised Rossby deformation length  p. 4  
Normalised friction parameter  p. 4  
σ  Tidal frequency  p. 4 
Rossby number of forced waves  p. 4  
Φ  Divergent displacement potential  p. 4 
Ψ  Rotational displacement streamfunction  p. 4 
n  Outward pointing unit vector defining the normal to the coast  p. 4 
dS  Infinitesimal surface element  p. 5 
𝒪  Domain occupied by the ocean basin 𝒪  p. 5 
ϑ𝒪  Boundary of the ocean basin  p. 5 
Conjugate of a complex number (z)  p. 5  
Normalised Laplacian operator  p. 5  
ALFs  Associated Legendre Functions of the first kind  p. 5 
Φ_{k}  kth eigenfunction of the set of SCHNs  p. 5 
λ_{k}  Eigenvalue associated with Φ_{k}  p. 5 
Ψ_{k}  kth eigenfunction of the set of SCHDs  p. 5 
v_{k}  Eigenvalue associated with Ψ_{k}  p. 5 
i  Imaginary number  p. 5 
SCHs  Spherical cap harmonics  p. 5 
SCHNs  SCHs with Neumann boundary conditions  p. 5 
SCHDs  SCHs with Dirichlet boundary conditions  p. 5 
δ_{k,j}  Kronecker delta function  p. 5 
p_{k}  Weighting coefficient of Φ_{k}  p. 5 
pk  Weighting coefficient of Ψ  p. 5 
SPHs  Spherical harmonics  p. 5 
β_{k,j}  Gyroscopic coefficient  p. 6 
Complex SPH associated with the coordinates of ℛ  p. 7  
l  Degree (or latitudinal wavenumber) of the SPHs  p. 7 
m  Order (or longitudinal wavenumber) of the SPHs  p. 7 
component of the normalised tidal gravitational potential  p. 7  
kth function of the set of SPHs associated with ℛ  p. 7  
kth component of the normalised tidal potential  p. 7  
σcomponent of the normalised tidal gravitational potential  p. 7  
σcomponent of the forcing term  p. 7  
component of  p. 7  
component of  p. 7  
,  Solid deformation factors  p. 7 
Degreel solid tidal gravitational Love number  p. 7  
Degreel solid tidal displacement Love number  p. 7  
Degreel solid loading tidal Love number  p. 7  
Degreel solid loading displacement Love number  p. 7  
Mean density of the solid regions  p. 7  
υ_{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 {φ_{j}}  p. 8 
N  Number of basis functions or truncation degree of the set {ψ_{j}}  p. 8 
K  Number of basis functions or truncation degree of the set {}  p. 8 
Φ  Vector defining φ in terms of the φ_{k}  p. 8 
Ψ  Vector defining ψ in terms of the ψ_{k}  p. 8 
Vector defining 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 
Α_{φ}, R_{φ}  Matrices accounting for the deformation of the solid regions  p. 8 
Transition matrix from the SPHs to the SCHs φ_{k}  p. 8  
Γ_{D}, Γ_{T}  Diagonal matrices of solid deformation factors  p. 8 
Λ, N  Diagonal matrices of eigenvalues  p. 8 
Β  Matrix of gyroscopic coefficients  p. 8 
B_{Φ,Φ}  First block matrix of B  p. 8 
B_{Φ,Ψ}  Second block matrix of B  p. 8 
B_{Ψ,Φ}  Third block matrix of B  p. 8 
Β_{Ψ,Ψ}  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 
Complex (σ, l, m)component of U_{D}  p. 9  
Complex (σ, l, m)component of the normalised surface elevation  p. 9  
Complex Love numbers  p. 9  
𝒯_{z}  Tidal torque exerted about the planet’s spin axis  p. 9 
dV  Infinitesimal volume parcel  p. 9 
δρ  Local tidal density variations  p. 9 
𝒱_{*}  Simply connected region that includes the planet but not the perturber  p. 9 
ϑ𝒱_{*}  Boundary of ϑ𝒱_{*}  p. 9 
r_{*}  Radius of the spherical volume 𝒱_{*}  p. 10 
a  Semimajor axis of the perturber  p. 10 
Ƥ_{T}  Timeaveraged power input by the tidal force  p. 10 
ρ  Local density  p. 10 
Ƥ_{T;oc}  Timeaveraged tidal input power in the ocean  p. 10 
Ƥ_{diss;oc}  Timeaveraged tidally dissipated power in the ocean  p. 10 
Ƥ_{diss}  Total timeaveraged tidally dissipated power  p. 11 
Ƥ_{T;sol}  Timeaveraged tidal input power in the solid regions  p. 11 
Ƥ_{diss;sol}  Timeaveraged tidally dissipated power in the solid regions  p. 11 
^{μ}  Effective shear modulus of the solid regions  p. 11 
τ_{M}  Maxwell time  p. 11 
τ_{A}  Andrade time  p. 11 
α  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}_{⦅}  Lunar mass  p. 12 
^{n}_{⦅}  Lunar mean motion  p. 12 
χ  Normalised semidiurnal tidal frequency  p. 13 
Ω_{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 function  p. 15 
*  Superscript referring to the values of the reference case defined in Table 2  p. 15 
Colatitude of the continental centre in the reference case  p. 15  
Angular radius of the continent in the reference case  p. 15  
H*  Oceanic depth in the reference case  p. 15 
Rayleigh drag frequency in the reference case  p. 15  
^{μ*}  Effective shear modulus in the reference case  p. 15 
Andrade time in the reference case  p. 15  
σ_{k}  Tidal frequency of the kth maximum  p. 17 
Δσ_{k}  kth frequency interval between two consecutive maxima  p. 17 
Mean frequency interval separating two consecutive relative maxima of the tidal torque  p. 18  
kth frequency interval normalised by the mean interval  p. 18  
CDF  Cumulative distribution function  p. 18 
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 , are the solutions of the equation (B.1)
which is found by separating the longitude (φ) and the colatitude (θ) in the eigenvalues problems associated with the Laplace equation in spherical coordinates (see Eq. (C.1)). The function h (θ) 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, nonintegral, or complex depending of 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 ϵ ℂ such that z < 1 and arbitrary complex constants, l and m, as (Abramowitz & Stegun 1972, Eq. (8.1.2)) (B.2)
In the above equation, Γ is the Gamma function introduced in Eq. (106), which is defined for ℜ (z) > 0 as (Abramowitz & Stegun 1972, Eq. (6.1.1)) (B.3)
The hypergeometric function_{2} F_{1} is defined as (Abramowitz & Stegun 1972, Eq. (15.1.1)) (B.4)
where (a)_{j} designates the Pochhammer symbol (Abramowitz & Stegun 1972, Eq. (6.1.22)), (B.5) (B.6)
The derivatives of the ALFs with respect to z are given by (B.7)
We note that the functions Γ and _{2}F_{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 which is defined as (e.g. Johansson 2016) (B.8)
For m ϵ ℤ, l ϵ ℝ and z = x ϵ ℝ such that −1 ≤ x ≤ 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) (B.9)
Here, the hypergeometric function is expressed as (B.10)
with the coefficients a_{j} recursively defined as (B.11) (B.12)
In the general case, the sum of Eq. (B.9) converges only for −1 < x ≤ 1, which corresponds to the inequality 0 ≤ θ < π in the change of coordinates x = cos θ. 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 θ = π in this case.
Since the degrees l tend to integral values (l = 0,1,2,…) as θ → π, the generalised ALFs given by Eqs. (B.2) and (B.9) continuously converge towards the wellknown ALFs of the SPHs, which are defined for −1 ≤ x ≤ 1, and l and m integers such that l ≥ m, as (Abramowitz & Stegun 1972; Arfken & Weber 2005) (B.13)
the P_{l} designating the Legendre polynomials, given by (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 ≈ −1, which occurs if the ocean is almost global. In practice, the evaluation is achieved by making use of the buildin hypergeometric function _{2}F_{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 θ_{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 zaxis 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 eigenvalueeigenfunction problem described by the Helmholtz equation (e.g. Haines 1985), (C.1)
where ∇^{2} designates the Laplace operator, defined for any function f (θ, φ) as (C.2)
the notation ϑ_{x} or ϑ_{xx} referring to the first or second partial derivatives with respect to the x coordinate. In both case, the variables θ and φ can be separated, which allows the problem to be solved as two decoupled eigenvalue problems.
Since the functions f (θ, φ) in Eq. (C.1) are defined for all φ, the boundary conditions on φ are those of continuity, (C.3) (C.4)
which restricts the values of m to be integral (m ϵ ℤ). 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 θ at the pole (θ = 0), (C.5) (C.6)
Fig. C.1 Geometry of the landocean 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 ℛ_{oc}: (𝒪, e_{x}, e_{y}, e_{z}) introduced in Sect. 2.1 is such that 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. 
Fig. C.2 Real degrees l_{k} of the SCHs for doublets (k, m) with 0 ≤ m ≤ k ≤ 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. 
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, 𝒮: (θ, φ) ∈ [0, π] × [0,2π[, with regularity boundary conditions at θ = π similar to those applied at θ = 0, (C.7) (C.8)
They are associated with integral degrees l ∈ ℕ such that l ≥ 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)) (C.9)
where the 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 (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 (C.12)
the notation dS referring to an infinitesimal surface element of the unit sphere, and to the conjugate of f. This scalar product allows the ℓ^{2}norm of the space of functions of 𝒮 to be defined, (C.13)
Hence, for any pair of SPHs and , (C.14)
The SCHs are obtained by solving the Helmholtz equation given by Eq. (C.1) on the domain filled by the spherical cap, 𝒞 : (θ, φ) e [0, θ_{0}] x [0,2π[, together with the boundary conditions specified at ∂𝒞, namely θ = θ_{0}. In the present model, the applied boundary conditions are either Neumann (fixed value for ∂_{θ}f at θ = θ_{0}) or Dirichlet (fixed value for the function f itself) conditions, which are respectively formulated as (C.15) (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).
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 (θ_{c}).
As the Helmholtz equation is the same for both the SCHs and the SPHs, solving the eigenfunctioneigenvalue 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 ∈ ℤ. As a consequence, the eigenfunctioneigenvalue problem becomes a rootfinding problem where one computes the series of real degrees l such that, alternately, (C.17) (C.18)
This can be achieved numerically using a standard NewtonRaphson 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 ∈ ℕ such that n ≥ m. These indices^{3} correspond to the integral degrees of the SPHs towards which the SCHs converge for the ℓ^{2}norm defined by Eq. (C.13) in the global ocean limit: l_{n} → n as θ_{0} → 180° 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 ≤ m ≤ 2 and various values of angular radius of the continent, θ_{c:} = 180 − θ_{0} (in degrees). They are plotted on Fig. C.2.
Eigenvalues of the SCHs with Neumann (SCHNs) or Dirichlet (SCHDs) boundary conditions for various angular radii of the supercontinent.
We note that l_{n} ≥ 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 θ_{0} ≤ 90° since there exists degrees l_{m} < m as long as θ_{0} > 90°. Asymptotic scaling laws may be used to estimate the values of the degrees l_{n} ≫ m as functions of the indices n and the angle θ_{0} (see e.g. Thébault et al. 2006, Eq. (40)). Particularly, these laws show that for θ_{0} ≪ 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,…, +∞ and l = l_{m}, l_{m+1},…, as (C.19)
where the normalisation coefficients α_{l,m} ∈ ℝ are to be specified. The corresponding eigenvalues λ_{l} are expressed as (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 ≤ m ≤ n ≤ 2 and various values of the continental angular radius (θ_{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 θ > θ_{0}, this scalar product reads (C.21)
For any pair of SCHs belonging to the same set, and , the property given by Eq. (C.14) for the SPHs, (C.22)
is enforced by defining the normalisation coefficients α_{l,m} ∈ ℝ as (C.23)
However, the functions in one set are not orthogonal to those in the other since (e.g. Haines 1985) (C.24)
for and taken in the sets of SCHNs and SCHDs, respectively. In the general case, and , which implies 〈, 〉 ≠ 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 𝒞.
In the case of real eigenfunctions defined as (C.25)
with φ_{m} = 0 if m ≥ 0 and φ_{m} = π/2 if m < 0, the orthogonality property given by Eq. (C.22) still holds and the normalisation coefficients defined in Eq. (C.23) become (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 θ_{0} increases, the ALFs smoothly evolve from the Bessel functions (θ_{0} → 0) to the standard ALFs of integral parameters (θ_{0} → π). Similarly, the SCHs evolve from the cylindrical harmonics, where there is no curvature effect (smallsea configuration), to the SPHs given by Eq. (C.9) (globalocean 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 sets of eigenfunctions. Besides, the limit θ_{0} → π is not strictly equivalent to the globalocean configuration since boundary conditions are still applied at ∂𝓒, 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, , both of them tending to zero as the size of the island decreases.
Appendix D Euler rotation matrix
We consider the frame of reference ℛ: (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 ℛ: (O, e_{X}, e_{Y}, e_{Z}) → ℛ_{rot}: (θ, 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 Zaxis through an angle α such that 0 ≤ α < 2π, (ii) a rotation about the new Y′axis through an angle β such that 0 ≤ β ≤ π, and (iii) a rotation about the new axis z through an angle γ such that 0 ≤ γ ≤ 2π. The frame of reference ℛ_{rot}: (o, e_{x}, e_{y}, e_{z}) is associated with the Cartesian coordinates (x, y, y). The three rotation angles, α, β and γ, are the socalled Euler angles.
Introducing the notations C_{ϕ} = cos _{ϕ} and S_{ϕ} = sin _{ϕ} for ϕ = α,β or γ, the standard Euler rotation matrix is defined in terms of the Euler angles as (e.g. Varshalovich et al. 1988, Eq. (54) p. 30) (D.1)
which is formed by the product of three rotation matrices, (D.2)
The inverse rotation matrix R^{−1} is given by (D.3)
The Euler rotation matrix and its inverse correspond to the change of basis matrices between the initial and rotated vectorial bases, (e_{X}, e_{Y}, e_{Z}) → (e_{x}, e_{y}, e_{z}), which is expressed as (D.4)
Besides, the rotation U → U_{rot} of any vector U through the Euler angles α,β and γ in the Cartesian coordinate system S {X, Y, Z} is expressed as U_{rot} = R_{α,β,γ}U.
In spherical coordinates, the transformation relations are deduced from Eq. (D .4) by making use of trigonometric identities. We denote by (, , ) the spherical coordinates associated with the initial frame of reference, ℛ, where designates the radial coordinate, the colatitude, and the longitude. The analogous coordinates in the rotated frame of reference, ℛ_{rot}, are denoted by (r, θ, φ), with r = . In the forward transformation, (, , ) → (r, θ, φ), the relations between angles (, ) and (θ, φ) are given by (Varshalovich et al. 1988, Sect. 1.4.1, Eqs (23)) (D.5) (D.6)
while the inverse relations are (D.7) (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 (D.9)
In addition, by considering the same vector as the image of e_{Z} through a rotation of Euler angles α, β, and γ, we can write e_{z} = R_{α,β,γ} (0,0,1)^{T}, where we made use of the standard Euler rotation matrix defined in Eq. (D.1). It follows (D.10)
By comparing Eqs. (D.9) and (D.10), we identify and . Considering the inverse transformation, ,weget (D.11)
which shows that e_{Z} does not depend on the rotation through the angle α in the oceanic frame of reference.
We note that the angle γ describes a rotation of coordinates around the axis defined by e_{Z}. Since the ocean is circular, the geometry remains unchanged through this rotation. The angle γ can therefore be set to γ = 0, which allows the position of the North pole in the ocean’s coordinate system to be simplified to e_{Z} = − sinβe_{x} + cosβe_{z}. The coordinates of the North pole in the coordinate system defined from the centre of the ocean basin are thus given by (θ_{NP}, φ_{ΝΡ}) = (β,π). At the end, the ocean’s coordinate system is simply obtained from the usual geocentric coordinate system through a rotation of Euler angles (α,β,γ) = (, , 0). We shall emphasise that the gyroscopic coefficients of the tidal theory only depend on the colatitude of the supercontinent’s centre since ß 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 α can be ignored owing to the periodicity of the tidal response across the longitudinal coordinate, , which implies that the longitude of the supercontinent’s centre can be set to without any loss of generality.
Appendix E Rotation of spherical harmonics
We consider a rotation ℛ: (O, e_{X}, e_{Y}, e_{Z}) → ℛ_{rot}: (O, e_{x}, e_{y}, e_{z}) described by the Euler angles introduced in Appendix D, which are denoted by α, β, and γ (see Fig. D.1). The initial frame of reference, ℛ, is associated with the spherical coordinate system (, , ), while the rotated frame of reference, ℛ_{rot}, is associated with the coordinate system (r, θ, φ). 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), (E.1)
and, conversely, the SPHs of the initial coordinate system are expressed as (E.2)
The notation refers to the Wigner Dfunctions (Varshalovich et al. 1988, Sect. 4.3, Eq. (1)), (E.3)
where is a real function expressed as (Varshalovich et al. 1988, Sect. 4.3.1, Eq. (2)) (E.4)
In the above equation, k runs over all integer values for which the factorial arguments are nonnegative, namely max (0, −q − m) ≤ k ≤ min (l − q, l − m).
Fig. E.1 Rotation matrix D of Eq. (E.5) for Euler angles (α,β,γ) = (120°, 60°, 270°). 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 ≤ m ≤ l) grouped together. 
The Wigner Dfunctions introduced in Eqs. (E.1) and (E.2) are the matrix elements of the rotation operator, (E.5)
where L is the chosen truncation degree, and is the matrix defined as (E.6)
Thus, the set of SPHs associated with the rotated coordinate system, denoted by Y ≡ [Υ_{1},…, Y_{K}]^{T}, is expressed as a function of the initial set, , through the algebraic relation (E.7)
which is analogous to the change of basis equation given by Eq. (D.4). The inverse transformation is given by (E.8)
We remark that , 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 (α,β,γ) = (120°, 60°, 270°) with a truncation degree set to L = 4. As previously observed in Eq. (E.1), the rotated SPHs are expressed in terms of the nonrotated 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, and , associated with the same coordinate system, ℛ, but not necessarily the same definition domains^{4}, denoted by and , respectively. Both and 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 and , where 𝓐 and g are assumed to have the same definition domains as and , respectively. The expressions of the standard Euler rotation matrices R_{𝓐} and R_{ℬ} describing these transformations are given by Eq. (D.1), and the corresponding Wigner rotation matrices for the functional bases, D_{𝓐} and D_{ℬ}, are given by Eq. (E.5). We note that the second Euler angle can only be set to β = 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 α and γ since these angles just introduce a phase lag in the component of the SCHs.
The initial and rotated functional basis vectors are denoted by , , Y_{𝓐}, and Y_{ℬ}, the vector Y_{𝓗} of any set 𝓗 being expressed as (F.1)
where the Y_{𝓗,k} designate the orthogonal functions of set 𝓗, and K_{𝓗} its total number of elements, which depends on the truncation degree LH introduced in Eq. (E.5). To lighten expressions, we assume that the setsof basis functions and 𝒜 contain M elements, and that and ℬ contain N elements, so that and . Using the above notations, we can write down the relations (F.2)
Fig. F.1 Transition matrices between the SCHNs, SCHDs and SPHs for an ocean basin of angular radius θ_{0} = 130°. 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,…, 3,4), with degrees (l ≥ m) grouped together. 
As a next step, we establish the relationship between the bases and , which share the same coordinate system. The functions and are expressed in terms of each other through the relations (F.3) (F.4)
In the above equations, the functions and are the residuals of the functions or when the latter are expanded in terms of the functions of set or , respectively. These residuals result both from the truncation of the sets and  which are mathematically defined as infinitedimension sets of basis functions for the definition domains and –, and from the fact that and are different from each other in the general case.
For instance, if , the functions of set g can be expressed as linear combinations of the elements of set on only. This is the case for SPHs, which cannot be approximated using SCHs while the opposite is true. Therefore, the function tends to be fully described by the functions of set as M → ∞ if , and consequently in that case. Conversely, for M = ∞ in the general case if . The residual function behaves in a similar way if the domains and are interchanged: as N → ∞ if , while if .
Since the functions of the two sets are assumed to be orthonormal by construction (see Appendix C), for all k ∈ {1,…, M} and j ∈ {1,…, N}. As a consequence, the expressions given by Eqs. (F.3) and (F.4) simplify to (F.5) (F.6)
which, introducing the vectors (F.7) (F.8)
Fig. F.2 Real part of complex SCH of order m = 1 and degree l_{3} = 6.806057 corresponding to an ocean basin of angular radius θ_{0} = 80° (i.e. = 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). 
In the above equations, the transition matrix is expressed as (F.10)
where , is the scalar product of the kth element of with the jth element of 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 and , namely .
With the used normalisations, the transition matrix satisfies the properties , , and , the notation I_{M} referring to the identity matrix of size M. However, one should pay attention to the fact that and 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 between the SCHNs, SCHDs, and SPHs for an ocean of angular radius θ_{0} = 130° (i.e. a continent of angular radius θ_{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 (F.11)
where the transition matrix P_{𝒜,ℬ} is given by (F.12)
A given function f defined on the sphere can thus be written both in terms of the basis functions of 𝒜 and those of ℬ as (F.13)
where we have introduced the residual functions ϵ_{𝒜,f} and ϵ_{ℬ,f}, and the coordinate vectors of f in the two sets of basis functions, y_{𝒜} and y_{ℬ}, defined as (F.14) (F.15)
The vectors y_{𝒜} and y_{ℬ} contain the coefficients of the expansion of fon the set of basis functions 𝒜 or ℬ, respectively. They are related to each other by the transition matrix P_{𝒜,ℬ} introduced in Eq. (F.12) through the relation equations (F.16)
Figure F.2 shows the three steps of the transition from SCHs to rotated SPHs through the example of the eigenfunction 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 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 θ = θ_{0} gives rise to ‘ringing’ when is expanded in series of SPHs. This effect also alters the expansions of the SCHNs due to their discontinuous transition at θ = θ_{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 β They are defined as (G.1) (G.2) (G.3) (G.4)
where the Φ_{k} and Ψ_{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 . 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 nonintegral 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, , and using trigonometric identities, we express cos as a function of the colatitude of the ocean’s centre () and the ocean’s coordinate system, (G.5) (G.6)
Considering the definitions of the SPHs and SCHs given by Eqs. (C.9) and (C.19), the functions Φ_{k} and Ψ_{k} are expressed in the general case as (G.7) (G.8)
where (cos θ) designates the ALF of the kth function of the set of SCHNs (Neumann conditions) used to expand the potential function, and (cos θ) the ALF of the kth 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 (G.9) (G.10) (G.11) (G.12)
Fig. G.1 Gyroscopic coefficients of an ocean basin of angular radius θ_{0} = 130° (i.e. θ_{c} = 50°) for various values of the colatitude of the ocean’s centre. Top: (polar ocean). Middle: (midlatitude ocean). Bottom: (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,…, 3, 4), with degrees (l ≥ 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. 
Interestingly, the coefficients β_{±k,± 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 (), the gyroscopic coefficients only mix up eigenmodes with the same order (m_{j} = m_{k}). In this configuration, . If the ocean is equatorial , 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 associated with eigenmodes having the same order, and one component in sin associated with eigenmodes having neighbouring orders. The gyroscopic coefficients that do not mix up potential and stream functions (β_{k,−j} and β_{−k,−j}) are pure imaginary numbers. Conversely, the coefficients that mix up the functions (β_{k,−j} and β_{−k,j}) are pure real numbers.
By comparing the expressions given by Eqs. (G.9) and (G.12), we remark that the β_{k,j} and β_{−k,−j} converge towards the same asymptotic values as θ_{0} → 180° since both the SCHNs and SCHDs converge towards the standard SPHs for the ℓ^{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 θ_{0} = 130° for three typical configuration: (polar ocean), (midlatitude ocean), and (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 υ_{k,q} introduced in Eq. (51) are evaluated similarly as the gyroscopic coefficients. These coefficients are defined as (G.13)
where the SPHs Y_{q} are written as products of functions of separated variables, (G.14)
After an integration by part, the coefficients are simply expressed as (G.15)
We note that the υ_{k,q} are zero if m_{k} = 0. Besides, if m_{k}  > 0, the Dirichlet condition satisfied by the Ψ_{k} at θ_{0} implies that G_{k} (θ_{0}) = 0, and the polar condition of Y_{q} implies that H_{q} (0) = 0. As a consequence, υ_{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 elastogravitational 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 degreel Love numbers of a homogeneous solid body of mass M_{core}, radius R_{core}, and uniform shear modulus μ, are expressed as (e.g. Munk & MacDonald 1960) (H.1)
In the above equation, the degreel dimensionless shearmodulus is defined as (H.2)
where designates the dimensionless viscoelastic rigidity accounting for the rheology of the material, and A_{l} the dimensionless parameter given by (CastilloRogez et al. 2011; Efroimsky 2012) (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} ≈ M_{p} and R_{core} ≈ R_{p}
The normalised rigidity 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 viscoelastic rheology (e.g. Correia et al. 2014), it is expressed as a simple function of the tidal frequency, σ. Following Bolmont et al. (2020), we opt for an Andrade rheology. Similarly as the Maxwell rheology, the Andrade model describes the viscoelastic deformation of the material. However, it additionally accounts for the unrecoverable creep forming the anelastic component of the response, which affects the frequencydependence of tidal dissipation in the highfrequency range (e.g. Efroimsky 2012). In the Andrade rheology, is expressed as (e.g. CastilloRogez et al. 2011) (H.4)
where τ_{Μ} designate the Maxwell time, τ_{Α} the Andrade time, Γ the Gamma function introduced in Eq. (B.3), and α the parameter determining the duration of the transient response in the primary creep (Castelnau et al. 2008; CastilloRogez et al. 2011).
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions, National Bureau of Standards: Applied Mathematics Series, 55 [Google Scholar]
 Arbic, B. K., & Garrett, C. 2010, Continent. Shelf Res., 30, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Arbic, B. K., Wallcraft, A. J., & Metzger, E. J. 2010, Ocean Model., 32, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Arfken, G. B., & Weber, H. J. 2005, Mathematical Methods for Physicists, 6th edn. (Elsevier Academic Press) [Google Scholar]
 AuclairDesrotour, P., Le PoncinLafitte, C., & Mathis, S. 2014, A & A, 561, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Mathis, S., & Le PoncinLafitte, C. 2015, A & A, 581, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Mathis, S., Laskar, J., & Leconte, J. 2018, A & A, 615, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Leconte, J., Bolmont, E., & Mathis, S. 2019a, A & A, 629, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Leconte, J., & Mergny, C. 2019b, A & A, 624, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuthe, M. 2016, Icarus, 280, 278 [CrossRef] [Google Scholar]
 Bills, B. G., & Ray, R. D. 1999, Geophys. Res. Lett., 26, 3045 [NASA ADS] [CrossRef] [Google Scholar]
 Biot, M. A. 1954, J. Appl. Phys., 25, 1385 [Google Scholar]
 Blackledge, B. W., Green, J. A. M., Barnes, R., & Way, M. J. 2020, Geophys. Res. Lett., 47, e85746 [NASA ADS] [CrossRef] [Google Scholar]
 Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 [NASA ADS] [CrossRef] [Google Scholar]
 Bolmont, E., Breton, S. N., Tobie, G., et al. 2020, A & A, 644, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Cadieux, C., Doyon, R., Plotnykov, M., et al. 2022, AJ, 164, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Carr, M. H., & Head, J. W. 2003, J. Geophys. Res. (Planets), 108, 5042 [NASA ADS] [Google Scholar]
 Carter, G. S., Merrifield, M., Becker, J. M., et al. 2008, J. Phys. Oceanogr., 38, 2205 [NASA ADS] [CrossRef] [Google Scholar]
 Castelnau, O., Duval, P., Montagnat, M., & Brenner, R. 2008, J. Geophys. Res. (Solid Earth), 113, B11203 [NASA ADS] [CrossRef] [Google Scholar]
 CastilloRogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. (Planets), 116, E09008 [CrossRef] [Google Scholar]
 Charbonneau, D., Berta, Z. K., Irwin, J., et al. 2009, Nature, 462, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodriguez, A. 2014, A & A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cowling, T. G. 1941, MNRAS, 101, 367 [NASA ADS] [Google Scholar]
 Daher, H., Arbic, B. K., Williams, J. G., et al. 2021, J. Geophys. Res. (Planets), 126, e06875 [NASA ADS] [Google Scholar]
 Dohm, J. M., Baker, V. R., Boynton, W. V., et al. 2009, Planet. Space Sci., 57, 664 [NASA ADS] [CrossRef] [Google Scholar]
 Eakins, B. W., & Sharman, G. F. 2010, NOAA National Geophysical Data Center, Boulder, CO [Google Scholar]
 Efroimsky, M. 2012, ApJ, 746, 150 [Google Scholar]
 Egbert, G. D., & Ray, R. D. 2001, J. Geophys. Res., 106, 22 [Google Scholar]
 Egbert, G. D., & Ray, R. D. 2003, Geophys. Res. Lett., 30, 1907 [NASA ADS] [CrossRef] [Google Scholar]
 Farhat, M., AuclairDesrotour, P., Boué, G., & Laskar, J. 2022a, A & A, 665, A1 [Google Scholar]
 Farhat, M., Laskar, J., & Boué, G. 2022b, J. Geophys. Res. (Solid Earth), 127, e2021JB023323 [NASA ADS] [CrossRef] [Google Scholar]
 FoxKemper, B., Ferrari, R., & Pedlosky, J. 2003, J. Phys. Oceanogr., 33, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Garrett, C. J. R., & Munk, W. H. 1971, Deep Sea Res. A, 18, 493 [NASA ADS] [Google Scholar]
 Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194 [Google Scholar]
 Gent, P. R., & McWilliams, J. C. 1983, Dyn. Atmos. Oceans, 7, 67 [CrossRef] [Google Scholar]
 Gerkema, T., & Zimmerman, J. 2008, Lecture Notes, Royal NIOZ, Texel [Google Scholar]
 Gordon, C., Webb, D., & Wolpert, S. 1992, Invent. Math., 110, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Green, J. A. M., & Huber, M. 2013, Geophys. Res. Lett., 40, 2707 [NASA ADS] [CrossRef] [Google Scholar]
 Green, J. A. M., Huber, M., Waltham, D., Buzan, J., & Wells, M. 2017, Earth Planet. Sci. Lett., 461, 46 [CrossRef] [Google Scholar]
 Grimm, S. L., Demory, B.O., Gillon, M., et al. 2018, A & A, 613, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haines, G. V. 1985, J. Geophys. Res., 90, 2583 [NASA ADS] [CrossRef] [Google Scholar]
 Han, L., & Huang, R. X. 2020, J. Phys. Oceanogr., 50, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Hendershott, M. C. 1972, Geophys. J., 29, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1981, A & A, 99, 126 [NASA ADS] [Google Scholar]
 Hwang, C., & Chen, S.K. 1997, Geophys. J. Int., 129, 450 [NASA ADS] [CrossRef] [Google Scholar]
 Johansson, F. 2016, ArXiv eprints [arXiv: 1606.06977] [Google Scholar]
 Kac, M. 1966, Am. Math. Monthly, 73, 1 [CrossRef] [Google Scholar]
 Kamata, S., Matsuyama, I., & Nimmo, F. 2015, J. Geophys. Res. (Planets), 120, 1528 [NASA ADS] [CrossRef] [Google Scholar]
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [Google Scholar]
 Kivelson, M. G., Khurana, K. K., Russell, C. T., et al. 2000, Science, 289, 1340 [CrossRef] [Google Scholar]
 Kodaira, T., Thompson, K. R., & Bernier, N. B. 2016, J. Geophys. Res. (Oceans), 121, 6159 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. K., Wolf, E. T., Arney, G., et al. 2017, ApJ, 845, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Lambeck, K. 1977, Philos. Trans. Roy. Soc. Lond. Ser. A, 287, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Laplace, P. S. 1798, Traité de mécanique céleste (Duprat J. B. M.), Tome 5, Livre XIII, 145 [Google Scholar]
 Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [Google Scholar]
 Léger, A., Selsis, F., Sotin, C., et al. 2004, Icarus, 169, 499 [Google Scholar]
 Le Provost, C., Genco, M. L., Lyard, F., Vincent, P., & Canceil, P. 1994, J. Geophys. Res., 99, 24, 777 [NASA ADS] [Google Scholar]
 Le Provost, C., Lyard, F., Molines, J. M., Genco, M. L., & Rabilloud, F. 1998, J. Geophys. Res., 103, 5513 [NASA ADS] [CrossRef] [Google Scholar]
 Levrard, B., Winisdoerffer, C., & Chabrier, G. 2009, ApJ, 692, L9 [NASA ADS] [CrossRef] [Google Scholar]
 LilloBox, J., Figueira, P., Leleu, A., et al. 2020, A & A, 642, A121 [CrossRef] [EDP Sciences] [Google Scholar]
 Lindzen, R. S., & Chapman, S. 1969, Space Sci. Rev., 10, 3 [NASA ADS] [CrossRef] [Google Scholar]
 LonguetHiggins, M. S. 1968, Philos. Trans. Roy. Soc. Lond. Ser. A, 262, 511 [NASA ADS] [CrossRef] [Google Scholar]
 LonguetHiggins, M. S., & Pond, G. S. 1970, Philos. Trans. Roy. Soc. Lond. Ser. A, 266, 193 [NASA ADS] [CrossRef] [Google Scholar]
 MacDonald, G. J. F. 1964, Rev. Geophys. Space Phys., 2, 467 [Google Scholar]
 Makarov, V. V., & Efroimsky, M. 2013, ApJ, 764, 27 [Google Scholar]
 Mamajek, E. E., Prsa, A., Torres, G., et al. 2015, ArXiv eprints [arXiv:1510.07674] [Google Scholar]
 Matsuyama, I. 2014, Icarus, 242, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuyama, I., Beuthe, M., Hay, H. C. F. C., Nimmo, F., & Kamata, S. 2018, Icarus, 312, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Merdith, A. S., Williams, S. E., Collins, A. S., et al. 2021, Earth Sci. Rev., 214, 103477 [NASA ADS] [CrossRef] [Google Scholar]
 Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics (McGrawHill Book Comp.) [Google Scholar]
 Motoyama, M., Tsunakawa, H., & Takahashi, F. 2020, Icarus, 335, 113382 [NASA ADS] [CrossRef] [Google Scholar]
 Munk, W. H., & MacDonald, G. J. F. 1960, J. Geophys. Res., 65, 2169 [NASA ADS] [CrossRef] [Google Scholar]
 Ogilvie, G. I. 2013, MNRAS, 429, 613 [Google Scholar]
 Ogilvie, G. I. 2014, ARA & A, 52, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Payne, M. J., & Lodato, G. 2007, MNRAS, 381, 1597 [Google Scholar]
 Peale, S. J., Cassen, P., & Reynolds, R. T. 1979, Science, 203, 892 [Google Scholar]
 Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 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) [Google Scholar]
 Proudman, J. 1920, Proc. Lond. Math. Soc., 2, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Provost, C. L., & Lyard, F. 1997, Progr. Oceanogr., 40, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, ApJ, 669, 606 [Google Scholar]
 Remus, F., Mathis, S., Zahn, J.P., & Lainey, V. 2012, A & A, 541, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scheller, E. L., Ehlmann, B. L., Hu, R., Adams, D. J., & Yung, Y. L. 2021, Science, 372, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Schwiderski, E. W. 1980, Rev. Geophys. Space Phys., 18, 243 [CrossRef] [Google Scholar]
 Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337 [Google Scholar]
 Strauss, W. A. 2007, Partial Differential Equations: An Introduction (John Wiley & Sons) [Google Scholar]
 Takeuchi, H., & Saito, M. 1972, Methods Computat. Phys., 11, 217 [Google Scholar]
 Thébault, E., Schott, J. J., Mandea, M., & Hoffbeck, J. P. 2004, Geophys. J. Int., 159, 83 [CrossRef] [Google Scholar]
 Thébault, E., Schott, J. J., & Mandea, M. 2006, J. Geophys. Res. (Solid Earth), 111, B01102 [Google Scholar]
 Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534 [Google Scholar]
 Torta, J. M. 2019, Surv. Geophys., 41, 201 [Google Scholar]
 Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A & A, 612, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tyler, R. H. 2008, Nature, 456, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Tyler, R. 2011, Icarus, 211, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Tyler, R. 2014, Icarus, 243, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Tyler, R. H. 2021, Planet. Sci. J., 2, 70 [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: University of Tokyo Press) [Google Scholar]
 Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics, 770 [Google Scholar]
 Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1988, Quantum Theory of Angular Momentum (World Scientific Publishing Company) [CrossRef] [Google Scholar]
 Watterson, I. G. 2001, J. Atmos. Ocean. Technol., 18, 691 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. J. 1973, Deep Sea Res. A, 20, 847 [NASA ADS] [Google Scholar]
 Webb, D. J. 1980, Geophys. J., 61, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. J. 1982, Geophys. J., 70, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Wunsch, C., Haidvogel, D. B., Iskandarani, M., & Hughes, R. 1997, Progr. Oceanogr., 40, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J. P. 1966, Ann. Astrophys., 29, 313 [NASA ADS] [Google Scholar]
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).
Other values of σ_{R} may be found in the literature. For example, Wunsch et al. (1997) uses a much smaller value to study the dynamics of the longperiod tides, σ_{R} = 2.5 × 10^{−7} s^{−1}, while Motoyama et al. (2020) set this value to σ_{R} ~ 4 × 10^{−6} s^{−1} following the prescription given by Schwiderski (1980).
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.
By definition domain, we refer to the domain where the functions are nonzero, the SCHs being defined on the full unit sphere similarly as the SPHs in the used convention (see Eq. (C.19)).
All Tables
Dimensionless control parameters determining the regime of the oceanic tidal response.
Evolution of the frequency spectrum of the quadrupolar tidal Love number as one parameter increases while the values of other parameters are fixed.
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 (θ_{c}).
Eigenvalues of the SCHs with Neumann (SCHNs) or Dirichlet (SCHDs) boundary conditions for various angular radii of the supercontinent.
All Figures
Fig. 1 Geometry of the studied system and associated parameters. The unit vector , which indicates the position of the continental centre, is defined as , with e_{z} pointing towards the centre of the ocean basin. In the diagram, the longitude of the continental centre is set to . 

In the text 
Fig. 2 Eigenfunctions describing the oceanic tidal response for a supercontinent of angular radius θ_{c} = 10° (i.e. θ_{0} = 170°), and associated tidal flows. Top: set {Φ_{k}} (SCHNs). Bottom: set {Ψ_{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,…, 3 (vertical axis) and −n ≤ m ≤ n (horizontal axis). Bright or dark colours designate positive or negative values of the eigenfunctions, respectively. Streamlines indicate the tidal flows corresponding to for the set {Φ_{k}} and to for the set {Ψ_{k}}. 

In the text 
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 χ = (Ω − n_{s}) /Ω_{Ε} (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 ( and 𝒯_{z;2} ≤ 0) and positive tidal input power (𝒫_{T} ≥ 0), while dashed lines indicate accelerating tidal torque and 𝒯_{z;2} > 0) and negative tidal input power (𝒫_{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. 

In the text 
Fig. 4 Frequency spectra of the imaginary part of the quadrupolar Love number associated with the semidiurnal tide, , 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 χ = (Ω − n_{s}) /Ω_{Ε} (horizontal axis), with Ω_{Ε} being the spin angular velocity of the actual Earth. The solid black line designates the reference case defined by Table 2: , , H* = 4.0 km, , μ* = 25.1189 GPa, and . 

In the text 
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: , , H* = 4.0 km, , μ* = 25.1189 GPa, and . 

In the text 
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: , , H^{*} = 4.0 km, , μ* = 25.1189 GPa, and . 

In the text 
Fig. C.1 Geometry of the landocean 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 ℛ_{oc}: (𝒪, e_{x}, e_{y}, e_{z}) introduced in Sect. 2.1 is such that 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. 

In the text 
Fig. C.2 Real degrees l_{k} of the SCHs for doublets (k, m) with 0 ≤ m ≤ k ≤ 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. 

In the text 
Fig. D.1 Euler angles corresponding to the standard Euler rotation matrix defined by Eq. (D.1). 

In the text 
Fig. E.1 Rotation matrix D of Eq. (E.5) for Euler angles (α,β,γ) = (120°, 60°, 270°). 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 ≤ m ≤ l) grouped together. 

In the text 
Fig. F.1 Transition matrices between the SCHNs, SCHDs and SPHs for an ocean basin of angular radius θ_{0} = 130°. 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,…, 3,4), with degrees (l ≥ m) grouped together. 

In the text 
Fig. F.2 Real part of complex SCH of order m = 1 and degree l_{3} = 6.806057 corresponding to an ocean basin of angular radius θ_{0} = 80° (i.e. = 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). 

In the text 
Fig. G.1 Gyroscopic coefficients of an ocean basin of angular radius θ_{0} = 130° (i.e. θ_{c} = 50°) for various values of the colatitude of the ocean’s centre. Top: (polar ocean). Middle: (midlatitude ocean). Bottom: (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,…, 3, 4), with degrees (l ≥ 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. 

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