The aim of this paper is to investigate the feasibility of well-posed lateral boundary conditions in a Fourier spectral semi-implicit semi-Lagrangian one-dimensional model. Two aspects are analyzed: (i) the complication of designing well-posed boundary conditions for a spectral semi-implicit scheme and (ii) the implications of such a lateral boundary treatment for the semi-Lagrangian trajectory computations at the lateral boundaries.
Straightforwardly imposing boundary conditions in the gridpoint-explicit part of the semi-implicit time-marching scheme leads to numerical instabilities for time steps that are relevant in today’s numerical weather prediction applications. It is shown that an iterative scheme is capable of curing these instabilities. This new iterative boundary treatment has been tested in the framework of the one-dimensional shallow-water equations leading to a significant improvement in terms of stability.
As far as the semi-Lagrangian part of the time scheme is concerned, the use of a trajectory truncation scheme has been found to be stable in experimental tests, even for large values of the advective Courant number. It is also demonstrated that a well-posed buffer zone can be successfully applied in this spectral context. A promising (but not easily implemented) alternative to these three above-referenced schemes has been tested and is also presented here.
In limited-area numerical weather prediction (NWP) models, the lateral borders of the domain are not physical boundaries for the flow, but artificial boundaries constructed in the aim that the solution of the governing equations over this limited area remain consistent with the one of the global atmosphere. This results in a so-called initial boundary value problem. To solve the equations describing the evolution of the atmosphere in that context, some information has to be supplied at these artificial boundaries throughout the course of the integration time. As shown by Oliger and Sundström (1978), a field must be supplied at the lateral boundary for each inward-pointing characteristic. This implies that only some subset of all the prognostic variables should be imposed there. If the number of imposed boundary fields has been correctly chosen in the above sense, the initial boundary value problem is said to be “well posed.”
The construction of such a well-posed lateral boundary condition (LBC) in NWP models is a delicate task. Nevertheless, some well-posed LBC’s treaments have been recently designed with in mind meteorological forecast models, but exclusively for discretization schemes based on local pointwise horizontal representations of the fields. Namely, for the High-Resolution Limited-Area Model (HIRLAM) model, McDonald (2000, 2002, 2003, 2005, 2006) has shown that it was possible to build various types of well-posed boundary treatments for semi-implicit (SI) semi-Lagrangian schemes (SL) based on a finite-difference gridpoint discretization with a staggered Arakawa C grid. Still in pointwise approaches, Lie (2001) has demonstrated the feasibility of well-posed boundaries with a semi-implicit semi-Lagrangian (SISL) model using an horizontal finite-element discretization. However, finite-difference and finite-element methods are not the only horizontal discretization schemes used in NWP models. Space discretization methods based on Fourier spectral horizontal representation of the fields are also commonly used in limited-area models, but until now the question of the feasibility of well-posed LBCs within this spectral approach has not yet been investigated, especially in conjunction with SISL algorithms.
In the NWP community, Fourier spectral limited-area models are common practice, as for example, the National Centers for Environmental Prediction (NCEP) model (Juang et al. 1997), the Aire Limitée Adaptation Dynamique Développement International (ALADIN) model (Bubnová et al. 1995), and the so-called spectral version of the HIRLAM model (Haugen and Machenhauer 1993). But, because of the contradiction between a local specification of LBCs and the global character of the Fourier spectral discretization, designing stable and accurate well-posed boundary conditions is a daunting task. In Fourier spectral limited-area models, one of the most difficult questions is to combine the periodicity constraint for the Fourier expansion with well-posed time-dependent LBCs. For this reason, and because of the success of the engineering solution of Davies (1976), all spectral limited-area NWP models have been so far formulated with a flow-relaxation scheme, as in the NCEP model (Juang and Kanamitsu 1994) and in the ALADIN model (Radnóti 1995). This scheme consists in overspecifying the boundaries and in damping the resulting noise in a relaxation buffer zone. This gives stable forecasts and leads to a very easy implementation. Despite this, the relaxation scheme has however some unavoidable weaknesses that might jeopardize the benefit that could be expected from high-resolution modeling and from sophisticated physical parameterizations (see Warner et al. 1997). Additional reasons to look for a well-posed lateral boundary strategy have been hightlighted by McDonald (1999), and will not be detailed here.
The aim of this paper is to examine if a well-posed formulation of LBCs is possible (and viable) in the framework of the spectral technique, since this is the method currently used in the ALADIN model and that might not be questioned in the short term.
In the finite-difference or in the finite-element SISL frameworks, the key point for designing well-posed LBCs consists in locally altering the gridpoint implicit operator induced by the semi-implicit scheme, with the boundary conditions. This leads to a space-dependent closed system that can be iteratively solved at each time step, thanks to efficient gridpoint solvers. For the Fourier spectral method to be attractive, the so-called implicit operator must remain horizontally homogeneous in such a way that it is readily inverted in spectral space. Hence, the above well-posed LBC strategy, implying a local modification of the implicit operator, is not suited for spectral models. The main issue is to find how well-posed LBCs can be incorporated into the Fourier spectral SISL scheme in a stable and accurate manner without jeopardizing the attractive aspects of the Fourier spectral method.
This paper explores these issues in the framework of a one-dimensional linearized shallow-water spectral model, by restricting the specification of LBCs to the gridpoint space. It is first shown that imposing the LBCs in a purely explicit way cannot lead to a sufficiently stable behavior for operational applications. Then an alternative method is examined: the LBCs are introduced in the explicit gridpoint part of the computations, but they are treated in an approximate implicit way through an iterative procedure, in the same spirit as the iterative centered implicit (ICI) schemes proposed in Bénard (2003).
Moreover, McDonald (2000) identified accuracy problems linked to semi-Lagrangian trajectories when the departure lies outside the domain. He proposed three solutions that are also examined in the present paper. However, in the context of McDonald’s study as well as in the present context, these proposals turn out to be either inaccurate in certain circumstances or too costly for future NWP applications. An alternative solution is introduced here, based on the substepping idea of Termonia and Voitus (2008), which circumvents these drawbacks.
This paper is organized as follows. The one dimensional two-time level (2TL) SISL spectral shallow-water model is presented in section 2. In section 3, we examine the feasibility of well-posed LBCs by imposing the characteristics boundary conditions, first with an explicit and then with an iterative implicit treatment. In section 4, the problem of the semi-Lagrangian trajectories originating from outside the domain is examined. Some numerical tests are presented in section 5. Finally, concluding remarks and a discussion are given in section 6.
2. The one-dimensional shallow-water model
Let us consider the linearized shallow-water equations for a one-dimensional flow on an f plane along the x direction in the absence of orography. Assuming that the basic state denoted by (u, , ) is constant and stationary, the evolution of small perturbations u, υ, and Φ around this basic state is examined, where Φ stands for the logarithm of the geopotential field gz; where z is the free surface height; and u and υ are, respectively, the x and y components of the horizontal wind. As a result, the linearized system is given by
where the basic-state advection velocity u, Coriolis parameter f, and geopotential height = are constants. The characteristic variables (or Riemann invariant variables) are given by υ, p = u + cΦ, and q = u − cΦ, with c = denoting the gravity wave phase velocity. These variables are, respectively, associated with the characteristic velocities u, (u + c), and (u − c). The wave solutions to (1)–(3) can be written as ψ = ψ̂ exp[î(kx − ωt)]. Substituting this general wave solution into (1)–(3) implies that ω admits three solutions: a slow advective solution ωa = ku, which corresponds to the geostrophic balance state (u = 0, and fυ = c2∂Φ/∂x), and two fast solutions, which correspond to the inertia–gravity waves adjustment modes with ω± = k(u ± ck), and ck =
In this paper, a spectral model version of the shallow-water equations in (1)–(3) with a structure following the proposal of Haugen and Machenhauer (1993) is formulated over the horizontal limited-area domain 0 ≤ x ≤ L (hereafter the “C zone”). This model is discretized on a collocation grid with the index i = (0, iL), where i = 0 represents the left artificial boundary of the physical domain and i = iL is the right artificial boundary (i.e., L = iLΔx). The C zone is then extended by a meteorologically meaningless extension zone (hereafter the “E zone”), which is labeled by i = (iL + 1, N − 1), making the whole domain (C + E) fully periodic (i.e., the value of all fields at points x = NΔx and x = 0 coincide). The length of the (C + E) area is denoted Lx (i.e., Lx = NΔx). In the E zone, all fields are supplemented by a periodic extension by means of a spline interpolation before applying the FFT transforms to compute the horizontal derivative terms in spectral space. The cubic splines used in this paper are exactly the same as the existing one in the ALADIN model. Alternative spectral model formulations with the Fourier extension, such as the one proposed by Boyd (2005) may be adopted, but they will not be investigated here: the work presented herein is restricted to the extension method used in the ALADIN model.
This spectral model version is used to solve the shallow-water system throughout the limited-area domain C zone, by supplying time-dependent boundary conditions only at the lateral boundaries of the C zone. The purpose in then to design the LBCs in such a way that the solution over this C zone is consistent with the exact analytical solution provided over the unbounded domain and considered as the “truth.” The evolution of the pronostic fields inside the E zone is not relevant, this zone is there only to allow the application of the Fourier spectral method.
A two-time-level semi-Lagrangian semi-implicit (2-TL SISL) discretization of (1)–(3) is performed. The set of shallow-water equations being linear, there are no additional sources of instabilities due to nonlinear residuals terms, the 2-TL SISL scheme is consequently unconditionnally stable and O(Δt2) accurate. The discretized equations of (1)–(3) can be written as
where the asterisk subscript denotes interpolations at the departure point of the SL trajectory x* = xi − uΔt, with xi = iΔx. The superscripts plus sign and 0 correspond to time levels t + Δt and t, respectively. X = (u, υ, Φ)T represents the state vector. Assuming that the initial state X0 is made periodic using the extension procedure, then transformed in spectral space, the detailed model structure in absence of any time-dependent lateral boundary conditions is organized as follows: the derivative of the (periodic) fields ∂xX0 are computed in spectral space after some suitable spectral truncation. An inverse Fourier transform is applied to X0 and to ∂xX0 in order to transfer these fields to gridpoint space. The right-hand side (rhs) explicit terms A0 = (AΦ, Au, Aυ)T are computed in gridpoint space at time t through (7)–(9), and then evaluated through cubic interpolation at the departure points to get the (periodic) fields R = (Ru, Rυ, RΦ)T over the whole (C + E) domain, which are then transformed into spectral space.
The left-hand side (lhs) implicit terms of (4)–(6) correspond to the implicit part of the SI scheme. This system can be put into the symbolic matrix form [𝗜 − (Δt/2)𝗟] X+ = R, where the definition of the linear differential operator [𝗜 − (Δt/2)𝗟] is directly derived from (4)–(6). This operator is inverted by solving a Helmholtz equation for the u field, which is obtained by substituting Φ by (4) in (5), and υ by (6) in (5):
In spectral space, the Helmholtz equation in (10) is solved by a simple algebraic division for each spectral mode:
where û+K, R̂uK, R̂υK, and R̂ΦK are the spectral coefficients associated to the wavenumber K for u+, Ru, Rυ, and RΦ, respectively, and 2πK/Lx the multiplicative factor of the derivative operator. Then, Φ̂+K and υ̂+K are easily obtained through back substitution in the spectral counterpart of (4) and (6). It is important to note that this solution in spectral space is allowed provided that the coefficients of the Helmholtz operator are spatially homogeneous.
The challenge is to insert well-posed boundary conditions into this SISL model structure. Sundström and Elvius (1979) have suggested to externally supply the well-posed incoming characteristic variables at the boundaries into the implicit part of the SI scheme, which leads to a local modification of the Helmholtz gridpoint solver at the boundaries. This idea has been successfully explored by McDonald (2000). In this gridpoint representation, it is always possible to impose boundary conditions anywhere locally. On the contrary, in spectral representation, the value at a single point cannot be modified by changing the coefficients of the spectral representation of the fields. This holds, a fortiori, in the spectral form (11) of the Helmholtz equation. Hence, local modifications can only be performed in gridpoint space computation of the spectral model structure. Besides, since the implicit part is evaluated in spectral space, it is impossible to impose a well-posed boundary condition in a purely implicit manner.
3. Well-posed boundary treatment
Oliger and Sundström (1978) have pointed out that well posedness implies that for each inward-pointing characteristic, one field must be supplied at the considered lateral boundary. In the context of the shallow-water system examined here, this means that if |u| < then two fields have to be supplied at the inflow lateral boundary and only one at the outflow. Inversely, if |u| > then three fields have to be supplied at the inflow boundary and none at the outflow. In this paper, the case |u| < will be deeply investigated. The other case can be treated in a similar way with some minor changes, but will not be examined here.
a. Explicit well-posed boundary conditions
Assuming that u > 0, without loss of generality, the left physical boundary (x = 0) becomes the inflow boundary and the right physical boundary (x = L = iLΔx) is the outflow boundary. Let us consider the case u < . As outlined above, two fields have to be supplied at x = 0 and only one at x = L. Thus, we impose p0,h0, υ0,h0, and q0,hiL at boundaries, where the superscript “h” denotes the external boundary data supplied by the “host” model following the terminology in McDonald (2002). In pratice, they are imposed through the Φ field and the υ field at the beginning of the gridpoint computations in such a way that Φ00 = (1/c)(p0,h0 − u00), υ00 = υ0,h0, and Φ0iL = −(1/c)(q0,hiL − u0iL).
The derivative terms involved in the rhs terms at boundaries i = 0 and i = iL are evaluated using a first-order off-centered finite-differente approximation in order to avoid the use of grid points from outside the C zone. For stability reasons, we used second-order centered finite-difference approximations to estimate the derivative terms at the immediately adjacent grid points from the boundaries i = 1 and i = iL − 1 (not to do this was found to lead to instabilities in our tests). Consequently, near the boundaries, the explicit rhs terms in (7)–(9) are evaluated by
where the tilde denotes a field modified by LBC specification. The same finite-difference approximations are employed to compute ∂u/∂x at the boundaries i = 0, 1, iL − 1, and iL.. In the interior of the physical domain (i.e., i = 2, . . . , iL − 2), the derivatives are computed spectrally.
Concerning SL trajectory computations at inflow boundaries, the so-called trajectory truncation method, as defined in McDonald (2000), is applied: when the departure point is outside the physical domain (i.e., x* < 0 or x* > iLΔx), the trajectory is simply truncated at the boundary of the physical domain C zone. Otherwise, when some points of the stencil of the cubic operator lay outside the domain [i.e., 0 < x* < Δx or (iL − 1) Δx < x* < iLΔx] the interpolations are performed only quadratically, see McDonald (2000) for more details. The resulting interpolated rhs terms are then periodically extended into the E zone before applying the FFT direct transform procedure and solving the implicit problem in spectral space.
To test the stability and the consistency of this explicit boundary option, we use the fastest wave solution corresponding to the frequency ω+. For this wave solution and for the wavenumber k = 2π/Lx, the variables Φ, u, and υ are given by
where a+ and λ+ are constants. This wave solution is periodic over the whole domain (C + E). The initial state is given by appropriate Φ(x, 0), u(x, 0), and υ(x, 0) fields. The initial and boundary conditions are chosen in such a way that this wave moves in the positive x direction with the velocity u + ck. The following parameters are used; Lx = 1200 km, L = 1000 km, Δx = 10 km, a+ = 10, λ+ = 0, f = 10−4 s−1, and ck = 300.4 m s−1. In the tests below, the length of the integration (T) is always chosen in such a way that the traveling wave exactly covers the distance Lx in the analytical solution. Thus, the initial and the final state should overlap, making any error “jump out” immediately. The characteristic boundary conditions are given by the analytical solution: p(0, t) = u(0, t) + cΦ(0, t), υ(0, t), and q(L, t) = u(L, t) − cΦ(L, t).
In a first test, the settings are chosen as follows: Δt = 49.94 s, c = 300 m s−1, u = 100 m s−1, so that the wave Courant number cΔt/Δx is 1.5, and the advective Courant number α = uΔt/Δx is 0.4994. This choice avoids severe SL trajectory truncation. The LBC variables are taken from the exact solution in (16)–(18) and supplied at all time steps. This is an idealization with respect to operational models where they come from a model of lower resolution and are interpolated in time. Complications due to the temporal interpolation (see, e.g., Warner et al. 1997; Termonia 2003) will not be treated here. Termonia (2004) proposed a method to monitor the quality of the interpolation operationally. A detailed study within the context of this work lies beyond the scope of this paper.
Figure 1 represents the fields u ± ckΦ and υ after 50 min of integration over the physical domain lying between 0 and 1000 km and the E zone between 1000 and 1200 km. As can be seen in this figure, these fields behave as expected at the boundaries and the numerical solution over of the whole domain (C + E) remains stable and accurate. In another test, the settings are now chosen as Δt = 100 s and u = 49.94 m s−1 in such a way that α remains unmodified, but now c(Δt/Δx) = 3. This test leads to a strongly unstable behavior (not shown). Further tests with various other settings indicate that this purely explicit boundary scheme becomes unstable when the wave Courant number cΔt/Δx is bigger than 2. Then a growing mode appears first at boundaries as a short wavenumber disturbance with a 2Δt time frequency, and progressively spreads out throughout the whole computational domain. This unstable behavior indicates that some kind of implicit treatment is needed.
b. Toward an iterative implicit well-posed LBC treatment
An appropriate treatment in order to stabilize the explicit LBC specification consists of introducing some additional implicit LBC corrections into the interpolated explicit rhs terms. These implicit corrections have to be constructed in such a way that they compensate the detrimental effect of the explicit LBC imposition in the derivatives and the Coriolis terms. By taking into account the form of the finite-difference approximations made at boundaries, these implicit LBC corrections can be symbolically written as
where superscript g denotes the guess solution obtained though the governing model equations, and lbc denotes the well-posed boundary conditions. As explained before, these terms cannot be incorporated into the spectral implicit part of the SI time stepping. In what follows, we propose an alternative solution for dealing with these additional implicit correction terms.
Here also, it is assumed that u > 0, without loss of generality. As in section 3a, the case u < is considered. The additional boundary corrections are designed for the characteristics LBC, hence Φlbc,+0 = (ph,+0 − u+0)/c, υlbc,+0 = υh,+0 at i = 0 (inflow boundary) and Φlbc,+ iL = −(qh,+iL − u+iL)/c at i = iL (outflow boundary). The guess boundary values Φg,+0, υg,+0, and Φg,+iL must be dictated by the 2-TL SISL discretized governing equations (4)–(6) at the boundaries. For consistency, let us use the simplest possible approximations: Φg,+0 = RΦ0 − (Δt/2Δx)(u1 − u0)+, υg,+0 = Rυ − (Δt/2)f u+0, and Φg,+iL = RΦiL − (Δt/2Δx) (uiL−uiL−1)+. These latter definitions of (Φlbc,+0, υlbc,+0, Φlbc,+iL), and (Φg,+0, υg,+0, Φg,+iL) are introduced into (19)–(22). Then, adding the implicit boundary correction to (5) for u field at the boundaries, the following system is obtained:
This equation has the same form as (5) with an additional term [δRu]+ that can be understood as a boundary correction term, given at the boundary grid point i = 0 and i = iL by
and at the adjacent boundary grid point i = 1 and i = iL − 1, the boundary correction terms are given by the relationships [δRu1]+ = 1/2[δRu0]+ and [δRuiL−1]+ = 1/2 [δRuiL]+. At the interior of the physical domain (i = 2, . . . , iL − 2), the boundary correction terms [δRui]+ are equal to zero. Notice that no boundary correction terms are required in (4) and (6); hence, [δRΦi]+ = [δRυi]+ = 0.
Equation (27) together with (4) and (6) form a closed system. It can be noticed that the boundary correction [δRu]+ has an implicit component due to the presence of the derivative term (∂Φ/∂x)+in its mathematical expression. This implies that the resulting operator has to be inverted over the whole domain, which in the spectral model is performed in spectral space. Unfortunately, because of the local form of this boundary correction term [δRu]+, the operator to be inverted has spatially variable coefficients, which makes the direct inversion in spectral space impossible in the simple usual form. However, this problem can still be solved through the spectral method, by making use of an iterative algorithm with a constant-coefficient preconditioner, in the same spirit as the ICI scheme proposed by Cullen (2001) and Bénard (2003). This iterative treatment will converge under suitable conditions. For example, the 2TL SISL scheme can be written as
In what follows, the quantities [δRϕ](k), [δRu](k), [δRυ](k) are called “the iterative LBC correction terms.” The superscript k is the integer index for iteration. For the sake of simplicity, the Coriolis terms are also treated iteratively. Thus, υh0 is supplied at t +Δt, at inflow directly by imposing [δRυ](k)0 = υh,+0 − Rυ0 + (fΔt/2)u(k)0. Iterative LBC correction terms for Φ in (4) are still not required, [δRΦ](k) = 0. The guess values Ũ+i are evaluated following an iterative relationships:
Finally, this yields the following expressions for LBC correction terms:
at the adjacent boundary grid point [δRu](k)1 = 1/2[δRu](k)0 and [δRu](k)iL−1 = 1/2[δRu](k)iL. This solves the implicitness of the derivative (∂Φ/∂x)+ and this also enables to circumvent the issue of the local structure of the boundary correction terms [δRu]+.
Please note that the convergence might depend on the specific form of the terms in (19)–(22). It is also important to emphasize that other options are possible. For instance, higher-order off-centered finite-difference approximations could be used to estimate (i) the derivatives (∂Φ/∂x) and (∂u/∂x) at boundaries in the explicit rhs terms (7)–(8), and (ii) the guess values Φg,+0 and Φg,+iL of the implicit correction terms in (19)–(22). However, this option would imply an increase in the number of unknown gridpoint values u+i in the system in (23)–(26), which leads to a much more complicated derivation of the boundary correction terms [δRu]+. An extensive study of such options lies beyond the scope of this paper.
In the case u > (not detailed herein), in which three fields have to be imposed at the inflow boundary (i = 0) and none at the outflow boundary (i = iL), we can proceed in a similar way as previously to design the iterative LBC correction terms [δRΦ](k), [δRu](k), and [δRυ](k) at the inflow boundary grid points i = 0 and i = 1 for (4)–(6). Although the spirit of the method is the same, the resulting details of the formulation significantly departs from previous ones in (35)–(38), and are not detailed here.
Let us now give a short description of the general time step organization of the iterative LBC correction procedure, regardless the values of u and . First let us assume that the explicit interpolated rhs terms are computed in gridpoint space as described in section 3a, henceforth designated by Rlbc. Then, the starting guess at the zeroth iteration is chosen as, X(0) = X0, and ∂xX(0) = ∂xX0. Thus, for k = 0, . . . , Niter − 1, the iterative LBC correction terms [δR](k) are computed using the known fields Rlbc, X(k), ∂xX(k), and the boundary data fields (denoted by Xh,+). Afterward, the final rhs terms R̃ = Rlbc + [δR](k) are periodically extended in the E zone and then transferred into spectral space by applying the FFT direct transform. The implicit problem [𝗜 − (Δt/2)𝗟]X̂(k+1) = is solved in spectral space as explained in section 2. The iterative process resumes by computing ∂xX(k+1) spectrally. To do this, an additional inverse FFT transform is needed at each iterations. The details of the data flow of the iterative LBC correction scheme is shown in Fig. 2. At the end of the iteration procedure the forecast is given by X+ = X(Niter).
Provided that only few iterations are required, this iterative LBC procedure can be easily incorporated into a general ICI scheme proposed by Cullen (2001) and more recently detailed by Bénard (2003). From a more practical standpoint, a specific ICI SL scheme, with Niter = 2 (predictor–corrector), has been implemented in the French NWP operational limited-area spectral models: the nonhydrostatic version of ALADIN (ALADIN-NH) and Application de la Recherche à l’Opérationnel pour la Méso-Echelle (AROME). In such a context the idea of an iterative well-posed LBC scheme is very attractive.
c. Numerical tests
The stability of this iterative LBC correction scheme is now discussed, using the same experimental environment as in section 3a. The main purpose is to address whether this new well-posed boundary scheme is stable for large wave Courant number as currently used in NWP application. The following experimental settings are used: Δt = 400 s, u = 12.5 m s−1, c = 300 m s−1, and Δx = 10 km; hence, α = 0.5 and the wave Courant number is equal to 12, which corresponds to a common target value in most of NWP pattern. For these experimental settings, numerical tests for different number of iteration Niter have shown that the solution remains stable provided that Niter ≥ 4. The mathematical reason of this stability criterion has not been investigated in this paper, we restrict our investigations to experimental observations only. This induces at least a factor of 4 in the computational cost of the implicit part of the dynamics, an overcost that is however similar to the one found in practice in the finite-difference context (see McDonald, 2000, p. 4049). Figure 3 represents the fields u ± ckΦ and υ after 50 min of integration for (Niter = 4), in the same format as Fig. 1. The forecast is stable with no sign of small-scale noise, and interestingly the numerical solution is consistent with the boundary conditions. Notice that a phase delay occurs between the analytical solution and the numerical solution. This can be accounted for by poor temporal representation of gravity waves by the SI scheme for large time steps.
These tests provide an experimental proof of the robustness of the iterative LBC correction scheme proposed herein. In particular, it has been shown that at least four iterations are required to maintain stability when an operational value of a wave Courant number is used [i.e., c(Δt/Δx) ≃ 12]. In what follows, we discuss the SL trajectory computations issue at inflow boundaries for large value of |α|.
4. SL trajectory at boundaries
When the departure point is outside the domain, a specific treatment is required. The simplest solution is the above-referenced “trajectory truncation” method, but as pointed out by McDonald (2000), this strategy is only accurate for |α| < 1. To address the limitation McDonald suggested two options that will be discussed here. Additionally, an alternative way for dealing with the SL trajectory issue at inflow boundaries is proposed. For the sake of simplicity, it is assumed that 0 < u < c without loss of generality. Hence in what follows, we consider that α > 1. The case α < −1 can be treated in a similar way, but is not examined here.
a. Time-interpolation scheme
This method consists in replacing the trajectory truncation by a time-interpolation scheme at the inflow boundary (x = 0), exactly as in section 2c of McDonald (2000). It has been implemented successfully in the present spectral model version, but as pointed out by McDonald (2000), the slowness of convergence for this method can be problematic. The number of iterations required to get an accurate forecast was found to be dramatically high, which renders this technique quite inefficient and therefore unattractive for NWP purposes (it will hence not be detailed furthermore here).
b. Well-posed buffer zone
A buffer zone, external to the inflow boundary, is constructed, so that if the departure point is located inside this buffer zone, fields can be interpolated in the same way as if they were in the interior. Since u > 0, a buffer zone is needed for the region x < 0. The ingoing characteristics p and υ are imposed at (i = 0), the third characteristic variable q is also needed to compute fields in the buffer zone, which, to maintain well posedness, has to be not externally imposed but extrapolated in a stable manner. As suggested by McDonald (2000), the equations for the characteristic fields can be used in conjunction with a semi-Lagrangian scheme to extrapolate the field corresponding to the outgoing characteristic q at the inflow boundary. The equations for the characteristic fields can be written as:
Equation (40) is solved by the semi-Lagrangian technique, by finding the value of the field at the departure point x*q = x0 − (u − c)Δt, it yields
Thus, q+0 at x = 0 are known, meaning that (∂q/∂t) and (∂2q/∂t2) can be approximated there. The characteristics values p0 and υ0 are imposed, thus (∂p/∂t), (∂υ/∂t), and (∂2p/∂t2), (∂2υ/∂t2), can also be estimated at x = 0. Thus, (39)–(40) can be used to compute (∂p/∂x), (∂q/∂x), which require a potentially problematic division by (u ± c), and (∂υ/∂x) can be obtained from (41). From these, (∂Φ/∂x) and (∂u/∂x) can be determined using the definitions of q and of p. Differentiating (39)–(41) with respect to x, the second-order terms (∂2Φ/∂x2), (∂2u/∂x2), and (∂2υ/∂x2) can be computed. Then the Taylor expansion is used to compute Φ, u, and υ in the buffer zone with a O(Δt2) level of accuracy. However, this method is no longer applicable when |u| ≃ |c|, which may occur in case of inertia–gravity wave propagation.
c. Substepping scheme at inflow
An alternative method based on the concept of extrinsic LBCs introduced by Termonia and Voitus (2008), is proposed here to address the SL trajectory issue. The key point of this method consists in replacing the trajectory truncation value of Ri, for i = 0, . . . , [α] at the inflow by a better approximation of the lhs terms [I − (Δt/2)L]X+, when α > 1. Obviously, to do this, a guess of the state vector X at time level t + Δt (denoted X̃+ hereafter) has to be estimated independently from the dynamical core of the SISL spectral model over a buffer zone1 of Nb = [α] + 1 internal grid points from the inflow boundary, (where [α] is the integer part of α). A substepping explicit time scheme between time level t and t + Δt with characteristics LBC is performed to yield X̃+ inside the buffer zone. For instance, in this paper, a characteristics-based backward/forward time discretization has been used as explicit substepping scheme to discretized the (39)–(41) such as
with τ = Δt/Nτ and Nτ = 1 + [(|u| + c) (Δt/Δx)]. The superscripts 0τ and +τ correspond to the time level t + nτ and t + (n + l)τ, respectively, where n is an integer between 0 and Nτ − 1. The subscripts *p, *q, and *υ denote interpolations to the departure points associated to each characteristic fields, which are given respectively by x*p = xi − (u + c)τ, x*q = xi − (u + c)τ, and x*υ = xi − uτ. To compute the substep at time t + (n + l)τ, the rhs terms Zpi, Zqi, and Zυi are computed at the points 0, . . . , (Nτ + Nb) − n. The values of p, q, and υ are then computed at the points [i, t + (n + l)τ] lying inside the buffer 0, . . . , (Nτ + Nb) − n − 1. Except at the boundary (i = 0), the choice of the time step τ ensures that trajectories for each characteristic are always contained inside the buffer 0, . . . , (Nτ + Nb) − n. This space–time organization of the substepping computations is illustrated in Fig. 4. At the inflow boundary (i = 0), P̃+τ0 = p+τ,h0, and υ̃+τ0 = υ+τ,h0, the value of the outgoing characteristic q0 is provided by (43). At the end of the substepping computation (n = Nτ − 1), the interpolated rhs terms R lying in the inflow buffer 0, . . . , Nb − 1 are replaced by a finite-difference approximation of [𝗜 − (Δt/2)𝗟]X̃+ using the values ũ+, Φ̃+, and υ̃+, that are deduced from the substepping computation of P̃+, q̃+, and υ̃+ at the points 0, . . . , Nb. In this implementation, first- and second-order finite-difference schemes have been used to approximate derivatives (∂Φ/∂x) and (∂u/∂x) depending on the position with respect to the edge.
5. Numerical testing
In this section, two numerical tests are performed. First, the efficiency of the two alternative methods to the trajectory truncation method (previously detailed in section 4) are discussed using a pure advection test with a large value of α. Second, the permeability of the boundaries to gravity wave is investigated using a radiation test.
a. Advection test: Bell-shaped perturbation entering and exiting the area
The results of tests with the slow solution of (1)–(3) for which u = 0, and fυ = (∂Φ/∂x) are discussed here. In that geostrophically balanced case, the analytical solution transports these fields with a velocity u without changing its shape. There, the analysis and the initialization produce a well-balanced initial state. On the boundary, meteorological fields should enter the domain without corruption from gravity waves; this can be modeled by imposing the slow solution at the boundaries. The analytical solution used here corresponds to a bell shape for Φ in geostrophic balance with υ and with zero divergence for u:
This mirrors the meteorological situation of a front entering the area through the western boundary and going out through the eastern boundary. The integration starts by placing the maximum of the bell-shaped Φ solution on the western boundary, xs = 0. Boundary conditions consistent with the bell shape moving in the positive direction with the velocity u are imposed. The experimental parameters are L = 1000 km, Δx = 10 km, f = 10−4 s−1, c = 300 m s−1, Δt = 416.667 s, u = 100 m s−1, thus α = 4.167. In these tests, the length of integration (T) will be always chosen so that in the analytical solution the bell curve moves exactly a distance L = iLΔx. The characteristics boundary conditions are imposed so that p(0, t) = u(0, t) + cΦ(0, t), υ(0, t), and p(L, t) = u(L, t) − cΦ(L, t) are given by (45). The initial state is also given by (45) with t = 0.
First, the forecast is run with trajectory truncation and four iterations are used for the iterative LBC correction scheme. The outcome is shown in Fig. 5. For clarity, the E zone has been dropped. It can be seen that the front has entered and gone out of the area quite accurately, there are neither sign of instabilities nor spurious reflections of (2Δx) shortwave noises form at boundaries. But, as pointed out by McDonald, the negative aspects are, of course, the slowing down of the Φ and υ solutions. The use of well-posed buffer zone or substepping scheme idea at inflow retain the positive aspects of the solution and almost completely eliminates the previous phase shift, as can be seen in Figs. 6 and 7. The well-posed buffer zone works in a satisfactory manner, but as outlined previously in section 4b this result holds only because in that case the gravity wave velocity significantly departs from the wind speed. Let us now consider the “substepping scheme at inflow.” The substepping parameters are given by Nτ = 17, Nbuf = 5, and, hence, τ = 24.51 s. The application of this scheme at inflow leads to a numerical solution, which is very similar to the analytical one; the phase shift has been removed, as seen in Fig. 7. Thereby, it appears as a good alternative for solving SL trajectory issue at inflow boundaries.
b. Radiation test: Adjustment case
Here, the behavior of the iterative implicit LBC scheme is examined when gravity waves leave the area of meteorological interest. It is demonstrated that imposing characteristics p and q on boundaries allows gravity waves to pass through with little reflection. The same adjustment experiment as the one presented in McDonald (2000) is reproduced here. Considering the steady-state solution to the Rossby adjustment problem, described in Gill (1982, section 7.2.2). The initial state is given by
with u(x, 0) = 0 and υ(x, 0) = 0 over the whole domain 0 ≤ x ≤ L. Experimental parameters are given by f = 10−4 s−1 and = 300 m s−1, then L = 30 000 km is 10 times the Rossby deformation radius (/f). Thus, the geostrophically balanced large-scale flow at the boundaries will remain steady with time. Constant mean wind velocity is set to u = 1 m s−1. The remaining parameters are chosen as Δx = 100 km, Δt = 600 s, and Φ0 = 10. This test is performed using the iterative implicit LBC treatment described in section 3b, moreover, since α is much less than unity, the trajectory truncation scheme can be used without significant accuracy problems. Concerning the boundary value, the ingoing characteristics are imposed so that: p(0, t) = Φ0, υ(0, t) = 0, and q(L, t) = −Φ0. The 10-day forecast is displayed in Fig. 8. It can be seen that the solution is not corrupted by gravity waves; there is no visible sign of reflection at boundaries, demonstrating the efficiency of the iterative implicit “characteristics” LBC scheme in terms of transparency.
6. Conclusions and discussion
The feasibility of a well-posed lateral boundary strategy for the Fourier spectral one-dimensional model using SISL time-marching schemes has been experimentally demonstrated. An iterative implicit LBC correction algorithm has been designed to control the detrimental effects in term of stability on the SI scheme, due to the specification of the well-posed characteristics LBC in the gridpoint explicit part of the model dynamics. This iterative approach is straightforward to implement when employing an iterative centered-implicit time (ICI) scheme instead of classical SI schemes. The salient result of the experimental tests is that stability for a significant value of the wave Courant number, typically c(Δt/Δx) ≃ 12, can be achieved provided that more than four iterations (Niter ≥ 4) are performed. The theoretical reason for this has not been elucidated here. Additional studies are needed to investigate some alternatives enabling the acceleration of the convergence of the above-proposed iterative approach.
As far as the SL computations are concerned, the so-called trajectory truncation scheme has been found to be stable in experimental tests, but quite inaccurate for large |α|. Both time-interpolation and well-posed buffer zone schemes overcome this accuracy limitation of trajectory truncation. However, these two schemes have both unavoidable weaknesses, highlighted by McDonald (2000), which still hold true in a spectral context. With the substepping scheme, proposed here, some of the main drawbacks of the preceding scheme can be circumvented. Substepping scheme does not require additional iteration, it can be applied in any meteorological pattern and provides very accurate result. Nevertheless, it has its own limitations: (i) it is not an easy-implementing alternative, especially when moving to two- (or three) dimensional systems; (ii) the substepping algorithm is also quite demanding in computing resources because of its need of rather large buffers near the inflow boundary.
This preliminary implementation identified the difficulties to build a stable and accurate well-posed boundary scheme for spectral SISL discretization for a simple one-dimensional linearized system. An important question has not been addressed by this paper: will the presence of nonlinear terms cause unexpected problems at boundaries? It is well known that the nonlinear residual terms are sources of instabilities for the SI time scheme, and thereby might jeopardize the stability of the iterative LBC scheme. Besides, Oliger and Sundström (1978), pointed out that “the important case where velocities change sign on the boundary and do not vanish in a neighbourhood of such a boundary point, is not covered by the well-posed initial-boundary value problem existing theory.”
In this preliminary work, no theory of stability of the LBC iterative scheme has been presented. The feasibility of a well-posed LBC strategy for the Fourier spectral SISL scheme has been experimentally demonstrated in the particular case of the one-dimensional linearized shallow-water system. As a result, there is no evidence that this iterative approach will work in more realistic cases (e.g., bidimensional, three-dimensional, or in nonlinear context). Therefore, “deeper works” are necessary.
As a next study, it is essential to test this iterative idea in nonlinear simplified systems. Furthermore, if we try to extend this study to two dimensions (x–y) the complications will increase. In particular, corners are known to induce instabilites if they are not carefully discretized, see Elvius and Sundström (1973) or more recently McDonald (2002, 2003). Thus, we plan to use McDonald’s approach and the iterative idea discussed in this paper to construct a two-dimensional spectral SISL shallow-water model with well-posed characteristics boundary conditions.
The approach of imposing the LBC conditions in an iterative procedure has been suggested independently by M. Hortal while this work was carried out. The authors thank Aidan McDonald and Mariano Hortal for fruitful discussions.
Corresponding author address: Fabrice Voitus, CNRM/GMAP, 42, Avenue G. Coriolis, F-31057 Toulouse CEDEX, France. Email: email@example.com