670 likes | 952 Views
Regional Coastal Ocean Modeling. Patrick Marchesiello Brest, 2005. The Coastal Ocean: A Challenging Environment. Geometrical constraints : irregular coastlines and highly variable bathymetry Forcing is internal (intrinsic), lateral and superficial: tides, winds, buoyancy
E N D
Regional Coastal Ocean Modeling Patrick Marchesiello Brest, 2005 Patrick Marchesiello IRD 2005
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
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
Coastal Model Inventory • POM • ROMS • MARS3D • SYMPHONIE • GHERM • HAMSOM • QUODDY • MOG3D • SEOM Finite-Difference Models Finite-Elements Models Patrick Marchesiello IRD 2005
Hydrodynamics Patrick Marchesiello IRD 2005
Momentum Primitive Equations:Hydrostatic, Incompressible,Boussinesq Tracer Hydrostatic Similar transport equations for other tracers: passive or actives Continuity Patrick Marchesiello IRD 2005
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
GENERALIZED -COORDINATE Stretching & condensing of vertical resolution • Ts=0, Tb=0 • Ts=8, Tb=0 • Ts=8, Tb=1 • Ts=5, Tb=0.4 Patrick Marchesiello IRD 2005
Horizontal Coordinate System • Orthogonal curvilinear coordinates Patrick Marchesiello IRD 2005
Primitive Equations in Curvilinear Coordinate Patrick Marchesiello IRD 2005
Simplified Equations • 2D barotropic • Tidal problems • 2D vertical • Upwelling • 1D vertical • Turbulent mixing problems (with boundary layer parameterization) Patrick Marchesiello IRD 2005
Barotropic Equations Patrick Marchesiello IRD 2005
Vertical Problems:Parameterization of Surface and Bottom Boundary Layers Patrick Marchesiello IRD 2005
w’T’ 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 • Interior conditions: • Current shear instability • Stratification Reynolds term: K theory Patrick Marchesiello IRD 2005
Bottom stress Drag Coefficient CD: γ1=3.10-4 m/s Linear γ2=2.5 10-3 Quadratic Surface and Bottom Forcing Wind stress Heat Flux Salt Flux Patrick Marchesiello IRD 2005
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
Simpson-Hunter criterium for tidal fronts position 1.5 < < 2 Application: Tidal Fronts ROMS Simulation in the Iroise Sea (Front d’Ouessant) H. Muller, 2004 Patrick Marchesiello IRD 2005
Bottom Shear Stress – Wave effect • Waves enhance bottom shear stress (Soulsby 1995): Patrick Marchesiello IRD 2005
Numerical Discretization Patrick Marchesiello IRD 2005
A Discrete Ocean Patrick Marchesiello IRD 2005
Structured / Unstructured GridsFinite 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
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
Horizontal and Vertical Grids Patrick Marchesiello IRD 2005
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
Exemple: Advection Equation Dx grid space Dt time step Dt Dx Patrick Marchesiello IRD 2005
Order of Accuracy First order Downstream Upstream 2nd order Centered 4th order Patrick Marchesiello IRD 2005
Numerical properties: stability, dispersion/diffusion • Leapfrog / Centered Tin+1 = Tin-1 - C (Ti+1n - Ti-1n) ; C = u0 dt / dx Conditionally stable: CFL condition C < 1 but dispersive (computational modes) • Euler / Centered Tin+1 = Tin - C (Ti+1n - Ti-1n) Unconditionally unstable • Upstream Tin+1 = Tin - C (Tin - Ti-1n) , C > 0 Tin+1 = Tin - 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) 2nd order approx to the modified equation: Patrick Marchesiello IRD 2005
Numerical Properties • A numerical scheme can be: • Dispersive: ripples, overshoot and extrema (centered) • Diffusive (upstream) • Unstable (Euler/centered) Patrick Marchesiello IRD 2005
Weakly Dispersive, Weakly Diffusive Schemes • Using high order upstream schemes: • 3rd order upstream biased • Using a right combination of a centered scheme and a diffusive upstream scheme • 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
Upstream Centered 2nd order flux limited 3rd order flux limited Durran, 2004 Patrick Marchesiello IRD 2005
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 4THORDER 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
Sensitivity to the Methods: Example ROMS – 0.25 deg OPA - 0.25 deg C. Blanc C. Blanc Patrick Marchesiello IRD 2005
Properties of Horizontal Grids Patrick Marchesiello IRD 2005
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
Horizontal Arakawa grids • B grid is prefered at coarse resolution: • Superior for poorly resolved inertia-gravity waves. • Good for Rossby waves: collocation of velocity points. • Bad for gravity waves: computational checkboard mode. • C grid is prefered at fine resolution: • 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
Arakawa-C Grid Patrick Marchesiello IRD 2005
Vertical Staggered Grid Patrick Marchesiello IRD 2005
Numerical Round-off Errors Patrick Marchesiello IRD 2005
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
Time Stepping Patrick Marchesiello IRD 2005
Time Stepping: Standard • Leapfrog: φin+1 = φin-1 + 2 Δt F(φin) • computational mode amplifies when applied to nonlinear equations (Burger, PE) • Leapfrog + Asselin-Robert filter: φin+1 = φfin-1 + 2 Δt F(φin) φfin = φin + 0.5 α(φin+1 - 2 φin + φfin-1) • reduction of accuracy to 1rst order depending on α (usually 0.1) Patrick Marchesiello IRD 2005
Time Stepping: Performance C = 0.5 C = 0.2 Kantha and Clayson (2000) after Durran (1991) Patrick Marchesiello IRD 2005
Time Stepping: New Standards • Multi-time level schemes: • Adams-Bashforth 3rd order (AB3) • Adams-Moulton 3rd order (AM3) • Multi-stage Predictor/Corrector scheme • Increase of robustness and stability range • LF-Trapezoidal, LF-AM3, Forward-Backward • Runge-Kutta 4: best but expensive Multi-time level scheme Multi-stage scheme Patrick Marchesiello IRD 2005
Barotropic Dynamicsand Time Splitting Patrick Marchesiello IRD 2005
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
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
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
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
Time SplittingExplicit free surface method Patrick Marchesiello IRD 2005