-SELFE Users Trainng Course- (2) SELFE numerical formulation Joseph Zhang, CMOP, OHSU

1. -SELFE Users Trainng Course-(2)SELFE numerical formulationJoseph Zhang, CMOP, OHSU

2. Dx x j j-1 j+1 N 1 Numerical methods for pde • Finite difference • Approximation of the differential operators using Taylor series • System of algebraic equations • Truncation error consistency, and order of convergence (Lax) • Stability • Higher-order approximation • Upwind • Finite element u Residual Basis/shape function Weighting function (satisfied at finite number of points) Galerkin FE Positive definite operator  positive definite matrix Symmetric mass matrix Least square Strong vs. weak formulation

3. x j j-1 j+1 N 1 x j j-1 j+1 N 1 j Shape function • Used to approximate the unknown function • Although usually used within each individual elements, shape function is global not local! • Must be sufficiently smooth to allow integration by part in weak formulation • Mapping between local and global coordinates • Assembly of global matrix Conformal Non-conformal

4. Semi-implicit Eulerian-Lagrangian Finite Element (SELFE) • A formal semi-implicit finite-element framework • A key step is to decouple the continuity and momentum equations through the bottom boundary layer • Unstructured grid in the horizontal • Hybrid SZ in the vertical • Eulerian-Lagrangian method (ELM) for momentum advection • Volume conservation is not enforced numerically (but good) • Numerical efficiency: semi-implicit time stepping; ELM • Mild stability constraint: horizontal viscosity, baroclinic pressure • Moderately more expensive than ELCIRC • Convergence rate: 2nd order for uniform grid • Optional higher-order numerics for momentum and transport • Treats inundation/drying in a natural way (validated for inundation benchmarks) • Bed deformation scheme incorporated (e.g., tsunami)

5. Horizontal grid (i,j)

6. Vertical grid (1) Pure S zone z SZ zone s zone S zone s=0 Nz hc hs hs S-levels kz s=-1 kz-1 Z-levels 2 1

7. S-coordinates qf=10-3, qb=1 qf=10, qb=1 hc hs qf=10, qb=0 qf=10, qb=0.7

8. z=-hs Vertical grid (2)

9. nk+1 k+1 k s2 z2 s1 z1 Vertical grid (3) • 3D computationa unit: uneven prisms • S-coordinates (Song and Haidvogel 1994) are ideal for shallow region • Z-coordinates are necessary to “stabilize” the deep region • s-coordinates are used where S-coordinates are invalid • Equations are not transformed into S- or s-coordinates, but solved in their original (Z) form • Interpolation mode: along Z or S (in pure S region) • Hydrostatic consistency: • pressure Jacobian with higher-order integration • Z-method (preferred)

10. Characteristic line t+Dt t Semi-implicit scheme Continuity • Implicit treatment of divergence and pressure gradient terms by-passes most severe CFL condition • ELM: takes advantage of both Lagrangian and Eulerian methods • Grid is fixed in time, and time step is not limited by Courant number condition • Advections are evaluated by following a particle that starts at certain point at time t and ends right at a pre-given point at time t+Dt. • The process of finding the starting point of the path (foot of characteristic line) is called backtracking, which is done by integrating dx/dt=u3 backward in time. • To better capture the particle movement, the backward integration is often carried out in small sub-time steps (Dt /N). • Simple backward Euler method as an option • 5th-order embedded R-K method as an alternative • Interpolation-ELM does not conserve mass; integration ELM does • Numerical diffusion vs. dispersion Momentum b.c. dispersion diffusion • Other considerations • Implicitness factor 0.5 ≤q≤1 to ensure stability • Explicit treatment: Coriolis, baroclinicity, horizontal viscosity, part of bottom velocity

11. ub db Finite-element formulation (1) • A Galerkin weighted residual statement for the continuity equation: • Vertical integration of the momentum equation: • Momentum equation applied to the bottom boundary layer: via z=-h

12. Finite-element formulation (cont’d) • Vertically integrated velocity becomes: • Finally, substituting this eq. back to continuity eq. we get one equation for elevations alone: I2 I3 I1 I4 I5 I6 • Notes: • When node i is on boundaries where essential b.c. is imposed, h is known and so no equations are needed, and so I2 and I6 need not be evaluated there • When node i is on boundaries where natural b.c. is imposed, the velocity is known and so the last term is known  only the first term is truly unknown!

13. Matrices: I1 Wj i i’ is local index of node I in Wj. reaches its max. when i’=l, and so does I1. so diagonal is dominant if the averaged friction-reduced depth ≥0 (can be relaxed) i’ Mass matrix is positive definite and symmetric! > (i’,2) (i’,1)

14. Matrices: I3 and I5 n j i Gij Flather b.c. (overbar denotes mean) j Similarly

15. Matrices: I4 Has most complex form We have decomposed vector f into two parts: one for averaging and the other for integration Most complex part is the baroclinicity at elements/sides

16. Baroclinicity: pressure Jacobian method • No longer used in the newest version • In S layers In Z layers, the derivative can be evaluated directly The vertical integration is done using higher-order rule (Song and Haidvogel) Advantage: no special treatment is needed near surface and bottom Disadvantage: pressure gradient errors

17. Baroclinicity: Z-coordinate method (preferred) Interpolate density along horizontal planes Advantage: alleviate pressure gradient errors (Fortunato and Baptista 1996) Disadvantage: near surface or bottom Method: compute derivatives along two direction at node i and then convert them back to x,y. 2 ELM transport 4 3 (j,k) • Density is defined at nodes and whole levels • Density gradient calculated at side centers and whole levels • Either constant extrapolation (shallow) or more conservative method (deep) is used near bottom • Cubic spline is used in all interpolation of r (and necessary S,T) • Mean density profile removed (optional) • Solve 2 eqs for the density gradient 1 (j,k) Horizontal boundary or dry interface: no flux b.c.

18. Baroclinicity: upwind/TVD transport 3 1 • Density is defined at prism centers (half levels) • Density gradient calculated at prism centers first (for continuity eq.); the values at face centers are averaged from those at prism centers (for momentum eq.) • Either constant extrapolation (shallow) or more conservative method (deep) is used near bottom to avoid spurious flow • Cubic spline is used in all interpolation of rexcept for the averaging for sides • Mean density profile removed (optional) • 3 eqs for the density gradient vs. 2 unknowns – averaging needed for 3 pairs; e.g. 0 i 2 • degenerate cases • Vertical integration using trapzoidal rule

19. Horizontal viscosity Evaluate it at 3 side centers, and integrate using linear quadrature side j The normal derivative is computed using average between adjacent elements Use linear shape function within each element and chain ruleto calculate the derivatives

20. Nz …….. kb+2 kb+1 db z=-h kb Momentum equation • (u,v) solved at side centers and whole levels • Easy of imposing b.c. • Staggering for stability Galerkin FE in the vertical g Solution at bottom from the b.c. there (u=0 for c≠0) The matrix is tri-diagonal

21. 3 II I 1 2 III i Velocity at nodes • Needed for ELM (interpolation) • Method 1: averaging around ball (most diffusive) • Method 2: use linear shape function (conformal ornon-conformal); optional averaging • Filtering of modes • Needs velocity b.c. • Linear interpolation at foot of char. Line • Conformal • Non-conformal – discontinuity along each side

22. 4 1 0 3 2 Shaprio filter Few ocean models are truly free from spurious numerical modes (myriads of processes) A simple filter to reduce sub-grid scale (unphysical) oscillations while leaving the physical signals intact (a=0.5 optimal) No filtering for boundary sides – need b.c. there especially for incoming flow Pre-processor to check the geometry (ipre=1 with 1 processor) Violations usually occur at boundary (fort.11) boundary internal extreme

23. x ELM with Kriging • Best linear unbiased estimator for a random function f(x) • “Exact” interpolator • Works well on unstructured grid Drift function Fluctuation K is called generalized covariance function Minimizing the variance of the fluctuation we get 2-tier Kriging cloud

24. nk+1 k+1 k Vertical velocity • Vertical velocity using Finite Volume Bottom b.c.

25. Inundation algorithms • Algorithm 1 • Update wet and dry elements, sides, and nodes at the end of each time step based on the newly computed elevations • Algorithm 2 (accurate inundation for wetting and drying with sufficient grid resolution) • Wet  add elements & extrapolate velocity • Dry  remove elements • Iterate until the new interface is found at step n+1 • Extrapolate elevation at final interface (smoothing effects) dry A Gn B wet Extrapolation of surface

26. Char. line Transport: ELM • Galerkin FE applied to both nodes and sides • no horizontal diffusion Apply mass lumping reduces FE to FD l+1 x l Tl l-1 l-2 • ELM: mass conservation not guaranteed • Linear with element splitting in the horizontal; cubic spline in the vertical • Quadratic in 3D (works only in estuary due to dispersion) • Impose bounds from surrounding nodes to alleviate dispersion • Kriging (coming up) • Vertical interpolation: placement of T and k • Upwind bias for quadratic interpolation

27. Transport: upwind (1) k S j S Ti,k k-1 S Finite volume discretization in a prism (i,k): mass conservation m≠n, Dt’ ≠Dt; Ti=Ti,k; uj is outward normal velocity Focus on the advection first From continuity equation S+: outflow faces; S -: inflow faces Upwinding

28. Transport: upwind (2) k S j Drop “m” for brevity S Ti,k k-1 S Max. principle is guaranteed if Courant number condition (3D case) Global time step for all prisms Implicit treatment of vertical advective fluxes to bypass vertical Courant number restriction in shallow depths (Modified Courant number condition) • Mass conservation • Budget – vertical and horizontal sums • movement of free surface: violation of volume conservation

29. Solar radiation total solar radiation Budget The source term in the upwind scheme (F.D.) to ensure total heat is conserved (black body) This can cause overheating in shallow depths: adjust albedo

30. Transport: TVD (1) Start from discretized advection upwind Now include downwind component Flux limiter function (depends only on face j); Downwind and central schemes; The key is to find an appropriate function that does not violate max principle upwind Sum of coefficients = 1; so max principle is satisfied if and only if the coefficient of Tiis non-negative Put the last term in similar form as the upwind term d

31. Transport: TVD (2) Max principle is satisfied if (Courant number condition) m Upwind ratio r S+ S- i j Since The Courant number condition is as compared to upwind

32. Transport: TVD (3) • TVD scheme • Mass conservative with max principle; second-order accuracy • Fully explicit scheme due to nonlinearity – vertical Courant number most severe • Boundaries and wet/dry interface – revert back to upwind • Other higher-order schemes: MUSCL; WENO Choice of flux limiting function (gradient detector) Super Bee (compressive) Minmod (diffusive) Osher Van Leer Second order TVD region (Sweby 1984)

33. Turbulence closure (1) Essential b.c. Natural b.c. • FE formulation • Use natural b.c. first and essential b.c. is used to overwrite the solution at boundary • Neglect advection • The production and buoyancy terms are treated explicitly or implicitly depending on the sum to enhance stability • k and y are defined at nodes and whole levels, and diffusivities at half levels and

34. Turbulence closure (2)

35. GOTM turbulence library • General Ocean Turbulence Model • One-dimensional water column model for marine and limnological applications • Coupled to a choice of traditional as well as state-of-the-art parameterizations for vertical turbulent mixing (including KPP) • Some perplexing problems on some platforms – robustness? • Bugs page • Division by 0? www.gotm.net

36. Numerical stability • Semi-implicitness circumvents CFL (most severe) (Casulli and Cattani 1994) • ELM bypasses Courant number condition for advection • Implicit transport scheme along the vertical bypasses Courant number condition • Explicit terms • Baroclinicity  internal Courant number restriction • Horizontal viscosity/diffusivity  diffusion number condition • Upwind and TVD transport – Courant number condition • Coriolis – stable but prone to “modes” 3 1 2

37. Lateral b.c. and nudging

38. Summary