© F. Debras et al. 2019

1. Introduction

The study of the influence and propagation of waves in planetary atmospheres is often performed using models with several simplifications, most notably the assumption of a zero or zonally–symmetric and constant initial zonal flow (e.g. Kasahara & Qian 2000), or restriction to a beta–plane solution (e.g. Lindzen 1967). Such simplifications allow analytical prediction of the key wave mechanisms, and in some cases a complete understanding of their structure (even in the mathematical sense; see Matsuno 1966).

Despite the simplifications, such studies have provided significant insight into atmospheric dynamics. Wheeler & Kiladis (1999) for example demonstrate that the propagation of waves can be linked to convective motions in the atmosphere of the Earth, and the resulting description matches analytical theories (e.g. Vallis 2006; Holton 1992) remarkably well. Baroclinic and barotropic instabilities are also known to have an impact on the circulations of planetary atmospheres (see Williams 2003).

However, for more complex mean flows or situations where individual terms in the hydrodynamical equations are not clearly dominant, analytical treatments rapidly become impractical. Additionally, some wave structures, or modes, are only supported by the full equations, being effectively “filtered out” by the simplifications. Notably, Wang & Mitchell (2014) numerically identify a Rossby–Kelvin wave mode for a planetary atmosphere that cannot be recovered in the quasi–geostrophic equations (see Gill 1980 or Vallis 2006 for further details).

The detection and characterisation of a specific class of exoplanets, hot Jupiters, has provided impetus to the study of non-axisymmetrically forced planetary atmospheres. Hot Jupiters are Jovian planets in short-period orbits close to their host star, and likely have synchronised rotational and orbital speeds such that a single hemisphere or day side faces the host star at all times (see Baraffe et al. 2010). The slow rotation (periods of ∼4 days) and Jovian radii suggest such atmospheres exist in a regime where the Rossby number is of order unity , where U is the characteristic flow velocity, L a characteristic length scale, and f the Coriolis parameter, meaning rotation is neither dominant nor negligible. Observational evidence has indicated the presence of fast zonal “jets” (zonally coherent flows) of a few km s⁻¹ (Louden & Wheatley 2015). The mechanism for accelerating these zonal flows has been explored by Showman & Polvani (2011), building on the linear studies of Matsuno (1966) and Gill (1980) and using a two-layer equivalent depth approach.

To further study the linear waves and instabilities present in these atmospheres, we developed a public code, ECLIPS3D¹ (Eigenvectors, Circulation and Linear Instabilities for Planetary Science in 3 Dimensions), which we benchmark in this work. More globally, this code can be used in the study of linear stable or unstable modes within any planetary atmosphere from an arbitrary initial steady state.

We expand upon Thuburn et al. (2002), who studied propagating wave modes in an axisymmetric atmosphere at rest to include linear modes in an atmosphere with a steady background flow in three-dimensional spherical coordinates. We detail the structure of ECLIPS3D including the equations solved and the process of obtaining a solution. We detail the different sets of equations implemented (axisymmetric 2D, 3D, two-layer equivalent depth) as well as the time-dependent or independent solutions (waves, instabilities, and standing circulation). Finally, we implemented a posteriori calculation of an energy equation for a given solution. These semi-analytical results allow the verification of the frequency of the modes (and the growth and damping rate for instabilities), as well as isolation of the dominant mechanism providing insight into the physical phenomena driving the instability.

In Sect. 2, we outline the equations implemented in ECLIPS3D (with the full equations detailed in Appendix A), and the procedure for solving them, alongside the method of calculating the energy equation. We then benchmark ECLIPS3D against a range of classical calculations of waves, instabilities, and circulations in Sect. 3, including a setup similar to Showman & Polvani (2011). Finally, we draw conclusions and comment on future developments and applications for ECLIPS3D in Sect. 4.

2. The algorithm

2.1. Linearised equations

Thuburn et al. (2002) show that even the simplest atmospheric waves exhibit behaviour that cannot be accurately expressed by separating variables (requiring a height-dependent shift in latitude; see Thuburn et al. 2002, for details). Therefore, in the general case, no assumption can be made on the mathematical expression of the wave regarding spatial coordinates. As our steady state is arbitrary, we linearise the full equations with no simplification. However, to more easily describe the main capabilities of ECLIPS3D we detail how ECLIPS3D solves the Euler equations, omitting diffusion or viscosity, although dissipative processes have been implemented (and discussed in the steady circulation case in Sect. 3.5 and Appendix A.5). The basic equation set is as follows.

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

(1g)

where u, v, and w are the components of velocity in the longitudinal (λ), latitudinal (ϕ), and vertical (r) directions, ρ is the density, p the pressure, T the temperature, θ the potential temperature, p₀ a reference pressure, R the gas constant divided by mean molecular weight, c_p the heat capacity, g the gravitational acceleration (and is a function of r, see Appendix A.1), Ω the rotation rate of the planet, r the radial distance from the centre of the planet, λ the longitude, ϕ the latitude, and finally Q is the heating rate (if present). Equations (1a)–(1c) represent momentum conservation, (1d) mass conservation, (1e) conservation of energy, (1f) is the equation of state (here an ideal gas), and (1g) defines potential temperature, closing the set. Thuburn et al. (2002) showed that potential temperature is more appropriate than normal temperature for studies of the linear modes. Finally, D/Dt is the Lagrangian or material derivative and t is time.

Solving for waves or instabilities then requires these equations to be linearised. We follow the definitions of Thuburn et al. (2002) for the perturbed variables (this scaling comes from Daley 1988), which greatly simplify the equations when the steady state is axisymmetric and at rest. This choice was made for easier comparison and benchmarking with Thuburn et al. (2002), but a user of the code can change the implemented equations easily without affecting the method of solution. Namely, we write

(2a)

(2b)

(2c)

(2d)

(2e)

where a prime denotes a linearised variable and an i subscript the initial steady state.

If the heating rate Q is non-zero, it has to be properly included in the linearised equations. When solving for waves and instabilities, we simply linearise Q and include it in the left-hand side of the equations. This is detailed in Appendix A.4. When looking for steady, linear circulation (and not free or forced waves) Q(r, ϕ, λ) is specified and considered small enough to only trigger a linear response. A dissipative mechanism must also be added in order to reach a linear steady state. This setup is similar to that of Showman & Polvani (2011), which is one of our benchmark cases, and is detailed further in Appendix A.5.

Thuburn et al. (2002) considered the linearised equations for the case where the initial atmospheric state is axisymmetric, in hydrostatic balance and at rest. In the more general case, linearisation of each of the terms from Eqs. (1a)–(1g) must be completed as shown in Appendix A, alongside the resulting final equation set, Eqs. (A.2)–(A.5). These final, linearised equations are then discretised and solved within ECLIPS3D as detailed in Sect. 2.4. Additionally, we implemented a two-layer model following Showman & Polvani (2011), based on their Eqs. (9) and (10) (linearised versions of which are given in Appendix A.3). Other equation sets (e.g. shallow water, anelastic, and so on) and geometries could be implemented within ECLIPS3D with relative ease if required.

2.2. Boundary conditions

For inviscid flows, there must be no normal flow at the limits of the domain in order to obtain a well posed problem with complete boundary conditions. Namely, we impose that v′cos(ϕ) tends to zero at the poles, and w′ is zero at the top of the atmosphere as a no escape condition. At the inner boundary, we impose a solid boundary with w′ equal to zero. For hot Jupiters, this requires the modelled domain to extend to high-enough pressures for the atmosphere to reach a quiescent region not involved in the acceleration of superrotation. This inner boundary condition can be easily changed if mass flows or energy transfer with the deep atmosphere need to be considered. If the density of the upper atmosphere is too low, unphysical velocities might arise. In the physical applications of ECLIPS3D so far, we have solved this problem by reducing the extent of the atmosphere but a smoothing of the higher atmosphere could be implemented (as done for example in GCMs, see Mayne et al. 2014a,b).

With the above choice of boundary conditions and implemented equations, the code will only recover standing waves in the vertical direction. However, Wu et al. (2001) expressed the importance of vertical wave propagation in the context of Matsuno-Gill structures (Matsuno 1966; Gill 1980), relevant for understanding the climate of Earth (see notably Sarachik & Cane 2010) but also regarding the spin-up of superrotation in hot Jupiters (Showman & Polvani 2011). Additionally, a numerical way to mimic evanescent waves is to impose a damping region at the top of the atmosphere, as is done for example for GCM studies of hot Jupiters (see Mayne et al. 2014a), which prevents the wave from reflecting but allows it to propagate.

In this paper, the boundary conditions have been chosen to enable comparison with previous studies, but any user of the code can easily apply different boundary conditions. Additionally, adapting the equations to include a damping layer in ECLIPS3D poses no theoretical or numerical issue. Vertically propagating waves can therefore be recovered with ECLIPS3D.

2.3. Energy equation

Following the method of Thuburn et al. (2002), we calculated an energy equation by combining the linearised Euler equations and integrating them over the whole atmospheric volume with appropriate boundary conditions. We derive this equation in the same context as Thuburn et al. (2002), with an initially axisymmetric, hydrostatically balanced atmosphere at rest. The general case is shown in Appendix B. We assume that the linearised variables X′ can be expressed as X′(r, ϕ, λ, t) = X(r, ϕ)e^{−i(σt + mλ)}, where the real part of σ is the mode frequency and its imaginary part the growth rate (if non-zero), X(r, ϕ)∈ℂ, as explained in Sect. 2.4, and m ∈ ℤ:

(3)

where the star symbol ^⋆ denotes the complex conjugate, with f = 2Ωcos(ϕ), F = 2Ωsin(ϕ), a the radius of the planet, H the height of the top of the atmosphere, and the sum of the kinetic, thermobaric, and elastic energies of the perturbation, with and the initial buoyancy frequency and sound speed, respectively. Using integration by parts we can express:

(4)

with V the volume. This is possible only if ∫∫∫_VEr²cos(ϕ) drdϕdλ ≠ 0 which is the case where . For this work, we assume a stably stratified atmosphere (). As stated by Thuburn et al. (2002), Eq. (4) shows that σ can only be real in this case and no instability can grow around an atmosphere at rest with no heating and the boundary conditions we have described.

Once the variables u′,v′,… are known Eq. (4) can then be integrated numerically in ECLIPS3D and be used to verify the obtained frequency and identify the dominant physical processes (e.g. in Sect. 3.1 we show that an acoustic wave is largely dominated by the terms involving the pressure, whereas a Rossby wave is dominated by the f and F terms).

When using the code, the use of this a posteriori energy equation is therefore a powerful tool both for diagnosis of the dominant physical mechanism, as well as a validation of the numerical results. It requires an interpolation of the output variables and their derivative on a common grid, which in the current version is performed with linear interpolation, but the consistency of the results with the energy equation confirm that a more sophisticated interpolation would not change the physical interpretation provided by this energy equation (see notably Sect. 3.2).

2.4. Method of solution

The derived linearised equations of motion can be expressed as:

(5)

where D is a differential linear operator. If we introduce the operator A being:

(6)

it is clear that A commutes with D as the latter is only time-dependent through the ∂/∂t terms. Therefore, the vector subspaces of D remain stable upon application of A. As A is diagonal the kernel of D can be decomposed on the eigenvectors e_σ of A. Such eigenvectors are well known: e_σ ∝ e^−iσt where σ ∈ ℂ (the −i term is simply the convention we choose; the real part of σ is therefore the frequency and the imaginary part the growth rate). Therefore, coupled with appropriate boundary conditions, we can then solve Eq. (5) by decomposing them as:

(7)

with and remembering that the evolution of the perturbed quantity is then the real part of the above expression. The actual solution is an infinite sum over all σ but the eigenmodes, and therefore the waves, are the individual projections to a single value. It is worth noting that, a priori, σ could take continuous values. For an atmosphere initially at rest, Matsuno (1966) show that only discrete values are solutions, but on the other hand baroclinic waves exhibit a continuous range of frequencies (see e.g. Charney 1947)

2.4.1. Time-dependent solution

From the discussion above, we can re-write our equations as a complex eigenvalue–eigenvector problem:

(8)

(with B = D − A).

Equation (8) can become difficult or even impossible to solve analytically. However, Thuburn et al. (2002) discretise this equation using a staggered grid of points, turning the analytical matrix B into a finite numerical matrix. This allows spatial derivatives to be calculated using finite differences. The staggered grid was selected carefully for precision and stability by Thuburn et al. (2002), where the 2D axisymmetric version is presented. We adopt a staggered grid in ECLIPS3D, shown in Fig. 1, which resembles the one used in Thuburn et al. (2002) but adapted to 3D and with different staggering of the u and v variables at the poles to simplify the boundary condition.

Fig. 1. Example cell of the 3D staggered grid adopted in ECLIPS3D, based on Thuburn et al. (2002). i discretises the longitudinal variable λ, j the latitude ϕ, and k the radial variable r. Here, u and v are staggered in latitude, with v running from the south to north pole and thereby having an additional latitude point. Also, p is staggered in height with w and θ with the latter two variables running from the bottom to the top of the atmosphere resulting in an additional height point.

We have N_tot = N_λ(2N_ϕN_r + (N_ϕ + 1)N_r + 2N_ϕ(N_r + 1)) points in our grid, with N_λ, N_ϕ, and N_r being the number of points in each coordinate, meaning the matrix B will be of size (5 * N_tot)² as there are five inter-dependent variables. However, each variable at a given point only depends on the values of all variables over the closest points in the grid. Therefore, B is an extremely sparse matrix.

Once the matrix B is filled with discretised values from Eq. (8), at the staggered grid points we must find the eigenvectors of this matrix. ECLIPS3D uses the ScaLAPACK² library for parallel linear algebra (Blackford et al. 1997). To express the eigenvectors of a complex matrix we first calculate the upper Hessenberg form of the matrix, then find the Schur decomposition before identifying the eigenvectors themselves³. Finally, the eigenvectors are returned to their original form via multiplication with the matrix of transformation.

However, this process of solving for the eigenvectors yields an eigenvector for each row in the matrix. From this set we must select those of interest, representative of physical modes in the atmosphere in question. To identify the interesting eigenvectors we employ two methods. Firstly, for instabilities with positive exponential growth rates, we assume that the modes that will lead the dynamical instability have the highest growth rate, and only select the fastest-growing modes. For the case of no instability we first filter out modes arising from numerical errors (e.g. extreme values at the poles), and then manually select modes from the solution set. Analytical expectations then determine the modes of interest. Globally, this selection process needs to be performed thoroughly and be based on analytical expectations of the modes to look for. In the current version of the code, the selection is performed independently of the matrix calculation, and therefore allows different eigenvectors to be isolated in a single ECLIPS3D run.

2.4.2. Time-independent solution

For the case of a steady circulation, without time-dependence and with constant heating rate, Eq. (5) can be expressed as:

(9)

where C is similar to B in Eq. (5) with the inclusion of a drag term if required (see Appendix A.5). Solving this problem is much easier than the time-dependent case, as we simply need to express C on the staggered grid and invert it to obtain the unique solution to these equations.

2.5. ECLIPS3D resolution

In order to achieve the highest possible resolution, we implemented two versions of ECLIPS3D: one that solves for the whole eigenvector spectrum and another one dedicated to solving a reduced number of selected eigenvectors.

For the case of numerically solving for all potential eigenvectors, the computational expense (both in computation time and memory) increases steeply with the number of points in the matrix, and can rapidly become inhibitive, limiting the resolution. Specifically, within our computational framework, ECLIPS3D can be used to solve 2D problems (axisymmetric, two layer, or barotropic equations; see following section) in cases with up to 100 × 100 points within a day of real time. The efficient parallelisation of eigenvector calculations is still an active area of research in the computer science community, and increasing the number of processors does not significantly accelerate calculations of this type. Higher-resolution problems therefore take much longer, with our benchmarking tests suggesting the time taken scales with the number of points as roughly or . Additionally, eventually with increasing numbers of processors, the communication between the processors becomes the primary overhead or limitation in the calculation, whereas using too few processors leads to saturation of the available memory (see the ScaLAPACK documentation for details). These computational limitations are amplified in 3D, where calculations at resolutions of 25 × 20 × 20 points on 64 processors require around four days. We are currently working on adapting ECLIPS3D to employ sparse matrix libraries (as discussed, the matrix we are solving is sparse) and are following the developments in computer science research regarding eigenvector calculations.

However, by solving for a reduce set of selected eigenvectors we can commensurately increase the resolution of the ECLIPS3D setup, whilst retaining a similar computational expense to the case where the target eigenvectors are not restricted. This approach can be taken when there are some existing or prior constraints on the frequency and/or growth rate of the eigenvectors, for example in the case of an instability where one requires the fastest growing mode. Solving for ten eigenvectors with a resolution of 25 × 20 × 20 points takes less than two hours to converge on 64 processors, although these calculations are still limited by the physical memory available on the processors. As mentioned, we are working on implementing sparse matrix-solving libraries in ECLIPS3D which will overcome this limitation.

Globally, the search for particular waves can be optimised by the combination of both a complete eigenvector, and specified or restricted eigenvector setup of ECLIPS3D. First, one would calculate a whole spectrum of low-resolution eigenvectors and identify the most interesting ones, before studying them in much higher resolution with the faster version of the code. However, for the time-independent solution, hence matrix inversion, the resolution can be much higher because inverting a matrix is a well parallelised and efficient process,. For example, for 25 × 25 × 25 points the matrix inversion takes less than an hour on 16 processors.

Finally, we stress that the symmetries of the problem can allow us to restrict our study to only one hemisphere, doubling the effective resolution with the same number of points. In this paper for example, the modes we present are always symmetric about the equator, but ECLIPS3D can solve for both symmetric and antisymmetric modes at the equator.

3. Benchmarking

We have applied ECLIPS3D to five well-studied cases from the literature to benchmark the code. These tests are explained in this section. First we reproduce the results of Thuburn et al. (2002) for a simple, hydrostatically balanced, and zonally symmetric atmosphere at rest (Sect. 3.1). This is followed by the case of an unstable jet providing a steady initial circulation as introduced by Wang & Mitchell (2014) (Sect. 3.2). We also present results for the baroclinic instability test of Jablonowski & Williamson (2006) and Ullrich et al. (2014) (Sect. 3.3). In order to implement longitudinal variation in the steady state, we study the stability of Rossby–Haurwitz waves (Haurwitz 1940), as done in Thuburn & Li (2000). Finally, ECLIPS3D is applied to the case of a linear steady-state circulation with atmospheric drag following Komacek & Showman (2016) (Sect. 3.5).

3.1. Initial atmospheric rest state

We first apply the 2D, axisymmetric version of ECLIPS3D to solve for the eigenmodes of an initially axisymmetric, isothermal, and hydrostatically balanced atmosphere at rest following Thuburn et al. (2002). Namely, the atmosphere is 80 km high with the bottom boundary at the Earth radius a = 6371 km, and the temperature is T  =  250 K corresponding to N²  =  3.83 × 10⁻⁴ s⁻². The value of the other parameters are R  =  287.05 J kg⁻¹ K⁻¹, c_p  =  1005.0 J kg⁻¹ K⁻¹, Ω  =  7.292 × 10⁻⁵ s⁻¹, and g  =  9.8062 m s⁻². This version of ECLIPS3D needs to assume an integer wavenumber m in longitude, as in Thuburn et al. (2002). Table 1 shows the frequencies of the modes from both this study and that of Thuburn et al. (2002) revealing agreement better than 3%; the discrepancies are due to slightly different initialisations and grid staggering. When matching their setup exactly we return matching results to within machine precision. Additionally, our resulting eigenfunctions have the same shape in height and latitude and global values as those found in Thuburn et al. (2002). For example, we isolate and present both an acoustic and Rossby wave recovered by ECLIPS3D in Figs. 2 and 3, to be compared to Figs. 1 and 2 of Thuburn et al. (2002). The acoustic mode shows a vertical compression mode, with little energy in the horizontal velocities and an opposite phase between the pressure and vertical velocity perturbations. The tilt in the zonal velocities close to the pole confirms the impossibility of obtaining solutions with separate functions in the latitudinal and vertical directions in spherical geometry, as noted by Thuburn et al. (2002). The Rossby modes on the other hand are dominated by pressure and horizontal velocity perturbation, and if projected on a latitude-longitude plane we would recover the rotating winds around pressure maxima or minima in the mid-latitudes. Once again the shape of the waves clearly confirms the impossibility to separate variables.

Fig. 2. Values of the five perturbed variables u′,v′,p′,w′, and θ′ obtained with ECLIPS3D for an acoustic wave with longitudinal wave number 1, with units proportional to their influence on the energy of the wave (as our solutions are from linear theory, all values are defined relative to an unknown proportionality value). This mode is to be compared to Fig. 2 of Thuburn et al. (2002).

Fig. 3. Same as Fig. 2 but for a Rossby wave, to be compared with Fig. 1 of Thuburn et al. (2002).

The difference in the forcing mechanisms between acoustic and Rossby waves, namely the pressure gradient and Coriolis force, respectively, mean we expect their global features to differ. The calculations from our energy equation, shown in Table 1, are in excellent agreement with the obtained numerical frequencies. These calculations also allow one to investigate the restoring force. For example, for the Rossby waves the f and Ff terms of Eq. (4) account for 90% of the value of σ, whereas they are negligible compared to the terms involving the pressure, and its derivative, for the acoustic waves.

Although this problem is axisymmetric, it can be used to test the 3D version of ECLIPS3D. In Fig. 4 we present a single mode from the 3D case, to be compared to Fig. 2. For this mode (and the other modes not presented explicitly here) we obtain the same height and latitude behaviour. For the additional dimension, longitude, we recover oscillatory modes with an arbitrary integer number, 2m, of zeros in longitude (corresponding to a wavenumber m). The frequency of the obtained modes again, as with the 2D version of ECLIPS3D, match those of Thuburn et al. (2002) to better than ∼3% (in this case errors are also introduced by the discretisation in longitude).

Fig. 4. Same as Fig. 2 but from the 3D version.

3.2. Unstable jet

The next benchmark case is an initial state which includes an initial velocity field. Here we follow Wang & Mitchell (2014) who identified an exponentially growing linear mode, bringing eastward momentum to the equator under axisymmetric forcing. This study essentially identifies unstable modes in an atmosphere similar to Thuburn et al. (2002) but including a mid–latitude unstable jet. Namely, the initial velocity is controlled by a given latitude ϕ₀ through:

(10)

where α controls the decay of the velocity field towards the pole. The value of α is not given in Wang & Mitchell (2014); here we chose α = 50 which mimics the shape of their initial velocity field in their Fig. 1.

Wang & Mitchell (2014) identify two instabilities, firstly a well–known baroclinic instability (such as studied in Sect. 3.3), and secondly a new instability not captured by analytical treatments under the β–plane approximation. This new mode results in the convergence of eastward momentum at the equator, and is related to the Rossby and Froude numbers. Wang & Mitchell (2014) term this new instability the Rossby–Kelvin instability as it emerges from interaction between the mid-latitude Rossby waves and the Kelvin wave (an equatorially confined gravity wave with zero meridional velocity). In Fig. 5 we show the characteristics of the mode identified by our own study using the 2D axisymmetric ECLIPS3D, which is to be compared with Fig. 1a of Wang & Mitchell (2014). Figure 5 demonstrates the excellent agreement of the structure of the mode found using both ECLIPS3D and that of Wang & Mitchell (2014).

Fig. 5. Pressure (colour scale) and wind (vector arrows) for the most unstable mode obtained with ECLIPS3D from the setup of Wang & Mitchell (2014), which is to be compared with their Fig. 1a.

Following Wang & Mitchell (2014) we explore the effect on the most unstable mode of varying the planetary parameters. Figure 6 shows results for different ϕ₀, a characteristic latitude of the initial flow (see Wang & Mitchell 2014, for definitions), and the Burger number Bu = ((N_iH)/(2Ωa))² where H is a characteristic height, revealing a change in the growth rate as H is altered (see Wang & Mitchell 2014, for more details). Figures 6a and b therefore represent the most unstable mode for a broad inital jet in latitude, up to 50°, the growth rate of which is about Ω/5 and the characteristic height a third of the total height. On the other hand, Figs. 6c and d represent the most unstable mode for a narrower inital jet, with ϕ₀ = 35°, the growth rate of which is about Ω/20 and characteristic height a fifth of the total height. The obtained growth rates are consistent with those presented in Fig. 2a of Wang & Mitchell (2014).

Fig. 6. Pressure perturbations (colour scale) and wind vectors (arrows, panels b and d only) for the case of a planet with the same radius and rotation rate as the Earth but an isothermal temperature pressure profile set at 500 K. Panels a and c: latitude against height at a longitude of 200°; panels b and d: longitude against latitude at a height of 25 000 m. We report the values of ϕ0 and Bu as defined in Wang & Mitchell (2014) and the growth rate σgrowth: (a) and (b): ϕ0 = 50°, Bu ∼ 0.2 and , and (c) and (d): ϕ0 = 35°, Bu ∼ 0.05 and . These results indicate that the growth rate of the most unstable mode is dependent on the characteristic height of the wave as found in Wang & Mitchell (2014).

3.3. Baroclinic instability

Jablonowski & Williamson (2006) detail a baroclinic instability test for GCMs using pressure as a vertical coordinate. Ullrich et al. (2014) adapted this test for height-based GCMs. In this test a perturbation to a steady longitudinal wind at mid-latitudes leads to a dynamical instability growing in a few days (for Earth-like conditions). In full 3D GCM simulations many phenomena act simultaneously, meaning that reproducing an instability with the exact same behaviour and time evolution is unlikely. However, we can expect to reproduce the most unstable modes, which will drive the evolution of the atmosphere in the simulations.

We implemented the initial state prescribed in Ullrich et al. (2014) in both the axisymmetric 2D and 3D versions of ECLIPS3D (without the prescribed perturbation of Ullrich et al. 2014, as ECLIPS3D intrinsically perturbs steady states). We only show the 2D results as, similarly to the first test case, 2D and 3D are in excellent agreement. In this test, ECLIPS3D identifies the stable modes of Thuburn et al. (2002), slightly modified by the mean flow and the angular dependency of pressure and temperature. For the unstable modes, a continuum in frequency is returned (discretised by the numerical precision of the algorithm, controlled by the number of points in the matrix B), as expected from analytical treatment. ECLIPS3D identifies the most unstable mode at m = 5 with a growth rate of 6.4 × 10⁻⁶ s⁻¹ (∼2 days) which is presented in Fig. 7. Figures 4 and 5 of Ullrich et al. (2014) demonstrate that the instability dominates the flow after 8 days, consistent with our growth rate. Additionally, the shape of our instability has similar features to the thermodynamic state of the atmosphere in Ullrich et al. (2014) after 8 days. Our instability indeed exhibits a tilt in the plot of height versus latitude (Fig. 7b) as can be seen in Fig. 6 of Ullrich et al. (2014). The pressure also exhibits a similar sharp decrease just above the surface. One must note here that our results include an uncontrolled phase in longitude coming from the axisymmetry of our setup.

Fig. 7. Pressure (colour scale, pascals) and horizontal winds (arrows panel a only) produced using ECLIPS3D, for the baroclinic instability setup of Ullrich et al. (2014). Panel a: near-surface pressure as a function of longitude and latitude. Panel b: pressure as a function of longitude and height at 50° latitude.

The only difference between 2D and 3D is the precision due to discretisation. In 3D, ECLIPS3D is limited as the size of the matrix to invert is much bigger than in 2D. Due to the shape of the staggered grid (Fig. 1), the first point in pressure is not the surface pressure and the sharp decrease is less obvious than in 2D (not shown).

Comparison between the semi-analytical calculation from our a posteriori energy equation and the numerical eigenvalue obtained for this mode reveals close agreement, on the real and imaginary part of σ, within a few percent. As explained in Appendix B, we can decompose this energy equation into three components (five in the case with meridional and vertical velocities) comprising the terms coming from the equations at rest with no angular dependency on the thermodynamic variables, the terms arising from the angular dependency of the steady state, or the terms coming from the initial zonal velocity. From analytical considerations, we expect the frequency to be dominated by the velocity terms, where the global mean flow excites the modes at specific phase velocity. However, the growth rate should be controlled by the angular-dependent terms as the baroclinic instability arises from horizontal gradients in the pressure and temperature (see e.g. Vallis 2006).

Our energy-based calculation gives a frequency of 1.09 × 10⁻⁵ s⁻¹ where ECLIPS3D finds 9.93 × 10⁻⁶ s⁻¹. In the calculation, the velocity terms account for more than 80% of the frequency, confirming the analytical predictions. The calculated growth rate is 6.63 × 10⁻⁶ s⁻¹, close to the ECLIPS3D value of 6.39 × 10⁻⁶ s⁻¹, with the angular terms accounting for 96%. These results show that our a posteriori energy equation can be a powerful tool to obtain insight into the physics of numerical eigenvectors.

3.4. Rossby–Haurwitz waves

The previous steady states we studied are axisymmetric and include no background velocity (Sect. 3.1) or only a zonal velocity (Sects. 3.2 and 3.3). Ideally, we would also benchmark ECLIPS3D using a test including meridional velocities as well as a dependency on longitude. Unfortunately, there are no such non-linear steady states in 3D, in which the analytical theory can provide us with predictions to compare with. We therefore consider a 2D non-axisymmetric problem, with steady zonal and meridional winds: Rossby–Haurwitz waves. Here we identify the most unstable modes around two steady configurations of this setup which we detail below.

Rossby–Haurwitz waves are analytical solutions of the non-linear barotropic vorticity equations. They were discovered by Haurwitz (1940) by perturbing the non-divergent equations and solving them non-linearly. If the flow remains incompressible at all times, these waves remain analytical solutions of the full equations and propagate without changing their form at constant speed. With an appropriate choice of parameters, this speed can be zero, and these waves become a stationary, steady solution of the non-divergent barotropic vorticity equations, hence another test case for ECLIPS3D.

The non-divergent barotropic equations are simply Eqs. (1a) and (1b) with w = 0 with an imposed null divergence:

(11)

This latter equation does not involve any time derivative, and therefore we have to slightly adapt the structure of the code for this setup. Instead of solving an eigenvector problem, we solve a generalised eigenvector problem where

(12)

where x is an eigenvector, B the linearised matrix of equations, and C a diagonal matrix with some zeros in the diagonal. The linearisation is straightforward, as the equations are similar to the full set of equations and the divergence equation is already linear.

The interest of this test lies in the stability analysis of these waves. Hoskins (1973) showed that Rossby–Haurwitz waves are stable for longitudinal wave number R <  5 and unstable for higher wave numbers. However, in the analysis of Hoskins some triad interactions were missing, as simplifications were required to treat the problem analytically. Inspired by Baines (1976), Thuburn & Li (2000) resolved the issue by showing that a Rossby–Haurwitz wave of wavenumber 4 is unstable, because of an interaction with wave numbers 1, 3, and 5 (see also Lynch 2009).

ECLIPS3D does not make any assumption on the shape of the perturbation needed to trigger an instability, nor on the instability itself. We therefore expect to obtain unstable modes around a steady R = 4 wave. Our setup is similar to the classical benchmarking test of Williamson et al. (1992) and Thuburn & Li (2000), with a vorticity ψ being:

(13)

where K and ω are constants and R is the longitudinal wavenumber. In order to obtain a steady, stationary wave we must also impose (see Haurwitz 1940):

(14)

With R = 4 and Ω = 7.29 × 10⁻⁵ s⁻¹, the Earth rotation rate, this leads to ω ≈ 5 × 10⁻⁶, close to the value chosen by Williamson et al. (1992) and Thuburn & Li (2000)ω ≈ 7.8 × 10⁻⁶. With R = 2, ω ≈ 1.5 × 10⁻⁵. In accordance to their setup, we impose arbitrarily K = ω.

We present the initial steady states for R = 4 and R = 2 in Figs. 8a and c. As expected from Thuburn & Li (2000), we obtain a linear instability for the R = 4 setup, displayed in Fig. 8b. This instability oscillates with a period of ≈3 days, with an exponential growth timescale of 6 days. The timescale for instability is globally coherent with the results of Thuburn & Li (2000), which find that the flow becomes significantly altered after day 20. For R = 2, we also find an instability with striking resemblance to the R = 4 instability in shape, a period of just under a day, and growth timescale of 4 days. This mode is shown in Fig. 8d, and this result is in contradiction with Hoskins (1973) but in accordance with Baines (1976) who shows analytically that all the R ≥ 2 (n ≥ 3 in their study, where we have R = n − 1 here) Rossby–Haurwitz waves can be unstable. Interestingly, we also found an instability for a R = 3 Haurwitz wave but with a growth timescale of more than a hundred days. Such an instability would probably be smoothed out by any source of dissipation or diffusion in a GCM. Globally, our code agrees well with the theoretical study of Rossby–Haurwitz waves, and demonstrates the proper treatment of longitudinal-dependent steady state and meridional velocities in ECLIPS3D.

Fig. 8. Pressure (colour scale, a and c in pascals and arbitrary units for b and d) and horizontal winds (arrows) for Rossby–Haurwitz waves and most unstable modes. Panels a and c: initial steady states for R = 4 and R = 2, respectively (see text). Panels b and d are the most unstable mode obtained with ECLIPS3D for the R = 4 and R = 2 setup, respectively.

3.5. Linear steady circulation with drag

As discussed in the introduction, our development of ECLIPS3D was largely driven by studies of the acceleration of zonal flows in hot Jupiter atmospheres. For these planets, analytical studies have shown linear steady states to be of vital importance (see Showman & Polvani 2011). Therefore, we have also implemented the capability to calculate linear steady states in ECLIPS3D (see Sect. 2), which is a much simpler process compared to the identification of eigen modes. To benchmark this section of the code, we compare our results to those obtained in the study of Komacek & Showman (2016), in particular the case they present in their Fig. 5 where they solve the full Navier-Stokes equations, but with such a low heating rate that only the linear terms contribute. This requires the addition of a linear drag in the linearised equations following the depth-dependent behaviour of that adopted by Komacek & Showman (2016). Additionally, the heating is performed via a Newtonian relaxation with a height-dependent radiative constant. The resulting equation set is shown in Appendix A.5. Figure 9 then presents the resulting linear circulations obtained using ECLIPS3D which show excellent qualitative agreement with the results of Komacek & Showman (2016; see their Fig. 5). Komacek & Showman (2016) do not present the vertical structure of their circulation.

Fig. 9. Pressure (colour scale, Pascals) and horizontal winds (arrows) as a function of latitude and longitude at a height of 5 × 106 m, for a steady-state circulation with a forcing of ΔT = 100 K (see Komacek & Showman 2016, for definition). Results to be compared to Fig. 5 of Komacek & Showman (2016). (a) τdrag = 105 s and τrad = 104 s. (b) τdrag = 103 s and τrad = 103 s. (c) τdrag = 106 s and τrad = 104 s. (d) τdrag = 104 s and τrad = 106 s.

4. Conclusion

In this paper, we introduce and benchmark ECLIPS3D: a parallel code for identifying linear instabilities, waves, and circulations around a steady state of the Navier–Stokes equations in planetary atmospheres. The linearised equations only omit viscosity and an a posteriori energy equation is used to identify contributions from each component. The time-dependent eigenvector solution or time-independent matrix inversion calculations are performed through discretisation onto a staggered grid and subsequently ScaLAPACK routines.

The benchmarks cover various well-studied wave and instability tests, namely a simple atmosphere at rest (Thuburn et al. 2002), a Rossby–Kelvin unstable jet (Wang & Mitchell 2014), a baroclinically unstable jet (Ullrich et al. 2014), an unstable Rossby–Haurwitz wave (Thuburn & Li 2000), and a linear circulation with atmospheric drag (Komacek & Showman 2016). For all these setups, ECLIPS3D is able to produce excellent qualitative agreement with the previous works. We demonstrate that our a posteriori energy equation is a viable tool to verify the results and identify the dominant terms. We are currently preparing a follow-up study to explore the momentum transfer in hot Jupiter atmospheres and explore the stability of the initial conditions for GCMs using ECLIPS3D (Debras et al., in prep.).

ECLIPS3D currently has several limitations, primarily its computational efficiency, leading to limitations on resolution, particularly for 3D cases. We are working on several methods to improve this issue for example using libraries adapted to sparse matrices, or splitting the eigenvector solution into several sub-matrices as opposed to a single large matrix (potentially useful as the time taken to solve this type of problem increases faster than linearly with matrix size). This splitting of the matrix may be particularly well suited to a spectral decomposition as we are searching for the most unstable mode, and are not necessarily trying to capture the entire “shape” of the mode. This could be done through spherical harmonics in the horizontal or Chebyshev’s spectral decomposition in the radial direction.

Despite its limitations, ECLIPS3D in its current version still represents a powerful resource which can be used to study instabilities in 2D situations under axisymmetry or cases where a two–layer model is applicable, or for low-resolution 3D problems. The code itself can easily be adapted to different situations, in spherical coordinates, with additional physics or alternative boundary conditions. As the structure of the code is independent of the underlying equations, meaning alternatives can easily be implemented in terms of symmetries and coordinate systems.

Finally, ECLIPS3D could be applied to a wide range of astrophysical problems. The most obvious one, for which ECLIPS3D was designed, is the study of instabilities and linear circulations for planetary atmospheres, but the range of applicability is greater. Asteroseismology for example requires the need to linearise the equations of motion and identify the leading modes, sometimes with complicated circulation or thermodynamic state inside the star. Adapted to cylindrical geometry, ECLIPS3D could be a powerful tool to identify the possible instabilities in protoplanetary disks, where instabilities creating pressure traps are proposed to be an efficient way of making planets through core accretion. The addition of a magnetic field in the equations implemented in the code would not pose any theoretical challenge either, which could provide a significant amount of information on the linear behaviour of astrophysical fluids in more general cases.

Acknowledgments

FD is indebted to the ScaLAPACK team for both the provision of their libraries and their responses to numerous questions. FD thanks the European Research Council (ERC) for funding under the H2020 research & innovation programme (grant agreement #740651 NewWorlds). NJM is part funded by a Leverhulme Trust Research Project Grant and partly supported by a Science and Technology Facilities Council Consolidated Grant (ST/R000395/1), both of which we gratefully acknowledge. This study uses material produced using Met Office Software. Additionally, used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National e-Infrastructure. Additionally, this research made use of the ISCA High Performance Computing Service at the University of Exeter. Finally, this work was also partly funded by the ERC grant No. 320478-TOFU.

References

Appendix A: Equations in ECLIPS3D

A.1. General 3D case

In this appendix we detail the derivation of the full linearised equations for the longitudinal component of the momentum equation, and in the interests of brevity provide only the final expressions for the remaining components. These equations assume a dependence of g on r as g ∝ 1/r², and every other quantity is dependent on the three spatial variables, longitude, latitude, and radial distance from the centre of the planet. Here, we consider that Q is equal to zero for simplicity, and relax this assumption in Appendix A.4.

The longitudinal equation of momentum is as follows.

(A.1)

Term by term we obtain (refer to Eq. (2a) for definition of perturbed variables)

Therefore, the final five perturbed equations are

(A.2)

(A.3)

(A.4)

(A.5)

(A.6)

with ∂X′/∂t = iσX′.

A.2. Two-dimensional axisymmetric case

For the axisymmetric case, the equations are directly obtained from the 3D case by choosing a longitudinal wavenumber m such that X′(t, r, ϕ, λ) = X′(r, ϕ)e^iσte^i(m/2π)λ.

A.3. Two-layer equivalent depth

The reference for this particular case can be found in Showman & Polvani (2011) or Vallis (2006). We consider a dynamic layer above a quiescent layer, reservoir of mass or energy, and study the horizontal winds u and v as well as the height of the layer h, depending on both x and y, the cartesian horizontal coordinates. We follow the definitions of Showman & Polvani (2011) for the variables, hence consider adimensional equations. There is consequently only three equations to be implemented:

(A.7)

(A.8)

(A.9)

where H = H(x, y) is the initial steady height.

A.4. Heating rate

Particular care must be taken when dealing with the heating rate. If we call the heating rate Q_i, as the initial state is steady the zeroth order term will cancel the terms involving Q_i in Eq. (1e). However, two situations must be considered: the θ/T factor in Eq. (1e) has to be linearised, and will be a source of additional terms. Additionally, if Q_i depends on the atmospheric state (for example with Newtonian heating, see next appendix), a perturbation in the atmosphere will be associated with a change Q′ in Q_i. Therefore, if we write equation Eq. (A.6) as

(A.10)

where the gρ_i/θ_i factor arises from the definition of θ′, the final equation involving the heating rate is

(A.11)

where Q′, if it exists, depends linearly on the linearised atmospheric variables.

Moreover, obtaining Eq. (A.5) implies to use Eq. (1e), and therefore additional terms also have to be included. More precisely, one could show that Eq. (A.5) can be written as

(A.12)

With the Q terms we obtain

(A.13)

These new terms in Eqs. (A.5) and (A.6) have to be implemented in the matrix from which we solve for eigenvectors, but do not lead to a change in the way of finding the eigenvectors.

A.5. Steady linear circulation

Following Showman & Polvani (2011) and subsequently Komacek & Showman (2016), we have implemented the possibility to solve for linear steady states instead of waves and instabilities. We therefore have to impose a heating of the atmosphere, associated to dissipative processes in order to reach a steady state.

This heating function is extremely different from the heating of Appendix A.4. In Appendix A.4, we linearised the heating term coming from the initial steady solution of Navier Stokes equations. Here, we prescribe a small forcing of the atmosphere that causes it to depart from its initial steady state, and seek the new steady state that the atmosphere will reach at the linear order (because the heating has a small amplitude). For simplicity reasons, we consider that the initial steady state was obtained without forcing of the atmosphere (hence Q in Eqs. (1e), (A.11) and (A.13) is identically null), and use Q_l to refer to the low-amplitude linear forcing we impose.

In that case, the perturbed variables are assumed to be constant with time (σ is taken to be zero). The dissipative effects will simply be linear drags in Eqs. (1a)–(1c) expressed as −v/τ_drag where τ_drag is a characteristic time for the drag, eventually dependent on the space coordinates (see Showman & Polvani 2011).

If Q_l is constant, then we simply have to modify Eqs. (A.5) and (A.6) in a similar way to in Appendix A.4:

which, as Q_l is order 1, simply consist in neglecting the second-order terms in Eqs. (A.13) and (A.11) and moving the constant heating terms to the right-hand side. A dissipative or diffusive process could also be added in the energy equation.

Additionally, a special case must be discussed: Newtonian Heating (see e.g. Mayne et al. 2014b). In that case, Q_l is not constant but depends on the thermodynamic state of the atmosphere. More precisely, using Q_N to denote the Newtonian heating rate:

(A.14)

where T_eq is a prescribed equilibrium temperature and τ_rad a characteristic radiative time, both depending on space variables.

For the linear forcing approximation to remain correct, T_eq must be sufficiently close to the initial temperature T_i, but then a small change in T_i will have an impact on Q_N of the same order of Q_N itself. With our choice of perturbed variables, it is easy to show that

(A.15)

and subsequently

(A.16)

Finally Eqs. (A.5) and (A.6) can be rewritten as (using )

(A.17)

(A.18)

To summarise, when looking for a steady linear circulation with Newtonian heating, we need to invert the matrix C as advertised in Sect. 2.4.2, where C arises from Eqs. (A.2)–(A.6) and includes the heating rates and dissipations expressed in the text of this appendix, and in Eqs. (A.17) and (A.18).

Appendix B: A posteriori energy equation

In order to obtain a semi analytical verification for the frequency, we integrate the energy of the modes over the whole volume, and express it as an a posteriori condition on the frequency σ. In this part, we assume that the bottom boundary condition is a no-escape condition (w′ = 0) and that the initial state is in the hydrostatic balance: , with no initial heating (see Appendix A.4). These assumptions could be relaxed, but would be sources of numerous additional terms whereas they are always verified in our setups.

In the 2D axisymmetric case at rest with no angular dependency in the initial variables, Thuburn et al. (2002) used as variables u′, −iv′, −iw′, p′ and θ′ because this simplifies greatly the calculation. In order to allow for easier verification of our equations, we adopt the same definition for the perturbed variables. However, for simplicity reason, we drop the primes in the following equation and use v and w, not iv′ and iw′. Therefore, one has to remember that the v and w expressed in the following equations are actually −iv′ and −iw′ where v′ and w′ are the solutions of Eq. (8). The other variables are not affected.

Denoting a complex conjugate by a star, we express the integral of energy as

(B.1)

where Ω is the whole volume, dV = r²cosϕdrdϕdλ, which is the infinitesimal volume, and (A.2) is the left-hand side of the complete Eq. (A.2) and so on.

The calculations are highly cumbersome, but present no particular difficulty. In order to have a physical insight in the leading mechanism from this a posteriori energy equation, we have decided to separate this integral into five parts:

• The first part involves only the thermodynamic initial state (no velocities) with a dependency on the radial variable r solely. An initial atmosphere at rest with no angular dependency would have contributions to the energy only from this part.

• The second part involves the terms coming from the angular dependency in the thermodynamic steady variables only.

• The third part comes from the steady zonal velocity u_i.

• The fourth part is generated by the steady meridional velocity v_i.

• And finally the last part is due to the initial steady vertical velocity w_i.

Denoting this decomposition of the energy integral as [1]−[5], and remembering ∂/∂t = −iσ we obtain an a posteriori equation on σ:

(B.2)

This is similar to Eq. (4). The 1/ρ_i factor might seem useless as it is already in Eq. (B.1) but is a necessary density weighting to obtain the appropriate equations. E is unchanged:

(B.3)

After sorting (real component, then imaginary then complex), the calculation gives

(B.4)

(B.5)

(B.6)

(B.7)

(B.8)

All Tables

All Figures

Fig. 1. Example cell of the 3D staggered grid adopted in ECLIPS3D, based on Thuburn et al. (2002). i discretises the longitudinal variable λ, j the latitude ϕ, and k the radial variable r. Here, u and v are staggered in latitude, with v running from the south to north pole and thereby having an additional latitude point. Also, p is staggered in height with w and θ with the latter two variables running from the bottom to the top of the atmosphere resulting in an additional height point.

In the text

	Fig. 2. Values of the five perturbed variables u′,v′,p′,w′, and θ′ obtained with ECLIPS3D for an acoustic wave with longitudinal wave number 1, with units proportional to their influence on the energy of the wave (as our solutions are from linear theory, all values are defined relative to an unknown proportionality value). This mode is to be compared to Fig. 2 of Thuburn et al. (2002).
In the text

	Fig. 3. Same as Fig. 2 but for a Rossby wave, to be compared with Fig. 1 of Thuburn et al. (2002).
In the text

	Fig. 4. Same as Fig. 2 but from the 3D version.
In the text

	Fig. 5. Pressure (colour scale) and wind (vector arrows) for the most unstable mode obtained with ECLIPS3D from the setup of Wang & Mitchell (2014), which is to be compared with their Fig. 1a.
In the text

Fig. 6. Pressure perturbations (colour scale) and wind vectors (arrows, panels b and d only) for the case of a planet with the same radius and rotation rate as the Earth but an isothermal temperature pressure profile set at 500 K. Panels a and c: latitude against height at a longitude of 200°; panels b and d: longitude against latitude at a height of 25 000 m. We report the values of ϕ0 and Bu as defined in Wang & Mitchell (2014) and the growth rate σgrowth: (a) and (b): ϕ0 = 50°, Bu ∼ 0.2 and , and (c) and (d): ϕ0 = 35°, Bu ∼ 0.05 and . These results indicate that the growth rate of the most unstable mode is dependent on the characteristic height of the wave as found in Wang & Mitchell (2014).

In the text

	Fig. 7. Pressure (colour scale, pascals) and horizontal winds (arrows panel a only) produced using ECLIPS3D, for the baroclinic instability setup of Ullrich et al. (2014). Panel a: near-surface pressure as a function of longitude and latitude. Panel b: pressure as a function of longitude and height at 50° latitude.
In the text

	Fig. 8. Pressure (colour scale, a and c in pascals and arbitrary units for b and d) and horizontal winds (arrows) for Rossby–Haurwitz waves and most unstable modes. Panels a and c: initial steady states for R = 4 and R = 2, respectively (see text). Panels b and d are the most unstable mode obtained with ECLIPS3D for the R = 4 and R = 2 setup, respectively.
In the text

Fig. 9. Pressure (colour scale, Pascals) and horizontal winds (arrows) as a function of latitude and longitude at a height of 5 × 106 m, for a steady-state circulation with a forcing of ΔT = 100 K (see Komacek & Showman 2016, for definition). Results to be compared to Fig. 5 of Komacek & Showman (2016). (a) τdrag = 105 s and τrad = 104 s. (b) τdrag = 103 s and τrad = 103 s. (c) τdrag = 106 s and τrad = 104 s. (d) τdrag = 104 s and τrad = 106 s.

In the text