A Global Nonhydrostatic Atmospheric Model with a Mass- and Energy-conserving Vertically Implicit Correction (VIC) Scheme

We start by describing our model Euler equations in 1D to present the derivation and formulation of our VIC scheme. Discretized using the simple forward-Euler time integration scheme, the governing equation can be given by

$$\begin{eqnarray}&&\displaystyle \frac{{{\boldsymbol{Q}}}_{i}^{n+1}-{{\boldsymbol{Q}}}_{i}^{n}}{{\rm{\Delta }}t}+\displaystyle \frac{{\sigma }_{i+1/2}{{\boldsymbol{F}}}_{i+1/2}^{n+1}-{\sigma }_{i-1/2}{{\boldsymbol{F}}}_{i-1/2}^{n+1}}{{\rm{\Delta }}{V}_{i}}={{\boldsymbol{X}}}_{i}^{n+1}+{{\boldsymbol{Y}}}_{i}^{n},\end{eqnarray} \tag{ 2 }$$

where ${{\boldsymbol{Q}}}_{i}^{n}$ represents the vector of volume-averaged conservative variables in the ith cell at the nth time step (i.e., current time step); ${{\boldsymbol{F}}}_{i-1/2}^{n+1}$ and ${{\boldsymbol{F}}}_{i+1/2}^{n+1}$ are numerical fluxes across the left (bottom) and the right (top) interfaces of the ith cell at the (n + 1)th time step (i.e., next time step), respectively; ${{\boldsymbol{Y}}}_{i}^{n}$ and ${{\boldsymbol{X}}}_{i}^{n+1}$ are explicitly and implicitly treated source terms, respectively. ΔV_i is the volume of the cell, which is constant in the Cartesian coordinate system but changes with latitude and radius in the spherical polar coordinate system. ${\sigma }_{i+1/2}$ and ${\sigma }_{i-1/2}$ are face areas of the ith cell in the vertical direction. On the right-hand side of the equation, most of the source terms, such as gravitational acceleration, Coriolis force, and centrifugal force, can be treated implicitly (i.e., included in ${{\boldsymbol{X}}}_{i}^{n+1}$) because they can be written analytically in terms of flux Jacobian, which is the first-order partial derivative of the forcing with respect to ${\boldsymbol{Q}}$. Radiative forcing is treated explicitly because it usually does not have an analytical Jacobian.

Once we know ${\rm{\Delta }}{{\boldsymbol{Q}}}_{i}={{\boldsymbol{Q}}}_{i}^{n+1}-{{\boldsymbol{Q}}}_{i}^{n}$, then we can update ${{\boldsymbol{Q}}}_{i}^{n+1}$ by ${{\boldsymbol{Q}}}_{i}^{n+1}={\rm{\Delta }}{{\boldsymbol{Q}}}_{i}+{{\boldsymbol{Q}}}_{i}^{n}$. Then, the problem becomes how to acquire the flux and source terms at the (n + 1)th time step, ${{\boldsymbol{F}}}_{i-1/2}^{n+1}$, ${{\boldsymbol{F}}}_{i+1/2}^{n+1}$, and ${{\boldsymbol{X}}}_{i}^{n+1}$, from the conservative variables at the nth step. Following the ideas in Fernandez (1988) and Viozat (1997), we use the Roe scheme to acquire the implicit fluxes across cell interfaces at the (n + 1)th time step (Roe 1981),

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{i+1/2}^{n+1}=\displaystyle \frac{1}{2}({{\boldsymbol{F}}}_{i}^{n+1}+{{\boldsymbol{F}}}_{i+1}^{n+1}-| {A}_{i+1/2}^{n+1}| ({{\boldsymbol{Q}}}_{i+1}^{n+1}-{{\boldsymbol{Q}}}_{i}^{n+1})),\end{eqnarray} \tag{ 3 }$$

where $| {A}_{i+1/2}^{n+1}|$ is the Roe-averaged flux Jacobian across the right interface of the ith cell (Roe 1981). This formulation can also be applied to the fluxes at the nth time level. $| {A}_{i+1/2}^{n}|$ performs as a "stabilization term" to stabilize the flux at the cell interface. It is different from the hyperviscosity that is implemented in many GCMs. It is commonly considered as a stabilization term in the Roe scheme. The detailed discussion of $| {A}_{i+1/2}^{n}|$ is well presented in Roe (1981), LeVeque et al. (2002, Chap. 4.14), and Toro (2013, Chap. 11). We provide the analytical derivation of $| {A}_{i+1/2}^{n}|$ in Appendix B.

Note that the Roe-averaged Jacobian at the (n + 1)th time step, $| {A}_{i+1/2}^{n+1}|$, is unknown at the nth time step. As a workaround, we consider the first-order approximation, $| {A}_{i+1/2}^{n}| =| {A}_{i+1/2}^{n+1}| +{ \mathcal O }({\rm{\Delta }}t)$, which can be justified by two reasons: (1) $| A|$ is dictated by the flow and acoustic speed, and they are locally uniform in the atmosphere, therefore the temporal evolution of $| A|$ remains trivial over a small time step Δt (Δt is usually much smaller than 10³ s); (2) the desired numerical stability of the Roe scheme is not severely compromised by the approximation of using $| {A}_{i+1/2}^{n}|$, since $| A|$ acts as the numerical diffusion term. We will demonstrate the validity of this approximation by benchmark simulations in Section 4.

To simplify the notation in this section, we define the flux Jacobian, ${{\boldsymbol{J}}}_{i}^{n}$, and the Jacobian of the source term, ${{\boldsymbol{X}}}_{i}^{{\prime} n}$, as

$$\begin{eqnarray}&&{{\boldsymbol{J}}}_{i}^{n}\equiv \displaystyle \frac{{\rm{\partial }}{\boldsymbol{F}}}{{\rm{\partial }}{\boldsymbol{Q}}}{| }_{i}^{n}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{i}^{{\rm{{\prime} }}n}\equiv \displaystyle \frac{{\rm{\partial }}{\boldsymbol{X}}}{{\rm{\partial }}{\boldsymbol{Q}}}{| }_{i}^{n}.\end{eqnarray} \tag{ 5 }$$

The numerical flux, ${{\boldsymbol{F}}}_{i}^{n+1}$, and the implicitly treated source terms, ${{\boldsymbol{X}}}_{i}^{n+1}$, at the (n + 1)th time step are unknown. To compute these values, we use a Taylor expansion to linearize the flux at the (n + 1)th time step. For example, we predict ${{\boldsymbol{F}}}_{i}^{n+1}$ and ${{\boldsymbol{X}}}_{i}^{n+1}$, at the (n + 1)th step from the nth time step by a second-order accurate Taylor expansion,

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{i}^{n+1}\approx {{\boldsymbol{F}}}_{i}^{n}+{{\boldsymbol{J}}}_{i}^{n}{\rm{\Delta }}{{\boldsymbol{Q}}}_{i}+{ \mathcal O }[{\left({\rm{\Delta }}{{\boldsymbol{Q}}}_{i}\right)}^{2}]\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{i}^{n+1}\approx {{\boldsymbol{X}}}_{i}^{n}+{{\boldsymbol{X}}}_{i}^{{\prime} n}{\rm{\Delta }}{{\boldsymbol{Q}}}_{i}+{ \mathcal O }[{\left({\rm{\Delta }}{{\boldsymbol{Q}}}_{i}\right)}^{2}],\end{eqnarray} \tag{ 7 }$$

where ${{\boldsymbol{F}}}_{i}^{n}$ and ${{\boldsymbol{X}}}_{i}^{n}$ are numerical flux and source terms at the nth time step, respectively, and ${\rm{\Delta }}{{\boldsymbol{Q}}}_{i}={{\boldsymbol{Q}}}_{i}^{n+1}-{{\boldsymbol{Q}}}_{i}^{n}$. As previously mentioned, our numerical scheme solves ${\rm{\Delta }}{{\boldsymbol{Q}}}_{i}$ to obtain ${{\boldsymbol{Q}}}_{i}^{n+1}$.

Note that Roe flux is not available at cell centers (i.e., ${{\boldsymbol{F}}}_{i}^{n}$ and ${{\boldsymbol{F}}}_{i+1}^{n}$ in Equation (6)) because fluxes are evaluated at the cell interfaces. These terms can be eliminated by Equation (3) at the nth time step to assemble the interface flux at time step n (i.e., ${{\boldsymbol{F}}}_{i+1/2}^{n}$),

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{i}^{n}+{{\boldsymbol{F}}}_{i+1}^{n}=2{{\boldsymbol{F}}}_{i+1/2}^{n}+| {A}_{i+1/2}^{n}| ({{\boldsymbol{Q}}}_{i+1}^{n}-{{\boldsymbol{Q}}}_{i}^{n}).\end{eqnarray} \tag{ 8 }$$

Then, we can assemble linear equations about ${\rm{\Delta }}{{\boldsymbol{Q}}}_{i}$ by combining Equations (6)–(8). As a consequence, the tridiagonal linear system can be rewritten as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{{\rm{\Delta }}{V}_{i}}\left[{c}_{i}{\rm{\Delta }}{{\boldsymbol{Q}}}_{i-1}+(\displaystyle \frac{{\boldsymbol{ \mathcal I }}}{{\rm{\Delta }}t}-{{\boldsymbol{X}}}_{i}^{{\prime} n}+{a}_{i}){\rm{\Delta }}{{\boldsymbol{Q}}}_{i}+{b}_{i}{\rm{\Delta }}{{\boldsymbol{Q}}}_{i+1}\right]\\ \ =\ -\displaystyle \frac{{\rm{\Delta }}{\sigma }_{i+1/2}{{\boldsymbol{F}}}_{i+1/2}^{n}-{\rm{\Delta }}{\sigma }_{i-1/2}{{\boldsymbol{F}}}_{i-1/2}^{n}}{{\rm{\Delta }}{V}_{i}}+{{\boldsymbol{X}}}_{i}^{n}+{{\boldsymbol{Y}}}_{i}^{n},\end{array}\end{eqnarray} \tag{ 9 }$$

with coefficients a_i, b_i, and c_i being

$$\begin{eqnarray}\begin{array}{rcl}{a}_{i} & = & \displaystyle \frac{1}{2}\left[[,({\sigma }_{i+1/2}-{\sigma }_{i-1/2}){{\boldsymbol{J}}}_{i}^{n}+{\sigma }_{i-1/2}| {A}_{i-1/2}^{n}| \right.\\ & & +\,\left.{\sigma }_{i+1/2}| {A}_{i+1/2}^{n}| ,]\right],\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{b}_{i}=\displaystyle \frac{{\sigma }_{i+1/2}}{2}\left({{\boldsymbol{J}}}_{i+1}^{n}-| {A}_{i+1/2}^{n}| \right),\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}&&{c}_{i}=-\displaystyle \frac{{\sigma }_{i-1/2}}{2}\left({{\boldsymbol{J}}}_{i-1}^{n}+| {A}_{i-1/2}^{n}| \right),\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{ \mathcal I }}$ is the identity matrix. When i = 1 and i = m, these equations shall be modified to satisfy different boundary conditions, which is detailed in Section 3 and Appendix C.

We use the fifth-order accurate weighted essentially non-oscillation (WENO5) reconstruction scheme and Riemann solvers in SNAP (Li & Chen 2019) to compute the Riemann state at each interface i ± 1/2, which are used to solve the Riemann problems at i ± 1/2 to calculate the Roe fluxes, ${{\boldsymbol{F}}}_{i+1/2}^{n}$ and ${{\boldsymbol{F}}}_{i-1/2}^{n}$. A Riemann solver has a good capability of capturing shocks, which might be helpful in studying climatology on exoplanets (Fromang et al. 2016). The left-hand side of Equation (9) is a tridiagonal linear system that can be solved via the simple Gaussian elimination method. Thus, all terms in Equation (9) can be evaluated at the nth time step. ${\rm{\Delta }}{{\boldsymbol{Q}}}_{i}$ is solved by inverting Equation (9), and finally the conservative variables at step n + 1 are

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{i}^{n+1}={\rm{\Delta }}{{\boldsymbol{Q}}}_{i}+{{\boldsymbol{Q}}}_{i}^{n}.\end{eqnarray} \tag{ 13 }$$

The previous section provided the formulation of the VIC scheme in a single-stage and forward-Euler time integration scheme. Here we discuss how to implement the VIC scheme in a multi-stage Runge–Kutta (RK)-type time integration scheme to achieve higher numerical stability and accuracy in time (Shu & Osher 1988). The formulation of this implicit scheme is slightly different from the total variation diminishing (TVD) method in the original paper (Li & Chen 2019). Here, we list algorithms of a third-order accurate Runge–Kutta time integration (RK3) method as an example. We also present both explicit and implicit algorithms for comparison. The original explicit RK3 time integration scheme is

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{Q}}}_{i}^{(1)} & = & {{\boldsymbol{Q}}}_{i}^{n}-\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}{V}_{i}}({\sigma }_{i+1/2}{{\boldsymbol{F}}}_{i+1/2}^{n}-{\sigma }_{i-1/2}{{\boldsymbol{F}}}_{i-1/2}^{n})\\ & & +\ {\rm{\Delta }}t({{\boldsymbol{X}}}_{i}^{n}+{{\boldsymbol{Y}}}_{i}^{n}),\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{ \mathcal Q }}}_{i}^{(2)} & = & \displaystyle \frac{3}{4}{{\boldsymbol{ \mathcal Q }}}_{i}^{n}+\displaystyle \frac{1}{4}\left[{{\boldsymbol{ \mathcal Q }}}_{i}^{(1)}-\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}{V}_{i}}\left({\sigma }_{i+1/2}{{\boldsymbol{ \mathcal F }}}_{i+1/2}^{(1)}\right.\right.\\ & & -\ \left.\left.{\sigma }_{i-1/2}{{\boldsymbol{ \mathcal F }}}_{i-1/2}^{(1)}\right)+{\rm{\Delta }}t({{\boldsymbol{ \mathcal X }}}_{i}^{(1)}+{{\boldsymbol{ \mathcal Y }}}_{i}^{(1)}),]\right],\end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{Q}}}_{i}^{n+1} & = & \displaystyle \frac{1}{3}{{\boldsymbol{Q}}}_{i}^{n}+\displaystyle \frac{2}{3}\left[{{\boldsymbol{Q}}}_{i}^{(2)}-\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}{V}_{i}}({\sigma }_{i+1/2}{{\boldsymbol{F}}}_{i+1/2}^{(2)}\right.\\ & & -\ \left.{\sigma }_{i-1/2}{{\boldsymbol{F}}}_{i-1/2}^{(2)})+{\rm{\Delta }}t({{\boldsymbol{X}}}_{i}^{(2)}+{{\boldsymbol{Y}}}_{i}^{(2)}),]\right],\end{array}\end{eqnarray} \tag{ 16 }$$

where the superscripts (1) and (2) represent the first and the second intermediate states, respectively. We use the intermediate states of conservative variables, ${{\boldsymbol{Q}}}_{i}^{(1)}$ or ${{\boldsymbol{Q}}}_{i}^{(2)}$, to update the flux at the intermediate states, ${{\boldsymbol{F}}}_{i+1/2}^{(1)}$ and ${{\boldsymbol{F}}}_{i+1/2}^{(2)}$.

The implicit scheme replaces the explicit fluxes ${\boldsymbol{F}}$ and forcing ${\boldsymbol{X}}$ by their implicit counterparts. For example, the first step of the RK3 integration scheme becomes

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{{\boldsymbol{Q}}}_{i}^{n}}{{\rm{\Delta }}t}+\displaystyle \frac{{\sigma }_{i+1/2}{{\boldsymbol{F}}}_{i+1/2}^{(1)}-{\sigma }_{i-1/2}{{\boldsymbol{F}}}_{i-1/2}^{(1)}}{{\rm{\Delta }}{V}_{i}}={{\boldsymbol{X}}}_{i}^{(1)}+{{\boldsymbol{Y}}}_{i}^{n}.\end{eqnarray} \tag{ 17 }$$

According to Equation (9), ${\rm{\Delta }}{{\boldsymbol{Q}}}_{i}^{n}$ is solved by inverting the tridiagonal system defined by Equation (17). Then, the conservative variables are updated by the following equation:

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{i}^{(1)}={\rm{\Delta }}{{\boldsymbol{Q}}}_{i}^{n}+{{\boldsymbol{Q}}}_{i}^{n}.\end{eqnarray} \tag{ 18 }$$

Similarly, the other two intermediate stages are updated sequentially:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{\Delta }}{{\boldsymbol{Q}}}_{i}^{(1)}}{{\rm{\Delta }}t/4}+\displaystyle \frac{{\sigma }_{i+1/2}{{\boldsymbol{F}}}_{i+1/2}^{(2)}-{\sigma }_{i-1/2}{{\boldsymbol{F}}}_{i-1/2}^{(2)}}{{\rm{\Delta }}{V}_{i}}={{\boldsymbol{X}}}_{i}^{(2)}+{{\boldsymbol{Y}}}_{i}^{(1)},\\ {{\boldsymbol{Q}}}_{i}^{(2)}={\rm{\Delta }}{{\boldsymbol{Q}}}_{i}^{(1)}+\displaystyle \frac{1}{4}{{\boldsymbol{Q}}}_{i}^{(1)}+\displaystyle \frac{3}{4}{{\boldsymbol{Q}}}_{i}^{n},\end{array}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{\Delta }}{{\boldsymbol{Q}}}_{i}^{(2)}}{2{\rm{\Delta }}t/3}+\displaystyle \frac{{\sigma }_{i+1/2}{{\boldsymbol{F}}}_{i+1/2}^{n+1}-{\sigma }_{i-1/2}{{\boldsymbol{F}}}_{i-1/2}^{n+1}}{{\rm{\Delta }}{V}_{i}}={{\boldsymbol{X}}}_{i}^{n+1}+{{\boldsymbol{Y}}}_{i}^{(2)},\\ {{\boldsymbol{Q}}}_{i}^{n+1}={\rm{\Delta }}{{\boldsymbol{Q}}}_{i}^{(2)}+\displaystyle \frac{2}{3}{{\boldsymbol{Q}}}_{i}^{(2)}+\displaystyle \frac{1}{3}{{\boldsymbol{Q}}}_{i}^{n}.\end{array}\end{eqnarray} \tag{ 20 }$$

In this section, we demonstrate how to combine the explicit scheme in the horizontal plane (x–y plane) with the implicit scheme in the vertical z-direction. In particular, ${\boldsymbol{H}}$ is the flux in the vertical direction and is treated implicitly. The explicit scheme is computationally efficient in one time step but subject to the conditional numerical stability. The computational efficiency of an explicit nonhydrostatic atmosphere model is generally limited by the acoustic speed and atmospheric vertical resolution. In contrast, the implicit scheme trades off the computational efficiency against numerical stability. Here we combine the advantages of the two numerical schemes to solve the 3D Euler equations.

Most of the existing atmosphere models with the vertically implicit scheme use a dimensionally split method to compute horizontal and vertical terms separately (e.g., Ullrich & Jablonowski 2012a, 2012b; Mendonça et al. 2016). Flux divergence in different directions is updated by different time integration schemes in some algorithms (e.g., Bao et al. 2015). Here, we use a dimensionally unsplit method to update conservative variables simultaneously in all x-, y-, and z-directions for each half time step and the final time step. The discretized equation using the dimensionally unsplit formulae for the 3D case is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{\boldsymbol{ \mathcal P }}}_{{ijk}}^{n+1}-{{\boldsymbol{ \mathcal P }}}_{{ijk}}^{n}}{{\rm{\Delta }}t}+\displaystyle \frac{{\sigma }_{i+1/2;;}{{\boldsymbol{ \mathcal F }}}_{i+1/2;;}^{n}-{\sigma }_{i-1/2;;}{{\boldsymbol{ \mathcal F }}}_{i-1/2;;}^{n}}{{\rm{\Delta }}{V}_{{ijk}}}\\ \,+\,\displaystyle \frac{{\sigma }_{;j+1/2;}{{\boldsymbol{ \mathcal G }}}_{;j+1/2;}^{n}-{\sigma }_{;j-1/2;}{{\boldsymbol{ \mathcal G }}}_{;j-1/2;}^{n}}{{\rm{\Delta }}{V}_{{ijk}}}\\ \,+\,\displaystyle \frac{{\sigma }_{;;k+1/2}{{\boldsymbol{ \mathcal H }}}_{;;k+1/2}^{n+1}-{\sigma }_{;;k-1/2}{{\boldsymbol{ \mathcal H }}}_{;;k-1/2}^{n+1}}{{\rm{\Delta }}{V}_{{ijk}}}\\ \,={{\boldsymbol{ \mathcal X }}}_{{ijk}}^{n+1}+{{\boldsymbol{ \mathcal Y }}}_{{ijk}}^{n},\end{array}\end{eqnarray} \tag{ 21 }$$

where ${\boldsymbol{F}}$, ${\boldsymbol{G}}$, and ${\boldsymbol{H}}$ are numerical fluxes in x-, y-, and z-directions (or θ, ϕ, and r in a spherical coordinate system), respectively. Note that only the vertical (or radial) flux gradient is treated implicitly. We simplified the notation of the subscript by using ";" to omit some cell-centered indices. For example, ${{\boldsymbol{Q}}}_{i+1/2;;}\equiv {{\boldsymbol{Q}}}_{i+1/2,j,k}$.

Substituting Equations (6) and (3) into Equation (21), the discretized Euler equations can be rewritten as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{{\rm{\Delta }}{V}_{{ijk}}}\left[{c}_{{ijk}}({{\boldsymbol{Q}}}_{;;k-1}^{n+1}-{{\boldsymbol{Q}}}_{;;k-1}^{n})+\left(\displaystyle \frac{{\boldsymbol{ \mathcal I }}}{{\rm{\Delta }}t}-{{\boldsymbol{X}}}_{{ijk}}^{{\prime} n}+{a}_{{ijk}}\right)\right.\\ \ \times \ \left.({{\boldsymbol{Q}}}_{{ijk}}^{n+1}-{{\boldsymbol{Q}}}_{{ijk}}^{n})+{b}_{{ijk}}({{\boldsymbol{Q}}}_{;;k+1}^{n+1}-{{\boldsymbol{Q}}}_{;;k+1}^{n}),]\right]={{\boldsymbol{R}}}_{{ijk}},\end{array}\end{eqnarray} \tag{ 22 }$$

where ${{\boldsymbol{R}}}_{{ijk}}$ is

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{R}}}_{{ijk}} & = & -\displaystyle \frac{{\sigma }_{i+1/2;;}{{\boldsymbol{F}}}_{i+1/2;;}^{n}-{\sigma }_{i-1/2;;}{{\boldsymbol{F}}}_{i-1/2;;}^{n}}{{\rm{\Delta }}{V}_{{ijk}}}\\ & & -\ \displaystyle \frac{{\sigma }_{;j+1/2;}{{\boldsymbol{G}}}_{;j+1/2;}^{n}-{\sigma }_{;j-1/2;}{{\boldsymbol{G}}}_{;j-1/2;}^{n}}{{\rm{\Delta }}{V}_{{ijk}}}\\ & & -\displaystyle \frac{{\sigma }_{;;k+1/2}{{\boldsymbol{H}}}_{;;k+1/2}^{n}-{\sigma }_{;;k-1/2}{{\boldsymbol{H}}}_{;;k-1/2}^{n}}{{\rm{\Delta }}{V}_{{ijk}}}+{{\boldsymbol{X}}}_{{ijk}}^{n}+{{\boldsymbol{Y}}}_{{ijk}}^{n}\\ & = & \displaystyle \frac{{{\boldsymbol{Q}}}_{{ijk}}^{(e)}-{{\boldsymbol{Q}}}_{{ijk}}^{n}}{{\rm{\Delta }}t}.\end{array}\end{eqnarray} \tag{ 23 }$$

${{\boldsymbol{Q}}}_{{ijk}}^{(e)}$ denotes the conservative variables predicted by the explicit formulation. The time integration scheme for the 3D case is similar to the 1D time integration scheme in Section 2.2. We observe that, for some 3D simulations, the horizontal Courant–Friedrichs–Lewy (CFL) number can reach as large as 1.6 (i.e., c_sΔt/min(Δx, Δy) ≈ 1.6) with the implicit RK3 time integration scheme. The reason for this "ultrastable" behavior is probably due to the unsplit nature of our implicit scheme, such that relaxing the theoretical CFL limit in one direction helps to improve numerical stability in the other directions. A rigorous numerical study of the VIC scheme will be undertaken in forthcoming studies. Finally, we summarize the whole implicit scheme in the following steps:

• 1.  
  Perform an explicit forward step and calculate ${{\boldsymbol{R}}}_{{ijk}}$ according to Equation (23).
• 2.  
  Calculate coefficients of the block tridiagonal matrix (a_ijk, b_ijk, c_ijk) according to Equations (10)–(12).
• 3.  
  Solve the block tridiagonal matrix column by column as defined by Equation (22).
• 4.  
  Update the conserved variables according to Equation (13).
• 5.  
  Repeat the preceding steps for each stage of a multi-stage time integration scheme similar to Equations (18)–(20).

Because steps 1, 4, and 5 are needed for any explicit integration scheme, an existing explicit model can employ our implicit scheme by simply adding additional implicit correction steps outlined in step 2 and step 3, which makes the scheme extremely flexible and versatile.

The VIC scheme relaxes the CFL limitation by damping fast acoustic waves. In most benchmark tests that are designed for atmospheric simulations, the vertical wind velocity is slower than the acoustic speed by more than two orders of magnitude, meaning that the time step is predominantly limited by the acoustic wave speed. The VIC scheme resolves this issue by imposing numerical diffusivity, which suppresses the accumulation of the spurious numerical noise.

We validate the model by simulating the propagation of the linearized acoustic wave. This test is very similar to the linear wave convergence test in Stone et al. (2008) but we focus on how numerical diffusivity damps the amplitude of the acoustic wave. Linear acoustic waves are launched by an initial pressure perturbation. The polytropic index, γ, is set as 5/3. Pressure, density, and velocity of the uniform background are initialized with 1.0 Pa, 0.6 ${\rm{g}}\,{\mathrm{cm}}^{-3}$, and 0 ${\rm{m}}\,{{\rm{s}}}^{-1}$, respectively. The corresponding acoustic speed, c_s, is 1 ${\rm{m}}\,{{\rm{s}}}^{-1}$. The wave amplitude is set as 10⁻³ Pa. The wave propagates in a 1D double-periodic tube whose wavelength is equal to the length of the domain. We simulate the propagation of acoustic waves with different time steps (i.e., CFL numbers) and collect the density profiles after 1 s in the simulation.

The left plot in Figure 1 shows the density profiles after one wave period with different CFL numbers but the same resolution. The numerical diffusivity damps the acoustic wave amplitude significantly for large CFL numbers. The usage of the large CFL number also cause the dispersion error for the acoustic waves. The right plot in Figure 1 shows that the damping rate grows almost linearly with increasing CFL number. The increase of the damping rate can be inferred by analyzing Equation (9). The numerical diffusivity is imposed by the off-diagonal elements in Equation (9). They are independent of the changing time step and depend only on the spatial resolution. When we fix the spatial resolution and increase Δt, the off-diagonal elements become more and more dominant because 1/Δt in diagonal terms reduces. Therefore, the damping rate increases linearly with increasing CFL number.

Zoom In Zoom Out Reset image size

We have shown that the numerical diffusivity imposed by the VIC can damp the acoustic wave amplitude from the linearized acoustic wave test. This helps the VIC to achieve numerical stability when a large time step is adopted. However, large numerical diffusivity may hinder the ability to resolve turbulence. It is important to necessary the VIC scheme's ability to correctly track turbulent motions. In this and the next section, we will focus on the VIC scheme's capability of resolving nonlinear density currents and turbulent motions.

The Straka sinking bubble test is a standard benchmark test, which is designed for the validation of nonhydrostatic atmospheric dynamical cores (Straka et al. 1993). This case simulates the fluid motion in a nonlinear density current generated by a sinking cold bubble. Several physical processes (i.e., Kelvin–Helmholtz instability) should be produced after the bubble is dropped on the surface (Straka et al. 1993). The simulation is carried out in a 25.6 km by 6.4 km closed domain. The initial background temperature structure is set as dry adiabat in the vertical direction. The reference surface pressure is set as 1 bar (10⁵ Pa). A cold bubble, whose center is 15 K colder than the local background temperature, is put aloft. The cold bubble is expected to drop to the surface and generate several nonlinear density currents. In the original test, the explicit diffusion is also applied to momentum and energy equations to guarantee convergence (Straka et al. 1993). In this work, however, we do not apply the diffusion terms in our equations in order to check the performance difference between our VIC scheme and the explicit scheme. The initial temperature is given by

$$\begin{eqnarray}&&\begin{array}{c}{\rm{\Delta }}T=\left\{\begin{array}{c}0,\,{\rm{i}}{\rm{f}}\,L\gt 1,\\ -15[\cos (\pi L)+1]/2,\,{\rm{i}}{\rm{f}}\,L\leqslant 1,\end{array}\right.\end{array}\end{eqnarray} \tag{ 26 }$$

where $L=\{{[(x-{x}_{c})/{x}_{r}]}^{2}+{[(z-{z}_{c})/{z}_{r}]}^{2}\}{}^{1/2}$, ${x}_{c}=0$ km, ${x}_{r}=4$ km, ${z}_{c}=3$ km, and ${z}_{r}=2$ km.

For the explicit test, we adopt the same numerical techniques as adopted in Li & Chen (2019) using the Low Mach Number Approximate Riemann Solver (LMARS), WENO5 reconstruction scheme, and RK3 time integration scheme, to ensure our result is comparable with the numerical solutions in Li & Chen (2019). We first investigate our VIC scheme's capability of capturing the position of the bubble debris and Kelvin–Helmholtz rotors. Then, we compare the potential temperature difference between the results of the VIC scheme and the explicit SNAP from Li & Chen (2019) to estimate the role of the numerical diffusivity.

Simulation results at t = 900 s with an aspect ratio of 1 (i.e., Δx/Δz = 1) are presented in Figure 2. The maximum CFL numbers are 0.5 and 1.0 for the explicit and implicit schemes, respectively. Here, we present and compare the results of using the explicit scheme with CFL  = 0.4, the VIC scheme with CFL  = 0.4, and the VIC scheme with CFL = 0.8 in Figure 2. Simulation results show that they generally converge to the same solution morphologically. Three Kelvin–Helmholtz rotors are correctly produced in all three cases. The position of the density currents' outer limb is about 15,300 m in all results. Slight potential temperature differences can be seen from the plots. They are caused by numerical diffusion. The potential temperature deviations of the coldest part of the bubble debris are roughly 11.9 K, 13.4 K, and 12.9 K for the three cases, respectively. For the explicit cases, the numerical diffusivity becomes smaller when the CFL number becomes larger, whereas the numerical diffusivity caused by the VIC scheme becomes larger when the time step becomes larger. This is expected because the non-diagonal terms, which act as the numerical diffusivity, are more dominant than the diagonal terms when the time step becomes larger.

Zoom In Zoom Out Reset image size

We also present two innovative simulation results with the large aspect ratio (i.e., Δx/Δz = 10) in Figure 3. The large aspect ratio of the vertical resolution and horizontal resolution allows us to use a much larger CFL number for the VIC scheme. The VIC scheme adopts a time step that is larger than the time step in the explicit case by a factor of 20. The dynamical evolution of the bubble is not affected by the usage of a much larger time step in the VIC scheme but the computational efficiency is significantly improved. It takes about 31 minutes for the explicit scheme to finish the computation on Pleiades with 32 CPU cores (Sandy Bridge processors with frequency 2.6 GHz). The VIC scheme, on the other hand, only takes about 3 minutes, illustrating that the computational efficiency is improved by one order of magnitude with no influence on dynamics.

Zoom In Zoom Out Reset image size

Figure 4 shows the evolution of total mass and energy as a function of time for cases using the large aspect ratio. This result is presented to show the performance of the conservation law on both the VIC and the explicit schemes. It shows that total mass linearly increases as a function of time on the machine precision level in both explicit and implicit cases. The simulation using the VIC scheme has a smaller mass error than the explicit one due to the smaller workload and fewer time steps. The fractional variation of total energy in both cases is about 10⁻⁹ of the initial value.

Zoom In Zoom Out Reset image size

In sum, the Straka sinking bubble test shows that, despite the small and dynamically trivial potential temperature differences among these cases, our VIC scheme can reproduce numerical results in the explicit integration. It also shows that the VIC scheme can well conserve total mass and energy in a closed system (e.g., with reflecting boundary conditions).

The third benchmark test is the Robert rising bubble test, which is designed to test the model performance under a weak forcing. The simulation results of the Straka sinking bubble test show that the VIC scheme can correctly track the fluid motion in the gravity-dominant regime. However, further validation of turbulence tracking is necessary for the regime with weak buoyant forcing. The Robert rising bubble test is a suitable benchmark test for our needs (Robert 1993). This benchmark test has some practical implications because there are ubiquitous convective air parcels in the moist convective layer (i.e., troposphere) and, unlike the situation in the Straka bubble test, atmospheric convection is usually triggered by a small temperature perturbation (i.e., less than 1 K). The simulation is initialized with a Gaussian-shaped warm bubble and finished with a very turbulent snapshot (Robert 1993).

Similar to the Straka benchmark test, the ambient atmospheric temperature is adiabatic with 303.15 K at the surface in a 1.5 km by 1 km closed box. The central temperature of the Gaussian-shaped bubble is set to be 0.5 K warmer than the background. The analytical formulation of the initial temperature is calculated from the potential temperature

$$\begin{eqnarray}&&\begin{array}{c}{\rm{\Delta }}\theta =\left\{\begin{array}{c}A,\,{\rm{i}}{\rm{f}}\,r\leqslant a,\\ A\,{e}^{-{\left(r-a\right)}^{2}/{s}^{2}},\,{\rm{i}}{\rm{f}}\,r\gt a,\end{array}\right.\end{array}\end{eqnarray} \tag{ 27 }$$

where ${r}^{2}={(x-{x}_{0})}^{2}+{(z-{z}_{0})}^{2}$, A = 0.5 K, a = 50 m, s = 100 m, ${x}_{0}=500$ m, and ${z}_{0}=260$ m. The warm bubble is expected to rise upward due to buoyancy. Robert (1993) shows that the uplifting motion stops at about 18 minutes (simulation time) and the Kelvin–Helmholtz instability is developed at the tail of the warm air parcel (Robert 1993). It is necessary to apply a fine resolution, at least 5 m in horizontal or vertical, to resolve the Kelvin–Helmholtz instability at 18 minutes (Robert 1993).

We present the simulation results with a large aspect ratio (i.e., ${\rm{\Delta }}x/{\rm{\Delta }}z=10$) in Figure 5. Similar to the Straka sinking bubble test, this setup allows us to test the performance of our VIC scheme in a time step an order of magnitude larger. The implicit numerical solutions are mostly identical to explicit results. Our VIC scheme can reproduce the same turbulent patterns and Kelvin–Helmholtz instability as in the explicit solution at 18 minutes. The predominant difference is still the potential temperature. The head of the bubble showing in the implicit result is also slightly larger than in the explicit one. In general, this simulation result agrees with the solutions reported in the previous studies (Robert 1993; Chen et al. 2013; Guerra & Ullrich 2016; Li & Chen 2019). Although our implicit scheme possesses a strong numerical diffusion for large CFL numbers, the numerical diffusivity specifically suppresses the growing numerical noise of acoustic waves without any significant influence on turbulent flows induced by the very weak forcing.

Zoom In Zoom Out Reset image size

The simulation efficiency is significantly improved by the VIC scheme for simulations with large aspect ratio. Both explicit and implicit simulations use 50 CPU cores on Pleiades (Sandy Bridge processors) for the computation. The VIC scheme allows a time step that is a factor of 20 larger than in the explicit case. The implicit scheme is faster than the explicit scheme by a factor of eight in this case.

It is important for an atmospheric dynamical core to correctly resolve the dispersion relation of gravity waves because it preserves the information on how much energy and momentum are conveyed by waves. For example, inertial gravity waves play important roles in the momentum budget in the middle atmosphere, such as generating the quasi-biennial oscillation in Earth's stratosphere (Baldwin et al. 2001) and the quasi-quadrennial oscillation in Jupiter's stratosphere (Cosentino et al. 2017). Here, we present simulation results of a gravity wave benchmark test, which were first introduced by Skamarock & Klemp (1994).

The model setup is similar to Skamarock & Klemp (1994) and Chen et al. (2013). We adopt a stable and hydrostatic atmosphere in a 300 km × 10 km domain with a small initial potential temperature perturbation to launch a train of gravity waves. The background atmosphere is initialized with a constant buoyancy frequency (i.e., Brunt–Väisälä frequency), $N={10}^{-2}$ s⁻¹. The surface pressure is set as 1 bar. The analytical potential temperature perturbation is

$$\begin{eqnarray}&&{\rm{\Delta }}\theta ={\rm{\Delta }}{\theta }_{0}\displaystyle \frac{\sin (\pi z/H)}{1+{\left(x-{x}_{c}\right)}^{2}/{a}^{2}},\end{eqnarray} \tag{ 28 }$$

where ${\rm{\Delta }}{\theta }_{0}=0.01$ K is the maximum potential temperature perturbation; ${x}_{c}=100$ km, a = 5 km, and H is the domain height, 10 km. A prescribed horizontal wind profile, u = 20 ${\rm{m}}\,{{\rm{s}}}^{-1}$, is initialized as the background wind for the wave propagation. Thus one set of gravity waves propagates eastward with respect to the reference wind field, while another set of waves is quasi-stationary.

The reflecting boundary condition is adopted for the lower and upper boundaries. The double-periodic condition is used in the horizontal direction. Our setup is different from the original test (Skamarock & Klemp 1994) but is close to the simulations in Chen et al. (2013) and Bao et al. (2015). We simulate nonhydrostatic wave activities with density variations instead of using an anelastic assumption in the original paper (Skamarock & Klemp 1994). Furthermore, similar to Chen et al. (2013), we do not include Coriolis forces in the simulation to simplify the problem.

The simulation results in a train of dispersive gravity waves at 3000 s as shown in Figure 6. Both explicit and implicit schemes can accurately resolve the phase and dispersion relation of the wave. This is crucial for simulating the propagation of momentum and energy. On the other hand, gravity waves with a small amplitude near $\sim 160$ km are not resolved using our VIC scheme, indicating that the VIC scheme has a difficulty resolving very weak perturbations (i.e., ${\rm{\Delta }}\theta \sim 0.001$ K). Furthermore, the explicit scheme can preserve the symmetrical structure of waves, but it seems that the performance of our VIC scheme is not as good as that of the explicit scheme under weak perturbations (i.e., at about $\sim$ 170 km in Figure 6(b)).

Zoom In Zoom Out Reset image size

The Held–Suarez experiment is designed for the intercomparison of GCM dynamical cores to produce the global-scale atmospheric features of an Earth-like planet (Held & Suarez 1994). This ideal climatology test focuses on the long-term and statistically averaged final state. The final quasi-steady state of this test is similar to an Earth-like atmosphere with a latitudinal temperature gradient and large-scale winds. The original work applied Newtonian cooling and Rayleigh drag schemes to simplify the radiative forcing and boundary effects, respectively. The fluid motion is quasi-geostrophic under the equator-to-pole temperature gradient and the Coriolis force with Earth's rotational rate. In the long-term averaged equilibrium state, thermal wind theory predicts two sub-tropical jets.

The initial condition of our simulation is set as a quasi-hydrostatic atmosphere with random and small temperature perturbations deviating from dry adiabats to break the initial symmetry. The thermal structure is relaxed to a reference profile via a Newtonian cooling scheme,

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+...\,=\,-{k}_{T}(\phi ,\sigma )\rho {c}_{v}[T-{T}_{\mathrm{eq}}(\phi ,p)]-\rho {wg},\end{eqnarray} \tag{ 29 }$$

where c_v is the isochoric specific heat; σ is the ratio of the local pressure, p, to the surface pressure, p_s; k_T is the temperature damping strength as a function of the latitude and $\sigma ;$ T_eq is the reference temperature profile as a function of latitude and σ. K_T and T_eq are adopted from the original paper (Held & Suarez 1994).

Rayleigh drags are applied to all three momentum equations,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}(\rho {\boldsymbol{u}})}{{\rm{\partial }}t}+...\,=\,-{k}_{v}(\sigma )\rho {\boldsymbol{u}}-2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{u}}+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 30 }$$

where k_v is Rayleigh drag strength and ${\boldsymbol{\Omega }}$ is the planetary rotation vector. Geometric terms shall be written on the left-hand side of the equation. We treat both Newtonian cooling and Rayleigh drag friction as the explicit body forcing (i.e., explicit source terms), using the same parameters as used in the original paper (Held & Suarez 1994).

The simulation domain is 25 km in height, latitude from −90°S to 90°N, and the longitude from $0^\circ$ to $360^\circ$. We adopt spherical polar coordinates with the latitude and longitude divided uniformly in degrees. The polar wedge scheme in Athena⁺⁺ (Stone et al. 2020) is applied to the polar region as the boundary conditions at the poles. Two solid-wall boundaries (i.e., reflecting or no-flux boundary condition) are applied to both the bottom and the top of the domain (Table 1). The spatial resolution is about 28 in latitude (64 cells), $2\buildrel{\circ}\over{.} 8$ in longitude (128 cells), and 625 meters in height (40 layers). A 5 km thick sponge layer (e.g., Rayleigh drag layer) is applied at the top of the atmosphere for both explicit and implicit schemes to absorb the vertically propagating acoustic and gravity waves. The CFL number for the explicit scheme is fixed as 0.3, which is calculated as CFL $=\,{c}_{s}{\rm{\Delta }}t/{\rm{\Delta }}z$. We adopt CFL = 0.5 for the VIC scheme, in which case the CFL number is limited by the horizontal resolution, satisfying CFL $=\,{c}_{s}{\rm{\Delta }}t/\min ({\rm{\Delta }}x,{\rm{\Delta }}y)$.

Similar to the original work, the simulation reaches the quasi-equilibrium state after about 200 simulation days given that the radiative cooling and Rayleigh drag timescales are about tens of days (Held & Suarez 1994). Simulation results are averaged from 200 to 1200 days.

First, we present the averaged state of zonal-mean zonal wind, meridional wind, and temperature in Figure 7. We also present the eddy analysis of this test in Figure 8 with the zonal-mean meridional eddy momentum flux and zonal-mean eddy kinetic energy. The maximum prograde jet speeds for explicit and implicit models are 30.97 ${\rm{m}}\,{{\rm{s}}}^{-1}$ and 30.84 ${\rm{m}}\,{{\rm{s}}}^{-1}$, respectively. The wind structures and jet speeds are very close to previously published simulations (e.g., Held & Suarez 1994; Lin 2004; Ullrich & Jablonowski 2012b; Mayne et al. 2014a; Mendonça et al. 2016).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Second, in Figure 9, we show the temporal evolution of total mass, energy, and AAM of both the explicit and implicit simulation results. The fractional variation of total mass increases linearly at the machine precision level (i.e., $\sim {10}^{-9}$ of the initial value) at a rate of $\sim 3\times {10}^{-9}$ yr⁻¹ and $\sim 3\times {5}^{-10}$ yr⁻¹ for explicit and implicit results, respectively. The changes in total mass are similar to the result of the Straka sinking bubble test, in which the VIC scheme results have a smaller change in total mass. This shows that the VIC scheme can conserve total mass just like the explicit scheme in a spherical polar coordinate system. The accumulated machine precision error of total mass is negligible given that the change in total mass is much smaller than ${10}^{-6}$ even if the model with the VIC scheme is integrated for thousands of years. The total energy varies at the start by about 0.6% and reaches a quasi-steady state after 200 days. The fractional variation of total energy is less than 10⁻³ after 200 days. The total AAM decreases from the initial value by about 0.1% and varies by about 0.2% after 200 days. Note that the total AAM is not expected to be rigorously conserved because we solve momentum equations instead of angular momentum equations. This result shows that our simulation does achieve the steady state after 200 days.

Zoom In Zoom Out Reset image size

In sum, our model can reproduce the benchmark result for Earth-like global-scale atmospheric dynamics. The VIC scheme achieves the mass conservation to the machine precision level. In the steady state, the total energy and AAM vary at a level below 1%.

Similar to previous tests, the VIC scheme can significantly improve computational efficiency. First, the time step increases by a factor of $\sim 450$ with the VIC scheme. Second, with the use of 512 CPU cores on Pleiades (Sandy Bridge processors), the computational efficiency is improved by more than two orders of magnitudes. The polar convergence issue of the spherical polar coordinate system still limits the spatial resolution in the polar region and therefore limits the time step and the computational efficiency, which can be improved by implementing a cubed-sphere coordinate system (Putman & Lin 2007) or using static mesh refinement (Zhu & Stone 2018) in our future studies.

Our final global-scale test is on the atmospheres of hot Jupiters, which are Jovian-size extrasolar planets very close to their host stars. These giant planets are likely to be synchronously orbiting around their host stars due to strong gravitational tides, meaning that the dayside of such a planet is always irradiated by its host star. Even though most hot Jupiters could still be rapid rotators (rotational period of about three days) like Jupiter, the strong, permanent day–night irradiation patterns distinguish the climate state on these emerging planetary populations from any previously known planetary atmospheres in the solar system. Observations and theories (e.g., Showman & Guillot 2002; Knutson et al. 2008; Showman et al. 2009) have shown that the equatorial super-rotating jet on hot Jupiters can be subsonic or sonic and the dayside and nightside temperature contrast can exceed 1000 K. The atmospheres of some ultrahot Jupiters can be hotter than 2500 K, so that hydrogen molecules on the dayside can be thermally dissociated (Bell & Cowan 2018; Tan & Komacek 2019). To prepare our model for future application to diverse exoplanetary atmospheres, here we validate our global GCM in this new regime of a synchronously rotating hot Jupiter. In particular, we validate our model against an experiment called the shallow hot Jupiter test. Several previous models have performed this hot Jupiter benchmark test and shown relatively good agreement (Menou & Rauscher 2009; Heng et al. 2011; Bending et al. 2012; Mayne et al. 2014b, 2017; Mendonça et al. 2016). But to date the published results were from either the model using primitive equations, hydrostatic models (Menou & Rauscher 2009; Heng et al. 2011; Bending et al. 2012) or a nonhydrostatic model with implicit and semi-implicit schemes (Mendonça et al. 2016; Mayne et al. 2017). In addition to testing our VIC scheme, here we also provide the first nonhydrostatic, explicit integration results for this case without damping the acoustic waves.

Similar to previous works, we simulate a canonical hot Jupiter with a planetary radius of ${10}^{5}$ km and a gravitational acceleration of 8 ${\rm{m}}\,{{\rm{s}}}^{-2}$. The simulation domain covers 4500 km in height with 40 layers, $-90^\circ {\rm{S}}$ to $90^\circ {\rm{N}}$ in latitude with 64 cells, and $0^\circ$ to $360^\circ$ in longitude with 128 cells. The bottom pressure and temperature are set as 1 bar and 1600 K, respectively. The initial temperature structure is isothermal (1600 K) from the surface to the top. The atmospheric temperature is relaxed to a reference temperature profile by a Newtonian cooling scheme, but no Rayleigh drag is applied. We adopt a 500 km thick sponge layer at the top of the domain to absorb vertically propagating waves. The same Newtonian cooling setup is used as in previous works. The reference temperature profile is given by

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}={T}_{\mathrm{vert}}+{\beta }_{\mathrm{trop}}{\rm{\Delta }}{T}_{E-P}\cos \lambda \cos \phi ,\end{eqnarray} \tag{ 31 }$$

where λ is the longitude, ϕ is the latitude, and T_vert and ${\beta }_{\mathrm{trop}}$ are parameters adopted from Menou & Rauscher (2009), Heng et al. (2011), and Mendonça et al. (2016). The temperature difference between equator and pole, ${T}_{E-P}$, is 300 K.

We choose the time steps corresponding to CFL = 0.3 for the explicit case and CFL = 0.9 for the implicit case.

Figure 10 shows the zonal-mean zonal wind, meridional wind, and temperature structure in the statistically averaged state from Day 200 to Day 1200. Figure 11 shows the zonal-mean meridional eddy momentum flux, vertical eddy momentum flux, and eddy kinetic energy. The zonal-mean zonal wind pattern exhibits a prograde equatorial super-rotating jet and retrograde jets in the mid-latitude. The zonal-mean temperature pattern in the lower atmosphere shows a strong latitudinal variation, with two sub-tropical maxima and two polar maxima. The results from the explicit and VIC cases are almost identical. The maximum zonal-mean prograde zonal wind speeds are 1052.95 ${\rm{m}}\,{{\rm{s}}}^{-1}$ and 1109.61 ${\rm{m}}\,{{\rm{s}}}^{-1}$ for explicit and implicit simulations, respectively. The location of the super-rotating jet ranges from about $-20^\circ {\rm{S}}$ to $20^\circ {\rm{N}}$. The difference in jet speed is also associated with a slight difference in the equatorial temperature structure shown in Figure 10. The maximum westward jet speeds are −588.41 ${\rm{m}}\,{{\rm{s}}}^{-1}$ and −581.59 ${\rm{m}}\,{{\rm{s}}}^{-1}$ in the explicit and implicit cases, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The zonal-mean zonal wind and zonal-mean temperature generally show similar structures to previous results. Both the explicit and VIC schemes produce a slightly slower equatorial jet than previous works (Menou & Rauscher 2009; Heng et al. 2011; Bending et al. 2012; Mayne et al. 2014b, 2017; Mendonça et al. 2016). The zonal-mean temperature structure at the equator is also slightly different from previous models (note that these differences also exist among previous models). These subtle differences might come from different capabilities of resolving eddies in different models. Although the agreement between our implicit and explicit model results demonstrates a good performance of our VIC scheme in exoplanet simulations, future nonhydrostatic models with explicit time integration schemes would be needed to validate our explicit model results.

It was also recently argued that simulations of the deep atmospheres of hot Jupiters might require a very long time integration to achieve the steady state (Mayne et al. 2017; Deitrick et al. 2020; Sainsbury-Martinez et al. 2019; Wang & Wordsworth 2020). The shallow hot Jupiter benchmark case has a relatively short radiative timescale and is expected to converge early. We also integrated our case for 12,000 days to ensure that simulations have achieved the steady state but did not find a significant change after 1200 days.

Figure 12 shows the change in total mass, energy, and AAM for the shallow hot Jupiter test in the first 1200 days. The explicit and implicit schemes show similar behavior (Figure 12). Similar to the Held–Suarez test, the total mass is well conserved at the machine precision level with a slight, linear increase with time. In the steady state, total energy varies at a level below 1% and total AAM varies within a few per cent. Both the total energy and AAM at the steady state are different from their initial values by a few per cent but show no trend of continuous loss or increase, even after long-term integration of 12,000 days.