1. Introduction
The accuracy and efficiency of weather and climate models has been greatly enhanced by the introduction of better numerical algorithms for the solution of the equations of motion. Two of the most notable advances were the development of the semi-implicit (SI) scheme for treating the gravity-wave adjustment process and the semi-Lagrangian scheme for the advection processes. For a general review of semi-Lagrangian schemes, see [
1,
2]; for a review of recent and emerging time integration schemes, see [
3].
Many modern operational NWP models use a semi-implicit scheme for time integration, increasing efficiency by enabling the use of a larger time step. However, this comes at a price: while stabilisation is achieved by slowing down the high-frequency gravity waves, the meteorologically significant components of the flow are also distorted by the time averaging of the SI scheme. For models that use a semi-Lagrangian advection scheme, the problem is magnified, as further increases in the time step result in greater errors.
It was pointed out in Lynch and Clancy [
4] that the Laplace transform (LT) method with analytic inversion gives an exact treatment of the linear modes. This is due to the fact that the LT scheme does not involve time-averaging of the linear terms. Harney and Lynch [
5] described a Laplace transform integration scheme in an Eulerian baroclinic model and showed that it yields more accurate forecasts than SI.
In this paper, we describe the implementation and performance of an integration scheme based on semi-Lagrangian advection and Laplace transform adjustment. The new scheme is incorporated in a spectral baroclinic atmospheric model
Peak, described in Ehrendorfer’s book [
6]. It is validated by comparison with a semi-Lagrangian semi-implicit scheme. Both schemes have been comprehensively verified against the Eulerian semi-implicit scheme originally implemented in the model and fully documented in [
6].
Table 1 shows four integration schemes. EuSI is the original Eulerian semi-implicit scheme described in [
6]. EuLT is the Eulerian Laplace transform scheme described in Harney and Lynch, 2019 [
5]. The two Lagrangian schemes are denoted LaSI and LaLT.
Numerical tests were carried out to confirm the correctness of the model implementations. A comparison of their performance in simulating a growing baroclinically unstable disturbance is described in
Section 4.1. The four simulations are quite similar, confirming the integrity of the codes.
The goal of this study is to demonstrate that Laplace transform adjustment outperforms the semi-implicit scheme for time steps typically used for semi-Lagrangian advection, and provides a means of producing more accurate weather forecasts and climate simulations. In
Section 2, we review the semi-implicit semi-Lagrangian schemes and outline the Laplace transform adjustment scheme. We perform a simple analysis to compare the accuracy of the two adjustment schemes. In
Section 3, we apply the Laplace transform scheme to the
Peak model equations. Treatment of orography and of diffusion are given particular attention. In
Section 4, a series of tests comparing simulations using the LaLT scheme and the semi-implicit scheme are described. Conclusions are presented in
Section 5. Two appendices are included, one on the calculation of the commutator, and one relating the Laplace transform scheme to exponential integration schemes.
2. Outline of the Semi-Implicit and Laplace Transform Schemes
Before considering the details of the
Peak model, we describe the general approach to integration using the LaSI and LaLT schemes. The original EuSI scheme is comprehensively documented by Ehrendorfer [
6] and the EuLT scheme by Harney and Lynch [
5].
2.1. Outline of the LaSI Scheme
We consider a general equation written in Lagrangian form
where
are the linear terms and
are the nonlinear terms. Before we convert this grid-point equation to the spectral domain, we discretise it in time using a three time-level Lagrangian scheme. For a trajectory from departure point
D at time
through arrival point
A at
, the equations are
where nonlinear terms are evaluated at the midpoint
M and superscripts −, 0 and + denote values at times
,
and
. We can write the solution at the advanced time as
where the right hand side may be computed from known quantities. We now convert Equation (
2) to spectral form, multiplying by
and integrating over the sphere to get
where
is the vector of spherical harmonic coefficients of
, and
is similarly related to
. The right hand side may be computed from known quantities. We note that the linear operator
commutes with the spectral transform.
2.2. Outline of the LaLT Scheme
We consider again the general equation in Lagrangian form (
1). Before converting this grid-point equation to the spectral domain, we take the Laplace transform
along the trajectory starting from the departure point
D at time
:
where
s is the complex variable conjugate to
t and
. As for the LaSI scheme, we assume that the nonlinear term varies slowly, and approximate it by its value at the mid-point of the trajectory from
to
.
To get an equation for
we have to change the order of the Laplace transform and the linear operator
; in general, these do not commute. We define the commutator
so that the term
becomes
, and write the transformed equation
For computation of the commutator, see
Appendix A.
We now convert to spectral form, integrating over the sphere, to get
where
is the vector of spherical harmonic coefficients of
,
are the coefficients of
, and
and
are the coefficients of
and
. The solution may be written
where, after inverse transformation to the time domain, the right hand side may be computed from known quantities.
2.3. Accuracy Analysis
An accuracy analysis of the Laplace transform (LT) and semi-implicit (SI) schemes was presented in [
5]. The main conclusions are reviewed here. Considering a simple oscillation equation
that represents a component of the full system used in numerical weather prediction, the analysis showed that LT is more accurate than SI for both the linear and nonlinear terms. Assuming that the nonlinear term varies slowly, we take it to be a constant
. The exact solution of (
3) at time
is then
where
is the solution at time
.
The SI approximation to (
3) is
Solving for the new value
, we have
Comparing this with the exact solution (
4), we find that both the linear and nonlinear components of the solution are misrepresented. For the exact solution (
4),
is multiplied by
, while for the SI solution (
5) the multiplier is
Thus, although there is no amplification, the phase error increases with
. For the nonlinear term, the SI scheme has both modulus and phase errors; for details, see Section 2 in [
5]. The errors in the SI scheme are significant when the time step is large.
We now apply the Laplace Transform
to Equation (
3), taking the origin of time to be
. The transform of (
3) is
The solution for
follows immediately:
The inverse Laplace transform
with time set to
yields
equal to the exact solution (
4) at time
. Thus, to the extent that the nonlinear term can be regarded as constant, the LT scheme is free from error.
2.4. Filtering with the LT Scheme
The LT scheme filters high frequency components by using a modified inversion operator
: this can be done numerically by distorting the Bromwich contour for the inversion integral to a closed curve excluding poles associated with the high frequencies, as in [
7,
8,
9,
10,
11]. In the present study, as in those of Lynch and Clancy [
4] and Harney and Lynch [
5], we invert the transform analytically, explicitly eliminating components with frequency greater than a specified cut-off frequency
.
Filtering may be done with a sharp cut-off at
, or with a smooth function such as a Butterworth filter having frequency response function
Formally, we can define the modified inversion operator as the composition of the filter and the inverse Laplace transform:
. In the numerical tests reported in
Section 4 we set
and choose the cut-off period
h and
. These values were chosen on the basis of many tests which showed that they yielded the best results for the given model resolution. Assuming that
, the inverse Laplace transform of (
6) at time
gives
This agrees with the exact analytic result (
4) for both the linear and nonlinear terms.
Equation (
8) is formally identical to Equation (
4) in Cox and Matthews [
12], for which the local truncation error is
. Assuming that the nonlinear term varies slowly, we have
, and the local truncation error is third order in time. We note that (
8) is also equivalent to a component of the vector Equation (
2) in Clancy and Pudykiewicz [
13]. For the relationship between Laplace transform integration and exponential integrators, see
Appendix B.
3. Laplace Transforming the Peak Equations
We now describe the application of the Laplace transform scheme to the
Peak model equations. We begin with Equations (9.38)–(9.41) in Section 9.3 of Ehrendorfer [
6], but written in Lagrangian form:
The notation is generally conventional and equations similar to these have appeared frequently, going back to Hoskins and Simmons [
14]. The dependent variables are vorticity (
), divergence (
), temperature (
) and log surface pressure
), where
Pa is the reference pressure. All variables are in the spectral domain but, for compactness, we suppress the spectral indices so that, for example,
is written simply as
. Bold-face variables are vectors with values at all
model levels, which are equally spaced in
-coordinates. The explicit expressions for the matrices
and
are given in Ehrendorfer [
6] and vectors are defined as follows:
,
and
. Linear damping with coefficient
may be applied to all variables except the surface pressure. The surface orography is represented by
.
Diffusion is generally required to control spurious oscillations. Explicit diffusion schemes may become unstable for large time steps. To avoid this problem, we may use a scheme that is fully implicit. However, when combined with the Laplace transform scheme, this leads to additional complications. To circumvent these, we employ an operator splitting method. In the first stage, the inviscid Equations (
9)–(
12) with
, are advanced to time
. This is then followed by a second stage in which the damping terms are integrated analytically (see
Section 3.8).
3.1. Lagrangian Departure Points and Values
The Cartesian coordinates of a grid point at latitude
and longitude
are
We compute
for each grid-point and store them; they are independent of the model level. We use a method described by McGregor [
15] to determine the departure points
. The inverse transformation is
(We use the Fortran function atan2 here). The inverse transform must be done every time step, as the departure points change. They also depend on the model level.
Values of the variables at the departure points are obtained by interpolation in
space. A bicubic interpolation scheme is used. To facilitate interpolation, the grid point arrays are extended by two rows or columns around the boundaries. In the east–west direction, the extension is cyclic. At the poles, we repeat the two boundary rows, reversing their order and shifting the longitude by
to allow for crossing the poles [
15].
3.2. Orography
Preliminary tests with real data produced instabilities that were insensitive to damping and to spectral truncation. The problem first became manifest over the Andes and Himalayas. As is well known, orography can lead to problems with semi-Lagrangian integration schemes. Ritchie and Tanguay [
16] proposed a modification that alleviates this problem. An orographic term is subtracted from the pressure, yielding a much smoother field that is more accurately interpolated to the departure points.
The method of Ritchie and Tanguay is used in the IFS model and described in the documentation [
17]. The variable
in (
12) is separated into a constant part involving the orography and a variable part independent of orography:
where
is, by definition, the value for an isothermal hydrostatic atmosphere. Thus,
. We can then write
The term
involving the gradient of orography is computed in an Eulerian manner, and the continuity Equation (
12) becomes
The advected variable is much smoother than the original variable since the underlying orography has been extracted.
3.3. Inviscid Stage
We integrate the equations using a split scheme, where the diffusion is omitted in the first stage and integrated analytically in the second. The nonlinear terms are evaluated at the central time and averaged between the departure and arrival points:
This averaging is similar to the approach adopted in the ECMWF IFS model [
17].
We now take the Laplace transform of (
10)–(
12), after applying (
15) to change the pressure variable. Since the Laplace operator
does not commute with the spacial Laplacian
, we introduce the commutator
The factor
is included here to ensure that
is independent of
s (see
Appendix A). The transformed equations may be written
Using the equations for
and
, these quantities are eliminated from the divergence equation to obtain an equation for a single variable,
:
We transfer all terms involving
to the left:
3.4. Transforming to Vertical Eigenmodes
We now define the vertical structure matrix
Since
is symmetric, its eigen-structure can be written
We now transform to vertical eigen-modes by multiplying (
21) by
:
We note that, for a specific spectral component with total wavenumber
ℓ, the Laplacian has a simple form
. Then, defining
, the matrix on the left hand side of (
22) is
Since this is a diagonal matrix, the equation falls into
separate scalar equations, one for each vertical mode. We group the right hand terms of (
22) according to powers of
s:
where the vectors
,
and
are
Multiplying (
22) by the inverse of the diagonal matrix
, the equation for the
k-th component is
We apply the operator
to (
24), noting that the vertical transform and Laplace transform commute. The terms can be inverted using standard results from Laplace transform theory [
18]. The value at time
is denoted by a + superscript:
The filter response function
was defined in Equation (
7) above.
3.5. Inverse Vertical Transformation
Let us define four diagonal matrices
(
will be needed below). Then (
25) can be written
We can now calculate the divergence at the advanced time,
For compactness, we define the propagation matrices:
(
will be used below). Then we can write the solution as
The -matrices can be pre-computed and stored, since they do not depend on the model variables.
3.6. Temperature and Pressure
We return to Equations (
19) and (
20):
Noting that
and
, dividing these equations by
s and applying the operator
at time
, we have
Both (
28) and (
29) require computation of
This term involves a convolution integral that may be approximated by the trapezoidal rule
where
is the average of old and new values of
. However, this method was found not to perform in a satisfactory manner. We therefore employed an alternative strategy: the term (
30) was computed by noting that
We divide (
24) by
s to give
We invert this using standard results for Laplace transforms to obtain
Then using the
-matrices, we can write
Noting (
31) and using the
-matrices, we can now write
Finally, using (
28) and (
29), the values of
and
are
3.7. Integrating the Vorticity
To complete the inviscid stage of the time step, we take the Laplace transform of the vorticity Equation (
9)
Dividing by
s and applying
at
yields
which is a standard centred Lagrangian step along the trajectory. As the poles are at
, the Laplace transform has no filtering effect here.
As is usual with the leapfrog model, a Robert-Asselin filter [
19] is applied to the prognostic variables to prevent separation of the solutions at odd and even time steps. The coefficient is fixed at
in all cases.
3.8. Diffusion Stage
The governing equation for a spectral component of any of the variables
,
or
may be written in the form
We split the right hand side into two parts and integrate them separately. We first integrate the inviscid equation
over a time interval
with initial condition
and denote the result as
. We then integrate the equation
analytically with the initial condition
to get
which is the required solution at time
. This second stage is applied to the divergence, vorticity and temperature; the surface pressure is not damped in this way.
The parameter
depends only upon the horizontal scale. The diffusion is assumed to be of the form
In the spectral domain, the damping coefficient becomes
We note that the spectral equations are unchanged in form by the addition of the hyper-diffusion term (); only the value of the coefficient is changed. The -term more strongly damps the smaller scales. Having applied diffusion, we have all the model variables at the advanced time, and a new time step can be taken.
4. Numerical Evaluation of the Integration Schemes
In this section, we describe a series of tests comparing simulations using the Laplace transform scheme (LaLT) and the semi-implicit scheme (LaSI). The LaSI and LaLT models were run with 20, 40 and 60 min time steps. For LaLT, a cut-off value h was set for all time steps. In most cases, the reference forecast was an integration of the LaSI model with a time step min.
Eulerian models are subject to a Courant-Friedrichs-Lewy stability condition. For an advection speed m/s, the non-dimensional stability ratio is unity for a time step s or 25 min. In fact, both the Eulerian models, EuSI and EuLT, were found to be unstable for a time step of 24 min. The Lagrangian models are not subject to this limitation.
The horizontal resolution of the model was at triangular truncation . The colocation grid corresponding to this has grid points, with a grid interval of approximately 150 km. In all cases, there were 20 vertical levels, uniformly spaced in -coordinates.
The default setting of the diffusion coefficient was m s : the damping of a component of total wavenumber ℓ is s . The default value implies an e-folding time of 2.2 h for the shortest waves represented at truncation T85. Sixth order diffusion ( m s ) was also applied for runs with a 60 min time step. We note that LaLT consistently required less explicit diffusion than LaSI.
4.1. Initial Validation Tests
Numerous tests were carried out to confirm the correct operation of the model codes. For short time steps, the Eulerian and Lagrangian advection produced similar results, as did the semi-implicit and Laplace Transform adjustment. As an example of the performance of the four models, simulations of a growing baroclinically unstable disturbance [
20] are shown in
Figure 1. The initial field is a zonally symmetric flow with a small perturbation on the Greenwich meridian. The four models were integrated for 12 days, each with a time step
min. The damping coefficient was
and all other parameter settings were equal for the four runs. The four simulations are quite similar, although we notice a slight damping for the Lagrangian runs, associated with the interpolations involved in the treatment of advection.
4.2. Kelvin Waves
Kelvin waves are eastward propagating waves that play an important role in atmospheric dynamics. Clancy and Lynch [
7] showed that the LT scheme had a significantly smaller phase error than the semi-implicit scheme for the integration of these waves. Exact Kelvin Wave initial conditions can be generated using the method of Kasahara [
21]. In this study, we use a simple analytical approximation described in [
22]. We examine the solutions for zonal wave numbers 1 and 4. The wave amplitude is 100 m in both cases. The theoretical period for the Kelvin wave with zonal wavenumber
is about 32 h and for
is about 8.3 h ([
21], Figure 9).
Figure 2 shows that the root mean square errors for wavenumber 1 for LaLT (blue lines) are significantly smaller than for LaSI (red lines). Forecasts were run with time steps of 20, 40 and 60 min (solid, dashed and dotted lines, respectively).
For the LaSI forecast of wave number 4 (
Figure 3) with a 20 min step, the propagation of the wave lags by a full wavelength by the end of the integration, reducing the rms error. The errors for the LaSI scheme with the larger time steps oscillate wildly as the wave moves into and out of phase with the reference solution. The errors for LaLT are much smaller and behave in a more realistic and steady manner.
4.3. The Five-Day Wave
The Five-day Wave RO(1,2) is the gravest symmetric rotational Hough mode of zonal wavenumber 1. It is closely related to the initial state chosen by Lewis Fry Richardson for their preliminary shallow-water forecast experiment ([
23], Section 4.1). The initial conditions for a three-dimensional Five-day Wave were implemented in
Peak. The initial pressure amplitude was 10 hPa, with mean pressure 1000 hPa. As no zonal mean flow was included, the wave has a period close to 5 days.
Figure 4 shows the root mean square errors in surface pressure for the LaSI scheme (red) and the LaLT scheme (blue), with time steps of 20 min (solid lines), 40 min (dashed lines) and 60 min (dotted lines). The reference is an SI forecast with time step of 10 min. The error level for LaLT is significantly less than that of the LaSI scheme, especially for the longer time step. The scores for vorticity at 250 hPa (
Figure 5) confirm the superior performance of LaLT.
4.4. Rossby-Haurwitz Wave
Rossby-Haurwitz (RH) waves are exact solutions of the nonlinear barotropic vorticity equation. While they are not eigenfunctions of the shallow water equations, they have frequently been used as test cases. Following Phillips [
24], the RH(4,5) wave was chosen as Test Case 6 by Williamson et al. [
25]. This test case has been extended to three dimensions; the initial vorticity field is as in the barotropic case, the divergence is zero and a vertical temperature profile and surface pressure field are defined; for details, see [
26].
The LaSI and LaLT schemes, with time steps of 20, 40 and 60 min, were compared to a reference run of LaSI using a time step of 10 s. No diffusion was used for the reference or LaLT runs, but the LaSI runs were unstable. This was overcome by applying horizontal diffusion with a coefficient
m
s
.
Figure 6 shows the root mean square errors in surface pressure for the LaSI scheme (red) and the LaLT scheme (blue). The error level for LaLT is substantially less than for the LaSI scheme.
Figure 7 shows scores for vorticity near the tropopause (250 hPa). The forecasts for LaSI and LALT have comparable errors, but there is a slight advantage for the Laplace transform scheme.
4.5. Flow over a Mountain
Test Case 5 of Williamson et al. [
25] treats a zonal flow over an isolated mountain. The mountain is centred at
with maximum height 2000 m. No analytic solution is known so, as usual, we take the LaSI run with
min as a reference.
Figure 8 shows the root mean square errors in surface pressure for the LaSI scheme (red) and LaLT scheme (blue), with time steps of 20 min (solid), 40 min (dashed) and 60 min (dotted lines). The error levels for the two schemes are very close in value.
4.6. A Baroclinically Unstable Wave
Polvani et al. [
20] devised a test case for baroclinic instability. The initial conditions consist of a non-divergent zonal flow with constant surface pressure. A small perturbation is added to the temperature to trigger the development of baroclinic instability. The test case of Polvani et al. was used by Ehrendorfer [
6] to validate the
Peak model. With a fixed value for the diffusion coefficient, the initial conditions are ‘numerically convergent’ as shown in [
20] using two different numerical models. A test case quite similar to that of Polvani et al. was constructed by Jablonowski and Williamson [
27].
We use the test case of Polvani et al. [
20] to show that the LaLT scheme can accurately simulate baroclinic development. Using this case, Harney and Lynch [
5] showed that the EuLT scheme can accurately simulate baroclinic development. In
Figure 1, above, we showed the 12-day forecasts for all four schemes all with a small time step
min. There was no substantive difference between the four schemes. They are also indistinguishable from the results plotted in ([
20],
Figure 4). Thus, all four schemes are capable of forecasting baroclinic development with high precision.
Quantitative scores confirm that the differences in performance are small:
Figure 9 shows the root mean square errors in surface pressure (hPa) for forecasts with the LaSI scheme (red) and LaLT scheme (blue), with time steps 20, 40 and 60 min. As usual, the reference forecast is LaSI with a 10 min time step. It is clear that the time truncation error grows with forecast range. It is also clear that the error is greater for larger time steps. The important point is that the errors for the two integration schemes are very similar in magnitude.
4.7. Real Data Test
The simple wave tests described above indicate superior accuracy for the LT scheme compared to the semi-implicit scheme. The ultimate conclusion on superiority of the scheme must involve comprehensive comparisons for a large range of meteorological conditions. As a first step, a single test using real atmospheric data is described here.
Data was retrieved from the European Centre for Medium-Range Weather Forecasts MARS archive. The date chosen was 00 UTC on 15 October 2017, the day before a major storm,
Ophelia, reached Ireland. This data comprised temperature, divergence and vorticity fields on 25 pressure levels, surface pressure and the relevant orography field. These fields where interpolated onto the 20 sigma levels and reduced to the spectral resolution T85 used for the
Peak forecasts. The process of interpolation introduced noise, which was removed by initialization, as described by Harney and Lynch [
5]. Using initialized data, six forecasts were performed using the LaSI and LaLT schemes with time steps of 10, 20 and 40 min.
Figure 10 shows the root mean square error for surface pressure. The reference is a forecast using LaSI with a time step of 5 min. The red curves are for LaSI and the blue ones for LaLT. For the 10 min step, the errors are comparable for the two models, although the error during the initial day is smaller for LaLT. For the larger time steps, the Laplace transform scheme is clearly superior to the semi-implicit scheme. Scores for mid-troposphere vorticity (
Figure 11) confirm the superior performance of LaLT.
5. Conclusions
An integration scheme using a Laplace transform for adjustment combined with a semi-Lagrangian advection scheme (LaLT) has been found to yield results at least as accurate as the popular semi-implicit semi-Lagrangian scheme. The numerical tests described in
Section 4 clearly indicate the superior accuracy for the LaLT scheme compared to the semi-implicit scheme LaSI. The single experiment with real data reinforces this advantage.
Ultimate conclusions on the superiority of the LaLT scheme require more comprehensive comparisons for a large range of meteorological conditions. The potential for operational implementation of LaLT would depend upon more exhaustive testing with higher spatial resolution and incorporating a full package of physical processes.
There are well-known advantages of using a two time level scheme for Lagrangian advection. There appear to be no difficulties in principle combining such a scheme with Laplace transform adjustment.
The efficient formulation of the LaLT scheme, with analytical inversion of the Laplace transform, is made possible through the use of a spectral model. An active debate on the future of spectral models has been ongoing for decades. The global simulation of an entire season of the Earth’s atmosphere, with a 1 km grid and upper boundary of 80 km [
28], suggests that spectral models will continue to be competitive in the future.
For practical reasons, the semi-Lagrangian method used in this study was applied only to the horizontal advection. However, there is no difficulty to include vertical advection in the scheme. The algorithmic complexity of the LaLT scheme is comparable to that of LaSI. Averaging over a range of integrations, the processing time for LaLT was just 6% more than for LaSI, so the Laplace transform scheme is computationally competitive.
In summary, the main result of this study is that the LaLT scheme is clearly superior in accuracy to the LaSI scheme for the large time steps typically used with semi-Lagrangian advection. The evidence presented gives a clear indication of the practical potential of the LaLT scheme for integrating weather and climate prediction models.