Abstract

In this paper, theoretical and numerical analyses of the properties of some complex semi-Lagrangian methods are performed to deal with the issues of the instability associated with the treatment of the nonlinear part of the forcing term. A class of semi-Lagrangian semi-implicit schemes is proposed using a modified TR-BDF2 method, which is the combination of the trapezoidal rule (TR) and the second-order backward differentiation formula (BDF2). The process used for the nonlinear term includes two stages as predictor and corrector in the trapezoidal method and one stage for the BDF2 method. For the treatment of the linear term, the implicit trapezoidal method is employed in the first step, the explicit trapezoidal method in the second step, and the implicit BDF2 method in the third step. The combination of these techniques leads to a family of schemes that has a large region of absolute stability, performs well for the purely oscillatory cases, and has good qualities in terms of accuracy and convergence. The use of the explicit method for the linear term in the second step makes the proposed class of schemes competitive in terms of efficiency compared to some well-known schemes that use two steps. Numerical experiments presented herein confirm that the proposed class of schemes performs well in terms of stability, accuracy, convergence, and efficiency in comparison with other, previously known, semi-Lagrangian semi-implicit schemes and semi-implicit predictor–corrector methods. The potential practical application of the proposed class of schemes to a weather prediction model or any other atmospheric model is not discussed and could be the subject of other forthcoming studies.

1. Introduction

Semi-Lagrangian (SL) semi-implicit integration methods, first proposed by Robert (1981), have been extensively studied and widely incorporated into atmospheric numerical models. A comprehensive review of the applications of SL methods in atmospheric problems was given by Staniforth and Côté (1991). Other extensive studies have also been conducted to examine SL methods (e.g., Bonaventura 2000; White and Dongarra 2011). A three-time-level SL method was implemented operationally at the European Centre for Medium-Range Weather Forecasts in 1991, and was documented and described in Ritchie et al. (1995). Tanguay et al. (1990) generalized the use of a semi-implicit algorithm in order to integrate the fully compressible nonhydrostatic equations. Ritchie (1991) and Ritchie and Beaudoin (1994) applied the semi-Lagrangian semi-implicit method to a multilevel spectral primitive equation model.

McDonald and Bates (1987) and Temperton and Staniforth (1987) proposed two-time-level SL methods. This motivated many authors to use two-time-level SL methods instead of three-time-level schemes, since they reduce the number and impact of computational modes (e.g., McDonald and Haugen 1992; Temperton et al. 2001; Hortal 2002). The use of more than two time levels leads in general to computational modes having comparable amplitudes to those of physical modes, which can influence the accuracy of the solution if no efficient technique is used to damp those modes.

Hortal (2002) proposed a two-time-level SL method called the Stable Extrapolation Two-Time-Level Scheme (SETTLS), which has a region with absolute stability independently of the Courant number. Gospodinov et al. (2001) proposed a family of second-order two-time-level SL schemes that contain an undetermined parameter α. Durran and Reinecke (2004) showed that the size of the absolute stability region of the family of second-order schemes proposed by Gospodinov et al. (2001) varies dramatically according to the undetermined parameter, and that the optimal region is obtained for α = ¼, which corresponds to SETTLS. Still, SETTLS is not a perfect choice because the points on the imaginary axis are outside the domain of absolute stability, except for the origin point. Therefore, SETTLS can generate growth parasites for a purely oscillatory nonlinear term depending on the size of the time step.

The combination of the trapezoidal rule (TR) and the second-order backward differentiation formula (BDF2) is more commonly known as the TR-BDF2 method. Based on the analysis done by Dharmaraja (2007), the TR-BDF2 method performs well in terms of stability for solving partial differential equations. The method is able to compute the solutions using a reasonable step size for stiff problems. This motivated us to use this method to develop the semi-Lagrangian schemes.

In this paper, theoretical and numerical analyses are performed to study the properties of more complex methods. The aim is to try to avoid the problems of instability associated with the treatment of the nonlinear part of the forcing term. We propose a class of semi-Lagrangian semi-implicit schemes using a modified TR-BDF2 method. We use two stages as predictor and corrector in the trapezoidal method and one stage for the BDF2 method. The family of second-order approximations proposed in Gospodinov et al. (2001) is used in the predictor of the TR method. For the linear term we use the implicit trapezoidal method in the first step, the explicit trapezoidal method in the second step, and the implicit BDF2 method in the third step. Following Hortal (2002), the equation of the semi-Lagrangian trajectory used in the SETTLS method is obtained using an average approximation of the acceleration between the departure point and the arrival point. Since the middle point is in the interval where the average approximation of the acceleration is considered, in the proposed method we use the same average of the acceleration to obtain the position of the middle point. An explicit equation is used for this point and the only iterative equation is the one used to obtain the departure point. The potential practical application of the proposed schemes to a weather prediction model or any other atmospheric model is not analyzed in the current paper. The analysis of the application of the proposed class of schemes to these models will be considered in future studies.

The remainder of the paper is organized as follows. Section 2 introduces the notation and explicitly states the problem and expectations of the paper. Section 3 gives details on the proposed schemes for the nonlinear term as well as their linear stability. In section 4, the proposed schemes for the linear term are presented as well as their stability analysis. In section 5 the full semi-Lagrangian semi-implicit schemes are presented, and their accuracy, efficiency, and convergence are analyzed and compared to other existing schemes. Finally, section 6 provides some concluding remarks.

2. Problem and notation

The forcing source term of the general equation within the framework of the SL method can be split into a fast linear term and a slower, residual nonlinear term as follows:

 
formula

We assume that the two terms depend on the field variable ψ(x, t), x, and time. At each time tn = nΔt, where Δt is the time step, the values of the parameters ψ and the velocity V are known for all times tp = pΔt with p = 0, 1, 2, …, n. The proposed method requires for each arrival point xj the evaluation of the position at time tn of the air parcel arriving at a point (xj, tn+1) and the position of the middle point on the semi-Lagrangian trajectory. The details of the computation of the departure point and the position of the middle point along the semi-Lagrangian trajectory are presented in  appendix B.

In general, for solving (1), explicit methods are used for the nonlinear term N, and the linear term L is calculated using implicit methods. The superscript plus signs (+), zeros (0), and minus signs (−) are used to denote the variables at time levels t + Δt, t, and t − Δt, respectively. The subscripts A and D are used to denote the parameters at the departure point xD = x(t) and the arrival point xA = x(t + Δt), respectively. For the discrete form of the schemes, we define as an estimate of the departure point of the fluid parcel at time tn that arrives at (xj, tn+1). To abbreviate the notation, for any step p, we use and as the values of the nonlinear term and the linear term, respectively, at the departure point of the trajectory originating at and arriving at (xj, tn+1). Following theorem 1 in Gospodinov et al. (2001), the following discrete approximations of (1),

 
formula

lead to second-order accuracy for any undetermined parameter α. We use the terminology α approximation for these second-order approximations of the temporal derivative of the function ψ along the interval of endpoints A and D of the SL trajectory. Durran and Reinecke (2004) presented a stability analysis of a class of explicit semi-Lagrangian schemes using (2), taking into consideration only the nonlinear term. They showed that the size of the region of absolute stability of the second-order schemes derived from (2) varies dramatically according to parameter α. SETTLS corresponds to the optimal region obtained by setting α = ¼. Following the stability analysis of this scheme using von Neumann’s method, the imaginary points are outside the domain of absolute stability, except for the origin. In this study, we propose a family of semi-implicit semi-Lagrangian schemes, which does not suffer from this disadvantage. We will use, as in the case of the SETTLS method, only one departure point, and we will examine several schemes by using the α approximation for the treatment of the nonlinear term. Three steps are used, where the nonlinear term is treated explicitly in all steps, and the linear term is treated implicitly in the first and third steps and explicitly in the second step. For simplicity of terminology, in the second step we will use the term trapezoidal method even if a decentering is applied. The proposed techniques for linear and nonlinear terms are combined to obtain the proposed class of semi-Lagrangian semi-implicit schemes, which provides interesting stability properties. Other advantages of this class of schemes will be demonstrated, such as the accuracy, convergence, and efficiency, which are compared to those obtained by using other existing semi-implicit predictor–corrector schemes and semi-Lagrangian semi-implicit methods.

3. Explicit scheme for the nonlinear term

a. The proposed explicit predictor corrector scheme

For the nonlinear term, we propose an explicit semi-Lagrangian scheme based on a modified TR-BDF2 method combining the trapezoidal rule and the second-order backward differentiation formula. Some details about the initial BDF2 method are given in  appendix A. In the modified TR-BDF2 method, two stages are used as predictor and corrector in the trapezoidal method, and one stage for the BDF2 method:

 
formula

The terms in the right-hand side of each stage correspond to the explicit form for the nonlinear operator. In the first stage, we calculate the predictor value ψ* by approximating the value of the nonlinear term Nn+1/4. This value is obtained by using a trapezoidal treatment of Nn and the α approximation of the nonlinear term denoted by at the midpoint of the SL trajectory originating at and arriving at (xj, tn+1). In the second stage, the corrector value of ψn+1/2 is calculated by using the trapezoidal method. We compute the value of N* by using ψ* obtained from the first stage, and apply a trapezoidal method using the two values Nn and N* to obtain the source term for this corrector stage. A large domain of absolute stability of the scheme is obtained and, as will be demonstrated below, the decentering parameter θ has significant ability to improve the stability of this scheme for purely oscillatory solutions. The range values of θ will be chosen in order to ensure good accuracy of the proposed scheme.

In the third stage, the BDF2 method is applied to obtain the next value ψn+1 by using an explicit approximation of the nonlinear term Nn+1 on the right-hand side of the third equation in (3), where Nn+1/2 is computed by using the value of ψn+1/2 obtained from the second stage. The value used in the predictor formula is estimated by using the α approximation as

 
formula

where, for example, N(x, tn):= N[ψ(x, tn), x, tn]

A family of SL schemes is obtained based on the values of the undetermined parameter α and the decentering parameter θ. As will be confirmed in the next section, the use of the corrector in the second stage improves the stability of the proposed schemes. We will demonstrate that α = ¼ is the optimal choice with the decentering θ in the interval [0.5, 0.75] to provide good accuracy.

b. Stability analysis of the proposed scheme

Following Durran (2010), the stability of the SL scheme defined by the undetermined parameter α and the decentering parameter θ will be analyzed by using the forced one-dimensional SL equation

 
formula

where σ = λ + , λ and ω are real constants. The stability will be analyzed using a constant advecting velocity U. Equation (5) becomes

 
formula

For any initial condition ψ(x, 0) = f(x), the analytical solution is

 
formula

We assume a nonamplifying solution by imposing the condition λ ≤ 0. The stability analysis of the SL schemes is examined in terms of the parameter α used in the α approximation of the predictor and the decentering parameter θ used in the second step. The analysis is performed using von Neumann’s method, which is applied for a single Fourier mode of the form for any point (jΔx, nΔt). When this mode is substituted into (5), letting z = λΔt + Δt, one obtains after simplification the following equation for the amplification factor Ak:

 
formula

where , , and s = xAxD.

The region of stability depends on the Courant number, which is proportional to the parameter s. The region on the λΔtωΔt plane where the scheme, for the case θ = 0.7, is absolutely stable is plotted in Fig. 1 for several values of the parameter of approximation α. Since the solutions of (8) are periodic in ks, the curves plotted in Fig. 1 are given for several discrete values of ks covering the whole periodic domain [0, 2π]. The regions of absolute stability are only slightly sensitive to the values of α, but the part of the imaginary axis inside the region of absolute stability is variable, and the most interesting is for α = ¼, which corresponds to ωΔt ∈ [−0.708, 0.708]. This has a positive influence for stability and, especially, for solutions that have a purely oscillatory character. The decentering acts directly on the size of the stable imaginary part, as shown in Fig. 2. Table 1 shows the values of (ωΔt)max of the segments [−(ωΔt)max, (ωΔt)max] on an imaginary axis, which are included in the regions of absolute stability depending on the values of the approximation parameter α and the decentering parameter θ. To obtain the proposed class of SL schemes, we set α = ¼, which improves stability. As shown in Table 1, the use of this value ensures a large stability zone along the imaginary axis for all non-high-decentering parameter θ.

Fig. 1.

Region of absolute stability for six values of the parameter α plotted for 15 Courant–Friedrichs–Lewy (CFL) number values. The central white region corresponds to the absolutely stable region independent of the CFL number for a decentering parameter θ = 0.7.

Fig. 1.

Region of absolute stability for six values of the parameter α plotted for 15 Courant–Friedrichs–Lewy (CFL) number values. The central white region corresponds to the absolutely stable region independent of the CFL number for a decentering parameter θ = 0.7.

Fig. 2.

Region of absolute stability of the proposed schemes (α = ¼) for four values of the decentering parameter θ plotted for 15 CFL number values. The central white region corresponds to the absolutely stable region independent of the CFL number.

Fig. 2.

Region of absolute stability of the proposed schemes (α = ¼) for four values of the decentering parameter θ plotted for 15 CFL number values. The central white region corresponds to the absolutely stable region independent of the CFL number.

Table 1.

Value of (ωΔt)max of the segment along an imaginary axis [−(ωΔt)max, (ωΔt)max] included in the region of absolute stability depending on the values of the approximation parameter α and the decentering parameter θ.

Value of (ωΔt)max of the segment along an imaginary axis [−(ωΔt)max, (ωΔt)max] included in the region of absolute stability depending on the values of the approximation parameter α and the decentering parameter θ.
Value of (ωΔt)max of the segment along an imaginary axis [−(ωΔt)max, (ωΔt)max] included in the region of absolute stability depending on the values of the approximation parameter α and the decentering parameter θ.

To conduct a comparison and in order to show the usefulness of using the trapezoidal corrector iteration, the regions of absolute stability are analyzed in the same way as the proposed class of schemes, but without using any corrector value for the variable ψn+1/2. We will use the two-step method and, first, we examine the case in which the value is considered in the BDF2 iteration to be an estimation of Nn+1/2. The same form as (8) is obtained using the parameters γ1 = −1 and γ2 = −z. For this case, as can be seen in Fig. 3, the results are similar to those obtained by Durran and Reinecke (2004) for the second-order two-time-level semi-Lagrangian schemes developed by Gospodinov et al. (2001) and in the particular case for the SETTLS method. The region of absolute stability is very sensitive and varies excessively depending on the value of parameter α. For the optimal region of absolute stability, the origin point is the only imaginary part.

Fig. 3.

Region of absolute stability of the scheme without a corrector step and using the value in the step employing the BDF2 method, for six values of the parameter α plotted for 15 CFL number values.

Fig. 3.

Region of absolute stability of the scheme without a corrector step and using the value in the step employing the BDF2 method, for six values of the parameter α plotted for 15 CFL number values.

In the case of the two-step method in which we consider the value of Nn+1/2 by using ψn+1/2 obtained from the first stage, a similar form of (8) is obtained with and . For this case, the region of absolute stability is still largely reduced compared to the region of absolute stability of the proposed class of schemes. This two-step method also presents the disadvantage of being stable on the imaginary axis of the λΔtωΔt plane at only a single origin point.

c. Accuracy of oscillatory solutions

To learn some basic information about the accuracy of the proposed class of schemes for oscillatory solutions, the following equation is analyzed in this section:

 
formula

where ω is a real number. In our analysis we will use the decentering parameter θ = 0.7 in the second step of (3). Note that the analytical solution of oscillatory equation (9) corresponding to one mode can be obtained by setting f(x) = eikx and σ = in (7), which leads to an analytical amplification factor of Aan = ei(ωΔtks). Therefore, for an accurate scheme, the module of the numerical amplification factor Ak and the relative phase change defined as should be close to unity. The polynomial equation for the numerical amplification factor can be obtained by setting z = Δt in (8). We use the ωΔtks plane to plot the parameters studied, where ωΔt ∈ [0, 0.7] and ks ∈ [0, 2π]. In Fig. 4 we plot the module of the errors |AkAan| of the amplification factor, the damping, and the relative phase change. The errors are negligible, as shown in Fig. 4a, which is also confirmed in Figs. 4b,c, where we observe negligible damping and phase errors.

Fig. 4.

Solutions to the oscillatory equation using the proposed schemes for the nonlinear term. (a) Module of the errors of the amplification factors |AkAan|. (b) Damping of the solution. (c) Relative phase change. (d) Module of the computational mode.

Fig. 4.

Solutions to the oscillatory equation using the proposed schemes for the nonlinear term. (a) Module of the errors of the amplification factors |AkAan|. (b) Damping of the solution. (c) Relative phase change. (d) Module of the computational mode.

The proposed class of schemes presents the advantage of having a damped computational mode as shown in Fig. 4d compared to the computational modes of the class of semi-implicit predictor–corrector schemes developed by Clancy and Pudykiewicz (2013). Note that more comparisons in terms of accuracy between those schemes and the proposed full semi-implicit semi-Lagrangian schemes will be presented in detail in section 5.

4. The proposed scheme for the linear term

a. Treatment of the linear term

We seek a numerical scheme for the linear term that is stable and accurate, in particular for the case of oscillatory solutions. This scheme should have some compatibility with the proposed scheme for the nonlinear term. Since three steps are used for the nonlinear term, in order to reduce the computational cost of the full scheme, we try to use two implicit steps and one explicit step for the linear term. We will use the implicit trapezoidal method in the first step, the explicit trapezoidal method in the second step, and the implicit BDF2 method in the third step:

 
formula

The scheme for the linear term, which is the subject of this section, and the scheme for the nonlinear term given in the previous section are considered separately. It is well known that when the two schemes are combined, they interact in a complex way. Following our numerical tests, the use of the same decentering parameter θ in the second step for the schemes used for linear and nonlinear terms is required to avoid a large impact on the accuracy. On the right-hand side of the third equation in (10), a family of second-order approximations of the linear term Ln+1 at time tn+1 is considered and the choice of the appropriate value of parameter μ is required.

b. Stability analysis and choice of parameter μ

First, we will consider the case without decentering (θ = 0.5) to determine the appropriate value of parameter μ for both stability and accuracy. In the same way as in the previous section, the stability analysis of the SL schemes in (10) is examined using von Neumann’s method, and a single Fourier mode of the form for any point (jΔx, nΔt) is considered. Substituting this mode into (5) and letting z = λΔt + Δt, one obtains the following equation for the numerical amplification factor Ak:

 
formula

For the proposed schemes, we will try to choose the value of parameter μ that guarantees stability and accuracy, in particular for the purely oscillatory cases. Note that in general the amplification factor of the form

 
formula

which has the module

 
formula

where c is a positive real number, leads to stable results (λ ≤ 0) and the scheme preserves perfectly the purely oscillatory solutions.

In our case we try to have the amplification factor Ak, obtained from (11), with the same form as (12). Therefore, the two second-degree polynomials in the numerator and denominator of (11) should have one common root and one root of the opposite sign. This condition leads to μ = 0.5 and we obtain

 
formula

which leads to a stable scheme that perfectly preserves the oscillatory solutions.

In Fig. 5a we plot the module of the amplification factor, which is obtained by using (11) for the oscillatory equation (9) and without decentering (θ = 0.5), as a function of the parameters μ and ωΔt. In this case the schemes are stable for μ ≤ 0.5 and we do not observe any damping in the case of oscillatory solutions for μ = 0.5. For general equation (5), in Fig. 5b we plot for the case θ = 0.5 the module of the amplification factor |Ak| of the solutions for the proposed schemes (μ = 0.5). The schemes are stable and the oscillatory solutions are perfectly preserved. The analytical solutions that correspond to |Aan| = eλΔt are represented in the same figure by the vertical lines.

Fig. 5.

Treatment of the linear term. (a) Module of the amplification factor, for (9) using θ = 0.5, plotted as function of μ and ωΔt. (b) Module of the amplification factor for (5) using μ = 0.5 and θ = 0.5 plotted along the λΔtωΔt plane.

Fig. 5.

Treatment of the linear term. (a) Module of the amplification factor, for (9) using θ = 0.5, plotted as function of μ and ωΔt. (b) Module of the amplification factor for (5) using μ = 0.5 and θ = 0.5 plotted along the λΔtωΔt plane.

Note that in the case in which decentering is considered, the value μ = 0.5 leads to an amplification factor Ak closer to (12) with the remainder of order |z|2, where |z| ≪ 1. For this case, if we consider the purely oscillatory solutions, the module of the amplification factor can be approximated for z = Δt and under the assumption |z| ≪ 1 by

 
formula

In Fig. 6a, we plot the module of the amplification factor for (5) for θ = 0.55, and the results confirm that the scheme is stable. In the same way, if we use θ = 0.45 for (5), the results are not stable, as shown in Fig. 6b. Finally, we use μ = 0.5 for the full semi-implicit semi-Lagrangian combination, and θ is the only variable parameter, with 0.5 ≤ θ ≤ 0.75.

Fig. 6.

Treatment of the linear term. Module of the amplification factor for (5), plotted along the λΔtωΔt plane. (a) Example of a stable case using μ = 0.5 and θ = 0.55. (b) Example of an unstable case using μ = 0.5 and θ = 0.45.

Fig. 6.

Treatment of the linear term. Module of the amplification factor for (5), plotted along the λΔtωΔt plane. (a) Example of a stable case using μ = 0.5 and θ = 0.55. (b) Example of an unstable case using μ = 0.5 and θ = 0.45.

5. Numerical experiments of the semi-Lagrangian semi-implicit combination

a. The full semi-Lagrangian semi-implicit scheme

Following the previous sections, we will consider the following family of semi-Lagrangian semi-implicit schemes, which are obtained by considering α = 0.25 and μ = 0.5:

 
formula

These schemes are examined, and we will demonstrate their advantages in terms of accuracy, efficiency, and convergence compared to other existing methods such as the semi-implicit predictor–corrector schemes proposed by Clancy and Pudykiewicz (2013) and the semi-Lagrangian semi-implicit scheme studied by Cullen (2001).

b. Linear stability and error analysis of the two-frequency system

We consider the following semi-Lagrangian equation:

 
formula

where Ω is the frequency of the fast waves and ω is the frequency of the slow waves, which is the remaining nonlinear term (|Ω| > |ω|). If we consider the conjugate of the two members in (17), we obtain a similar equation using the two frequencies −Ω and −ω, and the conjugate function is its solution. Therefore, in our analysis and without loss of generality, we will use the solutions of (17) under the assumption that Ω is positive and Ω > |ω|.

The analysis is performed using von Neumann’s method for a single Fourier mode of the form for any point (jΔx, nΔt). This mode is substituted into (17) to obtain the following equation for the amplification factor Ak:

 
formula

where

 
formula

and

 
formula

The complex numbers z1 and z2 are related to the fast and slow waves, respectively, and are given by z1 = iΩΔt and z2 = Δt.

In our analysis we will consider a parameter of decentering θ = 0.7 for the proposed schemes. In Fig. 7a we plot the module of the amplification factor Ak corresponding to the physical solution for five values of the parameter ks on the ωΔt − ΩΔt plane. The dashed lines are the limit defined by Ω = |ω| of the domain subject of our analysis. The results confirm the stability of the proposed schemes. In Figs. 7b, 7c, and 7d, we plot the module of the errors |AkAan| of the amplification factors, respectively, for the three values ks = 2π/3, ks = π, and ks = 4π/3. Figure 8 shows the module of the errors for the same system (17) for the four semi-implicit predictor–corrector schemes proposed in Clancy and Pudykiewicz (2013), which perform well for this test. The authors used the notations AM2*–LFT, AM2*–ABM, T–ABT, and AM2*–ABT for their schemes, where they used the leapfrog trapezoidal (LFT) method, the Adams–Bashforth trapezoidal (ABT) scheme, and the Adams–Bashforth–Moulton (ABM) method as explicit predictor–corrector schemes. For the implicit schemes, the authors used the trapezoidal method (T) and a method denoted by AM2*. The schemes proposed in this paper perform very well in terms of accuracy compared to the three schemes AM2*–LFT, AM2*–ABM, and AM2*–ABT, as shown in Figs. 7 and 8. Note that the best method in Clancy and Pudykiewicz (2013) for this test in terms of amplitude is the method T–ABT, but it presents some phase errors and it is less accurate than the proposed schemes, as shown in Figs. 7 and 8.

Fig. 7.

The proposed full semi-Lagrangian semi-implicit schemes using μ = 0.5 and θ = 0.7. (a) Module of the amplification factors for five CFL number values. (b) Module of the errors of the amplification factors |AkAan| for ks = 2π/3. (c) Module of the errors for ks = π. (d) Module of the errors for ks = 4π/3.

Fig. 7.

The proposed full semi-Lagrangian semi-implicit schemes using μ = 0.5 and θ = 0.7. (a) Module of the amplification factors for five CFL number values. (b) Module of the errors of the amplification factors |AkAan| for ks = 2π/3. (c) Module of the errors for ks = π. (d) Module of the errors for ks = 4π/3.

Fig. 8.

Module of the errors of the amplification factors |AkAan| for the schemes proposed by Clancy and Pudykiewicz (2013) for the (a) AM2–LFT, (b) AM2–ABM, (c) T–ABT, and (d) AM2–ABT methods.

Fig. 8.

Module of the errors of the amplification factors |AkAan| for the schemes proposed by Clancy and Pudykiewicz (2013) for the (a) AM2–LFT, (b) AM2–ABM, (c) T–ABT, and (d) AM2–ABT methods.

In the same way, we will compare the errors of the amplification factor for the semi-Lagrangian semi-implicit method studied by Cullen (2001) and those obtained by using the proposed schemes. The amplification factor of the scheme studied by Cullen (2001) applied to the same equation (17) is as follows:

 
formula

where, as previously defined, z1 = iΩΔt and z2 = Δt.

In Fig. 9a we plot the module of the errors |AkAan| of the amplification factor for the method studied by Cullen (2001). As can be seen in Figs. 9 and 7, for positive and large values of ωΔt the proposed schemes are more accurate than the method studied by Cullen (2001). Note that the errors of this method are mainly due to the phase errors, as shown in Fig. 9b.

Fig. 9.

The scheme studied by Cullen (2001). (a) Module of the errors of the amplification factors |AkAan|. (b) Relative phase change.

Fig. 9.

The scheme studied by Cullen (2001). (a) Module of the errors of the amplification factors |AkAan|. (b) Relative phase change.

c. Accuracy, efficiency, and convergence of the proposed class of schemes

In this section we will demonstrate the advantage of the proposed class of schemes in terms of accuracy, convergence, and efficiency. We will compare the results to those obtained by using the semi-Lagrangian semi-implicit scheme studied by Cullen (2001). We analyze the following system considered by Cullen (2001):

 
formula

where ϕ is the geopotential, D is the divergence term, and c is the gravity speed. We consider c0 as the reference value of c, which is used to split the source term c2D into two parts. The first part, , based on the reference value c0, is considered to be the linear term, which corresponds to the fast waves, and the second part is the residual term , which is considered to be the nonlinear term. Since this second term is responsible for the slow waves, we will consider the condition , which leads to δ ≥ 0.5, where δ = (c0/c)2. In the first equation of (19), the source term ∇2ϕ is considered to be the linear term. The error analysis and the numerical tests will be performed using one mode, as well as other solutions, which are obtained by combination of two modes with different phases.

1) Error analysis using one mode

In our analysis we consider a single wave of the form for any spatiotemporal point (jΔx, nΔt), where k is the wavenumber. We use the previous abbreviated notation for both linear and nonlinear terms in (16), and we consider the parameters Dp and ϕp on the trajectory at each step p. The time discretization for the first step of the proposed scheme (16) applied to system (19) can be written as

 
formula

The time discretization for the second step is as follows:

 
formula

The third and last step discretization is

 
formula

Following the expressions of (20), (21), and (22), respectively, for the three steps, the terms k2Δt, c2Δt, and δ can be considered to be dimensionless parameters. These parameters will be used to conduct a more general analysis of the characteristics of the proposed schemes and their comparison with those of the scheme studied by Cullen (2001). The three steps of the temporal discretization of the proposed schemes can be rewritten in matrix form using the vector variable Up = (Dp, ϕp)T at each step p. The vector value U* of the first step is substituted into the second step, and the result obtained for Un+1/2 is substituted into the third step to obtain a fourth-degree polynomial equation for the amplification factor Ak. In  appendix B we offer further explanations of the matrix form of the proposed schemes and the characteristics of the polynomial equation of the amplification factor. This polynomial has one root equal to zero. Therefore, it can be reduced to a three-degree polynomial, which has one real root corresponding to the computational mode. The two other roots are complex conjugate numbers and correspond to the inertio-gravity waves. The computational mode is largely damped, as shown in Fig. 10 for the two cases δ = 0.75 and δ = 1.25. Note that one of the advantages of the proposed schemes is that they have a damped computational mode, as opposed to the method studied by Cullen (2001), which has one meteorological mode that is undamped (equal to unity).

Fig. 10.

Module of the computational mode for the proposed schemes (θ = 0.7).

Fig. 10.

Module of the computational mode for the proposed schemes (θ = 0.7).

In Fig. 11, we plot the module of the errors |AkAan| of the amplification factor Ak in the k2Δtc2Δt plane for the cases δ = 0.5 and δ = 1. To compare between the accuracy of the proposed schemes and the method studied by Cullen (2001), the errors are integrated using

 
formula

where σ = k2Δt, = c2Δt, and is the area of the domain used in the integration.

Fig. 11.

Module of the errors of the amplification factors |AkAan| for the proposed schemes (θ = 0.7).

Fig. 11.

Module of the errors of the amplification factors |AkAan| for the proposed schemes (θ = 0.7).

Table 2 shows the values of the errors E corresponding to the amplification factors, which are obtained by using the method studied in Cullen (2001), and those obtained by using the proposed schemes with the decentering parameters θ = 0.7 and θ = 0.75 for five values of the dimensionless parameter δ. Following the results of this test, we conclude that the proposed schemes provide more accurate results compared to those obtained by using the method studied by Cullen (2001). To further verify the resolution quality, the efficiency, and the convergence of the proposed schemes, in the following section several tests are performed by using other solutions of (19), which are the combination of two different modes.

Table 2.

Errors of the amplification factors using the proposed schemes and the method studied by Cullen (2001).

Errors of the amplification factors using the proposed schemes and the method studied by Cullen (2001).
Errors of the amplification factors using the proposed schemes and the method studied by Cullen (2001).

2) Numerical test using two modes

In this section the numerical test is performed using an analytical solution that is the real part of the combination of two modes with different phases:

 
formula

where

 
formula

and ν = νR + νIi is a complex number with as its conjugate. This analytical solution will be used as the initial conditions to test the accuracy of the schemes. The numerical test is performed using the parameter c = 0.01 and the initial conditions with νR = 1, νI = −1, and k = 0.01. Figure 12 shows the evolution of the errors for the parameters D and ϕ until time T = 50 using the proposed schemes and the scheme studied by Cullen (2001) with a time step Δt = 0.03. Note that for this test both methods lead to accurate results and the proposed schemes are highly accurate compared to the method studied by Cullen (2001). In the following section we demonstrate that the high accuracy of the proposed schemes has a great advantage for using large time steps that meet the Lipschitz criterion in order to reduce the computational cost, which makes them competitive in terms of efficiency.

Fig. 12.

Evolution of the errors until time T = 50 using the proposed schemes and the method studied by Cullen (2001) for an initial condition that combined two modes (νR = 1, νI = −1) with Δt = 0.03. Shown are errors for parameters (a) D and (b) ϕ.

Fig. 12.

Evolution of the errors until time T = 50 using the proposed schemes and the method studied by Cullen (2001) for an initial condition that combined two modes (νR = 1, νI = −1) with Δt = 0.03. Shown are errors for parameters (a) D and (b) ϕ.

3) Efficiency of the proposed schemes

The most expensive parts in terms of computational cost are the dynamical part, which is related to the implicit resolution of the linear term, and the physical part, related to the explicit resolution of the nonlinear term. The computational cost of the explicit resolution of the linear term, which is the case in the second step in the proposed schemes, is negligible compared to the computational costs of the dynamical and physical parts.

For the proposed schemes we have one dynamical part and one physical part in the first step, one physical part in the second step, and one physical part and one dynamical part in the third step. The method studied by Cullen (2001) has one dynamical part and one physical part in each step. The efficiency of the proposed schemes and the efficiency of the method studied by Cullen (2001) will be compared by using the time steps satisfying the condition Δtm = 1.25Δtc, where Δtm and Δtc are the time steps used for the proposed schemes and the method studied by Cullen (2001), respectively. This condition is sufficient to give a reliable comparison of the efficiency, depending on the number of physical and dynamical parts in each method.

Equation (19) with the gravity speed c = 0.01 is solved, for two cases of analytic solutions of the form (24), using the proposed schemes and the method studied by Cullen (2001). We consider k = 0.01, and in the first case we use the initial conditions with νR = 1 and νI = 1 and in the second case we consider νR = 1 and νI = 4. In Fig. 13 we plot the evolution of the errors of the solutions (D, ϕ)T, which are obtained by using the proposed schemes with a time step Δt = 0.030 and the method studied by Cullen (2001) with a time step Δt = 0.024. The results of the proposed schemes remain more accurate in comparison to those obtained by using the method studied by Cullen (2001). Therefore we conclude that the proposed schemes perform well in terms of efficiency.

Fig. 13.

Evolution of the error of the solution (D, ϕ)T until time T = 50 using the proposed method with Δt = 0.030 and the method studied by Cullen (2001) with Δt = 0.024 for an initial condition that combined two modes, for cases (a) νR = 1 and νI = 1 and (b) νR = 1 and νI = 4.

Fig. 13.

Evolution of the error of the solution (D, ϕ)T until time T = 50 using the proposed method with Δt = 0.030 and the method studied by Cullen (2001) with Δt = 0.024 for an initial condition that combined two modes, for cases (a) νR = 1 and νI = 1 and (b) νR = 1 and νI = 4.

Our numerical tests show that the method used to compare the competitive efficiency by considering the time steps that satisfy the condition Δtm = 1.25Δtc is valid for the large interval of time steps that meet the Lipschitz criterion. In the following we give another numerical test to study the efficiency of the proposed method using the two-frequency system in (17). In this example we will consider the initial value of the trajectories x(0) ∈ [0, 1] and the analytical solution ψ(x, t) = (x + ρxt)e(iΩ+)t, where ρ = 0.8, ω = 1.0, Ω = 1.8, and the trajectories are given by x(t) = x(0)/(1 + ρt).

We use an approach similar to that in Smolarkiewicz and Pudykiewicz (1992) to satisfy the Lipschitz criterion. The necessary condition the time step Δt should satisfy to avoid the intersection of the computed trajectories is

 
formula

which leads to the restrictive condition for the time step Δt < 0.625.

For each value of Δt, we study the error defined by

 
formula

where , ψn is the numerical solution, and ψexact is the exact solution.

Figure 14 (right) shows the error as a function of the variable time step Δt for the proposed method. The efficiency of the proposed method and the efficiency of the method studied by Cullen (2001) are compared by using the time steps that satisfy the condition Δtm = 1.25Δtc. The comparison is performed using the large interval of the time steps that meet the Lipschitz criterion. Figure 14 (left) shows the error for the method studied by Cullen (2001) using a time step Δtc and the error for the proposed method using a time step Δtm = 1.25Δtc. The results of the proposed method remain accurate, which confirms that the proposed method performs well in terms of efficiency.

Fig. 14.

(a) The error function for the method studied by Cullen (2001) using the time step Δt and the error function for the proposed method using the time step 1.25Δt. (b) The error function for the proposed method using the time step Δt.

Fig. 14.

(a) The error function for the method studied by Cullen (2001) using the time step Δt and the error function for the proposed method using the time step 1.25Δt. (b) The error function for the proposed method using the time step Δt.

4) Convergence of the proposed schemes

In this section we will consider the variable U = (D, ϕ)T in our analysis. To study the convergence of the proposed schemes, we use the following error for each global time T:

 
formula

where Un and U(tn) are the numerical solution and the analytical solution, respectively, at time tn = nΔt. The time step is Δt = T/N, and N is the total number of steps necessary to reach the global time T. In the analysis, we seek for a small Δt an error model of the form E(T) = CΔtr, where E(T) is calculated using (28). The least squares method is used to find the values of the parameter C and the order of convergence r. The tests of convergence of the schemes are performed using the global time T = 50 for two cases. In the first case we use c = 0.1, and the initial conditions are defined by setting νR = 1, νI = 4, and k = 0.1. For this case we obtain C = 2.44 × 10−11 and r = 2.03 for the proposed schemes, and C = 2.46 × 10−11 and r = 2.02 for the method studied by Cullen (2001). The proposed schemes are better in terms of convergence compared to the method studied by Cullen (2001). In some cases we observe the same order of convergence, such as the second case in which we use c = 0.01, k = 0.01, νR = 1, and νI = 1. For this case we obtain C = 8.02 × 10−6 and r = 2 for the proposed schemes and C = 8.04 × 10−6 and r = 2 for the method studied by Cullen (2001). For this case, in which we have the same order of convergence, the coefficient C for the proposed schemes is less than the coefficient obtained for the method studied by Cullen (2001). This is justified by the fact that the proposed schemes are more accurate, as reported in the previous sections.

6. Conclusions

This paper presents theoretical and numerical analyses of the properties of more complex methods. We try to avoid the problems of instability associated with the treatment of the nonlinear part of the forcing term. A class of semi-Lagrangian semi-implicit schemes is proposed that uses a modified TR-BDF2 method based on the trapezoidal rule and the second-order backward differentiation formula. For the nonlinear term we used two stages as the predictor and corrector in the TR method and one stage for the BDF2 method. A family of second-order approximations derived by Gospodinov et al. (2001) is used in the first stage as the predictor of the TR method. For the linear term we used the implicit trapezoidal method in the first step, the explicit trapezoidal method in the second step, and the implicit BDF2 method in the third step. The use of the corrector stage improves the stability of the proposed schemes, and a large stable imaginary part is obtained. Following the numerical tests, the second-order approximation using α = ¼, which corresponds to the approximation of SETTLS, is the appropriate choice that guarantees the large domain of absolute stability and the optimal intersection of the region of absolute stability with the imaginary axis. This intersection is improved by using a decentering in the second step to obtain the schemes that perform well for purely oscillatory solutions. The numerical analysis demonstrates that the proposed class of semi-Lagrangian semi-implicit schemes performs well in terms of stability, accuracy, convergence, and efficiency. The potential practical application of the proposed family of schemes to a weather prediction model or any other atmospheric model is not analyzed and could be the subject of other forthcoming studies.

Acknowledgments

The authors are grateful to the anonymous reviewers for their valuable comments and suggestions, which helped to improve the quality of the paper significantly. The research was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) and Environment Canada.

Appendix A

TR-BDF2 Method

The composite second-order trapezoidal rule–backward differentiation method can be used to numerically solve the following ordinary differential equation:

 
formula

We use un to denote the computed approximation to the solution at time tn = nΔt. The original TR-BDF2 scheme is performed in two steps. In the first step, we apply the trapezoidal rule (TR) to advance our solution from tn to tn+γ,

 
formula

and then the second-order backward differentiation formula (BDF2) is used to advance the solution from tn+γ to tn+1:

 
formula

In our case we use γ = ½. Equations (A2) and (A3) can be solved by using implicit or explicit formulas for the source term.

APPENDIX B

Matrix Form of the Proposed Scheme and Trajectory Evaluation

a. Matrix form of the proposed scheme

The three steps given in section 5 by (20)(22) for temporal discretization of the proposed schemes applied to solve (19) can be rewritten in the following matrix form using the vector variable Up = (Dp, ϕp)T at each step p:

 
formula

The matrix j, j = 1, 2, …, 8 are given as

 
formula
 
formula
 
formula

The system (B1) can be rewritten in the following form:

 
formula

where and .

We obtain the equation of the amplification factor P(τ) = 0, where P is a fourth-degree polynomial given by

 
formula

where is the identity matrix 2 × 2. We have P(0) = det(−) = 0 since 3 is a singular matrix. If we exclude the zero root, we obtain a three-degree polynomial, which has one real root and the two other roots are complex conjugate numbers. The coefficients of matrices j can be expressed as a function of the parameters k2Δt, c2Δt, and δ = (c0/c)2, which are used in section 5 for the analysis.

b. Trajectory evaluation

We consider the following equation of the semi-Lagrangian trajectory:

 
formula

where x is the position vector of an air parcel and V is the corresponding three-dimensional velocity.

Following Hortal (2002), the approximation of the SETTLS method is obtained by using the Taylor expansion of the unknown position vector x around the departure point of the semi-Lagrangian trajectory:

 
formula

where AV indicates the approximation of the average value along the semi-Lagrangian trajectory between the arrival point A at time t + Δt and the corresponding departure point D at time t. Equation (B8) can be rewritten in the following form:

 
formula

Hortal (2002) proposed the following approximation of the average of acceleration along the semi-Lagrangian trajectory between the departure point D at time t and the arrival point A at time t + Δt,

 
formula

which leads to the following equation of the semi-Lagrangian trajectory:

 
formula

In the proposed method, for each arrival point xj at time t + Δt, the departure point is found by iteratively solving the previous equation.

Following Hortal (2002), SETTLS is based on the approximate value given by (B10) for the average of the acceleration . The proposed method involves the evaluation of terms at the middle point of the semi-Lagrangian trajectory. The middle point is located in the part of the SL trajectory where the average of acceleration is estimated. This value is used in the proposed method to obtain the position of the middle point, denoted by xtt/2, of the SL trajectory. We will consider the Taylor expansion, using Δt/2, of the unknown position vector x around the departure point of the SL trajectory:

 
formula

where the average of the acceleration is obtained using (B10). Therefore, the position of the middle point is obtained by using an explicit method, and the only iterative equation is the one used to obtain the departure point .

REFERENCES

REFERENCES
Bonaventura
,
L.
,
2000
:
A semi-implicit, semi-Lagrangian scheme using the height coordinate for a nonhydrostatic and fully elastic model of atmospheric flows
.
J. Comput. Phys.
,
158
,
186
213
, doi:.
Clancy
,
C.
, and
J. A.
Pudykiewicz
,
2013
:
A class of semi-implicit predictor-corrector schemes for the time integration of atmospheric models
.
J. Comput. Phys.
,
250
,
665
684
, doi:.
Cullen
,
M. J. P.
,
2001
:
Alternative implementations of the semi-Lagrangian semi-implicit schemes in the ECMWF model
.
Quart. J. Roy. Meteor. Soc.
,
127
,
2787
2802
, doi:.
Dharmaraja
,
S.
,
2007
:
An analysis of the TR-BDF2 integration scheme
. M.S. thesis, School of Engineering, Massachusetts Institute of Technology, 76 pp.
Durran
,
D. R.
,
2010
: Numerical Methods for Fluid Dynamics: With Applications to Geophysics. 2nd ed. Texts in Applied Mathematics, Vol. 32, Springer, 516 pp.
Durran
,
D. R.
, and
P. A.
Reinecke
,
2004
:
Instability in a class of explicit two-time-level semi-Lagrangian schemes
.
Quart. J. Roy. Meteor. Soc.
,
130
,
365
369
, doi:.
Gospodinov
,
I. G.
,
V. G.
Spiridonov
, and
J.-F.
Geleyn
,
2001
:
Second-order accuracy of two-time-level semi-Lagrangian schemes
.
Quart. J. Roy. Meteor. Soc.
,
127
,
1017
1033
, doi:.
Hortal
,
M.
,
2002
:
The development and testing of a new two-time-level semi-Lagrangian scheme (SETTLS) in the ECMWF forecast model
.
Quart. J. Roy. Meteor. Soc.
,
128
,
1671
1687
, doi:.
McDonald
,
A.
, and
J. R.
Bates
,
1987
:
Improving the estimate of the departure point in a two-time-level semi-Lagrangian and semi-implicit model
.
Mon. Wea. Rev.
,
115
,
737
739
, doi:.
McDonald
,
A.
, and
J.
Haugen
,
1992
:
A two-time-level, three-dimensional semi-Lagrangian, semi-implicit, limited-area gridpoint model of the primitive equations
.
Mon. Wea. Rev.
,
120
,
2603
2621
, doi:.
Ritchie
,
H.
,
1991
:
Application of the semi-Lagrangian method to a multilevel spectral primitive equations model
.
Quart. J. Roy. Meteor. Soc.
,
117
,
91
106
, doi:.
Ritchie
,
H.
, and
C.
Beaudoin
,
1994
:
Approximations and sensitivity experiments with a baroclinic semi-Lagrangian spectral model
.
Mon. Wea. Rev.
,
122
,
2391
2399
, doi:.
Ritchie
,
H.
,
C.
Temperton
,
A.
Simmons
,
M.
Hortal
,
T.
Davies
,
D.
Dent
, and
M.
Hamrud
,
1995
:
Implementation of the semi-Lagrangian method in a high-resolution version of the ECMWF forecast model
.
Mon. Wea. Rev.
,
123
,
489
514
, doi:.
Robert
,
A.
,
1981
:
A stable numerical integration scheme for the primitive meteorological equations
.
Atmos.–Ocean
,
19
,
35
46
, doi:.
Smolarkiewicz
,
P.
, and
J.
Pudykiewicz
,
1992
:
A class of semi-Lagrangian approximations for fluids
.
J. Atmos. Sci.
,
49
,
2082
2096
, doi:.
Staniforth
,
A.
, and
J.
Côté
,
1991
:
Semi-Lagrangian integration schemes for atmospheric models—A review
.
Mon. Wea. Rev.
,
119
,
2206
2223
, doi:.
Tanguay
,
M.
,
A.
Robert
, and
R.
Laprise
,
1990
:
A semi-implicit semi-Lagrangian fully compressible regional forecast model
.
Mon. Wea. Rev.
,
118
,
1970
1980
, doi:.
Temperton
,
C.
, and
A.
Staniforth
,
1987
:
An efficient two-time-level semi-Lagrangian semi-implicit integration scheme
.
Quart. J. Roy. Meteor. Soc.
,
113
,
1025
1039
, doi:.
Temperton
,
C.
,
M.
Hortal
, and
A.
Simmons
,
2001
:
A two-time-level semi-Lagrangian global spectral model
.
Quart. J. Roy. Meteor. Soc.
,
127
,
111
126
, doi:.
White
,
J. B.
, III
, and
J. J.
Dongarra
,
2011
:
High-performance high-resolution semi-Lagrangian tracer transport on a sphere
.
J. Comput. Phys.
,
230
,
6778
6799
, doi:.