. 86
( 136 .)


on a wide variety of scales including scales far below that which can be resolved by a global
model. Worse still, the topography has a discontinuous ¬rst derivative at the coasts where
the non-zero slope of land abruptly gives way to the zero slope of the oceans. (Of course,
the ocean is continually rippled by waves, but on the scale of a global forecasting model,

|| ,,
zz ||
,, y

Real Gibbs

Figure 18.13: Left: The real world. The lower boundary for an atmospheric model is the
jagged surface of the land (black) plus the ¬‚at surface of the ocean (confetti-shading). The
boundary has a discontinuous slope at the coast. Right: “Gibbs” world. When the lower
boundary of the atmosphere is approximated by a ¬nite sum of spherical harmonics, large
valleys are created in the ocean. The frowning ¬sh had better learn to breathe air!

these wind-induced ripples are negligible.) It follows that if the topography is expanded
as a series of spherical harmonics, the expansion will have only an algebraic rate of con-
vergence with large Gibbs ripples. (These ripples are sometimes called “spectral ringing”.)
Even at high resolution, the partial sums of the topography may have undulating valleys
one hundred meters deep in the oceans, as illustrated in the cartoon Fig. 18.13.
Both fronts and topography require smoothing or ¬ltering of spectral series. For fronts,
the ¬ltering is applied at every time step as an arti¬cial dissipation. Since front problems
are intimately connected with spectral blocking, aliasing instability, and other dif¬culties
described at length in Chapter 11, we shall concentrate on the dif¬culty peculiar to mete-
orology, which is the highly irregular lower boundary of the ¬‚uid. The topography needs
to be ¬ltered only once in a pre-processing step. However, the mechanics of making good
¬lters is similar for both once-only and every-step ¬ltering.

18.20.2 Mechanics of Filtering
If Q(», θ) denotes a function on the surface of a sphere with spherical harmonic coef¬cients
Qm , then its ¬ltrate QF can be represented in spectral space as

∞ n
QF (», θ) ≡ σnm Qm Yn (», θ)
n=0 m=’n

where the σnm are the ¬lter coef¬cients. QF can also be represented in physical space as an
2π π
QF (», θ) ≡ (18.69)
d» sin(θ )dθ w(» , θ ) Q(» , θ )
0 0

where the primed coordinates denote latitude and longitude in a coordinate system that
has been rotated so that the point (», θ) is the north pole of the rotated coordinate system.
The reason for this apparently complicated tactic of a different coordinate system inside
the integral for each point is justi¬ed by a theorem proved by Sardeshmuhk and Hoskins

(1984): if the weight function w(» , θ ) is “isotropic”, that is, a function only of the distance
from the point where Q is evaluated inside the integral to the point where QF is being
computed, then w is a function only of θ in the rotated coordinate system. Furthermore, if
the ¬lter is isotropic, they prove that the weights σnm must be a function of n only, that is,
must be the same for all spherical harmonics of the same subscript.
Sardeshmukh and Hoskins note that truncation of a spectral series is just a special case
of ¬ltering in which all the ¬lter weights are either one or zero. The triangular truncation
of order N , usually denoted as “TN ”, is the isotropic ¬lter

[triangular truncation] (18.70)
σnm =
The rhomboidal truncation, which keeps the same number of basis functions for each zonal
wavenumber [superscript] m, is not isotropic ” another reason why it has fallen from
For an isotropic ¬lter, the weight function varies only with θ in the integral representa-
tion of QF . It can therefore be expanded as a one-dimensional series of ordinary Legendre
polynomials (Sardeshmukh and Hoskins, 1984):

w(θ ) = σn (2n + 1)Pn (cos(θ )
4π n=0

In the limit that all σn = 1, this is just the Legendre series for the Dirac delta function.
When the series is truncated, the weight function is a ¬nite Legendre approximation to the
delta-function (Fig. 18.14).
Unfortunately, this approximation for ¬nite N is rather ugly. It resembles a Bessel func-
tion with a maximum at the origin and oscillations that decrease very slowly with radius.
Thus, values of Q far from the point (», θ) have a strong in¬‚uence on the values of QF at
that point. A weight function that decays rapidly and monontonically is much more reason-
able than one with wiggles and slow decay, as produced by truncation.
In this sense, truncating a series is a rather stupid idea: a ¬ltered approximation, if the
¬lter weights are chosen properly, is much more reasonable. But how are the weights to be

18.20.3 Spherical splines
Truncation of a spectral series generates an approximation which is the best possible in a
least-squares sense. The reason that the truncated series is unsatisfactory is that the ap-
proximation is very wiggly. Ocean valleys are both ugly and unphysical!
The remedy is to compute an approximation which minimizes a “cost” function that
not only penalizes least-square error, but also excessive wiggles. Spherical splines of order
k are a weighted sum of spherical harmonics which minimizes the cost function
J≡ d„¦ ( Q(», θ) ’ QN (», θ)) + ± 2k
d„¦ QN
„¦ „¦

where „¦ denotes the surface of the sphere and the approximation QN is
N n
am Yn (», θ)
QN = n
n=0 m=’n

When the parameter ± = 0, the cost function is the usual mean-square error and QN is
the un¬ltered truncation of the series. As ± > 0 increases, the cost function depends more

and more on the smoothness of the approximation since wiggles imply large coef¬cients
for high degree spherical harmonics. Recalling that

Yn = ’n(n + 1)Yn
2 m m

it follows that high degree coef¬cients contribute disproportionately to the second term
in the cost function. The cost can only be minimized by reducing the wiggles, that is, by
reducing the large n coef¬cients, even though this slightly increases the mean-square error.
One can prove that the minimum of the cost function is rigorously given by

am = (18.75)
n 2k (n + 1)2k
1 + ±n

where the bm are the usual un¬ltered spherical harmonic coef¬cients.

Lindberg and Broccoli(1996) have applied spherical splines to topography. Wahba(1990)
reviews the general theory. The method can be generalized by replacing a single power of
the Laplacian in the second term by a weighted sum of Laplacians to make the weights
in Eq.(18.75) almost anything that one pleases. Almost any reasonable choice ” Lindberg
and Broccoli set k = 1 so that the smoothness penalty was the mean-square norm of the
Laplacian ” will produce a smoother approximation than truncation.

Unsmoothed weight






-0.2 0.2
0 0
0.2 -0.2

Figure 18.14: Mesh plot of the weight function for an un¬ltered truncation of a spherical
harmonics series. (N = 100.)

18.20.4 Filter Order
It is convenient to describe the properties of a ¬lter in terms of a continuous function σ
such that the discrete weights are given by σ for arguments ranging between zero and one:
σn = σ , n = 0, 1, 2, . . . , N
N +1
De¬nition 41 (Filter Order) A ¬lter is said to be of ORDER “p” if
j 1
∼O N ’ ∞, ¬xed j (18.77)
σ ,
or equivalently, if 1 ’ σ is an analytic function that has a (p ’ 1)-st order zero at the origin.
The concept of the order is important for two reasons. First, if the function is in¬nitely
differentiable so that its spectral series is converging exponentially fast, then it is both nec-
essary and suf¬cient for the ¬lter order to be p so that the errors produced by the ¬ltering
are O(1/N p ) where N is the truncation. Put another way, it is desirable to use a ¬lter of
high order to avoid throwing away the advantages of a high order method when comput-
ing smooth solutions. Second, when the solution is not smooth and its un¬ltered series
converges only at an algebraic rate, it has been shown that, at least in one dimension, one
can recover the true solution to within O(1/N p ) by using a ¬lter of order p at any ¬xed
distance away from the singularity of the solution (Vandeven, 1991).
The bad news is that the proportionality constant in Vandeven™s theorem is propor-
tional to the p-th power of distance from the singularity, so that ¬ltering, at least ¬ltering
through weighted partial sums, is unable to avoid O(1) errors in the immediate neighbor-
hood of the front or shock. This is not too surprising because a slight shift in the location
of a jump discontinuity implies that a small region that should be on the low side of the
jump is placed on the high side, creating a pointwise error equal to the magnitude of the
jump. Away from the discontinuity where the solution is smooth, however, we can ¬lter
the slowly convergent series to recover the function to arbitrarily high order ” if and only
if we use a high order ¬lter.
Vandeven shows that to be truly p-th order, a ¬lter σ(ζ) must have also have a (p ’ 1)-st
order zero at ζ = 1. Spherical smoothing splines, even for large order k, technically ¬‚unk
these conditions. So also does the exponential ¬lter favored by Hoskins(1980),
σ (exponential) (ζ) ≡ exp(’constant(n(n + 1))k ζ) (18.78)
If, however the ¬lter weights are exponentially small for the high n coef¬cients, the ζ = 1
conditions are approximately satis¬ed. Spherical splines (for high powers k of the Laplace
operator) and the exponential ¬lter thus work as high order ¬lters in practice.
Vandeven himself devised a ¬lter which satis¬ed all the conditions of his theorem
exactly. Unfortunately, his “ultimate” ¬lter is an incomplete beta function. Boyd(1996c)
showed that there is a close relationship between Vandeven™s ¬lter and the two hundred-
year old series acceleration method of Euler. He also demonstrated that both can be ap-
proximated quite accurately by an analytical function that satis¬es Vandeven™s conditions
exactly, the “Erfc-Log” ¬lter:
± 
 2
’ log(1 ’ 4θ )
σErf c’Log (θ; p) ≡ erfc 2p1/2 θ (18.79)
 
2 4θ
θ ≡ |θ| ’ (18.80)

18.20.5 Filtering with Spatially-Variable Order
Boyd(1995e) shows that the optimum order p for a function with a discontinuity is an in-
creasing function of distance from the singularity. This is consistent with the philosophy of
Flux-Corrected Transport algorithms for ¬‚ows with shocks: low order method around the
shock, a higher order method farther away. For a Fourier series with a spatial period of 2π,
the optimal order is
poptimum (x) = 1 + N 2π

where d(x) is the distance from point x to the singularity. (This distance is calculated by
interpreting the interval x ∈ [0, 2π] as a circle.)
As a one-dimensional example, consider the “sawtooth” function de¬ned by
x ∈ [’π, π]
Sw(x) ≡ (18.82)


. 86
( 136 .)