9.1 2D Domains

9.1.1 2D Depth Averaged Equation of Motion

The equation of motion of a passive tracer \(\phi\) (as a concentration), in 2D depth averaged conservative formulation is:

\[\begin{equation} \frac{\partial(h \phi)}{\partial t} + \frac{\partial(h u \phi)}{\partial x} + \frac{\partial(h v \phi)}{\partial y} - \frac{\partial}{\partial x} \left( h D \frac{\partial \phi}{\partial x} \right) - \frac{\partial}{\partial y} \left( h D \frac{\partial \phi}{\partial y} \right) = h S_\phi \tag{9.1} \end{equation}\]

Where:

• \(\phi\) = Tracer concentration as mass (or mols) per unit volume
• \(h\) = Water depth
• \(u\) and \(v\) = Depth averaged velocity components in the x and y directions
• \(x\) and \(y\) = 2D spatial coordinates
• \(t\) = Time
• \(D\) = Isotropic dispersion plus turbulent diffusion coefficients values in the x and y directions
• \(S_\phi\) = Tracer source as mass (or mols) per unit volume per unit time

The first term in the equation can be expanded using the chain rule and then substituting the continuity equation to get:

\[\begin{equation} \frac{\partial(h\phi)}{\partial t} = h \frac{\partial \phi}{\partial t} + \phi \frac{\partial h}{\partial t} = h \frac{\partial \phi}{\partial t} - \phi \left( \frac{\partial (h u)}{\partial x} + \frac{\partial (h v)}{\partial y} \right) \tag{9.2} \end{equation}\]

Combining Equations (9.1) and (9.2) we get:

\[\begin{equation} \frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} + v \frac{\partial \phi}{\partial y} - \frac{1}{h} \left[ \frac{\partial}{\partial x} \left( h D \frac{\partial \phi}{\partial x} \right) + \frac{\partial}{\partial y} \left( h D \frac{\partial \phi}{\partial y} \right) \right] = S_\phi \tag{9.3} \end{equation}\]

9.1.2 Solution Method (TUFLOW Classic)

The 2D advection dispersion algorithm in TUFLOW Classic is based on the ULTIMATE QUICKEST method of Leonard (1991), Leonard & Niknafs (1991) and Leonard et al. (1993). The equation of motion is based on Equation (9.3), but reformulated to make use of the offset grid and the TVD interpolation scheme (i.e. it is better to first compute \(\phi\) at the cell faces and then compute \(\partial(u\phi)/\partial x\) rather than \(\partial\phi/\partial x\)):

\[\begin{equation} \frac{\partial \phi}{\partial t} + \frac{\partial (u \phi)}{\partial x} + \frac{\partial (v \phi)}{\partial y} - \frac{1}{h} \left[ \frac{\partial}{\partial x} \left(D_{x}\frac{\partial\phi}{\partial x}\right) + \frac{\partial}{\partial y} \left(D_{y}\frac{\partial\phi}{\partial y}\right) \right] = S_\phi + \phi \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) \tag{9.4} \end{equation}\]

TUFLOW AD applies the dispersion formulation described by Falconer et al. (2005). This formulation computes dispersion in the X and Y directions (\(D_x\) and \(D_y\) respectively, to suit the Cartesian computational grid) from user specification of longitudinal and transverse dispersion coefficients \(K_L\) and \(K_T\), respectively. Specifically, \(D_x\) and \(D_y\) are computed dynamically at each grid cell and timestep as follows:

X Direction dispersion plus turbulent diffusion:

\[\begin{equation} D_x = max \left( \frac{(K_L u^2 + K_T v^2) h \sqrt g}{|U| C}, D_{w} \right) \tag{9.5} \end{equation}\]

Y Direction dispersion plus turbulent diffusion:

\[\begin{equation} D_y = max \left( \frac{(K_T u^2 + K_L v^2) h \sqrt g}{|U| C}, D_{w} \right) \tag{9.6} \end{equation}\]

Where:

• \(K_L\) = User specified longitudinal dispersion coefficient
• \(K_T\) = User specified transverse diffusion coefficient
• \(g\) = Gravitational acceleration
• \(|U|\) = Velocity magnitude
• \(C\) = Chezy bed friction coefficient \(C = \frac{h^\frac{1}{6}}{n}\)
• \(D_w\) = User specified lower bound dispersion coefficient

The value of \(D_w\) can be specified as constant or spatially varied as required.

The source terms also includes sink terms such as settling (for particulate species) and decay. The scheme also includes representation of mixing due to sub-grid-scale turbulence and vertical shear via the dispersion formulation provided in Falconer et al. (2005).

The TUFLOW AD computational procedure used is an explicit scheme based on Leonard (1991). This contrasts to the TUFLOW Classic engine, which employs an implicit scheme. As such, TUFLOW AD is generally subject to stricter stability constraints than the Classic hydrodynamic engine. As such, the TUFLOW AD calculation takes the form of three steps within each timestep.

The first step involves calculation of the Courant-Friedrichs-Lewy (CFL) condition at all wet cells, where the CFL in 1 dimension is:

CFL Condition:

\[\begin{equation} CFL = \frac{u \bullet \Delta t}{\Delta x} \tag{9.7} \end{equation}\]

Where:

• \(u\) (or \(v\)) = Fluid velocity
  
• \(\Delta x\) = Grid scale
  
• \(\Delta t\) = Timestep

This condition is typically required to be less than 1.0 (additive for both X and Y directions) and has a broad physical interpretation requiring that the distance fluid is advected in one timestep (\(u∆t\)) is less than one grid cell (\(∆x\)).

The second step is the computation of a similar condition for the diffusive lengthscale (related to the Peclet number) that ensures that dispersion also does not cause instability at any grid cell. The CFL and dispersion dimensionless numbers are then added and the maximum sum at any given location within a timestep is used to compute the number of sub-stepping iterations required by TUFLOW AD to remain stable within one TUFLOW timestep.

The third step within each TUFLOW timestep is to execute the advection dispersion calculations for the required number of iterations, with a modified (smaller) \(∆t\).

The original ULTIMATE QUICKEST solution method has been enhanced and improved as applied in TUFLOW AD. For example, TUFLOW AD also employs adaptive computational stencil expansion where it identifies sharp constituent concentration gradients (Leonard & Niknafs, 1991). Where possible (i.e. away from boundaries and dry cells) and required, the ULTIMATE QUICKEST stencil is expanded from the standard third order scheme to a ninth order scheme, only on principle computational axes. Cross terms greater than third order are not included. If insufficient wet cells exist to switch to ninth order, then seventh and fifth order schemes are progressively tested (with commensurately decreasing stencils) until all required wet cells are located.

Application of the ULTIMATE limiter (Leonard, 1991) has been found to induce steady flow anisotropy when extended to multi-dimensional problems and the numerical cross terms associated with additional dimensions are included in calculations. Wu & Falconer (2000) developed a modification to the ULTIMATE limiter that reduces this anisotropy, and this has been applied within the TUFLOW AD computational engine.

9.1.3 Solution Method (TUFLOW HPC)

The TUFLOW HPC solver tracks areal density, \(b=h\phi\), of the passive tracer as the primary prognostic variable. With this change, equations (9.1) and (9.2) combine to become:

\[\begin{equation} \frac{\partial b}{\partial t} + \frac{\partial (ub)}{\partial x} + \frac{\partial (vb)}{\partial y} - \frac{\partial \phi_{bu}}{\partial x} - \frac{\partial \phi_{bv}}{\partial y} = S_b \tag{9.8} \end{equation}\]

Where:

• \(\phi_{bu}\) and \(\phi_{by}\) are the combined dispersive and diffusive unit fluxes of tracer in the x and y directions respectively
• \(S_b\) = source terms (areal density rate)

The diffusive fluxes of tracer across cell faces are comprised of two components. The first is dispersion (induced by shear in the vertical velocity profile) which acts in the direction of the flow, and the second is turbulent diffusion which acts in both longitudinal and transverse directions. Both calculations use the same formulation, but with different coefficients (Falconer et al., 2005):

\[\begin{equation} \left( \substack{ \phi_{bL} \\ \phi_{bT}} \right) = -\frac{|U|h\sqrt{g}}{C} \left[ \substack{K_L \:\:\: 0 \\ 0 \:\:\: K_T} \right] \left( \substack{ \frac{\partial \phi}{\partial L} \\ \frac{\partial \phi}{\partial T}} \right) \tag{9.9} \end{equation}\]

Where:

• \(L\) and \(T\) denote longitudinal and transverse directions

The concentration gradients in the longitudinal and transverse directions can be obtained from those in cartesian coordinates with a rotational transformation:

\[\begin{equation} \left( \substack{\frac{\partial \phi}{\partial L} \\ \frac{\partial \phi}{\partial T}} \right) = \left[ \substack{cos \theta \:\:\: -sin \theta \\ sin \theta \:\:\: cos \theta} \right] \left( \substack{\frac{\partial \phi}{\partial x} \\ \frac{\partial \phi}{\partial y}} \right) \tag{9.10} \end{equation}\]

Where

• \(\theta\) = angle of flow vector such that \(cos \theta = \frac{u}{|U|}\) and \(sin \theta = \frac{v}{|U|}\)

And similarly the dispersive and diffusive fluxes can be transformed back from flow direction coordinates to the cartesian coordinates with the reverse transform:

\[\begin{equation} \left( \substack{ \phi_{bu} \\ \phi_{bv}} \right) = \left[ \substack{cos \theta \:\:\: sin \theta \\ -sin \theta \:\:\: cos \theta} \right] \left( \substack{ \phi_{bL} \\ \phi_{bT}} \right) \tag{9.10} \end{equation}\]

Putting all of this together, and adding an isotropic diffusion term \(D_w\) (which can be used to represent wind/wave induced diffusion) yields:

\[\begin{equation} \left( \substack{\phi_{bu} \\ \phi_{by}} \right) = -hRDR^{-1} \left(\substack{\frac{\partial \phi}{\partial x} \\ \frac{\partial \phi}{\partial y}}\right) \tag{9.11} \end{equation}\]

\[\begin{equation} RDR^{-1} = \frac{|U|h\sqrt{g}}{C_z} \left[ \substack{ K_l cos^2 \theta + K_t sin^2 \theta \:\:\: (K_l - K_t) cos \theta sin \theta \\ (K_l - K_t) cos \theta sin \theta \:\:\: K_l sin^2 \theta + K_t cos^2 \theta} \right] + \left[ \substack {D_w \\ 0} \: \substack {0 \\ D_w} \right] \tag{9.12} \end{equation}\]

Note that this form is similar to that of Equations (9.5) and (9.6) but includes the off diagonal terms that arise in the matrix. Also note that the isotropic diffusion \(D_w\) in HPC is a global constant.

9.1.4 Local Constituent Transformation

In addition to pure advection and dispersion/diffusion, constituents simulated within TUFLOW AD (in both Classic and HPC) are modified by transient boundary conditions, and optional settling (for sediment constituents) and decay processes (dissolved and degradable constituents).

TUFLOW AD boundary conditions can be set to vary in time for each constituent, and can be applied to all TUFLOW boundaries that set water levels and flows (either user-specified or computed), such as HT or QT. TUFLOW AD also supports SA inflow boundaries, where flows and concentrations are used to compute mass loads that are delivered to the model domain, mixed with ambient water and then resultant concentrations computed, prior to execution of the advection routines. Consituents are not supported in rainfall boundaries, which are assumed to be pure water.

Settling of constituents to simulate removal of particulate matter from the water column has been included in the engine as a simple linear process. Once settled, constituents do not re-enter the computational domain. The local source contribution from settling is given by:

\[\begin{equation} S_\phi = -\frac{w_s}{h} \phi \: ; \: or \: S_b = -\frac{w_s}{h} b \tag{9.13} \end{equation}\]

Where:

• \(w_s\) is the settling velocity for the particles

TUFLOW AD also supports the decay of individual species (if positive decay rates are specified) and employs first order rate equations to do so. These equations draw on user defined decay rates:

\[\begin{equation} S_\phi = -\omega \phi \: ; \: or \: S_b = -\omega b \tag{9.14} \end{equation}\]

Where:

• \(\omega\) is the decay rate for the constituent

9.1.5 Groundwater

The 2023-03-AA release introduces tracking of constituents through the groundwater layers when run using the TUFLOW HPC solver, including return flow to surface water when horizontal hydraulic conductivity in included in the soils file. When this feature is used, groundwater return to the surface typically occurs at the bottom of slopes or adjacent to creeks, representative of catchment baseflows in reality. Refer to Section 9.4.6 for commands for setting initial concentrations for sub-surface layers. No dispersion is modelled in the groundwater layers. The 2d_bc “GT” type groundwater boundary does not support defined inflow values for constituents.

Advection-dispersion in groundwater is not supported in TUFLOW Classic or the TUFLOW HPC Quadtree Module.

9.1.6 Limitations

TUFLOW AD is designed to model dissolved and particulate constituent advection and dispersion in coastal waters, estuaries, rivers, and floodplains. This is achieved through solution of the 2D transport equation combined with sub-models for settling or decay.

Limitations to note include:

1. The dispersion scheme adopted by TUFLOW AD (Falconer et al., 2005) is such as to allow use of literature values for longitudinal and transverse dispersion coefficients (\(K_{L}\) and \(K_{T}\)). Users adopting literature values for these coefficients (as provided in (Falconer et al., 2005)), should however do so with extreme caution, as they are known to vary widely, and by up to several orders of magnitude. It is always preferable to use monitoring data to calibrate advection dispersion models (TUFLOW AD included) and this should be done whenever and wherever possible. If no such data is available, then literature values can be used for \(K_{L}\) and \(K_{T}\), however results need to be appropriately caveated, and TUFLOW AD predictions (as for any AD model) should be seen as qualitative or indicative at best.
   
2. Modelling predictions should also be cross checked with desktop calculations where possible. For example, this might include a hand calculation of expected salt masses in a given tidal system, with comparison made to TUFLOW AD outputs.
   
3. TUFLOW AD allows for specification and computation of large dispersion coefficients, and with the automatic substepping implementation it should generally remain stable. However, specification of large (i.e. greater than approximately 100-500) dispersion coefficients may lead to results that are not physically real or defensible. As such, (in conjunction with 2 and 3 above) results should always be sanity checked and correlated with measurements. Relying on uncalibrated model predictions is not recommended.