alexa Implementation of a Semi-implicit Time Integration Scheme in Non-Hydrostatic Euler Equations

ISSN: 2168-9679

Journal of Applied & Computational Mathematics

  • Research Article   
  • J Appl Computat Math 6:369, Vol 6(4)
  • DOI: 10.4172/2168-9679.1000369

Implementation of a Semi-implicit Time Integration Scheme in Non-Hydrostatic Euler Equations

Nam H* and Choi SJ
Numerical Model Team, Korea Institute of Atmospheric Prediction Systems, 4F, Hankuk Computer Building, 35 Boramae-Ro 5 GIL, Dongjak-Gu, Seoul 07071, South Korea
*Corresponding Author: Nam H, Numerical Model Team, Korea Institute of Atmospheric Prediction Systems, 4F, Hankuk Computer Building, 35 Boramae-Ro 5 GIL, Dongjak-Gu, Seoul 07071, South Korea, Tel: 02-6959-1600, Email: [email protected]

Received Date: Aug 16, 2017 / Accepted Date: Oct 13, 2017 / Published Date: Oct 27, 2017


A semi-implicit time integration scheme is implemented in a non-hydrostatic Euler problem on the cubed-sphere grid. The semi-implicit time integration scheme is a 3-stage additive RungeKutta method, which is an Implicit-Explicit (IMEX) multi-stage single-step scheme. The system treats acoustic and gravity waves implicitly and advection explicitly. In the implicit part, we compute the linear system by defining a proper linear operator. Then, the numerical results are presented to compare the semi-implicit time integration scheme with the explicit RungeKutta time integration scheme in non-hydrostatic Euler equations. In terms of accuracy, we demonstrate that the proposed method performs better than other schemes.

Keywords: Semi-implicit time integration scheme; Additive Runge- Kutta method (ARK); Explicit Runge-Kutta method; Cubed-sphere grid; Non-hydrostatic Euler equations


In general, time integration schemes may be categorized into two classes: one is an explicit time integration scheme, and the other is an implicit time integration scheme [1-4]. The advantages of the former class of schemes are simplicity and high parallel efficiency with minimal inter processor communications [5-7]. Giraldo [8] showed that the spectral element can linearly scale with an explicit leapfrog time integration scheme in an atmospheric model. The explicit time integration scheme is simple and efficient in paralleling. However, in a general linear system, the fastest waves are acoustic waves, which have little to no effect on large-scale processes [9-11]. In other words, acoustic waves are very fast but have little significance in terms of the accuracy of simulations. Explicit systems have to adhere to a very small time step restriction caused by a physical phenomenon that is essentially inconsequential. To overcome these issues, a time-split integration scheme was proposed [12,13-20]. In this scheme, fast acoustic waves use a smaller time step while slower waves use a larger time step [21-25]. The scheme is based on the Runge-Kutta methods. Choi and Hong [3] explored the time-split third-order Runge-Kutta explicit method described [12,18,19,25]. The scheme is associated with a mass vertical coordinate and a flux-form set of non-hydrostatic Euler equations. Additionally, each stage of the Runge-Kutta method is divided into a number of sub steps in which the time tendencies are updated using fast acoustic and gravity wave terms in the governing equations. The sub steps use a forward-backward time integration scheme in which the vertical coupling terms are treated implicitly. These are described in detail [3] and in the references therein. For the latter class of schemes, an implicit scheme requires expensive numerical solvers, because the non-linear system must be solved by computing the Jacobian matrix. Nevertheless, the time step used in implicit schemes is greater than that used in explicit schemes. St-Cyr and Neckels [21], the authors discuss fully implicit time integration schemes for discontinuous Galerkin methods applied to atmospheric flow. In particular, they proposed Jacobian-free Newton-Krylov solvers without the computation of the Jacobian matrix. Evans et al. [7] studied the fully implicit time integration scheme, which was implemented in shallow water equations on a sphere with a spectral finite element. From the numerical results, they determined the time step that would minimize the total computing time while maintaining sufficient accuracy for these problems.

Recently, more attention has been focused on the way in which implicit and explicit time integration schemes are combined to form semi-implicit schemes in the global atmosphere model [4-24]. Semiimplicit time integration schemes treat acoustic and gravity waves implicitly. Advection is usually semi-Lagrangian [5] but is sometimes explicit Eulerian [22]. In either case, the time step is not limited by the speed of sound, gravity wave speed or stratification [4]. With appropriate linearization, spatial discretization leads to a sparse, diagonally dominant matrix that can be solved using iterative solvers. Semi-implicit time integration schemes can therefore be very efficient due to their ability to use long time steps. Moreover, linearly implicit schemes, such as an additive Runge-Kutta scheme, are applied in nonhydrostatic models [9,10]. A new second-order additive Runge-Kutta scheme was introduced [9]. The idea behind the new additive Runge- Kutta method is the use of two different integrators for non-stiff and stiff terms, respectively. An implicit integrator will be used for the stiff portion that represents acoustic and gravity waves, whereas an explicit integrator will be used for the non-stiff portion that represents the advective terms.

Now, we are interested in the semi-implicit time integration scheme for a non-hydrostatic problem. In this paper, we attempt to implement the new second-order 3-stage additive Runge-Kutta scheme in the non-hydrostatic Euler problem. Various factors were considered in the selection of the semi-implicit scheme. It is important desirable to keep the spatial discretization unchanged and to retain the time integration scheme [3] to facilitate a clean comparison between the semi-implicit and explicit schemes.

This paper is organized as follows. In Section 2, we consider the non-hydrostatic Euler problem used as the governing equations. The mass vertical coordinate and a flux-form set of equations are presented in this section. In Section 3, we review the s-stage additive Runge-Kutta schemes. These are the semi-implicit time integration schemes. Then, we define the proper linear operator. Based on these equations, we apply the semi-implicit time integration scheme to the problem. Finally, some numerical experiments for the idealization test are presented in Section 4. The performances of the semi-implicit and explicit schemes are compared in terms of error estimations where the semi-implicit scheme has some advantages over the explicit scheme. This paper ends with the summary and concluding remarks.s

Governing Equations

In this section, we consider the non-hydrostatic Euler equations with a flux form [3,18]. The governing equations are first formulated using a terrain-following mass vertical coordinate [15]. In η-coordinate, we define


where μd is the mass of the dry air in the column. As usual, let pdh be the hydrostatic pressure of the dry atmosphere and let pdht be defined pdh at the top. The coupled variables are denoted as


where v=[u,v]T represents the velocities in the horizontal (zonal and meridional) directions and w=η/dt represents the velocity in the vertical (radial) direction. θ is the potential temperature.

Suppose that the reference state is in hydrostatic balance and is strictly only a function of height Equation. The reference state variables are denoted by:

Equation (1.1a)

Equation (1.1b)

Equation (1.1c)

Equation (1.1d)

Because the η-coordinate surfaces are not horizontal, the reference profiles Equation are functions of variables (x,y,η). New variables Equation from the perturbation forms eqn. (1.1) are obtained and Equation is defined with θ0=300(K).

By using the perturbation variables Equation the governing equations are given as follows:

Equation (1.2a)

Equation (1.2b)

Equation (1.2c)

Equation (1.2d)

Equation (1.2e)

Where Fv and FW represent the forces of Carioles and curvature, respectively. Equation is the absolute vertical vorticity with f=2Ωe Sin ψ and e=2Ωe cosψ. Here, Ωe is the angular rotation rate of the earth and ψ is the latitude. Let us denote by Equation the horizontal kinetic energy. In these equations, αd and α are the inverse density of the dry air and the full air including moisture (i.e., α=αd (1+qvapor+qcloud+qrain+···)-1, where q are the mixing ratios for hydrometeors), respectively. These are described in detail [3] and the values can be found therein. In the governing equations (1.2), differential operators are defined by


Where V=[U,V]T is the vector value and ϕ represents a generic variable.

Equation (1.3)

Also, the full pressure is

Equation (1.4)

Where γ=cp/cv: =1.4 is the ratio of the heat capacity for dry air, Rd is the gas constant for dry air, p0=105 Pa is a reference pressure and Equation.

The Euler problem in eqn. (1.2) employs the spectral element method for the horizontal discretization and the finite difference method for vertical discretization [3]. In addition, the time-split thirdorder Runge-Kutta explicit method [12,17-20,25], which we call SRK3 in this paper, is used as the time integration scheme. The method shows only second order accuracy for nonlinear equations. In the time integration scheme, we consider two kinds of time steps. Slow or lowfrequency modes are integrated using a third-order Runge-Kutta time integration scheme, as follows:


With larger time step Δt. The values Φn and Φn+1 are represented at time level n Δt and (n+1)Δt, respectively. The high-frequency acoustic modes are integrated over a smaller time step Δτ to maintain numerical stability. The horizontally propagating acoustic modes and gravity waves are integrated using a forward-backward time integration scheme and the vertically propagating acoustic modes and buoyancy oscillations are integrated using a vertically implicit scheme. The procedures is described in detail [3,18]. From the properties, SRK3 uses a larger time step than the original explicit Runge-Kutta scheme. Models such as the WRF model [18], MPAS [20], and NICAM [17] have used the same time integration scheme and have been shown to work effectively with Δt as the time step for low-frequency modes (the model time step). High frequency but meteorologically insignificant acoustic modes would severely limit the time step Δt of the third order Runge-Kutta scheme in the governing equations.

We consider the non-hydrostatic Euler problem in eqn. (1.2) as the governing equations and keep the spatial discretization, which is defined [3], unchanged. We simply replace SRK3 with a semi-implicit time integration scheme in the system. In the next section, we will explain the scheme.

Semi-implicit Time Integration Schemes

In this section, we will address the use of the additive Runge- Kutta method as the semi-implicit multi-stage scheme. The governing equations can be rewritten in the compact vector form as follows:

Equation (2.1)

where the vector value Equation and the right-hand side S(x) represents the remaining terms in the equations apart from the time derivatives in eqn. (1.2).

S-stage Additive Runge-Kutta (ARK) method

In order to construct the semi-implicit time integration scheme in eqn. (2.1), we have to define a linear operator L(x), which approximates S(x) and contains the terms responsible for the acoustic and gravity waves. Then, the equation (2.1) is rewritten as:

Equation (2.2)

and explicitly discretizes in time the terms in [S(x)-L(x)] while implicitly discretizing those in L(x).

As was done [1,9,11,16], we now consider the discretization of s-stage ARK method in eqn.(2.2) as the following:




For n=0,1,··· ,T, the vector xn+1 and xn represent the numerical solution at time level (n+1)Δt and nΔt (:=tn), respectively, where Δt is the time step and T is the final time. For each stage, we obtain the value Xi with the coefficients Equation for i, j=1,··· ,s, from the Butcher tableaux [2,9].

In the Butcher tableaux, the relation is as follows:


Where Equation represent the time when [S(x)-L(x)] and L(x) are evaluated, respectively. It means that at each stage the values are evaluated at Equation.

From the s-stage ARK method, we have to solve a few implicit parts for each step, which is equal to the cardinality of Equation. Additionally, the implicit part of s-stage ARK schemes can achieve Aand L-stability properties of an arbitrarily high order [9].

In fact, the scheme is intended to use explicit and implicit integrators for the non-stiff and stiff terms, respectively. The stiff part represents the acoustic and gravity waves, whereas the non-stiff part represents the adventive terms. From this idea, we will introduce a new second-order Runge-Kutta method as a semi-implicit time integration scheme that performs better than the other second order methods considered.

New 3-stage (ARK) method with linear operator

Giraldo [9] developed the original and new versions of the 3-stage ARK method. These are second-order method, and we call the new version of the 3-stage ARK method ARK2 in this paper. The main idea behind ARK2 is Equation. It means that L-stable, with minimal cost per step, has at least 3-stage with the first stage being explicit. The number of computing implicit part in ARK2 is smaller than in ARK. The majority of the computing time in this scheme is occupied by solving the implicit part in numerical experiments. These aspects are different from the original ARK method.

By applying the order conditions and stability constraints, the coefficients of ARK2 are obtained from the Butcher tableaux [2,9] as follows:


Since the diagonal components of A are zero, we computed explicitly in time the terms [S(x)−L(x)] in eqn. (2.3). For


has been computed only once in eqn. (2.3), where


To proceed, the specific times are Equation

In the main idea of the semi-implicit scheme, the linear operator part is computed implicitly and depends on the problem. If the correct operator L is not obtained, the semi-implicit method will not work because the operator L must be selected such that the fastest waves in the system are retained, albeit in their linearized forms. Deriving the linear operator L(·) in eqn. (2.2) is now straightforward. In this paper, we first implemented ARK2 in non-hydrostatic Euler problem. The system is based on a mass vertical coordinate.

Here, we assume that α=αd with dry air. From the governing equations, we define the linear operator L(·) as follows


where the new linearization forms of perturbation pressure Equation and the inverse of density Equation in eqn. (1.3-1.4) are given as follows:


Note that the full pressure in eqn. (1.4) cannot be written in perturbation form because of the exponent function in the expression. For small perturbation simulations, accuracy for perturbation variables can be maintained by linearizing the perturbation variables. The hydrostatic relation in the perturbation system becomes:


Recall that there are no linear parts for the potential temperature θ in eqn. (1.2).

High order schemes, specifically the third (4-stage) order ARK3, fourth (6-stage) order ARK4, and fifth (8-stage) order ARK5, were developed [11]. The schemes presented above are all explicit firststage, singly diagonal, second-stage order, L-stable methods. In our numerical results, we only attempt to apply the third-order ARK3.

Numerical Results

In this section, we present a numerical example of the Rossby- Haurwitz wave test for the no hydrostatic Euler problem on a cubedsphere grid. We use the ne30np4 resolution (approximately 110 km) in the experiments. For horizontal grid spacing, ne indicates the number of quadrilateral elements in each direction for each face of the cube, and np is the number of Gauss-LegendreLobatto quadrature points. Throughout the numerical study in this section, numerical simulations are run for 3600(s) with a time step of Δt.

We will compare the numerical results of the semi-implicit time integration scheme ARK2 and the time-split explicit integration scheme SRK3. As the reference values of the test, we use the third-order explicit Runge-Kutta method, which will be abbreviated as RK3, with a very small time step Δt=0.01(s). In general, RK3 is reasonably simple and robust. In these experiments of RK3, time steps smaller than 0.1(s) are used. If the time step is greater than 0.1(s), it is blow up. To apply the semi-implicit scheme, we have to solve the linear system using a numerical iterative solver. In this case, the generalized minimal residual method (GMRES) is applied without a preconditioned. The number of iterations in the numerical solver depends on the size of the matrix and its tolerance (tol). In this paper, we fixed tol=10−4. If we want to obtain a more accurate numerical solution by using the iterative solver, we can select a smaller tolerance than 10−4.

The numerical results in ARK3, which are zonal wind u, meridional wind v, surface pressure ps and temperature θ at 850 hPa with Δt=2.0(s), are presented in Figure 1. They are similar to the numerical results for RK3 (at Δt=0.1(s)) and SRK3 (at Δt=2.0(s)).


Figure 1: The numerical results u(top, left), v(top, right), ps(bottom, left) and θ(bottom, right) of ARK2 at 850 hPa are plotted with a final time T=3600(s) and Δt=2.0(s).

The differences between the numerical results u, v, ps and θ for SRK3 and those for ARK2 at 850 hPa with Δt=2.0(s) are presented in Figure 2. In figures, the L2-errors of different values for each variable are bounded by 10-5. SRK3 has faster computing times than ARK2 at the same time step because the linear system must be solved using the numerical iterative method in the semi-implicit scheme.


Figure 2: The differences in the numerical results u(top, left), v(top, right), ps(bottom, left) and θ(bottom, right) for SRK3 and ARK2 at 850 hPa are plotted with a final time T=3600(s) and Δt=2.0(s).

In this paper, the semi-implicit scheme ARK2 is compared with the time-split Runge-Kutta explicit scheme SRK3 and we emphasize the accuracy of the semi-implicit scheme. We will show that ARK2 is more accurate than SRK3 for L2-error estimation. In L2-error estimation, we compute only the L2-error of the surface pressure ps in ARK2 and SRK3 for each time step, as in the following equation:


where (ps)RK3, as the reference numerical solution, is the value of the surface pressure from RK3 at Δt=0.01(s). (ps)h, with h=SRK3 or ARK2, is the numerical solution of the surface pressure for each time step Δt=0.2, 0.4, ··· , 2.0(s).

Since SRK3 is the time-split scheme, each stage of the SRK3 is divided by a number of sub steps ns=2,4,8,16. ns is the ratio of the RK3 time step to the acoustic time step for the second and third RK3 sub steps. Figure 3 shows the L2-error of the surface pressure for each sub step ns. L2-errors for the surface pressure are decreased when ns increases. In addition, since ARK2’s values are bounded by 10−8 and SRK3’s values are bounded by 10−4, we cannot compare the values directly in Figure 3. Thus, the L2-error of surface pressure of only ARK2 for each time step is presented in the right-hand side of Figure 4. From these results, for the same time step, we know that the L2-error of ARK2 is much smaller than that of SRK3.


Figure 3: L2-error analysis of the surface pressure for each time step in SRK3 and ARK2. The L2-error estimations of SRK3 with the number of substeps ns=2,4,8,16, and ARK2 are plotted. The reference solution is the value of the surface pressure determined using RK3 at Δt=0.01(s).


Figure 4: L2-error analysis of the surface pressure for each time step in ARK2. The L2-error estimations of only ARK2 are plotted. The reference solution is the value of the surface pressure determined using RK3 at Δt=0.01(s).

We note that construction of the high-order ARK methods is complex and difficult. The third order ARK method with 4-stage, which we call ARK3, is only computed in this experiments. From Table 1, we observe that L2-errors for ARK2 and ARK3 are smaller than those for SRK3 and the values of ARK3 seem to be the smallest. However, since ARK3 has an additional stage, the total computing the time has to be increased. Thus, we know that ARK2 and ARK3 are more accurate than time split explicit scheme SRK3, while ARK2 and ARK3 are much slower than SRK3 in terms of computing time. If we consider the proper pre-conditioner in the iterative solver or transformation of the system with variable dependency, we can obtain a faster scheme. These are the topics of on-going research.

Ps 4.81e-06 2.26e-09 2.31e-10
u 2.15e-06 8.90e-07 2.01e-07
v 1.85e-06 9.56e-07 2.82e-07
q 6.95e-07 1.18e-10 8.33e-11

Table 1: At Δt =1.0(s), L2-error analysis of numerical results ps, u, v and θin SRK3, ARK2, and ARK3 as the high order semi-implicit scheme. The reference value is obtained using RK3 at Δt=0.01(s).


In this paper, we implement the 3-stage additive Runge-Kutta method, which is a second-order semi-implicit time integration scheme, on the non-hydrostatic Euler problem with mass vertical coordinates. In this scheme, we use a larger time step than that used in the third-order RungeKutta explicit scheme, as the computing time in the semi-implicit scheme is long. This is done with the aim of solving the linear system with a numerical iteration scheme. Moreover, the time step of the semi-implicit scheme is larger than that of the original explicit Runge-Kutta scheme, but smaller than that of the timesplit Runge-Kutta explicit scheme. However, the numerical results demonstrate that the s-stage additive Runge-Kutta semi-implicit time integration scheme is more accurate than other explicit schemes in the non-hydrostatic Euler problem.

If we use a higher order additive Runge-Kutta semi-implicit scheme in this problem, more accurate numerical solutions are obtained.


This work has been carried out through the R&D project on the development of global numerical weather prediction systems of Korea Institute of Atmospheric Prediction Systems (KIAPS) funded by Korea Meteorological Administration (KMA).


Citation: Nam H, Choi SJ (2017) Implementation of a Semi-implicit Time Integration Scheme in Non-Hydrostatic Euler Equations. J Appl Computat Math 6: 369. Doi: 10.4172/2168-9679.1000369

Copyright: © 2017 Nam H, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Select your language of interest to view the total content in your interested language

Post Your Comment Citation
Share This Article
Article Usage
  • Total views: 1345
  • [From(publication date): 0-2017 - Jun 20, 2018]
  • Breakdown by view type
  • HTML page views: 1269
  • PDF downloads: 76

Post your comment

captcha   Reload  Can't read the image? click here to refresh
Leave Your Message 24x7