Vertical coordinates

Schär, ETH Zürich

Vertical Coordinate Formulations for Atmospheric Models

Christoph Schär ETH Zürich schaer@env.ethz.ch

Supplement to Lecture Notes

“Numerical Modeling of Weather and Climate”

April 3, 2007

Schär, ETH Zürich

2

Vertical coordinates

z = const p = const

σ = 0

σ = 1 p = ps p = 0 σ = const

Terrain-following coordinates (

σ

-coordinates)

Terrain-following coordinates in ocean models

Height coordinates Pressure coordinates

Schär, ETH Zürich

Outline

Introduction to terrain-following coordinates

Coordinate transformations Pressure coordinates Isentropic coordinates

Sigma coordinates

Smooth terrain-following coordinates

Schär, ETH Zürich

4

Terrain-following Sigma coordinates

σ = 0

σ = 1 p = ps p = 0 σ = const

Surface pressure:

€

p

_s

= p

_s

(x,y,t)

p = 0

p = p

_s

σ = 0

σ = 1 Special case: upper boundary at p _t = 0

€

σ := p p

_s

σ-coordinate:

General case: upper boundary at p _t = const

€

σ := p − p

_t

p

_s

− p

_t

σ-coordinate:

Schär, ETH Zürich

Terrain-following Sigma coordinates

5 10 15 20 25 [km]

5 10 [km]

Schär, ETH Zürich

6

1 30 38 46 5 1

1 21 28 30 1 0.1

Terrain-following hybrid coordinates

p=const

σ=const Hybrid coordinates:

• transition between two different coordinate formulations here: σ-coordinates at low levels, p-coordinates at high levels

• finer descretization at lower levels

Schär, ETH Zürich

Example of hybrid coordinate

Standard coordinate setting in the aLMo (COSMO) model of MeteoSwiss and DWD σ-based hybrid coordinate

15

10

5

z [km]

Schär, ETH Zürich

8

Outline

Introduction to terrain-following coordinates Coordinate transformations

Pressure coordinates Isentropic coordinates

Sigma coordinates

Smooth terrain-following coordinates

Schär, ETH Zürich

Basic considerations and objectives

physical space

(1) Transformation should simplify structure of governing equations:

(Example: continuity equation)

height-coordinates:

pressure-coordinates:

€

∂u

∂ x



  

 

p

+ ∂v

∂y



  

 

p

+ ∂ω

∂ p = 0

€

∂ρ

∂t + ∂ ( ) ρu

∂ x + ∂ ( ) ρv

∂ y + ∂ ( ρw )

∂z = 0

(nonlinear, prognostic) (linear, diagnostic)

(2) Transformation should map computational domain on a simple grid:

physical space

complex computational

domain

computational space

simple

grid

Schär, ETH Zürich

10

Implementation of coordinate transformations

€

u_j+1,k−u_j−1,k

2Δx +ω_j,k+1−ω_j,k−1

2Δp =0

Example: two-dimensional continuity equation Step 1:

Transform equations

Step 2:

Descretize equations

€

∂ρ

∂t + ∂ ( ) ρ u

∂ x + ∂ ( ρ w )

∂ z =0 height-coordinates

z = const

€

∂u

∂ x



  

 

p

+ ∂ω

∂ p =0 pressure-coordinates

p = const

descretized p-coordinates

j+1 j j–1

k+1 k k–1

Schär, ETH Zürich

Transformation of wind

€

u := Dx

Dt , v : = Dy

Dt , w := Dz Dt Wind vector in (x,y,z)-coordinates

€

u ˜ := D x ˜

Dt , ˜ v : = D y ˜

Dt , ˜ w := D˜ z Dt Wind vector in generalized coordinates

€

( ˜ x , ˜ y , ˜ z )

€

ω : = D p D t

Vertical wind in (x,y,p)-coordinates

[hPa/s] Note: ascending ω < 0, descending ω > 0

€

θ

^.

:= Dθ D t

Vertical wind in (x,y,θ)-coordinates

[K/s] Note: Dθ/Dt=0 for adiabatic flows

Schär, ETH Zürich

12

Transformation of advection operator

Advection in (x,y,z)-coordinates

€

D D t = ∂

∂ t + D x D t u {

∂

∂ x



  

 

y,z

+ D y D t v {

∂

∂ y



  

 

x,z

+ D z D t w {

∂

∂ z



  

 

x,y

€

D D t = ∂

∂ t + D x ˜ D t

˜ u {

∂

∂ x ˜



  

 

y ˜ ,˜ z

+ D y ˜ D t

˜ v {

∂

∂ y ˜



  

 

x ˜ ,˜ z

+ D z ˜ D t

˜ w {

∂

∂ z ˜



  

 

x ˜ , ˜ y

€

D D t = ∂

∂ t + ˜ u ∂

∂ x ˜



  

 

y ˜ ,˜ z

+ ˜ v ∂

∂ y ˜



 



 

x ˜ ,˜ z

+ ˜ w ∂

∂ z ˜



  

 

x ˜ , ˜ y

Advection in generalized coordinates

€

( ˜ x , ˜ y , ˜ z )

The advection operator is defined as the total derivative

The transformed advection

operator has the same

formal structure as its

Cartesian counterpart

Schär, ETH Zürich

€

∂ χ

∂ s

∂ s

∂ x



  

 

z

Transformation of derivatives

Transformation of z-coordinates to s-coordinates:

New coordinate defined by

s=s(x,y,z) or z=z(x,y,s) Invertibility implies strictly monotone functions

€

∂χ

∂z = ∂ χ

∂s

∂ s

∂z

Transformation of ∂χ/∂z

z = const s = const

€

∂ χ

∂ x



 



 

z

€

∂ χ

∂ x



  

 

s

€

∂ χ

∂ x



 



 

z

= ∂ χ

∂ x



 



 

s

+ ∂ χ

∂ s

∂ s

∂ x



 



 

z

Transformation of (∂χ/∂x) _z

Schär, ETH Zürich

14

Outline

Introduction to terrain-following coordinates Coordinate transformations

Pressure coordinates

Isentropic coordinates Sigma coordinates

Smooth terrain-following coordinates

Schär, ETH Zürich

Pressure coordinates

In pressure coordinates, all variables are expressed as χ=χ(x,y,p,t). For instance, the pressure field p(x,y,z,t) is described by

z (x,y,p,t) height [m]

z-coordinates

pressure at 3150 m

p-coordinates

height of 700 hPa surface

The geopotential is the potential of gravity:

€

g = ∂φ ∂ z

The geopotential height is a gravity-adjusted height accounting for variations of g=g(x,y,z).

€

z

_g

= φ g

_o

= 1

g

_o

g(z)dz

z=0

z

∫ ^g

_o

^{=9.81 m/s}

²

or φ (x,y,p,t) geopotential [m

²

/s

²

] or z _g (x,y,p,t) geopotential height [m]

Schär, ETH Zürich

18

€

∂u

∂ x



 



 

p

+ ∂v

∂ y



 



 

p

+ ∂ω

∂ p =0

€

DT D t = ω

ρc

_p

+ H

€

∂φ

∂ p = – 1 ρ

€

D Dt = ∂

∂t + u ∂

∂x



  

 

p

+v ∂ dy



 



 

p

+ω ∂

∂p

€

ω = Dp Dt

€

Dv

Dt + fu = – ∂φ

∂y



  

 

p

+ F

^y

€

Du

Dt – fv = – ∂φ

∂x



  

 

p

+ F

^x

Hydrostatic system in pressure coordinates

Horizontal momentum equations

Hydrostatic equation Equation of state

Thermodynamic equation Continuity equation

p = ρ R T

with and

F^x, F^y, H

Externally specified force and heating rate

φ

=g

_oz

Geopotential

Schär, ETH Zürich

Outline

Introduction to terrain-following coordinates Coordinate transformations

Pressure coordinates Isentropic coordinates

Sigma coordinates

Smooth terrain-following coordinates

Schär, ETH Zürich

21

Isentropic coordinates

€

θ = T p

_ref

p



 



 

R c_p

p_ref 1000 hPa reference pressure R 287 J/(K kg) gas constant for dry air

c_p 1004 J/(K kg) specific heat of dry air at constant pressure

Definition: The potential (or isentropic) temperature θ of an air parcel at pressure p is the temperature that the parcel would acquire if adiabatically brought to some reference pressure p _ref .

In isentropic coordinates, all variables are expressed as χ=χ(x,y,θ,t).

Disdvantage: The invertibility condition (i.e., existence of an inverse coordinate transformation) implies that isentropic coordinates are restricted to flows with ∂θ/∂z>0.

€

θ

.

: = D θ D t = 0

Advantage: Isentropic coordinates are particularly attractive for adiabatic flows, as the

vertical velocity (=diabatic heating rate measured in K/s) vanishes, i.e.

Schär, ETH Zürich

Example 1: gravity waves

(Doyle and Smith, 2003, QJ)

Height [m]

Vertically propagating gravity wave

W [m/s]

1 2 3 4 5 6 7

Trapped gravity wave

8

280 289 298 307 316 325

280 289 298 307 316 325

Idealized simulations showing

θ

and w (U=20m/s, from left to right)

Trapped gravity waves are non-hydrostatic features and not represented in the hydrostatic isentropic

framework

Schär, ETH Zürich

23

Example 2: Cold front

(Neiman, Shapiro, Fedor, 1993, MWR, 121)

Schär, ETH Zürich

Isentropic coordinates

θ = θs θ = θt

θ = const

θ = const θ = const

Schär, ETH Zürich

25

F^x, F^y, H

Externally specified force and heating rate

φ = g_oz

Geopotential

Montgomery potential €

∂σ

∂t + ∂(σu)

∂ x



  

 

θ

+ ∂(σv)

∂ y



  

 

θ

+ ∂(σ θ)

^.

∂θ = 0

€

c

_p

T θ = ∂M

∂θ

€

D Dt = ∂

∂t + u ∂

∂x



  

 

θ

+ v ∂

∂y



 



 

θ

+ θ

^.

∂

∂θ

€

Dv

Dt + fu = − ∂M

∂y



  

 

θ

+ F

^y

€

Du

Dt − fv = − ∂M

∂x



  

 

θ

+ F

^x

Hydrostatic system in isentropic coordinates

Horizontal momentum equations

Hydrostatic equation Equation of state

Thermodynamic equation Continuity equation

p = ρ R T with

φ = gz with

M = φ + c

_p

T €

σ := − 1 g

∂p

∂θ

€

Dθ

D t = θ

^.

Schär, ETH Zürich

Grid structure in θ -coordinates

k = K+1/2

k = 1/2

vertical section

k = K–1/2 k = K–3/2

k = 3/2 k = K k = K–1

k = 1

horizontal section

j–2 j–1/2

j–1 j+1/2

j j+3/2

j+1

i–1 i–1/2

i i+1/2

i+1 i+3/2

i+2

u v v, σ, M u

M, σ, θ, θ, p

.

θ = θs+(k–1/2)Δθ

θ, θ, p

.

θ = θs

θt = θs+(k–1/2)Δθ

=θs+KΔθ

k = 5/2 k = 2

Schär, ETH Zürich

29

Simplified system in isentropic coordinates

Simplifications: ∂/∂y = 0, D/Dθ = 0, F = 0, f = 0

θ = θ_s θ = θ_t

θ = const

€

θ = θ

_t

p

_t

= const

θ = θ

_s

M

_s

= c

_p

T + gz

_s

Boundary conditions:

• upper boundary is horizontal isentropic surface

• lower boundary is isentropic surface

€

∂σ

∂t + ∂(σu)

∂ x



  

 

θ

= 0

€

Du Dt = − ∂M

∂x



  

 

θ

€

π = ∂M

∂θ

€

D Dt = ∂

∂t + u ∂

∂x



  

 

θ

Momentum equations

Continuity equation

Hydrostatic equation

with

with

€

σ = − 1 g

∂p

∂θ

€

π = c

_p

T

θ = c

_p

p p

_ref





 





 

R c_p

with

Exner function

Schär, ETH Zürich

Step 3: Diagnosis of Montgomery potential => M ⁿ⁺¹

Step 2: Diagnosis of pressure and Exner function => p ⁿ⁺¹ , π ⁿ⁺¹ Integration

Step 1: Integration of momentum and continuity equations => u ⁿ⁺¹ , σ ⁿ⁺¹

€

∂σ

∂t + ∂(σu)

∂ x



  

 

θ

= 0

€

Du Dt = − ∂M

∂x



  

 

θ

€

π = ∂M

∂θ

€

D Dt = ∂

∂t + u ∂

∂x



  

 

θ

Momentum equations

Continuity equation

Hydrostatic equation

with

with

€

σ = − 1 g

∂p

∂θ

€

π( p) = T

θ = c

_p

p p

_ref





 





 

R c_p

with

Simplified system in isentropic coordinates

Simplifications: ∂/∂y = 0, D/Dθ = 0, F = 0, f = 0

Schär, ETH Zürich

31

k = K+1/2

k = 1/2

vertical section

k = K–1/2

k = 3/2 k = K k = K–1

k = 1

u σ, M

θ = θ_s + (k–1/2)Δθ

p

θ = θs

θt = θs + KΔθ θ = θs + (K–1/2)Δθ

θ = θs + (1/2)Δθ θ = θ_s + (1/2)Δθ θ = θs + (K–1)Δθ

k

k+1/2 θ = θs + kΔθ

i i –1/2 i+1/2

Grid structure with simplified θ -coordinates

Primary model variables

σ_i,k (k=1,.. K)

isentropic density

M_i,k (k=1,.. K)

Montgomery potential

u_i+1/2,k (k=1,.. K)

horizontal velocity

p_i,k+1/2 (k=0,.. K)

pressure

Schär, ETH Zürich

Step 1: Integration of momentum and continuity equations

Introduce vertical descretization: σ

_k

, M

_k

, u

_k

For each level k, the resulting system is formally identical with the shallow-water equations.

€

∂u

_k

∂t + u

_k

∂u

_k

∂x = − ∂M

_k

∂x



  

 

θ

€

∂σ

_k

∂t + ∂ (σ

_k

u

_k

)

∂ x



  

 

θ

= 0

The isentropic system can be interpreted as a stack of shallow-water layers. The integration of the momentum and continuity equations employs the same methods as those in a shallow-water model. After a time-step, the variables u

_i,k

and

σ_i,k

will be known at time t+

Δt.

€

Du

Dt + ∂ ( h + H )

∂x = 0

€

∂H

∂t + ∂(uH )

∂x = 0

Isentropic system Shallow-water system

Schär, ETH Zürich

33

Step 2: Diagnosis of pressure and Exner function

Diagnose p

ⁿ⁺¹

from σ

ⁿ⁺¹

^:

€

σ = − 1 g

∂p

∂θ

€

σ

_iⁿ⁺¹_,k

= − 1 g

p

_{i,k+1 / 2}ⁿ⁺¹

− p

_{i,k−1 / 2}ⁿ⁺¹

Δθ

π = c

_p

T

θ = c

_p

p p

_ref





 





 

R c_p

Compute π

ⁿ⁺¹

^from p

ⁿ⁺¹

:

Integrate from top to bottom

€

p

_{i,k−1 / 2}ⁿ⁺¹

= p

_{i,k+1 / 2}ⁿ⁺¹

+ gΔθσ

_i,kⁿ⁺¹

σ

_k

p

_k+1/2

p

_k–1/2

Requires upper b.c.

p

_K+1/2

p

_K+1/2

Schär, ETH Zürich

€

π

_{i,k−1 2}

= M

_i,k

− M

_i,k−1

Δθ

€

π = ∂M

∂θ

Step 3: Diagnosis of Montgomery potential

Diagnose M

ⁿ⁺¹

from π

ⁿ⁺¹

:

Requires lower b.c.

M

_k=1

M

_k=1

€

M

_i,kⁿ⁺¹

= M

_i,k−1ⁿ⁺¹

+ Δθπ

_{i,k–1/ 2}ⁿ⁺¹ Integrate from bottom to top

π

_k–1/2

M

_k

M

_k–1

Schär, ETH Zürich

35

Outline

Introduction to terrain-following coordinates Coordinate transformations

Pressure coordinates Isentropic coordinates

Sigma coordinates

Smooth terrain-following coordinates

Schär, ETH Zürich

€

∂u

∂ x



 



 

p

+ ∂v

∂ y



 



 

p

+ ∂ω

∂ p =0

€

DT D t = ω

ρ c

_p

+ H

€

∂φ

∂ p = – 1 ρ

€

D Dt = ∂

∂t + u ∂

∂x



  

 

p

+v ∂ dy



  

 

p

+ω ∂

∂p

€

ω = Dp Dt

€

Dv

Dt + fu = – ∂φ

∂y



 



 

p

+ F

^y

€

Du

Dt – fv = – ∂φ

∂x



  

 

p

+ F

^x

Hydrostatic system in pressure coordinates

Horizontal momentum equations

Hydrostatic equation Equation of state

Thermodynamic equation Continuity equation

p = ρ R T

with and

€

p = 0 : ω = 0

€

p = p

_s

: ω = ∂p

_s

∂t +v

_s

⋅∇

_h

p

_s

Boundary conditions at top:

Boundary conditions at surface:

€

φ= g z

_s

Schär, ETH Zürich

37

€

∂ p

_s

∂t +∇

_σ

⋅ ( p

_s

v ) ^{+ p}

s

∂ σ ˙

∂σ =0

€

DT Dt = RT

c

_p

σp

_s

ω + H

€

σ p

_s

= ρ RT

€

∂φ

∂ σ = − RT σ

€

Du

Dt − fv = − ∂φ

∂x



  

 

σ

− RT p

_s

∂p

_s

∂x + F

^x

€

D Dt = ∂

∂t + u ∂

∂x



  

 

σ

+ v ∂ dy



  

 

σ

+ σ ˙ ∂

∂σ

€

σ = ˙ D σ Dt

Hydrostatic system in sigma coordinates

Horizontal momentum equations

Hydrostatic equation Equation of state

Thermodynamic equation Continuity equation

with and

€

ω = σ ∂p

_s

∂t + u ∂p

_s

∂x + v ∂p

_s

∂y



 



  + σ ˙ p

_s

€

p = 0 : σ = 0, ˙ σ = 0 p = p

_s

: σ =1, ˙ σ = 0 Boundary conditions at top:

Boundary conditions at surface: φ= g z

_s

Schär, ETH Zürich

k = 1/2

k = K+1/2 p = ps, σ = 1

p = 0, σ = 0

vertical section σ = p / p_s

k = 3/2 k = 5/2

k = K–1/2 k = 7/2

k =1 k =2 k =3

k = K horizontal section

j–2 j–1/2

j–1 j+1/2

j j+3/2

j+1

i–1 i–1/2

i i+1/2

i+1 i+3/2

i+2

u v T, v u

Τ, φ, σ

.

, ω _{φ, σ}

.

_{, ω}

σ = (k-1/2) Δσ

Grid structure in σ -coordinates

Schär, ETH Zürich

39

Vertical wind and surface pressure tendency

If the variables u, v, T, p

_s

are known, the vertical wind may be diagnosed Integrate continuity equation from σ’=0 to σ’=σ

Solve for vertical wind:

Prognosis of surface pressure tendency: Integrate from σ’=0 to σ’=1

€

σ ˙

€

σ ∂p

_s

∂t + ∇

_σ

⋅ ( p

_s

v )

0 σ

∫ ^{d σ + ′ p}

s

σ ˙

_σ

= 0

€

σ ˙ (σ ) = − σ p

_s

∂p

_s

∂t − 1

p

_s

∇

_σ

⋅ ( p

_s

v )

0 σ

∫ ^{d σ ′}

€

∂ p

_s

∂t = − ∇

_σ

⋅ ( p

_s

v )

0 1

∫ ^{d σ}

€

∂ p

_s

∂t + ∇

_σ

⋅ ( p

_s

v ) ^{+ p}

s

∂ σ ˙

∂σ =0 =>

Schär, ETH Zürich

Integration in sigma-coordinates

Variables u, v, T, p

_s

and φ known at time t

• Computation of surface pressure tendency ∂p

_s

/∂t

• Diagnosis of vertical wind

• Computation of tendencies ∂u/∂t, ∂v/∂t and ∂T/∂t

• Time step for u, v, T and p

_s

to time t+Δt

• Diagnosis of gepotential φ

Variables u, v, T, p

_s

and φ known at time t+Δt

€

σ ˙

Schär, ETH Zürich

41

Outline

Introduction to terrain-following coordinates Coordinate transformations

Pressure coordinates Isentropic coordinates

Sigma coordinates

Smooth terrain-following coordinates

Schär, ETH Zürich

Alpenpanorma (Bella Tola)

Real topography is extremely complex.

This implies awkward terrain-following grids.

How will this affect quality of numerical solutions?

=> Compare different coordinate formulations! <=

Schär, ETH Zürich

44

Sigma coordinates

5 10 15 20 25 [km]

5 10 [km]

Schär, ETH Zürich

Hybrid coordinates

5 10 15 20 25 [km]

5 10 [km]

Schär, ETH Zürich

46

Smooth Level Vertical (SLEVE) coordinate

5 10 15 20 25 [km]

5 10 [km]

Scale-dependent decay of terrain-features with height

Schär et al. 2003

Schär, ETH Zürich

Tested numerical schemes:

• Centered in space and time, 2nd order

• Centered in space and time, 4th order

• Upstream, 1st order

• MPDATA (Smolarkiewicz), 2nd order Numerics in conservative flux form

Idealized Advection Test

Tested coordinates:

• Sigma

• Hybrid

• SLEVE

• Height (reference)

8Δx Topography:

wave-length of 8Δx

Schär, ETH Zürich

48

Coordinates

Sigma Hybrid

SLEVE Reference

Schär, ETH Zürich

Error Fields of idealized advection test

Hybrid

SLEVE Reference

Sigma

Schär, ETH Zürich

52

Interpretation (for 2nd-order centered scheme)

Truncation error for 1D linear advection equation: ^∂ρ

∂t + ∂ ( ) u ρ

∂x = 0 (a) regular grid:

Δx x

(b) general grid:

Δx x

Δx

Schär, ETH Zürich

Interpretation (for 2nd-order centered scheme)

(a) regular grid:

Truncation error E

_reg

= Δx

²

6 u ∂

³

ρ

∂x

³

+ O ( Δx

³

) 2nd-order scheme

dominates in presence of small-scale anomalies

dominates in presence of small-scale grid deformations (b) general grid:

Transformed equation with Jacobian Truncation error

J = ∂X

∂x ≈ ΔX Δx

∂ρ

∂t + J ∂ ( ) u ρ

∂X = 0

E

_trans

= Δx

²

6 u ∂

³

ρ

∂x

³

+ J

²

Δx

²

12 u ∂

²

J

⁻²

∂x

²

∂ρ

∂x + 3 ∂J

⁻²

∂x

∂

²

ρ

∂x

²



 



  + O ( Δx

³

)

have same leading order

Truncation error for 1D linear advection equation: ^∂ρ

∂t + ∂ ( ) uρ

∂x = 0

Schär, ETH Zürich

54

Idealized Flow Test

Nonhydrostatic flow past small-amplitude topography Analytical solution

propagating hydrostatic wave

decaying nonhydrostatic wave

Schär, ETH Zürich

Idealized Flow Tests

SigmaSLEVE

MC2, N=0.01s

^-1

MC2, isothermal LM, N=0.01s

^-1

Schär et al. 2003, Klemp et al. 2003