Issue 
A&A
Volume 615, July 2018



Article Number  A23  
Number of page(s)  15  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201732249  
Published online  06 July 2018 
Oceanic tides from Earthlike to ocean planets
^{1}
IMCCE, Observatoire de Paris, CNRS UMR 8028, PSL,
77 Avenue DenfertRochereau,
75014 Paris, France
email: jacques.laskar@obspm.fr
^{2}
Laboratoire AIM ParisSaclay, CEA/DRF  CNRS  Université Paris Diderot, IRFU/DAp Centre de Saclay,
91191 GifsurYvette, France
^{3}
LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité,
5 place Jules Janssen,
92195 Meudon, France
email: stephane.mathis@cea.fr
^{4}
Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy SaintHilaire,
33615 Pessac, France
email: pierre.auclairdesrotour@ubordeaux.fr
Received:
7
November
2017
Accepted:
24
January
2018
Context. Oceanic tides are a major source of tidal dissipation. They drive the evolution of planetary systems and the rotational dynamics of planets. However, twodimensional (2D) models commonly used for the Earth cannot be applied to extrasolar telluric planets hosting potentially deep oceans because they ignore the threedimensional (3D) effects related to the ocean’s vertical structure.
Aims. Our goal is to investigate, in a consistant way, the importance of the contribution of internal gravity waves in the oceanic tidal response and to propose a modelling that allows one to treat a wide range of cases from shallow to deep oceans.
Methods. A 3D ab initio model is developed to study the dynamics of a global planetary ocean. This model takes into account compressibility, stratification, and sphericity terms, which are usually ignored in 2D approaches. An analytic solution is computed and used to study the dependence of the tidal response on the tidal frequency and on the ocean depth and stratification.
Results. In the 2D asymptotic limit, we recover the frequencyresonant behaviour due to surface inertialgravity waves identified by early studies. As the ocean depth and Brunt–Väisälä frequency increase, the contribution of internal gravity waves grows in importance and the tidal response becomes 3D. In the case of deep oceans, the stable stratification induces resonances that can increase the tidal dissipation rate by several orders of magnitude. It is thus able to significantly affect the evolution time scale of the planetary rotation.
Key words: hydrodynamics / planetstar interactions / planets and satellites: oceans / planets and satellites: terrestrial planets
© ESO 2018
1 Introduction
Even though the number of detected extrasolar planets located in the habitable zone of their host stars has kept growing continuously for the past two decades, two major discoveries recently aroused excitement within the community of exoplanet research. Firstly, in Summer, 2016, a telluric extrasolar planet was found orbiting the Sun’s closest stellar neighbour, the red dwarf Proxima Centauri (AngladaEscudé et al. 2016; Ribas et al. 2016). This planet, Proxima b, has a minimum mass of 1.3 M_{⊕} and an equilibrium temperature that is within the range where water could be liquid on its surface, which makes it resemble a potential twin sister of the Earth. A few months later, another star was in the spotlight, the ultracool dwarf star TRAPPIST1 (Gillon et al. 2017). The system hosted by TRAPPIST1 drew attention due to its remarkable architecture, composed of eight telluric planets of masses between 0.7 M_{⊕} and 1.2 M_{⊕} and radii between 0.1 R_{⊕} and 1.5 R_{⊕} (Gillon et al. 2017; Wang et al. 2017). Most of these planets (i.e. b, c, d, e, f, g, h) exhibit small densities suggesting that water stands for a large fraction of their masses (Wang et al. 2017). Hence, Proxima b and several of the telluric planets orbiting TRAPPIST1 could host deep oceans of liquid water.
Oceans play a crucial role in the evolution of planetary systems. Similarly to rocky cores, they are distorted by the tidal gravitational forcings of perturbers such as stars or satellites. The energy dissipated by the resulting oceanic tides can be of the same order of magnitude and even greater than that associated with the solid tide, as is seen on Earth. Observed as an exoplanet, the Earth appears as a very dry rocky planet with an oceanic layer of depth approximatively equal to 6 × 10^{−4} R_{⊕} (Eakins & Sharman 2010). Yet, the Earth’s ocean is today the main contributor to the total energy tidally dissipated by the EarthMoonSun system. It accounts for roughly 95% of the planetary dissipation rate of 2.54 TW generated by the Lunar semidiurnal tide (see Lambeck 1977; Ray et al. 2001; Egbert & Ray 2001, and references therein). Moreover, the oceanic dissipation rate can vary significantly owing to the frequencyresonant behaviour of the fluid layer (e.g. Webb 1980), thus leading to the presentday high transfer of angular momentum from the spin of the Earth to the orbit of the Moon (Lambeck 1980; Bills & Ray 1999). More generally, oceanic tides have an impact on the history of the spin rotation and states of equilibrium of the Earth (Neron de Surgy & Laskar 1997) and terrestrial exoplanets (Correia et al. 2008; Leconte et al. 2015; AuclairDesrotour et al. 2017b). Therefore, the oceanic tidal dissipation rate needs to be quantified in the case of extrasolar planets to study the history of observed planetary systems and constrain their physical properties.
Since Laplace’s pioneering works (Laplace 1798), the Earth’s oceanic tides have been examined through different approaches. The dissipation rate they induce is now well constrained thanks to the fit of the altimetric data provided by the TOPEX/Poseidon satellite with numerical simulations made using general circulation models (GCM; Egbert & Ray 2001, 2003). Moreover, numerous studies have characterized their dependence on the planet parameters by developing linear ab initio global models, such as the hemispherical ocean model (Proudman & Doodson 1936; Doodson 1938; LonguetHiggins & Pond 1970; Webb 1980). This approach has recently been used to estimate the tidal dissipation rate in the oceans of icy satellites orbiting Jupiter and Saturn (Tyler 2011, 2014; Chen et al. 2014; Matsuyama 2014).
Simplified ab initio models are very convenient to explore the domain of parameters because they involve a small number of control parameters and require a relatively small computational cost compared to numerical models taking into account local features (e.g. topography and mean flows). They are commonly based on a two dimensional (2D) spherical geometry and the thin shell hypothesis. In this approach, the tidal response of the ocean is controlled by the Laplace tidal equations (LTE), which describe its barotropic response, that is the component associated with the propagation of surfacegravity waves. This is typically the tidal response that the ocean will have if it is unstratified (uniform density). Nevertheless, solutions to the barotropic equations can be applied to the case of stably stratified oceans to describe the contribution of internal gravity waves, restored by the Archimedean force. As was shown a long time ago (e.g. Chapman & Lindzen 1970, in the case of the atmosphere), the tidal response of a stratified fluid layer in the classical tidal theory is characterized by a frequencydependent equivalent depth, which is analogous to the barotropic ocean depth (see e.g. Tyler 2011). Results obtained in the barotropic approximation can thus be reinterpreted by replacing the ocean depth by the equivalent depth associated with the studied mode. However, the complete solution controlling the tidal response of a deep stratified ocean, including its vertical component, has never been explicitly resolved nor analysed to our knowledge.
Therefore, following along the line of Webb (1980) and Tyler (2011), we develop here a linear global model of oceanic tides including threedimensional (3D) effects related to the ocean sphericity and vertical structure. This model, by taking into account the vertical stratification of the ocean with respect to convection (i.e. its convective stability), allows us to quantify in a consistent way the contribution of internal gravity waves to the oceanic tidal response through explicit solutions derived in simplified cases. The obtained results can be used as a tool to (1) characterize the tidal behaviour of deep global oceans, (2) constrain the structure and history of observed telluric planets from the architecture of their host planetary system, and (3) provide a simplified diagnosis of the oceanic tidal dissipation rate.
In Sect. 2, we establish the equations governing the dynamics of the 3D oceanic tidal response, identify the possible tidal regimes, and express the tidal Love numbers and torque exerted on the planet in the general case. In Sect. 3, we simplify the dynamics by assuming a uniform stratification as a first step and derive an analytic solution for the tidal response and the resulting dissipation rate. This solution is then applied to idealized Earth and TRAPPIST1 f planets in Sect. 4 in order to illustrate the difference between shallow and deep oceans. We also use this solution in Sect. 5 to explore the dependence of the dissipation rate on the ocean depth and stratification. Finally, we discuss the simplifications assumed in the modelling in Sect. 6 and give our conclusions in Sect. 7.
2 Dynamics of a gravitationally forced thick ocean
We establish in this section the equations that govern the dynamics of tides in a thick oceanic shell submitted to the gravitational tidal potential of a perturber. The considered system is a spherical telluric planet of external radius R uniformly covered by an oceanic layer of depth H (Fig. 1). The radius of the solid part is designated by R_{c} = R − H. We assume that the ocean rotates uniformly with the rocky part at the angular velocity Ω, the corresponding spin vector being denoted Ω. This allows us to write the fluid dynamics in the natural equatorial reference frame rotating with the planet, , where O is the centre of the planet, X_{E} and Y _{E} define the equatorial plane and Z_{E} = Ω∕Ω. We use the spherical basis and coordinates where r is the radial coordinate, θ the colatitude and φ the longitude (the position vector is r = r e_{r}). Finally, the time is denoted t.
Fig. 1 Quadrupolar semidiurnal oceanic tide in a rotating terrestrial planet gravitationally forced by a perturber (star of satellite). Reference frame and system of coordinates used in the modelling. The spin frequency of the planet and the orbital frequency of the perturber are denoted Ω and n_{orb}, respectively. The notation δ_{tide} designates the angular lag of the tidal bulge. 
2.1 Forced dynamics equations
The planet generates a gravity field, denoted g, oriented radially. The centrifugal acceleration due to rotation tends to make the effective acceleration depend on latitude. We ignore this dependence by considering that Ω ≪Ω_{c}, where represents the Keplerian critical rotation velocity (the rotation velocity of the planet is far slower than the critical velocity for which the distortion due to the centrifugal acceleration destroys the oceanic layer). We assume then that the ocean is stratified radially in pressure p and density ρ. These quantities are written (1)
where the superscript _{0} refers to the background distribution and δ to a small fluctuation at the vicinity of the equilibrium. Similarly, the velocity of the flows generated by the perturbation is denoted and the associated displacement ξ (such that V = ∂_{t} ξ, Unno et al. 1989). We note that, because of the solid rotation assumption, there is no meridional or zonal background circulation. We consider in the following that the tidal perturbation can be approximated by a linear model and therefore ignore terms of order greater than 1 with respect to ξ, V, δp and δρ. Moreover, we assume the Cowling approximation (Cowling 1941), that is, the selfgravitational effect of the variation of mass distribution is not taken into account. Finally, following early works (e.g. Webb 1980; Egbert & Ray 2001, 2003; Tyler 2011), we introduce dissipation by using a Rayleigh friction, σ_{R} V, where σ_{R} is a constant effective frequency associated with the damping (we note that Egbert & Ray 2001, 2003, consider two different models to describe friction, the Rayleigh friction model and a quadratic one depending on the total velocity of the flow and including all tidal components). As demonstrated by Ogilvie (2009), the results obtained with such a simplified approach are a reasonably good approximation of those obtained by taking into account the complete viscous force. It also seems to be a reasonable approach for situations such as exoplanets, where the solidfluid coupling is naturally unknown.
Under these assumptions, the NavierStokes equation describing the distortion of the ocean forced by the tidal potential U can be expressed in the frame rotating with the planet as (Gerkema & Zimmerman 2008) (2)
It is completed by the equation of mass conservation (3)
and the equation of buoyancy (4)
where designates the sound velocity, Γ being any state function and S_{0} the salinity of the fluid. It shall be noted here that we only consider isohaline and isentropic processes. Diabatic processes such as the thermal expansion of the ocean are ignored. The radial stratification of the oceanic layer is characterized by the Brunt–Väisälä frequency (e.g. Gerkema & Zimmerman 2008), (5)
where we have introduced the nondimensional reduced altitude x, defined by r = R_{c} + Hx. The three components of the momentum equation given by Eq. (2) are strongly coupled by the Coriolis acceleration. To make the problem tractable analytically, we assume the commonly used traditional approximation (e.g. Eckart 1960; Mathis et al. 2008; AuclairDesrotour et al. 2017a). This simplification consists in ignoring the latitudinal component of the rotation vector in the Coriolis acceleration (i.e. 2Ω sinθV _{r} and 2Ω sinθV _{φ} along e _{φ} and e _{r} respectively), giving (6)
thus making it possible to separate θ and r coordinates in solutions, as shown hereafter. As discussed in Sect. 6, the traditional approximation is usually considered to be satisfactory in stably stratified layers (2Ω ≪ N^{2}) and for tidal frequencies satisfying the condition 2Ω < σ ≪ N^{2}, which corresponds to the regime of superinertial waves (e.g. Mathis 2009; AuclairDesrotour et al. 2018).
Finally, the equations of mass conservation (3) and of the buoyancy (4) can be written respectively (7)
The tidal perturbation is periodic in time and longitude. Therefore, introducing the longitudinal wavenumber m and the tidal frequency σ, we expand the perturbed quantities in Fourier series of t and φ. Any quantity f is written (9)
where the f^{m,σ} are the spatial distributions of f associated with the doublet . In the case of a perfect fluid (σ_{R} = 0), the latitudinal structure of the perturbation is fully characterized by m and the real parameter ν = 2Ω∕σ, which is the socalled spin parameter (e.g. Lee & Saio 1997). Viscous friction slightly modifies the latitudinal structure by introducing a dependence on σ_{R}. Thus, we define the complex tidal frequency and spin parameter as (10)
and the latitudinal operators (11)
As mentioned above, the traditional approximation allows us to separate the coordinates r and θ in the Fourier coefficients of Eq. (9). Therefore, we write these functions as (13)
where n designates the latitudinal wavenumberof a component ( if ν ≤ 1; if ν >1), the are the socalled Hough functions (Hough 1898), defined on the interval , and the and stand for the latitudinal functions associated with the latitudinal and longitudinal velocities and displacements, respectively. Hough functions are the solutions of the eigenfunctionseigenvalues problem defined by the Laplace’s tidal equation (Laplace 1798), (14)
where Θ is a function, and the operator (15)
Hence, for a given n, is such that , the parameter being the associated eigenvalue. We note here that the only difference with the usual case where friction is not taken into account (e.g. Lee & Saio 1997) is the nature of which is a complex number and not a real one (see Eq. (10)). It follows that the and are complex functions and parameters in the general case (e.g. Volland 1974a,b). Looking at the expression of , one can identify two dissipation regimes:

A weakly frictional regime (σ≫ σ_{R}), where Hough functions and eigenvalues are real (they are denoted and ). This regime corresponds to the case treated usually in the theory of tides (Chapman & Lindzen 1970; Lee & Saio 1997; AuclairDesrotour et al. 2014). Tidal modes can be divided into two families: the socalled gravity modes (n ≥ 0), which are defined both in the regime of superinertial waves (ν≤ 1) and in the regime of subinertial waves (ν > 1), and the inertial modes (n < 0), defined only in the regime of subinertial waves (see Lee & Saio 1997). The Hough functions associated with gravity modes degenerate in the associated Legendre polynomials^{1} (with l = m + n) when ν → 0.

A strongly frictional regime (σ≲ σ_{R}), where friction significantly affects the horizontal structure of tidal waves. In this regime, the hierarchy of eigenvalues can change, as demonstrated by Volland (1974a) who notes that the corresponding critical points appear for σ~ σ_{R}. For σ_{R} ≳Ω, inertial modes converge towards gravity modes. Gravity and Rossby Hough functions merge asymptotically when σ_{R} → +∞, and converge to the associated Legendre polynomials. Thus, a strong friction makes the Coriolis effects negligible, as if the planet were not rotating. In this asymptotic regime, the angular lag between the tidal bulge and the direction of the perturber is given by the imaginary part of the vertical profiles of perturbed quantities. The behaviour of tidal modes in the dissipative regime is discussed thoroughly by Volland & Mayr (1972) and Volland (1974a).
Hough functions associated with symmetric modes of degree m = 2 are plotted in Fig. 2 for various values of the ratio σ_{R} ∕σ. The equations describing the vertical structure are obtained by substituting the expansions given by Eqs. (9) and (13) in Eqs. (6–8). In order to lighten expressions, the superscripts is omitted up to the end of Sect. 2.1, where no confusion will arise. Let us introduce the notation y_{n} = δp_{n}∕ρ_{0} and assume the ocean to be at the hydrostatic equilibrium (dp_{0}∕dx = −Hgρ_{0}), where x is the reduced altitude introduced in Eq. (5). After some manipulations, one obtains (16)
This system is then written as a single secondorder ordinary differential equation, (18)
Here, we have introduced the frequencydependent, dimensionless acoustic, and sphericity parameters (22)
The acoustic parameter can be written , where is the cutoff frequency of acoustic waves (also called the Lamb frequency) for the degreen mode when . The regime of acoustic waves thus corresponds to ε_{s;n}≳ 1. Typically, for an Earthlike planet of radius R_{⊕} = 6371 km and angular velocity Ω_{⊕} = 7.25 × 10^{−5} s^{−1} (NASA fact sheets), with a water oceanic layer characterized by c_{s} ≈ 1.545 km.s^{−1} (Gerkema & Zimmerman 2008), the gravest acoustic modes of the quadrupolar semidiurnal tide (m = 2, n = 0, ) appear for Ω ≳5.5 Ω_{⊕}. The other parameter, K_{°}, corresponds to the effect of the spherical geometry of the layer on the structure of tidal waves. It is usually ignored in 2D modellings where r = R. We see in the following section that K_{°} ∝ H∕R when H∕R ≪ 1. Eventually, using the change of variables , where is the radial function (23)
we write thevertical structure equation (Eq. (18)) as the Schrödingerlike equation (24)
in which represents the vertical wavenumber of the mode, defined as (25)
The vertical profiles of the other quantities are deduced from straightforwardly thanks to polarization relations. We get, for the velocity field
with the frequencydependent factors
Through equations describing the tidal dynamics, we identify the families of waves involved in the tidal response. First, because of rotation, inertial waves are generated in the frequency range σ < 2Ω. Then, we identify surfacegravity waves, which are restored by gravity and are characterized by the cutoff frequency . If the ocean is stably stratified (N > 0), internal gravity waves can propagate inthe frequency range σ < N. These wavesare restored by the Archimedean force. Finally, in the highfrequency range delimited by the acoustic cutoff frequency σ_{s} = c_{s}∕r, the tidal response is partly composed of horizontally propagating acoustic Lamb modes, restored by compressibility. The possible regimes of oceanic tides are determined by the hierarchy of the characteristic frequencies of the system, namely σ, σ_{R}, 2Ω, N, σ_{g} and σ_{s}. The frequency spectrum of waves composing the oceanic tidal response is summarized by Fig. 3.
Fig. 2 Real (top) and imaginary (bottom) parts of even Hough functions (n even) associated with a quadrupolar tidal perturbation (m = 2) and defined by with ν =1 (ordinary/gravity modes) and various positive values of γ = σ_{R}∕σ. From left to right, γ = 0 (nonfrictional regime), γ = 0.1 (weakly frictional regime), γ = 1 and γ = 10 (stronglyfrictional regime). Hough functions are plotted as functions of the latitude δ (degrees). 
2.2 Tidal potential, Love numbers, and tidal torque
The variation of mass distribution due to the tidal distortion modifies the selfgravitational potential of the planet. The tidal potential of the excitation is usually expanded in Fourier series and spherical harmonics (see for instance the Kaula’s multipolar expansion; Kaula 1962). It is thus convenient to write the fluctuations of the selfgravitational potential, denoted , in the same form, (36)
The radial profiles of the (m, σ)components are themselves written (37)
where we have introduced the contribution of the surface displacement and of the internal density variations . In the following, we use the subscripts ξ and ρ to make the distinction between the two contributions in a systematic way. At the planet surface (x = 1), the two components of the potential are expressed as
with representingthe universal gravity constant. Expanded in Hough functions, they write
In these expressions, the complex weighting coefficients are expressed as (42)
where the and designate the mutual projection coefficients of the Hough functions on normalized associated Legendre polynomials, defined respectively by the expansions
The coefficients quantify the coupling induced by the Coriolis effects in the dynamics of the tidal response. They are real in the nonfrictional case (σ_{R} = 0 and ), where and , the notation ⟨⋅, ⋅⟩ standing for the scalar product defined for any and as (45)
In the case of a nonrotating planet (Ω = 0 and ν = 0), there is no coupling. Therefore, if l = n + m = k, else . The notations and in Eqs. (40) and (41) stand for the vertical displacement and density variations due to the component of the forcing projected on the Hough function. The expressions of the components of the tidal potential (Eqs. (40) and (41)) allow us to compute complex Love numbers, which are commonly used to quantify tidal dissipation in celestial bodies (Tobie et al. 2005; Remus et al. 2012; Ogilvie 2014). Love numbers are defined as the ratio between the gravitational potential due to the distortion and the tidal gravitational potential taken at the external surface of the body (r = R). Thus, oceanic Love numbers associated with the triplet , denoted , are defined as (46)
Here, and stand for the tidal response of the ocean to the component of the excitating tidal potential. The are provided by the Kaula’s multipolar expansion of the tidal gravitational potential and are expressed as functions of the Keplerian elements of the planetperturber system (Kaula 1966; Mathis & Le PoncinLafitte 2009).
Becauseof internal dissipation, the variation of mass distribution due to the tidal perturbation generates a tidal torque. The torque exerted by the perturber on the oceanic layer with respect to the spin axis of the planet, denoted , is given by (47)
The components and of the mode are defined as
where ^{*} stands for the conjugate of a complex number and ℜ its real part (ℑ will be used here for the imaginary part), denotes the density at the surface ocean, the spatial domain filled by the oceanic shell at rest, its upper boundary, and dS and infinitesimal surface and volume parcels, respectively. Similarly to the tidal potential due to the distortion, these two components can be expanded on Hough functions and associated Legendre polynomials. For l ≥ m and k ≥ m, one thus obtains,
Fig. 3 Frequency spectrum of waves likely to be excited by the tidal perturbation and of the corresponding tidal regimes. The characteristic frequencies of the system are, from left to right, the drag frequency σ_{R}, the inertia frequency 2Ω, the surfacegravity waves cutoff frequency , the Brunt–Väisälä frequency N, and the acoustic cutoff frequency σ_{s} = c_{s}∕r. Their hierarchy can vary as a function of the physical and dynamical properties of the planet. In the case of the Earth’s ocean, σ_{R} ≈ 10^{−5} s^{−1} (Webb 1980), 2Ω = 1.46 × 10^{−4} s^{−1}, σ_{g} = 3.1 × 10^{−5} s^{−1}, N ≈ 2.2 × 10^{−3} s^{−1} (Gerkema & Zimmerman 2008), and σ_{s} = 6.3 × 10^{−3} s^{−1}. 
3 Tides in a thin stably stratified ocean
In some observed cases, the oceanic layer of terrestrial planets and satellites is thin compared to the radius of the body (for instance, the Earth’s ocean has an aspect ratio of 6 × 10^{−4}; Eakins & Sharman 2010). Therefore, in this section, we apply our general modelling to the simplified case of an oceanic layer of small depth H ≪ R, and compute an analytic solution of the oceanic tidal response, Love numbers, and tidal torque.
3.1 Tidal waves dynamics
The thin shell hypothesis allows us to simplify the physical setup. Given that the relative difference ofgravity between the lower and the upper boundaries is proportional to the ratio H∕R, g is supposed to be constant in the oceanic layer. Similarly, the radial variations of the tidal gravitational potential with x can be neglected (d U_{n}∕dx ≈ 0 and d ^{2} U_{n}∕dx^{2} ≈ 0). Furthermore, for simplicity, we assume uniform stratification and compressibility, that is, that the Brunt–Väisälä frequency (N) and the sound velocity (c_{s}) are constants. We note that the case of the Earth’s ocean does not well follow this assumption since N decays over several orders of magnitude from the pycnocline (the upper region of the ocean where the density is changing most rapidly) where its typical values are about 0.01 s^{−1}, to the abyssal ocean, where they fall to 0.001 s^{−1} (e.g. Gerkema & Zimmerman 2008). Although a more realistic model would be more appropriate in this case, the uniform stratification appears as a useful first step to solve the vertical structure equation and derive explicit solutions controlling the tidal response of a stratified ocean. As shown by Eq. (5), setting N^{2} to a constant implies a density decreasing exponentially with altitude, (52)
where τ designates the decreasing rate given by (53)
The acoustic and sphericity terms (Eq. (22)) are simplified into (54)
It follows that Φ_{n} (see Eq. (22)) writes (55)
Hence, we express the vertical structure equation (Eq. (24)) (56)
The vertical wavenumber to square , given by Eq. (25) in the general case, is now the constant (57)
In this expression, the first term is the vertical wavenumber of internal gravity waves; it dominates when , that is, for a stable stratification. Other terms correspond to acoustic and sphericity terms usually ignored in anelastic (c_{s} = +∞) and thinshell (H∕R → 0) approximations. Here, we will only ignore sphericity terms, that is, we assume that K_{°} = 0. We will keep acoustic terms because they do not bring supplementary mathematical complexities into the analytic treatment and can, furthermore, be comparable to those associated to stratification. For instance, in the case of the Earth’s ocean, density increases from 1022 kg m^{−3} at the surface to 1070 kg m^{−3} at 10 km depth (Gerkema & Zimmerman 2008), which gives the mean gradient d ρ_{0} ∕dz = −0.0048 kg m^{−2}. As , the two terms of Eq. (5) are comparable and must be retained. To solve the vertical structure equation, two boundary conditions are required. At x = 0, we set the impenetrable rigidwall condition ξ_{r} = 0. At the upper boundary, we apply the usual freesurface condition (e.g. Unno et al. 1989), δp = gρ_{0}ξ_{r}. It follows that (58)
where is the constant given by (59)
and and the dimensionless coefficients expressed as
3.2 Secondorder tidal Love number and tidal torque
By substituting Eq. (58) in Eqs. (29) and (33) we obtain the components of the variation of mass distribution intervening in the Love numbers and tidal torque (see Eqs. (50) and (51)) as explicit functions of the internal structure parameters, (62)
where the frequencydependent parameters and are expressed as
with γ = δ − τ and the dimensionless coefficients (65) (66)
Hence, stands for the intrinsic response of the planet due to the variations of the oceanic surface level, while corresponds to the contribution of internal inertialgravity waves. This will later be equal to zero in the case of an incompressible and neutrally stratified ocean (τ = 0). The contribution of internal gravity waves to the tidal response with respect to that of surfacegravity waves is weighted by the ratio . By using the expressions of the solution, we retrieve the weighting factor given in the literature (see e.g. Hendershott 1981, Eq. 10.40), that is, . The oceanic secondorder Love number is then deduced from Eqs. (40), (41) and (46) straightforwardly. In the quadrupolar approximation, where terms of degrees l > 2 are neglected, it writes (70)
where M_{oc} = 4πR^{2}Hρ_{s} designates the total mass of the ocean in the weak stratification approximation (τ ≪ 1). Similarly, the tidal quality factor (Mathis & Le PoncinLafitte 2009) is expressed as (71)
and the quadrupolar tidal torque, which corresponds to the l = m = 2 component, as (72)
where the quadrupolar component of the tidal potential is given by (73)
Expressed as a function of the degreetwo tidal Love number given by Eq. (70), the tidal torque is written (74)
which is the wellknown expression of the tidal torque associated with the quadrupolar semidiurnal tide (e.g. Efroimsky & Williams 2009; Makarov 2012; Correia et al. 2014).
Fig. 4 Vertical displacement due to the quadrupolar semidiurnal tide (, n_{orb} being the dynamical frequency) obtained using the analytic solution established for the thin layer approximation (Eq. (58)). The displacement, normalized by the quadrupolar tidal potential, is plotted in the equatorial plane of the planet as a function of longitude (horizontal axis) and normalized altitude (vertical axis). The position φ = 0 corresponds to the subperturber point; the oceanic floor and surface are located at x = 0 and x = 1, respectively. For this computation, the rotation period of the planet P_{spin} = 2π∕Ω and orbital period of the perturber P_{orb} = 2π∕n_{orb} are set to P_{spin} = 6.0 d and P_{orb} = 365.25 d. The used values of parameters are R = R_{⊕}, H = 100 km, g = 9.81 m s^{−2}, N = 10^{−3} s^{−1}, c_{s} = 1545 m s^{−1} and σ_{R}= 10^{−5} s^{−1}. 
3.3 Case of the neutrally stratified ocean (N = 0)
In analytic solutions describing oceanic tides, the stratification is usually not taken into account, the layer being 2D and the fluid assumed to be incompressible (e.g. Webb 1980; Tyler 2011, 2014; Chen et al. 2014; Matsuyama 2014). This reduces the oceanic tidal dynamics to horizontal flows. We show here that we recover the results given by this approach with our modelling, in which the shallowwater case is a particular case.
Following the early works mentioned above, let us set N = 0 (neutral stratification) and c_{s} = +∞ (incompressible fluid). The vertical wavenumber thus becomes (75)
As shown byEq. (75), is now directly proportional to the horizontal wavenumber . In the weakfriction approximation (σ_{R} ≪σ) and the regime of superinertial waves (), it does not depend on the tidal frequency. Indeed, in this case, the associated eigenvalues are with (the relation between the degree l of associated Legendre polynomials and the degree n has been defined as l = m + n). The frequency dependence of cannot be ignored any more in the regime of subinertial waves. Because of the thinlayer approximation, . It follows that and . Hence, in the case of the neutrally stratified incompressible ocean, the analytic solution of Eq. (58) simply reduces to (76)
where and are the frequencies of resonance of largescale () damped surface inertialgravity waves (see Fig. 3) associated with the gravity mode of degree n. These frequencies characterize the tidal response of the 2D incompressible ocean, where only surfacegravity waves can propagate (e.g. Webb 1980, Eq. (2.12)). They are expressed as (77)
We recognize in this expression the dispersion relation of largescale surfacegravity waves, (see e.g. Vallis 2006, in the adiabatic limit). Like the vertical wavenumber, and formally depend on σ through but this dependence can be neglected in the regime of superinertial waves if the weakfriction approximation is assumed. The above approximations imply τ = 0 and . Thus, the Love numbers and tidal torque (Eqs. (70) and (72)) are determined by the contribution of the surface displacement solely, which reads (78)
The expression of shows that, in the absence of stable stratification, the coupling induced by the Coriolis effect will lead to a set of resonances (one per Hough mode) enhancing the tidal dissipation and tidal torque at frequencies .
4 Application to illustrative cases
In this section, we apply the results established above to representative telluric planets by using the thinshell approximation of Sect. 3. We first treat the case of the Earth, which illustrates the thinlayer asymptotic behaviour. However, we shall bear in mind that the global ocean approximation is a rough hypothesis in the case of the Earth, where the interaction of the tidal perturbation with continents plays a central role (Webb 1980; Egbert & Ray 2001, 2003). We then apply the model to TRAPPIST1 f, which could be covered with a deep global ocean of liquid water owing to its small density (M_{p} ≈ 0.36 M_{⊕} and R ≈ 1.045 R_{⊕}, see Wang et al. 2017).
4.1 The Earth
The Sun and the Moon exert on the Earth gravitational forcings of comparable intensities. The resulting tidal perturbation is the sum of two contributions, the solid and oceanic tides, which drive the rotational evolution of the planet as well as the orbital evolution of the EarthMoon system. We focus here on the Lunar quadrupolar tide, which corresponds to the tidal frequency , where Ω and n_{orb} stand for the spin frequency of the Earth and orbital frequency of the Moon. We introduce the corresponding rotation (P_{rot} = 2π∕Ω) and orbital (P_{orb} = 2π∕n_{orb}) periods. Following Webb (1980), we use for the uniform oceanic depth the mean depth of the Earth’s ocean, that is H = 4 km (Webb 1980). The radius of the planetis set to R = R_{⊕} km (with R_{⊕} = 6378 km), the surface gravity and density to g = 9.81 m s^{−2} and ρ_{s} = 1022 kg m^{−3}. The Earth’s ocean is characterized by N ~ 10^{−4} − 10^{−2} s^{−1} and c_{s} ≈ 1545 m s^{−1} (Gerkema & Zimmerman 2008), which means that τ ≪ 1. Thus, the effects of stratification are negligible in this case. We set the Brunt–Väisälä frequency to the intermediate value N = 10^{−3} s^{−1}. All these quantities are summarized in Table 1.
The drag frequency characterizing the effective Rayleigh friction (σ_{R}) is more difficult to specify because it models the effects of several mechanisms, such as turbulent friction, viscous friction, friction with topography, and breaking of tidal waves (Garrett & Munk 1979; Garrett & Kunze 2007). Thus, we begin by studying the dependence of the oceanic tidal response on σ_{R}. The frequency spectra of the imaginary part of the tidal Love number (Eq. (70)) and of the associated tidal quality factor (Eq. (71)) are plotted in Fig. 5 (left column) as functions of the normalized frequency (Ω _{⊕} stands for the today Earth rotation rate) for various orders of magnitude of σ_{R}. We observe on these plots the resonances associated with surfacegravity modes modified by rotation. They correspond to the eigenfrequencies given by Eq. (77). As σ_{R} decreases, the variability of increases. Particularly, the level of the nonresonant background decreases proportionally to σ_{R}, while the peaks heights increase as , which is in good agreement with the scaling laws derived in Auclair Desrotour et al. (2015). In the asymptotic regime of strong friction, the behaviour of the ocean is regular. The imaginary part of the tidal Love number scales as in the zerofrequency limit (σ → 0) and as at σ → +∞.
The observed receding of the Moon is estimated to 3.82 cm per year (Dickey et al. 1994; Bills & Ray 1999). This corresponds to for the tidal Love number and for the tidal torque, which are the values obtained by setting the Rayleigh drag frequency to σ_{R} = 10^{−5} s^{−1} in our model (Fig. 5, right panel). We hence retrieve the ~ 30 h friction time scale estimated by Webb (1980) with a similar approach. This value of σ_{R} is used in Sect. 5 to explore the domain of parameters.
Fig. 5 Semidiurnal oceanic tide of the Earth (left column) and TRAPPIST1 f (right column). The logarithms of the imaginary part of the tidal Love number (left) and tidal quality factor (right) are plotted as functions of the normalized forcing frequency (where Ω_{⊕} designates the today rotation rate of the Earth) for various orders of magnitude of the drag parameter . These frequency spectra are computed using the expressions given by Eqs. (70) and (71). In each case, the parameter n_{orb} is assumed to be constant and the rotation rate of the planet varies with the tidal frequency following the formula Ω = n_{orb} + σ∕2. The resonances associated with surface inertialgravity modes are designated by black dashed lines in the top panels and numbers indicate the degree n of the corresponding Hough modes and the sign of eigenfrequencies given by Eq. (77). The values of parameters used for this evaluation are summarized in Table 1. 
4.2 TRAPPIST1 f
To illustrate the behaviour of a deep, stably stratified ocean, we consider the case of TRAPPIST1 f. TRAPPIST1 f is one of the eight telluric planets recently discovered in the vicinity of the star TRAPPIST1 (Gillon et al. 2017; Wang et al. 2017). Its radius is very close to that of the Earth (R = 1.045 M_{⊕}) whereas its mass is only M = 0.36 M_{⊕}. As a consequence, its mean density is equal to , which is far smaller than those of rocky planets (for instance, the mean density of the Earth is ; see NASA fact sheets). This means that water could stand for an important fraction of the planet mass. Typically, the water fraction is estimated to be between 25 and 100% of the planet mass. Moreover, with a greenhouse effect, TRAPPIST1 f is sufficiently irradiated by its host star to have a surface temperature compatible with liquid water (its black body equilibrium temperature is estimated to 219 K, see Wang et al. 2017). All these features argue for the possible presence of a deep oceanic layer on TRAPPIST1 f. In such a case, if the layer is stably stratified, the effect of stratification cannot be ignored any more because the variation rate of the density profile (τ defined by Eq. (53)) is not necessary very small with respect to 1. Thus, as in the case of the Earth, the tidal response is mainly composed of surfacegravity waves modified by rotation. But it is also composed of internal gravity waves, restored by the Archimedean force and inducing internal density variations (Fig. 3).
We investigate this complex behaviour by considering an idealized TRAPPIST1 f planet with a global ocean of depth H = 1000 km and a stratification similar to that of the Earth’s ocean, that is N = 10^{−3} s^{−1}. Assuming the orbital period of the planet (n_{orb}) to be constant, we study the frequency dependence of the quadrupolar component of the semidiurnal stellar tide (cf. Table 1) by letting Ω vary with σ. Similarly to the previous case, the corresponding spectra of and Q are plotted on Fig. 5 (right column) as a function of the forcing frequency. In the highfrequency range we retrieve the resonances associated with surfacegravity modes. The two peaks appearing around ω ≈±4.5 correspond to the Hough mode of degree n = 0, which is the surfacegravity mode of lowest eigenfrequency. Given that in the weakly frictional regime (see Eq. (77)), the resonances are translated towards the highfrequency range with respectto the case of the Earth (left column). In addition to surfacegravity modes, we observe resonances caused by the propagation of internal gravity waves in the lowfrequency range. These modes can exist for σ ≲ N (see Fig. 3), which corresponds to ω≲ 5.
5 From shallow to deep oceans
In the previous sections, we identified the depth of the ocean and the stability of its stratification as the key structure parameters characterizing the oceanic tidal response. The depth determines the eigenvalues of resonant surfacegravity modes and the Brunt–Väisälä frequency, the frequency range of internal gravity modes. Further, the contribution of these modes is weighted by the dimensionless number N^{2}H∕g, which compares the Archimedean force to the gravity force. Therefore, we examine in this section the sensitivity of tidal dissipation to H and N. We consider the case of idealized Earth sisters hosting a global ocean and submitted to the semidiurnal gravitational forcing of any perturber orbiting in their equatorial plane. These planets are characterized by various orders of magnitude of ocean depths (H = 10, 100, and 1000 km) and Brunt–Väisälä frequency (N = 10^{−4}, 10^{−3} s^{−1}), all the other parameters being unchanged (see Table 1). Following the discussion of the previous section, the frequency associated with the Rayleigh drag is set to σ_{R} = 10^{−5} s^{−1} in a first case (friction time scale of the Earth’s ocean). In a second case, it is set to σ_{R} = 10^{−6} s^{−1}, which corresponds to a weaker friction, to broaden the explored parameter space. Hence, for each planet, we plot in Fig. 6 the imaginary part of the tidal Love number as a function of the orbital period of the perturber and the rotation period of the planet. Then, assuming that the perturber is a Moonlike satellite (M_{pert} = 7.346 × 10^{22} kg), we plot the corresponding tidal torque exerted on planets in Fig. 7.
As it may be noted, the three maps plotted for H = 10 km (bottom panels of Figs. 6 and 7) have the same aspect whatever the value of N. This can be interpreted as the fact that the tidal response generated by the gravitational forcing exerted on the ocean is mainly composed of surfacegravity waves. In other words, the tidal response is dominated by its barotropic component, which causes the role played by the stratification to appear negligible in this case, as mentioned before. We identify the resonances due to surface inertialgravity modes observed in the case of the Earth in the previous section (Fig. 5, left column). They form a typical wings pattern, which is proper to these waves (see e.g. Fuller & Lai 2014, Fig. 7, in the case of white dwarves). However, we note that this pattern is not symmetrical with respect to the diagonal line designating synchronization. This is a consequence of the effects of the Coriolis acceleration. In the supersynchronous regime (area located below the diagonal line), they induce a distortion that couples the quadrupolar forcing to several Hough modes. These inertialgravity waves tend to be mixed with acoustic waves while P_{rot} → 0. In the subsynchronous regime (area located above the diagonal line), the forcing frequency is greater than the inertia frequency (superinertial regime) and the only mode coupled with the forcing is the gravity mode of degree n = 0. As a consequence, only one resonance appears in the superinertial regime. As the tidal potential scales as (see Eq. (73)), the tidal torque scales as (see Eq. (72)), which corresponds to the horizontal colour gradient of the maps plotted in Fig. 7.
By moving to deeper oceans, we observe that the resonances are translated towards the smallperiod range, following the scaling law identified above. For N = 10^{−3} s^{−1}, the stratification is sufficiently strong to allow internal gravity waves to propagate. Therefore, new resonances appear in the frequency range σ≲ N (see right panels of Figs. 6 and 7). This effect is enhanced for H = 1000 km. As the eigenfrequencies of internal gravity waves do not depend on Ω like those ofsurface gravity waves, their wings patterns are symmetrical with respect to the diagonal axis. Moreover, given that the contribution of internal waves is weighted by the dimensionless number N^{2} H∕g, the nonresonant background is amplified in the frequency range σ≲ N. The effects of the resonances associated with internal gravity waves are attenuated by the Rayleigh drag. When σ_{R} → 0, resonances tend to increase the dependence of the tidal torque and Love number on the forcing frequency (as expected from the general scaling laws derived in Auclair Desrotour et al. 2015).
6 Discussion
In this work, we have opted for a simplified linear analysis allowing us to widely explore the domain of parameters at a reasonable computational cost. This approach is convenient to examine the frequencyresonant behaviour of a global ocean and highlight the key parameters of the problem, which can be explained in more detail in a second phase. Particularly, we identify and quantify in a consistant way the contribution of internal gravity waves restored by the stable stratification. This contribution can be important in the case of deep oceans, where the combined effects of tides and stratification induce important density fluctuations. However, we shall discuss here the main simplifications assumed in the model:

Global ocean of uniform depth – The large scale topographical features of the oceanic floor and continents are not taken into account. Because of this spherical symmetry, the obtained oceanic tidal response is very regular and exhibits resonances corresponding to global spherical modes (see e.g. Fig. 5). The case of a hemispherical ocean centred at the equator was treated in early studies in the shallowwater approximation (Proudman & Doodson 1936; Doodson 1938; LonguetHiggins & Pond 1970; Webb 1980). Meridian continental shelves prevent largescale gravity waves from propagating. They thus couple the global gravitational forcing to modes determined by the scale of the basin. This significantly modifies the aspect of the frequency spectrum of the tidal response (see e.g. Webb 1980).

Rayleigh friction – All of the dissipative mechanisms damping the tidal response are reduced to a single parameter, the effective Rayleigh coefficient σ_{R}. This approximation is common in the literature (e.g. Webb 1980; Egbert & Ray 2001, 2003; Ogilvie 2009; Tyler 2011) and assumed for convenience. In addition to the traditional approximation, discussed in the following paragraph, it allows us to separate the x and θ coordinates in the dynamics. Rayleigh friction can be reasonably used to describe the drag caused by smallscale topographical features homogeneously distributed around the planet, as in the numerical model of Egbert & Ray (2001). It can also be considered a simplified approximation of viscous and turbulent frictions if a lengthscale L such that ΔV ~ V∕L^{2} is introduced. However, it does not at all model the local dissipation due the breaking of waves on continental shelves although this effect stands for the major part of the tidal energy in the case of the Earth.

Traditional approximation – This simplification is commonly used in the literature to solve the latitudinal and vertical structure of the tidal response separately (e.g. Unno et al. 1989). It consists in neglecting the latitudinal component of the rotation vector in the Coriolis acceleration. This means that we ignore both the radial component of the Coriolis force and the horizontal component of the Coriolis force associated with radial motions. Consequently, this approximation is appropriate to 2D models, where vertical motions are negligible relative to horizontal ones. In the case of deeper oceans, its domain of applicability depends on the hierarchy of the inertia (2Ω), Brunt–Väisälä (N), and forcing (σ) frequencies. In the superinertial asymptotic regime (2Ω≪σ), the forcing time scale is small compared to the rotation period, which makes the traditional approximation appropriate. However, when this condition is not satisfied, motions are fully tridimensional in the case of a neutrally stratified ocean (N ≈ 0). This requires one to use 3D (2D if solutions ∝ e^{imφ} are assumed) numerical or semianalytical methods (see e.g. Ogilvie & Lin 2004). The Archimedean force acts as a restoring force on fluid particles in the vertical direction. Therefore, if this restoring force is sufficiently strong, it prevents vertical motions from being comparable in amplitude to horizontal ones. This condition is expressed as 2Ω≪ N (e.g. AuclairDesrotour et al. 2017a). The validity of the traditional approximation has been examined in various fluid layers, from planetary oceans and atmospheres (e.g. Gerkema & Shrira 2005; Gerkema & Zimmerman 2008) to stellar interiors (e.g. Mathis 2009; Tort & Dubos 2014; Prat et al. 2017).

Solid rotation – As we consider in our model that the oceanic layer rotates uniformly with the planet at the angular velocity Ω, we ignore the complex interplay existing between tidal waves and mean flows (e.g. Favier et al. 2014; Guenel et al. 2016). By studying the propagation of internal gravity waves in a shear flow, Booker & Bretherton (1967) showed that the amplitude of oscillations is attenuated by a factor depending on the local Richardson number of the fluid, (where V_{0} stands for the velocity vector of the shear flow), as waves pass through a critical level at which their horizontal phase speed is equal to the zonal velocity of the flow. By this mechanism, tidal waves can transfer horizontal momentum to mean flows, which modifies the tidal response.

Freesurface boundary condition – In this work, a freesurface boundary condition is used to derive analytic solutions. This condition corresponds well to the case of external oceanic layers, where the upper surface can move freely. However, it is not appropriate to subsurface oceans such as those presumedly hosted by Europalike icy satellites (e.g. Khurana et al. 1998; Kivelson et al. 2002). In order to study these objects, results obtained in this work can easily be applied to the case of a rigid lid by replacing for instance the standard freesurface condition by a rigidwall condition.

Coupling with the solid part – In order to simplify the analysis, the solid part of the planet has been considered as a sphere of infinite rigidity. In reality, the solid part behaves as a viscoelastic body (e.g. Efroimsky 2012) forced both by the tidal potential and the loading induced by the tidal distortion of the ocean. This solidfluid coupling affects the frequencies given by Eq. (77) and generates an additional coupling between tidal modes (e.g. Matsuyama 2014). We simply ignored this coupling so that we did not need to consider the parameters of the solid response.
Fig. 6 Imaginary part of the tidal Love number associated with the quadrupolar oceanic tide. The logarithm of is plotted asa function of the orbital (horizontal axis) and rotation (vertical axis) periods in logarithmic scale by using Eq. (70) for σ_{R} = 10^{−6} s^{−1} (bottom) and σ_{R} = 10^{−5} s^{−1} (top), and various values of H and N. Horizontally, H = 10 km (left), H = 100 km (middle) and H = 1000 km (right). Vertically, N = 10^{−4} s^{−1} (bottom) and N = 10^{−5} s^{−1} (top). Colours correspond to logarithmic decades (colourbars on the right). The diagonal black dashed line corresponds to spinorbit synchronization. 
Fig. 7 Tidal torque associated with the Lunar quadrupolar oceanic tide. The logarithm of is plotted asa function of the orbital (horizontal axis) and rotation (vertical axis) periods in logarithmic scale by using Eq. (72) for σ_{R} = 10^{−6} s^{−1} (bottom) and σ_{R} = 10^{−5} s^{−1} (top), and various values of H and N. Horizontally, H = 10 km (left), H = 100 km (middle) and H = 1000 km (right). Vertically, N = 10^{−4} s^{−1} (bottom) and N = 10^{−3} s^{−1} (top). Colours correspond to logarithmic decades (colourbars on the right). The diagonal black dashed line corresponds to spinorbit synchronization. 
7 Conclusion
To better understand the tidal response of the deep global oceans potentially hosted by recently discovered extrasolar planets, we developed a linear model reducing the physics of tides to the essentials. Our approach evolves from the traditional 2D modelling inherited from the study of the Earth’s oceans to include 3D effects resulting from the properties of the oceanic vertical structure. Particularly, it introduces the contribution of internal gravity waves induced by the stable stratification and takes into account dissipative mechanisms with a Rayleigh drag, following the early work by Webb (1980). Hence, we wrote the equations describing the tidal response of a deep ocean in the general case, as well as the associated tidal torque, Love numbers, and quality factor. We then simplified these equations in the framework of the thinlayer approximation to compute an analytic solution. This solution was finally used to explore the domain of parameters. We first treated the cases of idealized Earth and TRAPPIST1 f planets to constrain the Rayleigh coefficient (σ_{R}) and illustrate asymptotic regimes. In a second phase, we examined in a systematic way the dependence of the tidal Love number and torque associated with the semidiurnal tide on the ocean depth (H) and Brunt–Väisälä frequency (N).
In the thin layer limit, we recover the behaviour identified by early studies (e.g. LonguetHiggins & Pond 1970; Webb 1980). The tidal response is composed of resonant surface inertialgravity modes due to the conjugated effects of gravity and rotation. In this case, the role played by stratification is negligible. The fluid compressibility can affect the tidal response in the highfrequency range, where horizontally propagating Lamb modes can be excited. The importance of the role played by the stratification grows with the ocean depth. The stable stratification allows internal gravity waves to propagate. It induces internal density fluctuations characterized by a frequencyresonant behaviour, similar to surface gravity waves. The contribution of this component is not negligible inthe case of deep oceans and sensitively increases the evolution timescales of the planetperturber system.
Although the linear analysis developed in this work is simplified with respect to numerical models, it provides a very convenient way to explore the broad domain of planetary parameters and unravel the complex dynamics of oceanic tides with a reasonable computational cost and few physical control parameters. These analytical results can be used in the future as benchmark validations for fully 3D GCM tidal studies. A next step should be to introducemeridian boundaries in the model in order to study the effect of blocking the progression of the tidal response and thus better characterize the role played by continental shelves in the oceanic tidal dissipation.
Acknowledgements
The authors acknowledge funding by the European Research Council through ERC grants SPIRE 647383 and WHIPLASH 679030, the Programme National de Planétologie (INSU/CNRS) and the CNRS PLATO grants at CEA/IRFU/DAp. They wish to thank the referee, Robert Tyler, for his helpful suggestions and remarks.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions, 1046 [Google Scholar]
 AngladaEscudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 AuclairDesrotour, P., Le PoncinLafitte, C., & Mathis, S. 2014, A&A, 561, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Auclair Desrotour, P., Mathis, S., & Le PoncinLafitte C. 2015, A&A, 581, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Laskar, J., & Mathis, S. 2017a, A&A, 603, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Laskar, J., Mathis, S., & Correia, A. C. M. 2017b, A&A, 603, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Mathis, S., & Laskar, J. 2018, A&A, 609, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bills, B. G., & Ray, R. D. 1999, Geophys. Res. Lett., 26, 3045 [NASA ADS] [CrossRef] [Google Scholar]
 Booker, J. R., & Bretherton, F. P. 1967, J. Fluid Mech., 27, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Chapman, S., & Lindzen, R. 1970, Atmospheric Tides. Thermal and Gravitational, 200 [Google Scholar]
 Chen, E. M. A., Nimmo, F., & Glatzmaier, G. A. 2014, Icarus, 229, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Levrard, B., & Laskar, J. 2008, A&A, 488, L63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cowling, T. G. 1941, MNRAS, 101, 367 [NASA ADS] [Google Scholar]
 Dickey, J. O., Bender, P. L., Faller, J. E., et al. 1994, Science, 265, 482 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Doodson, A. 1938, Phil. Trans. R. Soc. London, Ser. A, Mathematical and Physical Sciences, 237, 311 [NASA ADS] [CrossRef] [Google Scholar]
 Eakins, B. W., & Sharman, G. F. 2010, NOAA National Geophysical Data Center, Boulder, CO [Google Scholar]
 Eckart, C. 1960, Phys. Fluids, 3, 421 [Google Scholar]
 Efroimsky, M. 2012, ApJ, 746, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Efroimsky, M., & Williams, J. G. 2009, Celest. Mech. Dyn. Astron, 104, 257 [Google Scholar]
 Egbert, G. D., & Ray, R. D. 2001, J. Geophy. Res., 106, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Egbert, G. D., & Ray, R. D. 2003, Geophys. Res. Lett., 30, 1907 [NASA ADS] [CrossRef] [Google Scholar]
 Favier, B., Barker, A. J., Baruteau, C., & Ogilvie, G. I. 2014, MNRAS, 439, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J., & Lai, D. 2014, MNRAS, 444, 3488 [NASA ADS] [CrossRef] [Google Scholar]
 Garrett, C., & Kunze, E. 2007, Ann. Rev. Fluid Mech., 39, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Garrett, C., & Munk, W. 1979, Ann. Rev. Fluid Mech., 11, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Gerkema, T., & Shrira, V. I. 2005, J. Fluid Mech., 529, 195 [NASA ADS] [CrossRef] [Google Scholar]
 Gerkema, T., & Zimmerman, J. 2008, Lecture Notes, Royal NIOZ, Texel [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Guenel, M., Baruteau, C., Mathis, S., & Rieutord, M. 2016, A&A, 589, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hendershott, M. C. 1981, Evolution of physical oceanography, 292 [Google Scholar]
 Hough, S. S. 1898, Roy. Soc. London Philos. Trans. Ser. A, 191, 139 [Google Scholar]
 Kaula, W. M. 1962, AJ, 67, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Kaula, W. M. 1966, Theory of satellite geodesy. Applications of satellites to geodesy, 140 [Google Scholar]
 Khurana, K. K., Kivelson, M. G., Stevenson, D. J., et al. 1998, Nature, 395, 777 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kivelson, M. G., Khurana, K. K., & Volwerk, M. 2002, Icarus, 157, 507 [NASA ADS] [CrossRef] [Google Scholar]
 Lambeck, K. 1977, Phil. Trans. R. Soc. London, Ser. A, 287, 545 [Google Scholar]
 Lambeck, K. 1980, The earth’s variable rotation: Geophysical causes and consequences, 458 [Google Scholar]
 Laplace, P. S. 1798, Traité de mécanique céleste (Duprat J. B. M.) [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 [NASA ADS] [CrossRef] [Google Scholar]
 LonguetHiggins, M. S., & Pond, G. S. 1970, Phil. Trans. R. Soc. London, Ser. A, 266, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Makarov, V. V. 2012, ApJ, 752, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, S. 2009, A&A, 506, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Le PoncinLafitte C. 2009, A&A, 497, 889 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Talon, S., Pantillon, F.P., & Zahn, J.P. 2008, Sol. Phys., 251, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuyama, I. 2014, Icarus, 242, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
 Ogilvie, G. I. 2009, MNRAS, 396, 794 [NASA ADS] [CrossRef] [Google Scholar]
 Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
 Ogilvie, G. I., & Lin, D. N. C. 2004, ApJ, 610, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Prat, V., Mathis, S., Lignières, F., Ballot, J., & Culpin, P.M. 2017, A&A, 598, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Proudman, J.,& Doodson, A. 1936, Phil. Trans. R. Soc. London, Ser. A, Mathematical and Physical Sciences, 235, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Ray, R. D., Eanes, R. J., & Lemoine, F. G. 2001, Geophys. J. Int., 144, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Remus, F., Mathis, S., Zahn, J.P., & Lainey, V. 2012, A&A, 541, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ribas, I., Bolmont, E., Selsis, F., et al. 2016, A&A, 596, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Tort, M., & Dubos, T. 2014, Q. J. Royal Meteorol. Soc., 140, 2388 [Google Scholar]
 Tyler, R. 2011, Icarus, 211, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Tyler, R. 2014, Icarus, 243, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars, 420 [Google Scholar]
 Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics, 770 [Google Scholar]
 Volland, H. 1974a, J. Atm. Terr. Phys, 36, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Volland, H. 1974b, J. Atm. Terr. Phys, 36, 1975 [NASA ADS] [CrossRef] [Google Scholar]
 Volland, H., & Mayr, H. G. 1972, J. Atm. Terr. Phys, 34, 1769 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, S., Wu,D.H., Barclay, T., & Laughlin, G. P. 2017, ApJ, submitted, [arXiv:1704.04290] [Google Scholar]
 Webb, D. J. 1980, Geophys. J., 61, 573 [NASA ADS] [CrossRef] [Google Scholar]
The represent here the normalized associated Legendre polynomials, expressed as (Abramowitz & Stegun 1972).
All Tables
All Figures
Fig. 1 Quadrupolar semidiurnal oceanic tide in a rotating terrestrial planet gravitationally forced by a perturber (star of satellite). Reference frame and system of coordinates used in the modelling. The spin frequency of the planet and the orbital frequency of the perturber are denoted Ω and n_{orb}, respectively. The notation δ_{tide} designates the angular lag of the tidal bulge. 

In the text 
Fig. 2 Real (top) and imaginary (bottom) parts of even Hough functions (n even) associated with a quadrupolar tidal perturbation (m = 2) and defined by with ν =1 (ordinary/gravity modes) and various positive values of γ = σ_{R}∕σ. From left to right, γ = 0 (nonfrictional regime), γ = 0.1 (weakly frictional regime), γ = 1 and γ = 10 (stronglyfrictional regime). Hough functions are plotted as functions of the latitude δ (degrees). 

In the text 
Fig. 3 Frequency spectrum of waves likely to be excited by the tidal perturbation and of the corresponding tidal regimes. The characteristic frequencies of the system are, from left to right, the drag frequency σ_{R}, the inertia frequency 2Ω, the surfacegravity waves cutoff frequency , the Brunt–Väisälä frequency N, and the acoustic cutoff frequency σ_{s} = c_{s}∕r. Their hierarchy can vary as a function of the physical and dynamical properties of the planet. In the case of the Earth’s ocean, σ_{R} ≈ 10^{−5} s^{−1} (Webb 1980), 2Ω = 1.46 × 10^{−4} s^{−1}, σ_{g} = 3.1 × 10^{−5} s^{−1}, N ≈ 2.2 × 10^{−3} s^{−1} (Gerkema & Zimmerman 2008), and σ_{s} = 6.3 × 10^{−3} s^{−1}. 

In the text 
Fig. 4 Vertical displacement due to the quadrupolar semidiurnal tide (, n_{orb} being the dynamical frequency) obtained using the analytic solution established for the thin layer approximation (Eq. (58)). The displacement, normalized by the quadrupolar tidal potential, is plotted in the equatorial plane of the planet as a function of longitude (horizontal axis) and normalized altitude (vertical axis). The position φ = 0 corresponds to the subperturber point; the oceanic floor and surface are located at x = 0 and x = 1, respectively. For this computation, the rotation period of the planet P_{spin} = 2π∕Ω and orbital period of the perturber P_{orb} = 2π∕n_{orb} are set to P_{spin} = 6.0 d and P_{orb} = 365.25 d. The used values of parameters are R = R_{⊕}, H = 100 km, g = 9.81 m s^{−2}, N = 10^{−3} s^{−1}, c_{s} = 1545 m s^{−1} and σ_{R}= 10^{−5} s^{−1}. 

In the text 
Fig. 5 Semidiurnal oceanic tide of the Earth (left column) and TRAPPIST1 f (right column). The logarithms of the imaginary part of the tidal Love number (left) and tidal quality factor (right) are plotted as functions of the normalized forcing frequency (where Ω_{⊕} designates the today rotation rate of the Earth) for various orders of magnitude of the drag parameter . These frequency spectra are computed using the expressions given by Eqs. (70) and (71). In each case, the parameter n_{orb} is assumed to be constant and the rotation rate of the planet varies with the tidal frequency following the formula Ω = n_{orb} + σ∕2. The resonances associated with surface inertialgravity modes are designated by black dashed lines in the top panels and numbers indicate the degree n of the corresponding Hough modes and the sign of eigenfrequencies given by Eq. (77). The values of parameters used for this evaluation are summarized in Table 1. 

In the text 
Fig. 6 Imaginary part of the tidal Love number associated with the quadrupolar oceanic tide. The logarithm of is plotted asa function of the orbital (horizontal axis) and rotation (vertical axis) periods in logarithmic scale by using Eq. (70) for σ_{R} = 10^{−6} s^{−1} (bottom) and σ_{R} = 10^{−5} s^{−1} (top), and various values of H and N. Horizontally, H = 10 km (left), H = 100 km (middle) and H = 1000 km (right). Vertically, N = 10^{−4} s^{−1} (bottom) and N = 10^{−5} s^{−1} (top). Colours correspond to logarithmic decades (colourbars on the right). The diagonal black dashed line corresponds to spinorbit synchronization. 

In the text 
Fig. 7 Tidal torque associated with the Lunar quadrupolar oceanic tide. The logarithm of is plotted asa function of the orbital (horizontal axis) and rotation (vertical axis) periods in logarithmic scale by using Eq. (72) for σ_{R} = 10^{−6} s^{−1} (bottom) and σ_{R} = 10^{−5} s^{−1} (top), and various values of H and N. Horizontally, H = 10 km (left), H = 100 km (middle) and H = 1000 km (right). Vertically, N = 10^{−4} s^{−1} (bottom) and N = 10^{−3} s^{−1} (top). Colours correspond to logarithmic decades (colourbars on the right). The diagonal black dashed line corresponds to spinorbit synchronization. 

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.