Transcript Introduction to Numerical Weather Prediction

Effects of the Numerical
Approximations
26 September, 1 & 3 October 2012
Thematic Outline
• Linear stability criteria for finite difference methods to
solving the primitive equations
• Impacts of finite difference methods upon the phase and
group speeds of meteorological waves
• Aliasing of information across wavelengths
• The formulation and utility of numerical diffusion
• Considerations related to vertical coordinate selection
Linear Numerical Stability
• To this point, we have only touched on the idea of
numerical stability in terms of either…
– Stating (without proof) that a given term or method is
stable or potentially unstable.
– General evaluation of stability using the basic definitions
for the Courant number and CFL criterion.
• Now, we want to step through how the CFL criterion
can be derived and discuss the stability criteria for a
few different finite difference schemes.
Linear Numerical Stability
• What is meant by stability, anyway?
• If the atmospheric variables are wave-like, how do
the waves evolve through time within a model?
– Stable: waves evolve for meteorological reasons
– Unstable: waves grow exponentially for numerical / nonphysical reasons
• As we will see, different finite difference methods
have slightly different stability criteria.
Linear Numerical Stability
• Three modes of stability, as applied to finite
difference methods...
– Absolutely stable (ideal, but uncommon)
– Absolutely unstable (cannot use at all)
– Conditionally stable (most common; stable for specific
ranges of model params and/or meteorological conditions)
• Why is stability a concern? There are a lot of partial
derivatives in the primitive equations!
Linear Numerical Stability
Example: Horizontal Momentum Equations
Note the large numbers of partial derivatives present below!
u
u
u
u uv tan uw 1 p
 u
v w 


 2( w cos  v sin  )  Fx
t
x
y
z
a
a  x
time derivs. advection terms
pres. grad.
v
v
v
v u 2 tan uw 1 p
 u  v  w 


 2u sin   Fy
t
x
y
z
a
a  y
φ = latitude, a = radius of the Earth, Ω = rotational frequency of Earth, F = friction
Linear Numerical Stability
• All partial derivative terms impact numerical stability,
but the advection terms are the most troublesome.
• Recall: the stability of the advection terms partially
motivated semi-Lagrangian and spectral methods!
– In this discussion on stability, note that we are in physical
space utilizing an Eulerian reference frame.
• We first desire to consider linear stability in the
context of a horizontal advection term.
Linear Numerical Stability
• 1-Dimensional Advection eqn. for a Shallow Fluid:
h
t

j
h
 U
x

j
τ denotes the time step, j denotes location along x-axis
• Let us assume a wave-like (harmonic) structure for h,
i ( kxt )
ˆ
h  he
k = 2π/L, L = wavelength, ω = U*k = frequency, hˆ = amplitude
Linear Numerical Stability
• The frequency has both real and imaginary parts, i.e.,
  R  iI
• Substitute into our expression for h:
i ( kx(R iI ) t )
i ( kxRt iI t )
ˆ
ˆ
h  he
 he
• If we split the exponential function, we get:
i ( kxRt ) i ( iI t )
i ( kxRt ) I t
ˆ
ˆ
h  he
e
 he
e
Linear Numerical Stability
 t
• The e I term acts to add temporal dependence to
the amplitude hˆ . Since t is positive-definite,
– ωI > 0 gives exponential wave growth / amplification
– ωI < 0 gives exponential wave damping
– ωI = 0 imparts no change upon the amplitude
• The imaginary (non-physical) part of the frequency
determines the numerical stability of the solution.
– Stability criterion must relate e I t to the Courant number.
– The real portion of the frequency determines wave phase.
Linear Numerical Stability
• Example: forward differencing in time, backward
differencing in space
• Apply these methods to the 1-D advection equation:
hj 1  hj
t
 U
hj  hj 1
x
• Our stability criterion is defined in terms of
we rewrite the above in terms of it:
 1
hj
Ut 
 hj  
h j  hj 1
x

Ut
,
x
so
Linear Numerical Stability
• Let x = jΔx (j = location on axis, Δx = grid spacing) and
t = τΔt (τ = time step, Δt = time increment).
• Pose the wave solution in terms of the above:
i ( kjx  Rt )  It
ˆ
h j  he
e

• The wave solutions for τ+1 and j-1 follow from
substituting those in for τ and j above, respectively.
Linear Numerical Stability
• Plug the wave solution into the finite difference
approximation. Simplify by dividing through by
common exponential terms to obtain:
Ut
 I t i R t
e e
1  
(1  eikx )
x
• Use Euler’s relations to substitute for the imaginary
exponentials, then separate the equation into real
and imaginary components to obtain:
U t
(cos(kx )  1)
x
U t
sin( R t )  
sin(kx )
x
e I t cos( R t )  1 
real
e I t
imaginary
Linear Numerical Stability
• Square both sides of each equation, then add them.
• The terms with ωR go away via trig identity, where
sin2(ωRΔt) + cos2(ωRΔt) = 1. This allows us to obtain
an expression for e t , given by…
I
eI t  1  2
Ut
Ut
(coskx  1)(1 
)
x
x
• Note: since t = τΔt,
eI t  eIt  (eI t )
Linear Numerical Stability
• The value of ωIΔt determines how the amplitude of
the wave evolves through time, where,
– ωI > 0 gives exponential wave growth / amplification
– ωI < 0 gives exponential wave damping
– ωI = 0 imparts no change upon the amplitude
• We need to understand how eI t varies as a function
Ut
of the Courant number, x .
Linear Numerical Stability
ωI > 0
ωI = 0
ωI < 0
Linear Numerical Stability
Ut
 t
• For x = 1, e I = 1 and ωI = 0. The amplitude is
constant through time at hˆ and the solution is stable.
Ut
x
• For
> 1, eI t > 1 and ωI > 0. The amplitude thus
grows exponentially and the solution is unstable.
Ut
x
 I t
e
• For
< 1,
< 1 and ωI < 0. The amplitude thus is
damped with time and the solution is stable.
Linear Numerical Stability
• Numerical stability only requires that the solution
not grow exponentially with time.
• Thus, the stability criterion for these finite
Ut
1 .
differencing schemes is given by
x
– U is a function of the meteorology
– Δx is chosen depending on the scale of the features we
want to resolve with a minimum of truncation error
– Δt is left as the only “free” variable to be selected
Linear Numerical Stability
• The stability of this finite difference scheme also
varies for features of different wavelengths.
– Recall: k = 2π/L, where L is wavelength.
– For smaller L, damping for negative values of ωI is large.
– For larger L, damping for negative values of ωI is small.
• This isn’t necessarily a bad thing – if a small feature
isn’t resolved well, do we want it in the solution?
Linear Numerical Stability
• Similar exercises can be done for other finite
difference schemes…
– Forward in time, three-point centered in space:
 Ut 
2
1 
 (sin kx )
 x 
2
e
 I t

• Recall: stability only occurs when ωI ≤ 0 (i.e., wave has constant or
decaying amplitude with time).
• This is only true when the value under the radical is zero, or Δt = 0.
• Thus, this scheme guarantees exponential growth for all Δt and Δx
and is thus computationally unstable and entirely unusable.
Linear Numerical Stability
– Three-point centered in both time and space:
Ut
sin kx  1
x
• As kΔx approaches zero, so too does sin kΔx (small angle theorem).
• More generally, 0 ≤ sin kΔx ≤ 1.
• Thus, Ut  1
x
– Three-point centered in time, five-point centered in space:
Ut
 0.73
x
This is slightly different than the CFL criterion for other schemes.
Linear Numerical Stability
• Aside: even-ordered (2nd, 4th, etc.) centered in space
and time schemes do not damp for any stable value
of the Courant number.
• In other words, for all wavelengths (all L, all k), eI t  1
– We saw this for the three-point centered example on the
previous slide.
• Contrast to the wavelength dependence for forward
in time schemes described first (math & figure).
Linear Numerical Stability: Key Ideas
• Specific CFL criterion depends upon temporal and
spatial finite difference methods used.
• Not all combinations of methods are usable.
– Ex: forward in time, centered in space.
• Most methods are conditionally stable, but a few are
absolutely stable (e.g., implicit methods).
Linear Numerical Stability: Key Ideas
Ut
x
• As all stability criteria depend upon
careful time step selection is crucial!
in some way,
– U is a function of the meteorological conditions.
– Δx depends upon the features to be resolved with a
minimum truncation error for the differencing scheme.
• If an inappropriate Δt is selected, the exponential
amplitude growth term will eventually cause the
model to crash via numerical overflow.
Linear Numerical Stability: Key Ideas
• Practical recommendation: set Δt = 6Δx, where Δt is
in seconds and Δx is in kilometers.
– Ex: Δx = 10 km, so Δt = 60 sec.
– Proof: Ut  1 , so U (60s)  1 and thus U  166.67m s1
x
(10000m)
• Note that current-generation NWP models use more
sophisticated differencing schemes than considered
in our stability analysis.
– We have focused only on the simplest of schemes.
– Thus, either know the model or read the documentation!
Stability Visualization
• UΔt is the distance traveled by an advecting feature
in one time step.
• If UΔt > Δx (distance traveled > grid spacing), the
solution becomes unstable (i.e., Courant # > 1).
unstable
stable
• We want to ensure that the distance traveled can be
resolved within adjacent grid points (not beyond)!
Other Stability Considerations
• Δx can vary over a domain as a function of the
selected map projection or grid system. This means
that stability can vary as a function of location.
• The full primitive equations include partial
derivatives beyond just the advection terms.
– When all of these terms are accounted for, the relevant U
is truly the advective speed “U” plus the speed of the
fastest wave on the grid (cP).
– Because of this, for fixed Δx, Δt must be smaller than if the
cP term were not present.
Other Stability Considerations
• If acoustic waves are permitted, cP can be large!
– Example: cP = 300 m s-1, “U” = 100 m s-1, so U = 400 m s-1.
– For the standard CFL criterion, t  1 , or t  2.5x
x
400
– Helpful to know what the model does and doesn’t permit!
Other Stability Considerations
• We must also consider the vertical advection term in
our stability analysis.
• Analogous stability criterion:
W t
1
z
– Here, z is the vertical coordinate.
– Criterion is for a three-point centered difference scheme.
• Issues to keep in mind…
– W is the advective speed plus fastest-moving wave speed.
– Δz often varies with height and location on the grid.
Other Stability Considerations
• W can be large: 20-80 m s-1 in resolved convection or
~300 m s-1 if acoustic waves aren’t filtered.
• As Δz typically << Δx, this stability criterion can pose
a greater limit on Δt than the horizontal criterion.
–
–
–
–
Example: w = 20 m s-1, Δz = 500 m: Δt ≤ 25 s
Contrast: U = 100 m s-1, Δx = 4 km: Δt ≤ 40 s
With acoustic waves: Δt ≤ 1.67 s and 13.3 s, respectively
Thus, even more motivation to filter the fast waves or use
a scheme that handles them with fewer stability issues
(e.g., split-explicit and implicit differencing methods)!
Other Stability Considerations
• The most restrictive stability criterion always applies.
• A small buffer in the selected time step is often used
to account for sensitivities in the stability analysis.
Phase and Group Speed Errors
• First, a review of phase and group speeds…
– Phase speed Cp: speed of a wave’s phase lines along the
wave vector (propagation of the wave)
• Cp = ω/k, or frequency divided by wavenumber
– Group speed Cg: speed of wave energy propagation
• Cg = ∂ω/∂k, partial derivative of frequency w/r/t wavenumber
• If Cp = Cg, the frequency is not a function of
wavenumber and the wave is non-dispersive.
• Otherwise, the wave is dispersive.
Phase and Group Speed Errors
• Consider the 1-D advection equation from earlier:
h
h
 U
t
x
The advective speed of the wave in h is given by U.
• Forward in time, backward in space approximation:
 1
hj
Ut 
 hj  
h j  hj 1
x

The finite difference approximation causes the advective speed
to no longer be equal to U.
Phase and Group Speed Errors
• How to obtain the advective speed?
• Recall: assuming a wave-like solution for h, there are
real and imaginary portions to the approximation:
U t
(cos(kx )  1)
x
U t
sin( R t )  
sin(kx )
x
e I t cos( R t )  1 
real
e I t
imaginary
• We first divide the imaginary part by the real part.
Phase and Group Speed Errors
• This gives:
 I t
e
e I t
Ut
sin(kx)
sin( R t )

x

Ut
cos( R t )
1
(cos(kx)  1)
x
• Eliminate exponential terms and note that sin x/cos x = tan x:
Ut
sin(kx)
x
t an( R t )  
Ut
1
(cos(kx)  1)
x
• Because Cp = ω/k, ω = Cp*k:
Ut
sin(kx)

x
t an(C p kt )  
Ut
1
(cos(kx)  1)
x
Phase and Group Speed Errors
• Take the inverse tangent of both sides:
Ut


sin(kx)


1

x

C p kt  t an  
U

t
 1
(cos(kx)  1) 


x


• Solve for Cp:
Ut


sin(kx)


1
1
x

Cp 
t an  
kt
 1  Ut (cos(kx)  1) 


x


• Cp is a function of Courant number and wavenumber (and, by
extension, wavelength).
Phase and Group Speed Errors
• Dependence of Cp upon Courant # for varying wavelength:
Remember: valid for forward in time/backward in space scheme only!
Phase and Group Speed Errors
• For a Courant number of 0.5, Cp = U. As the Courant
number approaches 1, Cp approaches U.
• For 0 < Courant # < 0.5, Cp < U and all waves move
slower than the advective speed.
• For 0.5 < Courant # < 1, Cp > U and all waves move
faster than the advective speed.
• Problem is most noticeable at shorter wavelengths.
Phase and Group Speed Errors
• Because Cp depends upon k, this wave is dispersive.
• But, U does not depend upon k; it does not disperse
the wave in h.
• Thus, we have what is known as erroneous
numerical dispersion of the wave.
Phase and Group Speed Errors
• Repeat this exercise for the 3-pt centered in time and
space differencing scheme discussed earlier. Here,
Cp 
1
 Ut

sin 1  
sin kx 
kt
 x

• Cp is again a function of k, so this approximation
method is also numerically dispersive.
• There are two waves of interest above.
Phase and Group Speed Errors
• First wave: positive component
– Approximation to the physical wave.
– Moves in the correct direction but slower than U.
• Proof: U > 1, so ∆t/∆x < 1 for stability. Also, 0 < sin k∆x < 1, and so
∆t*sin(k∆x)/∆x < 1.
• Second wave: negative component
– Represents a computational / non-physical wave mode.
– Moves in the opposite direction with an amplitude much
smaller than that of the physical wave.
Phase and Group Speed Errors
• Dependence of Cp upon wavelength for varying Courant #:
Remember: valid for 3-pt centered in time/space scheme only!
Phase and Group Speed Errors
• Large wavelengths: Cp asymptotically approaches U.
• Short wavelengths: Cp < U or Cp << U.
• Larger values of the Courant number provide values
of Cp that are closer to U.
• Note that Cp ≠ Cg, as expected for dispersive solution.
Phase and Group Speed Errors
• How does all of this bear out in an actual model
solution? Consider Figure 3.26 (next slide).
• Courant # = 0.1, U = 10 m s-1, ∆x = 1 km, ∆t = 10 sec.
• The 1-D advection equation has been integrated
forward in time until the wave is back where it
started (note: domain is periodic in zonal direction).
Phase and Group Speed Errors
• Evolution of wave using 2nd order space / time differencing:
initial/final wave
Remember: valid for 3-pt centered in time/space scheme only!
Phase and Group Speed Errors
• The largest wave only slightly lags the actual wave.
• Its amplitude is less than that of the actual wave.
– The actual wave is comprised of a superposition of waves
at all wavenumbers.
– The finite difference approximation results in the
numerical dispersion of these waves, however.
– Recall: this isn’t damping because we said this scheme
does not damp for all stable values of the Courant number.
Phase and Group Speed Errors
• From Fig. 3.25, the largest (shortest) waves should
move fastest (slowest).
– The longest waves are close to where they should be.
– The shortest waves remain on the right side of the domain.
• Numerical diffusion helps to smooth things slightly.
– Higher order schemes: primarily impact short wavelengths.
– Lower order schemes: impact all wavelengths.
– Results in a broadening / “smoothing” of wave amplitude.
Phase and Group Speed Errors
• Dependence of numerical dispersion upon Courant #
Phase and Group Speed Errors
• From Fig. 3.25, Cp is reduced the most (least) for
small (large) stable Courant numbers.
• This is confirmed by Fig. 3.27.
– Waves at all wavelengths move faster for increased
Courant number, giving a more realistic solution.
• Aside: for fixed U and ∆x, larger Courant numbers
result from the use of a longer time step.
– Despite longer time step, a better solution is achieved due
to the dependence of phase speed upon Courant number.
Phase and Group Speed Errors
• The specific dependence of Cp upon Courant number
varies with the chosen finite difference scheme(s).
• Higher-order (more accurate) finite differencing
schemes produce less numerical dispersion.
• Nevertheless, all differencing schemes will produce
some numerical dispersion – thus, we want to select
an appropriate scheme to mitigate it!
Phase and Group Speed Errors
• 2nd Order, 3-pt Centered in Time & Space Solution:
Phase and Group Speed Errors
• 4th Order, 5-pt Centered in Space Solution:
Phase and Group Speed Errors
• 3rd Order Runge-Kutta in time, 6th Order Centered in
Space (used by the WRF-ARW model):
Phase and Group Speed Errors
• Now, consider a meteorological application.
• The initial fluid wave may be viewed as akin to a
sharp gradient (e.g., a field along a frontal boundary).
• Forecast models using finite differencing schemes in
which numerical dispersion is an issue will tend to
weaken sharp gradients because of this alone!
Phase and Group Speed Errors
• Last point: group velocity (wave energy propagation).
• Recall: advection produces a non-dispersive wave.
• But, in finite difference land, it is numerically
dispersive. Thus, Cp does not equal Cg.
• Let us examine how group velocity is impacted by
wavelength and Courant number.
Phase and Group Speed Errors
• For the 3-pt centered in time/space scheme,
Cp 
1
 Ut

sin 1  
sin kx 
kt
 x

• By definition,
Cg 
C p k 

 1
Ut

1 


sin

sin
k

x


k
k
k 

t

x



• Simplifying,
Cg 
U cos kx
 Ut

1 
sin kx 
 x

2
Phase and Group Speed Errors
• Dependence of Cg upon wavelength for varying Courant #:
Remember: valid for 3-pt centered in time/space scheme only!
Phase and Group Speed Errors
• For longer waves, Cg asymptotes toward U, similar to
Cp. These waves and their energy are less dispersive.
• For shorter waves, particularly below n = 6, Cg is at or
below zero, implying wave energy propagation
opposing that of the wave itself!
• The magnitude of the group velocity is closest to U
for increasing Courant number, similar to Cp’s value.
Phase & Group Speed Errors: Summary
• Finite differencing schemes cause non-physical
numerical dispersion of atmospheric fields.
• Numerical dispersion is largest (smallest) for lowerorder (higher-order) differencing schemes.
• Longer wavelengths generally are less dispersive.
• The Courant number modulates numerical dispersion
slightly in a way specific to the differencing scheme.
Aliasing
• Aliasing is a problem that arises when two waves
approximated on a grid interact with one another.
– Non-physical exchange of energy between wavenumbers
– Implications to numerical accuracy and stability
• Recall: the atmospheric (dependent) variables can be
represented as having wave-like structure.
• Multiple waves can interact when they each appear
in one term of the primitive equations.
Aliasing
• Example: 1-D non-linear advection term
u
u
x
– u can be written as the sum of an infinite number of waves:

u   am cos k m x
m 0
2m
where am = amplitude, km = L , m = # of waves in the
domain, and L is the length of the grid.
u
–
is obtained analytically from u:
x

u
   am km sin km x
x
m 0
– Thus, this advection term reflects the multiplication of two
waves over all m.
Aliasing
u  
 

u
   am coskm x   an kn sin kn x 
x  m0
 n1

where we have introduced n, analogous to m, for clarity. The
summation over n starts at 1 because sin(0) = 0.
Apply a trig identity:
Result:
sin( x  y )  sin( x  y )
sin x cos y 
2
a n am k n
sin(kn x  km x)  sin(kn x  km x)
2
aa k  
2 
2 

 n m n sin  (n  m)
x   sin  (n  m)
x 
2  
L 
L 


Aliasing
• When two wavenumbers interact, they produce two
sine waves, one with wavenumber n+m and one with
wavenumber n-m.
• In continuous space, where all n and m are possible,
this isn’t a problem.
• But, on a discrete grid, not all n and m are possible
by the very nature of discretization.
Aliasing
• Consider a discrete grid, such as in Figure 3.31…
• It takes at least three grid points to resolve a wave.
– Prove by crudely representing a wave structure above.
– Thus, the smallest resolvable wavelength is 2∆x.
• How is aliasing manifest on this grid for our example?
Aliasing
• Let L, the length of the grid, be equal to jmax*∆x,
where jmax is the maximum wavenumber and is even.
• If you have one wave, i.e., the length of the grid is
equal to the wave’s length, then the wavelength is
simply jmax*∆x.
• If you have two waves, i.e., the wave’s length is half
of that of the length of the grid, then the wavelength
is simply ½jmax*∆x.
Aliasing
(# of waves)
• In general, the wavelength will be equal to the length
of the grid divided by the number of waves:
L jmax
wavelength 
x
j
j
• If 2∆x is the smallest resolvable wavelength, what is
the corresponding value of j? j = ½jmax
Aliasing
• Thus, half of all wavenumbers j (from 1 ≤ j ≤ ½jmax)
can be resolved.
• However, the other half (from ½jmax ≤ j ≤ jmax) cannot
be resolved.
• All of this is because of the discretization of the grid
(i.e., having a finite number of points on which to
resolve features).
Aliasing
• Assume that both m and n describe resolvable
features.
• There are two waves defined by their interaction: an
n+m and an n-m wave.
– All n-m are resolvable because all n and m are resolvable.
– Some n+m are resolvable, but others are not (i.e., n+m >
½jmax).
Aliasing
• Let’s look at the n+m wave a bit more closely…
2

sin  (n  m)
L


x

• We defined L = jmax*∆x.
• x is simply some location along the grid, which can
be defined by the grid increment ∆x times the grid
point index j (j∆x).

2jx 

sin (n  m)
jmaxx 

Aliasing
• We want to consider only the problematic waves, but
more generally, we want to relate n+m to jmax.
• We accomplish this by stating n+m = jmax – s, where s
is less than ½jmax…


2jx 
2jx 
  sin  ( jmax  s )

sin  (n  m)
jmaxx 
jmaxx 


 2j  jmax  s  

2js 




 sin 
 sin  2j 

jmax
jmax 



Aliasing
• Apply another trig identity:
sin(x  y)  sin x cos y  cos x sin y

  2js 
  2js 
2js 
  sin 2j  cos
  cos2j sin

sin 2j 
jmax 

 jmax 
 jmax 
• By identity, sin(2πj) = 0 and cos(2πj) = 1. Thus…
  2js 
  2js 
 2js 





sin 2j  cos
 cos2j sin
 sin


 jmax 
 jmax 
 jmax 
Aliasing
• Note that L = jmax*∆x while x = j∆x.
• Substitute for ∆x in the expression for L to obtain:
j
jmax

x
L
• Thus, equivalently,
 2js 
2s 
  sin
sin
x
 L 
 jmax 
Aliasing
• This defines the unresolvable wave that arises from
the interaction of two resolvable waves.
• This wave has a wavenumber s, where s 
jmax  (n  m) .
• Now, let’s consider an example of two resolved
waves of wavelengths 2∆x and 4∆x.
Aliasing
• What are the corresponding wavenumbers?
L jmax

x
j
j
jmax
jmax
2 x 
x  n 
n
2
j
j
4x  max x  m  max
m
4
wavelength
• Consider the interacting n+m wave…
nm 
jmax jmax 3 jmax


2
4
4
Aliasing
• This wave is not resolvable!
• The wavelength can be obtained from the inverse of
the relationship on the prior slide, i.e.,
L jmax

x
j
j
j
4
wavelength max x  x
3 jmax
3
4
wavelength
• But, this unresolvable wave has energy! Where does
that energy go?
Aliasing
• Recall:
s  jmax  (n  m)
3 jmax jmax
s  jmax 

4
4
• The energy from the unresolvable wave thus is
contained in a wave with wavenumber j 4 .
max
• We know this wavenumber is associated with the
resolvable 4∆x wavelength.
– Thus, energy from unresolved waves is leaking into the
resolved waves! This is bad!
Aliasing
• Aliasing can be viewed in the context of Fig. 3.32…
• Energy is folded over the smallest resolvable
wavelength (2∆x) into the resolved wavelengths.
Aliasing
• How does this impact modelers or model users?
Energy
– Energy at wrong wavelengths can lead to errors in the
solution and is a possible source of instability.
– Example of energy issues: Figure 3.33…
Wavenumber
Aliasing
Energy
• Panel (a): model including numerical diffusion
– Model spectrum matches correct spectrum up to the
effective resolution (typically ~6∆x).
– Between the effective resolution and 2∆x, energy is
dampened due to diffusion.
Since we cannot resolve
these scales well, this is
desirable.
– Beyond 2∆x, there is no energy.
Wavenumber
Aliasing
Energy
• Panel (b): model without numerical diffusion
– Model spectrum resembles correct spectrum up to the
effective resolution (typically ~6∆x).
– Erroneous energy is found between the effective
resolution and 2∆x due to
aliasing.
– This is BAD. We cannot resolve
these scales well, plus these are
the scales with notable
phase/group speed issues!
Wavenumber
– Shows one benefit of diffusion.
Aliasing
• In this example, aliasing most strongly impacts the
2∆x to 4∆x wavelengths.
• How can we prove this?
• Consider Figure 3.32 once again. There are 42
combinations of n+m that result in aliasing…
n  12, m  1  12
n  11, m  2  11
n  10, m  3  10
n  9, m  4  9
n  8, m  5  8
n  7, m  6  7
Aliasing
• Next, consider the values of n+m…
– For n+m between 13 and 18, Figure 3.32 shows that the
energy goes into the 2∆x to 4∆x wavelengths.
– For n+m greater than 18, the energy goes into waves of
length greater than 4∆x.
– Of the 42 possible n+m values, 30 are between 13 and 18.
• The accumulation of energy in these wavelengths
can lead to non-linear instability and the failure of
the model simulation.
Aliasing
• How can we mitigate the impacts of aliasing?
– Use spectral or semi-Lagrangian methods such that spatial
variability is handled analytically.
• Spectral methods also explicitly handle unallowable wavenumber
interactions.
– Use a differencing scheme that selectively dampens the
shorter wavelengths.
• Forward in time, backward in space enabled damping.
• Even-ordered centered schemes do not dampen.
– Employ numerical diffusion to dampen the energy
contained at the shorter wavelengths (e.g., Figure 3.33).
Diffusion
• Diffusion acts to spatially spread or smooth features
in the heat, moisture, and momentum fields.
– This can act to dampen the amplitude of such features.
– Typically applied only to selected scales (selective filtering)
• Diffusion takes on two primary forms…
– Physical: manifest as mixing and turbulence in the
planetary boundary layer and free atmosphere
– Numerical: artificial; alleviates issues related to lateral
boundary noise, numerical wave dispersion, and aliasing
Diffusion
• Diffusion transfers higher values of a given field
down the gradient to where lower values are found.
• This acts to weaken the gradient of that field.
• Such gradients are physical or artificial in nature...
– The former are smoothed by turbulent mixing processes;
e.g., homogenization of a fluid.
– The latter are dampened to minimize spurious waves with
short wavelengths (numerical wave dispersion) and to
mitigate the potential for non-linear instability (aliasing).
Diffusion: Before
gradient small
gradient large
gradient small
resolvable wavelength unresolvable wavelength resolvable wavelength
10
12
14 16 22 28 30 32
34
36
1-D Frontal Example
∆T = 2°C
∆x = 10 km
Diffusion: After
gradient weaker
perhaps not very well-resolved, but better than before
22
10 12 14 16
28 30 32
34
1-D Frontal Example
∆T = 2°C
∆x = 10 km
36
Physical Diffusion
• Physical diffusion is typically parameterized.
– Boundary layer: contained within a planetary boundary
layer parameterization.
– Free atmosphere: contained within a turbulence
parameterization.
– Turbulence parameterizations are often encapsulated
within planetary boundary layer parameterizations.
• We’ll revisit both types of parameterizations in Ch. 4.
Numerical Diffusion
• Motivating Example #1: Numerical Wave Dispersion
• Finite difference approximations of varying accuracy
lead to varying degrees of wave dispersion.
• Diffusion operators of varying orders can be
developed to dampen selected wavelengths.
– Generally, higher order diffusion schemes focus their
damping properties on progressively smaller wavelengths.
Numerical Diffusion
• Evolution of wave using 2nd order space / time differencing:
initial/final wave
A progressive increase in damping of all wavelengths for lower order diffusion operators.
Numerical Diffusion
• Motivating Example #2: Aliasing
• Non-linear interaction of waves causes the erroneous
accumulation of energy at short wavelengths.
• If these wavelengths are selectively dampened,
energy at these wavelengths is reduced.
– Reduces effective resolution of model to that of the
smallest undamped wavelength.
– This is close to the wavelength that a grid can resolve well.
Numerical Diffusion
With Numerical Diffusion
Without Numerical Diffusion
Energy on y-axis
Wavenumber on x-axis
Numerical Diffusion
• It is a challenge to selectively or sufficiently dampen
the erroneous and poorly-resolved features while
not damping the physically-realistic features.
• Note that physically-realistic does not necessarily
imply accurate.
• How to represent numerical diffusion in the model?
Numerical Diffusion
• Typical numerical diffusion constructs take the form:
n
h
1
  12 K n  n h
t
h = any dependent variable, n = order of the diffusion operator (n = 0, 2, 4, 6),
and Kn = diffusion coefficient
• 0th and 2nd order diffusion constructs:
0
h
1
  12 K 0 0 h   K 0 h
t
2
h
1
  12 K 2 2 h  K 2 2 h
t
Numerical Diffusion
• Properties of the 0th order diffusion operator…
– Applies over all h; it is not scale-selective.
– Usage: mitigate noise along lateral and upper boundaries.
• Properties of the 2nd order diffusion operator…
– Diffusion depends upon the curvature of a field.
– Scale-selective, but less so than higher order schemes.
– Does not add new extrema; only acts to smooth curvature.
• Can repeat the above to get higher order schemes.
Numerical Diffusion: Example
• Recall: linear stability of finite differencing schemes
– The imaginary part of the frequency (ωI) determined
whether a solution would grow or dampen with time.
– While varying by finite differencing scheme, the value of ωI
typically depends upon wavelength in some fashion.
Example: forward-in-time,
backward-in-space
(graph is scheme-dependent)
Numerical Diffusion: Example
• Mathematically, as a harmonic solution for h,
I t i ( kxRt )
i ( kxt )
ˆ
ˆ
h  he
 he e
ωI = 0 describes a constant amplitude with time
ωI < 0 describes exponential amplitude damping
ωI > 0 describes exponential amplitude growth
I t
e
• We obtained expressions for
for several
combinations of finite differencing schemes.
– Earlier focus: when would solution exponentially grow?
– Now, we want to focus on when the solution will dampen.
Numerical Diffusion: An Aside
• Recall: selected differencing schemes intrinsically
damp short wavelengths for select Courant numbers.
– Example: forward-in-time, backward-in-space 2 slides back
– These schemes describe implicit diffusion.
• Diffusion applied in the form of an explicit diffusion
operator is known as explicit diffusion.
Numerical Diffusion: Example
• How does the stability analysis relate to diffusion?
– Diffusion operator is comprised of partial derivatives.
– It stands to follow, then, that a stability analysis may tell us
something about the wavelengths that it dampens.
• Consider the example of a forward-in-time,
2nd/4th/6th order centered-in-space 1-D diffusion,
n
h
nh
1
  12 K n n
t
x
h
 2h
 K2 2
t
x h
 4h
 K4 4
t
x h
 6h
 K6 6
t
x
Numerical Diffusion: Example
• Equation (3.63): e I t for each diffusion operator…
e I t


2  2 coskx  / x 2


4
6  8 coskx   2 cos2kx  / x 
 1  K n t 


6










20

30
cos
k

x

12
cos
2
k

x

2
cos
3
k

x
/

x


• Obtained by assuming harmonic structure for h and
complex structure for ω, plugging in to each
I t
e
equation, and manipulating the result to relate to
Numerical Diffusion: Example
• The value of ωI (and thus e I t ) depends upon…
– Damping coefficient K
– Time step ∆t
– Value of k∆x
• Since k = 2π/L, where L is wavelength – wavelength dependence.
• Damping occurs for e I t less than 1.
• For all “resolvable” wavelengths 2∆x to ∞, with K and
I t
e
∆t both positive,
will be less than or equal to 1.
Numerical Diffusion: Example
• Rewrite equation (3.63) for e I t in terms of L…
e I t



 2x  
2



2

2
cos

/

x





L








 2x 
 4x  
4

 1  K n t 
 6  8 cos
  2 cos
  / x 
 L 
 L 







2


x
4


x
6


x






6
 20  30 cos
  12 cos
  2 cos
  / x  



 L 
 L 
 L 
• The value of each cosine term depends upon the
horizontal grid spacing divided by the wavelength.
– Ratio is largest at small wavelengths and zero at infinity.
Numerical Diffusion: Example
• Small wavelength example: 2∆x
– L = 2∆x, so ∆x/L = ½.
– cos π = -1, cos 2π = 1, and cos 3π = -1.
e I t


 4 / x 2 
2  2 cos  / x 2




4
4
6  8 cos   2 cos2  / x 
 1  K n t 
  1  K n t 16 / x  


6
6












20

30
cos


12
cos
2


2
cos
3

/

x
64
/

x




• Long wavelength example: ∞
– L = ∞, so ∆x/L = 0 and cos 0 = 1.
I t
– In this case, the coefficients all add up to 0, and e = 1.
Numerical Diffusion: Example
Numerical Diffusion: Example
• In this example, the value of K is chosen so as to
completely dampen the 2∆x wave at every time step.
– Typically do not employ such strong damping.
– Does provide an elegant way of ‘normalizing’ the damping
for ease of interpretation, however.
• Slope of curve becomes steeper as order of diffusion
operator gets larger.
– This is desirable: shorter wavelengths get damped while
longer, better-resolved wavelengths are less damped.
Numerical Diffusion: Example
4∆x Wave
2nd order: eωit = 0.50
4th order: eωit = 0.75
6th order: eωit = 0.90
8∆x Wave
2nd order: eωit = 0.87
4th order: eωit = 0.97
6th order: eωit = 0.99
Numerical Diffusion: Example
2nd Order Diffusion
6th Order Diffusion
25∆x Square Wave
(Physically unrealistic, but sharp gradients help depict effects of diffusion.)
Run out 100 time steps; only term acting on wave is the diffusion operator.
Numerical Diffusion: Example
• 2nd order diffusion operator results in an amplitude
reduction of ~15% for this well-resolved wave.
• 6th order diffusion operator leaves the amplitude
largely unchanged while smoothing the gradient
(associated with a very short wave) at the corners.
– Note: 6th order operator can lead to local extrema (unlike
2nd order operator); must be accounted for in some way.
– Depiction on previous slide incorporates one such method.
Numerical Diffusion: Example
initial/final wave
Amplitude damping of waves at all wavelengths depends upon order of diffusion operator.
Numerical Diffusion: Considerations
• Diffusion must be calculated on a quasi-horizontal
surface and not a terrain-following surface.
• Example in text: isolated mountain feature…
– Typically, T(mtn) < T(sea level).
– Higher temperatures are diffused toward the top of the
mountain if diffusion applied on terrain-following surface.
– This artificially weakens the gradient and establishes a
non-physical vertical/overturning circulation!
Numerical Diffusion: Considerations
• Diffusion gives rise to the idea of effective resolution.
• Wavelengths of at least 2∆x can be ‘resolved’ on a
grid; longer waves are better resolved.
– If you dampen features on the smaller wavelengths, you
reduce your ability to ‘resolve’ those waves.
– Effective resolution is thus the smallest wavelength at
which approximately no damping is applied ( e  I t ~ 1).
Numerical Diffusion: Considerations
e I t  1
e
I t
1
7∆x is the shortest wavelength at which the simulated KE spectrum has an expected
structure; shorter wavelengths have been damped to reduce (spurious) energy.
Numerical Diffusion: Considerations
• Grid diffusion may also impact a model solution.
• We select an appropriate ∆t for our chosen ∆x in
order to maintain numerical stability.
– Information travels from one grid point to another at a rate
of speed given by ∆x/∆t.
– At every time step, a centered-in-space finite difference
depends upon the value of a variable ∆x away from it.
Numerical Diffusion: Considerations
• This information velocity is often faster than the
advective velocity!
– Example in text: ∆x = 25 km, ∆t = 350 s, ∆x/∆t = 71.43 m s-1
– …but, the advective velocity is only 5 m s-1!
– A parcel is stuck between grid points for 14.29 time steps.
• This can act to smooth fields somewhat, but it does
so in a non-physical, non-scale-selective manner.
Vertical Coordinates
• A wide range of vertical coordinate systems are used
in modern NWP applications.
• There exist three primary factors to consider when
selecting a vertical coordinate…
– Does it permit unevenly distributed vertical levels?
– How does it handle uneven terrain?
– Do any issues arise with transforming the derivatives in the
primitive equations into the chosen coordinate?
Vertical Coordinates: Considerations
• We earlier noted that vertical levels in a numerical
model are often not evenly distributed.
• Instead, more levels are typically desired where the
sharpest vertical gradients exist.
– e.g., in the boundary layer and near the tropopause
• Ideal: can unevenly distribute vertical levels based
upon the meteorological conditions without adding
unnecessary computational expense
Vertical Coordinates: Considerations
• The topography of the Earth is varied and is
associated with sharp horizontal gradients.
• Problem: in regions of sloping terrain, unless the
vertical coordinate follows the terrain, adjacent grid
points may be above and/or below ground.
• Ideal: must be able to handle sloping terrain in an
elegant yet numerically precise manner.
Vertical Coordinates: Considerations
• No matter what vertical coordinate is chosen, the
derivatives in the primitive equations must be
transformed into the new coordinate system.
• The transformation is fairly straightforward.
• Depending upon the chosen coordinate, however, it
can pose issues related to the computation of the
pressure gradient.
Vertical Coordinates: Possibilities
• Option #1: height above sea level (z)
• Two of three criteria met:
– Can be used with unevenly spaced vertical grids.
– Pressure gradient contains one term rather than two.
• Drawback: sloping terrain…
– Vertical levels can intersect terrain – not good.
– How to compute horizontal derivatives when one grid
point above ground is next to one below ground?
– How to efficiently handle ‘gaps’ in the grid?
Vertical Coordinates: Possibilities
z = 3000 m
z = 2000 m
z = 1000 m
Vertical Coordinates: Possibilities
• Option #2: pressure (p)
• Same two criteria met for p as for z.
– Added benefit: most observations are taken on pressure
levels, minimizing interpolation for model initialization.
• Same drawback related to terrain as with z as well.
– Added drawback: since p is not fixed to the terrain, the
points that are above/below ground can change with time!
– How to handle this temporally-changing grid structure?
Vertical Coordinates: Possibilities
p = 700 hPa
p = 850 hPa
p = 925 hPa
Vertical Coordinates: Possibilities
• Option #3: potential temperature (θ)
• θ is approximately conserved following the motion.
– Cross-isentrope flow is inherently small.
– In this coordinate, this equates to small vertical motions.
– Small vertical motion = small vertical advection and small
artificial grid diffusion in the vertical.
• Benefit: θ surfaces are naturally packed where
temperature gradients are large (i.e., fronts, etc.).
– Provides better resolution of features of interest.
Vertical Coordinates: Possibilities
pressure
potential temperature
Isentropes: solid
Isotachs: dashed
Isotachs: dashed
Vertical Coordinates: Possibilities
• Isentropes are tightly
packed in the middle
troposphere.
• Leads to very large
horizontal and vertical
gradients of θ and |v|
over a short distance.
Isentropes: solid
Isotachs: dashed
Vertical Coordinates: Possibilities
• Horizontal θ gradient is gone.
• Vertical θ gradient reduced.
• Each θ level follows along,
rather than cuts across, the
frontal boundary.
Isotachs: dashed
– Reduces magnitude of the
gradients of other fields.
Vertical Coordinates: Possibilities
• For this front, the isentropic coordinate gives us
improved vertical resolution in areas of interest.
– Matches one of our three criteria.
– Note the shading in Fig. 3.37 – can fit more grid points in
the box in (b) than in the box in (a)!
• Reduced gradient magnitudes improve the accuracy
of finite difference approximations for spatial and
vertical derivatives.
– Theoretically, this should improve forecast skill!
Vertical Coordinates: Possibilities
• Terrain: suffers from same issues as p and z.
• Pressure gradient: two terms result from transform
into isentropic coordinates.
–
–
–
–
p  M , where M is the Montgomery streamfunction.
M has two terms related to z and T.
An extra derivative is an extra source of truncation error.
This can lead to excessive non-physical errors being
introduced into the calculation of the pressure gradient.
Vertical Coordinates: Possibilities
• An additional issue: change of θ with height...
Potential temperature does not
always increase with increasing
height!
• Thus, one θ can equate to two
different levels!
• Must be accounted for artificially.
Vertical Coordinates: Possibilities
• Options #4 & 5: terrain-following sigma (σ-z, σ-p)
• Vertical coordinate based upon a normalized
pressure or normalized height…
p  pt

ps  pt
z  zt

z s  zt
s = surface, t = top of model
• Allowable values for σ range from 0 (p = pt, z = zt) to
1 (p = ps, z = zs).
Vertical Coordinates: Possibilities
Vertical Coordinates: Possibilities
• Sigma coordinates based off of pressure can change
with time (like p); height-based sigma preferred.
• The σ coordinate elegantly handles terrain issues.
What about the other two criteria?
– Non-uniform grid: can do this (see previous figure).
– Pressure gradient: split into two, based off of ps and φ.
• ps (surface pressure) and φ (geopotential) derivatives can be large
in regions of sloping terrain.
• With two derivatives, the same truncation error issues seen with
the θ coordinate are a concern here.
Vertical Coordinates: Possibilities
• σ or σ-like coordinates (e.g., the Eta coordinate
described in the text) are used in many popular
numerical models (WRF, MM5, GFS, etc.).
• It is also possible to define a hybrid coordinate, such
as the σ-θ coordinate described in the text, that
combines the best of each coordinate…
– In this example: packing of θ coordinate with the terrainfollowing capability of the σ coordinate.
– The σ (θ) coordinate would be used at lower (upper) levels.