Issue 
A&A
Volume 659, March 2022



Article Number  A43  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202141638  
Published online  03 March 2022 
Mass transport in a moist planetary climate model
National Space Institute, Technical University of Denmark,
Elektrovej,
2800,
Kgs. Lyngby, Denmark
email: joao.mendonca@space.dtu.dk
Received:
25
June
2021
Accepted:
6
October
2021
Planetary climate models (PCMs) are developed to explore planetary climates other than that of the Earth. Therefore, the methods implemented need to be suitable for a large diversity of conditions. Every planet with a significant atmosphere has condensible cycles (e.g. the hydrological cycle), which can play an essential role in the planet’s appearance and environment. We must accurately represent a condensible cycle in our planet simulations to build a powerful planetary climate predictor. OASIS is a 3D PCM capable of selfconsistently representing the main physical processes that drive a planet’s environment. In this work, we improve the representation of mass transport in OASIS, which is the first step towards a complete and flexible implementation of a condensible cycle. We implement an upwindbiased scheme on a piecewise linear approximation with a flux limiter to solve the mass transport equation. We first benchmark the new scheme on a 2D problem that confirms the superior properties of the new method over the central finitevolume method in terms of performance, accuracy, and shapepreserving mass distribution. Due to the new scheme’s less dispersive nature, we do not have to apply any unphysical diffusion to maintain the model stable. OASIS includes the new improved solver in the total mass and the tracer (e.g. clouds and individual gas chemical species) transport. We couple the new formulation with physical schemes and validate the new code on two 3D simulations of an ocean Earthlike planet and an ocean tidally locked planet. The new OASIS simulations are robust and do not show any known problems from the dynamicsphysics coupling. We show that the two simulations capture the main characteristics of ocean planet atmospheres and are easy to set up. We propose these two simulations as the first standard benchmark tests for models built to explore moist planetary environments.
Key words: hydrodynamics / planets and satellites: atmospheres / planets and satellites: terrestrial planets
© ESO 2022
1 Introduction
Planetary climate models (PCMs) are flexible numerical tools capable of representing vast planet conditions. Every detail on how equations are solved in the models has repercussions on the model output and modelled physics. In this work, we focus on the atmospheric mass transport problem and simple moist processes.
A gas can exist in multiple phases in a planet environment. It can, for example, condense and form clouds or rain and fall onto the planet’s surface. In this work, we define the general term condensible cycle as the ‘life’ cycle of a condensible gas in a planet’s environment. Any planet with a significant atmosphere is expected to host a least one condensible gas. The composition and spatial distribution of the condensible material are tightly connected with the context of the planetary environment. To have a deep understanding of how a climate of a particular planet works, we need to study the impact of all the possible condensible cycles for the specific planet. On Earth, it is the water cycle that plays a vital role in the planet’s climate. The condensible gases are expected to be diverse under different planet conditions, for example water on Earth, sulphuric acid on Venus (e.g. Taylor 2011), and methane on Titan (e.g. Hörst 2017). By having a robust representation of moist processes in the atmosphere, we improve our model predicting power in planetary climates. In this work we present the first steps towards a complete implementation of a condensible cycle in the OASIS model (Mendonça & Buchhave 2020).
Atmospheric circulation transports gas and clouds across the atmosphere. A robust and fast representation of the atmospheric transport in 3D simulations is essential to avoid compromising the physics in the simulations. Inadequate atmospheric transport models can result in a loss of mass conservation, the creation of artificial highfrequency waves, or a poor representation of how chemistry and clouds are redistributed across the atmosphere. Condensible constituents of the atmosphere interact in various ways with the radiation. The distribution of the components can have a substantial impact on the atmosphere’s energy budget (e.g. the planet’s bond albedo; Read et al. 2016). The heat absorbed or released during the different phase changes also impacts atmospheric thermodynamics and has an important role in the planet’s atmospheric circulation.
Some phenomena that involve complex representations of moist processes remain poorly explored by 3D simulations, such as atmospheric collapse (e.g. Wordsworth 2015) and runaway greenhouse phenomena (e.g. Ingersoll 1969). These phenomena require the representation of a complete condensible cycle, such as cloud microphysics and condensate transport, which can be the main component of the atmosphere. The thermodynamical equations in the model’s dynamical core need to capture the leading physics and avoid including assumptions such as those typically adopted from Earth climate models to facilitate the numerical schemes. One such assumption is the limit of small condensable content in the atmosphere compared to the main atmospheric component (Ullrich et al. 2017), such as water and nitrogen on Earth. The limit of the small condensable abundance simplifies the vertical integration of the equations in the dynamical core of Earth models since quantities such as gas constant and heat capacity remain fixed.
Clouds or hazes in terrestrial planets can mask the atmosphere below and challenge the detection of molecules in the atmosphere (e.g. Fauchez et al. 2019; Suissa et al. 2020; Komacek et al. 2020), which highlights the need to have a robust aerosol distribution in 3D numerical simulations to maximise the science output of observational campaigns. We need to have a better understanding of the limitations of our models. Poor representations of these atmospheric phenomena in the simulations can lead to a flawed interpretation of the observational data and analysis of unphysical phenomena from the 3D simulations. The details on how the clouds and chemistry are being transported in 3D models of planetary simulations are often left out of scientific papers, making it hard to compare our model results with other models applied to extrasolar planets. Every formulation has pros and cons, and 3D planetary climate users or developers should make more explicit the options adopted in their simulations. This work aims to give the first steps towards developing a complete moisture cycle in a 3D dynamical core, focusing on the problem of mass transport and latent heat release. We explain how the modeller can estimate the numerical schemes’ errors and identify the primary source of numerical errors, which is essential for future model intercomparison studies. By solving the total mass and tracer equations with the same solver, we maintain consistency between the atmospheric mass and any tracer in the atmosphere in terms of spatial distribution and accuracy, which becomes important for planetary conditions where the condensable gas has a large concentration in the atmosphere. As far as we know, in the planetary community only the model developed in Ge et al. (2020) can follow a similar approach due to the higherorder scheme implemented to solve the mass conservation equation.
In the next section we describe the new equations and numerical schemes implemented in the OASIS dynamical core (known as THOR); we also explore the accuracy of the new advection scheme. In Sect. 3 we present the simplified physics implemented in the 3D simulations, and in Sect. 4 we present the results of the 3D ocean planet benchmark tests. Concluding remarks are presented in Sect. 5.
2 New equations in the OASIS dynamical core
OASIS is a 3D virtuallab to study planetary environments that has been written from groundup (Mendonça & Buchhave 2020). This platform comprises different physics and chemistry modules that interact with each other to reproduce selfconsistent results. OASIS makes use of the flexible dynamical core, THOR, which has a 3D nonhydrostatic nature aimed at exploring a large diversity of planets (Mendonça et al. 2016 and Deitrick et al. 2020). One example of the OASIS robustness and success is the recent successful OASIS simulations of Venus (Mendonça & Buchhave 2020). The climate of Venus is very challenging to simulate due to its highly reflective clouds that entirely cover the planet and its massive atmosphere that requires tens of thousands of simulation days (e.g. Lebonnois et al. 2010 and Mendonça & Read 2016). Venus has strong winds that cause the atmosphere to rotate 60 times faster than the solid planet, in a phenomenon called ‘superrotation’ (e.g. Read 1986; SánchezLavega et al. 2017). The strong winds have an enormous impact on the planet’s climate. Using our stateoftheart PCM, OASIS, we were able to simulate the extreme conditions observed in the Venus cloud region (Mendonça & Buchhave 2020).
In this work, we mainly expand the scope of the atmospheric dynamics module (THOR) to represent the atmospheric transport in a more accurate and physically based way to couple the dynamical core with simple physical representations of the atmospheric water cycle.
2.1 Expanding THOR equations
THOR is the module that solves the resolved atmospheric fluid flow (Mendonça et al. 2016). THOR was built to be flexible and efficient in exploring planetary atmospheres, and the equations solved for dry atmospheres are presented in Mendonça et al. (2016) and Deitrick et al. (2020). The equations of motion are solved in a modified icosahedral grid (Tomita et al. 2001). In this work, we extend and modify the code THOR to include: (i) a tracer transport equation explicitly in the main set of equations of the dynamical core and (ii) improved consistency and accuracy of the mass transport representation.
The new set of equations solved in the dynamical core are the following:
where ρ is the atmospheric density, v is the velocity, P is the atmospheric pressure, g is the planet gravity, Ω is the planet’s rotation vector, θ is the potential temperature, T is the atmospheric temperature, c_{p} is the heat capacity of the atmosphere, Q is the heating or cooling from radiation, convection, condensation and evaporation, q_{i} is the concentration of a tracer with index i (e.g. clouds, chemical gases), and is the physical source or loss of the tracer i. In our new platform the prognostic variables are: ρu_{h}, ρw, ρ, ρθ, and ρq_{i}. In addition to the equations above, we assume that the ideal gas equation and the dry atmosphere’s thermodynamic equations are valid anywhere in the atmosphere. We have not further explored the impact of moist thermodynamic equations because we want to focus only on the impact of the atmospheric mass transport problem in this work. A complete coupling of the thermodynamic equations is planned for followup work. The atmosphere’s specific heat capacities and gas constants vary for composition, phase changes, and temperature variations in a moist atmosphere. The equations or methods have to be adapted to the needed atmospheric conditions. Recent Venus simulations (e.g. Lebonnois et al. 2010; Mendonça & Read 2016) are good examples on how variations in the specific heat capacity as a function of temperature can be incorporated in 3D PCMs. In addition, planetary conditions where the main atmospheric gas is condensing will require more complete equations than those used in Earth climate models.
As explained in Mendonça et al. (2016) and Deitrick et al. (2020), THOR discretises the equations horizontally in an Arakawa Agrid and the vertical levels on a nonstaggering grid. The spatial integration of the momentum and entropy equations are solved using the central finitevolumes methods from Satoh (2002), and Tomita & Satoh (2004). In this work, we upgrade the spatial integration of the density and tracer equations with a new upwindbiased scheme with a flux limiter (Miura 2007), as explained in the section below. The density and the tracers are solved with the same solver to maintain physical consistency.
The time integration is explained in Mendonça et al. (2016). Our numerical scheme involves splitting the time integration into short explicit steps for the terms associated with acoustic modes and large steps for the rest (e.g. Wicker & Skamarock 2002; Skamarock & Klemp 2008). The fast acoustic waves can compromise the stability of the model. The equations are also formulated to solve the vertical component of the momentum using an implicit method (HEVI method – e.g. Tomita & Satoh 2004). The implicit method in the vertical direction is necessary due to the often much higher spatial resolution in the vertical direction compared to the horizontal spacing.
2.2 Tracer transport
The tracer^{1} transport scheme is one of the most important routines in a climate model dynamical core. Equation (4) shown above represents the tracer transport. The right side term of the equation, S_{q}, is the physical source or sink of the tracer (e.g. condensation and evaporation). If S_{q} is zero, the tracer concentration is conserved. There are many ways of solving the tracer transport equation, and in this work, we improved the existing tracer transport in the dynamical core THOR (Mendonça et al. 2016). The original scheme implemented in THOR to represent the tracer transport is a central finitevolume scheme from Tomita et al. (2001). This scheme is well known to be dominated by numerical dispersion that creates highfrequency patterns in the solutions (e.g. Durran 1999). The nonphysical solutions could create maxima and minima artificially and hence produce nonphysical scenarios (e.g. rain produced by an overshoot in the cloud transport scheme). To overcome this numerical inaccuracy, we have in Mendonça et al. (2018) advected chemical species in the atmosphere of hot Jupiter planets with a hyperdiffuse term in the tracer transport equation. The diffuse term removes the highfrequency oscillations created artificially by the central finite volume. The hyperdiffuse term in the tracer transport helps to maintain the simulation stable, but the formulation is not physicalbased for mass transport. The new flux divergence operator implemented here was developed in Miura (2007). The new method is an upwindbiased advection scheme on a piecewise linear approximation, which we also complement with a flux limiter developed in Thuburn (1996). These two methods allow us to achieve two essential properties in the tracer transport scheme: mass conservation (from the finitevolume scheme) and shapepreserving mass distribution (from the fluxlimiter method). The fluxlimiter method removes, forexample, over or undershoot solutions, and is important to ensure the monotonicity of the gradient concentration.
To solve Eq. (4), we first apply the finitevolume discretisation from Miura (2007): (5)
where ρ is the atmospheric density, q the tracer concentration, A_{0} is the control cell area, Δt is the large time step in the dynamical core, l is the length of the pentagon or hexagon sides, v is the wind velocity, N_{s} is the number of sides of the pentagons or hexagons, and n is the normal vector at each side of the control volume. The indices 0 in Eq. (5) refers to quantities defined at the centre of the control volumes and R at cell face centres. Miura (2007) provides more details about the equations and discretisation. The upwind scheme alone does not maintain the monotonicity of the tracer distribution. To guarantee the monotonicity we apply the flux limiter from Thuburn (1996). The same procedure was used in Miura (2007). All the details on how to implement the flux limiter used in this work are explained in Thuburn (1996).
To guarantee that the new atmospheric transport solver is well implemented in THOR we perform the test case 1 described in Williamson et al. (1992). The idea of the test is simple, where a cosine bell shape distribution of a tracer is transported across the globe until its initial position. The simplicity of the test allows us to easily find mistakes in our implementation. The nondivergent wind that transports the tracer in our test is set to a solid body rotation,
where u and v are the zonal and meridional wind directions, respectively, ϕ is the longitude and λ is the latitude. The variable u_{0} is defined as: (8)
The constant R_{planet} is the radius of theplanet and P_{bell} the period that the cosine bell shape takes to do a complete revolution around the planet. We set these two constants to a radius of 6.371 × 10^{6} m and 12 days period. The cosine bell function was constructed following Williamson et al. (1992) work: (9)
where r is the great circle distance between any point and the centre of the cosine bell. The variable h represents the concentration of an arbitrary tracer. The initial centre for the cosine bell is (longitude, latitude) = (270°, 0°). Figure 1 shows a map of the initial cosine bell shape. Other constants in Eq. (9) were set to h_{0} = 1000 m and R = R_{planet}∕3. Other configurations of the nondivergent wind and initial position of the cosine bell distribution were explored with no significant impact in the final conclusions about the solvers accuracy.
For this simple test we implemented an explicit Euler time integration. The timesteps used for the different spatial resolutions and techniques are shown in Table 1. The maps in Fig. 2 show the difference between the initial cosine bell and after a complete revolution using different methods. In Fig. 2a, the wavy pattern due to the dispersive nature of the central finitevolume method is very clear, despite using a much shorter time step than the one used by the other methods. The central finitevolume developed in Tomita et al. (2001) is not an upwind solver, which makes the solutions less diffusive and more dispersive. To mitigate the instabilities created by the central finite volume, we could add explicit diffusion to smear the wavy pattern. However, solving the diffusion terms reduces the solver’s accuracy and makes it more timeconsuming since it implies solving Laplacians in an icosahedral grid. The large amplitude of the error shows that this method is not stable without a diffusion term. Maps 2b and c show the results for the case of the upwindbiased scheme with and without the flux limiter, respectively. If the flux limiter is not used, the amplitude of the error is smaller than using the flux limiter. However, the error is not monotonically decreasing from the centre of thecosine bell distribution, which is an indication that the shape of the distribution is not well preserved. In Fig. 2c, the flux limiter is applied, and despite the error being relatively large at the centre of the cosine bell shape, the shape of the distribution is preserved. The larger error in Fig. 2c is related to the intrinsic diffusion of the fluxlimiter scheme.
The plot in Fig. 3 shows the final equatorial tracer distribution for the simulations using the upwindbiased scheme withthe flux limiter for different spatial resolutions. We verify that by increasing the spatial resolution, the numerical diffusion becomes smaller and the numerical solutions are close to the exact solutions.
To measure quantitatively the accuracy of our schemes, we computed the error, l_{∞}, defined as (10)
where h and h_{exact} are the numerical integrated and the analytic solutions, respectively. l_{∞} is the largest error anywhere in our spherical grid. We consider l_{∞} to be the real error of the methods applied, which is not mitigated by an averaging method (e.g. root mean square error). In Fig. 4 we show that the upwindbiased scheme improves the accuracy of the numerical solver considerably without the need for extra diffusion components. The slopes of the upwind schemes are similar to the exact secondorder convergence. In some regions, the solver without using the flux limiter is even better than secondorder accuracy. For the central finite volume case we found as expected that the error increases with the spatial resolution, which is due to the dispersive nature of the solver. To improve the accuracy of using the central finite volume, as mentioned above, we can reduce the time step or apply explicit numerical diffusion.
Fig. 2 Latitudelongitude maps of the difference between the initial cosine bell and after a complete revolution across the sphere (spatial resolution glevel 6). The three maps refer to three different methods that used different time steps (Table 1): CFV (central finitevolume), UW without FL (upwind biased without flux limiter), and UW with FL (upwind with flux limiter). 
Time step and spatial resolutions for the cosine bell test cases using different transport schemes.
Fig. 3 Final tracer distributions at the equator from the cosine bell test case that uses the upwindbiased scheme with a flux limiter for different spatial resolutions. The solid black line is the analytical solution used for the initial condition (Eq. (9)). 
Fig. 4 l_{∞} convergence plot where glevel represents the different spatial resolutions. The various solid lines represent different methods: the central finite volume (red); the upwindbiased with flux limiter (green); and upwind biased without flux limiter (blue). The solid black line is a reference line with a slope of an exact secondorder convergence. 
2.3 3D transport and consistency in the atmospheric mass transport
The previous section described how the transport is calculated in the horizontal projection (parallel to the spherical surface). For the vertical integration, we used a secondorder central finite difference approach with the flux limiter from Thuburn (1996).
We implemented the tracer transport in the 3D dynamical core. The time integration was done in the large time step integration section (see Mendonça et al. 2016 for more information about THOR’s time integration). The tracer equation is integrated implicitly using the updated momenta values. We replaced the central finite volume withthe upwindbiased scheme to solve the mass continuity equation to maintain consistency between the tracers and the total atmospheric mass. This way, we maintain consistency in the spatial distribution of the various types of mass transported in the model. Other quantities, such as the momenta and potential temperature, are still integrated with the central finite volume followed with explicit hyperdiffusion as explained in Mendonça et al. (2016).
3 Benchmark tests
3.1 Simplified physics
To test our new 3D formulation, we explored first the ocean Earthlike test case proposed in Thatcher & Jablonowski (2016). The test is a modified version of the wellknown HeldSuarez test (Held & Suarez 1994). The new test includes simplified moist processes, boundarylayer mixing, latent heat exchange, and sensible heat between the planet surface and atmosphere. The moist physics schemes included have been shown to represent quantitatively well the climate state of simulations with more complex physical methods (Thatcher & Jablonowski 2016). In the section below, we summarise the main points of the physical schemes implemented from Thatcher & Jablonowski (2016) to make the reader more familiarised with the new routines implemented. More details on how to implement these schemes are available in Thatcher & Jablonowski (2016). Thatcher & Jablonowski (2016) provides supplementary FORTRAN routines with the simplified physics implemented, which we converted to the CUDA C programming language and to work on our specific grid indexing. The inclusion of simplified physics allows us to focus on the performance of the new methods implemented in OASIS to represent mass transport in a moist atmosphere.
It is well known that sudden adjustments of local atmospheric quantities due to moist processes can trigger undesirable largescale gravity waves (Thatcher & Jablonowski 2016). The 3D numerical experiments proposed in this work allow us to study if our numerical schemes are robust enough to prevent undesirable largescale gravity waves from being triggered under different planetary conditions (tidally and nontidally locked configurations).
Using a similar approach, we adapted the routines used to study the ocean Earthlike planet to explore planets under a tidally locked configuration. The two test cases (ocean Earthlike planet and ocean tidally locked planets) proposed in this work allowed us to verify that the new OASIS can reproduce robust results under two very distinct planet scenarios where moist processes play an important role in the atmospheric dynamics.
Below, we describe briefly the physical schemes implemented in OASIS to reproduce the two test cases. The values for the physical parameters used in the simulations are shown in Table 2.
3.1.1 Largescale precipitation
The physics modules for the benchmark test do not include a cloud scheme formulation. The gas that condenses in the atmosphere is removed instantaneously from the atmosphere without considering any atmospheric transport of the condensed material or reevaporation processes. The temperature and specific humidity trends during condensation are defined as (11) (12)
where L is the latent heat of vaporisationat 0 °C, c_{p} is the specific heat of the dry air, and C is the condensation rate. The two physical trends defined in these two equations, Eqs. (12) and (11), contribute to the physical trends in Eqs. (3)–(4) solved in the dynamical core. The condensation is defined by (13)
where Δt is the dynamical time step, q is the water vapour mass mixing ratio, T is the atmospheric temperature, R_{q} is the gas constant of the condensible gas (water), and q_{sat} is the saturationspecific humidity. The q_{sat} was formulated following Reed & Jablonowski (2012), (14)
where is the saturation vapour pressure at the control temperature T_{0}, ϵ is the ratio of the gas constant for dry air to that for water vapour.
If the relative humidity exceeds the saturation vapour pressure, latent heat is released and the condensate is immediately removed. This simple scheme allows us to calculate the largescale precipitation rate (P_{ls}, Thatcher & Jablonowski 2016) from (15)
where nv is the number of layers, k the layer index, and ρ_{q} is density of the condensible gas (water in our experiments) in the atmosphere.
Main parameters for the 3D simulations.
3.1.2 Boundary conditions
In both simulations explored, we prescribed the surface temperature of the ocean covered planets. The formulation for the surface temperaturein the ocean Earthlike simulation is given in Reed & Jablonowski (2012), (16)
where T_{min} is minimum temperature at the poles, Δ T is the temperature difference between equator and poles, θ is the latitude and Δ l is the width of the Gaussian function in Eq. (16).
For the tidally locked planet experiment, we modified Eq. (16) to represent the daynight contrast, (17)
where ϕ is the longitude. The parameters for the tidally locked planet were calibrated to the results of Yang et al. (2019).
3.1.3 Surface fluxes
To represent the mechanical interaction between the atmosphere and the surface we implemented two different schemes for the ocean Earthlike simulation and the tidally locked planet. For the ocean Earthlike simulation, we implemented the Rayleigh linear friction scheme exactly as done in the HeldSuarez experiment (Held & Suarez 1994), (18) (19)
where v_{h} is the 3D wind field projected parallel to the spherical surface, K_{v} is the strength of the frictional damping, σ_{b} sets the depth of the winds damped, and k_{f} is the strength of the surface dissipation. The values for the parameters used in this scheme are calibrated to do Earthlike simulations. Thatcher & Jablonowski (2016) used the same scheme and constant values. In the tidally locked planet, we did one further step in terms of complexity to allow the code to be more flexible but still trying to keep it simple in terms of implementation. We followed the implementation suggested in Reed & Jablonowski (2012). The time rate of each horizontal velocity component in the lowermost layer is defined as (20) (21)
where u_{a} and v_{a} are the horizontal components of the wind at the lowermost model layer, C_{d} is the drag coefficient and z_{a} is the altitude of the first layer (z_{a} is set to 50 m in both simulations). We are assuming that the wind velocities at the ocean surface are zero. The drag coefficient, C_{d}, depends on the magnitude of the horizontal wind at the lowermost layer as (22)
The constants C_{d0} and C_{d1} are the same as in Smith & Vogl (2008). The vertical velocity time rate is set to zero. This simple formulation represents the eddy turbulence caused by the interaction between ocean surface and atmosphere. The formulation for the kinematic eddy flux of water vapour and heat at the surface is the same for both simulations: (23) (24)
where C_{E} and C_{H} are the bulk transfer coefficients for water and heat, respectively. T_{s} is the surface temperature, and q_{sat,s} is the saturation specific humidity calculated at surface temperature. The bulk transfer coefficients were set to 0.0044, which are four times higher than as in Smith & Vogl (2008), as suggested in Thatcher & Jablonowski (2016). More details about the scheme are available in Reed & Jablonowski (2012).
3.1.4 Boundary layer
In the ocean Earthlike simulation, we applied the Rayleigh friction from Held & Suarez (1994) to represent the boundarylayer mixing of the horizontal winds. This simplification was applied to be consistent with the results from Thatcher & Jablonowski (2016). For the ocean tidally locked planet simulation, we used a simple diffusive boundary scheme. This scheme is an extension of the formulation described in the section above for the surface fluxes. The turbulent kinematic mixing is represented by (25) (26)
where the prime symbols represent the deviations from timeaveraged quantities, the overbar is the time average operator, w is the vertical velocity, and K_{m} is the eddy diffusivity coefficient for momentum. The eddy coefficient K_{m} is defined as (27)
The variable p_{top} sets the top of the boundary layer, which for simplicity is set at the same number for every column. For altitudes above p_{top}, the strength of K_{m} decreases exponentially with a rate set by p_{strato}.
The vertical turbulent mixing of potential temperature and specific humidity is the same in both simulations, (28) (29)
where K_{E} is the eddy diffusivity coefficient for energy defined as (30)
The formulation of K_{E} is identical to K_{m}. More details about this simple boundary layer scheme can be found in Reed & Jablonowski (2012) and Thatcher & Jablonowski (2016).
3.1.5 Radiation
A Newtonian temperature relaxation represents the impact of the radiation in the atmosphere. The scheme is very similar for the two planets explored. The diabatic heating for the ocean Earthlike experiment is almost identical to the HeldSuarez test (Held & Suarez 1994). The temperatures in the atmosphere are forced towards a radiativeconvective equilibrium temperature (T_{eq}) at a specific timescale (k_{T}): (31)
The equilibrium temperature is defined such as in Held & Suarez (1994), (32)
The values of each parameter are defined in Table 2. The differences between values used in our simulations (suggested in Thatcher & Jablonowski 2016) and the original HeldSuarez test (Held & Suarez 1994) are the values for T_{Equator} and . In Thatcher & Jablonowski (2016), T_{Equator} and are 315 K and 65 K, respectively, instead of 294 K and 60 K from Held & Suarez (1994).
The timescale for the temperature forcing is defined as (33)
The different parameters in the equation are constants defined in Table 2. To represent the radiativeconvective forcing in the tidally locked planet, we changed both Eqs. (32) and (33). The equations and the parameters were adapted to reproduce a similar temperature structure to the case of a tidally locked planet orbiting an Mstar from Yang et al. (2019). Other works have also used Newtonian temperature relaxation to explore the climate of terrestrial tidally locked planets, such as Heng et al. (2011) and Carone et al. (2015). Our modified equations for this experiment are (34) (35)
The relaxation time coefficient is set as a constant in the tidally locked planet simulation. In Eq. (34), the term cos θ_{z} refers to thecosine of the zenith angle. We again note that the value of each parameter is explicitly shown in Table 2.
4 Results from 3D simulations
This section presents the 3D simulations from the Earthlike planet and the tidally locked planet, where both have the surface covered with an ocean. As explained in the previous section, these simulations include simplified physics schemes and aim to reproduce conditions that allow us to evaluate the mass and moist transport accuracy in 3D simulations. Models containing poor representations of moist transport in the atmosphere will produce errors that propagate, for example, to the energy budget of the atmosphere, cloud cover, and atmospheric circulation, which can compromise the interpretation of observational data on exoplanets. To assess the simulations, we characterise the global circulation in both simulations and pay particular attention to possible sources of highfrequency waves associated with bad numerical implementations. For the ocean Earthlike simulation, Thatcher & Jablonowski (2016) showed that the simplified physics applied here can reproduce the results from models that used more complex physics formulations (e.g. including microphysics). In the section below, we compare our results to the results presented in Thatcher & Jablonowski (2016). We suggest that the results obtained in this work become the first step towards a complete benchmarking of future moist atmospheric simulators on terrestrial planets.
Fig. 5 Evolution of the total mass error in the ocean Earthlike simulation. 
4.1 Ocean Earthlike planet
We integrated our model for 1200 Earth days with a time step of 100 s. The first 200 Earth days are discarded as being part of the numerical spinup phase of the simulation.
The new solvers implemented in OASIS aim to improve the mass transport representation in the 3D simulations. Two of the main improvements are the shapepreserving of the transported mass distribution and mass conservation. Figure 5 shows how the total mass error (M_{error}) evolved during the simulation. The total mass error was calculated using the following equation, (36)
where M is the spatially integrated atmospheric mass for a particular instant and M_{i} is the spatially integrated atmospheric mass at the beginning of the simulations. The small error measured along the simulation is always below 5 × 10^{−13}%, and it does not present any constant drift, which highlights the robustness of the new formulation.
4.1.1 Zonal winds and temperature
The spinup of the simulated atmosphere is mainly driven by the thermal forcing and the rotational effects. In Fig. 6 we show the longitudinal and timeaveraged temperature and zonal^{2} winds. As expected, the temperatures are warmer at low latitudes, where the lapse rate^{3} is also larger. The tropopause is located at roughly 200 mbar. The cold and inactive stratosphere is located above the 200 mbar. Hotter layers in the atmosphere extend to higher altitudes compared to the traditional HeldSuarez Test (Held & Suarez 1994), which is driven by the surplus of heat provided from the condensation in the largescale precipitation scheme. The more significant latitudinal gradients in temperature drive the formation of stronger winds, as it is possible to derive from the approximated geostrophic thermal wind equation. The average longitudinal winds in our simulation reach a maximum of around 50–55 m s^{−1} at 200 mbar and 30° latitude. The winds decrease their magnitude in general with the decrease in altitude. Our winds and temperature results are quantitatively similar to the results presented in Thatcher & Jablonowski (2016). However, the midlatitude jets in our simulation are slightly weaker (roughly 5 m s^{−1} difference). The differences are likely caused by the different levels of diffusion applied in the various models or different model implementations, such as our model diagnostics being updated in the dynamical core instead of being done immediately in the physical core (Mendonça & Buchhave 2020; Deitrick et al. 2020).
Fig. 6 Final time and longitudinalaveraged zonal winds (m s^{−1}) and temperatures (K) for the ocean Earthlike experiment. The values were timeaveraged for 1000 Earth days. 
Fig. 7 Final time and longitudinalaveraged mass stream function function (10^{10} kg s^{−1}) for the ocean Earthlike experiment. The values were timeaveraged for 1000 Earth days. 
4.1.2 Mass stream function
In Fig. 7 we show the mass stream function, which we calculated from the following equation: (37)
where R_{planet} is the planet radius, ϕ is the latitude, g is the surface gravity, p is the atmospheric pressure and is the zonal and timeaveraged meridional wind^{4}. The results show the wellknown atmospheric cells: Hadley, Ferrel, and polar cells (latitudes higher than 50°. The latitudinal extension of the atmospheric cells is smaller in this test than in the HeldSuarez test, which assumes a dry atmosphere.The characteristics of the atmospheric cells simulated in this work are consistent with the results from, e.g. Lee et al. (2008) and Thatcher & Jablonowski (2016). The different positions of the atmospheric cells also impact the latitudinal position of the midlatitudinal jets that are moved into lower latitudes. The stronger upwelling motion by mean circulation is located inthe equatorial region at pressure levels around 800 mbar, where largescale precipitation is at its maximum.
4.1.3 Meridional eddy sensible heat flux and eddy kinetic energy
Wave activity has an important role in driving atmospheric circulation in the ocean Earthlike planet. In Fig. 8 we show the meridional eddy sensible heat flux and the eddy kinetic energy for the ocean Earthlike simulation. These two quantities are calculated from the eddy components^{5} of the temperature and wind field. In general, atmospheric waves transport sensible heat from the tropics to higher latitudes. The meridional eddy sensible heat flux () map shows a dominant poleward eddy heat transport below the pressure level 400 mbar. The poleward transport of eddy sensible heat flux becomes again more intense above 200 mbar. The maximum of meridional eddy sensible heat flux is located between 30° to 40° latitude at roughly 850 mbar. The eddy kinetic energy (0.5) shown in Fig. 8, reaches the maximum at the midlatitude jet’s location, where the largest dynamical waveactivity is present. Both meridional heat flux and eddy kinetic energy obtained here are very similar to the results from the dry atmosphere HeldSuarez experiment (Held & Suarez 1994; Mendonça et al. 2016).
4.1.4 Vertical velocity
Figure 9 shows the vertical velocity in the tropical region. The vertical velocity results are consistent with the massstream function shown above. The negative values near the equator represent the upward branch of the Hadley cell, which reverses its direction for latitudes above roughly 6°. The strength of the Hadley circulation is enforced by the condensation of the water vapour (release of latent heat), which is more abundant in the equatorial region.
In Fig. 10 we show a snapshot of the vertical velocity at a pressure level close to the surface. The predominant negative value in the equatorial region is related to the previously discussed upward branch of the Hadley cell. It is also clear from the patternsin the vertical velocity map that the wave activity propagates from the tropical region towards higher latitudes. The vertical winds near the surface are also important to check for any numerical noise associated with the underlying grid (grid imprinting). A small perturbation in the vertical wind field would be easily detected in a simulation due to its low magnitude compared toother diagnostic variables (Ullrich & Jablonowski 2012). In our model, we use the icosahedral grid. If the numerical noise becomes significant, we would see the geometrical structure of the grid, which is not the case in our simulations. In regions of a large gradient of the model diagnostics, our dispersive finitevolume methods, such as the one we apply to the momentum and entropy, can produce numerical noise, which would be easily seen as high gravity waves. However, our solver is robust and did not show any sign of numerical problems associated with the nature of the central finite volume.
Fig. 8 Final time and longitudinalaveraged meridional eddy sensible heat flux (m s^{−1} K) and eddy kinetic energy (m^{2} s^{−2} K) for the ocean Earthlike experiment. The values were timeaveraged for 1000 Earth days. 
Fig. 9 Final time and longitudinalaveraged vertical velocity (Pa s^{−1}) between 15° S and 15° N for the ocean Earthlike simulation. The values were timeaveraged for 1000 Earth days. 
Fig. 10 Instantaneous latitudelongitude map of vertical velocity (Pa s^{−1}) for the ocean Earthlike experiment. 
4.1.5 Kinetic energy spectrum
Kinetic energy spectra provide important information about how atmospheric motion is redistributed across different spatial scales. In Fig. 11 we show the rotational and divergence components of the kinetic energy as a function of spherical wavenumber (e.g. Augier & Lindborg 2013). Our results show that at large scales (low wavenumbers), the rotational component is on average orders of magnitude larger than the divergent component. The magnitude of the divergent component increases on average for smaller horizontal scales, where for example, diabatic heating can increase baroclinicity and intensify thermally direct circulation. The two components become comparable at wavenumbers close to the truncation scales, where the energy in the divergent component becomes similar to the energy of the rotational component. In Fig. 11, we include a reference line of k^{−3}, where k represents the spherical wavenumbers. The slope results from the downscale cascade of the enstrophy, which has been suggested by theoretical work and observations (e.g. Nastrom & Gage 1985). At higher resolution simulations, it would be possible to capture the transition to as a result of thedownscale cascade of energy (Lindborg & Cho 2000). Despite the simple physics represented in our simulation, our results are broadly consistent with previous studies based on observational or reanalysis datasets of Earth’s atmosphere (e.g. Koshyk et al. 1999). The model needs to remove enstrophy or energy down at the spectral truncation scale to reproduce physically consistent results. Otherwise, physical quantities would build up at the truncation resolution. It is common practice for largescale atmospheric models to remove kinetic energy near the truncation limit via, for example, explicit hyperdiffusion methods or numerical implicit schemes. Our spectrum in Fig. 11 shows a steeper slope than k^{−3} at higher wavenumbers as expected by the use of the fourthorder hyperdiffusion and divergence damping (Mendonça et al. 2016). In general, simulations should overestimate the dissipation near the truncation resolution to maintain good stability performance and avoid nonphysical artefacts arising from the sudden cut in resolution aliasing onto the wellresolved modes (smaller wavenumbers). Also, the physics at the truncation scales are not accurately solved by the atmospheric models. For an accurate representation of the largescale physical phenomena, we should mitigate the inaccuracies by removing the energy at these scales.
Fig. 11 Timeaveraged kinetic energy spectrum at 0.25 mbar for the ocean Earthlike simulation. The dashed and dotted lines are therotational and divergent components of the kinetic energy, respectively, and the solid line is the total (sum of the rotational and divergent components). 
4.1.6 Water vapour concentration, relative humidity, and precipitation
In Fig. 12a we show how the time and longitudinal average water vapour concentration is distributed across latitudes and pressures. The water vapour is evaporated from the ocean surface, predominantly from the equatorial region. The larger concentration of water vapour remains at low latitude, which is also the warmer atmospheric region, as we showed in Fig. 6. The elongated vertical feature in the equatorial region is driven by the accelerated vertical winds caused by the latent heat released by water condensation. Water vapour saturation is higher in the lower atmosphere,as is shown in Fig. 12b. The low saturation levels at low latitudes are caused by largescale precipitation.
In Fig. 13 we plot the simulated surface precipitation. As expected from the distribution of the water vapour concentration, the higher levels of precipitation are located in the equatorial region. Despite the complex distribution of rainfall across the surface, the pattern obtained (climate state) is persistent. At midlatitudes, instabilities grow by feeding on the available potential energy associated with the meridional temperature gradient (baroclinic instabilities). These instabilities transport higher levels of precipitation to higher latitudes. However, the region between 15 and 20° latitude, shows very low levels of precipitation. This complex system of precipitation is also obtained in simulations by complex Earth climate models (e.g. including cloud microphysics) on ocean Earthlike simulations (e.g. Thatcher & Jablonowski 2016).
The absence of any sign of numerical error that could compromise the physics of the simulations and the quantitatively very similar results between our simulation and Thatcher & Jablonowski (2016) demonstrates that our new OASIS has passed a critical ocean Earthlike benchmark test.
4.2 Ocean tidally locked planet
In this section we analyse the results from an ocean covered tidally locked planet experiment. The model parameters are estimated to reproduce quantitatively the atmospheric (plus ocean) temperature structure of the tidally locked planet orbiting an M star presented in Yang et al. (2019). The planet parameters for our simulations are shown in Table 2. The simulations in Yang et al. (2019) are used as a reference, and we did not try to reproduce the same results since we use simpler physical numerical schemes in our simulations (e.g. we do not include radiative feedbacks from clouds). Our primary focus is to build an easy setup simulation to be used in the future as a benchmark test for moist mass transport in tidally locked planets. Note that in Yang et al. (2019), different models with complex physics implemented obtained quantitatively different results on this particular planet setup.
Our simulation has been integrated for 1200 Earth days, similar to the ocean Earthlike simulation. The time step is 50 s, half of the value used for the ocean Earthlike simulation, due to numerical stability issues associated with the different atmosphericdynamics. Again, in this experiment, we discarded the first 200 Earth days as part of the numerical spinup phase of the simulation.
Figure 14 shows the total mass error for the ocean tidally locked planet simulation. As for the aqua Earthlike case, the error measured along the simulation is small and not larger than 7.5×10^{−13}%, which is again an example of the good properties of the newly implemented solver conserving atmospheric mass in the 3D simulations.
Fig. 12 Final time and longitudinalaveraged water vapour concentration (kg/kg) and relative humidity (%) for the ocean Earthlike simulation. The values were timeaveraged for 1000 Earth days. 
Fig. 13 Instantaneous latitudelongitude map of surface precipitation (mm h^{−1}) for the ocean Earthlike simulation. 
4.2.1 Zonal winds and temperature
The tidally locked planet is rotating five times slower than Earth, which strongly impacts atmospheric circulation. In Fig. 15a we present the averaged zonal winds. The winds obtained are, in general, slower than the winds obtained in the previous experiment. However, as expected from slowly rotating planets, the equatorial winds become stronger (e.g. Merlis & Schneider 2010). The weaker Coriolis force allows more efficient transport of heat from low to high latitudes. The more efficient poleward transport of heat reduces the latitudinal temperature gradients as shown in Fig. 15b compared to the ocean Earthlike planet case. Another significant difference between the previous experiment and the tidally locked planet case is the vertical gradient of the temperature. In the tidally locked experiment, the atmospheric temperatures reach a maximum at pressures of nearly 800 mbar due to the release of latent heat from water condensation. The latent heat released enhances the upward motion on the dayside of the planet. We note that the substellar point is at 180° longitude. The effect of the latent heat released produces a peaky feature in the temperature structure, as is seen in Fig. 16. The heat produced on the dayside is transported to the nightside by the prevailing eastward wind, which is driven by the rotation of the planet. The longitudinal heat transport is generally efficient enough to maintain the daynight differences very small, except for pressures higher than roughly 700 mbar. In the lower part of the atmosphere, the simulations showed an asymmetry in the temperature structure between the west and east side of the substellar point. The colder regions of the atmosphere near the surface are due to the very cold temperatures prescribed to reproduce the ocean conditions simulated inYang et al. (2019), which reach temperatures below the water freezing point in some regions of the planet.
Fig. 14 Evolution of the total mass error in the ocean tidally locked planet simulation. 
4.2.2 Mass stream function
The averaged mass stream function results are shown in Fig. 17. Two largescale circulation cells form in each hemisphere that extends from the equator to the poles. These cells make the heat transport from low latitudes to the poles more efficient than in the ocean Earthlike experiment. The planet size cells are formed due to the slow rotation of the planet that weakens the impact of the Coriolis force driving the circulation. Largescale cells such as the ones shown in Fig. 17 are expected in Venus (e.g. Mendonça & Read 2016) and Titan (Lebonnois et al. 2012). However, in Titan, the size and position of the largescale cells evolve in time due to the seasonality effects.
Fig. 15 Final time and longitudinalaveraged zonal winds (m s^{−1}) and temperatures (K) for the ocean tidally locked planet experiment. The values were timeaveraged for 1000 Earth days. 
Fig. 16 Final time and longitudinalaveraged temperature (K) across the equator for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 
4.2.3 Vertical velocity
The vertical component of the winds is shown in Fig. 18. The distribution of the vertical winds in the equatorial region is similar to the structure obtained in the ocean Earthlike planet experiment. However, the magnitude of the averaged upward winds is stronger in the tidally locked case. The vertical velocities are enhanced in the equatorial region due to the more intense latent heat released from water condensation. The stronger vertical values are located just above the warm region produced by the water condensation and centred at the pressure level 700 mbar.
Fig. 17 Final time and longitudinalaveraged mass stream function (10^{10} kg s^{−1}) for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 
4.2.4 Kinetic energy spectrum
We also computed the kinetic energy spectrum for the tidally locked planet simulation. The kinetic spectrum is shown in Fig. 19. Compared with the ocean Earthlike planet case, the presence of the stronger components at shorter wavenumbers is clear. As in the Earthlike case, the divergent component of the kinetic energy is much smaller than the rotational component at the smallest wave numbers. However, the two components become comparable at around the spherical wavenumber 10. The characteristics of the temperature forcing (tidally locked configuration) plus the slower rotation of the planet compared to the Earthlike experiment resulted in a stronger generation of divergent kinetic energy. The slope near the truncation limit becomes larger due to the hyperdiffusion and divergent damping processes applied. A buildup of energy indicates that models produce numerical noise that is not being correctly removed. This is not the case in our model, as inthe ocean Earthlike simulation, we also do not see a buildup of energy at the truncation limit in this experiment. The results show that our new methods are robust enough to simulate mass transport in moist atmospheres.
Fig. 18 Final time and longitudinalaveraged vertical velocity (Pa s^{−1}) between 15° S and 15° N for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 
Fig. 19 Timeaveraged kinetic energy spectrum at 0.25 mbar for the ocean tidally locked planet simulation. The dashed and dotted lines are the rotational and divergent components of the kinetic energy, respectively, and the solid line is the total (sum of the rotational and divergent components). 
4.2.5 Water vapour concentration, relative humidity, and precipitation
The averaged water concentration and relative humidity are shown in Fig. 20. The largest concentration of water vapour is located at low latitudes in the lower atmosphere. As expected from the temperature distribution, we also obtain an elongated vertical distribution at the equator due to the enhanced acceleration of the vertical winds by latent heat released from water condensation. Above the pressure level of 400 mbar, the concentrations of water vapour become very small. Figure 20 shows low levels of relative humidity in the lower atmosphere, which is caused by precipitation. Higher relative humidity values are located near the surface and in a pressure region between 50 and 200 mbar. The relatively high values of relative humidity at low pressures seem to indicate that the atmosphere is not efficient at removing water vapour at these altitudes since this region is not cold enough to condense water and vertical mixing is not efficient.
In Fig. 21 we show the surface precipitation. Our results show a permanent precipitation region on the dayside of the planet. The complete water cycle is only possible on the dayside, where a large amount of water is evaporated from the hot dayside surface regions, transported upwards and precipitates. This cycle does not extend to the nightside. This configuration is a characteristic of the parameter space explored, and it does not mean that permanent precipitation on terrestrial tidally locked planets will only happen on the dayside. A more thorough exploration of the parameter space (e.g. rotation rate and surface pressure) needs to be done in the future. The precipitation map does not register any possible indication of numerical error, which would be noticed by a random distribution of maxima in the precipitation. The precipitation structure forms a round shape structure with spiral arms, which resembles of tropical cyclones on Earth. However, the formation of the largescale precipitation structure onthe dayside of the tidally locked planet is different from Earth tropical cyclones, which are mainly driven by rotational effects. The tidally locked planet largescale precipitation region is formed mostly due to divergent processes triggered by the release of latent heat from the water condensation and the permanent dayside configuration.
5 Conclusions
The formulation presented in this work is part of the ambitious development of the PCM OASIS (Mendonça & Buchhave 2020) that has been built from scratch to explore a large diversity of planets. OASIS includes various physical and chemical schemes in a modular approach that allows planetary climate problems to be explored at different levels of complexity. Every module implemented in the OASIS platform is carefully assessed in terms of performance, accuracy, physical consistency, and flexibility to explore an extensive range of planetary conditions. A complete design of a condensible cycle in 3D models needs to include methods to represent, for example, mass transport, cloud formation, cloud radiative feedbacks, and moist atmospheric thermodynamics. In this work, we focus on the implementation of the first step towards the complete implementation of a robust condensible cycle in 3D simulations: the mass transport scheme.
We have implemented an upwindbiased scheme on a piece linear approximation with a flux limiter from Miura (2007) to represent the mass transport in the 3D simulations. The new numerical method has been successfully tested in problems with different complexity. We first explored the new solver on a 2D model where we transport a bellshaped cosine distribution across the planet. Compared to the old central finite volume scheme, the results of the new solver show a significantimprovement in terms of, for example, mass conservation, accuracy, shapepreserving properties, and avoiding negative solutions. The new solver has a secondorder accuracy and avoids using any fixer, such as hyperdiffusion applied to the mass variable, which improves the physical consistency in the model. We compared the 2D results for different horizontal resolutions and showed that for 1degreeresolution simulations the error is less than 1% for one complete revolution of the peak. The error needs to be kept as small as possible to avoid any unphysical feedback processes in the simulation. We recommend using at least 1degree resolution for moist atmospheric simulations to avoid impactful errors and a low degreeof implicit numerical diffusion.
Another goal of this work was to test our new formulation on 3D simulations that include a simple physical representation of the moist processes in the atmosphere. To accomplish this goal, we benchmarked our model against a test developed in the Earth climate community and propose a new tidally locked terrestrial planet benchmark test.
Our new model was capable of passing the ocean Earthlike benchmark test successfully. The results obtained with our new OASIS model are quantitatively similar to the results presented in Thatcher & Jablonowski (2016), who designed the test. The simulations produce robust representations of atmospheric mass transport and consistent water vapour distributions with more complex Earth climate models that include sophisticated schemes to represent cloud formation. It is well known that moist physics tendencies in 3D numerical simulation models can trigger undesirable largescale gravity waves (e.g. Thatcher & Jablonowski 2016), which were not detected in any of our simulations. Our results also showed that our new formulations do excellent work conserving the total atmospheric mass (error is less than 5 × 10^{−13}%).
We offer the exoplanet community a robust new benchmark test for 3D exoplanet models to simulate ocean tidally locked planets. The new test explores the challenges of simulating a moist atmosphere of a tidally locked planet. The bulk parameters of our tidally locked planet are based on the simulations studied in Yang et al. (2019), who explore a tidally locked planet orbiting an M star. Our results show a planet where water evaporates from the permanent dayside ocean and precipitates on the same side. Condensation occurs in the lower atmosphere and has a stronger impact on the temperature and wind field than in the ocean Earthlike case. In addition, the slower rotation of the planet drives the atmospheric circulation to produce stronger winds at the equatorial region, which creates permanent strong easterly winds that have a substantial impact on the transport of water vapour from the dayside to the nightside. Our simulations showed a good performance of the new solver in terms of accuracy and does not produce any signs of undesirable numerical highfrequency waves (e.g. numerical noise). Furthermore, the benchmark is easy to set up. We encourage other groups to carry out our proposed test as a first step in evaluating the methods used to explore the exotic planet environments present in ocean tidally locked terrestrial planets.
This work is the first part of a series of manuscripts that will focus on the full implementation of moist processes in 3D planetary models. We will present the complete formalism implemented in the OASIS model and propose benchmark tests with different levels of complexity to evaluate all the various steps of the numerical implementation.
Fig. 20 Final time and longitudinalaveraged water vapour concentration (kg/kg) and relative humidity (%) for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 
Fig. 21 Instantaneous latitudelongitude map of surface precipitation (mm h^{−1}) for the ocean tidally locked planet simulation. 
Acknowledgements
J.M.M. acknowledges financial support from the PRODEX project number 4000127377. The figures in this study have been done with the NCAR Command Language (Version 6.6.2) [Software. (2019). Boulder, Colorado: UCAR/NCAR/CISL/TDD, http://dx.doi.org/10.5065/D6WD3XH5.
References
 Augier, P., & Lindborg, E. 2013, J. Atm. Sci., 70, 2293 [NASA ADS] [CrossRef] [Google Scholar]
 Carone, L., Keppens, R., & Decin, L. 2015, MNRAS, 453, 2412 [Google Scholar]
 Deitrick, R., Mendonça, J. M., Schroffenegger, U., et al. 2020, ApJS, 248, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Durran, D. R. 1999, Numerical Methods for Wave Equations in Geophysical Fluid Dynamics SpringerVerlag, 465 [Google Scholar]
 Fauchez, T. J., Turbet, M., Villanueva, G. L., et al. 2019, ApJ, 887, 194 [Google Scholar]
 Ge, H., Li, C., Zhang, X., & Lee, D. 2020, ApJ, 898, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Held, I. M., & Suarez, M. J. 1994, Bull. Am. Meteorol. Soc., 75, 1825 [CrossRef] [Google Scholar]
 Heng, K., Frierson, D. M. W., & Phillipps, P. J. 2011, MNRAS, 418, 2669 [NASA ADS] [CrossRef] [Google Scholar]
 Hörst, S. M. 2017, J. Geophys. Res. Planets, 122, 432 [CrossRef] [Google Scholar]
 Ingersoll, A. P. 1969, J. Atm. Sci., 26, 1191 [NASA ADS] [CrossRef] [Google Scholar]
 Komacek, T. D., Fauchez, T. J., Wolf, E. T., & Abbot, D. S. 2020, ApJ, 888, L20 [Google Scholar]
 Koshyk, J. N., Boville, B. A., Hamilton, K., Manzini, E., & Shibata, K. 1999, J. Geophys. Res., 104, 27177 [NASA ADS] [CrossRef] [Google Scholar]
 Lebonnois, S., Hourdin, F., Eymet, V., et al. 2010, J. Geophys. Res. Planets, 115, E6 [Google Scholar]
 Lebonnois, S., Burgalat, J., Rannou, P., & Charnay, B. 2012, Icarus, 218, 707 [CrossRef] [Google Scholar]
 Lee, M.I., Suarez, M. J., Kang, I.S., Held, I. M., & Kim, D. 2008, J. Clim., 21, 4934 [NASA ADS] [CrossRef] [Google Scholar]
 Lindborg, E., & Cho, J. Y. N. 2000, Phys. Rev. Lett., 85, 5663 [NASA ADS] [CrossRef] [Google Scholar]
 Mendonça, J. M., & Buchhave, L. A. 2020, MNRAS, 496, 3512 [CrossRef] [Google Scholar]
 Mendonça, J. M., & Read, P. L. 2016, Planet. Space Sci., 134, 1 [CrossRef] [Google Scholar]
 Mendonça, J. M., Grimm, S. L., Grosheintz, L., & Heng, K. 2016, ApJ, 829, 115 [CrossRef] [Google Scholar]
 Mendonça, J. M., Tsai, S.m., Malik, M., Grimm, S. L., & Heng, K. 2018, ApJ, 869, 107 [CrossRef] [Google Scholar]
 Merlis, T. M.,& Schneider, T. 2010, J. Adv. Model. Earth Syst., 2, 13 [Google Scholar]
 Miura, H. 2007, Mon. Weather Rev., 135, 4038 [NASA ADS] [CrossRef] [Google Scholar]
 Nastrom, G. D., & Gage, K. S. 1985, J. Atm. Sci., 42, 950 [NASA ADS] [CrossRef] [Google Scholar]
 Read, P. L. 1986, Q. J. R. Meteorol. Soc., 112, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Read, P. L., Barstow, J., Charnay, B., et al. 2016, Q. J. R. Meteorol. Soc., 142, 703 [NASA ADS] [CrossRef] [Google Scholar]
 Reed, K. A., & Jablonowski, C. 2012, J. Adv. Model. Earth Syst., 4, M04001 [Google Scholar]
 SánchezLavega, A., Lebonnois, S., Imamura, T., Read, P., & Luz, D. 2017, Space Sci. Rev., 212, 1541 [CrossRef] [Google Scholar]
 Satoh, M. 2002, Mon. Weather Rev., 130, 1227 [NASA ADS] [CrossRef] [Google Scholar]
 Skamarock, W. C., & Klemp, J. B. 2008, J. Comput. Phys., 227, 3465 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. K., & Vogl, S. 2008, Q. J. R. Meteorol. Soc., 134, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Suissa, G., Mandell, A. M., Wolf, E. T., et al. 2020, ApJ, 891, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, F. W. 2011, Planet. Space Sci., 59, 889 [NASA ADS] [CrossRef] [Google Scholar]
 Thatcher, D. R., & Jablonowski, C. 2016, Geosci. Model Dev., 9, 1263 [NASA ADS] [CrossRef] [Google Scholar]
 Thuburn, J. 1996, J. Comput. Phys., 123, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Tomita, H., & Satoh, M. 2004, Fluid Dyn. Res., 34, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Tomita, H., Tsugawa, M., Satoh, M., & Goto, K. 2001, J. Comput. Phys., 174, 579 [NASA ADS] [CrossRef] [Google Scholar]
 Ullrich, P. A., & Jablonowski, C. 2012, J. Comput. Phys., 231, 5078 [NASA ADS] [CrossRef] [Google Scholar]
 Ullrich, P. A., Jablonowski, C., Kent, J., et al. 2017, Geosci. Model Dev., 10, 4477 [NASA ADS] [CrossRef] [Google Scholar]
 Wicker, L. J., & Skamarock, W. C. 2002, Mon. Weather Rev., 130, 2088 [NASA ADS] [CrossRef] [Google Scholar]
 Williamson, D. L., Drake, J. B., Hack, J. J., Jakob, R., & Swarztrauber, P. N. 1992, J. Comput. Phys., 102, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Wordsworth, R. 2015, ApJ, 806, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, J., Leconte, J., Wolf, E. T., et al. 2019, ApJ, 875, 46 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Time step and spatial resolutions for the cosine bell test cases using different transport schemes.
All Figures
Fig. 1 Initial cosine bell shape defined in Eq. (9). The colour map is in arbitrary units. 

In the text 
Fig. 2 Latitudelongitude maps of the difference between the initial cosine bell and after a complete revolution across the sphere (spatial resolution glevel 6). The three maps refer to three different methods that used different time steps (Table 1): CFV (central finitevolume), UW without FL (upwind biased without flux limiter), and UW with FL (upwind with flux limiter). 

In the text 
Fig. 3 Final tracer distributions at the equator from the cosine bell test case that uses the upwindbiased scheme with a flux limiter for different spatial resolutions. The solid black line is the analytical solution used for the initial condition (Eq. (9)). 

In the text 
Fig. 4 l_{∞} convergence plot where glevel represents the different spatial resolutions. The various solid lines represent different methods: the central finite volume (red); the upwindbiased with flux limiter (green); and upwind biased without flux limiter (blue). The solid black line is a reference line with a slope of an exact secondorder convergence. 

In the text 
Fig. 5 Evolution of the total mass error in the ocean Earthlike simulation. 

In the text 
Fig. 6 Final time and longitudinalaveraged zonal winds (m s^{−1}) and temperatures (K) for the ocean Earthlike experiment. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 7 Final time and longitudinalaveraged mass stream function function (10^{10} kg s^{−1}) for the ocean Earthlike experiment. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 8 Final time and longitudinalaveraged meridional eddy sensible heat flux (m s^{−1} K) and eddy kinetic energy (m^{2} s^{−2} K) for the ocean Earthlike experiment. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 9 Final time and longitudinalaveraged vertical velocity (Pa s^{−1}) between 15° S and 15° N for the ocean Earthlike simulation. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 10 Instantaneous latitudelongitude map of vertical velocity (Pa s^{−1}) for the ocean Earthlike experiment. 

In the text 
Fig. 11 Timeaveraged kinetic energy spectrum at 0.25 mbar for the ocean Earthlike simulation. The dashed and dotted lines are therotational and divergent components of the kinetic energy, respectively, and the solid line is the total (sum of the rotational and divergent components). 

In the text 
Fig. 12 Final time and longitudinalaveraged water vapour concentration (kg/kg) and relative humidity (%) for the ocean Earthlike simulation. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 13 Instantaneous latitudelongitude map of surface precipitation (mm h^{−1}) for the ocean Earthlike simulation. 

In the text 
Fig. 14 Evolution of the total mass error in the ocean tidally locked planet simulation. 

In the text 
Fig. 15 Final time and longitudinalaveraged zonal winds (m s^{−1}) and temperatures (K) for the ocean tidally locked planet experiment. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 16 Final time and longitudinalaveraged temperature (K) across the equator for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 17 Final time and longitudinalaveraged mass stream function (10^{10} kg s^{−1}) for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 18 Final time and longitudinalaveraged vertical velocity (Pa s^{−1}) between 15° S and 15° N for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 19 Timeaveraged kinetic energy spectrum at 0.25 mbar for the ocean tidally locked planet simulation. The dashed and dotted lines are the rotational and divergent components of the kinetic energy, respectively, and the solid line is the total (sum of the rotational and divergent components). 

In the text 
Fig. 20 Final time and longitudinalaveraged water vapour concentration (kg/kg) and relative humidity (%) for the ocean tidally locked planet simulation. The values were timeaveraged for 1000 Earth days. 

In the text 
Fig. 21 Instantaneous latitudelongitude map of surface precipitation (mm h^{−1}) for the ocean tidally locked planet simulation. 

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