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,

18.20. FRONTS AND TOPOGRAPHY: SMOOTHING/FILTERS 419

|| ,,

zz ||

yy

||

zz

,, y

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

∞ n

QF (», θ) ≡ σnm Qm Yn (», θ)

m

(18.68)

n

n=0 m=’n

where the σnm are the ¬lter coef¬cients. QF can also be represented in physical space as an

integral:

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

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

420

(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

n¤N

1

[triangular truncation] (18.70)

σnm =

otherwise

0

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

grace.

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

∞

1

(18.71)

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

chosen?

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

2

2

J≡ d„¦ ( Q(», θ) ’ QN (», θ)) + ± 2k

(18.72)

d„¦ QN

„¦ „¦

where „¦ denotes the surface of the sphere and the approximation QN is

N n

am Yn (», θ)

m

(18.73)

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

18.20. FRONTS AND TOPOGRAPHY: SMOOTHING/FILTERS 421

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

(18.74)

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

bm

n

am = (18.75)

n 2k (n + 1)2k

1 + ±n

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

n

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

800

600

400

200

0

-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.)

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

422

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

(18.76)

σ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)

σ ,

Np

N

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θ )

1

σErf c’Log (θ; p) ≡ erfc 2p1/2 θ (18.79)

2

2 4θ

where

1

θ ≡ |θ| ’ (18.80)

2

18.20. FRONTS AND TOPOGRAPHY: SMOOTHING/FILTERS 423

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

d(x)

(18.81)

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 ∈ [’π, π]

x/π

Sw(x) ≡ (18.82)