## 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:

We assume that the two terms depend on the field variable *ψ*(*x*, *t*), *x*, and time. At each time *t*_{n} = *n*Δ*t*, where Δ*t* is the time step, the values of the parameters *ψ* and the velocity **V** are known for all times *t*_{p} = *p*Δ*t* with *p* = 0, 1, 2, …, *n*. The proposed method requires for each arrival point *x*_{j} the evaluation of the position at time *t*_{n} of the air parcel arriving at a point (*x*_{j}, *t*_{n+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 *x*_{D} = *x*(*t*) and the arrival point *x*_{A} = *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 *t*_{n} that arrives at (*x*_{j}, *t*_{n+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 (*x*_{j}, *t*_{n+1}). Following theorem 1 in Gospodinov et al. (2001), the following discrete approximations of (1),

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:

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 *N*^{n+1/4}. This value is obtained by using a trapezoidal treatment of *N*^{n} and the *α* approximation of the nonlinear term denoted by at the midpoint of the SL trajectory originating at and arriving at (*x*_{j}, *t*_{n+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 *N*^{n} 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 *N*^{n+1} on the right-hand side of the third equation in (3), where *N*^{n+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

where, for example, *N*(*x*, *t*_{n}):= *N*[*ψ*(*x*, *t*_{n}), *x*, *t*_{n}]

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

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

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

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* + *iω*Δ*t*, one obtains after simplification the following equation for the amplification factor *A*_{k}:

where , , and *s* = *x*_{A} − *x*_{D}.

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 *θ*.

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 *N*^{n+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.

In the case of the two-step method in which we consider the value of *N*^{n+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:

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*) = *e*^{ikx} and *σ* = *iω* in (7), which leads to an analytical amplification factor of *A*_{an} = *e*^{i(ωΔt−ks)}. Therefore, for an accurate scheme, the module of the numerical amplification factor *A*_{k} 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* = *iω*Δ*t* in (8). We use the *ω*Δ*t* − *ks* 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 |*A*_{k} − *A*_{an}| 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.

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:

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 *L*^{n+1} at time *t*_{n+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* + *iω*Δ*t*, one obtains the following equation for the numerical amplification factor *A*_{k}:

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

which has the module

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 *A*_{k}, 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

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 |*A*_{k}| 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 |*A*_{an}| = *e*^{λΔt} are represented in the same figure by the vertical lines.

Note that in the case in which decentering is considered, the value *μ* = 0.5 leads to an amplification factor *A*_{k} 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* = *iω*Δ*t* and under the assumption |*z*| ≪ 1 by

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.

## 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:

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:

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 *A*_{k}:

where

and

The complex numbers *z*_{1} and *z*_{2} are related to the fast and slow waves, respectively, and are given by *z*_{1} = *i*ΩΔ*t* and *z*_{2} = *iω*Δ*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 *A*_{k} 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 |*A*_{k} − *A*_{an}| 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.

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:

where, as previously defined, *z*_{1} = *i*ΩΔ*t* and *z*_{2} = *iω*Δ*t*.

In Fig. 9a we plot the module of the errors |*A*_{k} − *A*_{an}| 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.

### 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):

where *ϕ* is the geopotential, *D* is the divergence term, and *c* is the gravity speed. We consider *c*_{0} as the reference value of *c*, which is used to split the source term *c*^{2}*D* into two parts. The first part, , based on the reference value *c*_{0}, 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 *δ* = (*c*_{0}/*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 *D*^{p} 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

The time discretization for the second step is as follows:

The third and last step discretization is

Following the expressions of (20), (21), and (22), respectively, for the three steps, the terms *k*^{2}Δ*t*, *c*^{2}Δ*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 **U**^{p} = (*D*^{p}, *ϕ*^{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 **U**^{n+1/2} is substituted into the third step to obtain a fourth-degree polynomial equation for the amplification factor *A*_{k}. 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).

In Fig. 11, we plot the module of the errors |*A*_{k} − *A*_{an}| of the amplification factor *A*_{k} in the *k*^{2}Δ*t* − *c*^{2}Δ*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

where *σ* = *k*^{2}Δ*t*, = *c*^{2}Δ*t*, and is the area of the domain used in the integration.

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.

#### 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:

where

and *ν* = *ν*_{R} + *ν*_{I}*i* 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.

#### 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 Δ*t*_{m} = 1.25Δ*t*_{c}, where Δ*t*_{m} and Δ*t*_{c} 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.

Our numerical tests show that the method used to compare the competitive efficiency by considering the time steps that satisfy the condition Δ*t*_{m} = 1.25Δ*t*_{c} 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Ω+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

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

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

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 Δ*t*_{m} = 1.25Δ*t*_{c}. 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 Δ*t*_{c} and the error for the proposed method using a time step Δ*t*_{m} = 1.25Δ*t*_{c}. The results of the proposed method remain accurate, which confirms that the proposed method performs well in terms of efficiency.

#### 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*:

where **U**_{n} and **U**(*t*_{n}) are the numerical solution and the analytical solution, respectively, at time *t*_{n} = *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*Δ*t*^{r}, 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:

We use **u**^{n} to denote the computed approximation to the solution at time *t*_{n} = *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 *t*_{n} to *t*_{n+γ},

and then the second-order backward differentiation formula (BDF2) is used to advance the solution from *t*_{n+γ} to *t*_{n+1}:

### 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 **U**^{p} = (*D*^{p}, *ϕ*^{p})^{T} at each step *p*:

The matrix _{j}, *j* = 1, 2, …, 8 are given as

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

where and .

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

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 *k*^{2}Δ*t*, *c*^{2}Δ*t*, and *δ* = (*c*_{0}/*c*)^{2}, which are used in section 5 for the analysis.

##### b. Trajectory evaluation

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

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:

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:

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*,

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

In the proposed method, for each arrival point *x*_{j} 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 *x*^{t+Δt/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:

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

*Numerical Methods for Fluid Dynamics: With Applications to Geophysics*. 2nd ed. Texts in Applied Mathematics, Vol. 32, Springer, 516 pp.