Oceanic tides from Earth-like to ocean planets
Pierre Auclair-Desrotour Stéphane Mathis Jacques Laskar Jérémy Leconte
astro-ph.EP 85-02
 
  
 
  

  
  IMCCE, Observatoire de Paris, CNRS UMR 8028, PSL, 77 Avenue Denfert-Rochereau, 75014 Paris, France

pierre.auclair-desrotour@u-bordeaux.fr, jacques.laskar@obspm.fr
Laboratoire AIM Paris-Saclay, CEA/DRF - CNRS - Université Paris Diderot, IRFU/DAp Centre de Saclay, F-91191 Gif-sur-Yvette, France
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, F-92195 Meudon, France

stephane.mathis@cea.fr
Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France 

jeremy.leconte@u-bordeaux.fr
    

  


  
  
   Oceanic tides are a major source of tidal dissipation. They drive the evolution of planetary systems and the rotational dynamics of planets. However, 2D models commonly used for the Earth cannot be applied to extrasolar telluric planets hosting potentially deep oceans because they ignore the three-dimensional effects related to the ocean vertical structure.
  
   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 modeling allowing to treat a wide range of cases from shallow to deep oceans. 
  
   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.
  
   In the 2D asymptotic limit, we recover the frequency-resonant behaviour due to surface inertial-gravity 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 become three-dimensional. 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 affect significantly the evolution time scale of the planetary rotation. 
  
   

  

Oceanic tides from Earth-like to ocean planets
    P. Auclair-Desrotour1,2,4
     S. Mathis2,3J. Laskar1
     J. Leconte4
          
    Received ...; accepted ...
=====================================================================================







§ 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 exoplanets. First, in Summer, 2016, a telluric extrasolar planet was found orbiting the Sun's closest stellar neighbour, the red dwarf Proxima Centauri <cit.>. This planet, Proxima b, has a minimum mass of 1.3 M_ and its equilibrium temperature is within the range where water could be liquid on its surface, which makes it appear as a potential twin sister of the Earth. A few months later, another star was in the spotlight, the ultra-cool dwarf star TRAPPIST-1 <cit.>. The system hosted by TRAPPIST-1 drew attention for 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_ <cit.>. 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 <cit.>. Hence, Proxima b and several of the telluric planets orbiting TRAPPIST-1 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 proved by the case of the 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_ <cit.>. Yet, the Earth's ocean is today the main contributor to the total energy tidally dissipated by the Earth-Moon-Sun system. It accounts for roughly 95 % of the planetary dissipation rate of 2.54 TW generated by the Lunar semidiurnal tide <cit.>. Moreover, the oceanic dissipation rate can vary significantly owing to the frequency-resonant behaviour of the fluid layer <cit.>, thus leading to the present-day high transfer of angular momentum from the spin of the Earth to the orbit of the Moon <cit.>. More generally, oceanic tides have an impact on the history of the spin rotation and states of equilibrium of the Earth <cit.> and terrestrial exoplanets <cit.>. 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 <cit.>, 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 (or GCM) <cit.>. Moreover, numerous studies have characterized their dependence on the planet parameters by developing linear ab initio global models, such as the hemispherical ocean model <cit.>. This approach has been used lately to estimate the tidal dissipation rate in the oceans of icy satellites orbiting Jupiter and Saturn <cit.>. 

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 computational cost relatively small compared to numerical models taking into account local features (e.g. topography, mean flows). They are commonly based on a two dimensional 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 surface gravity 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 it has been shown for long <cit.>, the tidal response of a stratified fluid layer in the classical tidal theory is characterized by a frequency-dependent equivalent depth, which is the analogous of the barotropic ocean depth <cit.>. 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 analyzed to our knowledge.



Therefore, following along the line of <cit.> and <cit.>, we develop here a linear global model of oceanic tides including three-dimensional 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 planet from the architecture of their host planetary system, and (3) provide a simplified diagnosis of the oceanic tidal dissipation rate.

In Section <ref>, 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 Section <ref>, 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 TRAPPIST-1 f planets in Section <ref> in order to illustrate the difference between shallow and deep oceans. We use it in Section <ref> to explore the dependence of the dissipation rate on the ocean depth and stratification. Finally, we discuss the simplifications assumed in the modeling in Section <ref> and give our conclusions in Section <ref>. 





§ DYNAMICS OF A THICK OCEAN FORCED GRAVITATIONALLY



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. <ref>). 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, ℛ_ E;T : { O, X_ E, Y_ E, Z_ E}, 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 ( e_r , e_θ, e_φ) and coordinates ( r, θ , φ) 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.  







 §.§ 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 Ω_ c = √(g / r) 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 


    [ p ( r , t ) = p_0 ( r ) + δ p ( r , t ),     ρ( r ,t ) = ρ_0 ( r ) + δρ( r , t ), ]
 

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 V( r , t ) = ( V_r , V_θ , V_φ) and the associated displacement ξ <cit.>. Note that, because of the solid rotation assumption, there is no meridional or zonal background circulation. We will 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 <cit.>, i.e. the self gravitational effect of the variation of mass distribution is not taken into account. Finally, following early works <cit.>, we introduce dissipation by using a Rayleigh friction, σ_ RV, where σ_ R is a constant effective frequency associated with the damping <cit.>. As demonstrated by <cit.>, 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 solid-fluid coupling is naturally unknown.

Under these assumptions, the Navier-Stokes equation describing the distortion of the ocean forced by the tidal potential U can be expressed in the frame rotating with the planet as <cit.>


    ∂_t V + 2 Ω∧V = - 1/ρ_0∇δ p - g/ρ_0δρe_r + ∇ U - σ_ RV.


It is completed by the equation of mass conservation


    ∂_t δρ + d ρ_0/dr V_r + ρ_0 ∇·V = 0,


and the equation of buoyancy


    ∂_t δρ + d ρ_0/dr V_r = 1/c_s^2( ∂_t δ p  + d p_0/dr V_r ),


where c_s ( r ) = Γ( p_0 , ρ_0 , S_0 ) 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 <cit.>,


    N^2 = - g/H[ d lnρ_0/dx + g H/c_s^2],
 

where we have introduced the non-dimensional reduced altitude x, defined by r = R_ c + H x. The three components of the momentum equation given by Eq. <ref> are strongly coupled by the Coriolis acceleration. To make the problem tractable analytically, we assume the commonly used traditional approximation <cit.>. 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


    [                        ∂_t V_θ - 2 Ωcosθ V_φ                                            =         - 1/r∂_θ( δ p/ρ_0 - U )  - σ_ R V_θ,;                        ∂_t V_φ + 2 Ωcosθ V_θ                                            =     - 1/r sinθ∂_φ( δ p/ρ_0 - U ) - σ_ R V_φ,;                                      ∂_t V_r                                            = - 1/ρ_0∂_r δ p - g/ρ_0δρ + ∂_r U - σ_ R V_r. ]


It thus makes possible the separation of the θ and r coordinates in solutions, as showed thereafter.  As discussed in Section <ref>, 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 super-inertial waves <cit.>. 

Finally, the equations of mass conservation (<ref>) and of the buoyancy (<ref>) can be written respectively


    ∂_t δρ + 1/r^2∂_r ( r^2 ρ_0 V_r ) = - ρ_0/r sinθ[ ∂_θ( sinθ V_θ) + ∂_φ V_φ],


and


    ∂_t δρ = 1/c_s^2∂_t δ p + ρ_0 N^2/g V_r.


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


    f = ∑_m,σ f^m,σ( r , θ) ^i ( σ t + m φ),



where the f^m,σ are the spatial distributions of f associated with the doublet ( m , σ). 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 so-called spin parameter <cit.>. Viscous friction modifies slightly the latitudinal structure by introducing a dependence on σ_ R. Thus, we define the complex tidal frequency and spin parameter as


    [ σ̃ = σ - i σ_ R                    ν̃ = 2 Ω/σ̃, ]


and the latitudinal operators 


    [           ℒ_θ^m,ν̃ = 1/1 - ν̃^2 cos^2 θ[ ∂_θ + m ν̃θ]                                                      ; ℒ_φ^m,ν̃ =  1/1 - ν̃^2 cos^2 θ[ ν̃cosθ∂_θ + m/sinθ] ,                                                       ]


such that


    [  V_θ^m,σ = i/σ̃ r ℒ_θ^m,ν̃( δ p^m,σ/ρ_0 - U^m,σ)                                                 ; V_φ^m,σ = - 1/σ̃rℒ_φ^m,ν̃( δ p^m,σ/ρ_0 - U^m,σ). ]


 

As mentioned above, the traditional approximation allows us to separate the coordinates r and θ in the Fourier coefficients of Eq. (<ref>). Therefore, we write these functions as


    [     ξ_r^m,σ = ∑_n ξ_r;n^m,σ( r ) Θ_n^m,ν̃( θ),     V_r^m,σ = ∑_n V_r;n^m,σ( r ) Θ_n^m,ν̃( θ),; ξ_θ^m,σ = ∑_n ξ_θ;n^m,σ( r ) Θ_θ ; n^m,ν̃( θ), V_θ^m,σ = ∑_n V_θ;n^m,σ( r ) Θ_θ ; n^m,ν̃( θ),; ξ_φ^m,σ = ∑_n ξ_φ;n^m,σ( r ) Θ_φ ; n^m,ν̃( θ), V_φ^m,σ = ∑_n V_φ;n^m,σ( r ) Θ_φ ; n^m,ν̃( θ),;     δ p^m,σ = ∑_n δ p_n^m,σ( r ) Θ_n^m,ν̃( θ),       δρ^m,σ = ∑_n δρ_n^m,σ( r ) Θ_n^m,ν̃( θ),;         U^m,σ = ∑_n U_n^m,σ( r ) Θ_n^m,ν̃( θ),                                                ]


where n designates the latitudinal wavenumber of a component (n ∈ℕ if | ν| ≤ 1 ; n ∈ℤ if | ν| > 1), the Θ_n^m,ν̃ are the so-called Hough functions <cit.>, defined on the interval θ∈[ 0 , π], and the Θ_θ ; n^m,ν̃( θ) = ℒ_θ^m,ν̃[ Θ_n^m,ν̃( θ) ] and Θ_φ ; n^m,ν̃( θ) = ℒ_φ^m,ν̃[ Θ_n^m,ν̃( θ) ] stand for the latitudinal functions associated with the latitudinal and longitudinal velocities and displacements respectively. Hough functions are the solutions of the eigenfunctions-eigenvalues problem defined by the Laplace's tidal equation <cit.>,


    ℒ^m,ν̃Θ = - ΛΘ,


where Θ is a function, Λ∈ℂ and ℒ^m,ν̃ the operator


    [                                                                   ℒ^m,ν̃ =                                   1/sinθd/d θ( sinθ/1 - ν̃^2 cos^2 θd/d θ);                                                                            - 1/1 - ν̃^2 cos^2 θ( m ν̃1 + ν̃^2 cos^2 θ/1- ν̃^2 cos^2 θ + m^2/sin^2 θ). ]




Hence, for a given n, Θ_n^m,ν̃ is such that ℒ^m,ν̃Θ_n^m,ν̃ = -  Λ_n^m,ν̃Θ_n^m,ν̃, the parameter Λ_n^m,ν̃ being the associated eigenvalue. We can note here that the only difference with the usual case where friction is not taken into account <cit.> is the nature of ν̃ which is a complex number and not a real one (see Eq. (<ref>)). It follows that the Θ_n^m,ν̃ and Λ_n^m,ν̃ are complex functions and parameters in the general case <cit.>. 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 Θ_n^m,ν and Λ_n^m,ν). This regime corresponds to the case treated usually in the theory of tides <cit.>. Tidal modes can be divided into two families: the so-called gravity modes (n ≥ 0), which are defined both in the regime of super-inertial waves (| ν| ≤ 1) and in the regime of sub-inertial waves (| ν| > 1), and the inertial modes (n < 0), defined only in the regime of sub-inertial waves <cit.>. The Hough functions associated with gravity modes degenerate in the associated Legendre polynomials[The P_l^m represent here the normalized associated Legendre polynomials, expressed as P_l^m   = [ ( - 1 )^m / ( 2^l l ! ) ] ( 1 - x^2 )^m/2( d^l+m / dx^l+m) ( x^2 - 1 )^l <cit.>.] P_l^m (with l = m + n) when ν→ 0.
   
  ∙ a strongly frictional regime (| σ| ≲σ_ R), where friction affects significantly the horizontal structure of tidal waves. In this regime, the hierarchy of eigenvalues can change, as demonstrated by <cit.> who notes that the corresponding critical points appears 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 <cit.> and <cit.>.
 
 
 Hough functions associated with symmetric modes of degree m = 2 are plotted on Fig. <ref> for various values of the ratio σ_ R / σ. The equations describing the vertical structure are obtained by substituting the expansions given by Eqs. (<ref>) and (<ref>) in Eqs. (<ref>-<ref>). In order to lighten expressions, the superscripts ( m , σ) will be omitted up to the end of Section <ref>, 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 = - H g ρ_0), where x is the reduced altitude introduced in Eq. (<ref>). After some manipulations, one gets
 

    [ d Y/dx =
    [  A_1 ( x )  B_1 ( x );  A_2 ( x ) B_2  ( x ) ]Y +
    [  C_1( x ); C_2 ( x ) ],                                                          Y = 
    [       y_n; r^2 ξ_r;n ] ]


and the coefficients


    [                       A_1 = N^2 H/g, A_2 = H ( Λ_n^m,ν̃/σσ̃ - r^2/c_s^2),;            B_1 = H/r^2( σσ̃ - N^2 ),                     B_2 = g H/c_s^2,;                     C_1 =  d U_n/dx,          C_2 = - H Λ_n^m,ν̃/σσ̃ U_n. ]



This system is then written as a single second-order ordinary differential equation,


    d^2 y_n/dx^2 + A ( x ) d y_n/dx  + B ( x ) y_n = C ( x ),


with the coefficients











    A ( x ) =     d lnρ_0/dx - K_∘ , 
     B ( x )  =     H^2/r^2Λ_n^m,ν̃( N^2/σσ̃ - 1 ) ( 1 - ε_s ; n)   
         - ( d/dx - gH/c_s^2 - K_∘) ( N^2 H/g), 
     C ( x )  =     H^2/r^2Λ_n^m,ν̃( N^2/σσ̃ - 1 ) U_n - ( g H/c_s^2 + K_∘) d U_n/dx + d^2 U_n/dx^2 .





Here, we have introduced the frequency-dependent, dimensionless acoustic and sphericity parameters


    [          ε_s;n = r^2 σσ̃/Λ_n^m,ν̃ c_s^2                                         K_∘ = r^2/N^2 - σσ̃d/dx(N^2 - σσ̃/r^2). ]


 

The acoustic parameter can be written ε_s;n = σσ̃ / σ_s;n^2, where σ_s;n = ( Λ_n^m,ν̃)^1/2 c_s /r is the cutoff frequency of acoustic waves (also called the Lamb frequency) for the degree-n mode when Λ_n^m,ν̃∈ℝ^+. The regime of acoustic waves thus corresponds to | ε_s;n| ≳ 1. Typically, for an Earth-like planet of radius R_ = 6371 km and angular velocity Ω_ = 7.25 × 10^-5 s^-1 (https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.htmlNASA fact sheets), with a water oceanic layer characterized by c_s ≈ 1.545   km.s^-1 <cit.>, the gravest acoustic modes of the quadrupolar semidiurnal tide (m =2, n = 0, Λ_0^2,ν̃∼ 11) 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 two-dimensional modelings where r = R. We will see in the next section that K_∘ ∝  H/R when H/R ≪ 1. Eventually, using the change of variables y_n = Φ_n^m,σΨ_n^m,σ, where Φ_n^m,σ is the radial function 




    Φ_n^m,σ( x ) =  exp[ - 1/2∫_0^x A ( x' )  x' ],


we write the vertical structure equation (Eq. <ref>) as the Schrödinger-like equation


    d^2 Ψ_n^m,σ/dr^2 + k̂_n^2 Ψ_n^m,σ =  ( Φ_n^m,σ)^-1 C ,


in which k̂_n represents the vertical wavenumber of the mode, defined as


    k̂_n^2 =  B - 1/2( d A/dx + A^2/2).


The vertical profiles of the other quantities are deduced from Ψ_n^m,σ straightforwardly thanks to polarization relations. We get, for the velocity field


    V_r ; n^m,σ = - i σ/H ( N^2 - σσ̃)[ Φ_n ( d Ψ_n/dx + 𝒜_n Ψ_n ) - d U_n/dx],



    V_θ ; n^m,σ = i/σ̃ r( Φ_n Ψ_n - U_n ),



    V_φ ; n^m,σ = - 1/σ̃ r( Φ_n Ψ_n - U_n ) ,


for the displacement


    ξ_r ; n^m,σ = - 1/H ( N^2 - σσ̃)[ Φ_n ( d Ψ_n/dx + 𝒜_n Ψ_n ) - d U_n/dx],



    ξ_θ ; n^m,σ = 1/σσ̃ r( Φ_n Ψ_n - U_n ),



    ξ_φ ; n^m,σ = i/σσ̃ r( Φ_n Ψ_n - U_n ),


and for scalar quantities


    δ p_n^m,σ =  ρ_0 Φ_n Ψ_n,



    δρ_n^m,σ = - ρ_0/gHN^2/N^2 - σσ̃[ Φ_n ( d Ψ_n/dx + ℬ_n Ψ_n ) - d U_n/dx],


with the frequency-dependent factors


    𝒜_n ( x )  =  1/2[  gH/c_s^2 - N^2 H/g + K_∘] , 
       ℬ_n ( x )  = 1/2[ Hg/c_s^2( 2 σσ̃/N^2 - 1 ) - H N^2/g + K_∘].


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 surface gravity waves, which are restored by gravity and are characterized by the cutoff frequency σ_ g = √(gH) / r. If the ocean is stably-stratified (N > 0), internal gravity waves can propagate in the frequency range | σ| < N. These waves are restored by the Archimedean force. Finally, in the high-frequency 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. <ref>. 




 §.§ Tidal potential, Love numbers and tidal torque


The variation of mass distribution due to the tidal distortion modifies the self-gravitational potential of the planet. The tidal potential of the excitation is usually expanded in Fourier series and spherical harmonics <cit.>. It is thus convenient to write the fluctuations of the self-gravitational potential, denoted 𝒰, in the same form,


    𝒰 = ∑_σ,l,m𝒰_l^m,σ( x ) P_l^m ( cosθ) ^i ( σ t + m φ).


The radial profiles of the ( m , σ)-components are themselves written
 

    𝒰_l^m,σ = 𝒰_ξ ; l^m,σ + 𝒰_ρ ; l^m,σ,
 

where we have introduced the contribution of the surface displacement 𝒰_ξ ; l^m,σ and of the internal density variations 𝒰_ρ ; l^m,σ. In the following, we will 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 


    𝒰_ξ ; l^m,σ( 1 ) =  4 π𝒢 R/2 l + 1ρ_0 ( 1 )  ξ_r ; l^m,σ( 1 ), 
       𝒰_ρ ; l^m,σ( 1 )  =  4 π𝒢 R/2 l + 1 H [ ∫_0^1 ( 1 + H/R x' )^l+2δρ_l^m,σ x' ],



 
with 𝒢 representing the universal gravity constant. Expanded in Hough functions, they write


    𝒰_ξ ; l^m,σ( 1  ) = 4 π𝒢 R/2 l + 1ρ_0 ( 1 ) ∑_n ∑_k ≥  m  C_l,n,k^m,ν̃ξ_r ; n,k^m,σ( 1 ) ,



    𝒰_ρ ; n^m,σ( 1 ) = 4 π𝒢 R/2 l + 1 H ∑_n ∑_k ≥ m C_l,n,k^m,ν̃∫_0^1 ( 1 + H/R x' )^l+2δρ_n,k^m,σ x'  .












In these expressions, the complex weighting coefficients C_l,n,k^m,ν̃ are expressed as 


    C_l,n,k^m,ν̃ = B_k,n^m,ν̃ A_n,l^m,ν̃,






where the A_n,l^m,ν̃ and B_k,n^m,ν̃ designate the mutual projection coefficients of the Hough functions on normalized associated Legendre polynomials, defined respectively by the expansions


    Θ_n^m,ν̃( θ) = ∑_l ≥ m  A_n,l^m,ν̃ P_l^m ( cosθ),



    P_k^m ( cosθ) = ∑_n B_k,n^m,ν̃Θ_n^m,ν̃( θ).


The C_l,n,k^m,ν̃ coefficients quantify the coupling induced by the Coriolis effects in the dynamics of the tidal response. They are real in the non-frictional case (σ_ R = 0 and ν̃ = ν), where A_n,l^m,ν = ⟨ P_l^m , Θ_n^m,ν⟩ and B_k,n^m,ν = ⟨Θ_n^m,ν , P_k^m ⟩ = A_n,k^m,ν, the notation ⟨· , ·⟩ standing for the scalar product defined for any P_l^m and Θ_n^m,ν as


    ⟨ P_l^m , Θ_n^m,ν⟩ = ∫_0^π P_l^m ( cosθ) Θ_n^m,ν( θ) sinθ θ.


In the case of a non-rotating planet (Ω = 0 and ν = 0), there is no coupling. Therefore, C_l,n,k^m,0 =  1 if l = n + m = k and C_l,n,k^m,0 =  0 else. The notations ξ_r ; n,k^m,σ = ξ_r ; n^m,σ( U_k^m,σ) and δρ_n,k^m,σ = δρ_n^m,σ( U_k^m,σ) in Eqs. (<ref>) and (<ref>) stand for the vertical displacement and density variations due to the P_k^m-component of the forcing projected on the ( n , m, ν̃)-Hough function. The expressions of the components of the tidal potential (Eqs. <ref> and <ref>) allow us to compute complex Love numbers, which are commonly used to quantify tidal dissipation in celestial bodies <cit.>. Love numbers are defined as the ratio between the gravitational potential due to the distortion and the excitating tidal gravitational potential taken at the external surface of the body (r = R). Thus, oceanic Love numbers associated with the triplet ( l , m , σ), denoted k_l^m,σ, are defined as  


    k_l^m,σ = . 𝒰_l^m,σ/U_l^m,σ|_x=1 .
 

Here, { k_l^m,σ} and  { k_l^m,σ} stand for the tidal response of the ocean to the ( l, m )-component of the excitating tidal potential. The U_l^m are provided by the Kaula's multipolar expansion of the tidal gravitational potential and are expressed as functions of the Keplerian elements of the planet-perturber system <cit.>. 

Because of 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


    𝒯 = ∑_m,σ( 𝒯_ξ^m,σ + 𝒯_ρ^m,σ).



The components 𝒯_ξ^m,σ and 𝒯_ρ^m,σ of the ( m , σ)-mode are defined as 


    𝒯_ξ^m,σ = {1/2ρ_ s∫_∂𝒱_0. ∂ U^m,σ/∂φ|_r=R[  ξ_r^m,σ( 1 ) ]^*  S  }, 
       𝒯_ρ^m,σ = {1/2∫_𝒱_0∂ U^m,σ/∂φ( δρ^m,σ)^* 𝒱},


where ^* stands for the conjugate of a complex number and  its real part ( will be used for the imaginary part), ρ_ s = ρ_0 ( 1 ) denotes the density at the surface ocean,  𝒱_0 the spatial domain filled by the oceanic shell at rest, ∂𝒱_0 its upper boundary, and S and 𝒱 infinitesimal surface and volume parcels respectively. Similarly as the tidal potential due to the distortion, these two components can be expanded on Hough functions and associated Legendre polynomials. One thus obtains, for l ≥ m and k ≥ m,


    𝒯_ξ^m,σ = - m π R^2 ρ_ s∑_l,n,k{( C_l,n,k^m,ν̃)^* U_l^σ,m( 1 ) [ ξ_r ; n,k^m,σ( 1 )  ]^*  },



    𝒯_ρ^m,σ = - m π R^2 H∑_l,n,k{( C_l,n,k^m,ν̃)^* ∫_0^1 U_l^m,σ[ δρ_n,k^m,σ]^*  x  }.





§ TIDES IN A THIN STABLY-STRATIFIED OCEAN



In some of observed cases, the oceanic layer of terrestrial planets and satellites is thin compared to the radius of the body <cit.>. Therefore, in this section, we apply our general modeling 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.



 §.§ Tidal waves dynamics



The thin shell hypothesis allows us to simplify the physical setup. Given that the relative difference of gravity between the lower and the upper boundary 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^2U_n /dx^2 ≈ 0). Furthermore, we assume for simplicity uniform stratification and compressibility, i.e. 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 <cit.>. 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. (<ref>), setting N^2 to a constant implies a density decreasing exponentially with altitude,


    ρ_0 ( x ) = ρ_ s^τ( 1 - x ),


where τ designates the decreasing rate given by 


    τ = H ( N^2/g + g/c_s^2).


The acoustic and sphericity terms (Eq. <ref>) are simplified into 


    [ ε_s ; n≈R^2 σ^2/Λ_n^m,ν̃ c_s^2                                                 K_∘≈ - 2 H/R. ]


It follows that Φ_n (see Eq. <ref>) writes 


    [   Φ_n ( x ) = ^δ x                    δ = 1/2( τ + K_∘). ]


Hence, we express the vertical structure equation (Eq. <ref>) 


    d^2 Ψ_n^m,σ/dx^2 + k̂_n^2 Ψ_n^m,σ =  H^2/R^2Λ_n^m,ν̃( N^2/σσ̃ - 1 ) U_n ^- δ x .


The vertical wavenumber to square k̂_n^2, given by Eq. (<ref>) in the general case, is now the constant 


    k̂_n^2 =  H^2/R^2Λ_n^m,ν̃( N^2/σσ̃ - 1 ) ( 1 - ε_s ; n) + ( g H/c_s^2 + K_∘) N^2 H/g - δ^2.
 



In this expression, the first term is the vertical wavenumber of internal gravity waves. It dominates when | σσ̃| ≪ N^2, i.e. for a stable stratification. Other terms correspond to acoustic and sphericity terms usually ignored in anelastic (c_ s = + ∞) and thin-shell (H/R → 0) approximations. Here, we will only ignore sphericity terms, i.e. assume that K_∘ = 0. We will keep acoustic terms because they do not bring supplementary mathematical complexities in 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 <cit.>, which gives the mean gradient d ρ_0 / dz = - 0.0048   kg.m^-2. As ρ_0 g / c_s^2 ≈ 0.0043   kg.m^-2, the two terms of Eq. (<ref>) are comparable and must be retained. To solve the vertical structure equation, two boundary conditions are required. At x = 0, we set the impenetrable rigid-wall condition ξ_r = 0. At the upper boundary, we apply the usual stress-free condition <cit.>, δ p = g ρ_0 ξ_r. It follows that



    [                                                       Ψ_n ( x ) =                                           Ψ_n^(0)/𝒟_n{𝒟_n ^- δ x.;                                                                   . +   ( 𝒞_n - δ) ^-δ[ 𝒜_n sin( k̂_n x ) - k̂_n cos( k̂_n x ) ]  .;                                                                                     . +  ( 𝒜_n - δ)  k̂_n  cos( k̂_n ( 1 - x ) )  .;                                                                                     . +  ( 𝒜_n - δ) 𝒞_n sin( k̂_n ( 1 - x ) )  }  , ]











where Ψ_n^(0) is the constant given by 


    Ψ_n^(0) =  H^2 Λ_n^m,ν̃/R^2 ( k̂_n^2 + δ^2 )( N^2/σσ̃ - 1  ) U_n,


and 𝒞_n and 𝒟_n the dimensionless coefficients expressed as


    𝒞_n    =   𝒜_n + H/g( N^2 - σσ̃), 
    𝒟_n    =   k̂_n ( 𝒞_n - 𝒜_n ) cos( k̂_n ) - ( 𝒜_n 𝒞_n + k̂_n^2 ) sin( k̂_n ).









 §.§ Second order tidal Love number and tidal torque


By substituting Eq. (<ref>) in Eqs. (<ref>) and (<ref>) we obtain the components of the variation of mass distribution intervening in the Love numbers and tidal torque (see Eqs. <ref> and <ref>) as explicit functions of the internal structure parameters,


    [            ξ_r ;n,k^m,σ( 1 ) = H 𝒬_ξ ; n^m,σ U_k^m,σ                                                      ∫_0^1 δρ_n,k^m,σ( x )   x = ρ_ s𝒬_ρ ; n^m,σ U_k^m,σ, ]


where the frequency-dependent parameters 𝒬_ξ ; n^m,σ and 𝒬_ρ ; n^m,σ are expressed as


    𝒬_ξ ; n^m,σ  =      - Λ_n^m,ν̃/R^2 σσ̃( k̂_n^2 + δ^2 ){𝒜_n - δ + 𝒟_n^-1[ ℰ_n + ℱ_n sin( k̂_n ) ] } , 
    𝒬_ρ ; n^m,σ =     - N^2 H Λ_n^m,ν̃/g R^2  σσ̃( k̂_n^2 + δ^2 ){τ^-1( ℬ_n - δ)  ( ^τ - 1 ) . 
       . + 𝒢_n + ℋ_n cos( k̂_n ) + 𝒦_n sin( k̂_n ) /𝒟_n (  γ^2 + k̂_n^2 ) },


with γ = δ - τ and the dimensionless coefficients




    ℰ_n  =    k̂_n ( 𝒜_n - δ) ( 𝒜_n - 𝒞_n ) ^δ,



    ℱ_n =    ( 𝒞_n - δ) ( 𝒜_n^2 + k̂_n^2 ) ,



    𝒢_n  =    k̂_n {(𝒜_n - δ) e^δ[ γ( ℬ_n - 𝒞_n ) + ℬ_n 𝒞_n + k̂_n^2 ]  . 
       . + ( 𝒞_n - δ) ^- γ[ γ( ℬ_n - 𝒜_n ) + 𝒜_n ℬ_n + k̂_n^2 ] },



    ℋ_n =     - k̂_n { e^τ( 𝒜_n - δ) [ γ( ℬ_n - 𝒞_n ) + ℬ_n 𝒞_n + k̂_n^2 ] . 
       . + ( 𝒞_n - δ) [ γ( ℬ_n - 𝒜_n ) + 𝒜_n ℬ_n + k̂_n^2 ]  },



    𝒦_n  =    ( 𝒞_n - δ) [ γ( 𝒜_n ℬ_n + k̂_n^2 ) + k̂_n^2  ( 𝒜_n - ℬ_n ) ] 
        - e^τ( 𝒜_n - δ) [ γ( ℬ_n 𝒞_n + k̂_n^2 ) + k̂_n^2 ( 𝒞_n - ℬ_n ) ].



Hence, 𝒬_ξ ; n^m,σ stands for the intrinsic response of the planet due to the variations of the oceanic surface level, while 𝒬_ρ ; n^m,σ corresponds to the contribution of internal inertial-gravity waves. This later will 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 surface gravity waves is weighted by the ratio | 𝒬_ρ ; n  / 𝒬_ξ ; n|. By using the expressions of the solution, we retrieve the weighting factor given in the literature <cit.>, that is | 𝒬_ρ ; n  / 𝒬_ξ ; n| ∼ N^2 H / g. The oceanic second order Love number k_2^2 is then deduced from Eqs. (<ref>), (<ref>) and (<ref>) straightforwardly. In the quadrupolar approximation, where terms of degrees l > 2 are neglected, it writes


    k_2^2 = 𝒢 M_ oc/5R∑_n ∈ℤ C_2,n,2^2,ν̃( 𝒬_ξ ; n^2,σ + 𝒬_ρ ; n^2,σ),


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 Q = |  k_2^2 / { k_2^2 }| <cit.> is expressed as


    Q = |  ∑_n  C_2,n,2^2,ν̃( 𝒬_ξ ; n^2,σ + 𝒬_ρ ; n^2,σ) /∑_n {  C_2,n,2^2,ν̃( 𝒬_ξ ; n^2,σ + 𝒬_ρ ; n^2,σ) }|,


and the quadrupolar tidal torque, which corresponds to the l = m =  2 component, as


    𝒯^2,σ = 1/2 M_ oc| U_2^2,σ|^2  ∑_n{ C_2,n,2^2,ν̃(  𝒬_ξ ; n^2,σ + 𝒬_ρ ; n^2,σ)  }.


where the quadrupolar component of the tidal potential U_2^2,σ is given by


    U_2^2,σ = √(3/5)( R/a)^2 𝒢 M_⋆/a.


Expressed as a function of the degree-2 tidal Love number given by Eq. (<ref>), the tidal torque is written 


    𝒯^2,σ = 3/2𝒢 M_⋆^2 R^5/a^6{ k_2^2 },


which is the well-known expression of the tidal torque associated with the quadrupolar semidiurnal tide <cit.>.





 §.§ Case of the neutrally stratified ocean (N = 0)



In the analytic modeling of oceanic tides, the stratification is usually not taken into account, the layer being two-dimensional and the fluid assumed to be incompressible <cit.>. This reduces the oceanic tidal dynamics to horizontal flows. We show here that we recover the results given by this approach with our modeling, in which the shallow water 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


    k̂_n = i H/R√(Λ_n^m,ν̃).


As shown by Eq. (<ref>), k̂_n is now directly proportional to the  horizontal wavenumber k̂_⊥ ; n = √(Λ_n^m,ν̃ / R^2 ). In the weak-friction approximation (σ_ R≪| σ|) and the regime of super-inertial waves (| ν̃| ≤ 1), it does not depend on the tidal frequency. Indeed, in this case, the associated eigenvalues are Λ_n^m,ν̃≈Λ_n^m,0 = ( m + n  ) ( m + n + 1 ) n ∈ℕ (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 Λ_n^m,ν̃ cannot be ignored any more in the regime of sub-inertial waves. Because of the thin-layer approximation, | k̂_n | ≪ 1. It follows that sin( k̂_n x ) ≈k̂_n x and cos( k̂_n x ) ≈ 1. Hence, in the case of the neutrally stratified incompressible ocean,  the analytic solution of Eq. (<ref>) simply reduces to


    Ψ_n^m,σ =  U_n  g H k̂_⊥ ; n^2 - H/Rσσ̃ x /( σ - σ_n^- ) ( σ - σ_n^+ ) ,


where σ_n^ - and σ_n^ + are the frequencies of resonance of large scale (| k̂_⊥ ; n R | ≪ 1) damped surface inertial-gravity waves (see Fig. <ref>) associated with the gravity mode of degree n. These frequencies characterize the tidal response of the two-dimensional incompressible ocean, where only surface gravity waves can propagate <cit.>. They are expressed as


    σ_n^± =  i  σ_ R/2±√(  g H k̂_⊥ ; n^2 - ( σ_ R/2)^2 ) .


We recognize in this expression the dispersion relation of large scale surface gravity waves, σ^2 = g H k̂_⊥^2 <cit.>. Like the vertical wavenumber, σ_n^ - and σ_n^ + formally depend on σ through Λ_n^m,ν̃ but this dependence can be neglected in the regime of super-inertial waves if the weak-friction approximation is assumed. The above approximations imply τ = 0 and 𝒬_ρ ; n^m,σ = 0. Thus, the Love numbers and tidal torque (Eqs. <ref> and <ref>) are determined by the contribution of the surface displacement solely, which reads


    𝒬_ξ ; n^m,σ = - Λ_n^m,ν̃/R^2 1/( σ - σ_n^- ) ( σ - σ_n^+ ).


The expression of Q_ξ,n^m,σ 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 σ_n^±. 



§ APPLICATION TO ILLUSTRATIVE CASES



In this section, we apply the results established previously to representative telluric planets by using the thin-shell approximation of Sect. <ref>. We first treat the case of the Earth, which illustrates the thin-layer 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 <cit.>. We then apply the model to Trappist-1 f, which could be covered with a deep global ocean of liquid water owing to its small density <cit.>.







 §.§ 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 Earth-Moon system. We focus here on the Lunar quadrupolar tide, which corresponds to the tidal frequency σ = 2 ( Ω - n_ orb), 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 <cit.>, we use for the uniform oceanic depth the mean depth of the Earth's ocean, that is H = 4 km <cit.>. The radius of the planet is 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 <cit.>, 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 <ref>.



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 <cit.>. 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. <ref>) and of the associated tidal quality factor (Eq. <ref>) are plotted in Fig. <ref> (left column) as functions of the normalized frequency ω = ( Ω - n_ orb) / Ω_ (Ω_ stands for the today Earth rotation rate) for various orders of magnitude of σ_ R. We observe on these plots the resonances associated with surface gravity modes modified by rotation. They correspond to the eigenfrequencies σ_n^± given by Eq. (<ref>). As σ_ R decreases, the variability of { k_2^2} increases. Particularly, the level of the non-resonant background decreases proportionally to σ_ R, while the peaks heights increase as σ_ R^-1, which is in good agreement with the scaling laws derived in <cit.>. In the asymptotic regime of strong friction, the behaviour of the ocean is regular. The imaginary part of the tidal Love number scales as { k_2^2 } ∝ σ in the zero-frequency limit (σ→ 0) and as { k_2^2 } ∝ σ^-3 at | σ| → + ∞. 


The observed receding of the Moon is estimated to 3.82 cm per year <cit.>. This corresponds to { k_2^2 }≈ 2.56 × 10^-2 for the tidal Love number and 𝒯≈ 4.50 × 10^16  N.m 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. <ref>, right panel). We hence retrieve the ∼ 30 hours friction time scale estimated by <cit.> with a similar approach. This value of σ_ R will be used in Section <ref> to explore the domain of parameters. 






 §.§ TRAPPIST-1 f


To illustrate the behaviour of a deep stably-stratified ocean, we consider the case of TRAPPIST-1 f. TRAPPIST-1 f is one of the eight telluric planets recently discovered in the vicinity of the star TRAPPIST-1 <cit.>. 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 ρ̅ = 1740   kg.m^-3, that is far smaller than those of rocky planets (for instance, the mean density of the Earth is ρ̅_ = 5514   kg.m^-3; see https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.htmlNASA fact sheets). This means that water could stand for an important fraction of the planet mass. Typically, the water fraction is estimated between 25 % and 100 % of the planet mass. Moreover, with a greenhouse effect, TRAPPIST-1 f is sufficiently irradiated by its host star to have a surface temperature compatible with liquid water <cit.>. All these features argue for the possible presence of a deep oceanic layer on TRAPPIST-1 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. <ref>) is not necessary very small with respect to 1. Thus, as in the case of the Earth, the tidal response is mainly composed of surface gravity waves modified by rotation. But it is also composed of internal gravity waves, restored by the Archimedean force and inducing internal density variations (Fig. <ref>). 

We investigate this complex behaviour by considering an idealized TRAPPIST-1 f planet with a global ocean of depth H = 1000 km and a stratification similar to that of the Earth's ocean, i.e. 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 <ref>) by letting Ω vary with σ. Similarly to the previous case, the corresponding spectra of { k_2^2 } and Q are plottted on Fig. <ref> (right column) as a function of the forcing frequency. We retrieve in the high-frequency range the resonances associated with surface gravity modes. The two peaks appearing around ω≈± 4.5 correspond to the Hough mode of degree n = 0, which is the surface gravity mode of lowest eigenfrequency. Given that σ_n^± ∝ √(gH) in the weakly frictional regime (see Eq. <ref>), the resonances are translated towards the high-frequency range with respect to the case of the Earth (left column). In addition to surface gravity modes, we observe resonances caused by the propagation of internal gravity waves in the low frequency range. These modes can exist for | σ| ≲ N (see Fig. <ref>), which corresponds to | ω| ≲ 5.








§ FROM SHALLOW TO DEEP OCEANS



In the previous sections, we have 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 surface gravity modes, the Brunt-Väisälä frequency the frequency range of internal gravity modes. Besides, 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 , 1000 km) and Brunt-Väisälä frequency (N =10^-4, 10^-3  s^-1), all the other parameters being unchanged (see Table <ref>). 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. <ref> 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 Moon-like satellite (M_ pert = 7.346 × 10^22 kg), we plot the corresponding tidal torque exerted on planets in Fig. <ref>. 



 As it may be noted, the three maps plotted for H = 10 km (bottom panels of Figs. <ref> and <ref>) 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 surface gravity waves. In other words, the tidal response is dominated by its barotropic component, which induces that the role played by the stratification appears to be negligible in this case, as mentioned before. We identify the resonances due to surface inertial-gravity modes observed in the case of the Earth in the previous section (Fig. <ref>, left column). They form a typical wings pattern, which is proper to these waves <cit.>. However, 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 super-synchronous regime (area located below the diagonal line), they induce a distortion that couples the quadrupolar forcing to several Hough modes. These inertial-gravity waves tend to be mixed with acoustic waves while P_ rot→ 0. In the sub-synchronous regime (area located above the diagonal line), the forcing frequency is greater than the inertia frequency (super-inertial 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 super-inertial regime. As the tidal potential scales as U_2^2,σ ∝  n_ orb^2 (see Eq. <ref>), the tidal torque scales as P_ orb^-4 (see Eq. <ref>), which corresponds to the horizontal color gradient of the maps plotted on Fig. <ref>. 

By moving to deeper oceans, we observe that the resonances are translated towards the small-period range, following the scaling law σ_n^± ∝ √(gH) 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. <ref> and <ref>). This effect is enhanced for H = 1000 km. As the eigenfrequencies of internal gravity waves do not depend on Ω like those of surface 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 non-resonant 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 the tidal torque and Love number on the forcing frequency <cit.>.



§ DISCUSSION



In this work, we have opted for a simplified linear analysis allowing us to explore widely the domain of parameters with a reasonable computational cost. This approach is convenient to examine the frequency-resonant behaviour of a global ocean and highlight the key parameters of the problem, which can be explained in more details 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. <ref>). The case of a hemispherical ocean centred at the equator was treated in early studies in the shallow water approximation <cit.>. Meridian continental shelves prevent large scales gravity waves to propagate. They thus couple the global gravitational forcing to modes determined by the scale of the basin. This modifies significantly the aspect of the frequency spectrum of the tidal response <cit.>. 



  ∙ 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 litterature <cit.> and assumed for convenience. In addition to the traditional approximation, discussed in the next 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 small scales topographical features homogeneously distributed around the planet, as in the numerical model of <cit.>. It can also be considered as a simplified approximation of viscous and turbulent frictions if a lengthscale L such that ΔV∼ V / L^2 is introduced. However, it does not model at all of the local dissipation due the breaking of waves on continental shelves whereas 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 <cit.>. 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 well appropriate to 2D models, where vertical motions are negligible relatively 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 super-inertial 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 to use 3D (2D if solutions ∝ ^i m φ are assumed) numerical or semi-analytical methods <cit.>. 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 to be comparable in amplitude to horizontal ones. This condition is expressed as | 2 Ω| ≪ N <cit.>. The validity of the traditional approximation has been examined in various fluid layers, from planetary oceans and atmospheres <cit.> to stellar interiors <cit.>. 



  ∙ 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 <cit.>. By studying the propagation of internal gravity waves in a shear flow, <cit.> showed that the amplitude of oscillations is attenuated by a factor depending on the local Richardson number of the fluid, R_ i = N^2  | d  V_0  / dr |^-2 (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 way, tidal waves can transfer horizontal moment to mean flows, which modifies the tidal response. 



  ∙ Free-surface boundary condition – In this work, a free-surface 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 Europa-like icy satellites <cit.>. 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 free-surface condition by a rigid wall 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 visco-elastic body <cit.> forced both by the tidal potential and the loading induced by the tidal distortion of the ocean. This solid-fluid coupling affects the frequencies given by Eq. (<ref>) and generates an additional coupling between tidal modes <cit.>. We simply ignored this coupling so that we do not need to consider the parameters of the solid response.






§ 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 leaves the traditional 2D modeling inherited from the study of the Earth's oceans to include three-dimensional 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 <cit.>. 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 thin-layer 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 TRAPPIST-1 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 <cit.>. The tidal response is composed of resonant surface inertial-gravity 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 high-frequency 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 frequency-resonant behaviour, like surface gravity waves. The contribution of this component is not negligible in the case of deep oceans and sensitively increases the evolution timescales of the planet-perturber 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. It should be considered as a next step to introduce meridian 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 continental shelves in the oceanic tidal dissipation.  





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. 





aa  
 



