Model help:ROMS: Difference between revisions
m (WikiSysop moved page Temp:ROMS to Model help:ROMS without leaving a redirect) |
m (Text replacement - "http://csdms.colorado.edu/wiki/" to "https://csdms.colorado.edu/wiki/") |
||
| Line 3: | Line 3: | ||
1) Log in to the wiki | 1) Log in to the wiki | ||
2) Create a new page for each model, by using the following URL: | 2) Create a new page for each model, by using the following URL: | ||
* | * https://csdms.colorado.edu/wiki/Model help:<modelname> | ||
* Replace <modelname> with the name of a model | * Replace <modelname> with the name of a model | ||
3) Than follow the link "edit this page" | 3) Than follow the link "edit this page" | ||
| Line 1,584: | Line 1,584: | ||
<span class="remove_this_tag">Follow the next steps to include images / movies of simulations:</span> | <span class="remove_this_tag">Follow the next steps to include images / movies of simulations:</span> | ||
* <span class="remove_this_tag">Upload file: | * <span class="remove_this_tag">Upload file: https://csdms.colorado.edu/wiki/Special:Upload</span> | ||
* <span class="remove_this_tag">Create link to the file on your page: <nowiki>[[Image:<file name>]]</nowiki>.</span> | * <span class="remove_this_tag">Create link to the file on your page: <nowiki>[[Image:<file name>]]</nowiki>.</span> | ||
Latest revision as of 17:15, 19 February 2018
ROMS
ROMS is a community model shared by a large user group around the world, with applications ranging from the study of entire ocean basins to coastal sub-regions.
Model introduction
ROMS incorporates advanced features and high-order numerics, allowing efficient and robust resolution of mesoscale dynamics in the oceanic and coastal domains. The model solves the free surface, hydrostatic, primitive equations of the fluid dynamics over variable topography using stretched; terrain-following coordinates in the vertical and orthogonal, curvilinear coordinates in the horizontal. This allows enhancement of the spatial resolution in the regions of interest.
Model parameters
Uses ports
This will be something that the CSDMS facility will add
Provides ports
This will be something that the CSDMS facility will add
Main equations
1) Equations of motion
a) Momentum balance in the x-directions
[math]\displaystyle{ {\frac{\partial u}{\partial t}} + \vec{v} \nabla u - fv = -{\frac{\partial \phi}{\partial x}} - {\frac{\partial}{\partial z}} \left ( \overline{u'w'} - \nu {\frac{\partial u}{\partial z}} \right ) + F_{u} + D_{u} }[/math] (1)
b) Momentum balance in the y-directions
[math]\displaystyle{ {\frac{\partial v}{\partial t}} + \vec{v} \nabla v + fu = -{\frac{\partial \phi}{\partial y}} - {\frac{\partial}{\partial z}} \left ( \overline{v'w'} - \nu {\frac{\partial v}{\partial z}} \right )+ F_{v} + D_{v} }[/math] (2)
c) Advective-diffusive equations for a scalar concentration field C(x,y,z,t) (e.g. salinity, temperature, or nutrients)
[math]\displaystyle{ {\frac{\partial C}{\partial t}} + \vec{v} \nabla C = {\frac{\partial}{\partial z}} \left ( \overline{C'w'} - \nu_{\theta} {\frac{\partial C}{\partial z}} \right )+ F_{C} + D_{C} }[/math] (3)
d) Equation of state
[math]\displaystyle{ \rho = \rho \left (T,S,P \right ) }[/math] (4)
e) Vertical momentum equation (assume that the vertical pressure gradient balances the buoyancy force)
[math]\displaystyle{ {\frac{\partial \phi}{\partial z}} = {\frac{- \rho g}{\rho_{o}}} }[/math] (5)
f) continuity equation for an incompressible fluid
[math]\displaystyle{ {\frac{\partial u}{\partial x}} + {\frac{\partial v}{\partial y}} + {\frac{\partial w}{\partial z}} = 0 }[/math] (6)
g) Reynolds stresses (an overbar represents a time average and a prime represents a fluctuation above the mean)
[math]\displaystyle{ \overline{u'w'} = - K_{M} {\frac{\partial u}{\partial z}} }[/math] (7)
[math]\displaystyle{ \overline{v'w'} = - K_{M} {\frac{\partial v}{\partial z}} }[/math] (8)
[math]\displaystyle{ \overline{C'w'} = - K_{C} {\frac{\partial C}{\partial z}} }[/math] (9)
2) Vertical boundary conditions
a) top boundary condition ( z = ζ (x,y,t ))
[math]\displaystyle{ K_{m} {\frac{\partial u}{\partial z}} = \tau_{S}^x \left (x,y,t \right ) }[/math] (10)
[math]\displaystyle{ K_{m} {\frac{\partial v}{\partial z}} = \tau_{S}^y \left (x,y,t \right ) }[/math] (11)
[math]\displaystyle{ K_{C} {\frac{\partial C}{\partial z}} = {\frac{Q_{C}}{\rho_{o} c_{P}}} }[/math] (12)
[math]\displaystyle{ w = {\frac{\partial \zeta}{\partial t}} }[/math] (13)
b) bottom boundary condition (z = -h(x,y))
[math]\displaystyle{ K_{m}{\frac{\partial u}{\partial z}} = \tau _{b}^x \left (x,y,t \right ) }[/math] (14)
[math]\displaystyle{ K_{m}{\frac{\partial v}{\partial z}} = \tau _{b}^y \left (x,y,t \right ) }[/math] (15)
[math]\displaystyle{ K_{C}{\frac{\partial C}{\partial z}} = 0 }[/math] (16)
[math]\displaystyle{ -w + \vec{v} \nabla h = 0 }[/math] (17)
c) Horizontal boundary conditions (Eastern and western boundary)(considering biharmonic friction)
[math]\displaystyle{ {\frac{\partial}{\partial x}} \left ( {\frac{h \nu}{mn}}{\frac{\partial^2 u}{\partial x^2}} \right ) = 0 }[/math] (18)
Horizontal boundary conditions (Northern and southern boundary)
[math]\displaystyle{ {\frac{\partial}{\partial y}} \left ( {\frac{h \nu}{mn}}{\frac{\partial^2 u}{\partial y^2}} \right ) = 0 }[/math] (19)
3) Vertical transformation equations
a) vertical coordinate transformation (been available in ROMS since 1999)
[math]\displaystyle{ z \left (x,y,\sigma, t \right ) = S \left (x,y,\sigma \right )+ \zeta \left (x,y,t \right ) [1 + {\frac{S \left (x,y,\sigma \right )}{h \left (x,y\right )}}] }[/math] (20)
[math]\displaystyle{ S \left (x,y,\sigma \right ) = h_{c} \sigma + [h \left (x,y\right ) - h_{c}] C \left (\sigma \right) }[/math] (21)
[math]\displaystyle{ S \left (x,y,\sigma \right ) =\left\{\begin{matrix} 0 & \sigma = 0, C \left (\sigma \right ) = 0 \\ - h \left (x,y \right ) & \sigma = -1, C \left (\sigma \right ) = -1 \end{matrix}\right. }[/math] (22)
b) vertical coordinate transformation (been available in ROMS since 2005)
[math]\displaystyle{ z \left (x,y,\sigma,t \right ) = \zeta \left (x,y,t\right ) + [\zeta \left (x,y,t\right ) + h \left (x,y \right )] S \left (x,y,\sigma \right ) }[/math] (23)
[math]\displaystyle{ S \left (x,y,\delta \right ) = {\frac{h_{c} \delta + h \left (x,y \right ) C \left ( \delta \right )}{h_{c} + h \left (x,y\right )}} }[/math] (24)
[math]\displaystyle{ S \left (x,y,\sigma \right ) =\left\{\begin{matrix} 0 & \sigma = 0, C \left (\sigma \right ) = 0 \\ -1 & \sigma = -1, C \left (\sigma \right ) = -1 \end{matrix}\right. }[/math] (25)
unperturbed depth:
[math]\displaystyle{ \hat{z} \left (x,y,\sigma \right ) \equiv h \left (x,y \right ) S \left (x,y,\sigma \right ) = h \left (x,y\right )[{\frac{h_{c}\sigma + h \left (x,y\right ) C \left (\sigma \right )}{h_{c} +h \left (x,y \right )}}] }[/math] (26)
[math]\displaystyle{ d \hat{z} = d \sigma h \left (x,y\right ) [{\frac{h_{c}}{h_{c} + h \left (x,y\right )}}] }[/math] (27)
4) Vertical stretching functions
a) Vstretching = 1 (Song and Haidvogel (1994))
[math]\displaystyle{ C \left (\sigma \right ) = \left (1-\theta_{B}\right ) {\frac{sinh \left (\theta_{S}\sigma \right )}{sinh \theta_{S}}} + \theta_{B}[{\frac{tanh[\theta_{S}\left (\sigma +{\frac{1}{2}}\right )]}{2 tanh \left ({\frac{1}{2}} \theta_{S}\right )}} - {\frac{1}{2}}] }[/math] (28)
b) Vstretching = 2 (A.Shchepetkin (2005))
[math]\displaystyle{ C \left (\sigma \right ) = \mu C_{sur}\left (\sigma\right ) + \left (1-\mu \right ) C_{bot} \left (\sigma \right ) }[/math] (29)
[math]\displaystyle{ C_{sur} \left (\sigma \right ) = {\frac{1 - cosh \left (\theta_{S} \sigma \right )}{cosh \left (\theta_{S} \right ) - 1}} }[/math] (30)
[math]\displaystyle{ C_{bot} \left (\sigma \right ) = {\frac{sinh[\theta_{B}\left (\sigma + 1 \right )]}{sinh \left (\theta_{B}\right )}} - 1 }[/math] (31)
[math]\displaystyle{ \mu = \left (\sigma +1 \right )^\alpha [1 + {\frac{\alpha}{\beta}} \left ( 1 - \left ( \sigma + 1 \right )^\beta \right )] }[/math] (32)
c) Vstretching = 3 (R.Geyer function for high bottom layer resolution in relatively shallow applications)
[math]\displaystyle{ C \left (\sigma \right ) = \mu C_{bot} \left (\sigma \right ) + \left (1 - \mu \right ) C_{sur} \left (\sigma \right ) }[/math] (33)
[math]\displaystyle{ C_{sur} \left (\sigma \right ) = - {\frac{log[cosh \left (\gamma \right ) abs \left (\sigma \right ) ^\left (\theta_{S}\right )]}{log [cosh \left (\gamma \right )]}} }[/math] (34)
[math]\displaystyle{ C_{bot} \left (\sigma \right ) = {\frac{log [cosh \left (\gamma \right ) \left (\sigma + 1 \right )^ \left (\theta_{B}\right )]}{log [cosh \left (\gamma \right )]}} - 1 }[/math] (35)
[math]\displaystyle{ \mu = {\frac{1}{2}} [1 - tanh \left (\gamma \left (\sigma + {\frac{1}{2}}\right )\right )] }[/math] (36)
d) Vstretching = 4 (A. Shchepetkin (2010) UCLA-ROMS current function) Surface refinement function
[math]\displaystyle{ \left\{\begin{matrix} C \left ( \sigma \right ) = {\frac{1 - cosh \left (\theta_{S}\sigma \right )}{cosh \left (\theta_{S}\right ) - 1}} & \theta_{S} \gt 0 \\ - \sigma ^2 & \theta_{S} \lt = 0 \end{matrix}\right. }[/math] (37)
Bottom refinement function (θB > 0)
[math]\displaystyle{ C \left (\sigma \right ) = {\frac{exp \left (\theta_{B} C\left (\sigma \right )\right ) - 1}{1 - exp \left (- \theta_{B}\right )}} }[/math] (38)
5) Terrain-following coordinate transformation (the dynamic equations and boundary conditions are transformed based on the following equations)
[math]\displaystyle{ \hat{x} = x }[/math] (39)
[math]\displaystyle{ \hat{y} = y }[/math] (40)
[math]\displaystyle{ \sigma = \sigma \left (x,y,z \right ) }[/math] (41)
[math]\displaystyle{ z = z \left (x,y,\sigma \right ) }[/math] (42)
[math]\displaystyle{ \hat{t} = t }[/math] (43)
[math]\displaystyle{ \left ({\frac{\partial}{\partial x}}\right )_{z} = \left ({\frac{\partial}{\partial x}}\right ) - \left ({\frac{1}{H_{z}}}\right ) \left ({\frac{\partial z}{\partial x}}\right )_{\sigma} {\frac{\partial}{\partial \sigma}} }[/math] (44)
[math]\displaystyle{ \left ({\frac{\partial}{\partial y}} \right )_{z} = \left ({\frac{\partial}{\partial y}}\right )_{\sigma} - \left ({\frac{1}{H_{z}}}\right ) \left ({\frac{\partial z}{\partial y}}\right )_{\sigma} {\frac{\partial}{\partial \sigma}} }[/math] (45)
[math]\displaystyle{ {\frac{\partial}{\partial z}} = \left ({\frac{\partial \sigma}{\partial z}}\right ) {\frac{\partial}{\partial \sigma}} = {\frac{1}{H_{z}}} {\frac{\partial}{\partial \sigma}} }[/math] (46)
6) Curvilinear coordinate transformation
relationship of horizontal arc length to the differential distance
[math]\displaystyle{ \left (ds \right )_{\epsilon} = \left ({\frac{1}{m}} \right ) d \epsilon }[/math] (47)
[math]\displaystyle{ \left (ds \right )_{\eta} = \left ({\frac{1}{n}}\right ) d \eta }[/math] (48)
velocity components in this coordinate system (equations of motion can be re-written with the following equaitons)
[math]\displaystyle{ \overrightarrow{v} \hat{\epsilon} = u }[/math] (49)
[math]\displaystyle{ \overrightarrow{v} \hat{\eta} = v }[/math] (50)
7) Time-stepping schemes
Assuming that multiple equations for any number of variables are time stepped synchronously
[math]\displaystyle{ {\frac{\partial \phi \left ( t \right )}{\partial t}} = F \left ( t \right ) }[/math] (51)
a) Euler time step
[math]\displaystyle{ {\frac{\phi \left( t + \Delta t \right ) - \phi \left (t\right )}{\Delta t}} = F \left (t\right ) }[/math] (52)
b) Leapfrog
[math]\displaystyle{ {\frac{\phi \left ( t + \Delta t \right ) - \phi \left ( t - \Delta t \right )}{2 \Delta t}} = F \left (t \right ) }[/math] (53)
[math]\displaystyle{ {\frac{\phi \left ( t + \Delta t \right ) - \phi \left ( t \right )}{\Delta t}} = {\frac{1}{2}} [F \left ( t \right ) + F ^* \left ( t + \Delta t \right )] }[/math] (54)
c) Third-order Adams-Bashforth (AB3)
[math]\displaystyle{ {\frac{\phi \left ( t + \Delta t \right ) - \phi \left ( t \right )}{\Delta t}} = \alpha F \left ( t \right ) + \beta F \left ( t - \Delta t \right ) + \gamma F \left ( t - 2 \Delta t \right ) }[/math] (55)
d) Forward-Backward
[math]\displaystyle{ {\frac{\partial \zeta}{\partial t}} = F \left (u \right ) }[/math] (56)
[math]\displaystyle{ {\frac{\partial u}{\partial t}} = G \left (\zeta \right ) }[/math] (57)
[math]\displaystyle{ \zeta^ \left (n + 1 \right ) = \zeta ^n + F \left ( u^n \right ) \Delta t }[/math] (58)
[math]\displaystyle{ u^ \left (n + 1 \right ) = u^n + G \left (\zeta^ \left (n + 1 \right ) \right ) \Delta t }[/math] (59)
e) Forward-Backward Feedback (RK2-FB) (one option for solving equation (56), (57))
[math]\displaystyle{ \zeta^ \left (n+1,* \right ) = \zeta^n + F \left (u^n \right ) \Delta t }[/math] (60)
[math]\displaystyle{ u^ \left (n + 1, * \right ) = u^n + [\beta G \left ( \zeta^ \left (n + 1, * \right )\right ) + \left ( 1 - \beta \right ) G \left (\zeta^n \right )] \Delta t }[/math] (61)
[math]\displaystyle{ \zeta^ \left (n+1 \right ) = \zeta^n + {\frac{1}{2}}[F \left (u^ \left (n+1, * \right ) \right ) + F \left ( u^n \right )] \Delta t }[/math] (62)
[math]\displaystyle{ u^ \left (n+1 \right ) = u^n + {\frac{1}{2}}[\epsilon G \left ( \zeta ^ \left (n+1 \right )\right ) + \left ( 1 - \epsilon \right ) G \left ( \zeta^ \left (n+1, *\right ) \right ) + G \left ( \zeta^n \right )] \Delta t }[/math] (63)
f) LF-TR and LF-AM3 with FB Feedback (another option to solve equation (56), (57))
[math]\displaystyle{ \zeta^ \left (n+1,* \right ) = \zeta^ \left (n-1\right ) + 2 F \left (u^n \right ) \Delta t }[/math] (64)
[math]\displaystyle{ u^ \left (n+1,*\right ) = u^\left (n-1\right ) + 2 \{ \left (1 - 2 \beta \right ) G \left (\zeta^n \right ) + \beta [G \left (\zeta^ \left (n+1,*\right )\right ) + G \left (\zeta^\left (n-1\right )\right )]\} \Delta t }[/math] (65)
[math]\displaystyle{ \zeta^\left (n+1\right ) = \zeta^n + [\left ({\frac{1}{2}} - \gamma \right ) F \left (u^ \left (n+1,*\right )\right ) + \left ({\frac{1}{2}} + 2 \gamma \right ) F \left (u^n \right ) - \gamma F \left (u^ \left (n-1\right )\right )] \Delta t }[/math] (66)
[math]\displaystyle{ u^ \left (n+1\right ) = u^n + \{ \left ({\frac{1}{2}} - \gamma \right ) [\epsilon G \left ( \zeta^ \left (n+1\right )\right ) + \left (1 - \epsilon \right ) G \left (\zeta^ \left (n+1,* \right ) \right )] + \left ({\frac{1}{2}} + 2 \gamma \right ) G \left (\zeta^n \right ) - \gamma G \left (\zeta^ \left (n-1\right )\right ) \} \Delta t }[/math] (67)
g) Generalized FB with an AB3-AM4 Step
[math]\displaystyle{ \zeta^ \left (n+1\right ) = \zeta^n + [\left ({\frac{3}{2}} + \beta \right ) F \left (u^n\right ) - \left ({\frac{1}{2}} + 2\beta \right ) F \left (u^ \left (n-1\right ) \right ) + \beta F \left (u^ \left (n-2\right )\right )] \Delta t }[/math] (68)
[math]\displaystyle{ u^ \left (n+1\right ) = u^n + [\left ({\frac{1}{2}} + \gamma + 2\epsilon \right ) G \left ( \zeta^ \left (n+1\right )\right ) + \left ( {\frac{1}{2}} - 2 \gamma - 3 \epsilon \right ) G \left (\zeta^n \right ) + \gamma G \left (\zeta^ \left (n-1\right )\right ) + \epsilon G \left (\zeta^ \left (n-2\right ) \right )] \Delta t }[/math] (69)
8) Conservation Properties
a) A tracer concentration equation in advective form
[math]\displaystyle{ {\frac{\partial C}{\partial t}} + \left ( u \Delta \right ) C = 0 }[/math] (70)
b) A tracer concentration equation in conservation form
[math]\displaystyle{ {\frac{\partial C}{\partial t}} + \Delta \left ( u C \right ) = 0 }[/math] (71)
c) Continuity equation
[math]\displaystyle{ \left ( \Delta u \right ) = 0 }[/math] (72)
d) Semi-discrete form of the tracer equation
[math]\displaystyle{ {\frac{\partial}{\partial t}} \left ({\frac{H_{z} C}{m n}}\right ) + \sigma_{xi} \left ({\frac{u \overline{H_{z}}^ \eta \overline{C}^ \eta}{\overline{n}^ \eta}}\right ) + \sigma_{\eta} \left ({\frac{v \overline{H_{z}}^ \eta \overline{C}^ \eta}{\overline{m}^ \eta}}\right ) + \sigma_{delta} \left ( \overline{C}^ \delta {\frac{H_{z} \Omega}{mn}}\right ) = {\frac{1}{mn}}{\frac{\partial}{\partial \delta}} \left ({\frac{K_{m}}{\Delta z}}{\frac{\partial C}{\partial \delta}} \right ) + D_{C} + F_{C} }[/math] (73)
e) Volume-averaged concentration over the grid box ΔV
[math]\displaystyle{ C = {\frac{mn}{H_{z}}} \int_{\Delta V} {\frac{H_{z}C}{mn}} \sigma_{\xi} \sigma_{\eta} \sigma{\delta} }[/math] (74)
9) Depth-Integrated Equations
a) Depth average of quantity A
[math]\displaystyle{ \overline{A} = {\frac{1}{D}} \int_{-1}^0 H_{z} A d \delta }[/math] (75)
b) total depth of the water volume
[math]\displaystyle{ D \equiv \zeta \left (\xi, \eta, t \right ) + h \left (\xi, \eta \right ) }[/math] (76)
10) Pressure gradient term in mode coupling
a) Pressure gradient
[math]\displaystyle{ -{\frac{g D}{n}}{\frac{\partial \zeta}{\partial \xi}} + [{\frac{gD}{n}} {\frac{\partial \zeta}{\partial \xi}} + F ] }[/math] (77)
b) Vertically integrated pressure gradient
[math]\displaystyle{ F = -{\frac{1}{\rho_{o} n}} \int_{-h}^ \zeta {\frac{\partial P}{\partial \xi}} dz = -{\frac{1}{\rho_{o}}}{\frac{g}{n}} \{{\frac{\partial}{\partial \xi}}\left ({\frac{\rho^* D^2}{2}}\right ) - \overline{\rho}D{\frac{\partial h}{\partial \xi}}\} = -{\frac{1}{\rho_{o}}}{\frac{g}{n}}D\{\rho^* {\frac{\partial \zeta}{\partial \xi}} + {\frac{D}{2}}{\frac{\partial \rho^*}{\partial \xi}} + \left (\rho^* - \overline{\rho} \right ) {\frac{\partial h}{\partial \xi}} \} }[/math] (78)
11) Horizontal and Vertical Advection
a) Advection of a tracer C
[math]\displaystyle{ {\frac{\partial}{\partial t}}{\frac{H_{z} C}{mn}} = -{\frac{\partial}{\partial \xi}}F^\xi - {\frac{\partial}{\partial \eta}}F^\eta - {\frac{\partial}{\partial \sigma}}F^\sigma }[/math] (79)
b) Advective fluxes (definition - Second-order centered and Fourth-order Centered - Third-order Upwind)
[math]\displaystyle{ F^\xi = {\frac{H_{z} u C}{n}} = {\frac{\overline{H_{z}}^\xi u \overline{C}^\xi}{\overline{n}^\xi}} = {\frac{\overline{H_{z}}^ \xi}{\overline{n}^\xi}}u \left (\overline{C}^ \xi - {\frac{1}{3}}{\frac{\partial G^\xi}{\partial \xi}}\right ) = {\frac{H_{z} u}{n}} \left (C - \gamma {\frac{\partial^2 C}{\partial \xi^2}} \right ) }[/math] (80)
[math]\displaystyle{ F^\eta = {\frac{H_{z} v C}{m}} = {\frac{\overline{H_{z}}^\eta v \overline{C}^\eta}{\overline{m}^\eta}} = {\frac{\overline{H_{z}}^\eta}{\overline{m}^\eta}} v \left ( \overline{C}^\eta - {\frac{1}{3}} {\frac{\partial G^\eta}{\partial \eta}}\right ) = {\frac{H_{z} v}{m}} \left ( C - \gamma {\frac{\partial^2 C}{\partial \eta^2}}\right ) }[/math] (81)
[math]\displaystyle{ F^\sigma = {\frac{H_{z} \Omega C}{mn}} = {\frac{\overline{H_{z}}^\delta \Omega \overline{C}^\delta}{mn}} = {\frac{\overline{H_{z}}^\delta}{mn}} \Omega \left (\overline{C}^\delta - {\frac{1}{3}}{\frac{\partial G^\delta}{\partial \delta}}\right ) }[/math] (82)
Centered fourth-order scheme in the vertical when the third-order upwind option is turned on
[math]\displaystyle{ F^s = {\frac{H_{z}w}{mn}}[-{\frac{1}{16}}C_{i,j,k} + {\frac{9}{16}}C_{i,j,k} + {\frac{9}{16}}C_{i,j,k+1} - {\frac{1}{16}}C_{i,j,k+2}] }[/math] (83)
Gradient terms (fourth-order centered and fourth-order Akima)
[math]\displaystyle{ G^\xi = \overline{\left ({\frac{\partial C}{\partial \xi}}\right )}^ \xi = 2 {\frac{\partial C}{\partial \xi_{i}}}{\frac{\partial C}{\partial \xi_{i+1}}} / \left ({\frac{\partial C}{\partial \xi_{i}}} + {\frac{\partial C}{\partial \xi_{i+1}}} \right ) }[/math] (84)
[math]\displaystyle{ G^\eta = \overline{\left ({\frac{\partial C}{\partial \eta}}\right )}^ \eta = 2 {\frac{\partial C}{\partial \eta_{j}}}{\frac{\partial C}{\partial \eta_{j+1}}} / \left ({\frac{\partial C}{\partial \eta_{j}}} + {\frac{\partial C}{\partial \eta_{j+1}}} \right ) }[/math] (85)
[math]\displaystyle{ G^\delta = \overline{\left ({\frac{\partial C}{\partial \delta}}\right )}^ \delta = 2 {\frac{\partial C}{\partial \delta_{k}}}{\frac{\partial C}{\partial \delta_{k-1}}} / \left ({\frac{\partial C}{\partial \delta_{k}}} + {\frac{\partial C}{\partial \delta_{k-1}}} \right ) }[/math] (86)
12) Horizontal Boundary Conditions
Radiation boundary condition
[math]\displaystyle{ {\frac{\partial \phi}{\partial t}} = - \left ( \phi_{xi} {\frac{\partial \phi}{\partial \xi}} + \phi_{\eta}{\frac{\partial \phi}{\partial \eta}} \right ) }[/math] (87)
[math]\displaystyle{ \phi_{\xi} = {\frac{F {\frac{\partial \phi}{\partial \xi}}}{\left ({\frac{\partial \phi}{\partial \xi}}\right )^2 + \left ({\frac{\partial \phi}{\partial \eta}}\right )^2}} }[/math] (88)
[math]\displaystyle{ \phi_{\eta} = {\frac{F {\frac{\partial \phi}{\partial \eta}}}{\left ({\frac{\partial \phi}{\partial \xi}}\right )^2 + \left ({\frac{\partial \phi}{\partial \eta}}\right )^2}} }[/math] (89)
[math]\displaystyle{ F = - {\frac{\partial \phi}{\partial t}} }[/math] (90)
13) Deviatory Stress Tensor
[math]\displaystyle{ F^u \equiv \hat{\xi} \left (\Delta \vec{\delta}\right ) = {\frac{mn}{H_{z}}} [{\frac{\partial}{\partial \xi}} \left ({\frac{H_{z} \sigma_{\xi \xi}}{n}} \right ) + {\frac{\partial}{\partial \eta}} \left ({\frac{H_{z} \delta_{\xi \eta}}{m}} \right ) + {\frac{\partial}{\partial s}} \left ({\frac{\delta_{\xi s}}{mn}}\right ) + H_{z} \delta_{\xi \eta} {\frac{\partial}{\partial \eta}}\left ({\frac{1}{m}}\right ) - H_{z} \delta_{\eta \eta}{\frac{\partial}{\partial \xi}} \left ({\frac{1}{n}}\right ) - {\frac{1}{n}} \delta_{ss} {\frac{\partial H_{z}}{\partial \xi}}] }[/math] (91)
[math]\displaystyle{ F^v \equiv \hat{\eta} \left (\Delta \vec{\delta}\right ) = {\frac{mn}{H_{z}}} [{\frac{\partial}{\partial \xi}} \left ({\frac{H_{z} \sigma_{\eta \xi}}{n}} \right ) + {\frac{\partial}{\partial \eta}} \left ({\frac{H_{z} \delta_{\eta \eta}}{m}} \right ) + {\frac{\partial}{\partial s}} \left ({\frac{\delta_{\eta s}}{mn}}\right ) + H_{z} \delta_{\eta \xi} {\frac{\partial}{\partial \xi}}\left ({\frac{1}{n}}\right ) - H_{z} \delta_{\xi \xi}{\frac{\partial}{\partial \eta}} \left ({\frac{1}{m}}\right ) - {\frac{1}{m}} \delta_{ss} {\frac{\partial H_{z}}{\partial \eta}}] }[/math] (92)
[math]\displaystyle{ \delta_{\xi \xi} = \left (A_{M} + \nu \right ) e_{\xi \xi} + \left ( \nu - A_{M} \right ) e_{\eta \eta} }[/math] (93)
[math]\displaystyle{ \delta_{\eta \eta} = \left (\nu - A_{M}\right ) e_{\xi \xi} + \left ( A_{M} + \nu \right ) e_{\eta \eta} }[/math] (94)
[math]\displaystyle{ \delta_{ss} = 2 \nu e_{ss} }[/math] (95)
[math]\displaystyle{ \delta_{\xi \eta} = \delta_{\eta \xi} = 2 A_{M} e_{\xi \eta} }[/math] (96)
[math]\displaystyle{ \delta_{\xi s} = 2 K_{M} e_{\xi s} }[/math] (97)
[math]\displaystyle{ \delta_{\eta s} = 2 K_{M} e_{\eta s} }[/math] (98)
[math]\displaystyle{ e_{\xi \xi} = m {\frac{\partial u}{\partial \xi}} + mnv {\frac{\partial}{\partial \eta}}\left ({\frac{1}{m}}\right ) }[/math] (99)
[math]\displaystyle{ e_{ss} = {\frac{1}{H_{z}}}{\frac{\partial \left (\omega H_{z}\right )}{\partial s}} + {\frac{m}{H_{z}}}u{\frac{\partial H_{z}}{\partial \xi}} + {\frac{n}{H_{z}}}v{\frac{\partial H_{z}}{\partial \eta}} }[/math] (100)
[math]\displaystyle{ 2e_{\xi \eta} = {\frac{m}{n}}{\frac{\partial \left (nv\right )}{\partial \xi}} + {\frac{n}{m}} {\frac{\partial \left (mu \right )}{\partial \eta}} }[/math] (101)
[math]\displaystyle{ 2e_{\xi s} = {\frac{1}{m H_{z}}}{\frac{\partial \left (mu\right )}{\partial s}} + m H_{z}{\frac{\partial \omega}{\partial \xi}} }[/math] (102)
[math]\displaystyle{ 2e_{\eta s} = {\frac{1}{n H_{z}}}{\frac{\partial \left (nv\right )}{\partial s}} + n H_{z}{\frac{\partial \omega}{\partial \eta}} }[/math] (103)
a) Transverse Stress Tensor
[math]\displaystyle{ H_{z}F^u = n^2 m {\frac{\partial}{\partial \xi}} \left ({\frac{H_{z} F^ \left (u \xi \right )}{n}}\right ) + m^2 n {\frac{\partial}{\partial \eta}} \left ({\frac{H_{z} F^ \left (u \eta \right )}{m}} \right ) }[/math] (104)
[math]\displaystyle{ H_{z}F^v = n^2 m {\frac{\partial}{\partial \xi}} \left ({\frac{H_{z} F^ \left (v \xi \right )}{n}}\right ) + m^2 n {\frac{\partial}{\partial \eta}} \left ({\frac{H_{z} F^ \left (v \eta \right )}{m}} \right ) }[/math] (105)
[math]\displaystyle{ F^ \left (u \xi \right ) = {\frac{1}{n}}A_{M}[{\frac{m}{n}}{\frac{\partial \left (nu\right )}{\partial \xi}} - {\frac{n}{m}}{\frac{\partial \left (mv\right )}{\partial \eta}}] }[/math] (106)
[math]\displaystyle{ F^ \left (u \eta \right ) = {\frac{1}{m}}A_{M}[{\frac{n}{m}}{\frac{\partial \left (mu\right )}{\partial \eta}} + {\frac{m}{n}}{\frac{\partial \left (nv\right )}{\partial \xi}}] }[/math] (107)
[math]\displaystyle{ F^ \left (v \xi \right ) = {\frac{1}{n}}A_{M}[{\frac{m}{n}}{\frac{\partial \left (nv\right )}{\partial \xi}} + {\frac{n}{m}}{\frac{\partial \left (mu\right )}{\partial \eta}}] }[/math] (108)
[math]\displaystyle{ F^ \left (v \eta \right ) = {\frac{1}{m}}A_{M}[{\frac{n}{m}}{\frac{\partial \left (mv\right )}{\partial \eta}} + {\frac{m}{n}}{\frac{\partial \left (nu\right )}{\partial \xi}}] }[/math] (109)
b) Rotated Transverse Stress Tensor
[math]\displaystyle{ H_{z}R^u = n^2 m {\frac{\partial}{\partial \xi}}\left ({\frac{H_{z} R^ \left (u \xi \right )}{n}}\right ) + m^2 n {\frac{\partial}{\partial \eta}}\left ({\frac{H_{z} R^ \left (u \eta \right )}{m}}\right ) + {\frac{\partial}{\partial s}}\left (R^ \left (us\right ) \right ) }[/math] (110)
[math]\displaystyle{ H_{z}R^v = n^2 m {\frac{\partial}{\partial \xi}}\left ({\frac{H_{z} R^ \left (v \xi \right )}{n}}\right ) + m^2 n {\frac{\partial}{\partial \eta}}\left ({\frac{H_{z} R^ \left (v \eta \right )}{m}}\right ) + {\frac{\partial}{\partial s}}\left (R^ \left (us\right ) \right ) }[/math] (111)
[math]\displaystyle{ R^ \left (u \xi \right ) = {\frac{1}{n}}A_{M}[{\frac{1}{n}} \left (m{\frac{\partial \left (nu\right )}{\partial \xi}} - m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nu \right )}{\partial s}}\right ) - {\frac{1}{m}} \left (n{\frac{\partial \left (mv\right )}{\partial \eta}} - n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mv\right )}{\partial s}}\right )] }[/math] (112)
[math]\displaystyle{ R^ \left (u \eta \right ) = {\frac{1}{m}}A_{M}[{\frac{1}{m}} \left (n{\frac{\partial \left (mu\right )}{\partial \eta}} - n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mu \right )}{\partial s}}\right ) + {\frac{1}{n}} \left (m{\frac{\partial \left (nv\right )}{\partial \xi}} - m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nv\right )}{\partial s}}\right )] }[/math] (113)
[math]\displaystyle{ R^ \left (us \right ) = m {\frac{\partial z}{\partial \xi}}A_{M}[{\frac{1}{n}}\left (m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nu \right )}{\partial s}} - m{\frac{\partial \left (nu\right )}{\partial \xi}}\right ) - {\frac{1}{m}}\left (n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mv\right )}{\partial s}} - n{\frac{\partial \left (mv\right )}{\partial \eta}}\right )] + n{\frac{\partial z}{\partial \eta}}A_{M}[{\frac{1}{m}}\left (n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mu\right )}{\partial s}} - n{\frac{\partial \left (mu \right )}{\partial \eta}}\right ) + {\frac{1}{n}}\left (m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nv \right )}{\partial s}} - m{\frac{\partial \left (nv\right )}{\partial \xi}}\right )] }[/math] (114)
[math]\displaystyle{ R^ \left (v \xi \right ) = {\frac{1}{n}}A_{M}[{\frac{1}{n}} \left (m{\frac{\partial \left (nv\right )}{\partial \xi}} - m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nv \right )}{\partial s}}\right ) + {\frac{1}{m}} \left (n{\frac{\partial \left (mu\right )}{\partial \eta}} - n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mu\right )}{\partial s}}\right )] }[/math] (115)
[math]\displaystyle{ R^ \left (v \eta \right ) = {\frac{1}{m}}A_{M}[{\frac{1}{m}} \left (n{\frac{\partial \left (mv\right )}{\partial \eta}} - n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mv \right )}{\partial s}}\right ) - {\frac{1}{n}} \left (m{\frac{\partial \left (nu\right )}{\partial \xi}} - m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nu\right )}{\partial s}}\right )] }[/math] (116)
[math]\displaystyle{ R^ \left (vs \right ) = m {\frac{\partial z}{\partial \xi}}A_{M}[{\frac{1}{n}}\left (m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nv \right )}{\partial s}} - m{\frac{\partial \left (nv\right )}{\partial \xi}}\right ) + {\frac{1}{m}}\left (n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mu\right )}{\partial s}} - n{\frac{\partial \left (mu\right )}{\partial \eta}}\right )] + n{\frac{\partial z}{\partial \eta}}A_{M}[{\frac{1}{m}}\left (n{\frac{\partial z}{\partial \eta}}{\frac{1}{H_{z}}}{\frac{\partial \left (mv\right )}{\partial s}} - n{\frac{\partial \left (mv \right )}{\partial \eta}}\right ) - {\frac{1}{n}}\left (m{\frac{\partial z}{\partial \xi}}{\frac{1}{H_{z}}}{\frac{\partial \left (nu \right )}{\partial s}} - m{\frac{\partial \left (nu\right )}{\partial \xi}}\right )] }[/math] (117)
14) Sediment model
a) bedload transport
[math]\displaystyle{ q_{bl} = \Phi \sqrt{\left ( s - 1\right ) g D_{50}^3} \rho_{s} }[/math] (118)
b) Seawater density
[math]\displaystyle{ \rho = \rho_{water} + \Sigma_{m=1}^ \left (N_{sed}\right ) {\frac{C_{m}}{\rho_{s,m}}} \left (\rho_{s,m} - \rho_{water}\right ) }[/math] (119)
c) Bottom Stress Reynolds-averaged Navier-Stokes equations
[math]\displaystyle{ K_{M}{\frac{\partial u}{\partial s}} = \tau_{bx} }[/math] (120)
[math]\displaystyle{ K_{M}{\frac{\partial v}{\partial s}} = \tau_{by} }[/math] (121)
Linear and quadratic drag-coefficient methods
[math]\displaystyle{ \tau_{bx} = \left ( \gamma_{1} + \gamma_{2} \sqrt{u^2 + v^2}\right ) u }[/math] (122)
[math]\displaystyle{ \tau_{by} = \left ( \gamma_{1} + \gamma_{2} \sqrt{u^2 + v^2}\right ) v }[/math] (123)
Classic vertical logarithmic profile
[math]\displaystyle{ |u| = {\frac{u_{*}}{\kappa}} ln \left ({\frac{z}{z_{0}}}\right ) }[/math] (124)
[math]\displaystyle{ |u| = \sqrt{u^2 + v^2} }[/math] (125)
[math]\displaystyle{ u_{*} = \sqrt{|\tau_{bx}| + |\tau_{by}|} }[/math] (126)
Kinematic stresses
[math]\displaystyle{ \tau_{bx} = {\frac{\kappa^2 u \sqrt{u^2 +v^2}}{ln^2 \left ({\frac{z}{z_{0}}}\right )}} }[/math] (127)
[math]\displaystyle{ \tau_{by} = {\frac{\kappa^2 v \sqrt{u^2 + v^2}}{ln^2 \left ({\frac{z}{z_{0}}}\right )}} }[/math] (128)
| Symbol | Description | Unit |
|---|---|---|
| Du,Dv,DT,DS | diffusive terms | - |
| Fu,Fv,FC | forcing terms | - |
| f(x,y) | coriolis parameter | - |
| g | acceleration of gravity | m/s2 |
| h(x,y) | bottom depth | m |
| ν,νθ | molecular viscosity and diffusivity | |
| Km,KC | vertical eddy viscosity and diffusivity | - |
| P | total pressure, approximately -ρo gz | |
| Φ(x,y,z,t) | dynamic pressure, Φ = (P/ρo) | |
| ρo + ρ(x,y,z,t) | total in situ density | |
| S(x,y,z,t) | salinity | |
| t | time | |
| T(x,y,z,t) | potential temperature | |
| u,v,w | the (x,y,z) components of vector velocity | - |
| x,y | horizontal coordinate | |
| z | vertical coordinate | |
| ζ(x,y,t) | surface elevation | |
| QC | surface concentration flux | - |
| τS x, τS y | surface wind stress | |
| τb x, τb y | bottom stress | - |
| S(x,y,δ) | a nonlinear vertical transformation functional (in the vertical transformation equations) | |
| ζ(x,y,t) | time-varying free-surface(in the vertical transformation equations) | |
| h(x,y) | unperturbed water column thickness (in the vertical transformation equations) | |
| σ | fractional vertical stretching coordinate (ranging from -1 <= δ <= 0) | |
| C(σ) | nondimensional, monotonic, vertical stretching function (ranging from -1 <= C(δ) <= 0) | |
| hc | a positive thickness controlling the stretching | |
| Hz | vertical grid thickness, equals to Δz / Δδ | |
| θS, θB | surface and bottom control parameters, have different values for varied Vstretching | |
| γ | a scale factor for all hyperbolic functions, used for Vstretching = 3 | |
| ε, η | curvilinear coordinates | |
| α, β, γ | coefficient used in Time-stepping Schemes: Third-order Adams-Bashforth (AB3) | |
| δξ, δη, δσ | simple centered finite-difference approximations to \partial / \partial ξ, \partial / \partial η, and \partial / \partial σ with the differences taken over the distances Δξ, Δη, and Δσ, respectively | |
| Δz | vertical distance from one ρ point to another | |
| D | total depth, equals to h + ζ | |
| bar{u}, bar{v} | depth-integrated horizontal velocities | |
| AM(ξ, η), KM(ξ, η, s) | spatially varying horizontal and vertical viscosity coefficients, respectively | |
| D50 | grain mean size | |
| ρs | grain density | |
| s | specific density in water, equals to ρs / ρ | |
| τc | critical shear stress | |
| Φ | non-dimensional transport rates | |
| qbl | dimensional bedload transport rates | |
| γ1 | linear drag coefficient for bottom stress calculation | |
| γ2 | quadratic drag coefficient for bottom stress calculation | |
| u* | friction velocity | |
| z0 | a constant bottom roughness length | |
| z | elevation above the bottom (vertical mid-elevation point of bottom cell) | |
| κ | von Karman's constant, equals to 0.41 | |
| zr | reference elevation | |
| ub | representative wave-orbital velocity amplitude | |
| T | wave period | |
| θ | wave propagation direction (clockwise from north) | degree |
| ωs | representative settling velocity | - |
| τb | bed stress associated with mean current above the wave-boundary layer | - |
| τw | wave motions | - |
Notes
1) Momentum equations
In the Boussinesq approximation, density variations are neglected in the momentum equations except in their contribution to the buoyancy force in the vertical momentum equation.
2) Vertical transformation equations
In the vertical transformation equation, if equation (a) is used, the value of hc and other input vertical transformation parameters must be strictly defined. Since the beginning of ROMS, hc = min (hmin,Tcline) where hmin = minval(h) and Tcline is the stretching parameter. It is not possible to build the vertical stretching coordinates when hc > hmin and using transformation (a) because it yields Δz / Δδ < 0. This has been a legacy issue in ROMS for years, and you either consistent with the previous set-up of your application or change the value of Tcline in all input NetCDF files. Therefore, in new applications you need to be sure that Tcline <= hmin when using the transformation (a). However, if transformation (b) is used, hc can be any positive value. It works for both hc <= hmin, or hc > hmin.
Compared to equation (a), Shchepetkin (personal communication) points out that equation (b) has several advantages:
- Regardless of the design of C(σ), it behaves like equally-spaced sigma-coordinates in shallow regions, where h(x,y) << hc. This is advantageous because it avoids excessive resolution and associated CFL limitation is such areas.
- Near-surface refinement behaves more or less like geopotential coordinates in deep regions (level thicknesses, Hz, do not depend or weakly depend on bathymetry), while near-bottom like sigma coordinates (Hz is roughly proportional to depth). This reduces the extreme r-factors near the bottom and reduces pressure gradient errors.
- The true sigma-coordinate system is recovered as hc is close to infinity. This is useful when configuring applications with flat bathymetry and uniform level thickness. Practically, you can achieve this by setting Tcline to 1.0d+16 in ocean.in. This will set hc=1.0*1016. Although not necessary, we also recommend that you set θS=0 and θB=0.
In an undisturbed ocean state, corresponding to zero free-surface, z=S(x,y,σ). Shchepetkin and McWilliams (2005) denotes transformation (a) as an unperturbed coordinate system since all the depths are not affected by the displacements of the free-surface. This ensures that the vertical mass fluxes generated by a purely barotropic motion will vanish at every interface. On the other hand, for transformation equation (b), in an undistrubed ocean state, it yields unperturbed depths as equation (26) and (27), as a consequence, the uppermost grid box retains very little dependency from bathymetry in deep areas, where h_c << h(x,y). For example, if h_c = 250m and h(x,y) changes from 2000 to 6000m, the uppermost grid box changes only by a factor of 1.08 (less than 10%).
Now the transformation equation (b) is the default value for ROMS.
3) Vertical Stretching functions
The vertical stretching functions have the following properties:
- C(σ) is a dimensionless, nonlinear, monotonic function;
- C(σ) is a continuous differentiable function, or a differentiable piecewise function with smooth transition;
- C(σ) must be discritized in terms of the fractional stretched vertical coordinate σ;
- C(σ) must be constrained by -1 <= C(σ) <= 0, that is C(σ)= 0 for σ = 0 at the free-surface, C(σ) = -1 for σ = -1 at the ocean bottom.
The 4 available stretching functions has been shown in equations (28)~(38), equation (28) shows the condition for Vstretching = 1, which has the following features:
- 0 <= θS <= 20, 0 <= θB <= 1, respectively.
- It is infinitely differentiable in σ.
- The larger the value of θS, the more resolution is kept above hc.
- For θB = 0, the resolution all goes to the surface as θS is increased.
- For θB = 1, the resolution goes to both the surface and the bottom equally as θS is increased.
- For θS is not equal to 0, there is a subtle mismatch in the discretization of the model equations, for instance in the horizontal viscosity term. We recommend that you stick with reasonable values of θS, say θS <= 8.
- Some applications turn out to be sensitive to the value of θS used.
Equation (29)~ (32) shows the condition for Vstretching = 2, which is similar in meaning to equation (28), but the hyperbolic function in Csur is cosh instead of sinh and that when σ is close to 0, ΔC/Δσ is also close to 0.
Equations (33)~(36) shows the condition for Vstretching = 3, which is designed for high bottom boundary layer resolution in relatively shallow applications, and the typical values for the control parameters are:
- γ = 3;
- θS = 0.65 for minimal increase of surface resolution;
- θS = 1.0 for significant surface amplification;
- θB = 0.58 for no bottom amplification;
- θB = 1.0 for significant bottom amplification;
- θB = 1.5 for good bottom resolution for gravity flows;
- θB = 3.0 for super-high bottom resolution.
Equations (37), (38) shows the condition for Vstretching = 4, where C(σ) is defined as a continuous, double stretching function. equation (38) is the second stretching of an already stretched transform (37). The resulting stretching function is continuous with respect to θS and θB as their values approach zero. The range of meaningful values for θS and θB are: 0 <= θS <= 10; 0 <= θB <= 4. However, users need to pay attention to extreme r-factor values near the bottom.
4) Coordinate transformation
The terrain-following coordinate transformation is used as from the point of view of the computational model, it is highly convenient to introduce a stretched vertical coordinate system which essentially "flattens out" the variable bottom at z = -h(x,y).
The Curvilinear coordinate transformation is used in many applications of interest (e.g., flow adjacent to a coastal boundary), the fluid may be confined horizontally within an irregular region. In such problems, a horizontal coordinate system which conforms to the irregular lateral boundaries is advantageous. It is often also true in many geophysical problems that the simulated flow fields have regions of enhanced structure (e.g., boundary currents or fronts) which occupy a relatively small fraction of the physical/computational domain. In these problems, added efficiency can be gained by placing more computational resolution in such region.
5) Time step scheme
The Euler time step is the simplest approximation, this method is accurate to first order in Δt, however, it is unconditionally unstable with respect to advection.
The Leapfrog time step is more accurate, but it is unconditionally unstable with respect to diffusion. Also, the even and odd time steps tend to diverge in a computational mode. This computational mode can be damped by taking correction steps. SCRUM's time step on the depth-integrated equations was a leapfrog step with a trapezoidal correction (LF-TR) on every step, which uses a leapfrog step to obtain an initial guess of Φ (t + Δt). This leapfrog-trapezoidal time step is stable with respect to diffusion and it strongly damps the computational mode. However, the right-hand-side terms are computed twice per time step.
The third-order Adams-Bashforth step is used for SCRUM's full 3-D fields. The coefficients for the equations are: α = 23/12, β = -4/3, γ = 5/12. This method requires one time level for the physical fields and three time levels of the right-hand-side information and requires special treatment on startup.
The Forward-Backward equation is a second-order accurate and is stable for longer time steps than many other schemes. It is however unstable for advection terms.
The Forward-Backward Feedback (RK2-FB)equation is one option for solving equations (56) and (57), if we set β=ε=0 in the equations, it becomes a standard second order Runge-Kutta scheme, which is unstable for a non-dissipative system. Adding the β and ε terms adds Forward-Backward feedback to this algorithm, and allows us to improve both its accuracy and stability. The choice of β and ε terms adds Forward-Backward feedback to this algorithm, and allows us to improve both its accuracy and stability. The choice of β = 1/3 and ε = 2/3 leads to a stable third-order scheme with αmax = 2.14093.
LF-TR and LF-AM3 with FB Feedback is another possibility for solving equations (56) and (57). In this equation, keeping γ = 1/12 maintains third-order accuracy. The most accurate choices for β and ε are β = 17/20 and ε = 11/20, which yields a fourth-order scheme with low numerical dispersion and diffusion and αmax = 1.851640.
One drawback of Forward-Backward Feedback and LF-TR and LF-AM3 with FB Feedback is the need to evaluate the right-hand-side terms twice per time step. The original forward-backward scheme manages to achieve αmax = 2 while only evaluating each right-hand-side term once. The Generalized FB with an AB3-AM4 Step is second-order accurate for any set of values for β, γ, and ε. It is third-order accurate if β = 5/12. However, it has a wider stability range for β = 0.281105. Shchepetkin and McWilliams (2009) choose to use this scheme for the barotropic mode in their model with β = 0.281105, γ = 0.0880, and ε = 0.013, to obtain αmax = 1.7802, close to the value of 2 for pure FB, but with better stability properties for the combination of waves, advection, and Coriolis terms.
6) Pressure Gradient Terms in Mode coupling
In the vertical pressure gradient equation, the first term has the disadvantage that after the barotropic time stepping is complete and the new free surface is substituted into the full baroclinic pressure gradient, its vertical integral will no longer be equal to the sum of the new surface slope term and the original coupling term based on the old free surface. It is one form of mode-splitting error which can lead to trouble because the vertically integrated pressure gradient is not in balance with the barotropic mass flux. The second equation is then used.
7) Horizontal and Vertical Advection
For the advection equations, the second-order centered scheme is know to have some unfortunate properties in the presence of strong gradient, such as large over- and under-shoots of tracers, leading to the need for large amounts of horizontal smoothing. ROMS provides alternative advection schemes with better behavior in many situation, but retains this one for comparison purposes.
The Third-order Upwind scheme includes both 1D and 2D. This scheme is known as Uniformly Third-Order Polynomial Interpolation Algorithm (UTOPIA).
8) Boundary Conditions
ROMS comes with a variety of boundary conditions, including open, closed, and periodic.
The Gradient boundary condition is extremely simple and consists of setting the gradient of a field to zero at the edge. The outside value is set to the closest interior value.
The Wall boundary condition is used in ROMS when no other boundary condition is chosen. It is a zero gradient condition for tracers and the surface elevation and zero flow for the normal velocity. For tangential velocities, the wall is treated as either no-slip or free-slip, depending on the value of gamma2 chosen by the user.
The Clamped boundary condition is almost as simple as is setting the boundary value to a known exterior value.
The Flather boundary condition is one option for the normal component of the barotropic velocity to radiate out deviations from exterior values at the speed of the external gravity waves, the exterior values are often used to provide tidal boundary conditions to the barotropic mode. However, there are times when only the tidal elevation is known; The Chapman boundary condition assumes all outgoing signals leave at the shallow-water wave speed of /sqrt(gD). This can be useful when using the Flather condition on the 2D momentum equations.
In realistic domains, open boundary conditions can be extremely difficult to get right. There are situations in which incoming flow and outgoing flow happen along the same boundary or even at different depths at the same horizontal location. Orlanski(1976) proposed a radiation scheme in which a local normal phase velocity is computed and used to radiate things out (if it indeed going out). This works well for a wave propagating normal to the boundary, but has problems when wave approach the boundary at angle. Raymond and Kuo (1984) have modified the scheme to account for progagation in all three directions, the Radiation boundary condition in ROMS only account for two horizontal direction. The radiation approach is appropriate for waves leaving the domain. A check is made to see which way the phrase velocity is headed. If it is entering the domain, a zero gradient condition is applied unless the next step is also specified.
ROMS has an option for providing radiation conditions on outflow and nudging to a known exterior value on inflow. This is implemented as a variation on radiation condition, requiring two timescales: the inflow nudging timescale and the outflow nudging timescale. These timescales are provided in the input to ROMS.
9) Sediment Model
The ROMS model represents sediment using separate cohesive and non-cohesive categories, each with an unlimited number of user-defined size classes. Each class has fixed attributes of grain diameter, density, settling velocity, critical shear stress for erosion, and erodibility constant. These properties are used to help determine the bulk properties of each bed layer.
a) Sediment bed
The sediment bed is represented by three-dimensional arrays with a fixed number of layers beneath each horizontal model cell. Each cell of each layer in the bed is initialized with a thickness, sediment-class distribution, porosity, and age. The mass of each sediment class in each cell can be determined from these values and the grain density. The bed framework also includes two-dimensional arrays that describe the evolving properties of the seabed, including bulk properties of the surface layer (active layer thickness, mean grain diameter, mean density, mean settling velocity, mean critical stress for erosion) and descriptions of the subgrid scale morphology (ripple height and wavelength). These properties are used to estimate bed roughness in the BBL formulations and feed into the bottom stress calculations. The bottom stresses are then used by the sediment routines to determine resuspension and transport, providing a feedback from the sediment dynamics to the hydrodynamics.
The bed layers are modified at each time step to account for erosion and deposition and track stratigraphy. At the beginning of each time step, an active layer thickness za is calculated based on the relation of Harris and Wiberg (1997). The thickness of the top bed layer has a minimum thickness equivalent to za. If the top layer is thicker than za, no action is required. If the top layer is less than za thick, then the top layer thickness is increased by entraining sediment mass from deeper layers until the top layer thickness equals za. If sediment from deeper than the second layer is mixed into the top layer, the bottom layer is split to enforce a constant number of layers and conservation of sediment mass. Each sediment class can be transported by suspended-load and/or bedload (described below). Suspended-load mass is exchanged vertically between the water column and the top bed layer. Mass of each sediment class available for transport is limited to the mass available in the active layer. Bedload mass is exchanged horizontally between the top layers of the bed. Mass of each sediment class available for transport is limited to the mass available in the top layer. Suspended-sediment that is deposited, or bedload that is transported into a computational cell, is added to the top bed layer. If continuous deposition results in a top layer thicker than a user-defined threshold, a new layer is provided to begin accumulation of depositing mass. The bottom two layers are then combined to conserve the number of layers. After erosion and deposition have been calculated, the active-layer thickness is recalculated and bed layers readjusted to accommodate it. This step mixes away any very thin layer (less than the active layer thickness) of newly deposited material. Finally the surficial sediment characteristics, such as D50, ripple geometry, etc., are updated and made available to the bottom stress calculations.
b) Suspended-sediment transport
Sediment suspended in the water column is transported, like other conservative tracers (e.g., temperature and salinity) by solving the advection–diffusion equation (73) with a source/sink term for vertical settling and erosion. The vertical advection algorithm includes a piece-wise parabolic method (Colella and Woodward, 1984) and a weighted essentially non-oscillatory (WENO) scheme (Liu et al., 1994). This method allows the integration bounds of depositional flux to use multiple grid boxes in the vertical direction, so it is not constrained by the CFL criterion. Zero-flux boundary conditions are imposed at the surface and bottom in the vertical diffusion equation. The source or sink term in the advection equation represents the net of upward flux of eroded material and downward settling.
c) Bedload transport
This version of ROMS implements two methods for computing bedload transport: 1) the Meyer-Peter Müeller (1948) formulation for unidirectional flow and 2) the formulae of Soulsby and Damgaard (2005) that accounts for combined effects of currents and waves. As shown in equation (118), the formulae depend on the characteristics of individual sediment classes including median size D50, grain density ρs, specific density in water s =ρs /ρ, and critical shear stress τc. Non-dimensional transport rates Φ are calculated for each sediment class and converted to dimensional bedload transport rates qb.
d) Morphology
The bed model accounts for changes in sea floor elevation resulting from convergence or divergence in sediment fluxes. These morphological changes can have significant influence on flow and transport when they are larger than a few percent of the water depth. The morphological changes are accounted for by equating the bottom boundary condition of the vertical velocity to the rate of change of elevation of the sea floor. This method is completely mass conserving and retains tracer constancy preservation.
A morphological scale factor is also provided to allow an increased rate of morphological change, which can be useful for simulating evolution over long time periods. Strategies for morphological updating are described by Roelvink (2006). In our implementation, bedload fluxes, erosion, and deposition rates are multiplied by a scale factor. A scale factor with a value of one has no effect, and values greater than one accelerate the bed response. For bedload transport, the scale factor is multiplied against the bedload transport rates. For suspended load transport, the scale factor multiplies the exchange of sediment (erosive or depositional flux) at the bed-water interface. The magnitude of sediment concentrations in the water column are not modified – just the exchange rate to and from the bed. For both bedload and suspended load, sediment is limited in availability as described previously, based on the true amount of sediment mass (not multiplied by the scale factor). This morphological scale factor method works well for systems with unlimited sediment in the bed. However, it can generate extra sediment in systems with limited supplies of bed sediment. This occurs when the amount of sediment to be eroded is limited by the amount available and application of the morphological scale factor cannot remove the scaled amount of sediment from the bed. Subsequent deposition does place a scaled amount of sediment on the bed thus creating new mass in the bed. Other approach (Lesser et al 2004) is to limit the amount of sediment fluxed to the water column in these situations. This gives unrealistically low sediment concentrations, but conserves bed sediment.
e) Bottom Stress
Reynolds stresses, production and dissipation of turbulent kinetic energy, and gradients in velocity and suspended-sediment concentrations vary over short vertical distances, especially near the bed, and can be difficult to resolve with the vertical grid spacing used in regional-scale applications. ROMS provides algorithms to parameterize some of these subgrid-scale processes in the water column and in the bottom boundary layer (BBL). Treatment of the BBL is important for the circulation model solution because it determines the stress exerted on the flow by the bottom. Determination of the BBL is even more important for the sediment-transport formulations because bottom stress determines the transport rate for bedload and the resuspension rate for suspended sediment.
ROMS implements either of two methods for representing BBL processes: (1) simple drag-coefficient expressions or (2) more complex formulations that represent the interactions of wave and currents over a moveable bed. The drag-coefficient methods implement formulae for linear bottom friction, quadratic bottom friction, or a logarithmic profile. The other, more complex methods, implement some of the many wave-current BBL models (e.g., Jonsson and Carlsen, 1976; Smith, 1977; Grant and Madsen, 1979; Madsen, 1994; Styles and Glenn, 2000) and couple them with calculations of bottom roughness. ROMS offers three methods that implement slightly different combinations of algorithms for the wave-current interactions and moveable bed roughness. The first method (SG_BBL) is based on the wave-current algorithm and the ripple geometry and moveable bed roughness of Styles and Glenn (2000, 2002). The second method (MB_BBL) uses efficient wave-current BBL computations developed by Soulsby (1995) in combination with sediment and bedform roughness estimates of Grant and Madsen (1982), Nielsen (1986) and Li and Amos (2001). These algorithms and an example of their use on the Southern California continental shelf are described by Blaas et al. (2005). The third method (SSW_BBL) implements either the wave-current BBL model of Madsen (1994) or that of Styles and Glenn (2000) along with moveable bed routines proposed by Wiberg and Harris (1994; Harris and Wiberg, 2001). The differences in approach among these routines are small, but they can produce significantly different results.
The linear and quadratic drag-coefficient methods is shown in equation (121) and (122), the bottom stresses computed from these formulae depend on the elevation of u and v (computed at the vertical mid-elevation of the bottom computational cell). Therefore, in this s-coordinate model, the same drag coefficient will be imposed throughout the domain even though the vertical location of the velocity is different.
The kinematic stresses calculation has the advantage that the velocity and the vertical elevation of that velocity are used in the equation. Because the vertical elevation of the velocity in the bottom computational cell will vary spatially and temporally, the inclusion of the elevation provides a more consistent formulation.
More complex routines are required to simulate BBL processes in the presence of waves and mobile sediment.The short (order 10-s) oscillatory shear of wave-induced motions in a thin (a few cm) wave-boundary layer produces turbulence and generates large instantaneous shear stresses. The turbulence enhances momentum transfer, effectively increasing the coupling between the flow and the bottom and increasing the frictional drag exerted on the mean flow, averaged over many wave periods. The large instantaneous shear stresses often dominate sediment resuspension and enhance bedload transport. Sediment transport can remold the bed into ripples and other bedforms, which present roughness elements to the flow. Bedload transport can also induce drag on the flow, because momentum is transferred to particles as they are removed from the bed and accelerated by the flow. Resuspended sediments can cause sediment-induced stratification and, at high concentrations, change the effective viscosity of the fluid.
Examples
An example run with input parameters, BLD files, as well as a figure / movie of the output
Follow the next steps to include images / movies of simulations:
- Upload file: https://csdms.colorado.edu/wiki/Special:Upload
- Create link to the file on your page: [[Image:<file name>]].
See also: Help:Images or Help:Movies
Developer(s)
- Hernan G. Arango
References
- Blaas, M., Dong, C., Marchesiello, P., McWilliams, J.C., Stolzenbach, K.D., 2005. Sediment-transport modeling on Southern Californian Shelves: a ROMS case study. Con- tinental Shelf Research 27, 832–853.
- Colella, P., and P. L. Woodward, 1984: The piecewise parabolic method (PPM) for gas-dynamical simulations. J. Comput. Phys., 54, 174–201.
- Grant, W.D., Madsen, O.S., 1979. Combined wave and current interaction with a rough bottom. Journal of Geophysical Research 84 (C4), 1797–1808.
- Grant, W.D., Madsen, O.S., 1982. Movable bed roughness in unsteady oscillatory flow. Journal Geophysical Research 87 (C1), 469–481.
- Harris, C.K., Wiberg, P.L., 2001. A two-dimensional, time- dependent model of suspended sediment transport and bed reworking for continental shelves. Computers & Geosciences 27, 675–690.
- Jonsson, I.C., Carlsen, N.A., 1976. Experimental and theoretical investigations in an oscillatory boundary layer. Journal of Hydraulic Research, ASCE 14 (1), 45–60.
- Lesser, G., J.A. Roelvink, J.A.T.M. Van Kester, and G.S. Stelling. 2004. Development and validation of a three-dimensional morphological model. Coastal Engineering, 51: 883–915.
- Li, M.Z., Amos, C.L., 2001. SEDTRANS96: the upgraded and better calibrated sediment-transport model for continental shelves. Computers & Geosciences 27, 619–645.
- Liu, X.-D., Osher S., and Chan T., Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994), pp. 200–212.
- Madsen, O.S., 1994. Spectral wave–current bottom boundary layer flows. In: Coastal Engineering 1994. Proceedings of the 24th International Conference on Coastal Engineering Research Council, Kobe, Japan, pp. 384–398.
- Meyer-Peter, E., and Müller, R. (1948). “Formulas for bed-load transport.” Proc., 2nd Meeting, IAHR, Stockholm, Sweden, 39-64.
- Nielsen, P., 1986. Suspended sediment concentrations under waves. Coastal Engineering 10, 23–31.
- Raymond, W. H. and H. L. Kuo, 1984: A radiation boundary condition for multi-dimensional flows, Quart. J. R. Met. Soc., 110, 535-551.
- Roelvink, J.A., 2006. Coastal morphodynamic evolution techniques. Coastal Engineering, 53: 277– 287.
- Shchepetkin, A. F., and J. C. McWilliams, 2005: The Regional Ocean Modeling System: A split-explicit, free-surface, topography following coordinates ocean model, Ocean Modelling, 9, 347-404.
- Shchepetkin, A. F., and J. C. McWilliams, 2009: Correction and Commentary for "Ocean forecasting in terrain-following coordinates: formulation and skill assessment of the Regional Ocean Modeling System" by Haidvogel et al., J. Comp. Phys., 228, 8985-9000.
- Smith, J.D., 1977. Modeling of sediment transport on continental shelves. In: Goldberg, E.D., McCave, I.N., O’Brien, J.J., Steele, J.H. (Eds.), The Sea, Vol. 6. Wiley-Interscience, New York, pp. 539–577.
- Soulsby, R.L., 1995. Bed shear-stresses due to combined waves and currents. In: Stive, M.J.F. (Ed.), Advances in Coastal Morphodynamics: An Overview of the G8-Coastal Morpho- dynamics Project, Co-Sponsored by the Commission of The European Communities Directorate General XII pp. 4.20–4.23.
- Soulsby, R. L. and Damgaard, J. S. 2005: Bedload sediment transport in coastal waters. Coastal Engineering, 52, 673-689.
- Styles, R., Glenn, S.M., 2000. Modeling stratified wave and current bottom boundary layers on the continental shelf. Journal of Geophysical Research 105 (C10), 24,119–24,139.
- Styles, R., Glenn, S.M., 2002. Modeling bottom roughness in the presence of wave-generated ripples. Journal of Geophysical Research 107 (C8), 24/1–24/15.
- Wiberg, P.L., Harris, C.K., 1994. Ripple geometry in wave- dominated environments. Journal of Geophysical Research 99 (C1), 775–789.
