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

SLIDE 1

1

SCHISM numerical formulation

Joseph Zhang

SLIDE 2

1 2 Ni 1

Ni

i

Horizontal grid: hybrid

3 1 2 3 Side 1 Side 2 Side 3 1 2 ys xs i 1 2 3 2 1 3 1 2 3 4 Side 1 Side 2 Side 3 Side 4 i 1 1 2 Ni 2 Ni 2 Counter-clockwise convention

Internal ‘ball’: starting elem is arbitrary
Boundary ball: starting elem is unique

SLIDE 3

3 s zone

Vertical grid (1): SZ

1 hs 2 Nz kz s=0 s=-1 S-levels Z-levels kz-1

z

(1 ) ( ) ( ) ( 1 0) tanh ( 1/ 2) tanh( / 2) sinh( ) ( ) (1 ) (0 1; 0 20) sinh 2tanh( / 2) min( , )

c c f f f b b b f f f s

z h h h C C h h h  s s s s  s   s s           



    



  



      

S zone SZ zone Terrain-following zone

hs hc

SLIDE 4

4

S-coordinates

f=10-3, any b f=10, b=1 f=10, b=0 f=10, b=0.7

hc

SLIDE 5

5

Vertical grid: SZ

z=-hs

The bottom cells are shaved – no staircases!

SLIDE 6

Vertical grid (2): LSC2 via VQS

Dukhovskoy et al. (2009)

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 Procedure for VQS

SLIDE 7

What to do with the staircases?

1 2 P Q R A B u u 1 2

LSC2= Localized Sigma Coordinates with Shaved Cells

(ivcor =1 in vgrid.in) *No masking of thin layers (c/o implicit scheme)

SLIDE 8

Vertical grid: LSC2

8

z (m)

Degenerate prism

SLIDE 9

9

Vertical grid (3)

3D computational unit: uneven prisms
Equations are not transformed into S- or s-coordinates, but solved in the original Z-space
Pressure gradient
Z-method with cubic spline (extrapolation on surface/bottom)

z

d



  



……. Nz Nz-1 kb+1 kb

nk

k (whole level) k-1

uj,k+1 vj,k+1 Ci,k  k-1 k …….

nk

k k-1

uj,k+1 vj,k+1 Ci,k  wi,k-1 wi,k-1

SLIDE 10

Bathymetry (m)

Polymorphism

Zhang et al. (2016)

SLIDE 11

Polymorphism

flood S (PSU)

x (m) z (m)

*Make sure there is no jump in Cd between 2D/3D zones

SLIDE 12

Semi-implicit scheme

1 1

(1 )

n n n n h h

dz dz t

 

   

 

 





 

 

u u

1 1 1 1

, at ; , at ,

n n n n w n n n n b

z z z h z    

   

            -    u τ u u

Continuity Momentum b.c. (3D)

| |

n n D b

C   u



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…

1 * *

; ( , , , , )

n

D x y z t t t Dt t

 -

    u u u u

≈ 12 Advection: 𝒗𝒐+𝟐 − 𝒗∗ Δ𝑢 = 𝒈 − 𝑕𝜄𝛼𝜃𝑜+1 − 𝑕(1 − 𝜄)𝛼𝜃𝑜 + 𝒏𝒜

𝒐+𝟐 − 𝛽 𝒗 𝒗𝒐+𝟐𝑀 𝑦, 𝑧, 𝑨

SLIDE 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 dx/dt=u3 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

Eulerian-Lagrangian Method (ELM) t+t t

   

Characteristic line

1 u t x   

diffusion dispersion advection

13 Flow

SLIDE 14

Operating range of SCHISM

…. 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)

14

𝐷𝐺𝑀 = ( 𝒗 + 𝑕ℎ)∆𝑢 ∆𝑦

SLIDE 15

15

ELM with high-order Kriging

Best linear unbiased estimator for a random function f(x)
“Exact” interpolator
Works well on unstructured grid (efficient)
Needs filter (ELAD) to reduce dispersion

x

2-tier Kriging cloud

1 2 3 1

( , ) (| |)

N i i i

f x y x y K    



   



x x 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

1 1 1

( , )

N i i N i i i N i i i i i i

x y f x y f   

  

   

  

Cubic spline

The numerical dispersion generated by kriging is reduced by a diffusion ELAD

filter (Zhang et al. 2016)

SLIDE 16

Semi-implicit scheme

1 1

(1 )

n n n n h h

dz dz t

 

   

 

 





 

 

u u

1 1 1 1

, at ; , at ,

n n n n w n n n n b

z z z h z    

   

            -    u τ u u

Continuity Momentum b.c.

| |

n n D b

C   u



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…

1 * *

; ( , , , , )

n

D x y z t t t Dt t

 -

    u u u u

≈ 16

𝒗𝒐+𝟐 − 𝒗∗ Δ𝑢 = 𝒈 − 𝑕𝜄𝛼𝜃𝑜+1 − 𝑕(1 − 𝜄)𝛼𝜃𝑜 + 𝒏𝒜

𝒐+𝟐 − 𝛽 𝒗 𝒗𝒐+𝟐𝑀 𝑦, 𝑧, 𝑨

SLIDE 17

17

Galerkin finite-element formulation

A Galerkin weighted residual statement for the continuity equation: Shape function (global)

න

Ω

𝜚𝑗 𝜃𝑜+1 − 𝜃𝑜 Δ𝑢 𝑒Ω + 𝜄 − න

Ω

𝛼𝜚𝑗 ∙ 𝑽𝑜+1𝑒Ω + න

Γ

𝜚𝑗𝑉𝑜

𝑜+1𝑒Γ + 1 − 𝜄

− න

Ω

𝛼𝜚𝑗 ∙ 𝑽𝑜𝑒Ω + න

Γ

𝜚𝑗𝑉𝑜

𝑜𝑒Γ = 0,

(𝑗 = 1, ⋯ , 𝑂𝑞) 𝜚𝑗: shape/weighting function; U: depth-integrated velocity; :implicitness factor

Need to eliminate Un+1 to get an equation for  alone. We’ll do this with the aid from momentum equation: 𝒗𝒐+𝟐 − 𝒗∗ Δ𝑢 = 𝒈 − 𝑕𝜄𝛼𝜃𝑜+1 − 𝑕(1 − 𝜄)𝛼𝜃𝑜 + 𝒏𝒜

𝒐+𝟐 − 𝛽 𝒗 𝒗𝒐+𝟐𝑀 𝑦, 𝑧, 𝑨

SLIDE 18

18

2D case

The unknown velocity can be easily found as:

𝑽𝑜+1 = ෙ 𝑯 − 𝑕𝜄Δ𝑢 𝐼2 ෩ 𝐼 𝛼𝜃𝑜+1

ෙ 𝑯 = 𝐼 ෩ 𝐼 𝑽∗ + ∆𝑢 𝑮 + 𝝊𝑥 − 𝑕(1 − 𝜄)𝐼𝛼𝜃𝑜

෩ 𝐼 = 𝐼 + (𝜓 + 𝛽|𝒗|𝐼)Δ𝑢

(explicit terms)

SLIDE 19

19

3D case

db

ub

z=-h

Vertical integration of the momentum equation Momentum equation applied to the bottom boundary layer 𝑽𝒐+𝟐 − 𝑽∗ Δ𝑢 = 𝑮 − 𝑕𝜄𝛼𝜃𝑜+1 − 𝑕 1 − 𝜄 𝛼𝜃𝑜 + 𝝊𝑥 − 𝜓𝒗𝑐

𝑜+1 − 𝛽 𝒗 𝑽𝜷

න

−ℎ 𝑨𝑤

|𝒗|𝒗𝑒𝑨 ≅ 𝒗 න

−𝒊 𝒜𝒘

𝒗𝑒𝑨 ≡ 𝒗 𝑽𝜷 Assuming: 𝒗𝒄

𝒐+𝟐 − 𝒗𝒄 ∗

Δ𝑢 = 𝒈𝒄 − 𝑕𝜄𝛼𝜃𝑜+1 − 𝑕 1 − 𝜄 𝛼𝜃𝑜 − 𝛽 𝒗𝒄 𝒗𝒄

𝒐+𝟐 + 𝜖

𝜖𝑨 𝜉 𝜖𝒗 𝜖𝑨 𝒗𝑐

𝑜+1 =

1 1 + 𝛽|𝒗𝑐|∆𝑢 𝒗𝑐

∗ + 𝒈𝒄∆𝑢 − 𝑕(1 − 𝜄)∆𝑢𝛼𝜃𝑜 −

𝑕𝜄∆𝑢 1 + 𝛽|𝒗𝑐|∆𝑢 𝛼𝜃𝑜+1 Therefore Now still need to find U

SLIDE 20

20

3D case: locally emergent vegetation

𝑽𝑜+1 = 𝑯𝟐 − 𝑕𝜄 ෡ 𝐼Δ𝑢 1 + 𝛽|𝒗|∆𝑢 𝛼𝜃𝑜+1 𝑽𝛽 = 𝑽n+1

We have So: 𝑯1 = 𝑽∗ + 𝑮 + 𝝊𝒙 ∆𝑢 − 𝑕 1 − 𝜄 ෡ 𝐼∆𝑢𝛼𝜃𝑜 − ෤ 𝜓∆𝑢(𝒗𝑐

∗ + 𝒈𝑐∆𝑢)

1 + 𝛽|𝒗|∆𝑢 ෤ 𝜓 = 𝜓 1 + 𝛽|𝒗𝒄|∆𝑢 ෡ 𝐼 = 𝐼 − ෤ 𝜓∆𝑢

SLIDE 21

21

3D case: locally submerged vegetation

Integrating from bottom to top of canopy: 𝑽𝛽 = 𝑽∗𝛽 + 𝑮𝜷∆𝑢 − 𝑕𝜄𝐼𝛽∆𝑢𝛼𝜃𝑜+1 − 𝑕(1 − 𝜄)𝐼𝛽∆𝑢𝛼𝜃𝑜 − 𝛽∆𝑢 𝒗 𝑽𝛽 + ∆𝑢𝜉 ቤ 𝜖𝒗 𝜖𝑨 −ℎ

𝑨𝑤

𝑽∗𝛽 = න

−ℎ 𝑨𝑤

𝒗∗𝑒𝑨, 𝑮𝜷 = න

−ℎ 𝑨𝑤

𝒈𝑒𝑨 𝜉 𝜖𝒗 𝜖𝑨 ≡ 𝑣′𝑥′ = 𝑺0𝑓𝛾2(𝑨−𝑨𝑤), (𝑨 ≤ 𝑨𝑤) Inside vegetation layer, the stress obeys exponential law (Shimizu and Tsujimoto 1994) 𝛾2 = 𝑂𝑤 𝐼𝛽 −0.32 − 0.85 log10 𝐼 − 𝐼𝛽 𝐼𝛽 𝐽 𝛾2 is a constant found from empirical formula: 𝐽 = 𝜓|𝒗𝑐| 𝑕𝐼 Therefore 𝜉 ቤ 𝜖𝒗 𝜖𝑨 −ℎ

𝑨𝑤

= 𝛾𝜓𝒗𝑐

𝑜+1

(𝛾 = 𝑓𝛾2 𝑨𝑤−𝑨𝑐 − 1) 𝑽𝛽 = 𝑯3 − 𝑕𝜄 ෡ 𝐼𝛽∆𝑢 1 + 𝛽|𝒗|∆𝑢 𝛼𝜃𝑜+1 Finally: 𝑯3 = 𝑽∗𝛽 + 𝑮𝜷∆𝑢 + 𝛾 ෤ 𝜓∆𝑢 𝒗𝑐

∗ + 𝒈𝑐∆𝑢 − 𝑕(1 − 𝜄) ෡

𝐼𝛽∆𝑢𝛼𝜃𝑜 1 + 𝛽|𝒗|∆𝑢 ෡ 𝐼𝛽 = 𝐼𝛽 + 𝛾 ෤ 𝜓∆𝑢 𝑽𝑜+1 = 𝑯𝟑 − 𝑕𝜄 ന 𝐼Δ𝑢𝛼𝜃𝑜+1 ന 𝐼 = 𝐼 − ෤ 𝜓Δ𝑢 − 𝑑 ෡ 𝐼𝛽 𝑑 = 𝛽|𝒗|∆𝑢 1 + 𝛽|𝒗|∆𝑢

SLIDE 22

22

General case

(Note the difference in ෙ 𝐼 between 2D/3D) 𝑽𝑜+1 = 𝑭 − 𝑕𝜄 ෱ 𝐼Δ𝑢𝛼𝜃𝑜+1 ෱ 𝐼 = 𝐼2 ෩ 𝐼 , locally 2D ෡ 𝐼 1 + 𝛽|𝒗|∆𝑢 , locally 3D emergent ന 𝐼, locally 3D submerged 𝑭 = ൞ ෱ 𝑯, 2D 𝑯𝟐, 3D emergent 𝑯𝟑, 3D submerged

In compact form:

SLIDE 23

23

Finite-element formulation (cont’d)

Finally, substituting this eq. back to continuity eq. we get one equation for elevations alone: Notes:

When node i is on boundaries where essential b.c. is imposed,  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
n LHS is known  only the first term is truly unknown!
b.c. is free lunch in FE model (cf. staggering in FV)
The inherent numerical dispersion in FE nicely balances out the numerical diffusion inherent in the implicit

time stepping scheme!

I3 I1 I4 I5 I6 sources (i=1,…,Np) න

Ω

𝜚𝑗𝜃𝑜+1 + 𝑕𝜄2Δ𝑢2 ෙ 𝐼𝛼𝜚𝑗 ∙ 𝛼𝜃𝑜+1 𝑒Ω = න

Ω

𝜚𝑗𝜃𝑜 + 𝜄Δ𝑢𝛼𝜚𝑗 ∙ 𝑭 + (1 − 𝜄)Δ𝑢𝛼𝜚𝑗 ∙ 𝐕𝑜 𝑒Ω − 𝜄Δ𝑢 න

Γ𝑤

𝜚𝑗 ෡ 𝑉𝑜

𝑜+1𝑒Γ𝑤 − (1 − 𝜄) Δ𝑢 න Γ

𝜚𝑗𝑉𝑜

𝑜𝑒Γ − 𝜄Δ𝑢 න ഥ Γ𝑤

𝜚𝑗𝑉𝑜

𝑜+1𝑒ത

Γ𝑤 + 𝐽7

SLIDE 24

24

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

j j+1 j-1

x

1 N j j+1 j-1

x

1 N

1

, ( )

N i i i j ij i

u u x   d



 



Conformal Non-conformal

1 1 1

0, 1,...,

m

N M N i i j i i j m i m i

L u d L u d

j N

   

   

              

 

    

(Elevation, velocity) (optional for velocity: indvel=0 -> less dissipation and needs further stabilization (MB* schemes)

SLIDE 25

25

Matrices: I1

1 2 2 1 1 ' ' 1 2 3 ', 1 , 1 1

ˆ ˆ ˆ ( ) 1 ˆ ' 12 4

i j i

N n n i i j j N i l n j l j j l j

I g t H d g t A Hi l A      d  

      

                  

  

reaches its max when i’=l, and so does I1.

' i l 

i’ is local index of node i in j.

3 ' 1

'

i





so diagonal is dominant if the averaged friction-reduced depth ≥0 (can be relaxed) Mass matrix is positive definite and symmetric! JCG or PETSc solvers work fine ' i i j i’ (i’,1) (i’,2) >

i34(j)

ෙ 𝐼 ෙ 𝐼

SLIDE 26

26

Matrices: I3 and I5

1 1 3 1 1 1 , 1 n , 1 n ,

ˆ ˆ ˆ ˆ ( ) ( ) 4

v ij z bs

n n i n v i n ij j N ij n n ij k ij k ij k j k k

I U d U d L z u u  

   



   

          

    

Flather b.c. (overbar denotes mean)

1 1 n n n 1 1 3

ˆ / ( ) ( ) 2( ) ( ) 2 6

n n ij ij n n ij i i j j j

u u g h gh U I L      

   




      




      



1 5 ' , 1 n , 1 n ,

ˆ ˆ ˆ ˆ ( ) ( ) 4

z bs ij

N ij n n n i n ij ij k ij k ij k j j k k

L I U d z u u 



  

        

   

Similarly i j ij j n

SLIDE 27

27

Matrices: I7

Related to source/sink Bi: ball of node i Um: elements that have source/sink Sj: volume discharge rate [m3/s]

SLIDE 28

28

Matrices: I4

3 ' 3 , ', ' 1

ˆ ˆ ˆ ˆ (1 ) (1 ) (1 ) 12

i

N j n n n n i i i j l i l j i j i j j l

A I t t d tA t       d    

 

     

 

        



           

  

U G U G

Has the most complex form Most complex part of E is the baroclinicity

c z

g d



    - 



f

at elements/sides

i34(j)

E

SLIDE 29

29

Baroclinicity: z-method

Interpolate density along horizontal planes
Advantage: alleviate pressure gradient errors (Fortunato and Baptista 1996)
Disadvantage: near surface or bottom
Density is defined at prism centers (half levels)
Re-construction method: compute directional derivatives along two direction at node i and

then convert them back to x,y

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.)

constant extrapolation (shallow) is used near bottom (under resolution)
Cubic spline is used in all interpolations of 
Mean density profile can be optionally removed
3 or 4 eqs for the density gradient vs. 2 unknowns – averaging needed for 3 pairs; e.g.
Degenerate cases: when the two vectors are co-linear (discard)
Vertical integration using trapzoidal rule

( ) z 

' ' ' ' , 1 , 1 1 ' ' ' ' , 2 , 2 2

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

x i k y i k x i k y i k

x x y y x x y y        





Element j

3 1 2 A B

SLIDE 30

30

Momentum equation: 2D case

𝒗𝒐+𝟐 = ෱ 𝐼 𝐼 𝒗∗ + (𝒈 + 𝝊𝒙/𝐼)∆𝑢 − 𝑕𝜄𝛼𝜃𝑜+1 − 𝑕(1 − 𝜄)𝛼𝜃𝑜 The depth-averaged velocity can be directly solved once unknown elevations are found:

SLIDE 31

31

Momentum equation: 3D case

Galerkin FE in the vertical …….. kb kb+1 Nz kb+2

z=-h db

1 1 1 1 1 1 1 1 1 1 1 1 1 1/ 2 1 1/ 2 , 1 1 1 1 1 1 , 1

( ) (2 ) ( 1 ) (2 ) 6 6 ( ) (2 ) ( 1 ) (2 6 6

b b z

n n n n n n n n n l l l l l l z l l l b l l l l k k l l n n l l l N w z l l b l

z z u u N l t k l t t z z z z t N l k l   d  d

          







     

    





 



                 



 

u

u u u u u u τ g g g A A + A A

1 1 )

( 1, )

n l b z

l k N





  g

The matrix is tri-diagonal (Neumann b.c. can be imposed easily)
Solution at bottom from the b.c. there (u=0 for ≠0); otherwise free slip
Vegetation terms are also treated implicitly

න

𝑨𝑐 𝜃

𝜒𝑚 𝑐𝒗𝑜+1 − Δ𝑢 𝜖 𝜖𝑨 𝜉 𝜖𝒗𝑜+1 𝜖𝑨 𝑒𝑨 = න

𝑨𝑐 𝜃

𝜒𝑚𝒉𝑒𝑨 , (𝑚 = 𝑙𝑐 + 1, ⋯ , 𝑂𝑨)

𝜒𝑚: vertical hat function; g: explicit terms 𝑐 = 1 + 𝛽|𝒗𝑜|Δ𝑢ℋ(𝑨𝑤 − 𝑨)

b b

SLIDE 32

32

Velocity at nodes

Needed for ELM (interpolation), and can serve as a way to further stabilize the

mom solver by adding some inherent numerical diffusion

Method 1: inverse-distance averaging around ball (indvel=1 or MA schemes) –

larger dissipation; no further stabilization needed

Method 2: use linear shape function (conformal) (indvel=0 or MB schemes)
Discontinuous across elements; averaging to get the final value at node
Generally needs velocity b.c. (for incoming info)
indvel=0 generally gives less dissipation but needs further stabilization in

the form of viscosity or Shapiro filter for stabilization (important for eddying regime)

1 II III I

 

u

u u u

i

1 2 3 I II III

Dispersion Diffusion Eternal struggle between diffusion and dispersion

SLIDE 33

33

Shaprio filter

A simple filter to reduce sub-grid scale (unphysical) oscillations while leaving the physical signals largely intact (=0.5

ptimal)

4 1

ˆ 4 , (0 0.5) 4

i i

 



   



    



u u u u

1 2 3 4

No filtering for boundary sides – so need b.c. there especially for incoming flow
Equivalent to Laplacian viscosity

SLIDE 34

34

Horizontal viscosity

Assuming uniformity

       





d n u A A u

II I

) (  

x 1 2 4 3 x0 1 2 3 4 5 6 I II I II P Q R S

) 4 ( 3 ) (

4 3 2 1

u u u u u A u

I



       

(like Shapiro filter)

 

4 4 3 3 2 2 1 1 4 3 2 1 2 2 3 2 3 2 2 2 1 2 3 4

20 ) ( 7 ) 4 ( u u u u u u u u u u u u u t u u u u u u

b a b a b a b a



    



     







 

Bi-harmonic viscosity (ideal in eddying regime) x 1 2 4 3 1a 1b 2a 2b 3a 3b 4a 4b Important for eddying regime for adding proper dissipation (Zhang et al. 2016)

SLIDE 35

35

Vertical velocity



Vertical velocity using Finite Volume

nk+1 k+1 k

ˆ S

ˆ P

1 1 1 1 1 1 1 1 1 1 1 , 1 1 , 3 1 1 ( , ) , ( , ), ( , ), 1 1

ˆ ˆ ( ) ( ) ˆ ˆ ˆ ( )/ 2 0, ( ,..., 1)

n x n y n z n x n y n z k k k k k i k k k k k k k i k k n n js i m i m js i m k js i m k b z m

S u n v n w n S u n v n w n P s q q k k N

                

 



    



3 1 1 1 1 1 1 1 1 , 1 ( , ) , ( , ), ( , ), 1 1 1 1 1 1 , 1 1 1

1 ˆ ˆ ˆ ˆ ˆ ( )/ 2 ( ) ( ) , ˆ ( ,..., 1)

n n n n x n y n x n y n z i k js i m i m js i m k js i m k k k k k k k k k k k i k k z m k k b z

w P s q q S u n v n S u n v n w n n S k k N

                 

  





     



, ,

ˆ j k

j k j

q   u n

Bottom b.c.

1 1 1 ,

, ( )

n x n y n z k k k k i k k b

b u n v n w n k k t

  

     

Surface b.c. not enforced (closure error)

i34(j)

෍

𝑘∈𝑇+

|𝑅𝑘| = ෍

𝑘∈𝑇−

|𝑅𝑘|

In compact form: This is the foundation of mass conservation and monotonicity in transport solvers

SLIDE 36

36

Inundation algorithms



Algorithm 1 (inunfl=0)



Update wet and dry elements, sides, and nodes at the end of each time step based on the newly computed elevations. Newly wetted vel., S,T etc calculated using average



Algorithm 2 (inunfl=1; 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)

Extrapolation of surface A wet n B dry

SLIDE 37

37

Transport: upwind (1)

1 1 1 1 1 5 , 1 , , , 1 * , , , 1 , 1/ 2 , 1/ 2 3 1

' ' ' ( ) , ( 1, , )

m m m m m m i k i k i k i k m i i i j j j i i k i i k i k j V i k i k m m j i h j j b z j ij

T T T T T T V u S T V Q A t t z z T T t S k k N    d

     







 



  



      



 

 

k k-1

Ti,k

3

ˆ ( ) ( )

h

T T T Q T t z z                       u

ˆ, 0, T T z z T z h n           - 

S Finite volume discretization in a prism (i,k): mass conservation Focus on the advection first

m≠n, t’ ≠ t; Ti=Ti,k; uj is outward normal velocity

S

1 *

'

m m i i j j n j V i

t T T Q T V

 

 



j j j

Q u S 

From continuity equation

| | | |

j j j j j S j S

Q Q Q







  

  

S+: outflow faces; S -: inflow faces

Upwinding

*

, ,

i j j

T j S T T j S





      

j (sources: bdy_frc, flx_sf and flx_bt)

^

i34(j)+2 i34(j)

Ti*

SLIDE 38

38

Transport: upwind (2)

1

' ' ' | | | | 1 | | | | (1 )

m i i j i j j j i j j n j S j S j S j S i i i i j

t t t T T Q T Q T Q T Q T V V V T T  






   

       




        



   

k k-1

Ti,k S S S j Drop “m” for brevity (and keep only advection terms)

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: small conservation error
Includes source/sink

1 −

∆𝑢′ 𝑊𝑗 σ𝑘∈𝑇− |𝑅𝑘| ≥ 0 ֜ ∆𝑢′ ≤ 𝑊𝑗 σ𝑘∈𝑇− |𝑅𝑘|

∆𝑢′ ≤ 𝑊

𝑗

σ𝑘∈𝑇−𝐼 |𝑅𝑘|

SLIDE 39

39

Solar radiation

1 2

/ /

1 ˆ , ( ) (0) Re (1 R)e

D d D d p

H Q H z H C z 



    



 

(0) ˆ

p

H Qdz C 







total solar radiation

Budget The source term in the upwind scheme (F.D.)

, , 1 , , 1 , ,

( ) ˆ

b

m m m m i k i k i i k i k m i i i p i k p m i k

H H A H H V V Q C z C H  



  

*Assume the heat is completely absorbed by the bed (black body) to prevent overheating in shallow depths

SLIDE 40

40

Transport: TVD2

) ,..., 1 ( , ˆ ) ( | |

1

M m Q V t C C Q V t C C

H H

S j m j j i m S j m i m j j i m n i m i

 



 

 

  



) ( , ) ( ) ~ ~ ( | | ~

1 z b S j j j j i S j i j j i M i i

,.., N k j Q V t C C Q V t C C

V V

    



 

 

  

)

( , ~

1 1 , 1 , 1 z b V n h i n k i n k i i i i n i

,.., N k k dV F V t z C z C V t A C C

i

                  



         









 

Solve the transport eq in 3 sub-steps: Horizontal advection (explicit TVD method or WENO) Vertical advection (implicit TVD2) Remaining (implicit diffusion)

1 1 1 1 2 / 1

] [ 2 ] [

    

  



       

n jup n jup n jup j j n jup n j

t C t C C C C r

1st step: explicit TVD or WENO3
2nd step: Taylor expansion in space and time

Space limiter Time limiter Duraisamy and Baeder (2007) IMPORTANT: built-in monotonicity constraint is very important for higher-order transport; use of additional diffusion would degrade order of accuracy

SLIDE 41

41

TVD2

) ( | | 2 1 ) ~ ~ ( 2 1 1 | | ~

1     









              



 

   





M i S j j S q q q j i S j j i j S p p p j i i

C s Q V t C C r Q V t C

V V V V

   



   









 



V

M q q q S p p M i i q V p i p S q i q q p

S q C C Q Q C C s S p C C Q C C Q r

V V

, ) ~ ( | | | | ) ~ ( , ) ~ ~ ( | | ) ~ ~ ( | |

1 1

Upwind ratio Downwind ratio

2 1 1         







 j S p p p

V r

 

) ( | | 2 1  





 







d  

V V

S j j S q q q j i

s Q V t

TVD condition: (iterative solver)

Iterative solver converges very fast
2nd-order accurate in space and time
Monotone
Alleviate the sub-cycling/small t for

transport that plagued explicit TVD

Can be used at any depths!

SLIDE 42

42

Turbulence closure (itur=3)

0,

r

, ,

k

k z h z n z h z l n z z l

    

                -     -    -  

2

1 ,

r

( ) ,

r

( ) ( ) ,

r

p m n

k z h c z l z h c k z h

 

          -     -    - u

Essential b.c. Natural b.c.

FE formulation

k,  and diffusivities are defined at nodes and whole levels
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 𝐸𝑙 𝐸𝑢 = 𝜖 𝜖𝑨 𝜉𝑙

𝜔 𝜖𝑙

𝜖𝑨 + 𝜉𝑁2 + 𝜆𝑂2 + 𝑑𝑔𝑙𝛽 𝒗 3ℋ 𝑨𝑤 − 𝑨 − 𝜗 𝐸𝜔 𝐸𝑢 = 𝜖 𝜖𝑨 𝜉𝜔 𝜖𝜔 𝜖𝑨 + 𝜔 𝑙 𝑑𝜔1𝜉𝑁2 + 𝑑𝜔3𝜆𝑂2 + 𝑑𝑔𝜔𝛽 𝒗 3ℋ(𝑨𝑤 − 𝑨) − 𝑑𝜔2𝜗𝐺𝑥𝑏𝑚𝑚

SLIDE 43

43

Turbulence closure (itur=3)

       

1 1 1 1 1 1 1 1 1/ 2 1 1 1 1 1 1 1 1 1/ 2 2 2 1 1/ 2 2 2 1 1/ 1/ 2

( ) 2 2 ( ) 6 ( ) 2 2 ( ) 6 2 6

n n n n n n l l l z l l l l k l l n n n n n n l l l b l l l l k l l n l l z l n l l

z k k l N k k k k t z z k k l k k k k k t z z M N l N t z M N k

 

     

             



   

  




        




        



  A A A( )

       

3 1/ 2 1 1 1 1 1 1/ 2 1 1 1 2 2 2 1/ 2 3 1/ 2 1 1 1/ 2 2 2 1 1 1 1/ 2 1/ 2

( ) (2 ) 6 (2 ) 2 ( ) ( ) (2 6 (2 ) 6

n n n l l l l n n n l l n l l n n l b l l l n n n l l l n l l

z c k l k k k k z M N z l k t c k l k k z M N k k k

 

   



      








         



                        




           A

1 1 ) n

             

( , , )

b z

l k N 

P P P P

SLIDE 44

44

Turbulence closure

     

1 1 1 1 1 1 1 1 1/ 2 1 1 1 1 1 1 1 1 1/ 2 2 2 1 1 3 1/ 2 1/ 2 1

( ) 2 2 ( ) 6 ( ) 2 2 ( ) 6 2 6

n n n n n n l l l z l l l l l l n n n n n n l l l b l l l l l l n n l l l l z

z l N t z z l k t z z c M c N k z l N t

   

                

             



  

  




        




               



A A A( )

       

2 2 1 1 1 3 1 1/ 2 1/ 2 3 1/ 2 1 1 1 1 2 1 1/ 2 2 2 1 3 1/ 2 1/ 2 2 2 1 3 1/ 2 1/ 2

(2 ) ( ) (2 ) 6 2 ( ) ( 6

n n n l l n l l n n n l wall l l l n n l l l n l b n l l

c M c N k z c c k l F z c M c N k z l k t c M c N k

       

          

    



   



           -                                



 A

 

1 1 1 3 1/ 2 1 1 1 2 1 1/ 2

2 ) ( ) (2 ) 6

n n l l n n n l wall l l l

z c c k l F

 

   

 







           -                     

( , , )

b z

l k N 

Q Q Q Q

SLIDE 45

45

GOTM turbulence library (itur=4)



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? Newer versions may have resolved this issue



Works better than itur=3 otherwise

www.gotm.net

SLIDE 46

Radiation stress

Sxx etc are functions of wave energy, which is integrated over all freq. and direction 46

SLIDE 47

WBBL (Grant-Madsen)

angle between bottom current and dominant wave direction

tb: current induced bottom stress

Uw: orbital vel. amplitude w: representative angular freq.

Ratio of shear stresses Wave-current friction factor The nonlinear eq. system is solved with a simple iterative scheme starting from =0 (pure wave), c=1. Strong convergence is observed for practical applications 47

SLIDE 48

48

Hydraulics

The 2 ‘upstream’ and ‘downstream’ reference nodes on each side of the structure are used to calculate the thru flow Q 1 2 Q Faces (b) 1 2 3 4 5 5’ 4’ 3’ 2’ 1’ 1-1’ etc are node pairs

SLIDE 49

49

Hydraulics

The differential-integral eq. for continuity eq. is modified as And for the momentum eq.

Add b as a horizontal boundary in backtracking
At b and also for all internal sides in b, the side vel. is imposed

from Q – like an open bnd (add to I3 and I5)

SLIDE 50

50 solver

Preparation (Mom+cont.)

backtracking

SCHISM flow chart

Initialization

r hot start

Completed? yes no

forcings

Turbulence closure w

Transport

utput

Momentum eq. Update levels & eqstate

end Most expensive part

SLIDE 51

51

Numerical stability



Semi-implicitness circumvents CFL (most severe) (Casulli and Cattani 1994)



ELM bypasses Courant number condition for advection (actually CFL>0.4)



Implicit TVD2 transport scheme along the vertical bypasses Courant number condition in the vertical



Explicit terms



Baroclinicity  internal Courant number restriction (usually not a problem)



Horizontal viscosity/diffusivity  diffusion number condition (mild)



Upwind, TVD2 transport or WEON3: horizontal Courant number condition  sub-cycling



Coriolis – “inertial modes” in eddying regime (Le Roux 2009; Danilov 2013); needs to be stabilized by viscosity/filter

2

0.5 t x    ' 1 g h t x   

1 1 1 1 2 3

2

n n n

u u u

  



SLIDE 52

52

Lateral b.c. and nudging

Variable Type 1 (*.th) Type 2: Type 3 Type 4 (*[23]D.th) Type 5 Type -1 Type -4, -5 (uv3D.th); nudging

Nudging/spon ge layer near boundary

 elev.th: Time history; uniform along bnd constant Tidal amp/phases elev2D.th.n c: time- and space- varying along bnd elev2D.th.n c: combination

f 3&4

Must =0 S&T, tracers

[MOD]_?.th

: relax to time history (uniform along bnd) for inflow Relax to constant for inflow Relax to i.c. for inflow [MOD]_3D. th.nc: relax to time- and space- varying values along bnd during inflow u,v flux.th: via discharge (<0 for inflow!) Via discharge (<0 for inflow!) Tides (uniform amp/phase along bnd) uv3D.th.nc: time- and space- varying along bnd (in lon/lat for ics=2) uv3D.th.nc: combination

f 3&4 (but

tidal amp/phases vary along bnd) Flather (‘0’ for ) Relax to uv3D.th.nc (2 separate relaxations for in &

utflow)

bctides.in param.in

SLIDE 53

4H’s of SCHISM

Hybrid FE-FV formulation
Hybrid horizontal grid
Hybrid vertical grid
Hybrid code (openMP-MPI)