SCHISM numerical formulation Joseph Zhang Horizontal grid: hybrid - PowerPoint PPT Presentation

1 SCHISM numerical formulation Joseph Zhang

Horizontal grid: hybrid 2 3 Side 2 4 3 Side 1 Side 2 Side 1 Side 3 Side 4 1 2 2 1 Side 3 3 2 1 Counter-clockwise convention 1 2 3 3 y s 2 2 x s 1 2 i i 1 N i 1 2 N i i 1 1 N i N i Internal ‘ball’: starting elem is arbitrary • Boundary ball: starting elem is unique •

Vertical grid (1): SZ 3 Terrain-following zone z SZ zone s zone S zone s =0 N z h c h s h s S -levels k z s =-1 k z - 1 Z -levels 2 1     s  s  - s -  s  (1 ) ( ) ( ) ( 1 0) z h h h C  c c    s  -    s tanh ( 1/ 2) tanh( / 2) sinh( )   s  -          f f f  ( ) (1 ) (0 1; 0 20) C   b b b f sinh 2tanh( / 2)  f f  min( , ) h h h s

S -coordinates 4  f =10 -3 , any  b  f =10,  b =1 h c  f =10,  b =0  f =10,  b =0.7

Vertical grid: SZ 5 z = - h s The bottom cells are shaved – no staircases!

Vertical grid (2): LSC 2 via VQS Dukhovskoy et al. (2009) Procedure for VQS Different # of levels are used at different depths • At a given depth, the # of levels is interpolated from the 2 adjacent • master grids Thin layers near the bottom are masked • Staircases appear near the bottom (like Z ), but otherwise terrain • following

What to do with the staircases? A u u P 1 2 B Q R 1 2 LSC 2 = Localized Sigma Coordinates with Shaved Cells ( ivcor =1 in vgrid.in) *No masking of thin layers (c/o implicit scheme)

Vertical grid: LSC 2 8 z (m) Degenerate prism

Vertical grid (3) 9 • 3D computational unit: uneven prisms • Equations are not transformed into S - or s -coordinates, but solved in the original Z -space  • Pressure gradient     d z • Z -method with cubic spline (extrapolation on surface/bottom) n k n k k N z k (whole level)  N z - 1 u j,k+1 v j,k+1  ……. u j,k+1 v j,k+1 k C i,k C i,k k- 1 ……. k -1 k -1 w i,k -1 k b + 1 w i,k- 1 k b

Polymorphism Zhang et al. (2016) Bathymetry (m)

Polymorphism z (m) flood x (m) S (PSU) *Make sure there is no jump in Cd between 2D/3D zones

Semi-implicit scheme 12   -  1 n n Continuity          -    1 n n (1 ) 0 u dz u dz  - - t h h 𝒗 𝒐+𝟐 − 𝒗 ∗ = 𝒈 − 𝑕𝜄𝛼𝜃 𝑜+1 − 𝑕(1 − 𝜄)𝛼𝜃 𝑜 + 𝒏 𝒜 𝒐+𝟐 − 𝛽 𝒗 𝒗 𝒐+𝟐 𝑀 𝑦, 𝑧, 𝑨 Momentum Δ𝑢    n 1 u      τ 1 n n n  , at ; z b.c. (3D)   w z    n 1  u      -   1 n n n n n , at , | | u z h C u    b D b z  - 1 n u u u D    * Advection: ; ( , , , , ) u x y z t t t ≈  * Dt t Implicit treatment of divergence, pressure gradient and SAV form drag terms by-passes  most severe stability constraints: 0.5 ≤  ≤ 1 Explicit treatment: Coriolis, baroclinicity , horizontal viscosity, radiation stress… 

Eulerian-Lagrangian Method (ELM) 13 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 +  t . The process of finding the starting point of the path (foot of characteristic line) is called  backtracking , which is done by integrating d x / dt = u 3 backward in 3D space and time. To better capture the particle movement, the backward integration is often carried out in small  sub-time steps  Simple backward Euler method  2nd-order R-K method  Interpolation-ELM does not conserve; integration ELM does Interpolation: numerical diffusion vs. dispersion    u t 1  x     Characteristic line t +  t advection t Flow dispersion diffusion

Operating range of SCHISM 14 …. for a fixed grid 𝐷𝐺𝑀 = ( 𝒗 + 𝑕ℎ)∆𝑢 ∆𝑦 ELM prefers larger  t (truncation error ~1/  t) • As a result, SCHISM requires CFL>0.4 • Convergence: CFL=const,  t  0 • For barotropic field applications: [100,450] sec • For baroclinic field applications: [100,200] sec (e.g., start with 150 sec) •

ELM with high-order Kriging 15 • Best linear unbiased estimator for a random function f ( x ) • “Exact” interpolator • Works well on unstructured grid (efficient) • Needs filter (ELAD) to reduce dispersion N          - ( , ) (| |) f x y x y K x x 1 2 3 i i  1 i Cubic spline Drift function Fluctuation K is called generalized covariance function  - 2 3 ( ) , ln , , K h h h h h Minimizing the variance of the fluctuation we get N    0 i  1 i 2-tier Kriging cloud N    0 x i i  1 i N    0 y i i  1 i x  ( , ) f x y f i i i • The numerical dispersion generated by kriging is reduced by a diffusion ELAD filter (Zhang et al. 2016)

Semi-implicit scheme 16   -  1 n n Continuity          -    1 n n (1 ) 0 u dz u dz  - - t h h 𝒗 𝒐+𝟐 − 𝒗 ∗ = 𝒈 − 𝑕𝜄𝛼𝜃 𝑜+1 − 𝑕(1 − 𝜄)𝛼𝜃 𝑜 + 𝒏 𝒜 𝒐+𝟐 − 𝛽 𝒗 𝒗 𝒐+𝟐 𝑀 𝑦, 𝑧, 𝑨 Momentum Δ𝑢    n 1 u      τ 1 n n n  , at ; z b.c.   w z    n 1  u      -   1 n n n n n , at , | | u z h C u    b D b z  - 1 n u u u D    * ; ( , , , , ) u x y z t t t ≈  * Dt t Implicit treatment of divergence and pressure gradient terms by-passes most severe CFL  condition : 0.5 ≤  ≤ 1 Explicit treatment: Coriolis, baroclinicity , horizontal viscosity, radiation stress… 

Galerkin finite-element formulation 17 A Galerkin weighted residual statement for the continuity equation: 𝜃 𝑜+1 − 𝜃 𝑜 𝛼𝜚 𝑗 ∙ 𝑽 𝑜+1 𝑒Ω + න 𝑜+1 𝑒Γ + 1 − 𝜄 𝛼𝜚 𝑗 ∙ 𝑽 𝑜 𝑒Ω + න 𝑜 𝑒Γ = 0, න 𝜚 𝑗 𝑒Ω + 𝜄 − න 𝜚 𝑗 𝑉 𝑜 − න 𝜚 𝑗 𝑉 𝑜 Δ𝑢 Ω Ω Γ Ω Γ (𝑗 = 1, ⋯ , 𝑂 𝑞 ) 𝜚 𝑗 : shape/weighting function; U : depth-integrated velocity;  :implicitness factor Need to eliminate U n+1 to get an equation for  alone. We’ll do this with the aid from momentum equation: 𝒗 𝒐+𝟐 − 𝒗 ∗ = 𝒈 − 𝑕𝜄𝛼𝜃 𝑜+1 − 𝑕(1 − 𝜄)𝛼𝜃 𝑜 + 𝒏 𝒜 𝒐+𝟐 − 𝛽 𝒗 𝒗 𝒐+𝟐 𝑀 𝑦, 𝑧, 𝑨 Δ𝑢 Shape function (global)

2D case 18 The unknown velocity can be easily found as: 𝑯 − 𝑕𝜄Δ𝑢 𝐼 2 𝑽 𝑜+1 = ෙ 𝐼 𝛼𝜃 𝑜+1 ෩ 𝑯 = 𝐼 𝐼 𝑽 ∗ + ∆𝑢 𝑮 + 𝝊 𝑥 − 𝑕(1 − 𝜄)𝐼𝛼𝜃 𝑜 (explicit terms) ෙ ෩ ෩ 𝐼 = 𝐼 + (𝜓 + 𝛽|𝒗|𝐼)Δ𝑢

3D case 19 Vertical integration of the momentum equation u b 𝑽 𝒐+𝟐 − 𝑽 ∗ 𝑜+1 − 𝛽 𝒗 𝑽 𝜷 = 𝑮 − 𝑕𝜄𝛼𝜃 𝑜+1 − 𝑕 1 − 𝜄 𝛼𝜃 𝑜 + 𝝊 𝑥 − 𝜓𝒗 𝑐 d b Δ𝑢 z = - h 𝑨 𝑤 𝒜 𝒘 Assuming: 𝒗𝑒𝑨 ≡ 𝒗 𝑽 𝜷 න |𝒗|𝒗𝑒𝑨 ≅ 𝒗 න −ℎ −𝒊 Momentum equation applied to the bottom boundary layer 𝒐+𝟐 − 𝒗 𝒄 ∗ 𝒗 𝒄 𝒐+𝟐 + 𝜖 𝜖𝑨 𝜉 𝜖𝒗 = 𝒈 𝒄 − 𝑕𝜄𝛼𝜃 𝑜+1 − 𝑕 1 − 𝜄 𝛼𝜃 𝑜 − 𝛽 𝒗 𝒄 𝒗 𝒄 Δ𝑢 𝜖𝑨 Therefore 1 𝑕𝜄∆𝑢 𝑜+1 = ∗ + 𝒈 𝒄 ∆𝑢 − 𝑕(1 − 𝜄)∆𝑢𝛼𝜃 𝑜 − 1 + 𝛽|𝒗 𝑐 |∆𝑢 𝛼𝜃 𝑜+1 𝒗 𝑐 1 + 𝛽|𝒗 𝑐 |∆𝑢 𝒗 𝑐 Now still need to find U 

3D case: locally emergent vegetation 20 We have 𝑽 𝛽 = 𝑽 n+1 So: 𝑕𝜄 ෡ 𝐼Δ𝑢 𝑽 𝑜+1 = 𝑯 𝟐 − 𝛼𝜃 𝑜+1 1 + 𝛽|𝒗|∆𝑢 ∗ + 𝒈 𝑐 ∆𝑢) 𝑯 1 = 𝑽 ∗ + 𝑮 + 𝝊 𝒙 ∆𝑢 − 𝑕 1 − 𝜄 ෡ 𝐼∆𝑢𝛼𝜃 𝑜 − ෤ 𝜓∆𝑢(𝒗 𝑐 1 + 𝛽|𝒗|∆𝑢 ෡ 𝐼 = 𝐼 − ෤ 𝜓∆𝑢 𝜓 𝜓 = ෤ 1 + 𝛽|𝒗 𝒄 |∆𝑢