SCHISM physical formulation Joseph Zhang Advection 2 = + - - PowerPoint PPT Presentation

1

SCHISM physical formulation

Joseph Zhang

Advection

𝐸𝑑 𝐸𝑢 = 𝜖𝑑 𝜖𝑢 + 𝒗 ∙ 𝛼𝑑 = 0 u : advective velocity (usually known)

2

Diffusion

𝜖𝑑 𝜖𝑢 = 𝜖 𝜖𝑨 𝜆 𝜖 𝜖𝑨 k : diffusivity [m2/s]

3

Advection-Diffusion

𝐸𝑑 𝐸𝑢 = 𝜖 𝜖𝑨 𝜆 𝜖 𝜖𝑨

4

Dispersion

No dispersion With dispersion

• Waves of different frequencies travel at different phase speeds
• Usually related to derivatives of odd order (of either physical or

numerical origin)

𝜒𝑢 + 𝜒𝑦𝑦𝑦 − 6𝜒𝜒𝑦 = 0

5

6



Continuity equation



Momentum equations

Governing equations: Reynolds-averged Navier-Stokes (with vegetation)

0, ( ( , )) w u v z        u u





     





h dz

t u

MSL z x East

1 ˆ + ; ( )

A z

D g g f g p d Dt z z



                                 



u u f f k u u

𝐸𝒗 𝑒𝑢 = 𝐠 − 𝑕𝛼𝜃 + 𝒏𝒜 − 𝛽 𝒗 𝒗𝑀 𝑦, 𝑧, 𝑨

Vertical b.c. (3D):

𝜉 𝜖𝒗 𝜖𝑨 = 𝝊𝑥, 𝑨 = 𝜃 𝜉 𝜖𝒗 𝜖𝑨 = 𝜓𝒗𝑐, 𝑨 = −ℎ 𝒏𝑨 = 𝜖 𝜖𝑨 𝜉 𝜖𝒗 𝜖𝑨 , 3𝐸 𝝊𝒙 − 𝜓𝒗 𝐼 , 2𝐸 𝑀 𝑦, 𝑧, 𝑨 = ቊℋ 𝑨𝑤 − 𝑨 , 3𝐸 1, 2𝐸 𝛽(𝑦, 𝑧) = 𝐸𝑤𝑂𝑤𝐷𝐸𝑤/2 (vegetation)

z=zv

vegetation

7



Equation of state



Transport of salt and temperature



Turbulence closure: Umlauf and Burchard 2003

Governing equations: Reynolds-averged Navier-Stokes (with SAV)

) , , ( T S p   

( ), ( , )

h

Dc c Q c c S T Dt z z k k                

shear stratification

 

,

n m p k

c 



  𝐸𝑙 𝐸𝑢 = 𝜖 𝜖𝑨 𝜉𝑙

𝜔 𝜖𝑙

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

vegetation

8

On spherical coordinates

We transform the coordinates instead of eqs There are 2 frames used: 1. Global 2. Lon/lat (local @ node/element/side) Changes to the code is mainly related to re-project vectors and to calculate distances O xg l f P l0 (x) f0 (y) yg zg z

• Avoids polar singularity easily
• Preserves all original matrix properties

9

Continuity equation

h h

dz w

     

 



u

( , , )

x y t

z x y t dz dx dy w u v dt dt x dt y t                     

Vertical boundary conditions: kinematic z surface bottom

( , )

x y

z h x y dz w uh vh dt      

Integrated continuity equation

( )

t x y x y h h h

dz u v uh vh dz dz h

  

   

  

               

  

u u u u

h

dz t







     



u

10

Hydrostatic model

Hydrostatic assumption

1

A z

P g P g d P z

 

         



Separation of horizontal and vertical scales

1 ( ) + ( )

A z

D d P f Dt z z

  

                      



u u u k u

Boussinesq assumption (=>incompressibility)

z z

d d

 

          

 

1 ˆ + ( )

A z

D g g f g p d Dt z z



                                



u u k u u

11

Momentum equation: vertical boundary condition (b.c.)

at =

w

z z      u τ | | , at

D b b

C z h z       u u u

Surface Bottom



Logarithmic law:



Reynolds stress:



Turbulence closure:



Reynolds stress (const.)



Drag coefficient:  

ln ( ) / , ( ) ln( / ) : bottom roughness

b b b

z h z z h z h z z         u u ( )ln( / )

b b

z z h z        u u

1/ 2 2 2 / 3 2 1

2 , , 1 | | 2 ( )

m m D b

s K l s g K B C l z h  k      u

1/ 2

| | , ( ) ln( / )

D b b b b

C z h z h z z k           u u u

2

1 ln

b D

C z  k



      

b

ub

z=h

12

Bottom drag

• Use of z0 seems most natural option, but tends to over-estimate CD in shallow area
• CD should vary with bottom layer thickness (i.e. vertical grid)
• Different from 2D, 3D results of elevation depend on vertical grid (with constant CD)

meters

13

Transport equation: source terms

( ) , S S E P z z k       

Precipitation and evaporation model E: evaporation rate (kg/m2/s) (calculated from heat exchange model) P: Precipitation rate (kg/m2/s) (measured) Heat exchange model

,

p

T H z z C k       ( ) H IR IR S E      

IR’s are down/upwelling infrared (LW) radiation at surface S is the turbulent flux of sensible heat (upwelling) E is the turbulent flux of latent heat (upwelling) SW is net downward solar radiation The solar radiation is penetrative (body force) with attenuation (which depends upon turbidity) acts as a heat source within the water.

1

p

SW Q C z    

14

Air-sea exchange

1. Momentum: near-surface winds apply wind-stress on surface (influences advection, location of density fronts)

• free surface height: variations in atmospheric pressure over the domain

have a direct impact upon free surface height; set up due to wind stresses 2. Heat: various components of heat fluxes (dependent upon many variables) determine surface heat budget

• shortwave radiation (solar) - penetrative
• longwave radiation (infrared)
• sensible heat flux (direct transfer of heat)
• latent heat flux (heating/cooling associated with

condensation/evaporation) 3. Mass: evaporation, condensation, precipitation act as sources/sinks of fresh water

1-2 are accomplished thru the heat/salt exchange model

15

Downwelling SW at the surface is forecast in NWP models - a function of time of year, time of day, weather conditions, latitude, etc [sflux_rad*.nc] Upwelling SW is a simple function of downwelling SW

Solar radiation

SW SW   

 is the albedo - typically depends on solar zenith angle and sea state Attenuation of SW radiation in the water column is a function of turbidity and depth D (Jerlov, 1968, 1976; Paulson and Clayson, 1977):

1 2

/ /

( ) (1 ) Re (1 R)e

D d D d

SW z SW 

 

        

R, d1 and d2 depend on water type; D is the distance from F.S.

Type R d1 (m) d2 (m) Jerlov I 0.58 0.35 23 Jerlov IA 0.62 0.60 20 Jerlov IB 0.67 1.00 17 Jerlov II 0.77 1.50 14 Jerlov III 0.78 1.40 7.9 Paulson and Simpson 0.62 1.50 20 Estuary 0.80 0.90 2.1

SW Depth (m)

I II III

‘6’ ‘7’

16

Infrared radiation

• Downwelling IR at the surface is forecast in NWP models - a function of air

temperature, cloud cover, humidity, etc [sflux_rad*.nc]

• Upwelling IR is can be approximated as either a broadband measurement solely

within the IR wavelengths (i.e., 4-50 μm), or more commonly the blackbody radiative flux e is the emissivity, ~1 s is the Stefan-Boltzmann constant Tsfc is the surface temperature

4 sfc

IR T es 

17

Turbulent Fluxes of Sensible and Latent Heat

In general, turbulent fluxes are a function of:

• Tsfc,, T

air

• near-surface wind speed
• surface atmospheric pressure
• near-surface humidity

Scales of motion responsible for these heat fluxes are much smaller than can be resolved by any operational model - they must be parameterized (i.e., bulk aerodynamic formulation)

where: u* is the friction scaling velocity T* is the temperature scaling parameter q* is the specific humidity scaling parameter a is the surface air density Cpa is the specific heat of air Le is the latent heat of vaporization The scaling parameters are defined using Monin-Obukhov similarity theory, and must be solved for iteratively (i.e., Zeng et al., 1998).

* * * *

' ' ' '

a pa a pa a e a e

S C w T C u T E L w q L u q            

18

Wind Shear Stress: Turbulent Flux of Momentum Calculation of shear stress follows naturally from calculation of turbulent heat fluxes Total shear stress of atmosphere upon surface:

2 * a w

u    

Alternatively, Pond and Picard’s formulation can be used as a simpler option

| |

a w ds w w

C    τ u u

' '

0.61 0.063 1000 max(6, min(50, ))

w ds w w

u C u u   

19

Atmospheric inputs to heat/momentum exchange model

• Forecasts of uw vw @ 10m, PA @ MSL,T

air and qair @2m

• used by the bulk aerodynamic model to calculate the scaling parameters
• Forecasts of downward IR and SW

data source supplying agency time period spatial resolution temporal resolution area of coverage data type MRF NCEP 04/01/2001-02/29/2004 1° x 1° 12-hour snapshots 129W-120W, 35N-51N forecast GFS NCEP 07/03/2003-present 1° x 1° 3-hour snapshots 180W-70W, 90S-90N forecast OSU-ARPS OSU 05/04/2001-02/25/2004 12 km 1-hour snapshots 128W-119W, 41N-47N (approx) forecast ETA/NAM NCEP 07/03/2003-present 12 km 3-hour snapshots Eta Grid 218, west of 100W forecast NCAR/NCEP Reanalysis (NARR) NCAR 01/01/1979-12/31/2006 12km 6-hour snapshots North America reanalysis

Atmospheric properties (“air”)

data source supplying agency time period spatial resolution temporal resolution area of coverage data type AVN (lo-res) NCEP 04/01/2001-10/28/2002 0.7° x 0.7° (approx) 3-hour averages 129W-120W, 35N-51N forecase AVN (hi-res) NCEP 10/29/2002-02/29/2004 0.5° x 0.5° (approx) 3-hour averages 129W-120W, 35N-51N forecast GFS NCEP 07/03/2003-present 0.5° x 0.5° (approx) 3-hour averages 180W-70W, 90S-90N forecast ETA/NAM NCEP 07/03/2003-present 12 km 3-hour snapshots Eta Grid 218, west of 100W forecast NCAR/NCEP Reanalysis (NARR) NCAR 01/01/1979-12/31/2006 12km 6-hour averages North America reanalysis

Heat fluxes (“rad”)

20

Turbulence closure: Umlauf and Burchard (2003)

• Use a generic length-scale variable to represent various closure schemes
• Mellor-Yamada-Galperin (with modification to wall-proximity function  k-kl)
• k-e (Rodi)
• k-w (Wilcox)
• KPP (Large et al.; Durski et al)

 

p m n

c k



 

 

3 3/ 2 1

c k



e





1/ 2

c k



 

' 1/ 2

c k



k 

k k  

  s 

 

  s 

Stability: Kantha and Clayson (smoothes Galperin’s stability function as it approaches max)

Diffusivity

'

2 , 2 0.4939 1 30.19

m h h h

c s c s s G

 

    0.392 17.07 1 6.127

h h m h

s G s G   

2 _ _ _ _ _ _

( ) 0.28 0.023 2 0.02

h u h u h c h h u h h c h c

G G G G G G G G         

Wall function

2 2 2 4

1

wall b s

l l F E E d d k k                𝐸𝑙 𝐸𝑢 = 𝜖 𝜖𝑨 𝜉𝑙

𝜔 𝜖𝑙

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

4

Hydraulics

Orifice equation

4

Hydraulics

Radial gate weir