B.4 Numerical Scheme
The system of equations is solved using a Finite Volume numerical scheme. Further details on Finite Volume methods for hyperbolic systems are provided in LeVeque (2002).
B.4.1 Discrete System
The spatial domain is discretised using contiguous, non-overlapping triangular and quadrilateral cells (or elements).
A cell-centred spatial discretisation is adopted for all NLSWE conserved variables. The discrete form of the equations for cell \(i\), with \(k=1,N_k\) cell-faces is:
\[\begin{equation} \frac{\partial \mathbf{U}_i}{\partial t} = -\frac{1}{A_i} \sum_{k=1}^{N_k} (\mathbf{F}_k \cdot \mathbf{n}_k) L_k + \mathbf{S}_i \tag{B.10} \end{equation}\]
In this discrete equation:
- \(\mathbf{U}_i\) represents the volume-average of the conserved variables in cell \(i\)
- \(A_i\) is the cross-sectional (plan) area of the cell
- \(\mathbf{S}_i\) is the volume-average source term
A first-order midpoint quadrature is used to evaluate the cell boundary flux integral, where:
- \(\mathbf{n}_k\) is the boundary/face unit normal vector for face \(k\)
- \(L_k\) is the corresponding face length
The discrete conserved variable field is assumed to be continuous within a cell but discontinuous at the cell faces.
The above has essentially converted the original system of partial differential equations to one of ordinary differential equations. This allows the solution of the conservation system of equations to be separated into a two-stage algorithm:
- The spatial integration of the discrete flux and source components (right hand side of Equation (B.10))
- The time integration of the discrete system of conservation equations
B.4.2 Spatial Order
The first-order form of the finite volume schemes assumes a piecewise constant \(\mathbf{U}_i\) within each model cell (LeVeque (2002)). Finite-volume schemes with higher-order spatial accuracy can be derived by re-constructing a piecewise continuous \(\mathbf{U}_i\left(x,y,z\right)\) within each model cell. For instance, a second-order spatial scheme can be derived by re-construction of a piecewise linear \(\mathbf{U}_i\left(x,y,z\right)\), while a third-order spatial scheme would require re-construction of a piecewise parabolic \(\mathbf{U}_i\left(x,y,z\right)\). It should be noted that the discrete \(\mathbf{U}_i\) remains discontinuous at cell-faces even for schemes higher than first-order (Hubbard (1999)).
The higher spatial orders can significantly reduce numerical diffusion where the physical system being solved includes large spatial gradients relative to the discrete mesh size. Numerical diffusion can also be reduced through selection of a finer mesh resolution; however the higher spatial order schemes will generally achieve this outcome with less increase in computational overhead.
In general, the solution will only benefit from higher spatial order when the spatial gradients become sufficiently large relative to the mesh size. This can only be determined by testing for improvements in the higher-order solution relative to the first-order solution. If the first-order and high-order solutions are more or less identical for the particular model purpose, then it is generally appropriate to adopt the first-order accuracy. However, if the solutions are significantly different this suggests that first-order numerical diffusion is substantial relative to the physical fluxes that are being resolved in the model. In this case the higher-order solution is probably of a higher quality, though care must be exercised with the higher order solutions to ensure that spurious overshoots at the cell faces are avoided by the reconstruction procedure.
The Total Variation Diminishing (TVD) property (and hence stability) of the higher-order scheme solution is achieved using a choice of gradient limiter schemes. A variety of gradient limiters are available in TUFLOW FV and are listed in order from least to most compressive:
- Horizontal (all after Batten et al. (1996)):
- Limited Central Difference (LCD)
- Maximum Limited Gradient (MLG)
- Vertical (all after Fringer et al. (2005))
- MINMOD
- Maximum Central
- Superbee
The most compressive schemes will maximise the resolution of sharp gradients but may do so at the expense of additional computational overhead. The most compressive gradient limiter schemes also increase the risk of generating spurious overshoots within the solution.
Within TUFLOW FV horizontal and vertical reconstructions are performed separately. A first-order horizontal reconstruction can be combined with a second-order vertical reconstruction, and vice-versa.
B.4.3 Mode Splitting
Efficient integration of the NLSWE is achieved through a mode-splitting scheme, where different components of the governing equations are updated using an appropriate timestep that considers both physical and numerical stability constraints (e.g. Shchepetkin & McWilliams (2005)).
A reduced set of equations comprising all terms other than the barotropic (or free-surface) pressure-gradients is initially partially solved. As part of this solution, an appropriate “internal mode” timestep is calculated that obeys both:
- Courant-Friedrichs-Levy (CFL) limits due to advective current speeds
- Peclet number (Pe) constraints imposed by diffusion terms
Prior to updating (or time-integrating) the solution, an external mode loop is entered, in which a 2D depth-averaged reduction of the 3D NLSWE is solved multiple times, using a timestep that obeys the barotropic Courant-Friedrich-Levy (CFL) constraint imposed by the shallow water wave speed:
\[\begin{equation} \bar{u} \pm \sqrt{gh} \tag{B.11} \end{equation}\]
where \(\bar{u}\) is the depth-averaged current speed. The external mode loop is repeated until the cumulative timestep is approximately equal to the internal mode timestep.
The depth-averaged inviscid fluxes from the external mode solution are then used to correct the internal mode inviscid fluxes, ensuring that they represent the total inviscid flux for the 3D solution. The corrected fluxes are then used to update the full 3D solution.
A stability constraint imposed by the baroclinic internal wave speed is not explicitly calculated and may not always be met by the mode splitting scheme. If oscillations in the pycnocline cause numerical instabilities, this can be addressed by reducing the upper-limiting timestep.
Viscous fluxes and both inviscid and viscous scalar transport fluxes are calculated only for the internal mode (outer) loop.
Mode splitting can be disabled for 2D simulations, and this configuration can be more computationally efficient for fast, shallow flow scenarios where the internal mode and external mode timesteps are similarly restrictive. Three dimensional simulations are only supported with mode splitting enabled.
B.4.4 Flux Terms
A key step in the Finite-Volume numerical scheme is the calculation of numerical fluxes across cell boundaries:
- Inviscid fluxes \((\mathbf{F}_x^I, \mathbf{F}_y^I, \mathbf{F}_z^I)\) represent the directly resolved flux of mass and momentum between adjacent cells
- Viscous fluxes \((\mathbf{F}_x^V, \mathbf{F}_y^V, \mathbf{F}_z^V)\) represent the mixing of mass and momentum that is not directly resolved as advection within the numerical model
B.4.4.1 Viscous Fluxes
Viscous flux terms are calculated using the traditional gradient-diffusion model, with a variety of options available for the calculation of eddy-viscosity and scalar diffusivity.
B.4.4.1.1 Horizontal Viscous Fluxes
The horizontal viscous fluxes \((\mathbf{F}_x^V, \mathbf{F}_y^V)\) are calculated according to Equations (B.5) (momentum) and (B.8) (scalars). The horizontal eddy-viscosity can be specified directly using a temporally-constant value or calculated using the Smagorinsky (1963) or J. Wu (1982) formulation.
The horizontal scalar-diffusivity tensor can also be specified directly (as an isotropic constant value), calculated using the Smagorinsky (1963) formulation or computed using the Elder formulation (Falconer et al. (2005)). The Elder model calculates a non-isotropic diffusivity tensor that accounts for velocity dispersion processes not resolved in 2D depth-averaged models:
\[\begin{equation} D_{xx} = \frac{(D_l u^2 + D_t v^2) h}{u_*}, \quad D_{yy} = \frac{(D_l v^2 + D_t u^2) h}{u_*} \tag{B.12} \end{equation}\]
\[\begin{equation} D_{xy} = D_{yx} = \frac{(D_l - D_t) uvh}{u_*} \tag{B.13} \end{equation}\]
where:
- \(D_l\) and \(D_t\) are the Elder coefficients in the directions lateral to and transverse to the local currents, respectively,
- \(u_* = \sqrt{\frac{|\tau_b|}{\rho}}\) is the friction velocity.
The observed range of values for \(D_l\) and \(D_t\) derived from measurements is discussed in Fischer et al. (1979). In 3D model simulations, the Smagorinsky formulation is generally more applicable. The Elder formulation is recommended only for 2D simulations.
B.4.4.1.2 Vertical Viscous Fluxes
The vertical viscous fluxes \(\mathbf{F}_z^V\) are calculated according to Equations (B.5) (momentum) and (B.8) (scalars). An unconditionally stable semi-implicit scheme is adopted in the discretization of \(\mathbf{F}_z^V\) to prevent timestep restrictions.
The vertical eddy viscosity, \(\nu_t\), can be directly specified as constant, computed using TUFLOW FV’s turbulence model, or calculated from the simple parametric model formulation including the Munk & Anderson (1948) stability function:
\[\begin{equation} \nu_{t0}=\kappa u_\ast z\left(c_1-c_2\frac{z}{h}\right) \tag{B.14} \end{equation}\]
\[\begin{equation} \nu_t=\nu_{t0} f(Ri) \tag{B.15} \end{equation}\]
where Ri is the gradient Richardson number defined as:
\[\begin{equation} Ri=\frac{N^2}{\left(\frac{\partial u}{\partial z}\right)^2} \tag{B.16} \end{equation}\]
and N is the Brunt-Vaisala frequency (or buoyancy frequency):
\[\begin{equation} N=\sqrt{-\frac{g}{\rho}\frac{\partial\rho}{\partial z}} \tag{B.17} \end{equation}\]
Scalar diffusivities are set at 0.74 * \(\nu_t'\).
B.4.4.2 Inviscid Fluxes
The inviscid fluxes \(\mathbf{F}_x^I, \mathbf{F}_y^I, \mathbf{F}_z^I\) represent the directly resolved flux of mass and momentum between adjacent cells. Inviscid fluxes are computed at each cell face based on the conserved variable state immediately on either side of the face.
For a first-order spatial scheme, these values are equivalent to the adjacent-cell averages. In higher-order schemes, the conserved variable state at the cell faces is reconstructed from the cell-averaged values.
B.4.4.2.1 Internal Mode
The internal mode inviscid flux calculations solve the full 3D NLSWE, excluding terms related to the free-surface pressure gradient. A centered scheme is used for the internal mode mass flux, while an upwind scheme is used for momentum flux terms.
Flux-like source terms originating from bed slope and baroclinic pressure gradients are included in the cell-face flux calculation, rather than being part of the volume-integrated source term in Equation (B.10) (see Section B.2).
The internal mode timestep is determined by a combination of internal advection CFL constraints and viscous flux Péclet constraints. A stable internal model timestep is selected before entering the external mode.
B.4.4.2.2 External Mode
The external mode inviscid flux calculations solve the 2D depth-averaged NLSWE. In 3D simulations, the external mode is initiated by computing depth-averages of the 3D conserved variable fields. Viscous fluxes and baroclinic pressure gradients are also depth-integrated at the start of the external mode loop.
The 2D depth-averaged NLSWE are solved using Roe’s approximate Riemann solver (Roe (1981)). Flux-like source terms, such as bed slope and depth-averaged baroclinic pressure gradients, are included in the cell-face flux calculation rather than being incorporated in the volume-integrated source term in Equation (B.10) (see Section B.2).
The external mode timestep is dictated by the surface gravity wave CFL constraint and depth-averaged viscous flux Péclet constraints (Murillo et al. (2005)). A stable timestep is selected for each external mode sub-timestep. Within the external mode loop, multiple sub-timesteps are executed before returning to the outer internal mode loop.
B.4.4.2.3 Flux Correction
The internal mode inviscid fluxes are corrected using the depth-averaged external mode fluxes, which have been integrated in time through the external mode loop.
When mode splitting is disabled, the full NLSWE (including free-surface pressure gradients) is solved directly using the same flux scheme as the External Mode calculation. This option is currently only available for 2D simulations and can be more computationally efficient than mode splitting for fast shallow flow conditions, for example urban floodplain flows.
B.4.4.3 Total Flux
The total flux vector is the sum of the corrected inviscid and viscous flux components.
B.4.4.4 Flux Spatial Integration
The first term on the right hand side of Equation (B.10) requires calculating the boundary integral of the total flux vector normal component, which is approximated using a midpoint quadrature rule.
For momentum flux terms, values are converted to momentum flux differences before integration. In spherical coordinates, the momentum and flux vectors are adjusted from face-centered to cell-centered using a parallel transport transformation. This accounts for rotational effects in the spherical coordinate system (Rossmanith et al. (2004)).
B.4.5 Computational Timestep
B.4.5.1 Stability Criterion
The nonlinear shallow water equations (NLSWE) are solved using a variable timestep scheme. As a simulation progresses, the solver uses the model state to compute a stable timestep by ensuring the following three stability criteria are met (Murillo et al. (2005)).
- Wave celerity criterion - Equation (B.18)
- Advective CFL - Equation (B.19)
- Peclet criterion for turbulent mixing stability - Equation (B.20)
\[\begin{equation} \frac{\left| \mathbf{u}\cdot\mathbf{n} \pm \sqrt{g\,h} \right|\, \Delta t}{L^*} \;\leq\; 1 \tag{B.18} \end{equation}\]
\[\begin{equation} \frac{\left| \mathbf{u}\cdot\mathbf{n} \right|\, \Delta t}{L^*} \;\leq\; 1 \tag{B.19} \end{equation}\]
\[\begin{equation} \frac{|\mathbf{D} \cdot \mathbf{n}| \, \Delta t}{(L^*)^2} \;\leq\; 1 \tag{B.20} \end{equation}\]
where:
- \(\mathbf{u}\) = depth averaged velocity vector (m/s)
- \(\mathbf{n}\) = outward unit normal vector (face normal)
- \(g\) = gravitational acceleration (m/s\(^2\))
- \(h\) = water depth (m)
- \(\Delta t\) = timestep (s)
- \(L^*\) = characteristic length scale (e.g., face length or cell size) (m)
- \(\mathbf{D}\) = calculated eddy viscosity (m\(^2\)/s)
The characteristic length scale is defined as:
\[\begin{equation} L^* = \frac{\min(A_i, A_j)}{L_k} \tag{B.21} \end{equation}\]
where:
- \(A_i, A_j\) = plan areas of the two adjacent cells (m\(^2\))
- \(L_k\) = length of the face separating the two cells (m)
B.4.5.2 Model Timestep 2D
In 2D, the variable timestep scheme ensures that the model satisfies each criterion across all cells while maintaining the largest possible stable timestep. The maximum allowable timestep is computed by rearranging Equations (B.18), (B.19), and (B.20) with respect to timestep.
\[\begin{equation} \Delta t_c \;\leq\; \frac{L^*}{\left| \mathbf{u}\cdot\mathbf{n} \pm \sqrt{g\,h} \right|} \tag{B.22} \end{equation}\]
\[\begin{equation} \Delta t_u \;\leq\; \frac{L^*}{\left| \mathbf{u}\cdot\mathbf{n} \right|} \tag{B.23} \end{equation}\]
\[\begin{equation} \Delta t_d \;\leq\; \frac{(L^*)^2}{|\mathbf{D} \cdot \mathbf{n}|} \tag{B.24} \end{equation}\]
where:
- \(\Delta t_c\) = timestep limited by wave celerity (s)
- \(\Delta t_u\) = timestep limited by advective velocity (s)
- \(\Delta t_d\) = timestep limited by turbulence diffusion (eddy viscosity) (s)
By default, TUFLOW FV executes the solution using the computational mode splitting approach described in Appendix B.4.3. This results in two model timesteps: an external timestep, which controls the update of the free surface and governs wave celerity stability, and an internal timestep, which constrains advection and diffusion within the 2D domain. This separation allows the solver to maintain stability for both surface wave propagation and internal flow processes while maintaining computational efficiency.
The maximum allowable external mode timestep is calculated using Equation (B.25).
\[\begin{equation} \Delta t_{\max\_external} \;=\; \min\!\left( \frac{1}{\tfrac{1}{\Delta t_c} + \tfrac{1}{\Delta t_d}},\; \Delta t_u \right) \tag{B.25} \end{equation}\]
The maximum allowable internal mode timestep is calculated using Equation (B.26).
\[\begin{equation} \Delta t_{\max\_internal} \;=\; \frac{1}{\tfrac{1}{\Delta t_u} + \tfrac{1}{\Delta t_d}} \tag{B.26} \end{equation}\]
By default, a single global user defined
\[\begin{equation} \Delta t_{\text{external}} \;=\; \mathrm{CFL}_{user} \cdot \Delta t_{\max\_external} \tag{B.27} \end{equation}\]
\[\begin{equation} \Delta t_{\text{internal}} \;=\; \mathrm{CFL}_{user} \cdot \Delta t_{\max\_internal} \tag{B.28} \end{equation}\]
Alternatively, the user may override the global CFL and specify separate scaling factors using
\[\begin{equation} \Delta t_{\text{external}} \;=\; \mathrm{CFL}_{user\_external} \cdot \Delta t_{\max\_external} \tag{B.29} \end{equation}\]
\[\begin{equation} \Delta t_{\text{internal}} \;=\; \mathrm{CFL}_{user\_internal} \cdot \Delta t_{\max\_internal} \tag{B.30} \end{equation}\]
The resultant internal and external mode timesteps can be further constrained to user defined minimum and maximum allowable timestep limits via the
An example of the calculated internal and external timesteps is shown in the TUFLOW FV Console Window (Figure B.1). The internal timestep is highlighted by the yellow box, the external by the red.
Figure B.1: Internal And External Model Timesteps
B.4.5.3 Model Timestep 3D
In 3D, the external mode timestep is determined using the same formulation as the 2D HD simulation class Equation (B.25). In evaluating the 3D external model timestep (B.31), the full water depth together with depth-averaged velocity and eddy viscosity are used to compute \(\Delta t_c\), \(\Delta t_u\), and \(\Delta t_d\).
\[\begin{equation} \Delta t_{\max\_external} \;=\; \min\!\left( \frac{1}{\tfrac{1}{\Delta t_c(2D)} + \tfrac{1}{\Delta t_d(2D)}},\; \Delta t_u(2D) \right) \tag{B.31} \end{equation}\]
Internal mode timestep calculations are completed on each 3D cell horizontal face (Equation (B.32)) and each vertical cell face (Equation (B.33)). This is completed across the model domain with the final 3D internal model timestep defined as the minimum of horizontal and vertical face timesteps (Equation (B.34)).
\[\begin{equation} \Delta t_{\max\_internal-h} \;=\; \min\left( \frac{1}{\tfrac{1}{\Delta t_u(i)} + \tfrac{1}{\Delta t_d(i)}} \right) \tag{B.32} \end{equation}\]
\[\begin{equation} \Delta t_{\max\_internal-v} \;=\; \min\left( \frac{\min\!\big(\Delta z^{+}_k,\;\Delta z^{-}_k\big)}{|w_k|} \right) \tag{B.33} \end{equation}\]
\[\begin{equation} \Delta t_{\max\_internal} \;=\; \min\!\left( \Delta t_{\max\_internal-h},\; \Delta t_{\max\_internal-v} \right) \tag{B.34} \end{equation}\]
where:
\(\Delta t_{\max\_internal-h}\) maximum allowable horizontal internal timestep (s)
\(\Delta t_{\max\_internal-v}\) maximum allowable vertical internal timestep (s)
\(\Delta t_u(i),\, \Delta t_d(i)\) are maximum allowable timestep limits at each 3D horizontal cell face \(i\) (s)
\(\Delta t_u(k)\) is the maximum allowable advective timestep limit at each 3D vertical cell face \(k\) (s)
\(\Delta z^{+}_k\) = thickness of the cell above face \(k\) (m)
\(\Delta z^{-}_k\) = thickness of the cell below face \(k\) (m)
\(w_k\) = vertical velocity at face \(k\) (m/s)
The resulting \(\Delta t_{\max\_external}\) and \(\Delta t_{\max\_internal}\) are applied in Equations (B.27) and (B.28) when using the global
B.4.5.4 General Usage Guidance
This section provides a suggested workflow for selecting initial upper and lower model
- Estimate the external mode timestep
Use Equation (B.22) with indicative values of flow velocity \(u\) (m/s), water depth \(h\) (m), and \(L^*\) (m). The shallow water wave celerity will generally control the timestep.
Example:
- \(u = 0.5\) m/s
- \(h = 10\) m
- \(L^* = 10\) m (from a 10 m x 10 m cell)
- \(g = 10\) m/s2
Therefore:
- Celerity: \(c = 10\) m/s
- Timestep: \(\Delta t_c = 0.95\) s
Estimate the upper timestep limit
Take 10 x the external mode timestep. For the above example: 9.5 s.Set timestep limits
Use the estimates from Steps 1 and 2 as the lower and upper arguments of theTimestep Limits command.Run the model
Review external and internal timesteps in the TUFLOW FV log file, console or CFL diagnostic files.Check timestep behaviour
- If timesteps remain within the initial estimates, proceed to Step 6.
- If the external timestep is equal to the minimum estimate, reduce the lower timestep limit and repeat Steps 2-5.
- If timesteps remain within the initial estimates, proceed to Step 6.
Optimse timestep
- Can the user defined lower timestep limit be safely increased after review of Step 5? For example, if the user specified lower timestep is 0.1 s, but the model is consistently running with a lower timestep of 0.5 s, raise the lower timestep to 0.4 s and the corresponding upper timestep to 4.0 s. Then repeat Step 4-5 to check the model behaviour. If not, proceed to Step 7.
- Review mesh performance
Identify cells controlling the minimum and mean timesteps. Adjust the mesh at these locations if possible. Diagnostic CFL outputs that assist in doing so are described at https://fvwiki.tuflow.com/A_Model_Runs_Slow.
B.4.6 Time Integration
This section describes time integration methods. Some material from Appendix B.4.5.2 and B.4.5.3 is necessarily repreated for completeness.
Both internal mode and external mode temporal integration is performed using an explicit Euler scheme. To maintain numerical stability, the time step must satisfy the Courant-Friedrich-Levy (CFL) criterion for wave propagation and advection, and the Péclet criterion for diffusion (Murillo et al. (2005)).
The external mode CFL criterion is:
\[\begin{equation} \left| u \cdot n \pm \sqrt{gh} \right| \frac{\Delta t}{L^*} \leq 1 \tag{B.35} \end{equation}\]
where:
- \(\Delta t\) is the integration timestep,
- \(L^*\) is a cell-size dependent length scale.
The internal mode CFL criterion is:
\[\begin{equation} \max (|u \cdot n|, c_{\text{baro}} ) \frac{\Delta t}{L^*} \leq 1 \tag{B.36} \end{equation}\]
where:
- \(c_{\text{baro}}\) is the baroclinic (internal) wave speed.
The Péclet criterion for diffusive terms is:
\[\begin{equation} \frac{|D \cdot n| \Delta t}{(L^*)^2} \leq 1 \tag{B.37} \end{equation}\]
The cell-size dependent length scale \(L^*\) is computed for each cell-face as:
\[\begin{equation} L^* = \frac{\min(A_i, A_j)}{L_k} \tag{B.38} \end{equation}\]
where:
- \(A_i, A_j\) are the areas of the adjacent cells,
- \(L_k\) is the face length.
A variable time step scheme is implemented to ensure that the CFL and Peclet criterion are satisfied at all points in the model with the largest possible time step. Outputs providing information relating to performance of the model with respect to the CFL criterion are provided to enable informed refinement of the model mesh.
In stratified flows the baroclinic wave speed may impose a constraint on the stable internal mode timestep. However, the internal mode timestep is not automatically adjusted to satisfy the baroclinic wave speed limit. Additionally, the mode splitting scheme stability may benefit from limiting the ratio between the internal and external mode timestep to around 10 or less.
Maximum and minimum timestep limits are specified by the user. The maximum limit should be used to limit the upper internal mode timestep. The minimum limit should be used to restrict the external mode timestep in the event of a model instability, as it is preferable to have the model violate the prescribed stability bounds than have the timestep decrease towards zero.
B.4.7 Wetting and Drying
In shallow regions, the momentum terms are dropped to maintain stability as the NLSWE approach the zero-depth singularity. Mass conservation is maintained both locally and globally to the limit of numerical precision across the entire numerical domain, including wetting and drying fronts. A conservative mass redistribution scheme ensures that negative depths are avoided at numerically challenging wetting and drying fronts, without requiring timestep adjustments (Brufau et al. (2004); Murillo et al. (2006)).
Regions of the model domain that become dry are automatically excluded from computations to improve numerical efficiency.
B.4.8 Source Terms
B.4.8.1 Bed Slope
Bed slope integral source terms are calculated using a face centred upwind flux correction within the internal and external mode numerical flux solvers.
\[\begin{equation} \int_{\Omega} -gh \nabla z_b \, d\Omega \approx \sum_{k=1}^{N_k} \beta^* (\Delta z_b )_k L_k \tag{B.39} \end{equation}\]
That is, the cell-face bed elevation jump \(\Delta z_b\) becomes a correction term \(\beta^* (\Delta z_b )\) to the cell-face numerical flux terms. This numerical approach provides consistent upwinding between flux and bed-slope source terms. This is essential to obtaining the required numerical balance between these terms, at for instance the quiescent state equilibrium.
Further details are provided in Hubbard & Garcia-Navarro (2000) and Murillo et al. (2006).
B.4.8.2 Coriolis Force
Coriolis forces due to Earth’s rotation are calculated as cell-averaged source terms in the momentum equation. The Coriolis coefficient \(f_c\) is calculated from:
\[\begin{equation} f_c = 2\Omega_r \sin\phi \tag{B.40} \end{equation}\]
where
- \(\Omega_r\) is the angular frequency of Earth’s rotation (rad/s) and
- \(\phi\) is the geographic latitude (radians).
Coriolis calculations require the
B.4.8.3 Wind Stress
The cell-averaged surface stress vector due to wind is calculated from:
\[\begin{equation} \tau_{sw} = \rho_a c_{dw} u_w |u_w| \tag{B.41} \end{equation}\]
where the wind drag coefficient is calculated using the empirical formula of J. Wu (1980) and J. Wu (1982):
\[\begin{equation} c_{dw} = \begin{cases} c_a, & w_{10} < w_a \\ c_a + \frac{(c_b - c_a)}{(w_b - w_a)} (w_{10} - w_a), & w_a \leq w_{10} < w_b \\ c_b, & w_{10} \geq w_b \end{cases} \tag{B.42} \end{equation}\]
with default parameters \((w_a, c_a, w_b, c_b) = (0.0 m/s, 0.8 \times 10^{-3}, 50.0 m/s, 4.05 \times 10^{-3})\).
B.4.8.4 Bed Friction
Bed friction momentum sink terms are calculated using a quadratic drag law:
\[\begin{equation} \tau_{bf} = \rho c_{db} u |u| \tag{B.43} \end{equation}\]
where the bottom drag coefficient can be calculated using a roughness-length relationship:
\[\begin{equation} c_{db} = \left( \frac{\kappa}{\ln \left( \frac{30 z'}{k_s} \right)} \right)^2 \tag{B.44} \end{equation}\]
The above relationship assumes a rough-turbulent logarithmic velocity profile in the lowest model layer, where \(\kappa\) is von Karman’s constant, \(k_s\) is the effective bed roughness length (equivalent Nikuradse roughness), and \(z'\) is the height of the bottom cell centroid above the bed.
Instead of specifying \(k_s\), Manning’s \(n\) roughness can be specified and is internally converted into an equivalent roughness length:
\[\begin{equation} k_s = 11 h \exp \left( -\frac{\kappa h^{1/6}}{\sqrt{g} n} \right) \tag{B.45} \end{equation}\]
Bed roughness values (\(k_s\) or Manning’s \(n\)) may be specified globally or be spatially varying via material specification.
The above bed friction formulations are applicable in both 2D (depth-averaged) and 3D configurations. In 2D situations, the Manning’s \(n\) formulation is equivalent to the following equation for the friction slope (Chow (1959)):
\[\begin{equation} S_f = \frac{\tau_{bf}}{\rho gh} = \frac{n^2 \bar{u} |\bar{u}|}{h^{4/3}} \tag{B.46} \end{equation}\]
When comparing 2D and 3D simulations using the same bed roughness parameters, calculated bed friction energy losses are typically not exactly equivalent except in the simplest fully-developed, uniform flow scenarios. This is because 2D models assume a logarithmic velocity profile extending over full depth, whereas 3D simulations resolve the vertical velocity profile, which may be non-logarithmic in more complex flow situations.
Integrated bed friction source terms are calculated using a semi-implicit discretisation in order to maintain unconditional numerical stability of these terms in high-velocity or shallow flows (Brufau et al. (2004)). Coupling of the internal and external modes is achieved by applying the internal mode (3D) bed friction as an explicit momentum sink/source term during external mode (2D) loop calculations.
B.4.8.5 Mean Sea Level and Baroclinic Pressure Gradients
Mean Sea Level Pressure and Baroclinic pressure gradient source terms are calculated as face-centred flux correction terms within the internal and external mode numerical flux solvers. That is, the pressure gradient terms are treated in a similar manner to the bed slope source terms as described in Section 4.7.1, i.e.,
\[\begin{equation} \int_{\Omega} -(\nabla P) \, d\Omega \cong \sum_{k=1}^{N_k} \eta^* (\Delta P)_k L_k \tag{B.47} \end{equation}\]
where
- \(\nabla P\) is the gradient of the combined atmospheric and baroclinic pressure fields and
- \(\eta^*\) is a face-centred flux correction due to the cell-face pressure jump \(\Delta P\).
This is analogous to converting the cell-volume source term integral into a cell-boundary source term integral using Gauss’s theorem.
B.4.8.6 Wave Radiation Stress
Wave fields are applied as spatially and temporally varying datasets on a 2D rectilinear/curvilinear grid.
Wave radiation stress gradients are calculated as cell-centred source terms:
\[\begin{equation} \int_{\Omega} \left( \frac{\partial s_{xx}}{\partial x} + \frac{\partial s_{xy}}{\partial y} \right) d\Omega \tag{B.48} \end{equation}\]
or as face-centred momentum flux source terms. The wave radiation stress gradients are distributed uniformly throughout the water column.
B.4.8.7 Scalar Decay
Tracer constituents can be specified with a linear scalar decay property (Equation (B.9)), where \(K_d\) is the constant linear decay coefficient. Scalar decay is discretised explicitly as a cell-centred integral source term. Numerical stability of this term is not guaranteed for large \(K_d\) or for large model timesteps.
B.4.8.8 Scalar Settling
Tracer and sediment constituents can be specified with a settling velocity \(w_s\) (Equation (B.9)).
Within the water column, the settling velocity contributes an additional (vertically downward) inviscid flux component. At the bed, the settling velocity contributes a sink from the water column and a source into the bed. In the case of sediment fractions and particulate water quality constituents, the mass transferred to the bed is subsequently tracked within the TUFLOW FV Sediment Transport and Water Quality Modules, respectively. The passive tracer constituent mass exiting the water column is no longer tracked.