© 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 high-frequency 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 extra-solar 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 inter-comparison 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 higher-order 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 virtual-lab to study planetary environments that has been written from ground-up (Mendonça & Buchhave 2020). This platform comprises different physics and chemistry modules that interact with each other to reproduce self-consistent results. OASIS makes use of the flexible dynamical core, THOR, which has a 3D non-hydrostatic 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 ‘super-rotation’ (e.g. Read 1986; Sánchez-Lavega et al. 2017). The strong winds have an enormous impact on the planet’s climate. Using our state-of-the-art 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: $$\begin{array}{l}\frac{\partial }{\partial t}(\rho )+\nabla \cdot (\rho v)=0,\\ \frac{\partial \rho v}{\partial }+\nabla \cdot (\rho v\otimes v)=-\nabla P-\rho g\widehat{r}-2\rho \$\Omega \$\times v,\\ \frac{\partial }{\partial t}(\rho \theta )+\nabla \cdot (\rho \theta v)=\frac{Q\rho \theta }{T{c}_p},\\ \frac{\partial }{\partial t}(\rho q_i)+\nabla \cdot (\rho q_{i}v)=S_{q_i},\end{array}$$

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 $S_q^i$ 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 follow-up 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 A-grid and the vertical levels on a non-staggering grid. The spatial integration of the momentum and entropy equations are solved using the central finite-volumes 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 upwind-biased 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 (HE-VI 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¹ 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 finite-volume scheme from Tomita et al. (2001). This scheme is well known to be dominated by numerical dispersion that creates high-frequency patterns in the solutions (e.g. Durran 1999). The non-physical solutions could create maxima and minima artificially and hence produce non-physical 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 hyper-diffuse term in the tracer transport equation. The diffuse term removes the high-frequency oscillations created artificially by the central finite volume. The hyper-diffuse term in the tracer transport helps to maintain the simulation stable, but the formulation is not physical-based for mass transport. The new flux divergence operator implemented here was developed in Miura (2007). The new method is an upwind-biased advection scheme on a piece-wise 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 finite-volume scheme) and shape-preserving mass distribution (from the flux-limiter method). The flux-limiter method removes, forexample, over- or under-shoot solutions, and is important to ensure the monotonicity of the gradient concentration.

To solve Eq. (4), we first apply the finite-volume discretisation from Miura (2007): $${\rho }_0^{t+\Delta t}{q}_0^{t+\Delta t}={\rho }_0^t{q}_0^t-\frac{\Delta t}{A_0}{\sum }_i^{N_s}(l_i{\tilde{\rho }}_{R_i}{\tilde{q}}_{R_i}{v}_{R_i}^{t+\Delta t/2}\cdot n_i),$$(5)

where ρ is the atmospheric density, q the tracer concentration, A₀ 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 non-divergent wind that transports the tracer in our test is set to a solid body rotation, $$\begin{array}{ll}u(\varphi ,\lambda )\hfill & =u_0\cos (\lambda )\hfill \\ v(\varphi ,\lambda )\hfill & =0,\hfill \end{array}$$

where u and v are the zonal and meridional wind directions, respectively, ϕ is the longitude and λ is the latitude. The variable u₀ is defined as: $$u_0=\frac{2\pi R_{\text{planet}}}{P_{\text{bell}}}.$$(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⁶ m and 12 days period. The cosine bell function was constructed following Williamson et al. (1992) work: $$h(r)=\{\begin{array}{ll}(h_0/2)[1+\cos \pi r/R]\hfill & \text{if}r<R\hfill \\ 0\hfill & \text{if}r\ge R,\hfill \end{array}$$(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₀ = 1000 m and R = R_planet∕3. Other configurations of the non-divergent 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 time-steps 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 finite-volume method is very clear, despite using a much shorter time step than the one used by the other methods. The central finite-volume 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 time-consuming 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 upwind-biased 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 flux-limiter scheme.

The plot in Fig. 3 shows the final equatorial tracer distribution for the simulations using the upwind-biased 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 $$l_{\infty }=\text{max}|h-h_{\text{exact}}|,$$(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 upwind-biased 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 second-order convergence. In some regions, the solver without using the flux limiter is even better than second-order 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. 1 Initial cosine bell shape defined in Eq. (9). The colour map is in arbitrary units.

Fig. 2 Latitude-longitude 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 finite-volume), UW without FL (upwind biased without flux limiter), and UW with FL (upwind with flux limiter).

Fig. 3 Final tracer distributions at the equator from the cosine bell test case that uses the upwind-biased 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 upwind-biased 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 second-order 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 second-order 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 upwind-biased 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 hyper-diffusion 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 Earth-like test case proposed in Thatcher & Jablonowski (2016). The test is a modified version of the well-known Held-Suarez test (Held & Suarez 1994). The new test includes simplified moist processes, boundary-layer 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 large-scale 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 large-scale gravity waves from being triggered under different planetary conditions (tidally and non-tidally locked configurations).

Using a similar approach, we adapted the routines used to study the ocean Earth-like planet to explore planets under a tidally locked configuration. The two test cases (ocean Earth-like 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 Large-scale 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 re-evaporation processes. The temperature and specific humidity trends during condensation are defined as $$\frac{\partial T}{\partial t}=\frac{L}{c_{\text{p}}}C$$(11) $$\frac{\partial q}{\partial t}=-C,$$(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 $$C=\frac{\text{d}{q}_{\text{sat}}}{\text{d}t}=\frac{1}{\Delta t}\left(\frac{q-q_{\text{sat}}(T)}{1+\frac{L}{c_p}\frac{L{q}_{\text{sat}}}{R_q{T}^2}}\right),$$(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 saturation-specific humidity. The q_sat was formulated following Reed & Jablonowski (2012), $$q_{\text{sat}}(T)=\frac{ϵ}{p}{e}_0^{\star }\exp \left[-\frac{L}{R_q}(\frac{1}{T}-\frac{1}{T_0})\right],$$(14)

where $e_0^{\star }$ is the saturation vapour pressure at the control temperature T₀, ϵ 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 large-scale precipitation rate (P_ls, Thatcher & Jablonowski 2016) from $$P_{\text{ls}}\approx \frac{1}{{\rho }_{q}g}{\sum }_{k=1}^{nv}{C}_k(p_{k+1/2}-p_{k-1/2}),$$(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.

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 Earth-like simulation is given in Reed & Jablonowski (2012), $$T_s=\Delta T\exp (-\frac{{\theta }^2}{2{(\Delta l)}^2})+T_{\text{min}},$$(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 day-night contrast, $$T_s=\Delta T\exp (-\frac{{\theta }^2}{2{(\Delta l)}^2})\exp (-\frac{{\varphi }^2}{2{(\Delta l)}^2})+T_{\text{min}},$$(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 Earth-like simulation and the tidally locked planet. For the ocean Earth-like simulation, we implemented the Rayleigh linear friction scheme exactly as done in the Held-Suarez experiment (Held & Suarez 1994), $$\frac{\partial v_h}{\partial t}=K_v(\sigma )v_h,$$(18) $$K_v(\sigma )=k_f\times \max (0.0,\frac{\sigma -{\sigma }_b}{1-{\sigma }_b}),$$(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 Earth-like 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 $$\frac{\partial u_a}{\partial t}=-\frac{C_d(\sqrt{u_a^2+v_a^2})u_a}{z_a},$$(20) $$\frac{\partial v_a}{\partial t}=-\frac{C_d(\sqrt{u_a^2+v_a^2})v_a}{z_a},$$(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 $$C_d=\{\begin{array}{ll}{C}_{d0}+C_{d1}(\sqrt{u_a^2+v_a^2})\hfill & \text{if}\sqrt{u_a^2+v_a^2}<20m{s}^{-1}\hfill \\ 0.002\hfill & \text{if}\sqrt{u_a^2+v_a^2}\ge 20m{s}^{-1}.\hfill \end{array}$$(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: $$\frac{\partial q_a}{\partial t}=-\frac{C_H(\sqrt{u_a^2+v_a^2})(T_s-T_a)}{z_a},$$(23) $$\frac{\partial T_a}{\partial t}=-\frac{C_E(\sqrt{u_a^2+v_a^2})(q_{\text{sat},\text{s}}-q_a)}{z_a},$$(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 Earth-like simulation, we applied the Rayleigh friction from Held & Suarez (1994) to represent the boundary-layer 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 $$\frac{\partial u}{\partial t}=-\frac{1}{\rho }\frac{\partial \rho \overline{w^{\prime }{u}^{\prime }}}{\partial z}=\frac{1}{\rho }\frac{\partial }{\partial z}\rho K_m\frac{\partial u}{\partial z}$$(25) $$\frac{\partial v}{\partial t}=-\frac{1}{\rho }\frac{\partial \rho \overline{w^{\prime }{v}^{\prime }}}{\partial z}=\frac{1}{\rho }\frac{\partial }{\partial z}\rho K_m\frac{\partial v}{\partial z},$$(26)

where the prime symbols represent the deviations from time-averaged 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 $$K_m=\{\begin{array}{ll}{C}_d(\sqrt{u_a^2+v_a^2})z_a\hfill & \text{if}p<p_{\text{top}}\hfill \\ C_d(\sqrt{u_a^2+v_a^2})z_a\exp (-{[\frac{p_{\text{top}}-p}{p_{\text{strato}}}]}^2)\hfill & \text{if}p\le p_{\text{top}}.\hfill \end{array}$$(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, $$\frac{\partial \Theta }{\partial t}=-\frac{1}{\rho }\frac{\partial \rho \overline{w^{\prime }{\Theta }^{\prime }}}{\partial z}=\frac{1}{\rho }\frac{\partial }{\partial z}\rho K_E\frac{\partial \Theta }{\partial z}$$(28) $$\frac{\partial q}{\partial t}=-\frac{1}{\rho }\frac{\partial \rho \overline{w^{\prime }{q}^{\prime }}}{\partial z}=\frac{1}{\rho }\frac{\partial }{\partial z}\rho K_E\frac{\partial q}{\partial z},$$(29)

where K_E is the eddy diffusivity coefficient for energy defined as $$K_E=\{\begin{array}{ll}{C}_E(\sqrt{u_a^2+v_a^2})z_a\hfill & \text{if}p<p_{\text{top}}\hfill \\ C_E(\sqrt{u_a^2+v_a^2})z_a\exp (-{[\frac{p_{\text{top}}-p}{p_{\text{strato}}}]}^2)\hfill & \text{if}p\le p_{\text{top}}.\hfill \end{array}$$(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 Earth-like experiment is almost identical to the Held-Suarez test (Held & Suarez 1994). The temperatures in the atmosphere are forced towards a radiative-convective equilibrium temperature (T_eq) at a specific timescale (k_T): $$\frac{\partial T}{\partial t}=-k_T(\varphi ,\sigma )[T-T_{\text{eq}}(\varphi ,\sigma )].$$(31)

The equilibrium temperature is defined such as in Held & Suarez (1994), $$\begin{array}{ll}{T}_{\text{eq}}(\varphi ,p)=\hfill & \text{max}\{200K,[T_{\text{Equator}}-{(\Delta T)}_y{\sin }^2\varphi \text{ vphantom}{\left(\frac{p}{p_{00}}\right)}^{\kappa }\hfill \\ \hfill & -{(\Delta \theta )}_z\log \left(\frac{p}{p_{00}}\right){\cos }^2\varphi ]{\left(\frac{p}{p_{00}}\right)}^{\kappa }\}.\hfill \end{array}$$(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 Held-Suarez test (Held & Suarez 1994) are the values for T_Equator and ${(\Delta T)}_y$. In Thatcher & Jablonowski (2016), T_Equator and ${(\Delta T)}_y$ 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 $$k_T(\varphi ,\sigma )=k_a+(k_s-k_a)\text{max}(0,\frac{\frac{p}{p_{00}}-{\sigma }_b}{1-{\sigma }_b}){\cos }^4\varphi .$$(33)

The different parameters in the equation are constants defined in Table 2. To represent the radiative-convective 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 M-star 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 $$\begin{array}{ll}{T}_{\text{eq}}(\varphi ,p)=\hfill & \max \{200K,[T_{\text{Equator}}-\cos {\theta }_z({(\Delta T)}_y\hfill \\ \hfill & -{(\Delta \theta )}_z\log (\frac{p}{p_{00}}))]{(\frac{p}{p_{00}})}^{\kappa }\},\hfill \end{array}$$(34) $$k_T(\varphi ,\sigma )=k_a.$$(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 Earth-like 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 high-frequency waves associated with bad numerical implementations. For the ocean Earth-like 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 Earth-like simulation.

4.1 Ocean Earth-like 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 spin-up 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 shape-preserving 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, $$M_{\text{error}}=\frac{M}{M_i}-1,$$(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⁻¹³%, and it does not present any constant drift, which highlights the robustness of the new formulation.

4.1.1 Zonal winds and temperature

The spin-up of the simulated atmosphere is mainly driven by the thermal forcing and the rotational effects. In Fig. 6 we show the longitudinal- and time-averaged temperature and zonal² winds. As expected, the temperatures are warmer at low latitudes, where the lapse rate³ 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 Held-Suarez Test (Held & Suarez 1994), which is driven by the surplus of heat provided from the condensation in the large-scale 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⁻¹ 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 mid-latitude jets in our simulation are slightly weaker (roughly 5 m s⁻¹ 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 longitudinal-averaged zonal winds (m s−1) and temperatures (K) for the ocean Earth-like experiment. The values were time-averaged for 1000 Earth days.

Fig. 7 Final time- and longitudinal-averaged mass stream function function (1010 kg s−1) for the ocean Earth-like experiment. The values were time-averaged 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: $$\Psi =\frac{2\pi R_{\text{planet}}\cos \varphi }{g}{\displaystyle {\int }_0^P\overline{[v]}}\text{d}p,$$(37)

where R_planet is the planet radius, ϕ is the latitude, g is the surface gravity, p is the atmospheric pressure and $\overline{[v]}$ is the zonal and time-averaged meridional wind⁴. The results show the well-known 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 Held-Suarez 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 mid-latitudinal jets that are moved into lower latitudes. The stronger up-welling motion by mean circulation is located inthe equatorial region at pressure levels around 800 mbar, where large-scale 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 Earth-like planet. In Fig. 8 we show the meridional eddy sensible heat flux and the eddy kinetic energy for the ocean Earth-like simulation. These two quantities are calculated from the eddy components⁵ 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 ($\overline{[v^{\prime }{T}^{\prime }]}$) 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$\overline{[u^{\prime }{u}^{\prime }+v^{\prime }{v}^{\prime }+w^{\prime }{w}^{\prime }]}$) shown in Fig. 8, reaches the maximum at the mid-latitude jet’s location, where the largest dynamical wave-activity is present. Both meridional heat flux and eddy kinetic energy obtained here are very similar to the results from the dry atmosphere Held-Suarez 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 mass-stream 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 finite-volume 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 longitudinal-averaged meridional eddy sensible heat flux (m s−1 K) and eddy kinetic energy (m2 s−2 K) for the ocean Earth-like experiment. The values were time-averaged for 1000 Earth days.

Fig. 9 Final time- and longitudinal-averaged vertical velocity (Pa s−1) between 15° S and 15° N for the ocean Earth-like simulation. The values were time-averaged for 1000 Earth days.

Fig. 10 Instantaneous latitude-longitude map of vertical velocity (Pa s−1) for the ocean Earth-like 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⁻³, 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 $k^{-\frac{5}{3}}$ 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 large-scale atmospheric models to remove kinetic energy near the truncation limit via, for example, explicit hyper-diffusion methods or numerical implicit schemes. Our spectrum in Fig. 11 shows a steeper slope than k⁻³ at higher wavenumbers as expected by the use of the fourth-order hyper-diffusion 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 non-physical artefacts arising from the sudden cut in resolution aliasing onto the well-resolved modes (smaller wavenumbers). Also, the physics at the truncation scales are not accurately solved by the atmospheric models. For an accurate representation of the large-scale physical phenomena, we should mitigate the inaccuracies by removing the energy at these scales.

Fig. 11 Time-averaged kinetic energy spectrum at 0.25 mbar for the ocean Earth-like 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 large-scale 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 mid-latitudes, 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 Earth-like 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 Earth-like 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 Earth-like simulation. The time step is 50 s, half of the value used for the ocean Earth-like 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 spin-up phase of the simulation.

Figure 14 shows the total mass error for the ocean tidally locked planet simulation. As for the aqua Earth-like case, the error measured along the simulation is small and not larger than 7.5×10⁻¹³%, 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 longitudinal-averaged water vapour concentration (kg/kg) and relative humidity (%) for the ocean Earth-like simulation. The values were time-averaged for 1000 Earth days.

Fig. 13 Instantaneous latitude-longitude map of surface precipitation (mm h−1) for the ocean Earth-like 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 Earth-like 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 sub-stellar 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 day-night 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 sub-stellar 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 large-scale 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 Earth-like 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. Large-scale 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 large-scale cells evolve in time due to the seasonality effects.

Fig. 15 Final time- and longitudinal-averaged zonal winds (m s−1) and temperatures (K) for the ocean tidally locked planet experiment. The values were time-averaged for 1000 Earth days.

Fig. 16 Final time- and longitudinal-averaged temperature (K) across the equator for the ocean tidally locked planet simulation. The values were time-averaged 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 Earth-like 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 longitudinal-averaged mass stream function (1010 kg s−1) for the ocean tidally locked planet simulation. The values were time-averaged 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 Earth-like planet case, the presence of the stronger components at shorter wavenumbers is clear. As in the Earth-like 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 Earth-like experiment resulted in a stronger generation of divergent kinetic energy. The slope near the truncation limit becomes larger due to the hyper-diffusion 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 Earth-like 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 longitudinal-averaged vertical velocity (Pa s−1) between 15° S and 15° N for the ocean tidally locked planet simulation. The values were time-averaged for 1000 Earth days.

Fig. 19 Time-averaged 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 large-scale 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 large-scale 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 upwind-biased 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 bell-shaped 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, shape-preserving properties, and avoiding negative solutions. The new solver has a second-order accuracy and avoids using any fixer, such as hyper-diffusion 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 1-degree-resolution 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 1-degree 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 Earth-like 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 large-scale 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⁻¹³%).

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 Earth-like 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 high-frequency 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 longitudinal-averaged water vapour concentration (kg/kg) and relative humidity (%) for the ocean tidally locked planet simulation. The values were time-averaged for 1000 Earth days.

Fig. 21 Instantaneous latitude-longitude 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

All Tables

All Figures

	Fig. 1 Initial cosine bell shape defined in Eq. (9). The colour map is in arbitrary units.
In the text

	Fig. 2 Latitude-longitude 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 finite-volume), 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 upwind-biased 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 upwind-biased 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 second-order convergence.
In the text

	Fig. 5 Evolution of the total mass error in the ocean Earth-like simulation.
In the text

	Fig. 6 Final time- and longitudinal-averaged zonal winds (m s−1) and temperatures (K) for the ocean Earth-like experiment. The values were time-averaged for 1000 Earth days.
In the text

	Fig. 7 Final time- and longitudinal-averaged mass stream function function (1010 kg s−1) for the ocean Earth-like experiment. The values were time-averaged for 1000 Earth days.
In the text

	Fig. 8 Final time- and longitudinal-averaged meridional eddy sensible heat flux (m s−1 K) and eddy kinetic energy (m2 s−2 K) for the ocean Earth-like experiment. The values were time-averaged for 1000 Earth days.
In the text

	Fig. 9 Final time- and longitudinal-averaged vertical velocity (Pa s−1) between 15° S and 15° N for the ocean Earth-like simulation. The values were time-averaged for 1000 Earth days.
In the text

	Fig. 10 Instantaneous latitude-longitude map of vertical velocity (Pa s−1) for the ocean Earth-like experiment.
In the text

	Fig. 11 Time-averaged kinetic energy spectrum at 0.25 mbar for the ocean Earth-like 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 longitudinal-averaged water vapour concentration (kg/kg) and relative humidity (%) for the ocean Earth-like simulation. The values were time-averaged for 1000 Earth days.
In the text

	Fig. 13 Instantaneous latitude-longitude map of surface precipitation (mm h−1) for the ocean Earth-like 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 longitudinal-averaged zonal winds (m s−1) and temperatures (K) for the ocean tidally locked planet experiment. The values were time-averaged for 1000 Earth days.
In the text

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

	Fig. 17 Final time- and longitudinal-averaged mass stream function (1010 kg s−1) for the ocean tidally locked planet simulation. The values were time-averaged for 1000 Earth days.
In the text

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

	Fig. 19 Time-averaged 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 longitudinal-averaged water vapour concentration (kg/kg) and relative humidity (%) for the ocean tidally locked planet simulation. The values were time-averaged for 1000 Earth days.
In the text

	Fig. 21 Instantaneous latitude-longitude map of surface precipitation (mm h−1) for the ocean tidally locked planet simulation.
In the text