Alternatives to the Robert-Asselin filter

Patrick Marsaleix, Francis Auclair, Thomas Duhaut, Claude Estournel, Cyril

Nguyen, Caroline Ulses

► To cite this version:

Patrick Marsaleix, Francis Auclair, Thomas Duhaut, Claude Estournel, Cyril Nguyen, et al. Alternatives to the Robert-Asselin filter. Ocean Modelling, 2012, 41, pp.53-66. 10.1016/j.ocemod.2011.11.002. hal-02110112

HAL Id: hal-02110112

https://hal.science/hal-02110112

Submitted on 28 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Alternatives to the Robert-Asselin filter

Patrick Marsaleix<sup>a,*</sup>, Francis Auclair<sup>a</sup>, Thomas Duhaut<sup>a</sup>, Claude Estournel<sup>a</sup>, Cyril Nguyen<sup>b</sup>,
Caroline Ulses<sup>a</sup>

The Leap Frog time stepping scheme (hereafter LF) partly loses its conservation properties when a Rob-
ert-Asselin filter (hereafter RA) is used to damp the computational mode. The LF + RA scheme actually
leads to a well-known long term attenuation of the physical mode. Besides, the stability of the LF, e.g.
the maximum permitted time step, is lowered by the use of the RA. Several methods, derived from the
Laplacian approach of Marsaleix et al. (2008), are presented as an alternative to the RA. It appears that
the physical mode is eventually much less impacted by higher order time filters. However, in some cases,
the stability of the time stepping scheme becomes worse than that of the LF + RA. A five points scheme
finally appears to preserve both the amplitude of the physical mode and the stability of the time stepping
scheme. The analysis of these filters is based on a triple approach: the kinetic energy balance, the ampli-
fication factors of the oscillation equation, numerical experiments performed with a 3D circulation ocean
model.

**1. Introduction**

The leapfrog time stepping scheme (thereafter the LF scheme) is
a second order, three time-level, and time centred scheme that has
been used for years in numerous General Circulation Models
(GCMs) (Mesinger and Arakawa, 1976).

There are principally two criticisms that are made about the LF
scheme. First of all, this scheme is less accurate than several other
methods, for instance the third-order Adams–Bashforth scheme
(Durran, 1991) or the predictor-corrector scheme (Shchepetkin
and McWilliams, 2005). The second criticism, upon which the pres-
ent paper focuses, concerns the Robert-Asselin time filter (Robert,
1966; Asselin, 1972). The latter is generally associated to the LF
scheme in order to counter the possible growth of the numerical
mode permitted by the Leapfrog scheme. This numerical mode is
characterised by the time splitting of the model solution into two
independent physical trajectories (respectively corresponding to
odd and even time steps). The Robert-Asselin time filter has how-
ever several drawbacks. First of all, it has a perceptible diffusive ef-
fect on the physical mode. Second, the properties of conservations

are not always clearly ensured (Leclair and Madec, 2009). Last, it
noticeably reduces the stability of the LF scheme (the computa-
tional cost is thus possibly increased by the requirement of a smal-
ler time step).

Despite evidence that other time stepping schemes can be more
efficient, the LF scheme is often used in weather and ocean models,
like for instance the widespread Princeton Ocean Model (Blumberg
and Mellor, 1987). As pointed out by Williams (2009), several rea-
sons can be invoked: the LF scheme is easy to implement and has
low run-time storage requirements. It is also computationally
inexpensive, notably because it requires only one evaluation of
the right-hand side of a prognostic equation.

Besides, on one hand, the accuracy of a time stepping scheme
largely depends on the ratio of the time step to the time scales of
the simulated processes. On the other hand, the stability of time
stepping methods is mostly constrained by the fastest gravity
waves, leading to a time step much smaller than the time scale
of the other processes. Thus the accuracy of the Leapfrog scheme
remains reasonably good over a wide range of ocean situations.

We may also note that complex numerical systems like modern
ocean models are globally coherent, e.g. a numerical scheme can
not be written independently from the other schemes of the model.
Amongst a lot of other possible examples, the time stepping

<footnote><sup>a</sup>Laboratoire d'Aérologie, CNRS and Toulouse University, 14 Avenue Edouard Belin, 31400 Toulouse, France</footnote>

<footnote><sup>b</sup>Observatoire Midi-Pyrénées, 14 Avenue Edouard Belin, 31400 Toulouse, France</footnote>

<footnote>* Corresponding author. Tel.: +33 (0)5 61 33 27 77.
E-mail address: Patrick.marsaleix@aero.obs-mip.fr (P. Marsaleix).</footnote>
scheme has a major impact on the way the internal and external modes are coupled in the POM and ROMS models (Ezer et al., 2002). Obviously, replacing a LF scheme by another time stepping method is not a trivial development.

For these reasons, and even though we may expect a growing diversity of numerical methods in the coming years, the Leapfrog will probably remain in use in several models. For the time being, the LF scheme combined to a Robert-Asselin filter is notably implemented in widespread ocean models like POM (Blumberg and Mellor, 1987) and NEMO (Leclair and Madec, 2009, 2011).

Our study has similarities with the papers of Williams (2009, 2011). Indeed, the objective of the present paper is clearly to focus on what we consider to be the principal weak point of the LF scheme, namely the Robert-Asselin filter used to damp the numerical mode permitted by the LF scheme. We thus propose to develop an alternative scheme to the actual Robert-Asselin filter. This new filter should bring significant improvements on the three shortcomings previously mentioned. Practically:

(1) It should let the physical mode unchanged.

(2) The tracer conservation should be unambiguously ensured.

(3) The stability should be higher than that of the LF scheme when combined with a Robert-Asselin filter.

Another expected requirement is that this new filter should be easy to implement in actual LF based models, e.g. it should not imperil the global consistency of these models. In other word, its implementation should not concern any other scheme than the time stepping scheme itself.

The paper is organised in the following manner. A modified formulation of the Robert-Asselin filter, namely the Laplacian scheme used in Marsaleix et al. (2008) (hereafter M08), is presented in Section 2. The contribution of the Laplacian and Robert-Asselin filters to the discrete energy balance of the physical and numerical modes is presented. Other filters, based on higher order schemes, are also examined. Amongst them, the FD filter is a “five points” scheme emerging as a convincing low-diffusive alternative to the Robert-Asselin filter. The amplification factor for the oscillation equation is analysed in Section 3. The different filters are finally tested on the basis of numerical experiments presented in Section 4.

## 2. A preliminary, physically based, examination of different filters

Several time filters derived from the Laplacian filter used in M08 will be considered in the following. We thus start to recall the differences between the Laplacian and Robert-Asselin filters.

### 2.1. Laplacian filter

Let us start with.

$$ \frac{\partial F}{\partial t} = \frac{\partial \psi}{\partial t} + \text{RHS}, \qquad (2.1) $$

where F represents any of the variables found in ocean models (e.g. velocities or tracers), and RHS the usual forcing terms of the momentum or tracer equations (Coriolis, pressure gradient, advective and diffusive terms). The laplacian time filter is $\frac{\partial \psi}{\partial t}$ with:

$$ \psi = K \frac{\partial F}{\partial t}, \qquad (2.2) $$

where K is a coefficient of diffusion (the word “diffusion” is used here for convenience, although K has the dimensions of time, which is not the usual dimensions of a diffusion coefficient). As the filter is the time derivative of $\psi$, the conserving properties discussed in

M08 are a consequence of the Green-Ostrogradski theorem. Indeed, the time integration of Eq. (2.1) leads to:

$$ F^t - F^0 = \psi^t - \psi^0 + \int_0^t (\text{RHS}) dt', \qquad (2.3) $$

where superscripts t and 0 refer to the present time and the initial time respectively. As we can expect $|\psi^t - \psi^0|$ to be much smaller than $|F^t - F^0|$, notably if we integrate over a long time, it happens that $F^t - F^0 \approx \int_0^t (\text{RHS}) dt'$. The conservation of F is not perfect (since $\psi^t - \psi^0 \neq 0$) but satisfying since the variation of F over time is quasi balanced by the RHS term.

### 2.2. Robert-Asselin filter

In M08 it is showed that the Robert-Asselin filter leads to transform (2.1) and (2.2) into:

$$ \frac{\partial F^t}{\partial t} = \text{RHS}^t + \frac{\Delta t}{2} \sum_{q=1,N-1} \left(\frac{v}{2}\right)^q \frac{\partial^2 F^{t-q\Delta t}}{\partial t^2}, \qquad (2.4) $$

where N is the number of iterations since the beginning of the simulation, $\Delta t$ is the time step increment, v is the Robert-Asselin coefficient and

$$ \frac{\partial^2 F^{t-q\Delta t}}{\partial t^2} = \frac{1}{\Delta t^2} \left( F^{t-(q-1)\Delta t} - 2F^{t-q\Delta t} + F^{t-(q+1)\Delta t} \right). \qquad (2.5) $$

The Robert-Asselin filter is thus equivalent to the discrete sum of decentred Laplacian type filters starting from the beginning of the simulation. Because of this, the conservation properties of the Robert-Asselin filter are possibly questionable since a cumulative effect of all these terms can not be excluded on the time filter balance. Some answers are provided by the discrete energy balance detailed in the following section. However, and realising that the Robert-Asselin coefficient v is a positive number generally much smaller than one (0.1 is a common value), we note that the terms $(\frac{v}{2})^q \frac{\partial^2 F^{t-q\Delta t}}{\partial t^2}$ of the discrete sum at the right hand side of (2.4) tend to vanish as the integer number q increases. If the Robert-Asselin coefficient is small enough, the discrete sum can be reasonably approximated to its first term (e.g. $\frac{\Delta t}{2} \sum_{q=1,N-1} (\frac{v}{2})^q \frac{\partial^2 F^{t-q\Delta t}}{\partial t^2} \approx \frac{v\Delta t}{4} \frac{\partial^2 F^{t-q\Delta t}}{\partial t^2}$) and becomes consequently almost equivalent to the Laplacian filter, provided that:

$$ K = \frac{v\Delta t}{4}. \qquad (2.6) $$

In this particular case, the conservation properties of the Laplacian and Robert-Asselin filters are obviously very close. We also retain that (2.6) gives a first (but possibly rough estimate if v is not small enough) estimate of an equivalent diffusion coefficient of the Robert-Asselin coefficient.

### 2.3. Kinetic energy balance

#### 2.3.1. Continuous time

In order to obtain a kinetic energy balance, Eq. (2.1) (where now F represents u, the velocity) is multiplied by the velocity and then integrated over time. For the sake of clarity the RHS term is arbitrarily dropped. This leads to:

$$ \int_{0}^{t} u \frac{\partial u}{\partial t} dt' = \int_{0}^{t} u \frac{\partial \psi}{\partial t} dt'. \qquad (2.7) $$

Supposing a zero velocity at t = 0, the left hand side of (2.7) is also:

$$ \int_{0}^{t} u \frac{\partial u}{\partial t} dt' = \frac{u(t)^2}{2}, \qquad (2.8) $$
and the right hand side of (2.7) leads to:

$$ 
\begin{aligned}
\int_0^t u \frac{\partial \psi}{\partial t} dt' &= \int_0^t \left( \frac{\partial u \psi}{\partial t} - \psi \frac{\partial u}{\partial t} \right) dt' \\
&= \int_0^t \frac{\partial u \psi}{\partial t} - K \left( \frac{\partial u}{\partial t} \right)^2 dt' \\
&= u^t \psi^t - \int_0^t K \left( \frac{\partial u}{\partial t} \right)^2 dt'.
\end{aligned} 
\quad (2.9) $$

As far as the two last terms at the right hand side of (2.9) are concerned, we note (provided that $K$ is not negligible) that $u\psi$ should be much smaller than $\int_0^t K(\frac{u\psi}{\partial t})^2 dt'$ since the latter is the integral of a quantity that is always positive. As a consequence, the energy balance of the time filter may be approximated by

$$ \int_0^t u \frac{\partial \psi}{\partial t} dt' \approx - \int_0^t K \left( \frac{\partial u}{\partial t} \right)^2 dt'. \quad (2.10) $$

The latter is systematically negative so that it is a sink of kinetic energy. It acts as a dissipation term.

### 2.3.2. Discrete time

We now consider the discrete time. Time is given by $t = n\Delta t$ where $n$ is an integer representing the number of iterations and $\Delta t$ is the time step increment. The time stepping scheme is the LF scheme. As far as (2.8) is concerned (and assuming a zero initial current) we now have:

$$ \int_0^t u \frac{\partial u}{\partial t} dt' = \sum_{0:N} \left( u^n \frac{u^{n+1} - u^{n-1}}{2\Delta t} \right) \Delta t = \frac{u^N u^{N+1}}{2}. \quad (2.11) $$

We first note that the last term at the right hand side of (2.11) is not exactly the square of a velocity, but the product of two consecutive velocities. It is consequently possible for the discrete kinetic energy to be negative in some cases. This occurs when the direction of the current is changing (e.g. when $u^N$ and $u^{N+1}$ have a different sign). This is far from being the most probable situation, except in (quite unrealistic) cases where the model solution is dominated by numerical instabilities characterised by a current reversal at even and odd iterations.

As the present paper focusses on one of the main drawbacks of the LF scheme, that is the possible time splitting of even and odd time steps, it can be suitable to separate the current into the additional contributions of a mean current (corresponding to the average of the solution at even and odd time steps) and a $2\Delta t$ periodic spurious oscillation. According to (2.11), we shall retain that the discrete energy of the physical modes should be positive when the discrete kinetic energy associated to the spurious oscillation (e.g. the divergence of even and odd time steps) is negative.

The numerical expression corresponding to (2.9) is:

$$ 
\begin{aligned}
\int_0^t u \frac{\partial \psi}{\partial t} dt' &= \sum_{0:N} u^n \frac{\psi^{n+1/2} - \psi^{n-1/2}}{\Delta t} \Delta t \\
&= u^N \psi^{N+1/2} - \sum_{0:N-1} \psi^{n+1/2} (u^{n+1} - u^n).
\end{aligned} 
\quad (2.12) $$

The last term at the right hand side of (2.12) corresponds to the dissipation of kinetic energy by the time filter (note that we assumed $u^0 = 0$). Following M08, the Laplacian filter involved in the computation of the velocity at time step $n+1$ is based on a scheme that is centred on time step $n-1$. This means that instead of the straightforward discretisation that would intuitively lead to:

$$ \psi^{n+1/2} = K(u^{n+1} - u^n)/\Delta t, \quad (2.13) $$

We rather use:

$$ \psi^{n+1/2} = K(u^n - u^{n-1})/\Delta t. \quad (2.14) $$

We will see that the discrete energy balance is, amongst others, a mean to explain why the choice of (2.14) is better than (2.13). Actually, we can compare the numerical expressions for the time filter balance respectively obtained with (2.13) and (2.14). We consider (2.13) first. The last term at the right hand side of (2.12) leads to:

$$ -\sum_{0:N-1} \psi^{n+1/2} (u^{n+1} - u^n) = -\sum_{0:N-1} K(u^{n+1} - u^n)^2 / \Delta t. \quad (2.15) $$

At first sight, (2.15) is in good agreement with the dissipation term obtained in the continuous space, given by (2.10), since both are always negative. On the other hand, using (2.14) leads to:

$$ -\sum_{0:N-1} \psi^{n+1/2} (u^{n+1} - u^n) = -\sum_{0:N-1} K(u^n - u^{n-1}) (u^{n+1} - u^n) / \Delta t. \quad (2.16) $$

We note that the right hand side of (2.16) should be negative (as expected) most of the time (since at physical frequencies ($u^n - u^{n-1}$) and ($u^{n+1} - u^n$) should be very close), but not in the case of the aforementioned $2\Delta t$ periodic spurious oscillation (since the latter would lead ($u^n - u^{n-1}$) and ($u^{n+1} - u^n$) to have opposite signs). This behaviour is in fact consistent with the discrete energy balance of the LF scheme. We indeed showed that the discrete kinetic energy associated to $2\Delta t$ oscillations is negative. The fact that (2.16) becomes positive simply permits to counter their growth. At low frequencies, (2.16) recovers the usual conclusions of the continuous space, in other words the time filter behaves as a dissipation term. On the other hand, the unconditionally negative scheme (2.15) will enhance the growth of $2\Delta t$ oscillations. The decentred form (2.14) is consequently better than the centred form (2.13) since the latter is potentially unstable. We now generalise these conclusions to other possibilities of centring in the numerical form of $\psi^{n+1/2}$. Practically the latter is rewritten:

$$ \psi^{n+1/2} = K(u^{n-m} - u^{n-1-m})/\Delta t. \quad (2.17) $$

where $m$ is an integer. The numerical forms (2.13) and (2.14) correspond to $m = -1$ and $m = 0$, respectively. Let us consider the case of $m = +1$. Instead of (2.15), (2.16), the energy balance of the time filter now leads to:

$$ -\sum_{0:N-1} \psi^{n+1/2} (u^{n+1} - u^n) = -\sum_{0:N-1} K(u^{n-1} - u^{n-2}) (u^{n+1} - u^n) / \Delta t. \quad (2.18) $$

As for the previous schemes (2.18) should be negative at low frequencies and thus should have the expected behaviour of a dissipation term. But in the case of a $2\Delta t$ oscillation, ($u^{n-1} - u^{n-2}$) and ($u^{n+1} - u^n$) have the same sign, leading (2.18) to recover the drawback of (2.15), practically to be unable to damp the oscillation.

These different results may be summarised as follows:

- For odd values of $m$ in (2.17), the Laplacian filter behaves as a dissipation term at low (physical) frequencies but favours the growth of $2\Delta t$ oscillations.

- For even values of $m$ in (2.17), the Laplacian filter behaves as a dissipation term at all frequencies.

### 2.4. A first step toward a hybrid filter preserving the low frequencies

We recall that one of the objectives of this paper is to build a time filter that prevents the time splitting of odd and even time steps and, at the same time, does not dissipate the physical modes. We will use the ambivalence of the Laplacian filter towards $2\Delta t$ oscillations, depending of the parity of $m$ in (2.17), to build hybrid schemes combining two Laplacian filters, one using $m = 0$ and the other one using $m = +1$ or $m = -1$, namely:
$$ \psi^{n+1/2} = +K(u^n - u^{n-1})/\Delta t - K(u^{n-m} - u^{n-1-m})/\Delta t \quad m = \pm 1. \quad (2.19) $$

Because of the opposite sign of its two terms, (2.19) should lead to the following properties:

- As far as 2Δt oscillations are concerned, the energy balance of the time filter will be positive. A damping effect is expected.

- As far as low frequencies are concerned, the second term of (2.19) is expected to balance the first one. The physical modes should remain unchanged.

We note that the use of (2.19) may be interpreted, in the continuous time, as a third derivative (hereafter TD) filter added at the right hand side of the model equations.

## 2.5. Dissipation of physical and numerical modes

In the following subsections, we examine the energy dissipated (or created) by different filters in two cases: a physical mode and a spurious 2Δt oscillation.

### 2.5.1. The Laplacian filter in M08

We recall that the Laplacian filter in M08 is based on the decented scheme (2.14), that is (2.17) with $m = 0$. The energy balance of the filter is given by (2.16). We first consider a case of 2Δt oscillation:

$$ u^n = u_0(-1)^n, \quad (2.20) $$

where $n$ is the number of iterations. Using (2.20) in (2.16), we obtain the energy balance of the filter, namely:

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx +Nu_0^2 4K/\Delta t. \quad (2.21) $$

Note that we made the hypothesis that the amplitude $u_0$ is constant over the time $NΔt$ of integration, in other words that $NΔt$ is not too long compared to the time scale of the dissipation effect. The right hand side of (2.21) is positive, as expected from the previous discussions in Section 2.3.2.

We now consider a case of low frequency oscillation:

$$ u^n = u_0 \cos(\omega n \Delta t + \phi). \quad (2.22) $$

Using (2.22) in (2.16) leads to:

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -\sum_{0,N-1} 2u_0^2 K \sin^2(\theta) [\cos(2\theta) - \cos(2\omega n \Delta t + 2\phi)]/\Delta t, \quad (2.23) $$

with

$$ \theta = \frac{\omega \Delta t}{2}. \quad (2.24) $$

As for the previous case, we assumed that $NΔt$ is short enough so that $u_0$ can be considered constant. The time step increment can also be defined such that the time of integration is a multiple of $\pi/\omega$ and simplify (2.23) according to:

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -2Nu_0^2 K \sin^2(\theta) \cos(2\theta)/\Delta t, \quad (2.25) $$

and last we assume that $Δt$ is much smaller than the period of the oscillation, in other words that $\theta \ll 1$, in order to approximate (2.25) to:

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -2Nu_0^2 K \theta^2/\Delta t. \quad (2.26) $$

We see that the kinetic energy dissipated by the filter principally depends on $\theta^2$.

### 2.5.2. The Laplacian filter with $m = \pm 1$ in (2.17)

In the case of the two schemes based on $m = -1$ or $m = +1$ in (2.17), and considering the spurious 2Δt oscillation (2.20), the energy balance of the filter is:

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -Nu_0^2 4K/\Delta t. \quad (2.27) $$

Compared to the result obtained with the Laplacian filter used in M08 (e.g. (2.21)), this energy balance is negative.

As far as the physical oscillation (2.22) is considered, the energy balance is:

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -\sum_{0,N-1} 2u_0^2 K \sin^2(\theta) \times [1 - \cos(2\omega n \Delta t + 2\theta + 2\phi)]/\Delta t, \quad (2.28) $$

in the case of $m = -1$ and

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -\sum_{0,N-1} 2u_0^2 K \sin^2(\theta) \times [\cos(4\theta) - \cos(2\omega n \Delta t - 2\theta + 2\phi)]/\Delta t \quad (2.29) $$

in the case of $m = +1$.

In both cases (e.g. $m = \pm 1$) and as expected from the discussions in Section 2.3.2, the filter dissipates the physical mode, but amplifies the spurious oscillation since (2.27) is negative.

### 2.5.3. The Robert-Asselin filter

As shown in Section 2.2, the Robert-Asselin filter is numerically equivalent to (2.4) and thus can be seen as a sequence of decented Laplacian filters cumulated from the beginning of the simulation. As far as current values of the Robert-Asselin coefficient are concerned (a typical value is $v = 0.1$), the terms of the discrete sum $\sum_{q=1,N-1} (\frac{v}{2})^q \frac{\partial^2 F^{t-q\Delta t}}{\partial t^2}$ become negligible when $q$ increases and we have seen that retaining the first term (e.g. $q = 1$) leads to the Laplacian filter used in M08. It is now convenient to approximate the Robert-Asselin scheme to the two first terms of the discrete sum, namely:

$$ \frac{\partial F^t}{\partial t} = RHS^t + \frac{\Delta t}{2} \left[ \frac{v}{2} \frac{\partial^2 F^{t-\Delta t}}{\partial t^2} + \left(\frac{v}{2}\right)^2 \frac{\partial^2 F^{t-2\Delta t}}{\partial t^2} \right], \quad (2.30) $$

that can be expressed similarly as (2.1) (e.g. $\frac{\partial F}{\partial t} = \frac{\partial \psi}{\partial t} + RHS$) provided that $\psi$ is suitably formulated, namely:

$$ \psi^{n+1/2} = \frac{\Delta t}{2} \left( \frac{v}{2} \frac{u^n - u^{n-1}}{\Delta t} + \left( \frac{v}{2} \right)^2 \frac{u^{n-1} - u^{n-2}}{\Delta t} \right). \quad (2.31) $$

Using (2.6), we obtain:

$$ \psi^{n+1/2} = K(u^n - u^{n-1})/\Delta t + 2K^2(u^{n-1} - u^{n-2})/\Delta t^2. \quad (2.32) $$

The latter corresponds to (2.17) with $m = 0$, plus $\frac{2K}{\Delta t}$ times (2.17) with $m = -1$. The dissipation of the spurious 2Δt oscillation (2.20) is thus easily deduced from the previous sections, e.g. by using the energy balances (2.21) and (2.26) (the latter is multiplied by $\frac{2K}{\Delta t}$ and then added to the former). We obtain:

$$ -\sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx +Nu_0^2 4K/\Delta t \left(1 - \frac{2K}{\Delta t}\right). \quad (2.33) $$

Since $\frac{2K}{\Delta t} = \frac{v}{2}$ is smaller than one, we have $0 < 1 - \frac{2K}{\Delta t} < 1$. This leads to the following two remarks:

- First, the condition for the damping of the spurious oscillation is satisfied since (2.33) is positive.

- Second, (2.33) is lower than (2.21) and thus the damping of the 2Δt spurious oscillation is more efficient in the case of the M08 Laplacian filter, or, alternatively, an equivalent damping effect is
obtained with the Robert-Asselin filter but with a stronger coefficient $v$. Supposing that $K^{RA}$ and $K^{LP}$ are the values of the diffusion coefficient in (2.33) and in (2.21) respectively, this is obtained provided that $K^{RA}(1 - \frac{2K^{RA}}{\Delta t}) = K^{LP}$. Using (2.6) in order to define $v^{RA} = 4K^{RA}/\Delta t$ and $v^{LP} = 4K^{LP}/\Delta t$, this can be alternatively formulated as:

$$v^{RA} \left(1 - \frac{v^{RA}}{2}\right) = v^{LP}, \qquad (2.34)$$

or alternatively $v^{RA} = 1 - \sqrt{1 - 2v^{LP}} \approx v^{LP} + \frac{v^{LP2}}{2} + \frac{v^{LP3}}{2} + \dots$. If the diffusion coefficient used in the M08 Laplacian filter is such that $v^{LP} = 0.1$, (2.34) says that an equivalent damping effect of the spurious $2\Delta t$ spurious oscillations is obtained with the Robert-Asselin filter provided that $v^{RA} \approx 0.1056$. We note that (2.6) and (2.34) are anyway equivalent for very small $v$.

In a similar way, the dissipation of the physical oscillation (2.22) can be deduced from (2.23) and (2.29). Compared to (2.26), the rough estimate of the filter energy balance obtained in the case of the M08 Laplacian filter (e.g. $-2Nu_0^2K\theta^2/\Delta t + \dots$), the Robert-Asselin filter leads to a slightly larger value, i.e. $-2Nu_0^2K\theta^2(1 + \frac{2K}{\Delta t})/\Delta t + \dots$. The physical mode is consequently more dissipated in the case of the Robert-Asselin filter.

### 2.5.4. TD filters

We first consider the TD filter based on (2.19) with $m = 1$, that is:

$$\psi^{n+1/2} = +K(u^n - u^{n-1})/\Delta t - K(u^{n-1} - u^{n-2})/\Delta t. \qquad (2.35)$$

As in Section 2.5.3, and now noting that (2.35) is simply (2.17) with $m = 0$ minus (2.17) with $m = 1$, the energy balance of this filter is simply deduced from previous calculus. In the case of the spurious $2\Delta t$ oscillation (2.20), the energy balance is given by (2.21) minus (2.27), that is:

$$- \sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx +Nu_0^2 8K/\Delta t. \qquad (2.36)$$

This is twice the energy balance of the M08 Laplacian filter. We thus deduce that a dissipation effect equivalent to that of the M08 Laplacian filter is obtained with a twice smaller coefficient $K$. In the case of the physical mode (2.22), the filter energy balance is deduced from the difference of (2.23) and (2.29), that is:

$$- \sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -2Nu_0^2 K \sin^2(\theta) (\cos(2\theta) - \cos(4\theta))/\Delta t. \qquad (2.37)$$

As in the previous sections, we suppose that $\theta$ (given by (2.24)) is much smaller than one and we thus approximate (2.37) to $-2Nu_0^2K(\theta^2 + \dots)(6\theta^2 + \dots)/\Delta t$ and obtain:

$$- \sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -12Nu_0^2 K \theta^4 / \Delta t. \qquad (2.38)$$

The latter depends on $\theta^4$ and is thus much smaller than the energy dissipated by the M08 Laplacian filter (e.g. $-2Nu_0^2K\theta^2/\Delta t$).

We now consider the TD filter based on (2.19) with $m = -1$, that is:

$$\psi^{n+1/2} = +K(u^n - u^{n-1})/\Delta t - K(u^{n+1} - u^n)/\Delta t. \qquad (2.39)$$

As far as the spurious $2\Delta t$ oscillation (2.20) is concerned, the filter energy balance is the same as for the $m = +1$ case, that is (2.36). In the case of the physical mode (2.22), the filter energy balance is deduced from the difference of (2.23) and (2.28), that is:

$$- \sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -2Nu_0^2 K \sin^2(\theta) (\cos(2\theta) - 1)/\Delta t \qquad (2.40)$$

when $\theta \ll 1$ (2.40) can be approximated to:

$$- \sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx +4Nu_0^2 K \theta^4 / \Delta t. \qquad (2.41)$$

As previously, (2.41) is several orders smaller than the energy balance of the M08 Laplacian filter. On the other hand, this filter no longer behaves like a dissipation term because of the positive sign in the right hand side of (2.41).

### 2.6. FD filter

At this stage, a combination of different TD schemes can be considered in order to improve the selective properties of the resulting filter. Practically, a linear combination of (2.35) and (2.39) leads to:

$$\psi^{n+1/2} = +K(u^n - u^{n-1})/\Delta t - (1 - \alpha)K(u^{n-1} - u^{n-2})/\Delta t + \alpha K(u^{n+1} - u^n)/\Delta t, \qquad (2.42)$$

with $0 \le \alpha \le 1$. We will see in a following section that this new filter can be seen as a fourth derivative (hereafter FD) filter in the continuous time. The special case $\alpha = 0$ corresponds to the TD filter with $m = 1$ (that is (2.35)) and the case $\alpha = 1$ to the TD filter with $m = -1$ (e.g. (2.39)). The corresponding energy balance is simply deduced from an analogous linear combination of the results given in the previous section. As far as the spurious $2\Delta t$ oscillation (2.20) is concerned, the energy balance remains unchanged (that is 2.36). In the case of the physical mode (2.22), the filter energy balance is $(1 - \alpha)$ times (2.37) plus $\alpha$ times (2.40), that is:

$$- \sum_{0,N-1} \psi^{n+1/2} (u^{n+1} - u^n) \approx -2Nu_0^2 K \sin^2(\theta) (\cos(2\theta) - \alpha) - (1 - \alpha) \cos(4\theta)/\Delta t. \qquad (2.43)$$

For $\theta$ small enough, a Taylor development leads to:

$$\cos(2\theta) - \alpha - (1 - \alpha) \cos(4\theta) = -\frac{4\theta^2}{2!}(4\alpha - 3) + \frac{\theta^4}{4!}(2^4 + 4^4(\alpha - 1)). \qquad (2.44)$$

The leading order term cancels for $\alpha = \frac{3}{4}$ for which the filter effect on the physical mode should be much smaller than that of the TD filters. On the other hand, using $\alpha = \frac{3}{4}$ and (2.44) in (2.43) leads to $-2Nu_0^2K\sin^2(\theta)(\frac{\theta^4}{4!}(2^4 - 4^3)/\Delta t)$. As the latter expression is positive, we can expect the time stepping scheme to be slightly unstable when $\alpha = \frac{3}{4}$. We will see in Section 3.6 that a smaller value of the tuning coefficient $\alpha$ is actually preferable, although $\alpha \ne \frac{3}{4}$ leads the FD accuracy to recover the $\theta^4$ dependency characterising the TD filter.

### 2.7. Partial conclusions

(1) About the Laplacian and Robert-Asselin filters: the diffusive impact on the physical mode is lower, and the filtering of the $2\Delta t$ oscillation is more efficient, in the case of the Laplacian filter. However, the two methods tend to become equivalent for small values of $\alpha$ and $K$.

(2) The generalised form of the Laplacian filter (based on (2.17)) dissipates the physical mode whatever the value of $m$ in (2.17). The $2\Delta t$ oscillation is dissipated if $m = 0$ but amplified if $m = \pm 1$. This ambivalent behaviour is used to build a hybrid (TD) filter that removes the spurious oscillations and leaves the physical modes unchanged.

(3) As far as the $2\Delta t$ oscillation is concerned: for a dissipation effect equivalent to that of the M08 Laplacian filter, $K$ is twice smaller in the case of the TD filter.
(4) As far as the physical mode is concerned: the impact of the TD scheme is several orders smaller than that of the Laplacian filter. It behaves as a dissipation term when the scheme is based on (2.19) with $m = +1$ and as an amplification term when $m = -1$.

(5) A linear combination of these two TD schemes enforces the selective properties of the filter. But the stability of the so-called FD filter is conditioned to the use of an appropriate tuning coefficient, still not defined at this stage of the paper.

### 3. The oscillation equation

#### 3.1. Amplification factor

The physically based approach of the time filters properties presented in Section 2 is particularly attractive since it provided several powerful conclusions with very little developments. On the other hand, as the calculus were based on what can be considered as rather restrictive hypothesis (the fact that $u_0$ is maintained constant over the time of integration of the energy balance, or the fact that we considered small time step increments), further analysis appears useful. For convenience (i.e. in order to easily distinguish Section 2 and Section 3 in the following) Section 3 is referred to as “mathematical” approach (meanwhile Section 2 is referred to as “physical approach”).

In the following, we examine the properties of various filters through the analysis of the oscillation equation used in Williams (2009). The latter is given by

$$ \frac{dF}{dt} = i\omega F. \qquad (3.1) $$

The discrete formulation of (3.1) based on the LF time stepping scheme combined to the M08 Laplacien (hereafter LP) filter is:

$$ F(t + \Delta t) = F(t - \Delta t) + \frac{v}{2}(F(t) - 2F(t - \Delta t) + F(t - 2\Delta t)) + 2i\omega\Delta tF(t). \qquad (3.2) $$

In order to make comparisons with Williams's study easier, we adopted similar notations. We notably used (2.6) (e.g. $v = \frac{4K}{\Delta t}$) in order to make $v$ appear in (3.2).

In the case of the TD scheme based on (2.19) with $m = +1$, and once again using Williams (2009) notations, we have:

$$ F(t + \Delta t) = F(t - \Delta t) + \frac{v}{2}(F(t) - 3F(t - \Delta t) + 3F(t - 2\Delta t) - F(t - 3\Delta t)) + 2i\omega\Delta tF(t), \qquad (3.3) $$

when $m = -1$ in (2.19), the TD filter leads to:

$$ F(t + \Delta t) = \frac{F(t - \Delta t) + \frac{v}{2}(3F(t) - 3F(t - \Delta t) + F(t - 2\Delta t)) + 2i\omega\Delta tF(t)}{1 + \frac{v}{2}}. \qquad (3.4) $$

We note that (3.4) is an implicit scheme when (3.3) is explicit. In the following, (3.4) is referred to as TDI and (3.3) as TDE.

As in Williams (2009) we calculate the amplification factor $A = F(t + \Delta t)/F(t)$ associated to the different schemes considered in the present study. As far as (3.2) (LP scheme) is concerned, the latter is obtained from:

$$ A^3 + \left(-\frac{v}{2} - 2i\omega\Delta t\right)A^2 + \left(-1 + v\right)A - \frac{v}{2} = 0. \qquad (3.5) $$

The polynomial corresponding to the Robert-Asselin filter can be found in Williams (2009) (see Eq. (11) with $\alpha = 1$ in the Williams's paper), namely:

$$ A^2 - (v + i\omega\Delta t)A + v - 1 + vi\omega\Delta t = 0. \qquad (3.6) $$

In the case of the TDE (3.3) scheme, the polynomial is:

$$ A^4 + \left(-\frac{v}{2} - 2i\omega\Delta t\right)A^3 + \left(-1 + \frac{3v}{2}\right)A^2 - \frac{3v}{2}A + \frac{v}{2} = 0, \qquad (3.7) $$

and in the TDI (3.4) case we obtain:

$$ \left(1 + \frac{v}{2}\right)A^3 + \left(-\frac{3v}{2} - 2i\omega\Delta t\right)A^2 + \left(-1 + \frac{3v}{2}\right)A - \frac{v}{2} = 0. \qquad (3.8) $$

There are 3 roots in the LP polynomial (3.5), 2 roots in the Asselin polynomial (3.6), 4 roots in the TDE polynomial (3.7) and 3 roots in the TDI polynomial (3.8). They are numerically computed using the “roots finder” routine of a mathematical toolbox. Fig. 1 presents the root corresponding to the physical mode (e.g. $A = 1$ for $\omega\Delta t = 0$) and the different “numerical mode” roots. As far as the diffusion coefficient is concerned, we chose a commonly used value (i.e. $v = 0.1$) for the LP filter. Taking (2.34) into account, a slightly higher value is actually taken for the Robert-Asselin filter, i.e. $v = 0.1056$. As far as the TDE and the TDI filters are concerned, $v = 0.05$, since we have shown that their diffusive effect is comparable to that of a LP filter but with a twice smaller coefficient.

First of all, Fig. 1 shows that the use of any of the considered filters (Fig. 1(b) corresponding to the Leap Frog scheme with no filter) leads to lower the stability of the LF scheme. As a matter of fact, in all cases, one of the roots is outside of the unity circle for some high enough values of $\omega\Delta t$ (within the range $[0,1]$), when the roots obtained in the case of the LF scheme with no time filter remain on the unity circle (as long as $\omega\Delta t \le 1$). Being on the unity circle does not mean that the LF scheme is exact but only that the amplitude of the oscillation is correctly represented. Note that If we extent the oscillation equation case to the case of the horizontal momentum equations, this property actually becomes synonymous with energy conservation (Marsaleix et al., 2008). On the other hand the LF introduces a well-known phase error (Durran, 1991).

#### 3.2. Damping of the $2\Delta t$ oscillations

As previously mentioned, a particular well-known feature of the LF scheme is to enable the existence of the “$2\Delta t$ numerical mode”. Fig. 1(b) evidences this property with a “numerical mode” root remaining on the unity circle whatever the time step within the range $0 \le \omega\Delta t \le 1$. All the time filter schemes considered here lead to a substantial damping of the $2\Delta t$ oscillation provided, in some cases, that $\omega\Delta t$ remains lower than a critical value smaller than one. As far as the TDE and LP filters are concerned, the critical value for $\omega\Delta t$ is 0.802 and 0.951, respectively. These values depend on $v$. They would be smaller if a higher coefficient had been chosen.

#### 3.3. Robert-Asselin versus Laplacian filters

As far as the LP and Robert-Asselin filters are concerned, the stability of the two corresponding time stepping schemes are quite similar since the $\omega\Delta t$ critical value over which one of the two roots crosses the unity circle is almost the same in both cases. We simply note that the LF + RA scheme (critical $\omega\Delta t = 0.9485$) is slightly less stable than the LF + LP scheme (critical $\omega\Delta t = 0.951$). As long as the physical mode is concerned, the physically based analysis (Section 2) suggests that the dissipative effect is lower in the case of the Laplacian filter. As Fig. 1 is not really appropriate to evidence small differences between schemes, Fig. 2 gives an alternative representation of the amplification factor around $A \approx 1$. The conclusions of Section 2 are confirmed by Fig. 2, at least for a wide range of relevant time steps (e.g. $\omega\Delta t < 0.86$). This is not true anymore for
$\omega\Delta t > 0.86$ but we note that such large time steps are not really representative of realistic applications.

### 3.4. TDE filter

Fig. 2 confirms that the dissipative effect on the physical mode is much smaller with the TDE filter than with the Robert-Asselin and the Laplacian filters. This rather satisfying behaviour is however tempered by the fact that the stability of the TDE filter is lower than that of the Robert-Asselin and Laplacian filters. Fig. 1 indeed shows that a "numerical mode" root crossing the unity circle for a relatively small time step ($\omega\Delta t = 0.802$).

### 3.5. TDI filter

On the other hand, the numerical modes are well damped by the TDI scheme whatever the time step value. Unfortunately, as far as the physical mode is concerned, the behaviour of the TDI filter is questionable. As expected from the physically based approach, the magnitude of the amplification factor is indeed bigger than one, whatever the time step. As the latter remains very close to one for a wide range of time steps, it is not clear whether the TDI filter should be discarded or not. We actually may wonder whether other possible forms of dissipation in the model (depending on the numerical scheme, advection is for instance potentially diffusive) can compensate for the inherently unstable

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Amplification factor of the oscillation equation problem as a function of the time step. The physical mode is represented by a solid line and the numerical modes with dashed lines. Labels indicate the value of the dimensionless time step $\omega\Delta t$. Exact solution: a; Leap Frog with no filter: b; Asselin filter: c; Laplacian filter: d; TDI filter: e; TDE filter: f; FD filter: g.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 1 (continued)</figcaption></figure>

nature of the LF + TDI time stepping scheme. We however note that the behaviour of the LF + TDI depicted by Fig. 1(e) is very similar to that of Fig. 3(c) in Williams (2009) obtained with a LF scheme associated to a modified Robert-Asselin filter preserving a three-time-level mean state. This property of conservation of the so-called RAW filter is ensured provided that a tuning coefficient, conveniently introduced by the author, is equal to 1/2. Williams (2009) underlines the unconditionally unstable nature of the RAW filter in this special case and recommends the use of a higher tuning coefficient in order to recover a conditionally stable filter (Williams, 2011; Amezcua et al., 2011). In the light of Williams's results, we can reasonably expect the TDI filter to be inappropriate in realistic simulations. This point will be actually assessed in Section 4.

### 3.6. FD filter

Fig. 1 shows some symmetries between the TDI and TDE shortcomings, suggesting that the latter could be eliminated, at least partially, by a combination of the two filters. The FD filter presented in Section 2.6 actually corresponds to the following linear combination:

$$ \alpha \times \text{TDI} + (1 - \alpha) \times \text{TDE}, \quad (3.8) $$

with $0 \le \alpha \le 1$. For the particular case $\alpha = 0.5$, we may note that this new filter can be seen as a discrete equivalent of:

$$ \frac{\partial F}{\partial t} = - \frac{K \Delta t^2}{2} \frac{\partial^4 F}{\partial t^4}. \quad (3.9) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Amplification factor around A ≈ 1 as a function of the dimensionless time step $\omega\Delta t$ and of the time filter scheme.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Modulus of the amplification error, $|1 - A|$, as a function of the dimensionless time step $\omega\Delta t$ and of the time filter scheme. The slope of the curves indicates the $\Delta t$ dependency order of the errors (note the logarithmic axis). For guidance, straight lines of slope 2 and 4 are also drawn.</figcaption></figure>

**Table 1**

Critical stability number $\omega\Delta t$ as a function of $v$ and $\alpha$ (underlined values corresponding to the best setting of the tuning coefficient $\alpha$).

<table><thead><tr><th></th><th>v = 0.05</th><th>v = 0.1</th><th>v = 0.2</th></tr></thead><tbody><tr><td>α = 0.5</td><td>0.975</td><td>0.950</td><td>0.900</td></tr><tr><td>α = 0.52</td><td><u>0.976</u></td><td>0.953</td><td>0.907</td></tr><tr><td>α = 0.54</td><td>0.975</td><td><u>0.954</u></td><td>0.913</td></tr><tr><td>α = 0.57</td><td>0.969</td><td>0.951</td><td><u>0.916</u></td></tr></tbody></table>

We can be tempted to find an analogy with the bilaplacian scheme often used as a spatial filter (Marchesiello et al., 2009). On the other hand, this analogy is limited by a somehow unexpected negative sign at the RHS of (3.9) and also by the fact that $\alpha = 0.5$ is a priori not the best setting since the energy analysis of the TDE and TDI filters (Section 2) actually suggests that a better compromise should be found with a higher weight given to TDI in (3.8).

As far as the oscillation equation is concerned, the general formulation of the LF associated to this hybrid filter leads to:
$$F(t + \Delta t) = \frac{F(t - \Delta t) + \frac{v}{2}((1 + 2\alpha)F(t) - 3F(t - \Delta t) + (3 - 2\alpha)F(t - 2\Delta t) + (\alpha - 1)F(t - 3\Delta t)) + 2i\omega\Delta tF(t)}{1 + \frac{v}{2}\alpha} \quad (3.10)$$

<figure><img src="image_6.png" /><figcaption>Fig. 4. Phase error as a function of the dimensionless time step $\omega\Delta t$ and of the time filter scheme. (a) Absolute phase error (°). (b) Relative phase error. The slope of the curves indicates of the $\Delta t$ dependency order of the errors (note the logarithmic axis). For guidance, a straight line of slope 2 is also drawn.</figcaption></figure>

with the related polynomial:

$$(1 + \frac{v}{2}\alpha)A^4 + \left(-\frac{v(1 + 2\alpha)}{2} - 2i\omega\Delta t\right)A^3 + \left(-1 + \frac{3v}{2}\right)A^2 - \frac{v(3 - 2\alpha)}{2}A + \frac{v}{2}(1 - \alpha) = 0. \quad (3.11)$$

Fig. 1(g) shows the amplification factor of the physical mode and of the computational modes corresponding to (3.11). Fig. 1(g) is obtained with the same diffusion coefficient than that previously used for the TDE and TDI filters, i.e. $v = 0.05$. Clearly, the physical mode is less impacted when the FD filter is used (see also Fig. 2). Beside, the computational modes are efficiently dissipated. Note that several values have been tested for the tuning coefficient $\alpha$. The highest
critical stability number, $\omega\Delta t = 0.976$, is obtained with $\alpha$ close to 0.52. This result actually depends on the diffusion coefficient. Indeed, when $\nu = 0.1$ (respectively, 0.2), the highest critical number, $\omega\Delta t = 0.954$ (respectively, 0.916), is obtained for $\alpha$ close to 0.54 (respectively, 0.57). These optimal values, as expected smaller than the theoretical value mentioned in Section 2.6, are close to the straightforward value $\alpha = 0.5$ for which the FD stability remains quite satisfying (see Table 1). In this latter case, the critical stability number is indeed equal to 0.975 (respectively, 0.950 or 0.900) for $\nu = 0.05$ (respectively, 0.1 or 0.2), that is very close to the optimised setting. Consequently, $\alpha = 0.5$ can be seen as an acceptable, v-independent, default tuning of the FD scheme.

### 3.7. Partial conclusions and technical considerations

The analysis of the oscillation equation confirms the principal conclusions of the physically based approach presented in Section 2. It notably confirms that the TDE and TDI filters have a much lower impact on the physical mode. Meanwhile the TDE and TDI filters successfully suppress the spurious computational mode. However, in the case of the TDE filter, this is true provided that the time step remains lower than a relatively small (i.e. compared to other filters) critical time step. We also note that both physical and mathematical approaches predict an unstable behaviour of the TDI filter concerning the physical mode. The latter is however relatively slight so that we will check (using a realistic simulation presented in Section 5) if it can be balanced by other forms of dissipation in the model. Finally, as none of the TDE and TDI filters is completely satisfying, the FD scheme ends up being a much better alternative. Last, Fig. 3 confirms the $\Delta t^4$ dependency of the amplitude error in the case of the TD and FD filters, as well as the $\Delta t^2$ dependency of the LP and Robert-Asselin filters, provided that $\Delta t$ remains small enough, (accuracy becoming sensitive to higher orders contributions for large time steps).

The present paper is focussed on the attenuation of the signal amplitude, leaving the phase errors problem for a future study. It is however worth noting that the phase errors are not aggravated by the different filters proposed here as low-diffusive alternatives to the Robert-Asselin filter. Fig. 4a, giving the absolute phase error as a function of the time step (parameters setting as for Fig. 1), shows that the phase accuracy is very similar for all the filter schemes within the first third of the range of stability. Noticeable differences only appear for larger time steps, the Robert-Asselin filter producing the largest errors and the TDE filter the smallest ones. The relative phase errors (Fig. 4(b)) are $\Delta t^2$ dependent, whatever the time filter, except for $\omega\Delta t$ close to one where the accuracy is dominated by higher orders.

According to (3.2) and (3.4) the computation of $F(t + \Delta t)$ using the Laplacian or the TDI filters is based on $F(t)$, $F(t - \Delta t)$, $F(t - 2\Delta t)$. On the other hand, the use of the TDE or FD schemes requires the storage of an additional state, e.g. $F(t - 3\Delta t)$.

The memory storage induced by all these filters is thus heavier than that of the Robert-Asselin and RAW filters, for which the computation of $F(t + \Delta t)$ only requires the knowledge of $F(t)$ and $F(t - \Delta t)$ (Williams, 2009). In order to estimate the relative impact on the overall model memory, we now look at the results of the compilation of a wide spread open access model, namely the pom98.f code (Blumberg and Mellor, 1987) available at http://www.aos.princeton.edu/WWWPUBLIC/htdocs.pom. This code is first compiled in its original version and then recompiled with the additional arrays required by the aforementioned filters. Results are summarised in Table 2. As all the ocean models do not use a turbulent closure scheme based on two prognostic variables (the POM model actually uses the Mellor and Yamada 1982, scheme), Table 2 considers two cases; case 1 where all the state variables (e.g. including the turbulent variables) are concerned by

<caption>**Table 2**</caption>

Relative increase of the overall memory requirement of the pom98.f code depending on the time filter scheme. Case 1: all the states variables are concerned, including the Mellor and Yamada (1982) turbulent scheme. Case 2: the turbulent variables are ignored.

<table><thead><tr><th></th><th>LP, TDI</th><th>TDE, FD</th></tr></thead><tbody><tr><td>Case 1</td><td>1.15</td><td>1.3</td></tr><tr><td>Case 2</td><td>1.1</td><td>1.2</td></tr></tbody></table>

the time filter and Case 2 ignoring the turbulent variables. The Laplacian and TDI filters lead to an overall memory representing 1.15 (case 1) or 1.1 (case 2) times that of the original code. These values become 1.3 (case 1) or 1.2 (case 2) when the TDE and FD filters are concerned. Let us note that for a given three dimensional grid, several models can have quiet different memory requirements. Ezer et al. (2002) for instance reported that the ROMS model is 46/29 times heavier than POM. As a matter of fact, we can expect models based on low order numerical schemes to require less storage than high order models for which the implementation of the FD filter would imply a more modest increase of the overall memory.

As far as the computation time is concerned, the most sophisticated filter considered in this study (e.g. the FD filter in its conservative form presented in the Appendix A) represents 1.8 times the computation expense of the Robert-Asselin filter, when the basic formulation of the FD filter (3.10) requires only 1.3 times that of the Robert-Asselin filter. The time filters considered here are anyway simple compared to the numerical complexity generally found amongst the other terms of the momentum and scalars equations (advection, pressure gradient schemes, turbulence closure, etc.). As a consequence, the computation expense induced by the time filter only represents a small ratio of the overall cost. Again considering the case of the pom98.f (we run the test case provided with the pom98.f code), we notice that the Robert-Asselin filter requires roughly 5% of the overall computation time. This cost would certainly be smaller in models using more sophisticated methods than those implemented in the POM code (Marchesiello et al., 2009; Shchepetkin and McWilliams, 2003; Chu and Fan, 1997).

## 4. Numerical experiments

The different filters are now tested using numerical simulations. The first test is based on the oscillation equation. Then we consider the case of internal waves propagating in a 2 dimensional vertical plan. Last we consider the case of a 3D realistic simulation. These tests are made with the 3D circulation ocean model described in Marsaleix et al. (2008, 2009).

### 4.1. Oscillation equation

The fully discrete form of the Oscillation equation analyzed in Section 3 is obtained by solving the equations for X and Y, the real and imaginary components of $F = X + iY$ (details in Williams (2009)), that is:

$$ 
\begin{aligned}
\frac{\partial X}{\partial t} &= -\omega Y, \\
\frac{\partial Y}{\partial t} &= +\omega X.
\end{aligned} 
\qquad (4.1) $$

The different discrete forms related to (4.1) depend on the time filter scheme. They are simply deduced from the discrete expressions given in Section 3 and Section 4, then replacing F by X or Y and $i\omega F$ by $-\omega Y$ or $+\omega X$.
The time step is such that $\omega\Delta t = \frac{2\pi}{100}$. If we were considering a semi-diurnal signal, $\Delta t$ would be of this order of 400 s, which could be considered as a possible realistic value for a circulation ocean model. On the other hand, the diffusion coefficient is set to $\nu = 0.1$ in order to remain consistent with the previous sections. The initial discrete state is given by the exact solution. Fig. 5 shows the different discrete solutions and the exact solution after 200 oscillations. The computational mode does not appear in the solution obtained with the LF alone. This is an expected result since the LF scheme neither amplifies nor damps the computational mode. The amplitude of the oscillation is perfectly conserved by the LF scheme alone. However, Fig. 5 evidences the well-known lack of accuracy of the LF scheme as long as the phase is concerned. Indeed, after 200 oscillations, the phase error is about 48° degrees, which is not negligible. On the other hand, as the leading order terms of the LF scheme error depend on $\Delta t^2$, decreasing the time step is an efficient way to reduce the errors. The dissipative impact of the low order time filter on the physical mode is obvious. The Laplacian and Robert-Asselin schemes almost suppressed roughly 90% of the physical mode. Between these two filters, the Laplacian scheme eventually appears as a “less worse” option. The benefit of the other filters is obvious. The LF solution is quasi unchanged by the use of the FD filter. The solutions obtained with TDE and the TDI filters are not that good but remain close to that of the LF scheme alone. We note that the TDI filter tends to aggravate the phase error of the LF scheme. On the other hand, the TDE scheme tends to diminish it. The reduction of the phase error is anyway too small to see in the TDE filter a way to overcome the lack of accuracy of the LF scheme.

## 4.2. Internal Gravity Waves in a 2DV plan

The following tests are based on simulations of internal waves propagating in an idealised, non rotating, flat, continuously stratified ocean. They are performed with the free surface sigma coordinate SYMPHONIE ocean model described in Marsaleix et al. (2006, 2008, 2009a). Applications of this model to realistic cases of gravity waves in the North-East Atlantic Ocean can be found

in Pairaud et al. (2008, 2010). The following simulations are inspired from Floor et al. (2011), and Marsaleix et al. (2009b). The horizontal resolution is 1 km. The vertical grid consists of 30 regularly spaced sigma levels. The density field is horizontally homogeneous and linear in $z$: $d\rho/dz = 10^{-3}$ kg m$^{-4}$. The background reference density is $\rho_0 = 1000$ kg m$^{-3}$. The corresponding Brunt-Väisälä Frequency is $N = (-g\rho_0^{-1}d\rho/dz)^{0.5} = 3.132 \times 10^{-3}$ s$^{-1}$. The constant sea-floor depth is $h = 1000$ m. The two dimensional vertical grid is 2000 km long. The first baroclinic mode of the analytical solution of the normal mode problem described in Gill (1982) is used as boundary conditions for the velocity at the incoming ($x=0$) lateral open boundary. Practically the forcing velocity is set to:

$$u = u_0 \cos(\pi z/h) \sin(\omega t - kx) \quad x = 0. \qquad (4.2)$$

The wave amplitude is $u_0 = 0.001$ ms$^{-1}$. This choice of small amplitude intends to limit the dissipation related to the turbulence closure and the bottom friction. Besides, centred, few diffusive, advection schemes are purposely used, so that dissipation is essentially caused by the different filters of the time stepping scheme. The period of waves is 12 h. The theoretical phase speed is $c = Nh/\pi \approx 0.996$ ms$^{-1}$ and the wave length is $l \approx 43$ km. The simulation is stopped before the waves reach the outgoing boundary in order to avoid any problem of spurious reflected waves. The theoretical maximum time step depends on the horizontal resolution and the

<caption>**Table 3**</caption>

First column: diffusion coefficient. Second column: approximated (since empirically determined) maximum allowed time step as a function of the time filter scheme. Third column: the latter is divided by the theoretical maximum time step permitted by the Leap Frog scheme with no time filter.

<table><thead><tr><th>Filter</th><th>ν</th><th>Δt<sub>max</sub></th><th>Δt<sub>max</sub>/Δt<sub>max</sub><sup>theory</sup></th></tr></thead><tbody><tr><td>Robert-Asselin</td><td>0.1058</td><td>476s</td><td>0.948</td></tr><tr><td>Laplacian</td><td>0.1</td><td>477s</td><td>0.950</td></tr><tr><td>TDE</td><td>0.05</td><td>406s</td><td>0.809</td></tr><tr><td>TDI</td><td>0.05</td><td></td><td></td></tr><tr><td>FD</td><td>0.05</td><td>488s</td><td>0.972</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Numerical solution of the oscillation equation (X in Eq. (4.1)) as a function of time (in number of oscillation periods) after about 200 oscillations.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Internal waves propagating in a two dimensional vertical plan. Surface velocity after 20 days around $x = 1000$ km. The time step is $\Delta t = 400$ s. The diffusion coefficient, $v$, is given in Table 3.</figcaption></figure>

phase speed. In the case of the Leap-Frog scheme with no time filter the latter is $\Delta t_{\text{max}}^{\text{theory}} = 0.5\Delta x/c \approx 502$ s (Nycander and Döös, 2003).

Table 3 gives the maximum permitted time step depending on the time filter scheme used in the simulation, empirically obtained by increasing the time step until the model blows up.

As far as the TDI filter is concerned, the simulation blows up whatever the value of the time step. As previously mentioned, we tried to control the expected growth of the physical mode through other dissipative processes in the model. We notably replaced the centred advection scheme of the momentum and tracer equations by an upwind scheme, since the latter is known for its strong diffusive properties. Our attempt was unfortunately unsuccessful. Actually, the diffusive effect of the upwind scheme preferentially acts at the grid mesh scale and is consequently rather powerless to prevent the growth of the physical mode. As a matter of fact, its wave length (43 km) is much larger than the grid mesh.

The maximum time step permitted by the TDE scheme is 406 s that is 80, 9% of the theoretical maximum time step. The stability of the Laplacian filter is a little bit better than that of the Robert-Asselin filter. The FD filter is the most stable scheme amongst the considered filters. Last, we note that Table 3 is consistent with the maximum time steps allowed by the oscillation equation for the different time filter schemes (see Fig. 1), as expected from previous theoretical studies (see for instance Eq. (7) in Shchepetkin and McWilliams, 2005).

Fig. 6 shows the surface velocity after 20 days, at 1000 km from the wave source (e.g. approximately 23 wave lengths). All the simulations use the same time step ($\Delta t = 400$ s). Fig. 6 confirms the comments made on Fig. 5. The LP and RA filters have a strong dissipative impact on the physical mode, the RA filter being the worse of the two filters. The amplitude of the physical mode is much better preserved with the FD and TDE filters, both filters giving very similar results. The phase lag reported in the previous section between the TDE and FD filter is too small to be detected on Fig. 6 but can be evidenced by a more careful examination (not shown). The latter is about 0.5° in favour of the TDE scheme. This difference is anyway far from being enough to compensate for the phase error. The latter is about 14°.

## 5. Conclusion

Four alternative schemes have been analysed with the aim of finding an efficient solution to the well-known shortcoming of the Robert-Asselin filter, namely the numerical dissipation of the physical mode.

Provided that their respective diffusion coefficients are tuned in order to obtain an equivalent damping effect on the computational mode, the Laplacian filter (Marsaleix et al., 2008) seems to be a better option than the Robert-Asselin filter. Indeed, the physical mode is less impacted and the stability of the time stepping scheme is stronger in the case of the Laplacian filter. Nevertheless, the behaviour of the two filters remains rather similar and for a sufficiently small value of the diffusion coefficient they can even be considered equivalent. The Laplacian filter is consequently not a really convincing alternative to the Robert-Asselin filter.

The analysis of the role of the time filter scheme in the discrete kinetic energy balance suggested the building of two higher order time filters. One is an explicit scheme based on 4 points (TDE) and the other one is a 4 points implicit scheme (TDI). Unfortunately, none of them is completely satisfying. Indeed, as long as the computational mode is concerned, the stability of the TDE scheme is sensibly lower than that of the other filters. On the other hand, the physical mode appears unstable, whatever the time step, in the case of the TDI filter.

A 5 points filter finally appeared as the best option amongst the schemes considered in the present study. Indeed, the FD filter let the physical mode almost unchanged and enables a larger time step compared to the other filters. On the other hand, the FD filter needs more computing memory storage than the RAW filter, a low diffusive alternative to the Asselin filter proposed by Williams (2009).

Interestingly, the TDE filter is the only scheme (amongst the schemes that have been considered here) able to reduce the phase error of the Leap-Frog scheme. The improvements (noticeable in our numerical tests) are however too small to really overcome the advantages of the FD filter. On the other hand the well-known problem of phase errors of the LF scheme remains unsolved by the FD filter.
Acknowledgements

This study was supported by the French programs ANR COM- ODO and GMMC. We used the 3D ocean circulation SYMPHONIE developed by the SIROCCO oceanography group of the Toulouse University and CNRS. We thank one anonymous reviewer and Paul Williams for their helpful comments as well as the Laboratoire d'Aérologie Computer Team, Serge Prieur, Laurent Cabanas, Jérémy Leclercq, Didier Gazen and Juan Escobar for their support.

Appendix A. The case of moving vertical levels

In sigma coordinate free surface ocean models, the depth of ver- tical levels is not constant in time. Considerations on the conserva- tion properties of the time stepping and advection schemes generally lead to a particular formulation of the model equations, practically:

$$
\frac{\partial \Delta z F}{\partial t} = \frac{\partial \Delta z \psi}{\partial t} + \Delta z \text{RHS}, \qquad (A1)
$$

where Δz is the thickness of the vertical layers. We note that the comments about the Green-Ostrogradski theorem made in Section 2.1 still apply to (A1). The discrete expression corresponding to (A1) is:

$$
\frac{\Delta z^{n+1} F^{n+1} - \Delta z^{n-1} F^{n-1}}{2\Delta t} = \frac{\Delta z^{n+1/2} \psi^{n+1/2} - \Delta z^{n-1/2} \psi^{n-1/2}}{\Delta t} + \Delta z^n \text{RHS}^n, \quad (A2)
$$

where Δz<sup>n+1/2</sup> = (Δz<sup>n</sup> + Δz<sup>n+1</sup>)/2. As Δz<sup>n+1</sup> is already determined at this stage, F<sup>n+1</sup> is simply obtained from:

$$
F^{n+1} = \left( \Delta z^{n-1} F^{n-1} + 2\Delta z^{n+1/2} \psi^{n+1/2} - 2\Delta z^{n-1/2} \psi^{n-1/2} + 2\Delta t \Delta z^n \text{RHS}^n \right) / \Delta z^{n+1}. \quad (A3)
$$

As far as the FD filter is concerned:

$$
\begin{align}
\psi^{n+1/2} &= \frac{K}{\Delta t} [F^n - F^{n-1} - \alpha(F^{n+1} - F^n) - (1-\alpha)(F^{n-1} - F^{n-2})] \nonumber \\
&= \frac{\nu}{4} [-\alpha F^{n+1} + (1+\alpha)F^n - (2-\alpha)F^{n-1} + (1-\alpha)F^{n-2}]. \tag{A4}
\end{align}
$$

Using (A4), (A1) leads to:

$$
F^{n+1} = \left\{ \Delta z^{n-1} F^{n-1} + \Delta z^{n+1/2} \frac{\nu}{2} \left[ (1+\alpha)F^n - (2-\alpha)F^{n-1} + (1-\alpha)F^{n-2} \right] \right. \\
\left. - \Delta z^{n-1/2} \frac{\nu}{2} \left[ -\alpha F^n + (1+\alpha)F^{n-1} - (2-\alpha)F^{n-2} + (1-\alpha)F^{n-3} \right] \right. \\
\left. + 2\Delta t \Delta z^n \text{RHS}^n \right\} / \left\{ \Delta z^{n+1} + \Delta z^{n+1/2} \frac{\alpha\nu}{2} \right\}. \quad (A5)
$$

The TDI and TDE filters are simply obtained by using (A5) with α = 1 or α = 0 respectively.

References

Amezcua, Javier, Kalnay, Eugenia, Williams, Paul D., 2011. The effects of the RAW Filter on the climatology and forecast skill of the SPEEDY Model. Mon. Weather Rev. 139, 608–619. doi:10.1175/2010MWR3530.1.

Asselin, R., 1972. Frequency filter for time integrations. Mon. Weather Rev. 105 (6), 487–490.

Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal circulation model. In: Heaps, N. (Ed.), Three-Dimensional Coastal Ocean Models, Coastal Estuarine Science, vol. 4. American Geophysical Union, pp. 1–16.

Chu, P.C., Fan, C., 1997. Sixth-order difference scheme for sigma coordinate ocean models. J. Phys. Oceanogr. 27, 2064–2071.

Durran, D.R., 1991. The third-order Adams–Bashforth method: an attractive alternative to leapfrog time differencing. Mon. Weather Rev. 119, 702–720.

Ezer, T., Arango, H., Shchepetkin, A.F., 2002. Developments in terrain-following ocean models: intercomparisons of numerical aspects. Ocean Modell. 4, 249–267.

Floor, J.W., Auclair, F., Marsaleix, P., 2011. Energy transfers in internal tide generation, propagation and dissipation in the deep ocean. Ocean Modell. 38, 22–40 <http://dx.doi.org/10.1016/j.ocemod.2011.01.009>.

Gill, A.E., 1982. Atmosphere Ocean Dynamics. Academic Press, p. 662.

Leclair, M., Madec, G., 2009. A conservative leapfrog time stepping method. Ocean Modell. 30, 88–94.

Leclair, M., Madec, G., 2011. ž-Coordinate, an arbitrary Lagrangian–Eulerian coordinate separating high and low frequency motions. Ocean Modell. 37, 139–152.

Marchesiello, P., Debreu, L., Couvelard, X., 2009. Spurious diapycnal mixing in terrain-following coordinate models: the problem and a solution. Ocean Modell. 26, 156–169.

Marsaleix, P., Auclair, F., Estournel, C., 2006. Considerations on open boundary conditions for regional and coastal ocean models. J. Atmos. Ocean. Technol. 23, 1604–1613 <http://dx.doi.org/10.1175/JTECH1930.1>.

Marsaleix, P., Auclair, F., Floor, J.W., Herrmann, M.J., Estournel, C., Pairaud, I., Uleses, C., 2008. Energy conservation issues in sigma-coordinate free-surface ocean models. Ocean Modell. 20, 61–89 <http://dx.doi.org/10.1016/
j.ocemod.2007.07.005>.

Marsaleix, P., Auclair, F., Estournel, C., 2009a. Low-order pressure gradient schemes in sigma coordinate models: the seamount test revisited. Ocean Modell. 30, 169–177 <http://dx.doi.org/10.1016/j.ocemod.2009.06.011>.

Marsaleix, P., Uleses, C., Pairaud, I., Herrmann, M.J., Floor, J.W., Estournel, C., Auclair, F., 2009b. Open boundary conditions for internal gravity wave modelling using polarization relations. Ocean Modell. 29, 27–42 <http://dx.doi.org/10.1016/
j.ocemod.2009.02.010>.

Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys. 20 (4), 851–875.

Mesinger, F., Arakawa, A., 1976. Numerical Methods Used in Atmospheric Models. GARP Publication Series No. 17.

Nycander, J., Döös, K., 2003. Open boundary conditions for barotropic waves. J. Geophys. Res. 108 (C5), 3168. doi:10.1029/2002JC001529.

Pairaud, I.L., Lyard, F., Auclair, F., Letellier, T., Marsaleix, P., 2008. Dynamics of the semi-diurnal and quarter-diurnal internal tides in the Bay of Biscay. Part 1: barotropic tides. Cont. Shelf Res. 28, 1294–1315. doi:10.1016/j.csr.2008.03.004.

Pairaud, I.L., Auclair, F., Marsaleix, P., Lyard, F., Pichon, A., 2010. Dynamics of the semi-diurnal and quarter-diurnal internal tides in the Bay of Biscay. Part 2: baroclinic tides. Cont. Shelf Res. 30, 253–269 <http://dx.doi.org/10.1016/
j.csr.2009.10.008>.

Robert, A.J., 1966. The integration of a low order spectral form of the primitive meteorological equations. J. Meteorol. Soc. Jpn. 44 (5), 237–245.

Shchepetkin, A.F., McWilliams, J.C., 2003. A method for computing horizontal pressure-gradient force in an oceanic model with non-aligned vertical coordinate. J. Geophys. Res. 108, 3090, 35.1–35.34.

Shchepetkin, A.F., McWilliams, J.C., 2005. The regional oceanic modeling system (ROMS): a split explicit, free surface, topography following coordinate oceanic model. Ocean Modell. 9, 347–404.

Williams, P.D., 2009. A proposed modification to the Robert–Asselin time filter. Mon. Weather Rev. 137, 2538–2546.

Williams, P.D., 2011. The RAW filter: an Improvement to the Robert–Asselin filter in semi-implicit integrations. Mon. Weather Rev. 139, 1996–2007. doi:10.1175/
2010MWR3601.1.