EDP Sciences
Free Access
Issue
A&A
Volume 603, July 2017
Article Number A107
Number of page(s) 44
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201628252
Published online 19 July 2017

© ESO, 2017

1. Introduction

Since the pioneering work of Kelvin (1863) on the gravitationally forced elongation of a solid spherical shell, tides have been recognized as a phenomenon of major importance in celestial mechanics. They are actually one of the key mechanisms involved in the evolution of planetary systems over long timescales. Because they affect both stars and planets, tides introduce a complex dynamical coupling between all the bodies that compose a system. This coupling results from the mutual interactions of bodies and is tightly bound to their internal physics. Solid tides have long been the subject of studies (see the reference works by Love 1911; Goldreich & Soter 1966). The corresponding tidal dissipation, which is caused by viscous friction, can have a strong impact on the orbital elements of a planet (e.g. Mignard 1979, 1980; Hut 1980, 1981; Neron de Surgy & Laskar 1997; Henning et al. 2009; Remus et al. 2012). This evolution depends smoothly on the tidal frequency (e.g. Efroimsky & Lainey 2007; Greenberg 2009). Fluid tides cause internal dissipation as well as solid tides, and the energy dissipated is even more dependent on the tidal frequency than the solid tides are because of a highly resonant behaviour (Ogilvie & Lin 2004; Ogilvie 2014; Auclair Desrotour et al. 2015). As a consequence, in both cases the phase shifted elongation of the shell induces a tidal torque which modifies the rotational dynamics of the body. This torque drives the evolution of the spin of planets and determines its possible states of equilibrium (Gold & Soter 1969; Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001; Correia et al. 2003, 2014; Correia & Laskar 2003; Arras & Socrates 2010).

The increased focus on exoplanets over the past two decades and the rapidly increasing number of discovered planetary systems(Mayor & Queloz 1995; Perryman 2011) have brought to light a huge diversity of existing orbital structures that can be very different from our own solar system (Fabrycky et al. 2012). Therefore, tidal effects now need to be studied in a general context to characterize these new systems, understand their orbital and rotational dynamics, and constrain the physics of exoplanets with the information provided by observational measurements.

In this work, we focus on Earth-like exoplanets, which are interesting for their complex internal structure and possible habitability. Indeed, a planet such as the Earth or Venus is composed of several solid and fluid layers with different physical properties, which implies that the tidal response of a layer cannot be simply correlated to its size. For instance, the Earth’s oceans, in spite of their very small depth compared to the Earth’s radius, are responsible for most of the tidal dissipation of the planet (Webb 1980; Egbert & Ray 2000; Ray et al. 2001). The atmospheric envelope of a terrestrial exoplanet sometimes represents a non-negligible fraction of its mass (about 20% in the case of small Neptune-like bodies). It is also a layer of great complexity because it is submitted to both tidal gravitational potentials due to other bodies and thermal forcing of the host star. Therefore, it is necessary to develop theoretical models of atmospheric tides with the smallest possible number of parameters. These models will be useful tools to explore the domain of parameters, and to quantify the tidal torque exerted on the spin axis of a planet and quantities used in the orbital evolution of planetary systems. One of these quantities is the second-order Love number (k2) introduced by Love (1911), which measures the quadrupolar hydrostatic elongation. The adiabatic Love number k2 has been estimated in the solar system (e.g. Konopliv & Yoder 1996; Lainey et al. 2007; Williams et al. 2014).

According to Wilkes (1949), the first theoretical global models of atmospheric tides were developed at the end of the eighteenth century by Laplace, who published them in Mécanique céleste (Laplace 1798). Laplace was interested in the case of a homogeneous atmosphere characterized by a constant pressure scale height and affected by tidal gravitational perturbations. He thus founded the classical theory of atmospheric tides. At the end of the nineteenth century, Lord Kelvin pointed out the unexpected importance of the Earth’s semidiurnal pressure oscillations in spite of the predominating diurnal perturbation (Kelvin 1882; Hagan et al. 2003; Covey et al. 2009, see also the measures of the daily variation of temperature and pressure, Fig. 13). This question motivated many further studies (e.g. Lamb 1911, 1932; Taylor 1929, 1936; Pekeris 1937) which contributed to enriching the classical theory of tides. In the 1960s, theoretical models were generalized in order to study other modes, such as diurnal oscillations (see Haurwitz 1965; Kato 1966; Lindzen 1966, 1967a,b, 1968), but still remained focused on the case of the Earth. They revealed the impact of thermal tides due to solar insolation by showing that these tides drive the Earth’s diurnal and semi-diurnal waves.Indeed, in the Earth atmosphere, the contribution of the tidal gravitational potential, which causes the elongation of solid layers, is negligible compared to the effect of the thermal forcing. A very detailed review of the general theory of thermal tides and of its history is given by Chapman & Lindzen (1970), hereafter CL70 (see also Wilkes 1949; Siebert 1961).

Leaning on this early work, here we generalize the existing theory to Earth-like planets and exoplanets. We introduce the mechanisms of radiative losses and thermal diffusion, which are usually not taken into account in geophysical models of atmospheric tides owing to their negligible effects on the Earth’s tidal response (radiation was, however, introduced with a Newtonian cooling by Lindzen & McKenzie 1967; Dickinson & Geller 1968, which inspired the present work). These mechanisms cannot be ignored because they play a crucial role in the case of planets like Venus, which are very close to spin-orbit synchronization. Models based on perfect fluid hydrodynamics, such as the one detailed by CL70, fail to describe the atmospheric tidal response of such planets because they are singular at synchronization. In these models, the amplitude of perturbed quantities like pressure, density, and temperature diverges at this point. The case of Venus was studied by Correia & Laskar (2001), who proposed an empirical bi-parameter model for the tidal torque describing two possible regimes: the regime of an atmosphere where dissipative processes can be ignored, far from synchronization, and its supposed regime at the vicinity of synchronization, where the perturbation is expected to annihilate itself. This behaviour was recently retrieved with global climate models (GCM) simulations by Leconte et al. (2015), who took into account dissipative mechanisms, thus showing that these mechanisms drive the response of planetary atmospheres in the neighbourhood of synchronization.

Our objective in this work is thus to propose a new theoretical global model controlled by a small number of physical parameters which could be used to characterize the response of any exoplanet atmosphere to a general tidal perturbation. We aim to answer the following questions:

  • How does the amplitude of pressure, density, temperature, andvelocity oscillations depend on the parameters of the system (i.e.rotation, stratification, radiative losses)?

  • What are the possible regimes for atmospheric tides?

  • How do Love numbers and tidal torques vary with tidal frequency?

First, in Sect. 2, we derive the general theory for thick atmospheres, which are characterized by a depth comparable to the radius of the planet. In this framework, we expand the perturbed quantities in Fourier series in time and longitude, and in series of Hough functions (Hough 1898) in latitude. This allows us to compute the equations giving the vertical structure of the atmospheric response to the tidal perturbation. Then, from the hydrodynamics we derive the analytic expressions of the variations of pressure, density, temperature, velocity field, and displacement for any mode. At the end of the section, the density oscillations obtained are used to compute the atmospheric tidal gravitational potential, the Love numbers, and the tidal torque exerted on the spin axis of the planet. In Sect. 3, we apply the equations from Sect. 2 to the case of a thin isothermal stably stratified atmosphere. In this simplified case, the terms associated with the curvature of the planet disappear and most of the parameters depending on the altitude, such as the pressure scale height of the atmosphere, become constant. In Sect. 4, we treat the case of slowly rotating convective atmospheres, such as near the ground in Venus (Marov et al. 1973; Seiff et al. 1980). In the next section, we derive the terms of thermal forcing from the Beer-Lambert law (Bouguer 1729; Klett et al. 1760; Beer 1852) and the theory of temperature oscillations at the planetary surface. In Sect. 6, the model is applied to the Earth and Venus. The corresponding torque and atmospheric Love numbers are computed from analytical equations. We end the study with a discussion and by giving our conclusions and prospects. To facilitate the reading, several technical issues have been deferred to the Appendix where an index of notations is also provided.

2. Dynamics of a thick atmosphere forced thermally and gravitationally

The reference book CL70 has set the basis of analytical approaches for atmospheric tides. This pioneering work focuses on fast rotating telluric planets covered by a thin atmosphere, such as the Earth.

The goal of the present section is to generalise this work. We establish the equations that govern the dynamics of tides in a thick fluid shell in the whole range of tidal frequencies, from synchronization to fast rotation. We consider a spherical telluric planet of radius R covered by a stratified atmospheric layer of typical thickness Hatm. The atmosphere rotates uniformly with the rocky part at the angular velocity Ω (the spin vector of the planet being denoted Ω). Therefore, the dynamics are written in the natural equatorial reference frame rotating with the planet, RE;T:{O,XE,YE,ZE}, where O is the centre of the planet, XE and YE define the equatorial plane and ZE = Ω/ |Ω|. We use the spherical basis and coordinates (r,θ,ϕ), where r is the radial coordinate, θ the colatitude, and ϕ the longitude. As usual, t stands for the time.

thumbnail Fig. 1

Spherical coordinate system associated with the equatorial reference frame of the planet, and the geometrical parameters of the system.

Open with DEXTER

The atmosphere is assumed to be a perfect gas homogeneous in composition, of molar mass M, and stratified radially. Its pressure, density, and temperature are denoted p, ρ, and T respectively. For the sake of simplicity, all dissipative mechanisms (e.g. viscous friction and heat diffusion) are ignored except for the radiative losses of the gas, which play an important role in the vicinity of synchronization. Tides are considered to be a small perturbation around a global static equilibrium state. The basic pressure p0, density ρ0, and temperature T0, are supposed to vary only with the radial coordinate. The tidal response is characterized by a combination of inertial waves, of typical frequency , due to Coriolis acceleration, and gravity waves which are restored by the stable stratification. The typical frequency of these waves, the Brunt-Väisälä frequency (e.g. Lighthill 1978; Gerkema & Zimmerman 2008), is given by (1)where Γ1 = (lnp0/lnρ0)S is the adiabatic exponent (the index S being the specific macroscopic entropy) and g the gravity. We assume that the fluid follows the perfect gas law, (2)where Rs = RGP/M is the specific gas constant of the atmosphere, RGP = 8.3144621 J mol-1 K-1 (Mohr et al. 2012) being the molar gas constant. In the following, the notations , , and are used to designate the real part, imaginary part, and conjugate of a complex number. The subscript symbols and stand for the Earth and Venus, respectively, which are taken as examples to illustrate the theory. All functions and fields displayed in the figures are computed with the algebraic manipulator TRIP (Gastineau & Laskar 2014) and plotted with gnuplot.

2.1. Forced dynamics equations

The atmosphere is affected by the tidal gravitational potential U and the heat power per mass unit J. In this approach, the equations of dynamics are linearized around the equilibrium state. The small variations in this state induced by the tidal perturbation are the velocity field V = Vθeθ + Vϕeϕ + Vrer; the corresponding displacement field ξ = ξθeθ + ξϕeϕ + ξrer defined as (3)and the pressure (δp), density (δρ), and temperature (δT) fluctuations where (4)This represents six unknown quantities to compute (Vr, Vθ, Vϕ, δp, δρ, δT), and we therefore solve a six-equation system. We first introduce the Navier-Stokes equation. We assume the Cowling approximation in which the self-gravitational potential fluctuations induced by tides are neglected (Cowling 1941). It is well adapted to waves characterized by a high vertical wavenumber as are found in the vicinity of synchronization, as detailed farther. In the equatorial reference frame RE;T, the linearized Navier-Stokes equation is thus written (5)which, projected on er, eθ, and eϕ, gives (6)(7)(8)Equations (6)–(8) are coupled by the Coriolis terms (characterized by the factor ) on the left-hand side. To simplify the equations, we assume the traditional approximation, i.e. to neglect the terms 2ΩsinθVr and 2ΩsinθVϕ in Eqs. (7) and (8), respectively. This hypothesis, commonly used in literature dealing with planetary atmospheres and star hydrodynamics (e.g. Eckart 1960; Mathis et al. 2008), can be applied to stratified fluids where the tidal flow satisfies two conditions: (i) the buoyancy (given by gδρ/ρ0er in Eq. (5)) is strong compared to the Coriolis term (2ΩV) and (ii) the tidal frequency (introduced below) is greater than the inertia frequency () and smaller than the Brunt-Väisälä frequency. These conditions are satisfied by the Earth’s atmosphere and fast rotating planets. The quantitative precision that can be expected with the traditional approximation may be lower in the vicinity of synchronization where the tidal frequency tends to zero; moreover, we note that it leads to issues with momentum and energy conservation in the case of deep atmospheres (for a discussion on the limitations of the traditional approximation, see Gerkema & Shrira 2005; White et al. 2005; Mathis 2009).

Hence, we obtain: (9)(10)(11)In strongly stratified fluids, the left-hand side of Eq. (11) is usually ignored because it is very small compared to the terms in δp and δρ in typical wave regimes. The radial acceleration will only play a role in regimes where the tidal frequency is close to or exceeds the Brunt-Väisälä frequency, which corresponds to fast rotators.

The second equation of our system is the conservation of mass, (12)which, in spherical coordinates, writes (13)The thermal forcing (J) appears on the right-hand side of the linearized heat transport equation (see CL70, Gerkema & Zimmerman 2008) given by (14)where and Jrad is the power per mass unit radiated by the atmosphere, supposed to behave as a grey body. We consider that JradδT. This hypothesis is known as “Newtonian cooling” and was used by Lindzen & McKenzie (1967) to introduce radiation analytically in the classical theory of atmospheric tides (see also Dickinson & Geller 1968). Physically, it corresponds to the case of an optically thin atmosphere in which the flux emitted by a layer propagates upwards or downwards without being absorbed by the other layers. In optically thick atmospheres, such as on Venus (Lacis 1975), this physical condition is not verified. Indeed, because of a stronger absorption, the power emitted by a layer is almost totally transmitted to the neighbourhood. Therefore, this significant thermal coupling within the atmospheric shell should ideally be taken into account in a rigorous way, which would lead to great mathematical difficulties (e.g. complex radiative transfers, Laplacian operators) in our analytical approach. However, recent numerical simulations of thermal tides in optically thick atmospheres by Leconte et al. (2015) show behaviour of the flow that is in good agreement with a model that uses radiative cooling. Therefore, in this work we assumeNewtonian cooling as a first model of the action of radiation on atmospheric tides. Newtonian cooling brings a new characteristic frequency, denoted σ0, which we call radiative frequency and which depends on the thermal capacity of the atmosphere. The radiative power per unit mass is thus written (15)Like the basic fields p0, ρ0, and T0, the radiative frequency varies with r and defines the transition between the dynamical regime, where the radiative losses can be ignored, and the radiative regime, where they predominate in the heat transport equation. Assuming that the radiative emission of the gas is proportional to the local molar concentration C0 = ρ0/M, it is possible to express Jrad and σ0 as functions of the physical parameters of the fluid (Appendix D), (16)and (17)the parameter ϵa being an effective molar emissivity coefficient of the gas and S= 5.670373 × 10-8 W m-2 K-4 the Stefan-Boltzmann constant (Mohr et al. 2012). The substitution of Eq. (15) in Eq. (14) yields (18)Finally, the system is closed by the perfect gas law (19)Substituting Eq. (19) in Eq. (18), we eliminate the unknown δT, and obtain (20)Because of the rotating motion of the perturber in the equatorial frame (RE;T), a tidal perturbation is supposed to be periodic in time (t) and longitude (ϕ). So, any perturbed quantity f of our model can be expanded in Fourier series of t and ϕ(21)the parameter σ being the tidal frequency of a Fourier component and m its longitudinal degree1. We also introduce the spin parameter (22)which defines the possible regimes of tidal gravito-inertial waves:

  • |ν| ≤ 1 corresponds to super-inertial waves, for which the tidal frequency is greater than the inertia frequency;

  • |ν| > 1 corresponds to sub-inertial waves, for which the tidal frequency is lower than the inertia frequency.

The parameters m and ν determine the horizontal operators

Here the operator is associated with and , while is associated with and (Lee & Saio 1997; Mathis et al. 2008). Hence, the traditional approximation allows us to write the horizontal component of the velocity field only as a function of the variations in the pressure and tidal gravitational potential

where ym,σ = δpm,σ/ρ0. Since Vm,σ = tξm,σ, we also get the corresponding horizontal component of the displacement vector from Eqs. (25) and (26)

thumbnail Fig. 2

Normalized Hough functions , , and corresponding to the Earth’s semidiurnal tide (m = 2 and ν = 1.0030). We note that ν = 2Ω /σ. Left: first gravity modes; 0 ≤ n ≤ 5. Right: first Rossby modes: − 6 ≤ n ≤ −1. In this case, |ν| is very close to 1 and the Rossby modes are trapped at the poles. The Earth’s semidiurnal tide is mainly described by gravity modes. For more details about the behaviour of Hough functions, see Fig. A.2.

Open with DEXTER

The substitution of Eqs. (25) and (26) in Eq. (13) yields (29)where m,ν is the Laplace operator, parametrized by ν and m and is expressed by (30)with ν = 2Ω /σ (Eq. 22). When ν = 0 (i.e. Ω = 0), m,ν is reduced to the horizontal Laplacian. Since the horizontal velocity field is not coupled with the other variables, only three unknowns remain. So, the new system is written

Under the traditional approximation, solutions of Eqs. (31)(33) can be obtainedunder the form of series of functions of separated coordinates,

(34)where the functions are the eigenvectors of the Laplace operator, the corresponding eigenvalues being denoted (with ν = 2Ω /σ). Laplace’s tidal equation, expressed as (35)describes the horizontal structure of the perturbation. The solutions of Eq. (35), are called Hough functions (Hough 1898)2. They form a complete set of continuous functions on the interval θ ∈ [0], parametrized by m and ν. They also verify the same boundary conditions as the associated Legendre polynomials (denoted ), which are solutions of Eq. (35) for ν = 0 with . In Figs. 2 and 3, the functions , , and are plotted for the Earth and Venus in the case of the solar semidiurnal tide defined by m = 2, σ = 2(Ω−norb), and ν = Ω/(Ω−norb), where norb stands for the mean motion of the planet. In the case of the Earth, the set is mainly composed of gravity modes. Since the spin parameter is slightly greater than 1, Rossby modes exist but are trapped at the poles. In the case of Venus, Hough functions are composed of gravity modes only because |ν| < 1. These functions are very similar to the associated Legendre polynomials , with l ≥ 2.

thumbnail Fig. 3

First gravity modes (0 ≤ n ≤ 5) of normalized Hough functions , , and corresponding to the semidiurnal tide of Venus-like exoplanets (m = 2 and ν = 0.4804). We note that ν = 2Ω /σ. For more details about the behaviour of Hough functions, see Figs. A.2.

Open with DEXTER

The system composed of Eqs. (31) to (33) can then be reduced to the pair of first-order differential equations (36)with the coefficients (37)and (38)where we introduce the sound velocity (39)The system of Eq. (36) can itself be reduced to a single equation in alone (40)with the set of coefficients (41)Equation (40) gives us the vertical structure of tidal waves generated in the fluid shell by both gravitational and thermal forcings. We have: (42)The last term, (dB1/ dr) /B1, is associated with curvature and can be ignored in the context of the “shallow atmosphere” approximation, as discussed in the next section. It is expressed as (see Eq. (50)) (43)We now transform Eq. (40) in a Schrödinger-like equation by introducing the radial function (44)and the variable change (45)where the volume is the wave function of tidal waves associated with the mode of degree n. We note that fixing σ0 = 0 allows us to retrieve the usual variable change, (Press 1981; Zahn et al. 1997; Mathis et al. 2008). With the general variable change, Eq. (40) becomes (46)which gives us the vertical profile of the perturbation3. The form of this equation is very common in the literature because it describes the radial propagation of waves in a spherical shell. The typical wavenumber of these waves, denoted , is given by (47)We introduce the corresponding normalized length scale of variations of the mode of degree n, (48)This parameter is the typical height over which spatial variations of the mode can be observed.

2.2. Polarization relations

The perturbed quantities are readily deduced from Ψ. Before going further, let us introduce the Lamb frequency (the cut-off frequency of acoustic waves) associated with the mode (ν,m,n),(49)and the ratio (50)The parameter |εs,n| measures the relevance of the anelastic approximation, which consists in neglecting all the terms that correspond to an acoustic perturbation in Eq. (36). When |εs,n| ≪ 1 (i.e. |σ| ≪ |σs;n|), this hypothesis can be assumed. However, particular caution should be used when ignoring some terms. For instance, since the first gravity mode of the Earth’s semidiurnal tide is characterized by |σ| ≈ |σs;0| ≈ 10-4 s-1, the anelastic approximation cannot be done. In order to be able to study cases of the same kind, the polarization relations that follow will be given in their most general form. Hence, by substituting the series of Eq. (34) in Eqs. (25)–(28), and (31)–(33), we obtain the radial profiles , , , , , , , , , and as functions of and its first derivative. We obtain for the Lagrangian displacement

for the velocity

and for scalar quantities (57)

The coefficients in these expressions are given by

(60)where we have introduced the curvature term (61)which will be ignored in the thin atmosphere approximation. These relations and Eq. (46) define a linear operator, denoted Υσ, which determines the vertical structure of the atmospheric tidal response. Let be the output vector of the model and the input vector. Eqs. (51) to (59) can be reduced to (62)For any quantity f, the component of describing the vertical profile of f is denoted ( for the pressure, for the density, etc.).

Polarization relations can be naturally divided into two distinct contributions. The first contribution, which is called the “horizontal part” in this work, is a linear combination of the forcings and . It corresponds to the component of the response that does not depend on the tidal vertical displacement. The length scale of this horizontal component is the typical thickness of the atmosphere (Hatm). The other part, expressed as a function of and called “vertical part”, thus results from the tidal fluid motion along the vertical direction. In the following, we use the subscripts H and V to designate the horizontal and vertical components.

As can be noted, polarization relations point out the necessity of taking into account dissipative processes to study tidal regimes close to synchronization. Indeed, without dissipation (σ0 = 0), synchronization (σ = 0) is characterized by a singularity. The terms associated with the horizontal component, which are directly proportional to and , tend to infinity at σ → 0. Terms in and behave similarly, as demonstrated in the next section. These terms are highly oscillating along the vertical direction over the whole depth of the atmosphere in this frequency range because when σ → 0. Hence, radiation regularizes the behaviour around σ = 0 and damps the vertical oscillations of waves associated with the vertical component. Diffusion acts in the same way by flattening the oscillations of smallest length scales (LV;n) (Press 1981; Zahn et al. 1997; Mathis et al. 2008). We return to this point when we solve analytically the case of the isothermal atmosphere with constant profiles for the forcings (Sect. 3).

The whole spectrum of possible tidal regimes is represented in Fig. 4. The radiative regime (blue) characterizes a flow where thermal losses due to the linear Newtonian cooling (see the heat transport equation, Eq. 18) predominate. In the dynamic regime (orange), dissipation can be neglected because the system is governed by the Coriolis acceleration. Finally, the mixed-acoustic regime (red) corresponds to high tidal frequencies comparable to the Lamb frequency, given by Eq. (49). This regime marks the limit of the traditional approximation.

thumbnail Fig. 4

Frequency spectrum of the atmospheric tidal response. The identified regimes are given as functions of the tidal frequency (σ).

Open with DEXTER

2.3. Tidal potential, Love numbers, and tidal torque

The new mass distribution resulting from tidal waves generates a tidal gravitational potential , which is a linear perturbation of the spherical gravitational potential of the planet. This atmospheric potential presents disymmetries with respect to the direction of the perturber. The tidal torque thus induced will affect the rotationnal dynamics of the planet over secular timescales and determine the possible equilibrium states of the spin, as demonstrated by Correia & Laskar (2001). This torque is deduced from the Poisson’s equation expressed in the co-rotating frame (63)the notation G designating the gravitational constant. Like δρ, the potential is expanded in series of functions of separated variables (64)where the (with l ∈N such as l ≥ |m|) are the normalized associated Legendre polynomials. Decomposing δρm,σ on this basis, we get the equation describing the vertical structure of , (65)For the upper boundary condition, the tidal potential is required to remain bounded at r → ∞. At r = 0, we impose the same condition. Therefore, the solution of Eq. (65) is (66)where and are the functions (67)Love numbers are defined as the ratio between the gravitational potential due to the tidal response of the atmosphere and the forcing potential at the upper boundary of the layer. Let us denote Ratm = R + Hatm this upper boundary. Then, at the upper boundary, the atmospheric tidal potential given by Eq. (66) is simply expressed: (68)Considering Eq. (72) and using the linearity property of , the density component can be expanded (69)with = and the change-of-basis coefficients, (70)which allows us to write the atmospheric tidal potential (71)Let us treat the case of a simplified planet-star system, where the planet P circularly orbits around its host star, of mass MS. The semi-major axis, obliquity, and mean motion of the planet are denoted a, ε, and norb. In this case, the modes of the gravitational and thermal forcings applied on the atmosphere can be written (Appendix C) (72)Hence the complex Love numbers can be written generically (73)In celestial dynamics, the second-order Love number (l = 2) is commonly used to quantify the tidal response of a body. Therefore, we illustrate the expression given by Eq. (73) by computing this coefficient for m = 2. For the sake of simplicity, we set the obliquity ε = 0 (the spin of the planet is supposed to be aligned with its orbital angular momentum). The tidal frequency is σSD = 2(Ω−norb) and the series of Eq. (C.14) are reduced to the terms characterized by (l,m,j,p,q) = (2,2,2,0,0)4(74)The tidal potential becomes (75)If the star can be considered a point-mass perturber (Ra), then rapidly decays with k. So, the terms of orders k> 2 are only generated by the thermal forcing. We also assume that the horizontal pattern of J is well represented by the function , which allows us to ignore other terms. Finally, we consider the case |ν| ~ 1 (|norb| ≪ |Ω|, see Eq. (22)), where and for n> 0 (). In this simplified framework, the second-order Love number can be approximated in magnitude order by (76)To conclude this section, we compute the tidal torque exerted on the atmosphere with respect to the spin axis. The monochromatic torque associated with the frequency σ writes (Zahn 1966) (77)where the notation V designates the volume of the atmospheric shell. Denoting the phase lag between the forcing and the response, we obtain (78)and the total tidal torque, (79)As previously done for the Love numbers, is expanded using Hough modes (80)with l ≥ |m|, n ∈Z, and k ≥ |m|. Indeed, we note here that is the component of density variations represented by caused by the component of the excitation represented by and projected on Θn. Hence, the parameter is the phase difference between the component of degree of the response and the component of degree l of the excitation. When the tidal gravitational potential is quadrupolar (Ra), the terms of orders higher than l = 2 can be neglected. Denoting and the quantities and for l = 2 (which should not be confused with the projections and of the forcings on the set of Hough functions), it follows (81)If the forcing is a positive real function, then the previous expression becomes (82)Substituting the polarization relation of density in Eq. (82), we obtain the torque as a sum of two contributions (83)where (84)(85)with and (86)Assuming |εs;n| ≪ 1 and noticing that , we establish that does not depend on the eigenvalues of Laplace’s tidal equation and can be written (87)The torque induced by this component is expanded in series of Hough functions and characterized by the solutions of the vertical structure equation. The typical scales of the vertical component are given by the vertical wavenumbers (see Eq. (47)) and can be very small compared to the scale of the horizontal ones. For example, in the vicinity of synchronization, if the radiative damping is ignored (CL70), which means that the characteristic scales of variations rapidly decay when σ → 0. Finally, we note that the series of modes in Eq. (81) can often be reduced to one dominating term. For instance, in the case of the Earth’s semidiurnal tide, and the expression above (Eq. (81)) can be simplified (88)

3. Tides in a thin stably stratified isothermal atmosphere

We establish in this section the analytical equations that describe the tidal perturbation in the case where HatmR. This theoretical approach has been exhaustively formalized for the Earth in “atmospheric tides” (CL70) and corresponds to the case of thin atmospheres. The formalism developed in the previous section allows us to go beyond this pioneer work, and more particularly to study regimes near synchronization thanks to the inclusion of radiative losses.

3.1. Equilibrium state

For a thin atmosphere (rR), the equilibrium structure can be determined analytically. We introduce the altitude z, such that r = R + z, which is here a more appropriate coordinate than r.

The fluid is supposed to be at hydrostatic equilibrium, which leads to (89)The local acceleration in RE;T, denoted aRE;T, is not equal to g because of the rotating motion. The total acceleration can be written (90)where ac is the centrifugal acceleration: (91)Here, Coriolis acceleration does not intervene because the global circulation of the atmosphere is ignored. The centrifugal acceleration (Eq. (90)) can be neglected with respect to the gravity if the rotation rate of the planet satisfies the condition Ω ≪ Ωc, where is the critical Keplerian angular velocity. In the case of an Earth-like planet, with g ~ 10 m s-2 and R ~ 6.0 × 103 km, Ωc ~ 20 rotations per day. From Eqs. (2) and (89), a typical depth can be identified, (92)which is the local characteristic pressure height scale. The parameter H represents the vertical scale of variations of basic fields. It is of the same order of magnitude as Hatm and allows us to introduce the reduced altitude (93)which is used in the following. Assuming that p0, ρ0, and T0 only depends on x, we deduce their expressions from Eqs. (2) and (89) (94)In the fluid shell, the gravitational acceleration does not vary very much (95)which means that g can be considered as a constant. Assuming that H also does not vary with the altitude, we obtain a constant profile of the temperature and exponentially decaying profiles of the pressure and density: (96)In this case, which corresponds to an isothermal atmosphere, 95% of the mass of the gas is contained within the interval x ∈ [0,3]. Therefore, we can write Hatm ≈ 3H. For the Earth, H ~ 8 km (CL70) and Hatm ~ 24 km.

3.2. Wave equation and polarization relations

The constant H hypothesis is very useful in order to simplify the expressions of Sect. 2. Indeed, with this approximation, the typical frequencies of the system no longer vary with the radial coordinate. The Brunt-Väisälä frequency simply writes (97)and the acoustic frequency associated with the mode , (98)with . Moreover, the horizontal structure of tidal waves does not change because Laplace’s tidal equation is not modified. So, the vertical structure equation and the radial polarization relations only are affected by the constant H hypothesis. The coefficients defined by Eq. (41) become (99)where (100)The parameter hn, which writes (101)is usually called the “equivalent depth” in the literature (e.g. Taylor 1936; Chapman & Lindzen 1970) and represents a typical length scale associated with the mode of degree n. We note that the curvature term in Eq. (42) vanishes in the shallow atmosphere approximation because the coefficient B1 (Eq. (37)) is a constant in this case. This allows us to simplify the vertical structure equation (Eq. (46)). Unlike the vertical wavenumber in the thick shell (Eq. (47)), the vertical wavenumber does not vary with x in this case. It is expressed (102)with the scale ratio . It can also be written (103)The term ςn(εs;n−1) comes from the radial acceleration in the Navier-Stokes equation (Eq. (11)) and is negligible for usual values of (|εs;n| ≪ 1 and ςnH2/R2 ≪ 1). Therefore, by considering the case where σ0 ≪ |σ|, we recover the wavenumber obtained by CL70 for the Earth’s semidiurnal tide. In the general case, the expression of Eq. (103) can be approximated by (104)which gives for the real and imaginary parts, introducing the characteristic depth hc = 4κH and fixing , (105)with (106)

thumbnail Fig. 5

Real and imaginary part of the vertical wavenumber as functions of the tidal (σ) and critical (σc;n) frequencies. Left: , with defined in Eq. (104), as a function of f(σ/σ0) (horizontal axis) and (vertical axis), where f is the function defined by . The colour of the map, denoted c, is defined by . Hence, the purple region stands for slowly oscillating waves. The luminous (dark) regions correspond to strongly oscillating waves propagating downwards (upwards).Right: as a function of the same parameters. The luminous (dark) regions correspond to strongly (weakly) evanescent waves. The positions of the planets (Earth , Venus ) on the map are determined by the gravity mode of degree 0 in the case of the semidiurnal tide (Λ;0 = 11.1596, Λ;0 = 7.1485). The tidal frequencies of the planets are given by σ/σ0 ≈ 194 and σ/σ0 = 0.43, with σ0 = 7.5 × 10-7 s-1.

Open with DEXTER

The behaviour of the vertical wavenumber defines the possible regimes of the perturbation. These regimes are represented by the map in Fig. 5, which shows the real and imaginary part of the vertical wavenumber as colour functions of the tidal frequency and critical frequency, σc;n, given by (the frequency σc;n is proportional to the Lamb frequency defined by Eq. (98)). Propagative modes are characterized by (right panel, dark regions), evanescent modes by (yellow regions). The transition between the radiative regime and the dynamic regime corresponds to |σ| ≈ σ0.

In the case of the thin isothermal atmosphere, the equation giving the vertical profiles (Eq. (46)) writes (107)where C is given by Eq. (99). Moreover, if the condition (108)is satisfied, the forcing in the right-hand member is dominated by thermal tides and the contribution of the gravitational tidal potential can be ignored in Eq. (107), which is rewritten (109)The polarization relations of the thin atmospheric shell are deduced from Eqs. (51) to (59)

with (119)and where , n, and are the constant coefficients of Eq. (60)5, (120)

3.3. Boundary conditions

Solving Eq. (107) requires that we choose two boundary conditions. Following CL70, we fix ξr = 0 at the ground. This corresponds to a smooth rigid wall and is equivalent to at x = 0. In the case of Earth-like exoplanets where the telluric surface is well defined, this condition is relevant. Thus, given the form of Eq. (107), can be written (121)where A and B are the complex functions (122)The parameter K ∈C is an integration constant which is fixed by the upper boundary condition.

The upper border should be treated very carefully, as has been discussed by Green (1965). The fluid envelopes of stars and giant gaseous planets, which are considered as bounded fluid shells, are usually treated with a free-surface condition applied at the upper limit x = xatm (e.g. Unno et al. 1989), i.e. (123)Hence, (124)with the complex coefficients (125)The atmosphere then behaves as a wave-guide of typical thickness xatm (Appendix E). In this case, atmospheric tides are analogous to ocean tides (see Webb 1980), the amplitude of the perturbation being highly frequency-resonant. However, this condition appears to be inappropriate for the thin isothermal atmosphere because the fluid is not homogeneous and there is no discontinuous interface in this case. Setting the free-surface condition would give birth to unrealistic resonances depending on xatm (see Fig. E.1 in Appendix F). Therefore, we have to choose a condition consistent with the exponentially decaying density and pressure. Following Shen & Zhang (1990), we consider that there is no material escape when x → + ∞, i.e. Vr|x → ∞ = 0. This means that the amplitude of oscillations decrease with the altitude and, consequently, that we have to eliminate the diverging term in Eq. (121)6.

Let us assume . As shown by CL70 (Eq. (88)), the above condition requires that (126)at the upper boundary (x = xatm). To make this condition relevant, xatm must be chosen so that xatm/xcrit ≫ 1, where xcrit is the typical damping depth computed from the case where C is a constant, and given by (127)we usually fix xatm = 10 in computations. The integration constant of Eq. (122) is then obtained as (128)

3.4. Analytical expressions of the semidiurnal tidal torque, lag angle, and amplitude of the pressure bulge

If and are assumed to be constant (with ), C does not vary with x and the expression of Eq. (121) can be simplified significantly. Indeed, in this case it can be written (129)We note that we use the convention , which implies at infinity. Substituting this simplified solution in polarization relations (Eqs. (110)to (118)), we compute three-dimensional variations of the velocity field, displacement, pressure, density, and temperature caused by the semidiurnal tide as explicit functions of the physical parameters of the atmosphere and tidal frequency. Let us introduce the function of the altitude, parametrized by (which can be , n, or ), and defined by (130)Then, polarization relations write

with (140)From the analytical solution and the polarization relation of density given by Eq. (138), we then compute the second-order Love number and the tidal torque. Integrating δρ over the whole thickness of the atmosphere, we obtain (141)where the expressions of n and are given by Eqs. (120) and (100), respectively. The semidiurnal tide is supposed to be represented by spatial forcings of the form and with the tidal frequency σ = 2(Ω−norb). The second-order Love number, deduced from Eq. (75), is thus expressed (142)Assuming |εs;n| ≪ 1, we identify the different contributions. It follows (143)where the subscripts therm and grav indicate the origin of a contribution (thermally or gravitationally forced) and H and V its horizontal and vertical components. The terms of this expansion are expressed as functions of the parameters of the atmosphere and tidal frequency

the parameter being a dimensionless constant given by (148)Contrary to gravitational Love numbers ( and ) which are intrinsic to the planet, thermal Love numbers ( and ) are proportional to the forcings ratio J2/U2 and thus depend on the properties of the whole star-planet system, particularly the semi-major axis and the stellar luminosity (e.g. Correia et al. 2008).

In the same way as we obtained Love numbers, the torque exerted on the atmosphere may be computed by substituting Eq. (141) in Eq. (80). This torque writes (149)Finally, assuming |εs;n| ≪ 1 and introducing the factor (150)allows us to write under the form (151)with

Let us consider the case of the pure thermal tide (U2 ≈ 0) at the vicinity of synchronization. If σ → 0, then , |ℬn| → + ∞, and . As a consequence, and . The tidal torque exerted on the stably stratified isothermal atmosphere is thus very small compared to the horizontal thermal component,(156)where we have assumed that the thermal forcing is in phase with the perturber (J2 ∈R). This behaviour had been identified before, for instance by Arras & Socrates (2010) who studied thermal tides in hot Jupiters. According to Arras & Socrates (2010), at the vicinity of synchronization, the vertical displacement of fluid due to the restoring effect of the Archimedean force (N2Vr term in the heat transport equation, Eq. (14)) exactly compensates the local density variations generated by the thermal forcing, which annihilates the quadrupolar tidal torque. We discuss this effect in Sect. 6 when considering the isothermal and stably stratified atmospheres of the Earth and a Venus-like planet (see Fig. 14). For a stably stratified structure, the tidal torque exerted on the atmosphere is weak and cannot balance the solid torque, which leads the planet’s spin to the synchronization configuration. We note here that the importance of stable stratification has been pointed out for the case of Jupiter-like planets by Ioannou & Lindzen (1993b) (see also Ioannou & Lindzen 1993a, 1994), who showed that the atmospheric tidal response was very sensitive to variations in the Brunt-Väisälä frequency (N).

3.5. Comparison with CL70

The vertical wavenumber is a key parameter of the tidal perturbation. Hence, to highlight the interest of taking into account dissipative mechanisms such as radiation, we consider the relative difference between the given by Eq. (103) and the value established via the classical theory of tides without dissipation (e.g. Wilkes 1949, CL70), denoted in reference to CL70 and given by (157)

thumbnail Fig. 6

Relative difference between the complex wavenumber obtained in this work (Eq. (103)) on the one hand, and the real one given by CL70 for Earth atmospheric tides () on the other hand. The difference η, given by Eq. (158), is plotted in logarithmic scale as a function of the normalized tidal frequency σ/σ0 and of the ratio hc/hn. The colour of the map, denoted c, is defined by the expression c = log (η). The horizontal and vertical axis represent f(σ/σ0) and , where f is the function defined by with X ∈R. The positions of the planets are determined by the gravity mode of degree 0 in the case of the semidiurnal tide. Dark (luminous) regions correspond to the dynamic (dissipative) regime. The vertical wavenumbers of Venus are strongly affected by dissipative processes, while those of the Earth are not.

Open with DEXTER

This difference, written (158)is plotted in Fig. 6 as a function of the reduced tidal frequency σ/σ0, and of the height ratio hc/hn. The supplementary term of Newtonian cooling in the heat transport equation induces an additional regime with respect to CL70. At low frequencies, around the radiation frequency, the vertical profiles are damped, while they are highly oscillating in the case where σ0 = 0. This thermal regime corresponds to the middle white-yellow region of the map. The two wavenumbers become similar for |σ| ≫ σ0. This is the regime of fast rotating planets studied by CL70. Beyond the threshold corresponding to |σ| ~ N, the traditional approximation assumed in Sect. 2 is no longer valid and the coupling between the three components of the Navier-Stokes equation is considered. Finally, the area located at hc/hn ≈ 1 represents the discontinuous transition between propagative and evanescent waves in CL70, this transition being regular in our model.

4. Tides in a slowly rotating convective atmosphere

In this section, we treat the asymptotic case of a slowly rotating convective atmosphere for which the equations of tidal hydrodynamics can be simplified drastically. This situation typically corresponds to Venus-like planets. Indeed, the rotation period of Venus, 243 days, is of the same order of magnitude as its orbital period. It follows that |ν| < 1 for the semidiurnal tide of Venus. As a consequence, the Coriolis acceleration can be neglected as the first step. The atmospheric layer located below 60 km is characterized by a strongly negative temperature gradient (Seiff et al. 1980). In this layer, the temperature decreases from 750 K to 250 K, which make it subject to convective instability (Baker et al. 2000). Therefore, the Brunt-Väisälä frequency of Venus is far lower than the value given by the isothermal approximation. As this frequency represents the “stiffness” of the Archimedean force which restores gravity waves, these waves cannot propagate in the unstable region above the surface.

4.1. Equilibrium distributions of pressure, density, and temperature

We assume Ω = 0 and N = 0. According to Eq. (1), the Brunt-Väisälä frequency can be expressed as(159)Therefore, the case N2 = 0 corresponds to the adiabatic temperature gradient, given by(160)For this structure, the reduced altitude is expressed by (161)and we thus obtain the distributions of the basic pressure height, temperature, pressure, and density, respectively (162)with .

4.2. Tidal response

For the sake of simplicity, we neglect the contribution of the gravitational forcing in the Navier-Stokes equation and examine the case of pure thermal tides. When σ → 0, the tidal response tends to its hydrostatic component, the equilibrium tide (Zahn 1966). We thus assume the hydrostatic approximation. This allows us to reduce the radial projection of the Navier-Stokes equation (Eq. (8)) to(163)Now, by expanding the perturbed quantities in Fourier series and substituting Eq. (163) in the heat transport equation (Eq. (20)), we get the simplified vertical structure equation of pressure (164)represents the complex damping rate of the perturbation. Following Dobrovolskis & Ingersoll (1980) and Shen & Zhang (1990), we consider that thermal tides are generated by an increase in heat at the ground and apply the profile of thermal forcing (165)where 1 /τJ ≪ 1 is the thickness of the heated region (see Eq. (213)for a physical expression of this parameter) and J0 the mean absorbed power per unit mass. At the upper boundary xatm ≫ 1, we apply the free-surface condition δp = 0 and get the pressure profile(166)from which we deduce the vertical profiles of the Lagrangian horizontal displacements (167)(168)of the horizontal velocities(169)(170)and of the density and temperature (171)(172)

4.3. Second-order Love number and tidal torque

The previous results enable us to compute the tidal Love numbers and torque associated with the atmospheric tidal response. Owing to the hydrostatic approximation, these quantities are directly proportional to δp(x = 0). We consider that the semidiurnal thermal tide can be reduced to its quadrupolar component, given by and . Hence, using the expressions Eq. (76) for the second-order Love number and Eq. (82) for the tidal torque, we obtain (173)and(174)

4.4. Comparison with previous models

The tidal torque plotted in Fig. 7 is identical to the horizontal component of the torque exerted on the stably stratified isothermal atmosphere given by Eq. (156) because convection eliminates the component resulting from the fluid vertical displacement. In agreement with the early result of Ingersoll & Dobrovolskis (1978) (Eq. (4)), it is of the same form as the result given by the Maxwell model (Correia et al. 2014), where we identify the Maxwell time . It also corresponds to the torque computed by Leconte et al. (2015) for Venus-like planets with numerical simulations using GCM, which suggests that the slow rotation and convective instability approximations are appropriate for planets of this kind. This torque is stronger than that computed in the case of the stably stratified atmosphere. Consequently, it can lead to non-synchronized rotation states of equilibrium by counterbalancing the solid tidal torque. Finally, it must be compared to the model introduced by Correia & Laskar (2001) in early theoretical studies of the tidal torque exerted on Venus atmosphere. In this model, the tidal torque is given by (175)where a and b are two empirical positive real parameters. As can be observed in Fig 7 where (174) and (175) are plotted as functions of the tidal frequency, looks like from Eq. (174) with amplitude maxima located at (176)This strengthens the statement of Correia & Laskar (2001), who expected that “further studies of atmospheres of extra solar synchronous planets [ ...] may provide an accurate solution for the case σ ≈ 0”. In Fig. 14, the expression of Eq. (175) is plotted for the Earth and Venus in addition to the analytical results of the model presented here with parameters a and b adjusted numerically.

thumbnail Fig. 7

Tidal torques exerted on the atmosphere in the case of the slowly rotating convective atmosphere (red line with the label “Conv”, Eq. (174)) and in the model by Correia & Laskar (2001) (green dashed line with the label “CL01”, Eq. (175)) as functions of the tidal frequency. Both functions are normalized by their respective maximum value, reached at σ = σ0; the tidal frequency is normalized by the Newtonian cooling frequency (σ0). The tidal forcing is quadrupolar, with U and J given by Sect. 4.3. The parameters a and b of Eq. (175) are adjusted so that the peaks of the two functions are superposed.

Open with DEXTER

5. Physical description of heat source terms

On the left-hand side of Eq. (46) and in polarization relations, the perturbation is driven by and . The tidal gravitational potential has been expanded in spherical harmonics for many years and there is no need to come back to it7. Nevertheless, it is necessary to establish physical expressions of the heat power per unit mass to compute a solution. The different contributions must be clearly identified. We detail them in the case of the optically thin atmosphere. The most obvious heat source, but not necessarily the one with the largest amplitude, is the absorption of the light flux coming from the star. The non-absorbed part of this flux reaches the ground and causes temperature oscillations at x = 0. This implies two contributions from the ground: a radiative emission in infrared and diffusive heating due to turbulence within the surface boundary layer.

5.1. Insolation

The heating fluxes coming from the ground are all oriented radially and proportional to the local insolation flux at x = 0. Therefore, they can be easily written as a product of separable functions: . This is not the case of the insolation flux, which goes through the spherical shell along the star-planet direction (see Fig. 10). Particularly, the atmosphere is partly enlightened by stellar rays on the dark side. We note IS, the bolometric flux transmitted to the atmosphere. If we assume that the star radiates like a black body, this flux is given by the integral of the spectral radiance IS;λ over a wide range of wavelengths (177)where IS;λ is given by Plancks’s law (Planck 1901): (178)the parameter c being the speed of light in vacuum, h the Planck constant, kB the Boltzmann constant, d the star-planet distance, RS the radius of the star, and TS its surface temperature. The flux is partly reflected by the atmosphere with the albedo Aatm;λ for the wavelength λ, which gives the effective heating flux (179)To describe the absorption of I0 by the atmosphere, we use the Beer-Lambert law (Bouguer 1729; Klett et al. 1760; Beer 1852). Light rays propagate along a straight line defined by the linear spatial coordinate l. So, the absorption of the incident spectral density by a gas G of molar concentration CG, illustrated by Fig. 8, can be written (180)the parameter εG;λ being the molar extinction density coefficient of the gas, which is supposed to depend on the light wavelength (λ) alone. This dependence should be taken into account because εG;λ can vary over several orders of magnitude, typically, εG;λλ-4 in the regime of Rayleigh diffusion (molecules are small compared to the wavelength). From Eq. (180), we get (181)and the incident flux (182)

thumbnail Fig. 8

Left: absorption of the flux emitted by the star along the star-planet direction. Right: absorption of the flux emitted by the ground.

Open with DEXTER

Then, performing the derivation of Eq. (181), the heat power per unit volume is obtained, (183)and the heat power per unit mass comes straightforwardly (184)We retrieve the formulation of the thermal forcing proposed by CL70. Considering that the gas is homogeneous in composition leads us to CG = C0 = ρ0/M. Therefore, the incident power flux and heat power per unit mass become (185)Near the planet-star axis, they are simply expressed as (186)where τλ is the dimensionless optical depth for the radiance of wavelength λ, (187)The parameter 1 /τλ corresponds to an absorption depth of the flux in a homogeneous fluid normalized by the pressure height scale (H). If τλ ≪ 1, the heating power of a wavelength λ is almost totally transmitted to the ground. On the contrary, if τλ ≫ 1, the flux does not reach the surface because it is entirely absorbed by the atmosphere (see Fig. 9). In this case, the atmosphere is not optically thin, and it cannot be assumed that its radiative losses are proportional to δT as was done in Sect. 2, and thermal diffusion has to be considered. Therefore, the case treated by the present work corresponds to τλ ≪ 1.

thumbnail Fig. 9

Vertical profiles of radiative fluxes (I), power per unit volume (J), and power per unit mass (J) for various values of the dimensionless optical depth τλ (Eq. (187)). As usual, the reduced altitude, x = z/H, is on the vertical axis. Left: incident flux, defined by Eq. (186), at ψ = 0. If τλ ≪ 1 the ground receives almost 100% of the total power. On the contrary, if τλ ≫ 1, all the power is absorbed by the atmosphere. Right: radiative flux emitted by the ground, defined by Eq. (206). For small τλ, the flux goes through the atmosphere without being absorbed. For high values of the parameter, it is absorbed by the lowest layers.

Open with DEXTER

The proportion of flux density transmitted to the ground αλ is deduced from Eq. (186) taken at x = 0, (188)Like the atmosphere, the ground has an albedo denoted Agr;λ, which implies that the upcoming flux is partly reflected. Therefore, we can write the absorption equation (189)and compute the reflected surface power (190)the coordinate ψ being the angle between the star-planet direction and the position vector of the current point. The heat power per unit mass and wavelength provided by the reflected flux is derived from the previous equation (191)So, the contribution of the reflected flux is finally obtained (192)

thumbnail Fig. 10

Left, from top to bottom: power per unit mass (Jinc, Eq. (184)) and power per unit volume (Jinc, Eq. 183) of a monochromatic incident flux normalized by their maxima and plotted as functions of the normalized spatial coordinates X (the direction of the star is along the horizontal axis) and Z (the spin axis is along the vertical axis). The planet is represented by a green disk of unit radius. The flux, of wavelength λ, propagates in an atmosphere assumed to be homogeneous in composition (CG = C0) and characterized by the ratio H/R = 0.1. The optical depth is given by τλ = 0.2. Right, from top to bottom: power per mass unit (Jrad, Eq. (206)) and power per volume unit (Jrad) of a monochromatic flux radiated by the ground, normalized by their maxima, and plotted as functions of X and Z, with τλ = 0.5.

Open with DEXTER

5.2. Radiative heating from the ground

The part of the incident flux which is not reflected causes surface temperature oscillations. This effect was studied in meteorological works throughout the twentieth century. We will follow here the theoretical approach proposed by Bernard (1962) for its simplicity and the physical landmarks that it brings. The ground is considered to be a black body of temperature Tgr, emitting the power flux Igr;BB given by the Stefan-Boltzmann law, (193)the parameter ϵgr being the emissivity of the ground and S= 5.670373 × 10-8 W m-2 K-4 the Stefan-Boltzmann constant (Mohr et al. 2012). Moreover, the atmosphere emits a counter radiation to the surface, which implies that the effective flux is less than Igr;BB. According to Bernard (1962), this counter-radiation can be assumed to be proportional to Igr;BB, and expressed by the semi-empirical formula (194)which allows us to take it into account without studying the whole coupled system ground-atmosphere. The factor ϵatm corresponds to an effective emissivity. So, introducing the spectral radiance of the atmosphere Iatm;λ, the albedo of the ground in the infrared Agr;IR can be defined by the relationship (195)and the effective flux emitted by the telluric surface can be written (196)Given that , Eq. (196) becomes (197)where ϵ = ϵgr(1−ϵatm) ≤ 1 is the effective emissivity of the ground. The terrestrial atmosphere being rather thin and transparent, ϵ ≈ 1 for the Earth. This expression will be used further to establish the heating by turbulent diffusion. We also introduce here the corresponding spectral radiance of the ground, (198)The radiative flux emitted by the ground, denoted Irad, is the solution of the absorption equation (similar to the one used to compute Iinc), (199)with the boundary condition at x = 0. Therefore, (200)and finally, (201)The temperature of the ground Tgr depends on the power reaching the surface, (202)Both quantities are linearized near the equilibrium, and Tgr = Tgr;0 + δTgr, and the perturbation is expanded in series (203)The spatial functions and are related by a transfer function of the tidal frequency, denoted Bgr, that is given explicitly in the next subsection (204)So, the perturbation of the spectral radiance emitted by the ground can be written (205)and the flux and heat power per mass unit (206)The partial derivative of is explicitly given by (207)Here, we note that the response of the ground to the tidal forcing induces a dependence on the tidal frequency through the oscillations of Tgr and the transfer function . Thus, the spatial functions describing the contribution of the ground are parametrized by σ contrary to those associated with the incident flux.

5.3. Boundary layer turbulent heat diffusion

Near the ground, heat transfers are dominated by turbulent mechanisms and diffusion. This is known as the planetary boundary layer and it is illustrated in Fig. 12. Thus, for x ≪ 1, thermal diffusion cannot be ignored as it is in Sect. 2 and drives the oscillations of temperature. However, we note that it is only valid at the vicinity of the surface, where friction with the ground is important. After CL70 and Siebert (1961), we neglect as a first step the heat transport through the troposphere due to advection to provide a first analytical treatment of the problem. At the interface, the diffusive fluxes in the ground and in the atmosphere write (208)the parameters kgr and katm representing the thermal conductivities of the ground and of the lowest layer of the atmosphere. In addition, denoting the associated thermal capacities cgr and catm, we introduce the thermal diffusivities (209)the parameter Katm being the vertical turbulent thermal diffusivity. Therefore, temperature variations near the surface are described by the equations of heat transport (e.g. CL70), (210)where only the diffusive terms are taken into account. The power balance of the surface is deduced from Eqs. (208) and (197), (211)and provides the temperature Tgr;0 of the ground at equilibrium. Linearized with variables of the form given by Eq. (21), Eq. (210) can be written (212)where and are the tidal frequency-dependent skin thicknesses of diffusive heat transport in the solid and fluid parts, respectively. Their expressions are given by (213)It should be noted that and both decay when the tidal frequency increases. Since corresponds to the typical depth of the boundary layer, it has to satisfy , otherwise the diffusive effects cannot be ignored in the dynamical core of the model (Sect. 2). Let us denote the radiative impedance of the ground considered as a perfect black body and submitted to a perturbation in temperature. The transfer function introduced in Eq. (204) can be deduced from the linearized power balance, (214)in which the expressions of Eq. (212) are substituted. We finally get (see Fig. 11) (215)

thumbnail Fig. 11

Nyquist diagram of the boundary layer transfer function in the frequency range − 103σ/σbl ≤ 103 (σ being the tidal frequency and σbl the frequency characterizing the boundary layer introduced in Eq. (217)). The parametrized complex transfer function (Eq. (215)) reduced by its value at σ = 0 is plotted as a function of the tidal frequency (σ) in the complex plane, where the horizontal axis corresponds to the real part of the function, , and the vertical axis to the imaginary part . The arrow indicates the direction of increasing σ. The gain of the transfert function is given by the distance of a point P of the curve to the origin O of coordinates , and its phase by the angle of the vector OP to the horizontal axis. This plot shows that the response of the ground is always delayed with respect to the thermal forcing except at synchronization, where and the gain is maximum. It also highlights the two identified regimes of the ground response: the lag tends to 0 (response in phase with the perturbation) and the gain increases at σ → 0, while the lag tends to ± π/ 4 and the gain to 0 at σ → ± ∞, the critical transition frequency being σbl.

Open with DEXTER

an expression that can be decomposed in its modulus and argument (216)The frequency σbl, which parametrizes , is a characteristic frequency reflecting the thermal properties of the diffusive boundary layer, and is expressed (217)where the parameters βgr and βatm are the thermal conductive capacities of the ground and of the atmosphere near the telluric surface (218)According to Bernard (1962), turbulent diffusion plays a more important role than thermal conduction in the ground, which means βatmβgr. However, the turbulent diffusivity Katm may vary over a wide range of orders of magnitude, taking extremal values above oceans (typically βatm ~ 104 J m-2 K-1 s− 1/2). For these reasons, Bernard (1962) prescribes for the Earth the effective mean value Katm ~ 10 m2 s-1, which had been prescribed before by Wilkes (1949) and which leads to σbl ~ 10-6 s-1.

thumbnail Fig. 12

Bichromatic case. Top: heat transfers between the atmosphere, space and the ground. Bottom: power balance of the interface ground-atmosphere.

Open with DEXTER

As shown by Eq. (216), the frequency ratio |σ| /σbl determines the angular delay of the ground temperature variations. If |σ| ≪ σbl, the response is in phase with the excitation. It corresponds to low conductive capacities. On the contrary, if |σ| ≫ σbl, then the delay tends to its asymptotical limit, π/ 4. By using the frequency σbl ~ 10-6 s-1 derived from the prescriptions of Wilkes (1949) and Bernard (1962), we note that the Earth’s solar diurnal tide belongs to the second category. This explains the position of the diurnal peak observed in surface temperature oscillations, three hours late with respect to the subsolar (midday) point. The gain of is damped by the conductive capacities of the material at the interface. Therefore, a very diffusive interface tends to attenuate the temperature oscillations of the ground, as expected.

We now consider the profile of temperature variations in the boundary layer obtained in Eq. (212) (219)where is the vertical damping ratio of the diffusive perturbation. Chapman & Lindzen (1970) shows that a heat power due to diffusion Jdiff can be deduced from the previous expression. Indeed, by writing the heat transport equation as (220)and linearizing it, we compute the expression of the tidal diffusive forcing term as (221)

5.4. Simplified bichromatic case

The expressions of the thermal forcings obtained in the previous section may appear of some complexity because of the dependence of albedos and optical depths on the wavelength. They can be simplified notably if each source is assumed to be monochromatic. In this case, the star emits in the frequency range of visible light, and the planet in the range of infrared. Consequently, the wavelength-dependent Agr;λ, Aatm;λ, ελ, ϵλ, αλ, and τλ are reduced to the constant coefficients Agr;V, Aatm;V, εV, ϵV, αV, and τV in the range of visible light and Agr;IR, Aatm;IR, εIR, ϵIR, αIR, and τIR in the infrared. Hence, the bolometric flux is given by (222)where TS is the temperature of the star. This flux is partly reflected back to space by the atmosphere so that the effective flux I0 getting through the layer can be written (223)The power of the incident flux Iinc decays along the path of a light beam, absorbed by the gas. The power absorbed per unit mass is denoted Jinc and the expressions of Eq. (185) become (224)From Iinc, we deduce the flux reflected by the ground Iref and the power absorbed per mass unit Jref, given by (225)Finally, the expressions of Eq. (206) are simply replaced by (226)The expressions of the powers per unit mass of the incident and radiated flux (Jinc and Jrad), and the corresponding powers per volume unit (Jinc and Jrad), are plotted in Fig. 10. The heating due to the turbulent boundary layer remains unchanged (see Eq. (221)).

6. Application to the Earth and Venus

To illustrate the regimes identified and the dependence of the atmospheric response on the forcing frequency, we apply the model to the cases of the semidiurnal thermal tides on Earth and Venus assuming first for both planets isothermal stably stratified atmospheres, as described by the equations in Sect. 3. We note here that the atmosphere of Venus is not stably stratified but convective in regions close to the surface (Seiff et al. 1980; Baker et al. 2000) where the tidal mass redistribution is the most important (Dobrovolskis & Ingersoll 1980; Shen & Zhang 1990). Nevertheless, examining the case of a stably stratified isothermal atmosphere for a Venus-like planet will allow us to better understand in an academic framework the role of stratification and radiative losses in the tidal response. To make a clear difference between Venus and the studied stably stratified Venus-like planet, we call the latter “VenusX”. For the sake of simplicity, we assume that planets orbit circularly in their equatorial plane. Hence, the tidal frequency8 is given by σ = 2(Ω−norb). One of the most crucial parameters of the model is the Newtonian cooling frequency (σ0). This parameter varies over several orders of magnitude with the altitude (Pollack & Young 1975). However, our purpose here is not to compute a quantitative tidal perturbation, but to illustrate qualitatively the non-dissipative (Earth) and dissipative (Venus) regimes. Therefore, following Lindzen & McKenzie (1967), we consider the interesting case where the radial profile σ0 is assumed to be constant. For both planets, we arbitrarily set σ0 = 7.5 × 10-7 s-1, which is the effective value of σ0 computed by Leconte et al. (2015) with a GCM for Venus. This value is such that σ ~ σ0 for Venus and σσ0 for the Earth. In the case of an optically thin atmosphere, the Newtonian cooling frequency can be estimated using Eq. (17). Moreover, we impose academic quadrupolar perturbations of the form (227)where U2 and J2 are fixed constant radial profiles. The tidal potential U2 is computed from the Kaula multipole expansion detailed in Appendix D, and explicitly given by the expression of Eq. (74). A zero-order approximation of the thermal power J2 can be quantified by writing the total surface power absorbed by the atmosphere as a function of θ and ϕ, which is then expanded in Fourier series of longitude and associated Legendre polynomials. Given that atmospheric tides are considered in this model as linear perturbations around the equilibrium state, the amplitudes of the perturbed quantities are proportional to the forcing.

thumbnail Fig. 13

Average daily variation in temperature in degrees (top) and pressure in mbar (bottom) over one year, plotted over solar time (in hours), the subsolar point corresponding to 12 h. The data was recorded every minute over the full year 2013 with a Vantage Pro 2 weather station at latitude 48.363°N. The time is converted to solar time and the data averaged over an exact count of 365 days. The average values of temperature (10.8973° C) and pressure (1015.83 mbar) have been removed to display only the variable part. In both figures, the raw data (in green) has been decomposed in a Fourier series over the 24 h period, limited to two harmonics (blue), or five harmonics (red). Temperature is dominated by the diurnal component, with a maximum value around 14h26mn. The largest harmonics of the pressure variations is the semidiurnal term, followed by the diurnal term. Two pressure maxima can be observed, at 9h38mn and 22h07mn.

Open with DEXTER

thumbnail Fig. 14

Tidal torque applied to the atmosphere of Earth (left panel) and Venus (right panel) by the thermal semidiurnal tide and its components as functions of the reduced frequency χ = (Ω−norb) /norb (horizontal axis). The forcings consist of academic quadrupolar perturbations expressed as and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3 for the Earth and U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3 for Venus. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively (we note that ). The label CL01 refers to the model by Correia & Laskar (2001), given by Eq. (175).

Open with DEXTER

The horizontal structure equation is solved using the spectral method described in Appendix B. The vertical structure equation is integrated numerically on the domain x ∈ [0,xatm] using a regular mesh with element of size δx. As pointed out by CL70, the number of points in the mesh, denoted N, must be sufficiently large to obtain vertical profiles of the perturbation with a good accuracy (see Appendix G). To fix N correctly in the case of the Earth, CL70 suggests using a criterion that we adapt for other configurations, i.e. (228)considering that the scale of the vertical variations is defined by the module of the complex wavenumber given in Eq. (107).

For each planet, we give here the spatial distribution of perturbed quantities and the evolution of the tidal torque with the tidal frequency. The numerical values used in these simulations are summarized in Table 1.

Table 1

Values of physical parameters used in simulations.

thumbnail Fig. 15

Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3. Left, from top to bottom: amplitudes of δp, δρ, δT, Vθ, Vϕ, and Vr in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (|δf|) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)].

Open with DEXTER

thumbnail Fig. 16

Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3. Left, from top to bottom: real parts of δp, δρ, δT, Vθ, and Vϕ taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}.

Open with DEXTER

thumbnail Fig. 17

Tidal surface oscillations of pressure, density, temperature, and their components at the equator of the Earth as functions of the reduced tidal frequency χ = (Ω−norb) /norb (horizontal axis). The perturbation is defined by the academic quadrupolar forcings and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively.

Open with DEXTER

6.1. The Earth

The Earth’s semidiurnal thermal tide corresponds to one of the cases detailed in CL70. Here, the radiation frequency is clearly negligible compared to the tidal frequency σ ≈ 2Ω. Therefore, the tidal response of the Earth’s atmosphere is driven by dynamical effects only and the radiative losses added in the heat transport equation can be ignored. This case is typical of a pure dynamic behaviour (see Figs. 5 and 6). Given that ν is slightly larger than 1, the horizontal structure of the perturbations is essentially described by gravity modes (Fig. 2). Rossby modes are trapped at the poles. Figures 15 and 16 show maps of the perturbed quantities in latitude (denoted δ), longitude (ϕ), and altitude (x), the subsolar point being indicated by the coordinates (δ,ϕ) = (0,0). In these figures, we can observe that the pressure and density variations present very similar spatial distributions. Particularly, they are in phase with each other and the phase lag is approximately the same in the whole layer, which causes the tidal bulge to materialize. The lag angle γ = π/ 4 of the bulge with respect to the subsolar point is given by Eq. (137) taken at x = 0 for the gravity mode of lowest degree (n = 0) and is not very sensitive to temperature structure (Lindzen 1968). For the Earth, these semidiurnal pressure peaks are measured around 9h38mn and 22h07mn (see Fig. 13).

The patterns of the velocities and are also explained well by the first horizontal Hough functions, and (see Fig. 2, middle and bottom left). In contrast, the behaviour of the temperature involves other modes.

We then compute the evolution of the tidal torque with the tidal frequency by changing the value of the spin frequency (Ω) in simulations. Hence, the reduced tidal frequency χ = (Ω−norb) /norb varies within the interval [− 15,15]. The torque is computed using the analytical formulae of Eqs. (141) and (149). We also compute the horizontal and vertical components isolated in Sect. 3. The results are plotted in Fig. 14.

thumbnail Fig. 18

Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3. Left, from top to bottom: amplitudes of δp, δρ, δT, Vθ, Vϕ, and Vr in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (|δf|) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)].

Open with DEXTER

thumbnail Fig. 19

Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3. Left, from top to bottom: real parts of δp, δρ, δT, Vθ, and Vϕ taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}.

Open with DEXTER

thumbnail Fig. 20

Tidal surface oscillations of pressure, density, temperature, and their components at the equator of Venus as functions of the reduced tidal frequency χ = (Ω−norb) /norb (horizontal axis). The perturbation is defined by the academic quadrupolar forcing and with U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively.

Open with DEXTER

At first, we note that the horizontal and vertical components of the torque are of the same order of magnitude and in opposition to one another, as shown by Eq. (141) taken at the limit σ → 0. This behaviour was identified previously in Sect. 3.4. It results from the fact that vertical displacement of fluid induced by the stable stratification compensates the local density variations in this frequency regime and strongly flattens the total tidal torque. As demonstrated in Sect. 4, this is not the case in general and it depends on the strength of the stratification. In the case of a slowly rotating convective atmosphere (N2 ≈ 0) forced by a heating located at the ground, only the horizontal component remains. This leads to a much stronger tidal torque, similar to those obtained by Leconte et al. (2015) with a GCM (given by Eq. (174)) and Correia & Laskar (2001) (see Eq. (175)).

In Fig. 17, we plot the surface variations of pressure, density, temperature, and their components as functions of the tidal frequency. For each quantity, the horizontal component varies smoothly while the vertical one is discontinuous at the synchronization. This discontinuity is a consequence of the dissymmetry between prograde and retrograde Hough modes due to the Coriolis acceleration (Fig. A.1), particularly as regards the gravity mode of the lowest degree, which is the most important one.

Assuming that the tidal torque is entirely transmitted to the telluric core of the planet, it is possible to estimate the timescale of the spin evolution induced by the atmospheric semidiurnal tide. Introducing the moment of inertia of the planet, , we write (229)Thus, the timescale of the spin evolution is given by (230)For the Earth’s semidiurnal tide, χ ≈ 365 and the horizontal component of the corresponding torque is . With ℐ = 8.02 × 1037 kg m2 (NASA fact sheets10), T ~ 180 Gyr. Therefore, the spin frequency of the Earth is not affected by atmospheric tides.

6.2. Venus

Geometrically speaking, Venus can be seen as an Earth-like planet, with scales of the same magnitude order for the telluric core and the thickness of the atmosphere. However, there are some important differences compared to the case of the Earth. The Venusian atmosphere does not have the same properties as the Earth’s atmosphere. First, the fluid layer is a hundred times more massive (denoting M the mass of the atmosphere, M ~4.8 × 1020 kg and M ~ 5.1 × 1018 kg). Thus, the surface pressure is around 93 bar. Second, it is optically thicker. The major part of the heating flux is thus deposited at high altitudes. The absorbed flux near the planet surface at the subsolar point Fabs is estimated at Fabs ~ 100−200 W m-2 (Avduevsky et al. 1970; Lacis 1975; Dobrovolskis & Ingersoll 1980), which implies J2 ~ 10-4 W kg-1. Finally, layers below 60 km are characterized by a strongly negative temperature gradient (Seiff et al. 1980) and are, therefore, weakly stratified or convective. As discussed in Sects. 3 and 4, these properties have a strong impact on the nature of the tidal response. Particularly, the tidal torque may be much stronger in the convective case than in the stably stratified one. In the present section, we study a Venus-like planet with a stably stratified isothermal atmosphere (VenusX) in order to understand the effects of stratification and radiative losses in a general context.

The other important difference between the Earth and Venus is in their rotational dynamics. While the spin frequency of the Earth is far higherthan its mean motion and the radiation frequency of its atmosphere, these three frequencies are of the same order of magnitude for Venus. Therefore, the Venusian semidiurnal thermal tide is characterized by a tidal frequency |σ| ~ σ0, which is typical of the thermal regime. The term describing radiative losses in the heat transport equation (Eq. (18)) plays an important role in this regime by damping the profiles in altitude of the tidal response; see Fig. 5. This figure shows that the imaginary part of , which is responsible for the damping, is always comparable to the real part in this frequency range.

Contrary to the case of the Earth, the semidiurnal thermal tide of VenusX belongs to the family of super-inertial waves, with ν ≈ 0.48, which means that the horizontal component of the tidal response is composed of gravity modes only. We note that the surface variations of perturbed quantities plotted in Fig. 19 are very well represented by the first mode. The patterns of the Hough functions , , and (see Fig. 3, red curves) clearly appear on the maps. Moreover, the lags of δp and δρ are different from that observed in the Earth’s semidiurnal tide. In particular, it is interesting to note that the pressure and density peaks are not superposed in this case. The effect of the damping caused by radiative losses can be observed in vertical cross-sections in Fig. 18, where the variations of pressure and density are located close to the ground. Given that the real and imaginary part of the vertical wavenumber are both very high, vertical component is both strongly oscillating and strongly damped.

As was previously done for the Earth, we draw the variations of the tidal torque (Fig. 14) and perturbed quantities at the ground (Fig. 20) around the synchronization for χ ∈ [−15,15]. The behaviour of the torque and its components is the same as previously observed for the Earth. The maxima of horizontal and vertical components are located at χ = ± χ with χ<χ owing to the difference of mean motions between the two planets. Since the tidal torque is proportional to the surface density of the atmosphere, it is a hundred times stronger in the case of Venus than in the case of the Earth for a given thermal forcing. This is also verified by the surface variations of the perturbed quantities, which – except for their amplitudes – are comparable to those of the Earth.

To conclude this section, we compute the timescale of the spin evolution. The Venusian semidiurnal tide is characterized by the frequency χ ≈ −1.92. With the moment of inertia ℐ = 5.88 × 1037 kg m2 (NASA fact sheets) and the horizontal component of the torque given by Fig. 14, i.e. , we obtain T ~ 10 Myr. The thermal forcing corresponds approximately to the mean level of the heating profile used by Dobrovolskis & Ingersoll (1980) (J2 ~ 10-4 W kg-1) so the obtained torque is of the same order of magnitude as the value given by this early work, . We note here that the amplitude of the response is directly proportional to the amplitude of the excitation if one of its two components, J or U, can be neglected (Sect. 2).

6.3. Discussion of physical assumptions

We have developed a new analytic model for atmospheric tides by generalizing the classical tidal theory detailed in CL70 to dissipative atmospheres.

This work leans on several assumptions which contribute to make computations simpler. The most important of them are listed here.

  • The traditional approximation. This hypothesis, com-monly used in the literature of geophysics, is very usefulbecause it allows us to separate the horizontal and verti-cal dependences in the Navier-Stokes equation. Hence,the horizontal and vertical structure can be computedseparately. This approximation is relevant in stronglystably stratified fluids where 2Ω < |σ| ≪ N, but is less appropriate if|σ| ~ 2Ω, as discussed in Mathis (2009). Its validity domain, not well defined, remains essentially limited to the regime of super-inertial waves (Prat et al. 2016). It is represented in Fig. 21, which shows the cases that can be treated analytically (regions a and b). Concerning the importance of the neglected Coriolis terms, we refer the reader to Gerkema & Shrira (2005), who explored the physics of gravito-inertial waves without the traditional approximation (see also Tort & Dubos 2014). The traditional approximation appears to be a relatively strong hypothesis because it does not take into consideration the case of rotating convective atmospheres (N2 ≈ 0, e.g. Ogilvie & Lin 2004). Convection is an important feature of the Earth’s nearly adiabatic layer (altitudes below 2 km) as well as the Venus’ atmosphere (Linkin et al. 1986). It has strong implications on the propagations of waves: gravity waves cannot exist in convective layers, which may significantly modify the tidal response, particularly if the convective turnover time is close to the period of the tidal forcing. In the general case, the effect of convection on the structure of waves can be studied either numerically with GCM (Fig. 21, d) or analytically with simplified models (e.g. based on a Cartesian geometry). The only exception for a global analytical study of a convective atmosphere is the case treated in Sect. 4 where stable stratification and rotation are too weak to be neglected (Fig. 21, b).

  • The linearization around the equilibrium. In Sects. 2 and 3, the equations have been linearized around the equilibrium. This method is only relevant in cases where the tidal response of the atmosphere can be considered as a small perturbation with respect to the equilibrium state. This condition is verified for the Earth and Venus.

  • The solid rotation approximation. We assumed that the atmosphere rotates with the solid core uniformly. Thus, the atmospheric circulation (both meridional and zonal) is ignored. This approximation is justified by the fact that the most important pressure and density variations are concentrated in the densest layers, near the ground. Given that these layers are bound to the solid core by viscous friction, they rotate at the same spin frequency. However, tidal waves might interact with the mean flow in cases where σ ≈ 0 (Goldreich & Nicholson 1989) and it will be interesting to take into account the global atmospheric circulation in further works, as is done in GCM simulations.

  • The use of simplified atmospheric vertical structures. In order to lighten calculations, we considered the cases of uniform background temperature (isothermal stably stratified atmosphere) and temperature gradient (convective atmosphere), and supposed that some parameters, such as σ0, do not vary with altitude. These simplified structures are a first step towards more sophisticated ones. Owing to the importance of the stratification in the tidal response, it will be necessary to treat the cases of observed structures in future works, for instance those of the atmospheres of Earth, Venus, and Mars.

  • The optically thin layer hypothesis. Radiative processes are introduced in the dynamics of atmospheric tides through Newtonian cooling. The linear dependence of the corresponding radiative term on the temperature perturbation can be interpreted as an optically thin layer approximation. Indeed, this term corresponds to the thermal behaviour of an atmosphere that would emit grey body radiation towards space and the planet surface. Actually, the flux emitted is partly absorbed by the atmosphere itself, which induces a complex thermal coupling between the different layers of the atmosphere. This is true particularly for optically thick atmospheres. In these cases, the Newtonian cooling should be seen as a convenient approximation of thermal radiation for an analytical approach where the dissipative processes cannot be treated without approximations. This approach seems to be well validated by numerical simulations using GCM (Leconte et al. 2015).

thumbnail Fig. 21

Asymptotic domains of parameters where the analytical approach presented in this work can be applied (green) and where the numerical approach is necessary (red) as functions of the frequency ratios N2/σ2 (horizontal axis) and (2Ω)2/σ2 (vertical axis). Green regions: stably stratified atmospheres satisfying the hierarchy of frequencies (2Ω)2N2 necessary to apply the traditional approximation a); weakly stratified or convective atmospheres of slowly rotating planets b). Red regions: stably stratified atmospheres where the hierarchy of frequencies (2Ω)2N2 is not satisfied c); super-adiabatic rotating convective atmospheres d).

Open with DEXTER

thumbnail Fig. 22

Predictive aspects related to the dependence of the atmospheric tidal torque on stratification. The neutral convective case implies a strong torque, which can lead to non-synchronized states of equilibrium. The stably stratified case implies a weak torque, which lets the planet evolve towards spin-orbit synchronization.

Open with DEXTER

7. Conclusions

We have set the bases of a new theoretical framework to study the atmospheric tides of Earth-like exoplanets analytically. This work is complementary to numerical simulations. The main objective is to explain the qualitative behaviour of an atmosphere submitted to both thermal and gravitational tidal perturbations and to explore the domain of key physical parameters. The parameters highlighted by the model are

  • the characteristic frequencies of the system (e.g. inertia,Brunt-Väisälä, radiative, boundary layer, tidal frequencies);

  • the involved characteristic length scales (e.g. planetary radius, pressure height scale of the atmosphere, length scales of the vertical variations of tidal waves);

  • the physical parameters of the planet (e.g. the molar mass of the atmosphere, the adiabatic exponent, the surface gravity, the surface pressure).

Using the equations of hydrodynamics, we have adapted the theoretical approach detailed in CL70 to thick dissipative atmospheric shells. We computed the analytic expressions of the horizontal and vertical structure equations, of the perturbed quantities (pressure, density, temperature, velocity, displacement), and of Love numbers and tidal torques exerted on the atmosphere. These expressions have been used to identify the possible regimes of atmospheric tides: a non-dissipative Earth-like dynamical regime characterizing fast rotating bodies, and a dissipative Venus-like thermal regime for planets near synchronization. The behaviours involved in the two cases are very different. The dynamical regime corresponds to the case studied by CL70. In this regime, tidal waves are either oscillating or evanescent depending of the equivalent depth (hn) associated with the mode. In the radiative regime, they are oscillating and evanescent, and the perturbation can be strongly damped.

The tidal response of the atmosphere will also be affected by the nature of the tidal excitation. A predominating thermal forcing generates thermal tides. On the contrary, if the thermal energy brought to the fluid layer is negligible compared to the effect of the tidal gravitational potential, the perturbation belongs to the family of gravitational tides. This case has not been studied in this work, but the vertical structure equation and polarization relations of the model clearly attest to behaviour that is qualitatively different from that observed for thermal tides. Finally, the tidal response of the atmospheric shell is obviously affected by lower and upper boundary conditions. We showed for instance that a free-surface condition and a no-material-escape condition would not give the same vertical profiles of the perturbation. Similarly, the ground may affect temperature oscillations by introducing a delayed thermal forcing locally.

Finally, the results obtained for the tidal torque are in agreement with the GCM simulations by Leconte et al. (2015). This torque appears to be strongly dependent both on the tidal frequency and on the atmospheric structure (Fig. 22). In the case of a stably stratified atmosphere, the tidal bulge vanishes at the vicinity of synchronization. As a consequence, the tidal torque exerted on the atmosphere is very weak and cannot balance the solid tidal torque, which should lead the planet to synchronization. On the contrary, if the stable stratification is sufficiently weak (weakly stratified or convective atmosphere), the atmospheric tidal torque can be strong. In this case we recover the Maxwellian behaviour obtained numerically by Leconte et al. (2015) and described by the early model of Correia & Laskar (2001). This would lead planets to non-synchronized states. These results highlight the predictive interest of the modelling of atmospheric tides by showing that it is possible to constrain the atmospheric structure of terrestrial planets from their rotation states of equilibrium.


1

In a binary star-planet system where the planet orbits circularly in the equatorial plane defined by XE and YE at the orbital frequency norb, the semidiurnal tide corresponds to m = 2, σ = 2(Ω−norb), and ν = Ω/(Ω−norb).

2

The properties of Hough functions and the method used to compute them in this work are detailed in Appendix A.

3

To solve the vertical structure equation numerically in any case, we use the method described in Appendix B.

4

The parameters j, p, and q are the usual indexes of the Kaula’s expansion of the tidal gravitational potential; see Appendix C.

5

These coefficients are the same as those given by Eq. (60) to an H factor.

6

In the case of the semidiurnal tide treated by CL70 (no dissipative processes), this condition cannot be applied because . To address this particular situation, CL70 proposes applying a radiation condition at the upper boundary, i.e. considering that the energy is only propagating upwards and eliminating the terms of the solution corresponding to an energy propagating downwards. However, in the present work, when σ0> 0, , so we do not have to consider this case.

7

For instance, see the Kaula expansion of the gravitational tidal potential (Kaula 1964) detailed in Appendix C and Mathis & Le Poncin-Lafitte (2009).

8

This frequency corresponds to the term defined by (l,m,j,p,q) = (2,2,2,0,0) in the multipole expansion of the forcings detailed in Appendix C.

Acknowledgments

We thank the anonymous referee for comments that helped to improve the paper. We thank Richard Lindzen for the very interesting discussions we had about the analytic modelling of the Earth’s atmospheric tides. We are also very grateful to Umin Lee for his help in understanding how to compute Hough functions efficiently. Finally, we thank Jeremy Leconte, who shared with us his helpful approach to the physical modelling of super-Earth atmospheres from the point of view of numerical simulations. P. Auclair-Desrotour and S. Mathis acknowledge funding by the European Research Council through ERC grant SPIRE 647383. This work was also supported by the Programme National de Planétologie (CNRS/INSU) and CoRoT/Kepler and PLATO CNES grant at CEA-Saclay.

References

Appendix A: Hough functions

Hough functions, named for Hough’s initial work (Hough 1898), are the solutions of Eq. (35) integrated as an eigenvalue problem of with appropriate boundary conditions (Longuet-Higgins 1968; Chapman & Lindzen 1970; Lee & Saio 1997). It is possible to compute these functions in several ways. Here, we present two possible approaches: the finite differences method proposed in the Numerical recipes (Press et al. 1986) and detailed in Lee & Saio (1997), and a spectral method developed from the works of Hough (1898), Longuet-Higgins (1968), and Chapman & Lindzen (1970).

Appendix A.1: Finite differences method

The operator m,ν of Eq. (30) writes (A.1)where μ = cosθ. The new Laplacian tidal equation is not modified when μ is replaced by μ. Therefore, the domain over which this equation is integrated can be reduced to 0 ≤ μ ≤ 1. Like the associated Legendre polynomials, the solutions are either even or odd with respect to the equatorial plane μ = 0. The even solutions satisfy , the odd solutions . For the sake of simplicity, the parity of the subscript used in the notations corresponds to the parity of the associated function. The solutions belonging to one of these two families are all defined by the boundary condition applied at μ = 0: for even solutions and for odd solutions. The operator m,ν has two regular singular points, a permanent one at μ = 1 and another one at μ = 1/|ν| in the sub-inertial regime, characterized by |ν| > 1 and, therefore, depending on the spin parameter. Both are sources of numerical difficulties. To tackle the problem of the singular point at μ = 1, we introduce the variable y(μ) defined by Press et al. (1986) and Lee & Saio (1997):

thumbnail Fig. A.1

Top: eigenvalues () of the gravity and Rossby modes of lowest degrees (− 6 ≤ n ≤ 5) as functions of the spin parameter (ν). The horizontal axis represents the function of ν given by , and the vertical axis the function of given by . Bottom: projection coefficients of the normalized Legendre polynomial on normalized Hough functions associated with these modes as functions of the spin parameter (see Eq. (70)).

Open with DEXTER

(A.2)Laplace’s tidal equation becomes (A.3)where (A.4)The boundary condition applied at μ = 0 depends on the parity of the sought solution, i.e. dy/ dμ = 0 at the equator for even modes and y = 0 for odd modes. We also apply at μ = 0 the normalization condition y = 1 for even modes and dy/ dμ = 1 for odd modes. Since and are not defined at μ = 1, the boundary condition at the poles is applied at μ = 1−ϵ, where 0 <ϵ ≪ 1. This condition is derived from the Sturm-Liouville theory, which requires that should be regular at μ = 1, and is given explicitly by (see Lee & Saio 1997).

In order to solve Eq. (A.3) numerically, the domain is discretized. Hence, the equation is reduced to a linear system of dimension NS = 2N, where N is the number of elements of the mesh. In the regime of sub-inertial waves, the singular point at μ = 1/|ν| corresponds to a turning point, in the vicinity of which the resolution of the mesh must be sufficiently high. Lee & Saio (1997) showed that a regular solution could be locally expanded in series of the form , where x = μ−1/|ν|, and satisfying the condition (A.5)This allows us to compute the solution properly in the sub-inertial regime. To construct the set of Hough functions, we use as predictors the asymptotic results given by Longuet-Higgins (1968) and Townsend (2003). Two types of modes can be identified.

The modes of the first type, called gravity modes, are defined for ν ∈R and are characterized by positive eigenvalues (Fig. A.1, top panel). In the non-rotating case, where ν = 0, they correspond to the associated Legendre polynomials (Fig. A.1, bottom panel). The gravity modes of Hough functions are ordered with positive n, by analogy with this particular case, for which the degrees l of Legendre polynomials are given by l = |m| + n. In the asymptotic cases where |ν| → + ∞, analytic solutions using Hermite polynomials can be computed (Townsend 2003). In the sub-inertial regime, the oscillations of gravity modes are trapped in an equatorial belt in the interval given by 0 ≤ μ ≤ 1/|ν| (see Fig. A.2). This explains why these modes are not well represented by associated Legendre polynomials for high values of |ν|.

The modes of the second type have various denominations in literature (see Longuet-Higgins 1968). We call them Rossby waves in reference to the early studies by Rossby et al. (1939). These modes are defined for |ν| > 1 only. Most of them are characterized by negative and are identified in this work by subscripts n< 0, e.g. Λ-1 ≥ Λ-2 ≥ Λ-3..., the parity of n corresponding to the parity of the associated function. Longuet-Higgins (1968) established the analytic expressions of the functions in the vicinity of the boundaries ν = ± 1 using Laguerre polynomials. In contrast with gravity modes, the oscillations of Rossby modes are concentrated around the poles in a region defined by 0 ≤ μ ≤ 1/|ν|. As pointed out by Lindzen (1966), Rossby modes have to be taken into account in the regime of sub-inertial waves, otherwise the set of eingenfunctions is not complete.

In Fig. (A.1), we can observe that the eigenvalues () strongly depend on the spin parameter. This dependence will affect the propagation of tidal waves along the vertical direction through their vertical wavenumbers. In addition, the Coriolis acceleration induces a dissymmetry between prograde (ν< 0) and retrograde (ν> 0) modes. Particularly, the gravity mode of lowest degree present two different asymptotic behaviours in the sub-inertial regime. The projection of the associated Legendre polynomial , which is representative of the semidiurnal tidal forcing, is also directly related to ν. We retrieve for the dissymmetry observed for the eigenvalues and note that the variations of the decomposition with ν are not simple to describe. In the present study, this decomposition is computed numerically with the spectral method presented below.

thumbnail Fig. A.2

Normalized Hough functions computed as functions of the colatitude (θ) for 0 ≤ n ≤ 5 (gravity modes) and various values of ν = (2Ω) /σ. Top-left: ν = 0; top-right: ν = 0.8; middle-left: ν = 1.2; middle-right: ν = 2; bottom-left: ν = 4; bottom-right: ν = 12. The case ν = 0 corresponds to the normalized associated Legendre polynomials . The two plots at the top represent the regime of super-inertial waves (|ν| ≤ 1) where Hough functions look like Legendre polynomials. In the regime of sub-inertial waves (|ν| > 1), functions become more and more evanescent when |ν| increases, and oscillations are trapped around the equator in the interval [θνθν], where .

Open with DEXTER

thumbnail Fig. A.3

Normalized Hough functions computed as functions of the colatitude (θ) for − 6 ≤ n ≤ −1 (Rossby modes) and various values of ν = (2Ω) /σ. Top-left: ν = 1.2; top-right: ν = 2; bottom-left: ν = 4; bottom-right: ν = 12. In the regime of sub-inertial waves (|ν| > 1), Rossby modes are trapped around the poles when |ν| → 1.

Open with DEXTER

Appendix A.2: Spectral method

The early studies by Longuet-Higgins (1968) and CL70 provide a detailed theoretical approach of Hough functions and present the very powerful general method introduced by Hough (1898) to compute the eigenvalues and eigenvectors of m,ν. In this method, which is not detailed here, Hough functions are expanded in series of normalized associated Legendre polynomials (Abramowitz & Stegun 1972). In the previous subsection, Eq. (A.1) highlights the symmetry between the two parameters defining these functions, m and ν. Hence, for m> 0 and ν ∈R, (A.6)Therefore, Hough functions and eigenvalues for m< 0 are readily deduced from those obtained for m> 0. This allows us to write (A.7)where s = |m|. A first possible way to compute the coefficients and eigenvalues , described by CL70, consists in finding the roots of a high-degree polynomial (typically, the degree N of this polynomial corresponds to the number of coefficients of the series given by Eq. (A.7)). With the current computing resources, N = 500 can be reached easily. However, this method rapidly becomes tricky to apply beyond N ~ 500 because of very large polynomial coefficients. Since the size of the range depends on N, this can limit the number of gravity modes or Rossby modes found. For instance, in the case of the Earth’s semidiurnal tide, where ν is very close to 1, the first Rossby modes are characterized by strongly negative eigenvalues. Therefore, only a few functions of this type can be computed. To address this issue, the previous problem may be written as an explicit eigenvalue problem. Hence, following Hough (1898) and CL70, we expand the non-normalized Hough functions in non-normalized associated Legendre Polynomials, denoted : (A.8)To lighten the expressions, the are simply denoted in the following analytic development. Let us introduce X = 1/ and the coefficients (A.9)Hough (1898) demonstrated the following equation, which corresponds to Eq. (60) in CL70 (A.10)with (A.11)The system of equations then reduces to the two eigenvalue problems (A.12)and (A.13)with the matrices(A.14)(A.15)They are solved numerically by linear algebra algorithms with great efficiency. For example, the case N = 1000 is computed on a laptop instantaneously. To make the comparison with Lee & Saio (1997) easier (see the previous subsection), we set at the equatorial point (θ = π/ 2) for even solutions and for odd solutions. Hence, the coefficients of Eq. (A.7) are explicitly given by (A.16)where (A.17)The solutions and projection matrices given by the spectral method are plotted in Figs. 2, 3, A.2, A.3. In the super-inertial regime (|ν| ≤ 1), the Hough functions basis is composed of gravity modes alone. Therefore, the projection matrix is a band-matrix. In this regime Hough functions are very similar to the associated Legendre polynomials. In the case ν = 0 (top left), the bases are exactly the same and the projection matrix is diagonal.

The sub-inertial regime is characterized by modes of both types, the transition being identified by Hough functions of the lowest degrees: Rossby modes on the left of the map (lowest eigenvalues) and gravity modes on the right (highest eigenvalues). The maps reveal the side effect inherent to the spectral method. Hough functions of the highest degrees are obviously not well represented by the family of associated Legendre Polynomials; more polynomials are necessary. This is explained by the asymptotic behaviour of Hough functions. At |ν| → 1+, Rossby modes are trapped around the poles and are thus well represented by highly oscillating associated Legendre polynomials (Fig. A.3). At |ν| → + ∞, gravity modes behave similarly and are trapped around the equatorial belt (Figs. 2, A.2).

Since the spectral method has been proved to be very efficient at computing a large number of eigenvalues and eigenfunctions at once, it is used for all calculations in this work. We note that some particular cases cannot be treated using Eq. (A.10). These cases have been discussed by CL70 and will not be developed in this work.

Appendix B: Numerical scheme used to integrate spatial differential equations

We detail here the numerical scheme used to integrate Laplace’s tidal equation (with the finite differences method) and the vertical structure equation. This method, introduced by Bruce et al. (1953) and described by Chapman & Lindzen (1970) for a particular case, is very convenient when solving second-order differential equations with non-constant coefficients. Those equations are written (B.1)where P, Q, and R are given functions of x. The domain on which Eq. (B.1) is integrated is divided into N intervals, i.e. N + 1 points, identified by the index n. Hence, for a regular mesh with elements of size δx, the first- and second-order derivative of a solution f are expressed (B.2)with n ∈ ˜1,N−1¨. Substituting Eq. (B.2) in Eq. (B.1), we obtain the numerical relations (B.3)with the coefficients (B.4)the Pn, Qn, and Rn coefficients being the components of the numerical vectors corresponding to P, Q, and R. Then, we introduce the coefficients αn and βn, such that (B.5)These coefficients are defined straightforwardly by the recursive relations (B.6)The first terms (α0 and β0) are deduced from the boundary condition at x = xinf, which can be written (B.7)where a0, b0, and c0 are real or complex coefficients. Thus, we obtain (B.8)and, step by step, all the following terms. At x = xsup, we apply the condition (B.9)with , which gives (B.10)The solution is finally obtained by computing the fn backwards using Eq. (B.5).

Appendix C: Decomposition of the thermal forcing in Keplerian elements

We establish here the general multipole expansion in spherical harmonics of the thermal power given in Eq. (72) using the notations of Mathis & Le Poncin-Lafitte (2009). We consider a planet (A) orbiting circularly around its star (B). Three reference frames, all centred on the centre of mass of the planet and represented in Fig. C.1, must be defined:

  • RR:{OA,XR,YR,ZR}, time independent inertial frame, with ZR the direction of the total angular momentum of the system.

  • RO:{OA,XO,YO,ZO}, the orbital frame linked to RR by the Euler angles:

    • IB, the inclination of the orbital frame with respect to ;

    • ωB, the argument of the pericentre;

    • , the longitude of the ascending node.

  • RE;T:{OA,XE,YE,ZE}, the spin equatorial frame rotating with the angular velocity ΩA, and linked to RR by three Euler angles:

    • εA, the obliquity, i.e. the inclination of the equatorial plane with respect to ;

    • ΘA, the mean sidereal angle such that ΩA = dΘA/ dt, which is the angle between the minimal axis of inertia and the straight line due to the intersection of the planes and ;

    • φA, the general precession angle.

thumbnail Fig. C.1

Inertial reference, orbital and equatorial rotating frames, and the associated Euler’s angles of orientation. Figure taken from Mathis & Le Poncin-Lafitte (2009).

Open with DEXTER

Finally, we denote a the semi-major axis and the mean anomaly with , nB being the mean motion. Neglecting the action of the centrifugal acceleration, the hydrostatic structure of the atmosphere presents a spherical symmetry, p0, ρ0, T0 being functions only of the radial coordinate. Hence, any heat forcing caused by the star can be written (C.1)where h is a function of r and ψ, the angle between the points M(r,θ,ϕ) and B(a,θAB,ϕAB). We decompose J in normalized Legendre polynomials (C.2)with . Using the addition theorem (C.3)with (C.4)we obtain (C.5)Kaula’s transform is finally introduced to express as a function of Keplerian elements (Kaula 1962; Mathis & Le Poncin-Lafitte 2009). Since the motion of the planet is circular (e = 0), we get (C.6)with (j,p,q) ∈ ˜−l,l¨ × ˜0,l¨ ×Z, where the κl,j coefficients are given by (C.7)The functions denoted , Fl,j,p, and Gl,p,q are called the obliquity, inclination, and eccentricity functions (for more details, see Mathis & Le Poncin-Lafitte 2009). The angle Ψl,m,j,p,q can be written (C.8)the parameter σl,m,p,q = mΩA−(l−2p + q)nB being the tidal frequency of the mode (l,m,p,q) and the phase lag. Substituting Eq. (C.6) in Eq. (C.5), we finally get(C.9)where (C.10)with ϕl,m,j,p,q = ψl,m,j,p,q/m.

With the same notations, the gravitational potential is written (C.11)The tidal frequency, written σ = mΩ−snorb, is defined by the doublet (m,s) ∈ Z2.

The perturbing tidal potential and thermal power finally write (C.12)where the spatial functions Um,σ and Jm,σ are expressed (C.13)with (C.14)the notation δ standing for the Kronecker symbol.

Appendix D: Radiated power as a function of physical parameters

The power per unit mass emitted by an optically thin atmosphere (Jrad) can be expressed as an analytic function of physical parameters in our simplified framework. Indeed, assuming that the flux radiated by a layer of thickness dx, denoted δI, propagates along the vertical direction without being absorbed by the other layers, we write (D.1)where φrad is the total power emitted per unit mass. Then, the atmosphere is treated as a grey body of molar emissivity per unit mass ϵa, which – using the molar concentration C0 introduced in Sect. 2 – gives (D.2)Combining Eqs. (D.1) and (D.2) allows us to get (D.3)Finally, we derive the perturbed radiative loss Jrad associated with the tidal perturbation as a function of δT by linearizing the previous expression around the equilibrium state, (D.4)which gives the radiation frequency introduced in Eq. (17): (D.5)

Appendix E: Analytic solution for constant profiles and upper boundary condition

We show here that the solution of the equation describing the vertical structure strongly depends on the choice of the upper boundary condition. As discussed in Sect. 3, (E.1)at the ground level x = 0.

Instead of the energy-bounded condition prescribed in the model, at the upper limit (xatm) we set the commonly used stress-free condition, (E.2)We obtain, at x = xatm, (E.3)with the complex coefficients (E.4)The parameter xatm must be fixed so that , which corresponds to xatm ≫ 1. A simple analytic solution can be computed assuming that and are constant profiles. The second member of the vertical equation, C, given by Eq. (99), becomes (E.5)the solution is written (E.6)with and , and the complex integration constant (E.7)Now, the condition of Eq. (108) characterizing thermal tides, i.e. (E.8)is assumed. The second-order Love number reduced to the main contribution of the thermal semidiurnal tide, given by Eq. (76), can be expressed (E.9)

thumbnail Fig. E.1

Tidal torques exerted on the atmosphere (in peta-Newton metres (PN.m)) induced by the first gravity mode of the Earth’s semidiurnal tide (n = 0) as a function of the reduced tidal frequency (χ) obtained with two different conditions at the upper boundary: the energy-bounded condition chosen for the model (red line) and a stress-free condition (green dashed line). The equation for the structure in the vertical (resp. horizontal) direction has been integrated with 10 000 mesh points (resp. 400 associated Legendre polynomials).

Open with DEXTER

for a thin isothermal atmosphere. As a zero-order approximation, we consider that the forcings are the constant profiles U2 = U2,2,2,0,0 and J2 = J2,2,2,0,0 used in Sect. 2 (with J2 ∈R). Therefore, introducing the complex impedance Ξatm, defined by (E.10)and substituting the analytic solution (Eq. E.6) in Eq. (117), we get (E.11)where Ξatm is explicitly given by (E.12)with (E.13)Similarly, the torque of Eq. (88) associated with the -contribution of the semidiurnal tide can be written (E.14)To illustrate the difference induced by the choice of the upper boundary condition, we compute the variation of the tidal torque generated by the Earth’s semidiurnal tide with the reduced tidal frequency (χ = (Ω−norb) /norb) for two conditions: the energy-bounded condition used in this work and the free-surface condition studied here. To simplify calculations, we focus on the gravity mode of lowest degree (n = 0). Results are plotted in Fig. E.1. In the first case, we obtain the smooth evolution observed in Fig. 14 (with slight differences in the vicinity of the synchronization because of the absence of Rossby modes), but the free-surface condition obviously leads to a highly resonant behaviour. This can be explained by the form of the solution given by Eq. (E.6), where the integration constant depends on and xatm. By fixing xatm, we implicitly define values of for which reaches a maximum. On the contrary, the integration constant obtained by applying the energy-bounded condition, given by Eq. (128), does not depend on the upper limit far beyond the critical altitude given by Eq. (127).

Appendix F: Supplementary materials regarding the dependence of the tidal response on the tidal frequency

thumbnail Fig. F.1

Tidal response of the Earth’s atmosphere to a quadrupolar forcing ( and ) as a function of the reduced tidal frequency χ = (Ω−norb)/norb (horizontal axis). For each quantity qn associated with the Hough mode of degree n, except the coefficients (top left) and the scales of variations LV;n, the function f(χ) = sign(qn)log (1 + |qn|) is plotted (vertical axis) for n ∈ ˜−6,5¨. Left, from top to bottom: a) coefficients (Eq. (70)); b) eigenvalues of Laplace’s tidal equation ( defined by Eq. (35)); c) normalized variation scales of the vertical structure of the vertical component, LV;n = 2π/ (Eq. (48)), plotted in logarithmic scale. Right, from top to bottom: d) real part of the vertical wavelengths (, Eq. (103)); e) imaginary part of the vertical wavelengths; f) spin parameter ν = 2Ω /σ (Eq. (23)).

Open with DEXTER

As a complement to the tidal torques (Figs. 14) and perturbed pressure, density, and temperature (Figs. 17) resulting from the Earth’s solar semidiurnal tide, we plot several quantities as a function of the tidal frequency that can be helpful for the diagnostic of results (Fig. F.1). These quantities are the coefficients introduced in Eq. (81) which weight the contribution of Hough modes to the tidal torque; the eigenvalues of Laplace’s tidal equation (); the length scale of vertical variations associated to Hough modes (LV;n); the spin parameter (ν); and the real and imaginary parts of the complex vertical wavenumber ( denoted kn in plots), which determines the nature (propagative or evanescent) of tidal waves.

These plots highlight the difference between prograde (ν< 0) and retrograde (ν> 0) modes. Indeed, it can be noted that the projection of the associated Legendre polynomial on Hough functions is not the same before and after the synchronization: the weight of the gravity mode of degree 0, which is the most representative one, varies and Rossby modes are predominating in a very short interval of positive tidal frequencies. This dissymmetry is retrieved in the variations of the horizontal eigenvalues, where Rossby modes only appear for σ> 0.

As predicted by its analytic expression (Eq. (102)), the real and imaginary part of the vertical wavenumber diverge at the synchronization. Therefore, Hough modes become both highly oscillating and strongly damped in this region. Their spatial variation scale along the vertical direction decays, which corresponds to what is observed in the variations of LV;n. Actually, below a given typical length scale, these modes are damped by dissipative processes such as thermal diffusion. Finally, thegraph of ν shows that this parameter strongly varies around the synchronization. The transition at σ = 0 is discontinuous, ν switching from − ∞ to + ∞.

Appendix G: Consequences of an insufficient spatial resolution on the numerical computation of the vertical structure

As pointed out by CL70 and detailed in Sect. 5, solving numerically the vertical structure equation (Eq. (107)) with the numerical scheme of Appendix B requires using a sufficiently high spatial resolution. This condition is mathematically expressed by the criterion given by Eq. (228) for a regular mesh and depends on the vertical wavenumbers that characterize Hough modes. To compute the response of a mode, it is necessary to use a spatial step δx lower than the corresponding variations length scale (LV;n, defined by Eq. (48)), otherwise the signal does not appear. We illustrate this point in here by considering the simplified forced wave equation (G.1)where kV ∈C is such that ℑ{kV} > 0, with the previously chosen boundary conditions, (G.2)We arbitrarily set kV = 100(1 + i) in order to satisfy the condition LV ≪ 1 (keeping in mind that xsup = 10). This corresponds to a strongly damped tidal mode similar to those predicted by the model in a frequency range close to synchronization. The critical number of points necessary to compute the response with our scheme is then Ncrit ~ 400. We solve Eq. (G.1) for various resolutions. The corresponding results are plotted in Fig. G.1 where they are also compared to the analytical solution given by Eq. (129). The flattening observed for N<Ncrit gives an idea of the kind of numerical difficulties that need to be addressed to compute tidal low-frequency waves with a GCM. Particularly, this shows that if we compute the tidal torque numerically by using the scheme of Appendix B with an insufficient resolution, we will only be able to get the large-scale contribution associated with the horizontal component (see Figs. 14).

thumbnail Fig. G.1

Real part of the solution of Eq. (G.1) as a function of x in logarithmic scale computed numerically for various mesh resolutions. The number of points of the spatial discretization is set at N = 20 (green dashed line), N = 200 (blue dashed line), N = 2000 (pink dashed line), and N = 20 000 (cyan dashed line). These numerical solutions are compared to the analytical solution given by Eq. (129) (red line).

Open with DEXTER

Appendix H: Nomenclature

Table H.1

List of the symbols used throughout this work and their designations in order of introduction.

All Tables

Table 1

Values of physical parameters used in simulations.

Table H.1

List of the symbols used throughout this work and their designations in order of introduction.

All Figures

thumbnail Fig. 1

Spherical coordinate system associated with the equatorial reference frame of the planet, and the geometrical parameters of the system.

Open with DEXTER
In the text
thumbnail Fig. 2

Normalized Hough functions , , and corresponding to the Earth’s semidiurnal tide (m = 2 and ν = 1.0030). We note that ν = 2Ω /σ. Left: first gravity modes; 0 ≤ n ≤ 5. Right: first Rossby modes: − 6 ≤ n ≤ −1. In this case, |ν| is very close to 1 and the Rossby modes are trapped at the poles. The Earth’s semidiurnal tide is mainly described by gravity modes. For more details about the behaviour of Hough functions, see Fig. A.2.

Open with DEXTER
In the text
thumbnail Fig. 3

First gravity modes (0 ≤ n ≤ 5) of normalized Hough functions , , and corresponding to the semidiurnal tide of Venus-like exoplanets (m = 2 and ν = 0.4804). We note that ν = 2Ω /σ. For more details about the behaviour of Hough functions, see Figs. A.2.

Open with DEXTER
In the text
thumbnail Fig. 4

Frequency spectrum of the atmospheric tidal response. The identified regimes are given as functions of the tidal frequency (σ).

Open with DEXTER
In the text
thumbnail Fig. 5

Real and imaginary part of the vertical wavenumber as functions of the tidal (σ) and critical (σc;n) frequencies. Left: , with defined in Eq. (104), as a function of f(σ/σ0) (horizontal axis) and (vertical axis), where f is the function defined by . The colour of the map, denoted c, is defined by . Hence, the purple region stands for slowly oscillating waves. The luminous (dark) regions correspond to strongly oscillating waves propagating downwards (upwards).Right: as a function of the same parameters. The luminous (dark) regions correspond to strongly (weakly) evanescent waves. The positions of the planets (Earth , Venus ) on the map are determined by the gravity mode of degree 0 in the case of the semidiurnal tide (Λ;0 = 11.1596, Λ;0 = 7.1485). The tidal frequencies of the planets are given by σ/σ0 ≈ 194 and σ/σ0 = 0.43, with σ0 = 7.5 × 10-7 s-1.

Open with DEXTER
In the text
thumbnail Fig. 6

Relative difference between the complex wavenumber obtained in this work (Eq. (103)) on the one hand, and the real one given by CL70 for Earth atmospheric tides () on the other hand. The difference η, given by Eq. (158), is plotted in logarithmic scale as a function of the normalized tidal frequency σ/σ0 and of the ratio hc/hn. The colour of the map, denoted c, is defined by the expression c = log (η). The horizontal and vertical axis represent f(σ/σ0) and , where f is the function defined by with X ∈R. The positions of the planets are determined by the gravity mode of degree 0 in the case of the semidiurnal tide. Dark (luminous) regions correspond to the dynamic (dissipative) regime. The vertical wavenumbers of Venus are strongly affected by dissipative processes, while those of the Earth are not.

Open with DEXTER
In the text
thumbnail Fig. 7

Tidal torques exerted on the atmosphere in the case of the slowly rotating convective atmosphere (red line with the label “Conv”, Eq. (174)) and in the model by Correia & Laskar (2001) (green dashed line with the label “CL01”, Eq. (175)) as functions of the tidal frequency. Both functions are normalized by their respective maximum value, reached at σ = σ0; the tidal frequency is normalized by the Newtonian cooling frequency (σ0). The tidal forcing is quadrupolar, with U and J given by Sect. 4.3. The parameters a and b of Eq. (175) are adjusted so that the peaks of the two functions are superposed.

Open with DEXTER
In the text
thumbnail Fig. 8

Left: absorption of the flux emitted by the star along the star-planet direction. Right: absorption of the flux emitted by the ground.

Open with DEXTER
In the text
thumbnail Fig. 9

Vertical profiles of radiative fluxes (I), power per unit volume (J), and power per unit mass (J) for various values of the dimensionless optical depth τλ (Eq. (187)). As usual, the reduced altitude, x = z/H, is on the vertical axis. Left: incident flux, defined by Eq. (186), at ψ = 0. If τλ ≪ 1 the ground receives almost 100% of the total power. On the contrary, if τλ ≫ 1, all the power is absorbed by the atmosphere. Right: radiative flux emitted by the ground, defined by Eq. (206). For small τλ, the flux goes through the atmosphere without being absorbed. For high values of the parameter, it is absorbed by the lowest layers.

Open with DEXTER
In the text
thumbnail Fig. 10

Left, from top to bottom: power per unit mass (Jinc, Eq. (184)) and power per unit volume (Jinc, Eq. 183) of a monochromatic incident flux normalized by their maxima and plotted as functions of the normalized spatial coordinates X (the direction of the star is along the horizontal axis) and Z (the spin axis is along the vertical axis). The planet is represented by a green disk of unit radius. The flux, of wavelength λ, propagates in an atmosphere assumed to be homogeneous in composition (CG = C0) and characterized by the ratio H/R = 0.1. The optical depth is given by τλ = 0.2. Right, from top to bottom: power per mass unit (Jrad, Eq. (206)) and power per volume unit (Jrad) of a monochromatic flux radiated by the ground, normalized by their maxima, and plotted as functions of X and Z, with τλ = 0.5.

Open with DEXTER
In the text
thumbnail Fig. 11

Nyquist diagram of the boundary layer transfer function in the frequency range − 103σ/σbl ≤ 103 (σ being the tidal frequency and σbl the frequency characterizing the boundary layer introduced in Eq. (217)). The parametrized complex transfer function (Eq. (215)) reduced by its value at σ = 0 is plotted as a function of the tidal frequency (σ) in the complex plane, where the horizontal axis corresponds to the real part of the function, , and the vertical axis to the imaginary part . The arrow indicates the direction of increasing σ. The gain of the transfert function is given by the distance of a point P of the curve to the origin O of coordinates , and its phase by the angle of the vector OP to the horizontal axis. This plot shows that the response of the ground is always delayed with respect to the thermal forcing except at synchronization, where and the gain is maximum. It also highlights the two identified regimes of the ground response: the lag tends to 0 (response in phase with the perturbation) and the gain increases at σ → 0, while the lag tends to ± π/ 4 and the gain to 0 at σ → ± ∞, the critical transition frequency being σbl.

Open with DEXTER
In the text
thumbnail Fig. 12

Bichromatic case. Top: heat transfers between the atmosphere, space and the ground. Bottom: power balance of the interface ground-atmosphere.

Open with DEXTER
In the text
thumbnail Fig. 13

Average daily variation in temperature in degrees (top) and pressure in mbar (bottom) over one year, plotted over solar time (in hours), the subsolar point corresponding to 12 h. The data was recorded every minute over the full year 2013 with a Vantage Pro 2 weather station at latitude 48.363°N. The time is converted to solar time and the data averaged over an exact count of 365 days. The average values of temperature (10.8973° C) and pressure (1015.83 mbar) have been removed to display only the variable part. In both figures, the raw data (in green) has been decomposed in a Fourier series over the 24 h period, limited to two harmonics (blue), or five harmonics (red). Temperature is dominated by the diurnal component, with a maximum value around 14h26mn. The largest harmonics of the pressure variations is the semidiurnal term, followed by the diurnal term. Two pressure maxima can be observed, at 9h38mn and 22h07mn.

Open with DEXTER
In the text
thumbnail Fig. 14

Tidal torque applied to the atmosphere of Earth (left panel) and Venus (right panel) by the thermal semidiurnal tide and its components as functions of the reduced frequency χ = (Ω−norb) /norb (horizontal axis). The forcings consist of academic quadrupolar perturbations expressed as and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3 for the Earth and U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3 for Venus. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively (we note that ). The label CL01 refers to the model by Correia & Laskar (2001), given by Eq. (175).

Open with DEXTER
In the text
thumbnail Fig. 15

Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3. Left, from top to bottom: amplitudes of δp, δρ, δT, Vθ, Vϕ, and Vr in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (|δf|) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)].

Open with DEXTER
In the text
thumbnail Fig. 16

Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3. Left, from top to bottom: real parts of δp, δρ, δT, Vθ, and Vϕ taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}.

Open with DEXTER
In the text
thumbnail Fig. 17

Tidal surface oscillations of pressure, density, temperature, and their components at the equator of the Earth as functions of the reduced tidal frequency χ = (Ω−norb) /norb (horizontal axis). The perturbation is defined by the academic quadrupolar forcings and with U2 = −0.985 m2 s-2 and J2 = 1.0 × 10-2 m2 s-3. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively.

Open with DEXTER
In the text
thumbnail Fig. 18

Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3. Left, from top to bottom: amplitudes of δp, δρ, δT, Vθ, Vϕ, and Vr in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (|δf|) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)].

Open with DEXTER
In the text
thumbnail Fig. 19

Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3. Left, from top to bottom: real parts of δp, δρ, δT, Vθ, and Vϕ taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}.

Open with DEXTER
In the text
thumbnail Fig. 20

Tidal surface oscillations of pressure, density, temperature, and their components at the equator of Venus as functions of the reduced tidal frequency χ = (Ω−norb) /norb (horizontal axis). The perturbation is defined by the academic quadrupolar forcing and with U2 = −2.349 m2 s-2 and J2 = 1.0 × 10-4 m2 s-3. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively.

Open with DEXTER
In the text
thumbnail Fig. 21

Asymptotic domains of parameters where the analytical approach presented in this work can be applied (green) and where the numerical approach is necessary (red) as functions of the frequency ratios N2/σ2 (horizontal axis) and (2Ω)2/σ2 (vertical axis). Green regions: stably stratified atmospheres satisfying the hierarchy of frequencies (2Ω)2N2 necessary to apply the traditional approximation a); weakly stratified or convective atmospheres of slowly rotating planets b). Red regions: stably stratified atmospheres where the hierarchy of frequencies (2Ω)2N2 is not satisfied c); super-adiabatic rotating convective atmospheres d).

Open with DEXTER
In the text
thumbnail Fig. 22

Predictive aspects related to the dependence of the atmospheric tidal torque on stratification. The neutral convective case implies a strong torque, which can lead to non-synchronized states of equilibrium. The stably stratified case implies a weak torque, which lets the planet evolve towards spin-orbit synchronization.

Open with DEXTER
In the text
thumbnail Fig. A.1

Top: eigenvalues () of the gravity and Rossby modes of lowest degrees (− 6 ≤ n ≤ 5) as functions of the spin parameter (ν). The horizontal axis represents the function of ν given by , and the vertical axis the function of given by . Bottom: projection coefficients of the normalized Legendre polynomial on normalized Hough functions associated with these modes as functions of the spin parameter (see Eq. (70)).

Open with DEXTER
In the text
thumbnail Fig. A.2

Normalized Hough functions computed as functions of the colatitude (θ) for 0 ≤ n ≤ 5 (gravity modes) and various values of ν = (2Ω) /σ. Top-left: ν = 0; top-right: ν = 0.8; middle-left: ν = 1.2; middle-right: ν = 2; bottom-left: ν = 4; bottom-right: ν = 12. The case ν = 0 corresponds to the normalized associated Legendre polynomials . The two plots at the top represent the regime of super-inertial waves (|ν| ≤ 1) where Hough functions look like Legendre polynomials. In the regime of sub-inertial waves (|ν| > 1), functions become more and more evanescent when |ν| increases, and oscillations are trapped around the equator in the interval [θνθν], where .

Open with DEXTER
In the text
thumbnail Fig. A.3

Normalized Hough functions computed as functions of the colatitude (θ) for − 6 ≤ n ≤ −1 (Rossby modes) and various values of ν = (2Ω) /σ. Top-left: ν = 1.2; top-right: ν = 2; bottom-left: ν = 4; bottom-right: ν = 12. In the regime of sub-inertial waves (|ν| > 1), Rossby modes are trapped around the poles when |ν| → 1.

Open with DEXTER
In the text
thumbnail Fig. C.1

Inertial reference, orbital and equatorial rotating frames, and the associated Euler’s angles of orientation. Figure taken from Mathis & Le Poncin-Lafitte (2009).

Open with DEXTER
In the text
thumbnail Fig. E.1

Tidal torques exerted on the atmosphere (in peta-Newton metres (PN.m)) induced by the first gravity mode of the Earth’s semidiurnal tide (n = 0) as a function of the reduced tidal frequency (χ) obtained with two different conditions at the upper boundary: the energy-bounded condition chosen for the model (red line) and a stress-free condition (green dashed line). The equation for the structure in the vertical (resp. horizontal) direction has been integrated with 10 000 mesh points (resp. 400 associated Legendre polynomials).

Open with DEXTER
In the text
thumbnail Fig. F.1

Tidal response of the Earth’s atmosphere to a quadrupolar forcing ( and ) as a function of the reduced tidal frequency χ = (Ω−norb)/norb (horizontal axis). For each quantity qn associated with the Hough mode of degree n, except the coefficients (top left) and the scales of variations LV;n, the function f(χ) = sign(qn)log (1 + |qn|) is plotted (vertical axis) for n ∈ ˜−6,5¨. Left, from top to bottom: a) coefficients (Eq. (70)); b) eigenvalues of Laplace’s tidal equation ( defined by Eq. (35)); c) normalized variation scales of the vertical structure of the vertical component, LV;n = 2π/ (Eq. (48)), plotted in logarithmic scale. Right, from top to bottom: d) real part of the vertical wavelengths (, Eq. (103)); e) imaginary part of the vertical wavelengths; f) spin parameter ν = 2Ω /σ (Eq. (23)).

Open with DEXTER
In the text
thumbnail Fig. G.1

Real part of the solution of Eq. (G.1) as a function of x in logarithmic scale computed numerically for various mesh resolutions. The number of points of the spatial discretization is set at N = 20 (green dashed line), N = 200 (blue dashed line), N = 2000 (pink dashed line), and N = 20 000 (cyan dashed line). These numerical solutions are compared to the analytical solution given by Eq. (129) (red line).

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.