Coastal Modeling

Report
Regional Coastal Ocean
Modeling
Patrick Marchesiello
Brest, 2005
Patrick Marchesiello IRD 2005
1
The Coastal Ocean:
A Challenging Environment





Geometrical constraints: irregular coastlines and highly
variable bathymetry
Forcing is internal (intrinsic), lateral and superficial: tides,
winds, buoyancy
Broad range of space/time scales of coastal structures
and dynamics: fronts, intense currents, coastal trapped
waves, (sub)mesoscale variability, turbulent mixing in
surface and bottom boundary layers
Heterogeneity of regional and local characteristics:
eastern/western boundary systems; regions can be
dominated by tides, opened/closed to deep ocean
Complexe Physical-biogeochemical interactions
Patrick Marchesiello IRD 2005
2
Numerical Modeling



Require highly optimized models of
significant dynamical complexity
In the past: simplified models due to limited
computer resources
In recent years: based on fully nonlinear
stratified Primitive Equations
Patrick Marchesiello IRD 2005
3
Coastal Model Inventory









POM
ROMS
MARS3D
SYMPHONIE
GHERM
HAMSOM
QUODDY
MOG3D
SEOM
Finite-Difference Models
Finite-Elements Models
Patrick Marchesiello IRD 2005
4
Patrick Marchesiello IRD 2005
5
Hydrodynamics
Patrick Marchesiello IRD 2005
6
Primitive
Equations:
Hydrostatic,
Incompressible,
Boussinesq
Momentum
Tracer
Hydrostatic
Similar transport equations
for other tracers: passive or
actives
Patrick Marchesiello IRD 2005
Continuity
7
Vertical Coordinate System

Bottom following
coordinate (sigma):
best representation of
bottom dynamics:

but subject to pressure
gradient errors on steep
bathymetry
Patrick Marchesiello IRD 2005
8
GENERALIZED
-COORDINATE
Stretching &
condensing of
vertical resolution
(a)
(b)
(c)
(d)
Ts=0, Tb=0
Ts=8, Tb=0
Ts=8, Tb=1
Ts=5, Tb=0.4
Patrick Marchesiello IRD 2005
9
Horizontal Coordinate System

Orthogonal curvilinear
coordinates
Patrick Marchesiello IRD 2005
10
Primitive
Equations in
Curvilinear
Coordinate
Patrick Marchesiello IRD 2005
11
Simplified Equations

2D barotropic


2D vertical


Tidal problems
Upwelling
1D vertical

Turbulent mixing problems (with
boundary layer parameterization)
Patrick Marchesiello IRD 2005
12
Barotropic Equations
Patrick Marchesiello IRD 2005
13
Vertical Problems:
Parameterization of Surface and
Bottom Boundary Layers
Patrick Marchesiello IRD 2005
14
Boundary Layer
Parameterization

Boundary layers are
characterized by strong
turbulent mixing

Turbulent Mixing depends on:

Surface/bottom forcing:



Wind / bottom-shear stress stirring
Stable/unstable buoyancy forcing
Reynolds term:
w’T’
Interior conditions:


Current shear instability
Stratification
Patrick Marchesiello IRD 2005
15
Surface and Bottom Forcing
Wind stress
Heat Flux
Salt Flux
Bottom stress
Drag Coefficient CD:
γ1=3.10-4 m/s Linear
γ2=2.5 10-3 Quadratic
Patrick Marchesiello IRD 2005
16
Boundary Layer Parameterization


All mixed layer schemes are based on
one-dimentional « column physics »
Boundary layer parameterizations are
based either on:



Turbulent closure (Mellor-Yamada, TKE)
K profile (KPP)
Note: Hydrostatic stability may require
large vertical diffusivities:


implicit numerical methods are best suited.
convective adjustment methods (infinite
diffusivity) for explicit methods
Patrick Marchesiello IRD 2005
17
Application: Tidal Fronts
ROMS
Simulation in
the Iroise
Sea (Front
d’Ouessant)
Simpson-Hunter
criterium for tidal
fronts position
1.5 <
<2
H. Muller, 2004
Patrick Marchesiello IRD 2005
18
Bottom Shear Stress – Wave effect
 c   /ln(za /z0 ) u2z z
2
 w  12 f w ub2 ;

a
f w  1.39(ub /z0 )0.52
Waves enhance bottom shear stress (Soulsby 1995):
3.2

  w  

 cw   c 1 1.2

 c   w  

Patrick Marchesiello IRD 2005
19
Numerical Discretization
Patrick Marchesiello IRD 2005
20
A Discrete Ocean
Patrick Marchesiello IRD 2005
21
Structured / Unstructured Grids
Finite Differences / Elements


Structured grids: the grid cells have the same
number of sides
Unstructured grids: the domain is tiled using more
general geometrical shapes (triangles, …) pieced
together to optimally fit details of the geometry


Good for tidal modeling, engineering applications
Problems: geostrophic balance accuracy, wave scattering
by non-uniform grids, conservation and positivity
properties, …
Patrick Marchesiello IRD 2005
22
Finite Difference (Grid Point) Method

If we know:



The ocean state at time t (u,v,w,T,S, …)
Boundary conditions (surface, bottom, lateral
sides)
We can compute the ocean state at t+dt
using numerical approximations of Primitive
Equations
Patrick Marchesiello IRD 2005
23
Horizontal and Vertical Grids
Patrick Marchesiello IRD 2005
24
Consistent Schemes:
Taylor series expansion, truncation errors


We need to find an consistent approximation
for the equations derivatives
Taylor series expansion of f at point x:
Truncation
error
Patrick Marchesiello IRD 2005
25
Exemple: Advection Equation
Dx
Dt
grid space
time step
Dt
Dx
Patrick Marchesiello IRD 2005
26
Order of Accuracy
First order
Downstream
Upstream
2nd order
Centered
4th order
Patrick Marchesiello IRD 2005
27
Numerical properties:
stability, dispersion/diffusion
Leapfrog / Centered
Tin+1 = Ti n-1 - C (Ti+1n - Ti-1n) ; C = u0 dt / dx
Conditionally stable: CFL condition C <
1 but dispersive (computational modes)

Euler / Centered
Tin+1 = Ti n - C (Ti+1n - Ti-1n)
Unconditionally unstable

Upstream
Tin+1 = Ti n - C (Tin - Ti-1n) , C > 0
Tin+1 = Ti n - C (Ti+1n - Tin) , C < 0
Conditionally stable,
not dispersive but diffusive
(monotone linear scheme)
Advection equation:
should be non-dispersive:the
phase speed ω/k and group
speed δω/δk are equal and
constant (uo)

Patrick Marchesiello IRD 2005
2nd order approx to the
modified equation:
28
Numerical Properties
A numerical scheme can be:
• Dispersive: ripples, overshoot and
extrema (centered)
• Diffusive (upstream)
• Unstable (Euler/centered)
Patrick Marchesiello IRD 2005
29
Weakly Dispersive, Weakly Diffusive Schemes

Using high order upstream schemes:


Using a right combination of a centered scheme
and a diffusive upstream scheme



3rd order upstream biased
TVD, FCT, QUICK, MPDATA, UTOPIA, PPM
Using flux limiters to build nolinear monotone
schemes and guarantee positivity and
monotonicity for tracers and avoid false extrema
(FCT, TVD)
Note: order of accuracy does not reduce
dispersion of shorter waves
Patrick Marchesiello IRD 2005
30
Upstream
Centered
2nd order
flux limited
3rd order
flux limited
Durran, 2004
Patrick Marchesiello IRD 2005
31
Accuracy
Numerical dispersion
2nd order
High order accurate methods:
optimal choice (lower cost for a given
accuracy) for general ocean circulation
models is 3RD OR 4TH ORDER accurate
methods (Sanderson, 1998)
With special care to:
• dispersion / diffusion
• monotonicity and positivity
• Combination of methods
4th order
2nd order
double resolution
Spectral method
Patrick Marchesiello IRD 2005
32
Sensitivity to the Methods: Example
OPA - 0.25 deg
ROMS – 0.25 deg
C. Blanc
Patrick Marchesiello IRD 2005
C. Blanc
33
Properties of Horizontal
Grids
Patrick Marchesiello IRD 2005
34
Arakawa Staggered Grids
Linear shallow water equation:
A staggered difference is 4 times
more accurate than non-staggered
and improves the dispersion
relation because of reduced use of
averaging operators

Patrick Marchesiello IRD 2005
35
Horizontal Arakawa grids

B grid is prefered at coarse resolution:




C grid is prefered at fine resolution:




Superior for poorly resolved inertia-gravity waves.
Good for Rossby waves: collocation of velocity points.
Bad for gravity waves: computational checkboard mode.
Superior for gravity waves.
Good for well resolved inertia-gravity waves.
Bad for poorly resolved waves: Rossby waves
(computational checkboard mode) and inertia-gravity
waves due to averaging the Coriolis force.
Combinations can also be used (A + C)
Patrick Marchesiello IRD 2005
36
Arakawa-C Grid
Patrick Marchesiello IRD 2005
37
Vertical
Staggered
Grid
Patrick Marchesiello IRD 2005
38
Numerical Round-off Errors
Patrick Marchesiello IRD 2005
39
Round-off Errors







Round-off errors result from inability of computers to represent a
floating point number to infinite precision.
Round-off errors tend to accumulate but little control on the
magnitude of cumulative errors is possible.
1byte=8bits, ex:10100100
Simple precision machine (32-bit):
1 word=4 bytes, 6 significant digits
Double precision machine (64-bit):
1 word=8 bytes, 15 significant digits
Accuracy depends on word length and fractions assigned to
mantissa and exponent.
Double precision is possible on a machine of any given basic
precision (using software instructions), but penalty is: slowdown
in computation.
Patrick Marchesiello IRD 2005
40
Time Stepping
Patrick Marchesiello IRD 2005
41
Time Stepping: Standard

Leapfrog: φin+1 = φi n-1 + 2 Δt F(φin)


computational mode amplifies when applied to
nonlinear equations (Burger, PE)
Leapfrog + Asselin-Robert filter:
φin+1 = φfi n-1 + 2 Δt F(φin)
φfi n = φi n + 0.5 α (φin+1 - 2 φin + φfin-1)
 reduction of accuracy to 1rst order depending on
α (usually 0.1)
Patrick Marchesiello IRD 2005
42
Time Stepping: Performance
C = 0.5
C = 0.2
Kantha and Clayson (2000) after Durran (1991)
Patrick Marchesiello IRD 2005
43
Time Stepping: New Standards

Multi-time level schemes:
Multi-time level
scheme



Adams-Bashforth 3rd order (AB3)
Adams-Moulton 3rd order (AM3)
Multi-stage Predictor/Corrector scheme
Multi-stage
scheme



Increase of robustness and stability range
LF-Trapezoidal, LF-AM3, Forward-Backward
Runge-Kutta 4: best but expensive
Patrick Marchesiello IRD 2005
44
Barotropic Dynamics
and Time Splitting
Patrick Marchesiello IRD 2005
45
Time step restrictions




The Courant-Friedrichs-Levy CFL stability condition on the
barotropic (external) fast mode limits the time step:
Δtext < Δx / Cext where Cext = √gH + Uemax
ex: H =4000 m, Cext = 200 m/s (700 km/h)
Δx = 1 km, Δtext < 5 s
Baroclinic (internal) slow mode:
Cin ~ 2 m/s + Uimax (internal gravity wave phase speed + max
advective velocity)
Δx = 1 km, Δtext < 8 mn
Δtin / Δtext ~ 60-100 !
Additional diffusion and rotational conditions:
Δtin < Δx2 / 2 Ah and Δtin < 1 / f
Patrick Marchesiello IRD 2005
46
Barotropic Dynamics


The fastest mode (barotropic) imposes a
short time step
3 methods for releasing the time-step
constraint:




Rigid-lid approximation
Implicit time-stepping
Explicit time-spitting of barotropic and baroclinic
modes
Note: depth-averaged flow is an
approximation of the fast mode (exactly true
only for gravity waves in a flat bottom ocean)
Patrick Marchesiello IRD 2005
47
Rigid-lid Streamfunction Method


Advantage: fast mode is properly filtered
Disadvantages:


Preclude direct incorporation of tidal processes, storm
surges, surface gravity waves.
Elliptic problem to solve:

convergence is difficult with complexe geometry; numerical
instabilities near regions of steep slope (smoothing required)
Matrix inversion (expensive for large matrices); Bad
scaling properties on parallel machines
Fresh water input difficult
Distorts dispersion relation for Rossby waves



Patrick Marchesiello IRD 2005
48
Implicit Free Surface Method


Numerical damping to supress barotropic
waves
Disadvantanges:


Not really adapted to tidal processes unless Δt is
reduced, then optimality is lost
Involves an elliptic problem


matrix inversion
Bad parallelization performances
Patrick Marchesiello IRD 2005
49
Time Splitting
Explicit free surface method
Patrick Marchesiello IRD 2005
50
Barotropic Dynamics:Time Splitting

Direct integration of barotropic equations, only few
assumptions; competitive with previous methods at
high resolution (avoid penalty on elliptic solver);
good parallelization performances

Disadvantages: potential instability issues involving
difficulty of cleanly separating fast and slow modes

Solution:


time averaging over the barotropic sub-cycle
finer mode coupling
Patrick Marchesiello IRD 2005
51
Time Splitting: Averaging
Averaging
weights
ROMS
Patrick Marchesiello IRD 2005
52
Time Splitting: Coupling terms
Coupling terms: advection (dispersion) + baroclinic PGF
Patrick Marchesiello IRD 2005
53
Internal
mode
Flow
Diagram
of POM
External
mode
Forcing terms
of external mode
Replace barotropic part
in internal mode
Patrick Marchesiello IRD 2005
54
Vertical Diffusion
Patrick Marchesiello IRD 2005
55
Vertical
Diffusion
Semiimplicit
CrankNicholson
scheme
Patrick Marchesiello IRD 2005
56
Pressure Gradient Force
Patrick Marchesiello IRD 2005
57
PGF Problem

Truncation errors are made
from calculating the baroclinic
pressure gradients across
sharp topographic changes
such as the continental slope

Difference between 2 large
terms

Errors can appear in the
unforced flat stratification
experiment
Patrick Marchesiello IRD 2005
58
Reducing PGF Truncation Errors

Smoothing the topography using a
nonlinear filter and a criterium:

Using a density formulation

Using high order schemes to reduce
the truncation error (4th order,
McCalpin, 1994)

Gary, 1973: substracting a reference
horizontal averaged value from
density (ρ’= ρ - ρa) before computing
pressure gradient
Rewritting Equation of State: reduce
passive compressibility effects on
pressure gradient

Patrick Marchesiello IRD 2005
r = Δh / h < 0.2
59
Equation of State
Full UNESCO EOS:
30% of total CPU!
Jackett &
McDougall,
1995: 10%
of CPU
Patrick Marchesiello IRD 2005
Linearization
(ROMS):
reduces PGF
errors
60
Smoothing methods


r = Δh / h is the slope of the logarithm of h
One method (ROMS) consists of smoothing ln(h)
until r < rmax
Res: 1 km
r < 0.25
Res: 5 km
r < 0.25
Senegal Bathymetry Profil
Patrick Marchesiello IRD 2005
61
Smoothing method and resolution
Standard Deviation [m]
Bathymetry Smoothing Error off Senegal
Convergence at
~ 4 km resolution
Grid Resolution [deg]
Patrick Marchesiello IRD 2005
62
Errors in Bathymetry data compilations
Gebco1 compilation
Etopo2: Satellite observations
Shelf errors
(noise)
Patrick Marchesiello IRD 2005
63
Wetting and Drying Schemes
Patrick Marchesiello IRD 2005
64
Wetting and Drying: Principles

Application:



Principles:



mask/unmask
drying/wetting
areas at every time
step
Criterium based on
a minimum depth
Requirements

Patrick Marchesiello IRD 2005
Intertidal zone
Storm surges
Conservation
properties
65

similar documents