7  2D Domains

This chapter of the Manual discusses features specifically related to 2D model domains. 1D domain features are discussed separately in Chapter 5 and 1D/2D linking is discussed in Chapter 10.

7.1 Solvers

TUFLOW solves the depth averaged 2D Shallow Water Equations (SWE). The SWE are the equations of fluid motion used for modelling long waves such as floods, urban stormwater inundation, dam failure hydraulics, ocean tides and storm surges. They are derived using the assumption of vertically uniform horizontal velocity and negligible vertical acceleration (i.e. a hydrostatic pressure distribution). These assumptions are valid where the wave length is much greater than the depth of water. In the case of the ocean tide, the wavelength is long, and even at the deepest parts of the ocean the SWE are applicable.

For the original TUFLOW Classic solver, the mass continuity and momentum conservation equations are available in the Section 6.3 of 2018 TUFLOW Manual. For TUFLOW HPC, the 2D SWE are described by the following partial differential equations of mass continuity and momentum conservation, in the X and Y directions for an in-plane Cartesian coordinate frame of reference. The equations are:

2D Continuity:

\[ \frac{\partial h}{\partial t} + \frac{\partial(hu)}{\partial x} + \frac{\partial(hv)}{\partial y} = S \tag{7.1}\]

X Momentum:

\[ \begin{split} \frac{\partial (hu)}{\partial t} + \frac{\partial (huu)}{\partial x} + \frac{\partial (hvu)}{\partial y} - \frac{\partial \left( h \nu_t \frac{\partial u}{\partial x} \right)}{\partial x} - \frac{\partial \left( h \nu_t \frac{\partial u}{\partial y} \right)}{\partial y} + gh\frac{\partial (z + h)}{\partial x} \\ + (\frac{gn^2}{h^{4/3}} + \frac{k}{2 \mathrm{\Delta}x})h\sqrt{u^2+v^2}u = S_u - \frac{h}{\rho}\frac{\partial P_a}{\partial x} + \frac{\tau_{x,wind}}{\rho} \end{split} \tag{7.2}\]

Y Momentum:

\[ \begin{split} \frac{\partial (hv)}{\partial t} + \frac{\partial (huv)}{\partial x} + \frac{\partial (hvv)}{\partial y} - \frac{\partial \left( h \nu_t \frac{\partial v}{\partial x} \right)}{\partial x} - \frac{\partial \left( h \nu_t \frac{\partial v}{\partial y} \right)}{\partial y} + gh\frac{\partial (z + h)}{\partial y} \\ + (\frac{gn^2}{h^{4/3}}{\rho} + \frac{k}{2 \mathrm{\Delta}y})h\sqrt{u^2+v^2}v = S_v - \frac{h}{\rho}\frac{\partial P_a}{\partial y} + \frac{\tau_{y,wind}}{\rho} \end{split} \tag{7.3}\]

Where:

• \(h\) = Depth of water
• \(u\) and \(v\) = Depth averaged velocity components in X and Y directions
• \(t\) = Time
• \(x\) and \(y\) = Distance in X and Y directions
• \(z\) = Bed elevation
• \(n\) = Manning’s bed friction coefficient
• \(k\) = Energy loss coefficient, used for local energy losses
• \(\mathrm{\Delta}x\) and \(\mathrm{\Delta}y\) = 2D cell size
• \(\nu_t\) = Turbulent kinematic viscosity
• \(P_a\) = Atmospheric pressure
• \(\rho\) = Density of water
• \(S\) = Areal volume source (volume per unit time per unit area), used for rain and infiltration
• \(S_u\) and \(S_v\) = Areal momentum source terms (force per unit area per unit fluid density), used for local energy losses

The terms of the SWE can be attributed to different physical phenomena. These are:

• Propagation of the wave due to gravitational forces.
• Transport of momentum by advection.
• Horizontal diffusion of momentum or sub-grid scale turbulence (see Section 7.1.4).
• External forces such as bed friction, rotation of the earth, wind, wave radiation stresses, and barometric pressure.

For further discussion relating to the SWE, please see Section 3.3.1.

TUFLOW Classic and TUFLOW HPC use different solution schemes to solve the SWE. Both approaches are discussed in the following sections.

7.1.1 TUFLOW HPC 2D Solver

TUFLOW HPC is an explicit solver using a finite volume method. It computes the volume flow across cell boundaries. Volume cannot leave one cell without being placed in its neighbour. As a result, 2D volume is conserved and 2D mass error is 0%. The transfer of momentum across cell boundaries is computed in the same way and once external forces are considered (bed slope, bottom friction, and wet perimeters of non-uniform depth) momentum is conserved.

The explicit finite volume scheme applies the conservation of mass over the cell for calculating the rate of change of cell depth. In Figure 7.1, the cell centre (for the cell in question) is given the notation cc, while the surrounding neighbours are given the notation n1..n4. The u velocity at the left and right faces are notated u1 and u2, and the v velocities at the bottom and top faces are notated v3 and v4. The cell width and height are Δx and Δy respectively.

TUFLOW HPC uses an automatic adaptive timestep to achieve unconditional stability, as mentioned in Section 3.4.4. It solves the 2D SWE on a uniform Cartesian grid. Water depth/level is calculated at the cell centres, and velocity components at the cell mid-sides or faces in the same manner as TUFLOW Classic.

Figure 7.1: TUFLOW HPC 2D SWE Finite Volume Scheme Approach

The time rate of change for the cell averaged depth is shown in Equation 7.4.

\[ \frac{A\partial h}{\partial t} = \Phi_{1} - \Phi_{2} + \Phi_{3} - \Phi_{4} + S_{Q} \tag{7.4}\]

Where:

• \(A\) = cell area
• \(h\) = depth
• \(t\) = time
• \(\Phi_{1}\) to \(\Phi_{4}\) = the four face fluxes
• \(S_Q\) = sources (rain and infiltration)

The volume fluxes across the four cell sides and the net volume from source boundaries determine the rate of volume change and the change in depth. Source boundaries include SA, ST and RF boundaries, soil infiltration, evaporation, and any flow linkages to 1D elements via SX links. By computing the face fluxes for all model faces, and referencing these when computing the depth derivative for each cell, volume conservation is guaranteed to numerical precision. At “Head” boundaries, the above equation does not apply. Rather, the defined head (level) is directly used to calculate the water volume and associated depth in the cell.

The calculation of the cell side volume fluxes is available in either 1st or 2nd order spatially. For the 1st order scheme, this uses depth of the upstream cell (often referred to as upwinding), bounded to be greater than or equal to 0, and less than or equal to the surface elevation of the upstream cell less the bed elevation at the cell side mid-point. For the 2nd order scheme the depth at the face is computed as the average of the two cell averaged depths, however, this method in its simplest form is not total variation diminishing (TVD) and is known to be unstable. A hybrid method is implemented in which the depth at the cell face transitions from interpolated depth to the upstream depth (1st order upwinding) when the solution shows short scale reversal or upstream controlled supercritical flow.

The solution of the cell side fluxes includes the inertia and sub-grid scale turbulence (eddy or kinematic turbulent viscosity) term. Turbulence is detailed further in Section 7.1.4. The cell side fluxes may also be factored down by flow constriction factors where sub-grid-scale obstructions exist.

Due to the explicit scheme, the calculation of flux for one cell face may be performed independently of the other faces, and likewise the summation of flux for each cell volume may be performed independently of the other cell volumes. Applying the same algorithm to millions of data elements is well suited to modern multi-core CPUs,and particularly suited to GPU hardware acceleration.

The 1st order approach can experience numerical diffusion, like all 1st order schemes, and does not resolve strongly two-dimensional hydraulics (e.g. flow expansion downstream of a constriction) as well as a 2nd order solution. The 2nd order solution demonstrates no discernible numerical diffusion, and resolves complex 2D hydraulics, including hydraulic jumps as demonstrated using the UK EA Benchmark Flume Test 6A Collecutt & Syme (2017). When running HPC solver the 2nd order spatial scheme is the default and recommended approach.

For further details on the scheme, refer to Collecutt & Syme (2017). Note that at the time this paper was written, the scheme utilised cell centre definitions for velocity, which was prone to a zero-energy ‘checker-board’ mode in the solution. Subsequent to the paper, cell mid-side points were adopted for the definition of u and v velocities which has eliminated the checker-board mode with only very minor changes to the results.

7.1.2 TUFLOW Classic 2D Solver

The scheme is unlike the TUFLOW HPC solution in that it solves the same shallow water equations implicitly using matrices, allowing much larger timesteps (hence why it is important to monitor mass error in implicit schemes such as TUFLOW Classic to check that the solution is converging).

An Alternating Direction Semi-Implicit (ADI) finite difference method is used for TUFLOW Classic’s computational procedure. It was originally based on the work of Stelling (1984). The method involves two stages per timestep, each having two steps, giving four steps overall. Two of the steps sweeps through the model domain solving tri-diagonal matrices, hence the “semi-implicit”. Due to the implicit solution, TUFLOW Classic is not suited to being parallelised for multi-core CPUs or for GPU acceleration.

Step 1 solves the momentum equation in the Y-direction for the Y-velocities. The equation is solved using a predictor/corrector method, which involves two sweeps. For the first sweep, the calculation proceeds column by column in the Y-direction. If the signs of all velocities in the X-direction are the same the second sweep is not necessary, otherwise the calculation is repeated sweeping in the opposite direction.

The second step of Stage 1 solves for the water levels and X-direction velocities by solving the equations of mass continuity and of momentum in the X-direction. Tri-diagonal equations are generated and solved across the 2D domain by substituting the momentum equation into the mass equation and eliminating the X-velocity. The water levels are calculated and back substituted into the momentum equation to calculate the X-velocities. This process is repeated for a recommended two iterations. Testing on a number of models showed there to be little benefit in using more than two iterations unless there are rapid changes in the hydraulic conditions per timestep as may occur with modelling inundation from a dam break.

Stage 2 proceeds in a similar manner to Stage 1 with the first step using the X-direction momentum equation and the second step using the mass equation and the Y-direction momentum equation.

The solution as formulated by Stelling has been enhanced and improved to provide much more robust wetting and drying of elements, upstream controlled flow regimes (e.g. supercritical flow and upstream controlled weir flow), modifications to cells to model structure obverts (e.g. bridge decks) and additional energy losses due to fine-scale features such as bridge piers.

TUFLOW Classic solves the 2D SWE on the same uniform Cartesian grid as used by TUFLOW HPC. Water depth/level is calculated at the cell centres, and velocity components at the cell mid-sides or faces.

7.1.3 Cell Schematisation

TUFLOW 2D models discretise the real world as a grid of connected square cells. A fixed grid 2D domain is defined by a bounding rectangle in the same manner a computer screen or digital photo is made up of a grid of pixels.

The physical properties (i.e. ground elevations, surface roughness, etc) of a 2D cell are defined as a minimum at the cell’s centre, mid-sides and corners as described in Section 7.2.1. High-resolution topographic detail for the cell’s storage and conveyance across the cell sides (faces) can be incorporated by using the Sub-Grid Sampling (SGS) feature available in TUFLOW HPC, as documented in Section 7.3.3. SGS has substantial accuracy benefits. It enables the use of much larger cell sizes without loss of hydraulic conveyance accuracy. It also results in less output sensitivity to the grid orientation.

7.1.3.1 Computational Points

To fully understand how TUFLOW functions, understanding of the computational role of the 2D cell and its four sides (also referred to as faces) is important. In the following sections, these model components are referred to as the ZC, ZU, ZV and ZH points, as shown in Figure 7.2. The description of these points are:

• ZC – the cell centre with ZC being the elevation used computationally within the cell.
• ZU – middle of the cell side of the Y-axis face with ZU being the elevation used computationally along the cell side.
• ZV – middle of the cell side of the X-axis face with ZV being the elevation used computationally along the cell side.
• ZH – corner of the cell. The ZH point plays no role in the computational hydraulics but is used for some output formats that require output at the cell corners.

Figure 7.2: Location of Zpts and Computation Points

7.1.3.1.1 ZC Point

The cell centre or ZC point:

• Defines the location where water levels are computed based on mass balance equation. Simply put, the net volume of water entering (or leaving) the cell across the four cell sides (or faces) must equal the change in volume of the cell over a timestep.
• The volume of a cell can either be simply represented as a square, flat bucket and calculated as the cell centre depth of water multiplied by the cell area (see Section 7.2.4.1) or, if Sub-Grid Sampling (SGS) (Section 7.3.3) is enabled, a curve of volume versus depth.
• The ZC value is typically the elevation at the cell centre or, if using SGS, elevations are sampled across the cell, as outlined in Section 7.3.3. The ZC value plus the Cell Wet/Dry Depth is the water surface elevation that controls when a cell becomes wet or dry (note that cell sides can also wet and dry).
• The ZC value also determines the bed slope when testing for the upstream controlled flow regime in TUFLOW Classic (see Section 7.4.2) and TUFLOW HPC (Section 7.3.4).

7.1.3.1.2 ZU and ZV Points

The cell sides or faces:

• Control how water is conveyed from one cell to another using the momentum equation, or when upstream controlled flow occurs, the relevant flow equation (eg. weir equation).
• The flow area across a cell face can either be simply represented as a rectangular section (i.e. cell side centre depth times the cell width) (see Section 7.2.4.1) or, if Sub-Grid Sampling (SGS) (Section 7.3.3) is enabled, a curve of flow area versus depth.
• Are deactivated if the whole of the cell has dried based on the ZC elevation as described above.
• Can wet and dry when the whole cell is wet (see Cell Wet/Dry Depth). This allows for the modelling of “thin” obstructions such as fences and thin embankments relative to the cell size (e.g. a concrete levee).
• Surface roughness information (e.g. Manning’s n) is inspected at the ZU and ZV points. See Bed Resistance Cell Sides and Section 7.2.4.1.

7.1.3.1.3 ZH Point

The cell corners:

• Play no role in the computational hydraulic calculations.
• Are used by the more common map output formats to map results due to mesh based outputs requiring elevations and results at the element (cell) vertices.

7.1.4 Turbulence

Turbulence within rivers plays a significant role in determining the mean flow velocity field and is integral to the overall energy loss mechanism. During a significant flood event, the majority of the flow is carried within, or near to, the river or open channel system (the primary flow path), and the water levels on the surrounding floodplains are principally determined by the energy losses along this flow path. Manning’s equation is accurate for straight or slowly varying channels but does not calculate energy losses due to sudden changes in flow direction or velocity, and it does not capture super-elevation at bends. To automatically and accurately capture these additional energy losses, a scheme must be physics based with a turbulence closure model. The turbulence closure model must be thoroughly bench-marked against appropriate test cases at a range of resolution scales. It is important that the turbulence closure model is cell size and timestep independent (i.e. the turbulence closure model parameters do not need to be adjusted if cell size or timestep changes). Collecutt et al. (2020) provides benchmarking of the TUFLOW turbulence schemes.

There are three approaches available to model losses associated with turbulence when using TUFLOW HPC, these are discussed in Section 7.3.2. The default scheme for TUFLOW HPC is the Wu approach (Section 7.3.2.4). There are two approaches available to model turbulence when using TUFLOW Classic, these are discussed in Section 7.4.1. The default scheme for TUFLOW Classic is a hybrid approach that includes elements of both the Constant and Smagorinsky approaches (Section 7.4.1.2).

7.1.4.1 Dry Wall Treatment

All of TUFLOW’s viscosity approaches feature an enhanced treatment of the viscosity term at dry boundaries. The enhanced boundary treatment corrects for unrealistic flow separation that would otherwise occur at the wet/dry interface. Figure 7.3 presents an example of a benchmark model test with and without the enhanced boundary treatment. The top image is the result without the enhanced boundary treatment. Non-physical flow separation is observed along the oblique wet/dry boundary. The enhanced approach is shown in the bottom image. It does not suffer from the flow separation, and the velocity increases gradually from left to right as the water depth gradually shallows (the model has a horizontal bed). The correct water surface slope is produced when compared with theory.

All of the viscosity approaches are cell-centred. This guarantees symmetry is achieved in hydraulic results when using a symmetrical model.

The top image in Figure 7.3 shows flow separation along a dry oblique boundary without enhanced treatment of viscosity term. The bottom image shows a correct velocity distribution using enhanced treatment of viscosity term at dry boundaries.

Figure 7.3: Effect of Enhanced Dry Boundary Viscosity Term Treatment

7.2 Common Functionality

The following sections are the same (or mostly so) between TUFLOW HPC and TUFLOW Classic.

7.2.1 Defining the Domain

Each fixed grid 2D domain is treated as a rectangle at any orientation. The orientation and dimensions of each 2D domain are defined using .tgc file commands. For the orientation it is recommended that the X-axis falls between 90° and –90° of East due to some post-processing software only operating within this range.

Several options are available for setting the 2D domains grid location and orientation. The options are:

• Using a four-sided polygon in a GIS layer to define the 2D grid orientation and dimensions (see Read GIS Location).
• Using a line (two vertices only) in a 2d_loc GIS layer (attributes shown in Table 7.1) to define the origin (first point of line) and orientation (based on the second point of the line) of the X-axis (see Read GIS Location), and Grid Size (N,M) or Grid Size (X,Y) to set the 2D grid X and Y dimensions.
• Using Origin, Orientation or Orientation Angle, and Grid Size (N,M) or Grid Size (X,Y). No GIS layers are required for this option. This option is used in TUFLOW Tutorial Module 1.
• Using a DEM to set the size and location of a 2D domain (see Read Grid Location). This option is useful where the model extent is the same as the DEM.

It is recommended during the initial model build process, to view the _dom_check file to ensure that the model domain is set up as intended.

Note, if using the TUFLOW HPC Quadtree functionality, the 2D domain is set up in the Quadtree Control File (.qcf), see Section 7.3.1.2.

Table 7.1: Location (2d_loc) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Comment	Optional field for entering comments. Not used.	Char(250)

After establishing the domain’s origin, orientation and extent, the .tgc command Cell Size is used to define the domain’s fixed grid resolution. For example, the below command in the .tgc sets the cell size to 10m (or 10ft is using Units == US Customary):

Cell Size == 10

Note

It is not necessary to specify domain dimensions that are an exact multiple of the domain’s Cell Size (e.g. when using the Grid Size (X,Y) command).

Within this computational domain, 2D cells can be set to be active or inactive, as described in Section 7.2.2. The redundant (inactive) areas around the edges of the domain’s bounding rectangle are automatically removed from the computation. However, when allocating memory for processing of the arrays this is still accounted for. To determine the redundant area, search the TUFLOW log file (.tlf) for “redundant”, as shown in Figure 7.4. If the redundant value is large, revising the model domain dimensions will reduce the memory usage during startup.

Figure 7.4: Redundant Perimeter Sections Reported in the TLF

7.2.2 Active / Inactive Cells

Each 2D cell is assigned a code to indicate its role. The available code types are listed in Table 7.2. The default value is one (1) for active (i.e. the cell can be wet or dry during a simulation).

Commands used to modify the cell codes are:

• Set Code in the .tgc file. This option sets the cell codes for the entire domain.
• Read GIS Code in the .tgc file. This option sets the cell codes based on the value in 2d_code input polygons, see Table 7.3.
• Read Grid Code in the .tgc file. This option sets the cell codes based on the values of an input raster grid.
• Read GIS Code BC in the .tgc file. This option sets the cell codes from the 2d_bc input polygon. The Type attribute must be set to “CD” and the code value is taken from the 2d_bc f attribute, see Table 8.6.

Note

When using the Read GIS Code command, code values are extracted from objects in a 2d_code layer (see Table 7.3). When using the Read GIS Code BC, code values are extracted from objects in a 2d_bc layer (see Table 8.6). Confusing the two GIS layer types and commands will result in WARNING 2320 message being issued.

A typical approach is to set all the cells to be inactive using Set Code == 0 and then read a GIS layer defining the active cells with Read GIS Code == <layer>. The GIS layer regions would have the attribute set to 1 for active. For example:

Set Code == 0                                     ! Set all cells to inactive (i.e. Code 0)
Read GIS Code == ..\model\gis\2d_code_M01_002.shp ! Reads in a polygon that sets cells to active (i.e. Code 1)

Boundary Cells are automatically set to 2 along external boundaries and links (e.g. HX, HS, HT, QT boundaries). This code value of 2 (see Table 7.2), is output in the TUFLOW grid check file and a value of 2 should not be used in input commands. The _grd_check file can be used to view the active cells (and their code values e.g. 1, 2 or -1) within the domain.

For TUFLOW Classic models containing multiple 2D domains (see Section 10.7.2), a useful option for setting the cell codes is the INVERT flag (see Read GIS Code Invert or Read GIS Code BC Invert), which allows the same code layer/polygons to be used for activate and deactivate regions in multiple 2D domains using a single GIS input.

TUFLOW automatically strips any redundant rows and columns around the active area of the model to reduce simulation times. This is described in Section 7.2.1.

Table 7.2: Cell Codes

Type	Code	Description
Inactive Cells	0	Inactive cells are cells that are totally removed from the computation. Maximising the area of inactive cells reduces computation time and output file sizes. For a fixed grid 2D domain, inactive cells still consume memory during a simulation, as these domains are stored on a row column rectangular grid. For a quadtree grid, inactive cells will consume temporary memory (CPU RAM) during model startup, but do not consume memory (CPU or GPU) during the simulation.
Active Cells	1	Active cells are cells that can wet and dry during a simulation.
Boundary Cells	2	Boundary cells indicate cells that are an external boundary (including some types of 1D/2D dynamic links). There should be an active cell on one side and an inactive or null cell on the other at an external boundary. If an external boundary is digitised inside an active area (i.e. not along the active area boundary), water can flow in both directions either side of the boundary line (only in rare situations would the modeller require this). Note: It is not necessary to manually specify boundary cells. Boundary lines are digitised in GIS layer(s) and TUFLOW automatically assigns the boundary code to the cells (see Section 8.4.1 and Read GIS BC).
Null Cells	-1	Inactive cells used to deactivate cells within the active domain. Null cells may be preferred to inactive cells as they are not excluded from the output mesh structure (eg. not excluded from the XMDF mesh). For two simulations to be compared in some post-processers (eg. SMS), they must have exactly the same mesh. For example, if an area in a model is removed (e.g. filling part of a floodplain), use null cells or raise the ground elevations in preference to using inactive cells so that the two simulations can be compared. Setting null cells is not supported when using the TUFLOW HPC Quadtree (see Section 7.3.1) functionality.

Table 7.3: 2D Code (2d_code) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Code	The code value (see Table 7.2) to be assigned to cells falling on or within the object.	Integer

7.2.3 Data Layering

TUFLOW models are set up using data layering. This is an important, and very useful concept.

Commands are applied in sequential order; and it is possible to override previous information with new data to modify the model in selected areas. For topography modifications, this is extremely useful where a base dataset exists, over which areas need to be modified to represent other scenarios such as a proposed development. This eliminates or minimises data duplication. This concept of layering datasets may also be applied to other GIS , including (but not limited to) 2d_mat, 2d_code and 2d_soil layers. For example, setting all cell codes inactive, then reading in a GIS polygon covering the active cells (as discussed in Section 7.2.2).

Figure 7.5 shows a visual representation of how TUFLOW interprets the following set of .tgc commands:

Read Grid Zpts == grid\DEM_5m.tif
Read Grid Zpts == grid\Bathy_1m.tif
Read GIS Z Shape == gis\2d_zsh_building.shp

The command lower in the TGC will take precedence over a command above where data overlaps and the command has the same functional purpose (eg. topography definition). The DEM_5m layer is read in first. The elevations are then replaced by the overlapping Bathy_1m layer. A 2D Z Shape (2d_zsh_buildings) is then read in which raises all elevations inside the polygon by 5m. The result is shown in the DEM_Z. The DEM_Z is a check file representing the elevation data that TUFLOW has used for the hydraulic calculations. For more information on check files see Section 14.6.

Figure 7.5: Visual Representation of Data Layering in TUFLOW

Data layering is discussed in further detail for elevation updates in Section 7.2.6.

7.2.4 Sampling of Data Sets

To setup the boundary of the 2D grid and assign data values to the cells from DEMs, TINs, land-use polygons, soil polygons and so forth, these data layers need to sampled, interrogated or interpolated. There are two approaches to sampling datasets in TUFLOW:

• Traditional approach (Section 7.2.4.1)
• Sub-Grid Sampling (SGS) approach - recommended (Section 7.3.3)

The traditional approach by 2D solvers is to sample a single value at the 2D cell centre (ZC). See Section 7.1.3 for a description of 2D cell control points. TUFLOW also supports sampling values, for some data types, at the centre of each cell side (ZV and ZU) and the corner of the cell (ZH). The cell side sampling for example, enables a cell side to be set to the height of a levee that is narrower in width than the 2D cell. Subsequently controlling when water can move from adjacent cells depending on the water level relative to the levee crest elevation. The traditional sampling approach is explained in Section 7.2.4.1.

As of the 2020-01 TUFLOW release, TUFLOW HPC supports the sampling of elevations at sub-grid locations, referred to as Sub-Grid Sampling (SGS). The TUFLOW Classic solver does not support SGS due to the nature of the implicit solution and preferred use of a fixed timestep. The many benefits of using SGS are detailed in Section 3.2.3, implementing the SGS approach is described in Section 7.3.3.

7.2.4.1 Traditional Sampling Approach

Topographically, the cell is treated as having a flat horizontal bed at a height set to the cell centre (ZC) elevation, and the cell faces as having flat, rectangular shaped sides at a height set to the cell mid-side (ZU and ZV) sampled elevations as illustrated in Figure 7.6. Using the traditional sampling approach, the topography of a cell is treated as follows:

• The cell volume is represented as a square bucket with a flat horizontal bed, and simply calculated as the cell centre depth times the cell area.
• The flow area across a cell face used for the momentum and mass balance equations is simply represented as a rectangular section (i.e. cell side centre depth times the cell width).
• The bed resistance term in the momentum equation (e.g. Manning’s equation) uses the Resistance Radius approach (i.e. the cell face radius value as used in Manning’s equation is set to the depth and the cell width or cell size is used for the wetted perimeter).

Figure 7.6: 2D Cell Topography - Traditional Approach

7.2.5 Cell Selection Approach

Input GIS data with line geometries, such as 2d_bc boundaries (Section 8.4.1) and 2d_zsh topographic modifiers (Section 7.2.6.2), must select or intersect with 2D cells in order for those cells to be modified. This selection is based on the cross-hair approach, or a modified spider web approach (for Quadtree - see Section 7.3.1), as described below.

7.2.5.1 Cross-Hair Approach

Cross-hairs are imaginary lines that extend from cell mid-side to mid-side. Cells are selected only if an input GIS line intersects these imaginary cross-hairs, as demonstrated in Figure 7.7. This means that if a line starts or crosses inside a cell, and does not intersect with the cross-hairs, that cell will not be selected and modified, or utilised by the input.

Figure 7.7: Cell Cross-Hair Selection Approach

7.2.5.2 Spider Web Approach

Since the 2020-10-AA TUFLOW release, a modified spider web cell selection approach has been adopted for HPC solver Quadtree models (see Section 7.3.1). Instead of using cross-hairs, the spider web approach draws an imaginary line from the cell centre to each connected cell face, as shown in Figure 7.8. In areas of consistent cell size, this will look identical to the cross-hair approach, however at transition zones, there will be a spider web line connecting to each nested grid cell, creating a continuous intersection line and removing the potential for gaps in the cell selection.

At external boundary lines where the spider web test may fail (since there are no neighbouring cells on one or more sides), TUFLOW will automatically buffer refined areas to fill gaps.

Figure 7.8: Cell Spider Web Selection Approach

7.2.6 Elevations

2D domain elevations are defined in the .tgc file. As mentioned in Section 7.2.3, a powerful feature of TUFLOW is its capacity to build the 2D elevations from any number of GIS layers and/or TINs. The typical approach adopted is as follows:

1. All elevations in the model start with an un-initialised value of 99999. TUFLOW will output an error if any elevations of 99999 (or higher) occur after the processing of the elevations. This is to indicate that elevations in active cells have not been initialised or set.
2. A default elevation is specified first in the .tgc file using Set Zpt. An inundation free elevation is usually specified. If your elevation data sets are intended to cover your entire active area this step can be omitted so as to trigger the error described above in case elevation data are unintentionally missing (e.g. unintended null areas in a DEM or accidental omission of a layer).
3. The elevations read directly from a DEM using Read Grid Zpts, or read from TIN in the SMS, 12D or LandXML TIN formats (see Read TIN Zpts).
4. If there are areas where data is either missing or erroneous, these can easily be corrected via interpolation using Read GIS Z Shape regions or polygons. This is achieved by digitising a polygon around the missing or erroneous data. TUFLOW interpolates elevations across the region based on the existing Zpts around the perimeter of the polygon. An example is shown in Section 7.2.6.2. Missing or erroneous data often occurs with aerial surveys where the sampling is in areas of thick vegetation, water or where post processing has poorly filtered the removal of elevated objects (e.g. buildings).
5. Sometimes there is a need to remove cut and fill works from the topography. For example, model calibration often requires the removal of infrastructure, such as levees, road embankments and developments that may not have existed at the time of the historic event to back date the topography within the model. This is easily done using Read GIS Z Shape by simply digitising polygons around the various features.
6. The base elevations set up in the previous step(s) can be modified to represent hydraulic controls, proposed works, failure of a flood defence wall, etc. Some examples are:
   • The crests of road/rail embankments, levees, fences and other solid obstructions are easily inserted using Read GIS Z Line, Read GIS Z Shape or Read GIS Z HX Line. Read GIS Z Shape is particularly powerful as 3D lines can be given a thickness making it very easy to quickly raise, or lower, elevations along a road alignment or a diversion channel where the width of the embankment/channel is wider than the 2D cell size.
   • The proposed cut and fill for a development or other works can be incorporated using Read GIS Zpts, Read GIS Z Shape, Read TIN Zpts or Create TIN Zpts. These powerful commands can set elevations based on regions. Within regions TINs can be generated from points and lines, and the perimeter of the TIN can be automatically merged with the existing Zpts. A TIN of the cut and fill produced by SMS or 12D can be read directly into TUFLOW using Read GIS Zpts.
7. If there is a need to simulate the failure of a flood defence wall or road/rail embankment, or the collapse of fences, Read GIS Variable Z Shape can be used to control the collapse of the embankment or fences over time. The collapse can be triggered to occur at a specified time, when a water level reached somewhere within the model, or based on the water level difference between two locations.

A 2D domain’s Zpts are built up using one or more of the commands shown in Table 7.4 and discussed in the following sections.

Table 7.4: 2D Zpt Commands

Command	Description
Set Zpt	Sets all Zpts over the whole 2D domain to the same value. Useful for providing an initial elevation prior to other commands as some Zpts in inactive (land) parts of the model may not receive a value. The default value for all Zpts is 99999. Every Zpt must be assigned a value, essentially making this command mandatory.
Read Grid Zpts	Directly interrogates a raster grid (e.g. tif, .flt or .asc) to define Zpt elevations. This command is similar to Read TIN Zpts but works on a grid rather than a TIN.
Read TIN Zpts	Reading of a TIN to set the Zpt values within the TIN.
Read GIS Z Shape	Powerful command to modify Zpt values using points, lines and polygons. Lines can vary in width from just the cell sides (thin), whole cells (thick) or be assigned a width (thickness). TINs are created within the polygons and incorporate elevations from points and lines that fall within the polygon. The perimeter of the polygon can be merged with the current Zpt values in part, or in its entirety. Read GIS Z Shape and Create TIN Zpts are excellent for removing bad data areas and for filling in null areas where the aerial survey has provided poor or no data. Another example is to remove buildings from a DEM.
Read GIS Variable Z Shape	Allows the user to define the eroded 3D shape of a section of the 2D domain, specify the period for collapse, and how the collapse is triggered (i.e. at a specified time or when a water level is reached or when a water level difference is exceeded). Raising of the Zpts over time is also permitted (e.g. to model the influence of flood defences during an event or a landslide filling a river).
Create TIN Zpts	Useful for having TUFLOW create a TIN using polygons for TIN boundaries and points and lines within the polygons to create the TIN. If no points are snapped to the perimeter vertices of a polygon, the elevations around the polygon’s perimeter are merged with the current Zpt values. The resulting TIN can optionally be exported to SMS and 12D TIN formats.
Read GIS Z Line	This is a legacy feature and it is recommended to use the Z Shape functionality (Section 7.2.6.2) instead.
Read GIS Z HX Line	Similar to Read GIS Z Line, but uses HX lines and ZP points in a 2d_bc layer (see Table 8.6) to adjust the 2D cell elevations along HX lines.
Read GIS Zpts	Typically used for simple modifications of sections of the topography. Examples are filling an area (defined by a region or polygon object) to the same elevation or dredging (lowering) a section of river using Read GIS Zpts ADD. Read GIS Z Shape offers greater functionality and may be preferable to using Read GIS Zpts to modify Zpts.
Maximum Points Maximum Vertices	Use these commands to increase the maximum number of elevation points or maximum number of vertices in a single line or polygon.
Default Land Z	Now rarely used in lieu of Set Zpt.
Interpolate ZC Interpolate ZHC Interpolate ZUV Interpolate ZUVC Interpolate ZUVH	Allows the interpolation of Zpts from other types of Zpts. Now rarely used as nearly all models assign values directly to all the Zpts. The original TUFLOW code only required input of ZH points, and Interpolate ZUVC provided a tool for interpolating the other Zpts. Models with “bumpy” terrain, such as that from airborne laser surveys, might benefit from using Interpolate ZHC or Interpolate ZUV. Models through urban areas where the DTM includes the buildings may benefit from using Interpolate ZC ALL LOWER, which reduces the amount of cells that become blocked out due to high ZC elevations from buildings.
ZC == MIN(ZU,ZV)	Rarely used. Sets the ZC (cell centre elevation) to the lower of the ZU and ZV (cell sides).

7.2.6.1 Direct Reading of DEM Grids

The use of the .tgc command Read Grid Zpts allows TUFLOW to directly interrogate (point inspect) a DEM to set the Zpt elevations. This command is similar to Read TIN Zpts but works on a grid rather than a TIN.

Grid formats currently supported are:

• GeoTIFF (.tif)
• Geopackage raster (.gpkg)
• Binary grid format (.flt/.hdr)
• ESRI ASCII Grid format (.asc)

Nearly all GIS software will support one of the above grid formats.

The GeoTIFF format is the default output format and preferred input format. It is faster to read and write compared to the .flt and .asc formats, and also supports compression. The format of grids can be converted (e.g. from .flt to .tif), the TUFLOW Wiki Raster Format Conversion details how to do this.

Like other .tgc commands, the command Read Grid Zpts can be specified multiple times. An option to specify ADD, MIN or MAX in the same way as for other similar commands is also available.

Clip regions can be specified as a second argument in the command Read Grid Zpts (and also Read TIN Zpts) by reading in a GIS layer containing one or more polygons to clip the area of Zpts to be inspected. For example, the command below will only assign elevations to Zpts that lie inside polygons within the 2d_clip_DEM layer:

Read Grid Zpts == grid\DEM_M01.tif | DEM\2d_clip_DEM.shp

The attributes of the clip layer are not used, and only polygons are processed. As such, any of the TUFLOW empty files can be used as the template for this layer. Polygons can have holes in them if required.

This is particularly useful for clipping out a TIN or DEM due to unwanted or irregular triangulation around the periphery, especially for secondary TINs/DEMs of proposed developments lying within the primary TIN/DEM.

Note

For your base Zpts from the primary DEM or TIN, do not clip this with your active 2d_code layer as this will cause problems with Zpts along any external 2d_bc boundaries. If no clip layer is specified, the Read Grid Zpts or Read TIN Zpts commands assign all Zpts falling within the Grid / TIN an elevation irrespective of whether a cell is active or inactive.

When reading grids into TUFLOW, the same interpolation approach occurs for all formats. The procedure TUFLOW follows is:

• The DEM elevation is assumed to be in the centre of the DEM grid cell (pixel), as shown by the four red circles in the picture below.
• A midpoint vertex (green circle) is defined in the middle of the four DEM elevations.
• The elevation of the midpoint vertex (green circle) is equal to the average of the four DEM elevations (red circles).
• The DEM grid cell centre points (red circles) and midpoint vertex (green circle) are used to create a TIN.
• Elevations within the TIN are interpolated from the associated TIN points. For example, a Zpt is included in the image below as a yellow point. The elevation assigned to that Zpt is the planar (linear) interpolation of the surrounding three elevation points associated with the top TIN triangle in the image.

Figure 7.9

7.2.6.2 Z Shape Layers (2d_zsh)

Read GIS Z Shape provides flexibility for manipulating Zpt values by combining various 2d_zsh GIS layer geometries and attribute options. This includes the ability to represent 3D breaklines and TINs.

For example, Z shape layers allow breakline features such as roads, railways, levees, creeks and drains to be applied quickly and easily. This is particularly useful where TUFLOW’s fixed grid discretisation may not capture crest or invert elevations from the DEM. Footprints of historical or proposed development (e.g. buildings or earthworks) can also be represented with ease. In many cases, Z shape layers remove the need to modify the DEM raster directly.

Table 7.5, Table 7.6 and Table 7.7 describe the attributes and options available for each geometry (point, line and polygon). Example use cases are provided in Section 7.2.6.2.2.

Zpt values are primarily modified using 2d_zsh polygon and line geometries. Additional flexibility is available when they are grouped with points, lines and other attribute options. 2d_zsh points do not directly modify Zpt values. Instead, they assign elevations to lines and polygons or help define TINs when used with a polygon. Line behaviour depends on classification. Lines are considered Thin, Thick or Wide based on the Shape_Width_or_dMax attribute value (see Table 7.6).

The following describes how each geometry modifies Zpt values:

1. Polygons:
   • All polygons are triangulated using the process described for Create TIN Zpts in Section 7.2.6.4. However, there are two key differences when applying this process to Z shape (2d_zsh) layers:
     • If a line is to be used for TIN generation, “TIN” must be specified for the Shape_Options attribute (see Table 7.6).
     • Merge polygon handling differs slightly. For Create TIN Zpts, polygon vertices without defined elevations from snapped points are automatically merged. For Read GIS Z Shape polygon layers, vertices can be assigned a special value of -99999 via a snapped point, allowing elevations to be interpolated from neighbouring defined vertices (see Table 7.5).
   • As shown in Figure 7.10:
     • The processed TIN triangulation polygons can be viewed in the _sh_obj_check layer.
     • Elevations of all control Zpts (ZC, ZH, ZU and ZV) within the polygon are modified based on the TIN triangulation.
   • When Sub-Grid Sampling (SGS) is enabled (with TUFLOW HPC) – see Section 7.3.3.1.2:
     • The 2d_zsh polygon modifies all SGS sample points within the polygon extent, based on the TIN triangulation (see Section 7.3.3.1.2).
     • If a crest is included in the modified topography and is not aligned exactly with a cell face, lower elevation values may be sampled in the SGS cell face width vs depth curve, potentially allowing water to leak at a lower elevation (see Section 7.3.3.1.3 – SGS Breakline Detection Delta). To avoid this, 2d_zsh breaklines with RIDGE or MAX specified for the Shape_Options attribute (see Table 7.6) should be added along the polygon or TIN crest to enforce the elevations.
   • When Quadtree is enabled (with TUFLOW HPC) – see Section 7.3.1.3:
     • For 2d_zsh polygons with merged elevations at undefined vertices extending across multiple quadtree nested layers, it is recommended to review the _zsh_zpt check file to ensure appropriate representation. CHECK/WARNING 2934 will be written in this instance.

Figure 7.10: Example of 2d_zsh polygon with Zpt modification

2. Thin and Thick Lines:
   • Line type is determined by the Shape_Width_or_dMax attribute value (see Table 7.6):
     • Thin lines have a value of zero (0).
     • Thick lines have a value greater than zero (0) and less than or equal to 1.5 times the 2D Cell Size.
   • As shown in Figure 7.11:
     • Thin lines modify cell faces (ZH, ZU and ZV Zpts) where the line intersects the ends of the model cell’s imaginary cross-hairs that extend from cell mid-sides to mid-side (see Section 7.2.5).
     • Thick lines modify entire cells (ZC, ZH, ZU and ZV control Zpts) along the line. The cells are selected if the thick line intersects any part of a cell’s imaginary cross-hairs that extend from cell mid-sides to mid-sides (see Section 7.2.5).
     • The MIN or GULLY Shape_Options (see Table 7.6) override the line type and define a path through the cell face and cell centre (ZU, ZV and ZC control Zpts only), lowering Zpt values only if the new value is lower than the current value to ensure a connected flow path. Note: The GULLY option is not intended to accurately represent the shape of a waterway and is not recommended for use with SGS, as it can overestimate conveyance. See Section 7.3.3 and the TUFLOW Wiki Quadtree and SGS FAQ page for more details.
     • Elevations assigned to control Zpts selected by either a thin or thick line are based on the perpendicular intersection of that Zpt with the line. If there is no perpendicular intersection, the nearest vertex on the line is used.
   • From the 2023-03 release onwards, TUFLOW HPC applies the full weir equation in the 2D domain along thin and thick breaklines using the advanced weir formulation (see Section 5.7.4.3). Default weir flow parameters can be adjusted separately for thin and thick breaklines – see Section 7.3.4.
   • When Sub-Grid Sampling (SGS) is enabled (with TUFLOW HPC) – see Section 7.3.3.1.2:
     • The 2d_zsh line modifies elevations based on the cell’s control Zpts (unlike polygons, which directly modify all SGS sampled elevations), as described in Section 7.3.3.1.3 – SGS Z Shape Line Approach.
     • By default, thin lines apply a gradient along the cell face using side and corner Zpts (ZH, ZU and ZV), while thick lines apply a sloping cell area for the cell volumes using ZC, ZU, ZV and ZH Zpts.
     • Use of the MIN or GULLY Shape_Options (see Table 7.6) is not recommended, as it can overestimate conveyance – see Section 7.3.4.
     • It is strongly recommended to add MAX or RIDGE lines (see Table 7.6) for every hydraulic control (e.g. levees (artificial or natural), road or rail embankments) to reduce “leakiness” caused by SGS sampling. See Section 7.3.3.1.3 – SGS Breakline Detection Delta.
   • When Quadtree is enabled (with TUFLOW HPC):
     • A modified cell cross-hair selection approach, referred to as the “spider web” approach, is applied. Instead of using cross-hairs, this method draws imaginary lines from the cell centre to each connected cell face, ensuring continuous cell selection. See Section 7.2.5.2.

Figure 7.11: Example of 2d_zsh GULLY, thin and thick line Zpt modification

3. Wide Lines:
   • Wide lines have a Shape_Width_or_dMax (see Table 7.6) value greater than 1.5 times the 2D Cell Size.
   • As shown in Figure 7.12, a buffer polygon is created around the line, extending to a distance of half the Shape_Width_or_dMax value (on either side of the line). All control Zpts (ZC, ZH, ZU and ZV) within the buffer polygon are assigned elevations based on a perpendicular intersection with the line and the interpolated line elevation at that intersection. This buffer polygon is written to the _sh_obj_check layer.
   • When Sub-Grid Sampling (SGS) is enabled (with TUFLOW HPC) – see Section 7.3.3.1.2:
     • The 2d_zsh line will modify elevations based on the cell’s control Zpts (rather than directly modifying all sampled SGS elevations as for polygons), as described in Section 7.3.3.1.2 and Section 7.3.3.1.3 – SGS Z Shape Line Approach.
     • By default, wide lines apply a gradient along the cell face using the side and corner Zpts (ZH, ZU and ZV), and a sloping cell area for the cell volumes (using ZC, ZU, ZV and ZH Zpts). SGS points are sampled from these.
     • For partially selected cells at the buffer’s boundary, a sloping cell area is determined based on the cell’s modified Zpts within the buffer and the unmodified Zpts outside the buffer, and all SGS sample points will be updated within the partially selected cells.
     • It is recommended to consider alternative TIN based modification methods via Read GIS Z Shape 2d_zsh polygons or Create TIN Zpts to improve representation.

Figure 7.12: Example of 2d_zsh wide line Zpt modification

Table 7.5: 2D Z Shape (2d_zsh) Point Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Z	Elevation of the point. Points are only used to assign elevations along lines and polygon perimeters, or for the creation of TINs within polygons. They do not modify Zpts directly. An elevation of -99999 has a special meaning when the point is snapped to a vertex of a polygon. The -99999 indicates to ignore the elevation at that vertex and of any automatically inserted vertices between that vertex and the two neighbouring vertices. Instead the elevations are based on the elevations of the neighbouring vertices. If a neighbouring vertex also has a -99999 point snapped to it, the next vertex is used, and so on. This feature is very useful, as illustrated in the example above. Note: This special meaning does not apply to points snapped to lines, only to polygons. For lines that have at least one point snapped to a vertex, any remaining vertices where no points are snapped will be interpolated by default.	Float
2	dZ	Change in elevation of the Z attribute of the 2d_zsh object at the point (does not change the Zpt elevations directly). To leave Z unchanged, set to zero. Useful for easily adjusting the height of the shape object without having to change the original Z value. For example, if the Z attribute has been set to 5mAHD, a dZ value of 0.1m will raise the point’s Z elevation to 5.1mAHD.	Float
3	Shape_Width_or_dMax	Not used.	Float
4	Shape_Options	Not used.	Char(20)

Table 7.6: 2D Z Shape (2d_zsh) Line Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Z	Elevation of the line if no points are snapped to the line. If the ADD option is specified, the value entered is used to increase (positive Z values) or decrease (negative Z values) the elevation of the Zpt values by the amount specified (i.e. a value of 0.5 will raise existing Zpt values by 0.5m). Otherwise ignored.	Float
2	dZ	Ignored if the Z attribute is specified. If point(s) (with the Z attribute set) are snapped to the line, this attribute behaves as described in Table 7.5.	Float
3	Shape_Width_or_dMax	No TIN: Provided TIN is not specified for the Shape_Options attribute, this attribute specifies the width (thickness) of the breakline in metres. If equal to zero, only the elevations along the cell sides and corners are adjusted (i.e. a thin breakline). If less than or equal to 1.5 times the 2D Cell Size, a line of whole cells is adjusted (i.e. a thick breakline). Otherwise, for wide lines, all Zpts within a distance of half Shape_Width_or_dMax are adjusted. TIN: If TIN is specified for the Shape_Options attribute, this attribute controls the maximum distance between automatically added intermediate vertices. If set to zero, half the 2D domain’s Cell Size is used. If less than zero no intermediate vertices are inserted.	Float
4	Shape_Options	ADD: Add the shape’s Z attribute value to the current Zpts. If ADD is specified, any automatic merging around the region perimeter is ignored (to merge the perimeter with existing Zpts when using the ADD option, specify a value of zero for the Z and dZ attributes for the region). MAX, RIDGE or RAISE: Only changes a Zpt elevation if the Z Shape elevation at the Zpt is higher. MIN, GULLY or LOWER: Only changes a Zpt elevation if the Z Shape elevation at the Zpt is lower. TIN: Indicates the line is only to be used for of generation of TINS from polygons (only sections of TIN lines that fall within the polygon(s) are used).	Char(20)

Table 7.7: 2D Z Shape (2d_zsh) Polygon Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Z	If the NO MERGE option is specified, used to set the elevation of the polygon perimeter if there are no points snapped to the perimeter. If the ADD option is specified, the value entered is used to increase (positive Z values) or decrease (negative Z values) the elevation of the Zpt values by the amount specified. Otherwise ignored.	Float
2	dZ	If the Z attribute is specified and the NO MERGE option is specified, this attribute behaves as described in Table 7.5. If the Z attribute is not specified and the NO MERGE option is specified, this attribute behaves as described in Table 7.5, however Z is assumed to be zero (0). Otherwise ignored.	Float
3	Shape_Width_or_dMax	The maximum distance between intermediate points inserted around the polygon’s perimeter to interpolate elevations from the current Zpts. If set to zero, half the 2D Cell Size is used. If less than zero no intermediate vertices are inserted. Not used if the NO MERGE or ADD options are specified.	Float
4	Shape_Options	ADD: Add the shape’s Z attribute value to the current Zpts. If ADD is specified, any automatic merging around the region perimeter is ignored (to merge the perimeter with existing Zpts when using the ADD option, specify a value of zero for the Z and dZ attributes for the region). MAX, RIDGE or RAISE: Only changes a Zpt elevation if the Z Shape elevation at the Zpt is higher. MIN, GULLY or LOWER: Only changes a Zpt elevation if the Z Shape elevation at the Zpt is lower. MERGE (default): The elevations at perimeter vertices that do not have an elevation point snapped to them are merged with the current Zpt values. This option will be applied if MERGE ALL or NO MERGE are not specified. MERGE ALL: Ignores elevations from any points snapped to the perimeter and merges all perimeter vertices with the current Zpt values. NO MERGE: When there are no points snapped to the perimeter, does not merge the perimeter elevations with the current Zpt values, and sets the elevation from the Z attribute.	Char(20)

7.2.6.2.1 Order of Processing

When a mixture of different geometries and attribute options occurs within the same 2d_zsh GIS layer to control how Zpt values are modified, the order in which they are processed follows the below protocols:

1. Polygons are processed first, in the order that they occur in the 2d_zsh GIS layer.
2. Wide lines are processed second, in the order that they occur in the 2d_zsh GIS layer. If the buffer polygons of two wide lines overlap each other, the latter one prevails (i.e. the modifications do not compound). In this situation, it is recommended to separate the two wide lines into two different 2d_zsh GIS layers.
3. Thin and thick lines are processed last, depending on the Shape_Options attribute (see Table 7.6), in the following order:
   • All ADD lines are processed first.
   • Next, lines with no Shape_Options specified (i.e. NULL) – these will modify all Zpts affected by the line.
   • Then, lines with GULLY, LOWER or MIN.
   • Finally, lines with RIDGE, RAISE or MAX.
   • If there are multiple lines with the same Shape_Options, they are processed in the order that they appear in the 2d_zsh GIS layer.

The priority can be further controlled by using different 2d_zsh layers and controlling the order in which layers are listed and subsequently processed in the .tgc file (with latter occurrences of the Read GIS Z Shape command building on any previous modifications).

The same Zpts may be modified numerous times through the application of separate 2d_zsh layers or GIS objects that may overlap, with the latter occurrence prevailing based on the above processing protocols. Each sequential modification builds upon the last change. Where this occurs, the _2d_zpt_check will include multiple separate objects stacked at the same Zpt location, which detail each sequential change from the previous modified elevation assigned from all objects in all 2d_zsh layers.

7.2.6.2.2 Examples

Some examples of using Read GIS Z Shape are given below. Models set up using these topography update features are provided in the Topography Features Example Model Dataset on the TUFLOW Wiki.

Note 7.1: Example 1: Triangulating Elevations over a Null Area.

The image below shows an example of a DEM that is missing data over a small area within the 2D domain. Gaps in DEM coverage can sometimes occur over water bodies and occasionally between the tiles of ALS or LiDAR data received from a third-party provider. An example of using the MERGE topography update feature is provided in the TUFLOW Tutorial Module 2

Table 7.8: Triangulating Elevations over a Null Area

	The .tgc command Set Zpt may be used to quickly and easily assign elevations to Zpts within these areas, however the limitations of the command mean the same elevation will be assigned to all null Zpts across the entire 2D domain. This may not be suitable in situations where there are multiple gaps in coverage or where the gap is located on steep terrain. Read GIS Z Shape with 2d_zsh polygons may instead be used to triangulate Zpt values based on the Zpt elevations of the polygon perimeter.
	Import in an empty 2d_zsh GIS layer, and digitise a polygon around the gap in coverage as shown. Ensure there is a reasonable buffer around the null area. The attributes of the polygon may be left blank. Alternatively, a value may be entered in the Shape_Width_or_dMax attribute to control the maximum distance between intermediate points inserted around the polygon’s perimeter to interpolate elevations. When left blank, this distance is half the 2D cell’s size.
	This image shows the resulting _DEM_Z check file. The _zsh_zpt_check layer can be used to view the final Zpt elevation assigned.

Note 7.2: Example 2: Use of the NO MERGE and ADD Shape Options.

The 2d_zsh NO MERGE option can be used to assign a single elevation to all Zpts within the 2d_zsh polygon (read by the Read GIS Z Shape command). An example of using the NO MERGE topography update feature is provided in TUFLOW Tutorial Module 2. This may be useful to set the elevation of a polygon to a known finished floor level of a proposed development. Digitise a polygon within an empty 2d_zsh GIS layer, and populate the Z attribute of the 2d_zsh polygon with the desired elevation. Set the Shape_Options attribute to NO MERGE. The example below will assign an elevation of 45.0mAHD to all Zpts located within the 2d_zsh polygon.

Figure 7.13: 2d_zsh Example - NO MERGE

Note that if the NO MERGE option is omitted and no points are snapped to the perimeter of the polygon, the Z attribute will be ignored and the Zpt elevations will be triangulated based on the Zpt elevations of the polygon perimeter.

Alternatively, to raise the polygon by a fixed value (i.e. to represent the slab height of a building) enter this value in the Z attribute and set the Shape_Options to ADD. An example of using the ADD topography update feature is also provided in TUFLOW Tutorial Module 2.

Figure 7.14: 2d_zsh Example - ADD

TUFLOW will add the value entered in the Z attribute to the existing Zpt elevations within the polygon. The entry within the figure will raise Zpt elevations by 1m. The use of a negative value will lower the Zpt elevations by the value of the Z attribute. The _zsh_zpt_check layer can be used to view the elevation points (Zpts) that have been modified.

Note 7.3: Example 3: Raising an Embankment.

The 2d_zsh MAX option can be used to enforce the highest elevations from the 2d_zsh or the underlying topography. A simple example of when this may be desirable is enforcing an embankment across a creek to a target elevation, which merges into the banks on either side. The 2d_zsh line is used to represent the embankment.

As shown in the image below, the Z attribute contains the target elevation, the Shape_Width_or_dMax attribute has been set to the model cell size (5m) and the Shape_Options attribute is specified as MAX. As such, it is a thick breakline (refer to Section 7.2.6.2), modifying whole cells. The MAX option will adjust the underlying Zpt elevations only when the target Z shape elevation at the Zpt is higher, i.e. the underlying Zpt elevations will only be raised, never lowered. For the below example, this means the higher creek bank elevations above the level of the enforced embankment are left unadjusted.

The _zsh_zpt_check file can be used to view the final Zpt elevation assigned. This is shown in the example below.

Another example of using the MAX topography update feature is provided in TUFLOW Tutorial Module 2.

Figure 7.15: 2d_zsh Example - MAX

Alternatively, if a 3D breakline has been digitised, the dZ attribute on either the line itself, or the snapped points, may be used to raise or lower the embankment. The elevation increases or decreases the Zpt by the amount specified in the dZ attribute. This may be useful for enforcing the crest of an embankment based on surveyed data and then adjusting the elevation to represent different development scenarios.

In the example below, a 3D breakline (i.e. a line of varying elevation) has been created by snapping points to either end of the line. The elevations along the line are determined by a linear interpolation of the Z attribute of the points layer. In this example, a positive dZ attribute of 0.2 has been applied to the points layer, so the elevations along the line will raise by 0.2m. The elevations are then interpolated based on these revised values.

The _zsh_zpt_check file can be used to view the final Zpts elevation assigned (input Z + input dZ). This is shown in the example below.

Figure 7.16: 2d_zsh Example - 3D breakline

Note 7.4: Example 4: Removing ridges from a poorly triangulated DEM.

Table 7.9: Removing ridges from a poorly triangulated DEM

	The image shows two false ridges indicated by the H letters. These were caused by a poor triangulation of the TIN software used to create the DEM. These ridges caused unrealistic flow patterns, as shown by the velocity vectors. Note, the blue cells at the bottom are the downstream boundary cells.
	This image is of the DEM_Z and _zpt_check check file layers from the TUFLOW simulation. This is how TUFLOW interprets the DEM data. The false ridges are clearly shown.
	To remove the ridges, import an empty 2d_zsh layer, and digitise a polygon around the ridges as shown. By default (i.e. using the default attribute values), the elevations assigned around the perimeter of the polygon are interpolated from the current Zpt values. One problem with this approach is that the elevations along the right-hand side (i.e. along the edge of the floodplain between Locations A and B) are interpolated from the high Zpts along this boundary.
	To solve this problem, digitise points either into the 2d_zsh layer or into another 2d_zsh…_P (points) layer that snap to the vertices of the polygon where the high elevations occur. Assign a Z attribute of -99999 to each point, as shown in the image. The -99999 indicates to not interpolate an elevation from the existing Zpts. Instead, the elevations at Locations A and B are used to interpolate elevations at vertices where -99999 has been assigned.
	Use Read GIS Z Shape to process the polygon and points and generate the TIN as shown in the image. The TIN can be viewed by importing the _sh_obj_check GIS layer.
	This image shows the new DEM_Z check file after the above 2d_zsh layer has been applied. As can be seen, the ridges have been removed and the flow patterns are now realistic.

Note 7.5: Example 5: Highway Embankment Removal Example.

The figures below demonstrate a case where a new highway, present in the DEM, needed to be removed because the calibration flood events occurred before the highway was built. To remove the highway, a 2d_zsh polygon was digitised around the highway. All attributes were left as their defaults and no elevation points were required.

To reflect the topography of the original (pre-highway) terrain, a 2d_zsh line and point feature were added (using the method described in Example 3 above). Either Read GIS Z Shape or Create TIN Zpts can be used to apply this topographic modification.

Figure 7.17: 2d_zsh Example - Highway Removal

7.2.6.3 Variable Z Shape Layer (2d_vzsh)

TUFLOW 2D model topography can be varied over time to simulate breaching of embankments, raising of flood defences during an event, or the filling of a river due to a landslide, by using Read GIS Variable Z Shape. The 2d_vzsh layer is used to define the final topographic shape at the end of the topography change period. As summary of the GIS layer attributes is provided in Table 7.10. The first four attributes of 2d_vzsh are the same as the 2d_zsh layer. They are used to define the finished state of the variable geometry. Additional attributes have been added to allow the user to define how/when the breach commences and for how long. The breach/fill can be triggered using a number of methods:

• At a specified time (example provided in TUFLOW Tutorial Module 10 - Part 1);
• When the water level reaches a specified height at a specified (trigger) location (example provided in TUFLOW Tutorial Module 10 - Part 2); or
• When the water level difference between two triggers exceeds a specified amount.

The 2d_vzsh interpolates the breach height vertically between the initial and final elevations at each cell position along the breach. This occurs linearly with time by default, but can use a temporal cubic transition with “CUBIC” specified in the ‘Shape_Options’ (see Table 7.10). This cubic transition is calculated using Equation 7.5. This differs from 1D breaches (see Section 5.9.8.2) that allow specification of a rectangular or trapezoidal shape, with the opening size assumed to increase linearly with time both horizontally and vertically.

\[ z(t) = (1 - a)z_{0} + a z_1 \tag{7.5}\]

\[ a = (3 - 2x)x^2 \tag{7.6}\]

\[ x = \frac{t - t_0}{t_1 - t_0} \tag{7.7}\]

Important

Overlapping 2d_vzsh polygons are not recommended.

Variable Z Shapes can be restored once or repeatedly. Examples would be a breach of a flood defence wall or levee that is reinstated 6 hours later, or a sand bar of a creek entrance that repeatedly opens and closes. To use the restore feature, two additional attributes Restore_Interval and Restore_Period are required as described below in Table 7.10. For a single restoration event, only these two additional attributes are required. To restore repeatedly, “REPEAT” must be specified in the Shape Options in column 4. Repeated restoration is only possible for the water level and water level difference trigger methods, as a time trigger will not be able to be reached on a second occasion.

Note, this variable geometry feature should be used instead of the 2d_bc VG option (unless the rate of change of the erosion/fill is non-linear) as 2d_vzsh layers are easier to define and manage.

Table 7.10: 2D Variable Z-Shape (2d_vzsh) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Z	Point: Same as for Table 7.5. Line: Same as for Table 7.6. Polygon: Same as for Table 7.7.	Float
2	dZ	Point: Same as for Table 7.5. Line: Same as for Table 7.6. The exception is for a thin line, for which the final elevations along the line can be set to dZ metres above the current Zpt values. For example, if a thin line represents a fence that is being collapsed, and the finished height of the collapsed fence is to be 0.2m above existing ground levels, specify a dZ value of 0.2. Note this only applies if the NO MERGE Shape_Options has not been specified. Polygon: Same as for Table 7.7.	Float
3	Shape_Width_or_dMax	Point: Same as for Table 7.5. Line: Same as for Table 7.6. Polygon: Same as for Table 7.7.	Float
4	Shape_Options	Point: TRIGGER or TRIGGER 1D: Indicates the point is not an elevation point, but a trigger location. The trigger must be given a name using the Trigger_1 attribute. TRIGGER 1D is required if a 1D node water level is used to trigger a 2D variable Z-Shape. The trigger point must be snapped to the 1D node or channel end in the 2d_vzsh layer to achieve this. Line or Polygon: Same as for Table 7.6 (line) or Table 7.6 (polygon). Except the ADD option is not supported. REPEAT: Specify this option for the variable Z shape to repeatedly function indefinitely based on the below trigger and restore attributes. CUBIC: Variable Z shapes can be failed using a temporal cubic transition (see Equation 7.5). This may be desirable for some dam, levee, or embankment failures, depending on the construction material or design of the failing structure. When CUBIC is specified in the Shape_Options attribute for a 2d_vzsh line or polygon, the Z-values transition following this non-linear cubic curve. This option has no effect on point geometries (used as triggers). Thin Line: NO MERGE: For thin lines (Shape_Width_or_dMax = 0), the final elevations along the line are as specified. If NO MERGE is not specified for a thin line, the final elevations are set to be the same as the current Zpt values plus the dZ value. REPEAT: Specify this option for the variable Z shape to repeatedly function indefinitely based on the below trigger and restore attributes.	Char(20)
5	Trigger_1	Point: If Shape_Options is set to TRIGGER or TRIGGER 1D, enter the name of the trigger location. The name can contain any characters and can include spaces. Otherwise not used. Line or Polygon: To commence the failure at a specified time leave blank. To commence failure based on reaching a water level elsewhere in the model, enter the name of the trigger location. Thin Line: For thin lines there are two special options as follows. Specify DEPTH to have the failure commence once the depth of water adjacent to the cell side exceeds the amount specified for Trigger_Value. Specify DEPTH DIFF to have the failure commence once the difference in water level across the cell side exceeds the amount specified for Trigger_Value.	Char(20)
6	Trigger_2	Point: Not used. Line or Polygon: The name of a second trigger location (only needed if the breach is to be initiated on a water level difference between two trigger locations).	Char(20)
7	Trigger_Value	Point: Not used. Line or Polygon: If Trigger_1 is blank, the simulation time in hours that the breach is to commence. If Trigger_1 is specified and Trigger_2 is left blank, the water level at Trigger_1 that needs to be reached to trigger the failure. If both Trigger_1 and Trigger_2 are specified, the water level difference between Trigger_1 and Trigger_2 that needs to be exceeded to trigger the failure. The water level difference is taken as the absolute of the difference between Trigger_1 and Trigger_2, so there is no need to specify a negative value. Thin Line: If “DEPTH” is specified for Trigger_1, the depth in metres adjacent to the cell side that needs to be exceeded to trigger the failure at the cell side. If “DEPTH DIFF” is specified for Trigger_1, the water level difference in metres across the cell side that needs to be exceeded to trigger the failure. For all of the above options the length units are metres if modelling in SI units, or feet, if using Units == US Customary.	Float
8	Period	Point: Not used. Line or Polygon: Time in hours over which the variation in Zpt elevations occurs.	Float
9	Restore_Interval	The time in hours between when the variable Z shape has finished altering the geometry and when to start restoring the Zpts back to their original values. Note: “REPEAT” must be specified in Column 4 to allow repeated triggering and restoration of the Variable Z Shape (for the water level and water level difference options), otherwise restoration will only occur once.	Float
10	Restore_Period	Time in hours over which the variation in Zpt elevations occurs to restore the Zpts back to their original values.	Float

Note 7.6: Variable Z Shape Example: Breaching of an Embankment

The image below shows an example of a 2d_vzsh layer. The solid magenta line is the polygon, the magenta dashed lines are lines used to enforce TIN breaklines, and the four magenta points all have an elevation of 41.0m. The Shape_Options attribute for the polygon was set to MIN this means that the Zpt elevation can only be lowered (i.e. eroded), and the dashed lines have Shape_Options of TIN (to indicate that they are to be used for TIN generation, and not for Z lines). The vertices of the polygon that do not have a point snapped to them will be automatically assigned an elevation based on the existing Zpt values. The polygon vertices with the points snapped to them are assigned the elevation of the point (in this case, all at 41.0m). The elevation of the dashed lines will be constant at 41m as they are snapped to the 41m points.

The only other object in the layer is the yellow pin point labelled A. This is a trigger point. Its only attribute values are: TRIGGER for the Shape_Options attribute; and A for the Trigger_1 attribute. This sets the point as a trigger point and the “A” is the name of the trigger point. The magenta polygon also has attribute values of: A for Trigger_1 (this indicates that the erosion trigger is based on the water level at Trigger Location A); a Trigger_Value of 42.0 (i.e. when the water level at location A reaches 42m, start the erosion); and a Period value of 1.0 indicating that the erosion takes one hour to complete.

Figure 7.18: Example 2d_vzsh layer

The final eroded Zpt values are based on the TIN created by TUFLOW (the grey triangles in the image below). The central section will be horizontal at 41m, sloping up either end to elevations based on the road level. The _vzsh_zpt_check layer is useful to view the Zpts affected by the variable Z shape. This layer is also shown in the image below. The green triangles indicate that the Zpt level is to be eroded, and the crosses indicate no change (this is because of the MIN Shape_Options). The final eroded Zpt values are labelled in the image below. Other useful attributes are also available in this layer.

Figure 7.19: Example _vzsh_zpt_check layer

The images below show the modelled breach which occurs using the example above. Each image is in half hour intervals. The colour shading is of the elevations (specify ZH as a Map Output Data Types to view the changes in ground level over time).

Figure 7.20: Modelled breach using 2d_vzsh layer

7.2.6.4 3D TIN Layers (2d_ztin)

TINs (triangulations) of elevation points and 3D lines within a polygon can be carried out using Create TIN Zpts. This is particularly useful for modifying the Zpt elevations where there have been, or are proposed, changes to the base DEM Zpt values.

Note, for large datasets it is likely to be much more efficient to use GIS or 3D surface modelling software to triangulate the data, and read triangulated data in a supported format (see Read TIN Zpts) or to convert to raster (see Read Grid Zpts).

The protocols applied to the Create TIN Zpts command are:

1. A TIN is created for each polygon in the 2d_ztin layer.
2. Any points found within a polygon are used when generating the TIN.
3. Any lines are converted to points, and those points falling within the polygon are used for the TIN creation. Lines are converted to points as follows:
   1. All vertices (nodes) of the line are converted to points.
   2. The dMax attribute is used to insert additional vertices between the line’s vertices. For example, if dMax is set to 10, then additional intermediate vertices are inserted at least every 10 metres (or feet if using Units == US Customary) between the existing vertices where the distance between the existing vertices exceeds this value. If the dMax attribute does not exist or is zero, half the 2D domain’s Cell Size is used as the dMax value.
   3. If there are any points snapped to the line’s vertices, the elevations of these points are used to set the elevations at all the vertices generated along the line. In this way, a 3D breakline effect can be produced within the TIN. If there are no points snapped to the line, the line’s Z attribute elevation is used giving the effect of a horizontal line.
4. The perimeter of the polygon/TIN can either be merged with the current Zpt values or have its own values as follows:
   1. If there are no points snapped to the perimeter of the polygon, the elevations of the polygon’s perimeter vertices, and of any automatically inserted vertices, are based on the current Zpt values (i.e. the Zpt values assigned by any prior commands).
   2. If there are one or more points snapped to the polygon’s perimeter vertices, the perimeter is not merged with the Zpt values, and the elevations of the snapped points are used to assign elevations to the perimeter vertices and any automatically inserted vertices.
   3. The frequency of any automatically inserted points around the perimeter is controlled by the dMax attribute. If the dMax attribute does not exist or is zero, half the 2D domain’s Cell Size is used.

The 2025 release introduced a new method “Method B” as beta functionality for triangulating TINs that applies to any of the commands that build a TIN, such as Create TIN Zpts and Read GIS Z Shape commands. The new method can be applied using the TIN Triangulation Approach command. “Method A” is the original triangulation approach. Method B uses the same inputs as Method A and creates the same input points and elevations, however, creates the triangles differently. While Method A adds additional vertices to lines to encourage edges to follow breaklines, these are strongly enforced using Method B. Method B also generates more balanced triangles avoiding long skinny triangles or clumps of triangles.

Models using the Create TIN Zpts functionality are provided in the Static Topography Updates Example Model Dataset on the TUFLOW Wiki.

A useful quality control option of Create TIN Zpts is the WRITE TIN option. This will write a TIN to either the SMS TIN or SMS 2DM format depending on the setting of TIN Output Format. If this option is specified, a SMS .tin or .2dm file is written for each TIN generated, and the triangles are written to the _sh_obj_check layer. This means that the TIN can be cross-checked, viewed, and edited and modified in SMS (tin or 2dm) or QGIS (2dm). Read TIN Zpts can be used to assign Zpt elevations from the modified SMS TIN.

A second argument to specify a GIS layer containing one or more polygons to clip the area of Zpts to be inspected can be used with the Read TIN Zpts command. Refer to Section 7.2.6.1 for more information.

Table 7.11: 2D Tin (2d_ztin) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Z	Point: Elevation of the point. Line: Elevation of the line. Ignored if there are any points snapped to the line’s vertices. Polygon: Not used.	Float
2	dMax (optional)	Point: Not used. Line: Maximum distance between automatically created intermediate vertices. If set to zero or this attribute does not exist, half the 2D domain’s Cell Size is used. If less than zero no intermediate vertices are inserted. Polygon: Same as for Line above.	Float

7.2.6.5 Z Line Layers (2d_zln)

This is a legacy feature and it is recommended to use the Z Shape functionality (Section 7.2.6.2) instead. For details on this feature see Read GIS Z Line in Appendix C.

7.2.6.6 Zpt Layers (2d_zpt)

Zpt (2d_zpt) layers are associated with the following legacy commands, which only remain for version backward compatibility:

• Read RowCol Zpts - for setting base elevations
• Read GIS Zpts - for simple topographic modifications

It is recommended to instead use Read Grid Zpts (Section 7.2.6.1) for setting base elevations and Read GIS Z Shape with 2d_zsh layers (Section 7.2.6.2) for modifying base topography.

7.2.6.7 Using Multiple Layers and Points Layers

GIS layers used for Read GIS Z Shape, Read GIS Variable Z Shape, Create TIN Zpts, Read GIS Layered FC Shape, Read GIS Z Line (legacy), Read GIS Z HX Line (legacy) and Read GIS FC Shape (legacy) can be split into more than one layer to better manage the variety of data these commands sometimes require.

For example, one layer may contain the elevation points, another the TIN lines and polygons and another the 3D Z lines. This is useful in terms of managing the data, and especially when interrogating and/or viewing the data in GIS. It is a requirement of the shapefile format that the different geometries (points, lines and regions) are in separate shapefiles. The TUFLOW empty template files include the following filename suffixes to differentiate which files are suitable for point, line or region features.

• _P for point features (e.g. 2d_zsh_M03_002_P.shp)
• _L for line features (e.g. 2d_zsh_M03_002_L.shp)
• _R for region or polygon features (e.g. 2d_zsh_M03_002_R.shp)

This is optional for MapInfo users; the different geometries can occur in the same MapInfo file or can be separated if preferred.

A maximum of nine (9) layers per command line is allowed, and each layer is separated by a vertical bar (“|”). For example, to read a Z Shape layer which has both line and points, the command may be:

Read GIS Z Shape == gis\2d_zsh_M03_002_L.shp | gis\2d_zsh_M03_002_P.shp

If using the GeoPackage format && is also used to specify more than one layer from a common database in the same command line. An example is provided below. See Section 4.3.3 for further details.

Read GIS Z Shape == gis\2d_zsh_R.gpkg | gis\2d_zsh_002.gpkg >> 2d_zsh_L && 2d_zsh_P

7.2.6.7.1 Point Only Layers

To minimise the number of attributes, some/all points may optionally be placed into a separate layer with less attributes as discussed below. This simplifies the datasets making them easier to manage and interrogate.

A layer is treated as a separate points layer if:

1. It has less attributes than the minimum required for the command. For most commands there is only one attribute for the points layer (i.e. Elevation or Z) as described in Table 7.12. The exception is Read GIS FC Shape, which requires the first two attributes. This option requires that the points layer be defined within the command line syntax. For example: Read GIS Z Shape == gis\2d_zsh_M03_002_L.shp | gis\2d_zsh_M03_002_P.shp
2. The points file uses the same filename as the associated line or region file with the addition of “_pts” as a suffix to its filename (for example 2d_zsh_M03_002_pts.shp will be automatically associated with the line file 2d_zsh_M03_002.shp). This option is supported for backward compatibility; however, it’s recommended that this option not be used (it is preferable to enter the filename of the second layer so that it is clear as to which layers are being used). Read GIS Z Shape == gis\2d_zsh_M03_002.shp

The data processing logic for points layers is outlined below:

1. The specified layer, 2d_zsh_M03_002.shp, is opened. This layer may or may not contain elevation points. If any elevation points exist they are used.
2. A separate points layer can optionally be used to specify additional points or all of the points. The layer can be specified in one of two ways:
   1. Entering the pathname of the points layer after the main layer. A “|” must be used to separate the two layers. The points layer must be the second layer specified. For example: Read GIS Z Shape == gis\2d_zsh_M03_002_L.shp | gis\2d_zsh_M03_002_P.shp
   2. Alternatively, name the points layer the same as the main layer, but with a “_pts” extension. If a layer exists with the “_pts” extension, TUFLOW automatically assumes this layer is associated with the main layer and includes all points within this layer when applying the above commands. For this example, the layer would be named 2d_zsh_M03_002_L_pts.shp.
3. The first approach (i) above prevails over the second (ii) if both apply.
4. If neither (i) or (ii) apply, TUFLOW assumes there is no separate points layer.

Multiple points layers can be specified. The points layer can be referenced in any location except for the first layer within the command line entry. For example, the below syntax will produce an error due to the points file being the first entry.

Incorrect:

Read GIS Z Shape == gis\2d_zsh_M03_002_P.shp | gis\2d_zsh_M03_002_L.shp

Correct:

Read GIS Z Shape == gis\2d_zsh_M03_002_L.shp | gis\2d_zsh_M03_002_P.shp

Table 7.12: 2D Z (2d_z_) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Z	Elevation (or change in elevation for ADD option) of the point.	Float

7.2.7 Land Use (Materials)

7.2.7.1 Bed Resistance

The default bed resistance formulation is the Manning’s equation, where bed resistance is calculated via Equation 7.8:

\[ \tau_{bed} = \frac{\rho g V^2 n^2}{y^{\frac{1}{3}}} \qquad (N/m^2) \tag{7.8}\]

or, if using English Units (Units == ENGLISH):

\[ \tau_{bed} = \frac{\rho g V^2 n^2}{2.208 y^{\frac{1}{3}}} \qquad (lbf/ft^2) \tag{7.9}\]

Where:

• \(\rho\) = density of water
• \(g\) = acceleration due to gravity
• \(V\) = depth averaged velocity
• \(n\) = Manning’s bed friction coefficient
• \(y\) = depth of water

Manning’s n values can be varied with depth (as user specified curve or using the Log Law formula (see Section 7.2.7.3.1)) or varied with velocity-depth product (VxD).

For TUFLOW Classic, bed resistance can also be set to use either Manning’s M values (1/n) or Chezy coefficients using the Bed Resistance Values command in the .tcf file. As this is a TUFLOW Classic feature only, it is further discussed in Section 7.4.3.

The bed resistance values for 2D domains are created using GIS layer polygons or rasters of different bed resistance zones. The recommended approach is to use materials to define how the bed roughness varies over the model. Each material is defined by a positive integer ID which represents a different roughness category. GIS layers of land-use or vegetation often make excellent material layers. Examples of different material categories are river in-bank, bank vegetation, pasture, maintained grass, roads, buildings, forest, mangroves, etc. Each material is assigned a constant Manning’s n value, depth or VxD varying Manning’s n. The material layer can also be used to set rainfall losses (if using direct rainfall - see Section 8.4.3), fraction impervious, storage area and land-use hazard categories (see Section 7.2.7.3.1).

Material/roughness values are used by TUFLOW during conveyance calculations at the cell mid-sides (refer to Section 7.1 and Section 7.1.3). However, rainfall losses and fraction impervious are applied to the cell and not cell sides, and materials ID values are sampled at both cell mid-sides and cell centres.

The most common approach is to digitise one or more 2d_mat materials layers (see Table 7.13) and assign Manning’s n values to the materials using Read Materials File. This approach allows the easy adjustment of Manning’s n values, for example during model calibration or sensitivity testing.

When creating the base 2d_mat layer, it is common practice to not digitise the most common or the most difficult to digitise material, and instead use the following data layering of commands in the .tgc file (see Section 4.1.8).

• Use Set Mat to set the most common material to all cells in a 2D domain.
• Use Read GIS Mat or Read Grid Mat to allocate the remaining material values.

The Read GIS Mat and Read Grid Mat commands may be used as many times as required to further modify the materials in parts of a 2D domain. Each subsequent dataset will overwrite the preceding assigned material value, as described in Section 7.2.3.

The default material value is zero. As a material value of zero is not allowed, every cell and cell-side must be assigned a material value using Set Mat, Read GIS Mat and/or Read Grid Mat in the .tgc file (it is good practice to always set a default materials value using the Set Mat as the first material command in the .tgc file). The assigned material ID values do not need to be continuous but must be within the range 1 to 32,767.

Table 7.13: 2D Materials (2d_mat) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Material	The material ID value referenced within a Materials File (see Section 7.2.7.3).	Integer

7.2.7.2 Log Law Depth Varying Bed Resistance

At very shallow depths the Manning’s n value and/or equation may not be a reliable estimate of bed resistance. The Log Law or “Law of the Wall” approach offers a theoretically based derivation of resistance based on a bed shear analysis. This relationship along with benchmarking against flume test results was used by Boyte (2014) to derive the following equation that varies Manning’s n with depth based on the roughness height of the surface:

\[ n = \frac{\kappa y^{\frac{1}{6}}}{\sqrt{g}\left[\ln\left( \frac{y}{z_{0}} \right) - 1 \right]} \tag{7.10}\]

\[ z_{0} = \ \frac{k_{s}}{30} + \frac{0.11\nu}{U_{f}} \tag{7.11}\]

Where:

• \(k_s\) is the roughness height in m
• \(\kappa\) is typically in the range 0.38 to 0.42 (recommend 0.4)
• \(y\) is depth
• \(\nu\) is the kinematic viscosity and is set to 10\(^-6\) m\(^2\)/s
• \(U_f\) is the friction velocity defined as \(\sqrt{Sgy}\) where \(y\) approximates \(A/P\) and \(S\) is the water surface slope

Equation 7.10 becomes singular as \(\ln\left( \frac{y}{z_{0}} \right)\) approaches 1, and asymptotes to zero as depth becomes large. It is necessary to apply sensible limits for the resulting \(n\) under these circumstances. For depths such that \(\ln\left( \frac{y}{z_{0}} \right) < 1.000001\) the Manning’s n is set to \(n = 10\). For all other cases, the resulting value for \(n\) is bounded to be in range \(n_{limit} < n < 10.\)

Figure 7.21 illustrates how the equivalent Manning’s n varies with depth using the log law for a roughness height of 10mm (0.01m) that would be applicable to a small pebble bed. The different series are the variations in the slope, \(S\), where 0.001 is 0.1% slope, 0.02 is 2% slope and 0.1 is a 10% slope. As can be seen there is a significant variation in Manning’s n below 2cm (0.02m) and a trend to a n value of around 0.018, with only a minor variation due to slope.

In terms of applying a limiting n value, if, for example, \(n_{limit}\) was set to 0.02, then the Manning’s n value would not fall below 0.02.

Figure 7.22 shows a comparison using the Log Law versus a constant Manning’s n value (Boyte, 2014). The thesis investigated the use of the Law of the Wall for direct rainfall modelling using TUFLOW. The flow depths in this example range from 4 to 20cm and the roughness height, \(k_s\), was 3cm.

The Log Law Depth Varying Bed Resistance is activated by entering special characters in the 2nd column of the Materials File. See Section 7.2.7.3 and Table 7.14.

Figure 7.21: Example of Log Law Variation of Manning’s n with Depth

Figure 7.22: Example of Log Law versus Constant Manning’s n with Depth

7.2.7.3 Materials File

Materials file(s) contain information on a material’s roughness and, optionally, rainfall losses if using direct rainfall. The file is referenced within the .tcf file using Read Materials File and can be in one of two formats (.csv or .tmf). The .csv format is the recommended of the two options. It supports all functionality. The .tmf format does not. For example, the log law bed resistance option is only available via the .csv format. The .csv format also supports curves of Manning’s n versus depth.

More than one materials file may be specified by repeat occurrences of the command Read Materials File however, most models will use only a single materials file. Any combination of .tmf and .csv files can be used and up to 1,000 materials are allowed in total.

If a second argument is provided with Read Materials File, this value is used to factor all Manning’s n values. For example, the following command increases all Manning’s n values by 10%:

Read Materials File == My_Materials.tmf | 1.1

7.2.7.3.1 .csv Format

The .csv format materials file is a comma delimited text file containing Manning’s n and other information for different materials (e.g. land-uses). The format is intended to be generated from an Excel file database of materials and associated data in a similar manner to BC databases (with the option of using the Excel TUFLOW Macro .xlam macros to export to the .csv format - the Excel TUFLOW Tools.xlam can be downloaded from here). The .csv can also be written from text editor if preferred.

The format of the materials.csv file is described in Table 7.14.

Tip

The .csv format offers access to all materials features, whereas the .tmf format does not.

Table 7.14: Materials .csv File Format

No.	Description
1	Mat (Material ID) number, which must be an integer.
2	Contains information on the bed resistance values (usually Manning’s n). The options available are: Single n value (the n value is assumed constant at all water depths). Two pairs of depth and n values (see the y\(_1\), n\(_1\), y\(_2\), n\(_2\) option described in Table 7.15 for the .tmf format (i.e. Columns 5 to 8). The values must be separated by commas. A .csv file containing n vs depth values. The data is linearly interpolated between values. The first and last n values are applied as constants for depths below the first depth value and above the last depth value, respectively. TUFLOW Classic is limited to 100 data points, while TUFLOW HPC is limited to 8. By default, when the n vs depth values are read they are assumed to be in the first two columns of the .csv file containing the curve’s data. However, you can specify the columns to be used by placing the column label after the .csv filename using a vertical bar | to separate the column label from the .csv filename. If only one column label is given, this is assumed to be the depth column and the next column must contain the Manning’s n values. As for boundary data, TUFLOW searches down the rows until it finds the first numeric values and will start reading from this row until a row is found with no more numeric values. To view how the n values vary with depth over time, use the “n” option for Map Output Data Types. If the first four characters are “VxD:” the Manning’s n values for the material will be varied according the velocity times depth product, rather than depth. The VxD versus n values are provided in the same manner as for depths described above. Note that the Manning’s n value never increases during the simulation, it can only remain the same or decrease. This feature was provided to allow the user to reduce Manning’s n with increasing VxD to represent the “stripping” of river banks in high VxD areas, and then keeping the reduced Manning’s n as the flood recedes, hence why the n value cannot increase once decreased during the simulation. VxD was used in preference for simply V to avoid reducing n values at very shallow depths (but high velocities) when the flood first inundates the 2D cells. If the first four characters are “Log:” the bed resistance as defined by the theoretically derived log law for shallow flows is applied. Three numbers are required after the “Log:” in space or comma delimited format. These are: Ks (roughness height), Kappa (0.3 to 0.4) and the limiting Manning’s n value. See Section 7.2.7.2 for more information.
3	Sets the rainfall loss parameters using the initial loss/continuing loss option. The initial/continuing loss is entered as two comma delimited numbers in a similar manner to the third and fourth column values in the .tmf format. See Table 7.15. Refer to Section 7.2.7.4.
4	Reserved.
5	Defines the Storage Reduction Factor (SRF) value. If no fifth column entry exists, no SRF is applied. The default is an SRF of 0 (i.e. no change in storage). See Section 7.2.10.1 for more information.
6	Defines the Fraction Impervious of the overlying material type. The value entered should be a number from 0.0 to 1.0 where 0.0 is fully pervious and 1.0 is fully impervious. The default is a value of 0.0, assuming that the overlying material is 100% pervious. This feature is used to influence the amount of water that is infiltrated into the ground with the soil infiltration feature. Refer to Section 7.2.8 for more information. Note: This option works with the Soli Infiltration feature (see Section 7.2.8). It does not apply to materials rainfall losses (Column No 3 above) when applying direct rainfall (see Section 8.4.3).

An example of a materials.csv is provided in Figure 7.23. To give a description of the material, this must be done after all inputs for that material and must be preceded by a “!” or “#”.

In the example shown in Figure 7.23:

• Material 1 has a constant n value of 0.03 and no rainfall loss parameters.
• Material 5 varies n with depth using the four y\(_1\), n\(_1\), y\(_2\), n\(_2\) values (as per .tmf format approach) and no rainfall loss parameters.
• Materials 11 and 12 use depth vaying n sourced from a file called Grass.csv. As only one column label has been specified, the y values must occur under that label and the n values must occur in the next adjoining column (see Figure 7.24). An IL of 10mm and CL of 2mm/h for both materials will be used for any direct rainfall.
• Materials 21 and 22 use depth vaying n sourced from a file called Trees.csv. As two column labels have been specified, the y and n values must occur under the specified labels (see Figure 7.24). An IL of 20mm and CL of 2mm/h for Material 21, and 25 and 2 for Material 22, will be used.
• Material 31 uses the log law approach (Section 7.2.7.2) with a roughness height of 0.01m, Kappa value of 0.4 and limiting Manning’s n of 0.02. IL = 0mm and CL = 0mm. The fraction impervious value is set to 1.0 (fully impervious).
• Material 32 use depth vaying n sourced from a file called Buildings.csv. In Buildings.csv Column A is used for y values and Column B for n values. IL = 0 and CL = 0. The storage reduction factor is set to 0.3 (30% reduction in storage) and the fraction impervious value is set to 1 (fully impervious).

Figure 7.23: Example of Materials .csv File Format

Figure 7.24: Example of the Grass.csv file read into the Materials.csv

Figure 7.25: Example of the Trees.csv file read into the Materials.csv

7.2.7.3.2 .tmf Format

The .tmf format is a text file containing Manning’s n and other information for different materials (e.g. land-uses). The file can contain comments using the “#” and/or “!” comment characters at any location. When comment characters are specified, the remainder of the line is ignored. The format of the materials .tmf file is described in Table 7.15. The first two columns are mandatory and must be specified. All other columns are optional. A maximum of 100 different materials can be specified in this format.

Tip

The .tmf format does not offer all materials functionality, whereas the .csv format in the above section does.

Table 7.15: Materials .tmf File Format

No.	Description
1	Mat (Material ID) number, which must be an integer.
2	Manning’s n value. Note, if the four values in columns 5 to 8 are specified, the Manning’s n value in this column is ignored and not used.
3	Sets the initial loss if using a direct rainfall boundary. Refer to Section 7.2.7.4 and Section 8.4.3. The units are mm, or if Units == US Customary inches.
4	Sets the continuing loss rate if using a direct rainfall boundary. Refer to Section 7.2.7.4 and Section 8.4.3. The units are mm/hr, or if Units == US Customary inches/hr.
5	y\(_1\) - The depth below which the Manning’s n value n1 (column 6) is applied.
6	n\(_1\) -The Manning’s n value applied below depth y\(_1\) (column 5).
7	y\(_2\) - The depth above which the Manning’s n value n\(_2\) (column 8) is applied.
8	n\(_2\) -The Manning’s n value applied above depth y\(_2\) (column 7). Between y\(_1\) and y\(_2\), the Manning’s n value is interpolated between n\(_1\) and n\(_2\) according to Bed Resistance Depth Interpolation. When specifying values for columns 5 to 8, initial and continuing loss values must be specified in columns 3 and 4 as described above (use zero values if not using direct rainfall).
9	Reserved.
10	Defines the Storage Reduction Factor (SRF) value. The default is an SRF of 0. Enter 0 to ensure there is no change in 2D cell storage for the material type. See Section 7.2.10.1 for more information.
11	Defines the fraction impervious of the overlying material type. The value entered should be a number from 0.0 to 1.0 where 0.0 is fully pervious and 1.0 is fully impervious. The default is a value of 0.0, assuming that the overlying material is 100% pervious. This feature is used to influence the amount of water that is infiltrated into the ground with the soil infiltration feature. Refer to Section 7.2.8 for more information. Note: This option works with the Soli Infiltration feature (see Section 7.2.8). It does not apply to materials rainfall losses (Column No 3 above) when applying direct rainfall (see Section 8.4.3).

The .tmf file format is shown in the examples below. See Set Mat, Read GIS Mat and Read Grid Mat for assigning the material IDs to the 2D domains. These material values may also be used to define bed resistance values across 1D XZ cross-sections (see Section 5.6.1.1.2).

! Comments and blank lines are allowed in this file
! First value is the Mat value, Second is the Manning's n value
1, 0.03      ! waterways
2, 0.08      ! river banks
11, 0.06     ! grazing land
12, 0.04     ! parks and gardens
13, 0.15     ! sugar cane
14, 0.12     ! natural forest
15, 0.02     ! roads

To include the initial loss (mm) and the continuing loss rate (mm/h) optionally enter a third and fourth value as shown below. If an IL is specified, a CL must also be specified otherwise an ERROR occurs. Both can be omitted, in which case, they are both set to zero.

1, 0.03             ! waterways
2, 0.08, 20, 2      ! river banks
11, 0.06, 20, 2     ! grazing land
12, 0.04, 5, 1.5    ! parks and gardens
13, 0.15, 10, 2     ! sugar cane
14, 0.12, 30, 2.5   ! natural forest
15, 0.02, 2, 0      ! roads

To vary n values with depth (m) using two pairs of values optionally enter values in the fifth to eighth columns as shown below in lines 2 and 3 (Materials 2 and 11) below. IL and CL values must be entered (use zero if not relevant). If depth varying manning’s n is not used, column 5 to 8 should be left blank, not set to zero.

1, 0.03                                 ! waterways
2, 0.08, 20, 2, 0.3, 0.15, 0.5, 0.08    ! river banks with long grass 0.3m high
11, 0.06, 20, 2, 0.1, 0.1, 0.15, 0.06   ! grazing land
12, 0.04, 5, 1.5                        ! parks and gardens
13, 0.15, 10, 2                         ! sugar cane
14, 0.12, 30, 2.5                       ! natural forest
15, 0.02, 2, 0                          ! roads

To reduce the storage of cells, enter a SRF value in the tenth column as shown in line 7 for Material 15 below. The example shown reduces the storage of all cells with Material ID of 15 by 20%. Note that a value of -1 has been entered for the ninth column as this field has been reserved. If not used, specify 0 to ensure no change in storage for the material.

1, 0.03                                 ! waterways
2, 0.08, 20, 2, 0.3, 0.15, 0.5, 0.08    ! river banks with long grass 0.3m high
11, 0.06, 20, 2, 0.1, 0.1, 0.15, 0.06   ! grazing land
12, 0.04, 5, 1.5                        ! parks and gardens
13, 0.15, 10, 2                         ! sugar cane
14, 0.12, 30, 2.5                       ! natural forest
15, 0.02, 2, 0, , , , , -1, 0.2         ! roads with width reduced by 20%

To specify a fraction impervious, enter a value between 0.0 and 1.0 in the eleventh column as shown in line 4 for Material 12 below. The example shown partially restricts the rate of infiltration by applying a fraction impervious of 0.1 or 10%.

1, 0.03                                 ! waterways
2, 0.08, 20, 2, 0.3, 0.15, 0.5, 0.08    ! river banks with long grass 0.3m high
11, 0.06, 20, 2, 0.1, 0.1, 0.15, 0.06   ! grazing land
12, 0.04, 5, 1.5, , , , , -1, 0, 0.1    ! parks and gardens (10% impervious)
13, 0.15, 10, 2                         ! sugar cane
14, 0.12, 30, 2.5                       ! natural forest
15, 0.02, 2, 0, , , , , -1, 0.2         ! roads with width reduced by 20%

7.2.7.4 Rainfall Losses

Rainfall losses applied through the Materials file (Section 7.2.7.3) remove the loss depth from the rainfall before it is applied as a boundary on the 2D cells (Section 8.4.3). Rainfall losses are ideal for modelling situations in which water is prevented from reaching the ground, such as through the interception by trees.

Important

The ILCL rainfall losses are different to the ILCL infiltration losses that can be applied using the .tsoilf file (refer to Section 7.2.8.2).

The ILCL soil infiltration losses will infiltrate ponded water into the ground. It is possible to use both methods in the same simulation – for example, rainfall that doesn’t reach the ground would be modelled as a material IL/CL (applied as a loss to the rainfall) and infiltration into the ground as IL/CL via soil types. The log file (see Section 14.4.1) will report the material and soil properties separately as shown in the example below:

Example Material Properties:

#4 - Material 4:
Fixed Manning’s n = 0.030
IL = 1.0mm, CL = 0.0mm/h
Landuse Hazard ID not set.
SRF (Storage Reduction Factor) = 0.
Fraction Impervious = 0.

Example Soil Properties:

#1 - Soil 1 [based on pre-defined soil type SAND]:
Suction = 49.5 mm
HydCond = 117.8 mm/hr
Porosity = 0.417
Initial Moisture = 0.2
Soil Capacity = 0.217

Specifying the “fraction impervious” on the material allows the materials and the soils to be independent (i.e. the same soil can be present under both road and forest). This fraction impervious only applies to the infiltration into the soil, not to the Materials file rainfall losses.

7.2.8 Infiltration

Four methods are available to infiltrate water from the surface into the ground:

• Green-Ampt (GA) - Section 7.2.8.1.1;
• Horton (HO) - Section 7.2.8.1.2;
• Initial Loss/Continuing Loss (ILCL) - Section 7.2.8.1.3; and
• SCS Curve Number (SCS) (HPC only) - see Section 7.3.5.1.1.

When using TUFLOW HPC, it is possible to model multiple vertical groundwater layers, as well as enable horizontal flow of groundwater. As these features are only available in TUFLOW HPC, they are discussed in Section 7.3.5.2.

The first three methods mentioned above (GA, HO, and ILCL) monitor the amount of water infiltrated and if the ground becomes saturated, infiltration ceases. With these methods the amount of water that can be infiltrated depends on:

• The infiltration approach and soil parameters used for the top most (infiltration) soil layer.
• The number, thickness, porosities, and initial moisture of the groundwater layers.
• The fraction impervious value of the overlying material layer.

The cumulative infiltration data is stored and computed at cell centres. As such all soil data and groundwater layer geometry is sampled from the GIS / Grid layers at cell centres. Likewise the material “fraction impervious” data is sampled at cell centres (unlike the Manning’s values which are sampled at face centres).

The following .tgc commands are used to set the ID for the soil type referenced in the TUFLOW soils file (.tsoilf) (Section 7.2.8.2):

• Set Soil;
• Read GIS Soil; and
• Read Grid Soil.

Each of these commands has the option of specifying which soil layer(s) the command applies to. Note that 2D infiltration is activated by the occurrence of one of these commands. If none of the commands exists for a 2D domain, soil infiltration does not occur for that domain.

If soils are specified, the soil ID for each cell are written to the _grd_check layer. The .tlf file contains the parameters for each Soil ID.

7.2.8.1 Infiltration Layer

Only wet 2D cells can infiltrate water into the ground. If using soil infiltration with a single groundwater layer, the default groundwater depth is infinite. The depth (thickness) of the groundwater layer can be set globally using the following commands:

• Set Soil Thickness
• Set Soil Base Elevation

They can also be set spatially with the following GIS and Grid commands:

• Read GIS Soil Thickness
• Read GIS Soil Base Elevation
• Read Grid Soil Thickness
• Read Grid Soil Base Elevation

The soil thickness sets the layer depth from the surface elevation. The soil base elevation sets the absolute elevation of the bottom of the layer. If both methods are specified for a given grid cell, the highest of the two will be adopted. The input units should be in metres or feet. To set an initial groundwater level within the soil layer, see Section 8.8.2.

In releases prior to 2023-03, the groundwater level (or depth) was set to represent a water table. For backward compatibility, it is still possible to the set a water table, using these commands:

• Set GWD or Set GWL
• Read GIS GWD or Read GIS GWL
• Read Grid GWD or Read Grid GWL

The amount of water that enters the soil is also dependent on the fraction impervious value of the overlying material layer. The default is that the overlying material is 100% pervious (i.e. 0% impervious). However, if, for example, a concrete parking lot overlies a sandy soil, the imperviousness of the parking lot can be specified as 100% (i.e. totally), or 90% (i.e. partially) restrict the rate of infiltration. This is described in the materials file in Section 7.2.7.

If using TUFLOW HPC’s SGS functionality (Section 7.3.3), cells can be considered “partially wet”. This necessitated options for factoring down the rate of infiltration into the topmost groundwater layer, the methods are discussed in Section 7.3.3.1.3.

7.2.8.1.1 Green-Ampt (GA)

The Green-Ampt approach varies the rate of infiltration over time based on the soil’s hydraulic conductivity, suction, porosity and initial moisture content. The method assumes that as water begins to infiltrate into the soil, a line develops differentiating between the “dry” soil (with moisture content\(\ \theta_{i}\)) and the “wet” soil (with moisture content equal to the porosity of the soil η). As the infiltrated water continues to move through the soil profile in a vertical direction, the soil moisture changes instantly from the initial content to a saturated state. This concept is presented in Figure 7.26.

Figure 7.26: Green-Ampt Model Concept¹

The basic form of the Green-Ampt equation is expressed as follows:

\[ f(t) = K\left( 1 + \ \frac{\mathrm{\Delta}\theta\left( \varphi + h_{0} \right)}{F(t)} \right) \tag{7.12}\]

Where:

• \(t\) = time
• \(K\) = saturated hydraulic conductivity
• \(\mathrm{\Delta}\theta\) = defined as the soil capacity (the difference between the saturated and initial moisture content)
• \(\varphi\) = soil suction head
• \(h_{0}\) = depth of ponded water
• \(F(t)\) = cumulative infiltration calculated from:

\[ F(t) - \ \mathrm{\Delta}\theta\left( \varphi + h_{0} \right)\ln\left( 1\ + \ \frac{F(t)}{\mathrm{\Delta}\theta\left( \varphi + h_{0} \right)} \right) = Kt \tag{7.13}\]

The TUFLOW Wiki Green-Ampt Infiltration Parameters provides a detail description of the various Green-Ampt parameters, including how each influences the soil infiltration behaviour.

Models using the Green-Ampt infiltration loss approach are provided in the Loss Options Example Model Dataset on the TUFLOW Wiki.

United States Department of Agriculture (USDA) soil types have been hardwired into TUFLOW and are presented in Table 7.16, along with the soil parameters. Alternatively, it is possible to define a customised soil type by specifying user defined values as shown in Table 7.17.

Table 7.16: USDA Soil types for the Green-Ampt Infiltration Method

USDA Soil Type	Suction (mm)	Suction (inches)	Hydraulic Conductivity (mm/hr)	Hydraulic Conductivity (in/hr)	Porosity (Fraction)
Clay	316.3	12.453	0.3	0.012	0.385
Silty Clay	292.2	11.504	0.5	0.02	0.423
Sandy Clay	239.0	9.409	0.6	0.024	0.321
Clay Loam	208.8	8.22	1.0	0.039	0.309
Silty Clay Loam	273.0	10.748	1.0	0.039	0.432
Sandy Clay Loam	218.5	8.602	1.5	0.059	0.33
Silt Loam	166.8	6.567	3.4	0.134	0.486
Loam	88.9	3.5	7.6	0.299	0.434
Sandy Loam	110.1	4.335	10.9	0.429	0.412
Loamy Sand	61.3	2.413	29.9	1.177	0.401
Sand	49.5	1.949	117.8	4.638	0.417

The Green-Ampt approach is appropriate for simulating single rainfall events where the evapotranspiration and gravity-driven drainage is not significant. The 2023-03 release introduced functionality in TUFLOW HPC to allow for horizontal movement of soil water (Section 7.3.5.2.1), to extend TUFLOW’s capacity to conduct long term simulations with multiple rainfall events. This means, ideally, the initial soil moisture and the cumulative infiltration must be ‘reset’ before each rainfall event. The 2023-03 release implemented a simple change to the Green-Ampt equation to consider the changing initial soil moisture and the cumulative infiltration.

\[ f(t) = K\left( 1 + \ \frac{\phi\left( \varphi + h_{0} \right)}{h_{soil}} \right) \tag{7.14}\]

Where:

• \(\phi\) = porosity of the soil
• \(h_{soil}\) = Soil water depth.

The initial moisture is no longer used to set \(\mathrm{\Delta}\theta\) (the difference between the saturated and initial moisture content) in Equation 7.12, and the cumulative infiltration is replaced by the soil water depth at the time rainfall event.

These two approaches produce slightly different results for the infiltration rate with the Green-Ampt method. If the previous approach is required, and the model has only one soil layer and no horizontal advection, then it can be selected by using the Defaults == Pre 2023 command.

7.2.8.1.2 Horton (HO)

The Horton approach to infiltration uses the following equation:

\[ {f = f}_{c} + \left( f_{0} - f_{c} \right)e^{- kt} \tag{7.15}\]

Where:

• \(f_0\) is the initial infiltration rate in mm/hr or inches/hr (if using Units == US Customary)
• \(f_c\) is the final (indefinite) infiltration rate
• \(t\) is time in hours (period of time that the cell is wet)
• \(k\) is the Horton decay rate.

If an initial loss (IL) is specified, the initial loss is applied first, followed by the Horton infiltration. Figure 7.27 shows an example of how the infiltration rate varies over time for \(f_0\) equal to 3, \(f_c\) equal to 1 and \(k\) equal to 0.1.

Models using the Horton infiltration loss approach are provided in the Loss Options Example Model Dataset on the TUFLOW Wiki.

Figure 7.27: Example of Horton Infiltration Rate over Time

7.2.8.1.3 Initial Loss/Continuing Loss (ILCL)

The Initial Loss/Continuing Loss (IL/CL) method is a more simplistic approach compared to the Green-Ampt and Horton infiltration methods. The IL/CL method infiltrates water based on an initial volume (at any rate) then transitions to a constant rate after the initial loss volume is exceeded.

Note

The IL/CL infiltration is separate to the IL/CL materials values used to generate excess rainfall for direct rainfall simulations (refer to Section 7.2.7.4 for further information).

Models using the IL/CL infiltration loss approach are provided in the Loss Options Example Model Dataset on the TUFLOW Wiki.

7.2.8.2 Soils File (.tsoilf)

Soils infiltration is applied to the model by defining a soils (.tsoilf) file, which is read into the .tcf using the Read Soils File command. Table 7.17 presents the parameters of the .tsoilf file and Figure 7.28 shows an example of a completed file. A number of example models demonstrating the various soil options are available in the Loss Options Example Model Dataset on the TUFLOW Wiki. The soils (.tsoilf) file is similar to the materials file where you assign a positive integer ID to each soil, define the infiltration method (options are “NONE”, “GA”, “HO” and “ILCL”) followed by the soil parameters as the remaining values. The Porosity (saturated moisture content), Initial Moisture (fraction of the soil that is initially wet), Max Ponding Depth and Horizontal Hydraulic Conductivity are all optional with default values of 1.0, 0.0, 0.0 and 0.0 respectively.

Note

Table 7.17 does not include the additional attributes in regards to horizontal advection or the two addiitonal soil types available if using TUFLOW HPC (SCS and convection). For these attributes, refer to Table 7.30 instead.

Table 7.17: Soil File (.tsoilf) Parameters

Column No.	No Infiltration	Green-Ampt	Green-Ampt.1	Horton	Initial/Continuing Loss
1	Soil ID	Soil ID	Soil ID	Soil ID	Soil ID
2	NONE	GA	GA	HO	ILCL
3		USDA Soil Type (see Table 7.16) (mandatory)	Suction (mm or inches) (mandatory)	Initial Loss (mm or inches) (mandatory)	Initial Loss (mm or inches) (mandatory)
4		Initial Moisture (Fraction) (optional)	Hydraulic Conductivity (mm/h or in/h) (mandatory)	Initial Loss Rate (\(f_0\)) (mm/h or in/h) (mandatory)	Continuing Loss (mm/h or in/h) (mandatory)
5		Max Ponding Depth (m or ft) (optional)	Porosity (Fraction) (mandatory)	Final Loss Rate (\(f_c\)) (mm/h or in/h) (mandatory)	Porosity (Fraction) (optional)
6			Initial Moisture (Fraction) (optional)	Exponential Decay Rate (k) (\(h^{-1}\)) (mandatory)	Initial Moisture (Fraction) (optional)
7			Max Ponding Depth (m or ft) (optional)	Porosity (Fraction) (optional)
8				Initial Moisture (Fraction) (optional)

Note for the Green-Ampt method:

• the initial moisture and porosity values in Table 7.17 above are fractions. The soil capacity is defined as the difference between the saturated moisture content (porosity) and the initial moisture content, hence the initial moisture should not exceed the porosity otherwise the soil capacity is set to zero and no infiltration will occur for that soil type. A WARNING 2508 is issued if this occurs.
• the max ponding depth value is an optional value that can be used, if desired, to set a limit for the depth of ponded water (\(h_0\)) value used in the Green-Ampt Equation 7.12. The minimum of the water depth and the max ponding depth value is used as the \(h_0\) value. The default max ponding depth value is 0, to be consistent with the basic form of the Green-Ampt equation, as hydrology models do not necessarily have a depth calculated at cells.

Figure 7.28: Example Soils .tsoilf File Format

One or more soils need to be specified globally and/or via GIS layers/raster grids to activate the infiltration feature. The first attribute of the GIS layer/s must be the Soil ID referenced within the .tsoilf file, in the same way that a 2d_mat layer references a Material ID stored within the materials .tmf or .csv file.

Table 7.18: 2D Soil (2d_soil) GIS Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Soil	The soil ID value referenced within a Soils File (see Table 7.17).	Integer

Each soil type can have a different infiltration method and infiltration parameters assigned to it, including a no infiltration option, as shown in Table 7.17. Other parameters that can optionally be set are:

• The imperviousness of the surface (see Fraction Impervious parameter within the materials definition in Table 7.15 or Table 7.14); and
• Groundwater or impervious level beneath the ground surface (refer to Section 8.8.2).

The current limit to the number of soils types is 1,000 for TUFLOW Classic simulations and 255 for TUFLOW HPC simulations.

7.2.9 Hydraulic Structures

7.2.9.1 Introduction

Bridges, box culverts and other structures that constrict flow can be modelled in 2D directly, rather than using 1D elements, provided the flow width of the structure is of a similar or larger size than the 2D cell size. Cells are modified in their height (invert and obvert) and width. For bridges, additional losses associated with piers and flow reaching the underside of the deck are specified. For box culverts, the additional resistance for vertical walls is specified. Additional form losses (energy head losses) can be specified for all flow constrictions.

Weir flow (across levees and other embankments) is modelled in 2D domains by default, this is handled differently between TUFLOW Classic and TUFLOW HPC. Refer to Section 7.4.2 and Section 7.3.4 respectively.

Modelling hydraulic structures in 2D domains must be carried out with a good understanding of the limitations of different approaches and the different flow regimes possible. The modeller must understand why and where the energy losses occur when assigning form losses to a 2D cell or contraction and expansion losses to a 1D element (Syme, 2001).

It is important to note that contraction and expansion losses associated with structures are modelled differently in 1D and 2D schemes. 1D schemes rely on applying form loss coefficients, as they cannot simulate the horizontal or vertical changes in velocity direction and speed. 2D schemes model these horizontal changes and do not require the introduction of form losses to the same extent as that required for 1D schemes. 2D schemes still however require the introduction of additional form losses since they do not model losses in the vertical or fine-scale horizontal effects (such as around a bridge pier). See Syme (2001) for further details.

The following webinars by Bill Syme and Greg Collecutt (two TUFLOW Developers) discuss the theory behind the energy losses and affluxes modelling associated with hydraulic structures.

• Modelling Energy Losses at Structures
• 1D, 2D & 3D Hydraulic Modelling of Bridges

Additionally, the 2D Hydralic Structures Wiki Page contains information regarding 2D hydraulic structures, including theory, setup and Frequently Asked Questions (FAQs).

It is recommended that the losses through a structure be validated through:

• Calibration to recorded information (if available).
• Cross-checked using desktop calculations based on theory and/or standard publications (e.g. “Hydraulics of Bridge Waterways” (Bradley, 1978) or “Guide to Bridge Technology Part 8, Hydraulic Design of Waterway Structures” (Austroads, 2018).
• Crosschecked with results using other hydraulic software.

To validate structure flows and energy losses:

• Specify time-series output (PO) lines of flow (Q_) and flow area (QA) across the structure (see Section 11.3.2). Upstream and downstream water levels may also be specified using PO points or extracted from the map (e.g. XMDF) output.
• Using the upstream and downstream water levels, determine whether flow is upstream or downstream controlled and estimate the flow using theoretical equations or other methods.
• Using publications such as “Hydraulics of Bridge Waterways” (Bradley, 1978) or “Guide to Bridge Technology Part 8, Hydraulic Design of Waterway Structures” (Austroads, 2018), determine the energy loss coefficient and compare this with the total energy loss calculated in the model. The total energy loss is given below. Clearly, any energy losses associated with bed resistance (e.g. Manning’s equation) need to be allowed for by subtracting this term from the calculated head difference \((h_{1} - h_{2})\).

\[ \zeta_{total} = ({h}_{1} - h_{2})\frac{2g}{V^{2}} \tag{7.16}\]

Where:

• \(\zeta_{total}\) = total energy loss
• \(h_{1}\) = upstream head
• \(h_{2}\) = downstream head
• \(V\) = depth and width averaged velocity (i.e. \(\frac{Q\_}{QA}\))

Table 7.19 lists the recommended approaches for modelling 1D and 2D structures using TUFLOW. This is discussed further in the following paragraphs.

Table 7.19: Hydraulic Structure Modelling Approaches

Structure	1D Approach	2D Approach
Box Culvert (For culverts with a steep slope, use a 1D element)	OK	OK
Circular Culvert	OK	N/A
Bridge	OK	OK
Weirs	OK	OK

1D Approach

1D structures are discussed in Section 5.7. 1D modelling is the preferred approach where the total structure width is less than one or two 2D cells. Entry and/or exit losses are defined for each structure. Testing has shown that these losses may need to be reduced where the structure width is significant compared with the cell size (Syme, 2001).

1D structures can be linked to the 2D domain using either an SX or HX connection, as outlined in (see Section 10.2). The influence of these connection types on the modelled flow behaviour is shown in Figure 7.31.

SX Link: Momentum is not transferred into or out of the 1D element to/from the 2D domain. “Suppressed” flow patterns in the 2D domain occur at the structure outlet when using 1D elements, especially if the structure width is significant compared with the cell size. The water tends to spread, rather than jet out, as there is no inertia across the link. The effect of this is illustrated in Figure 7.30, which shows the effect on flow patterns and the preservation of inertia across 1D/2D links when modelling a structure. Figure 7.30 (red velocity arrows) is that using SX links, whilst Figure 7.29 (green arrows) is that using a fully 2D solution. As can be seen, using a SX link the water tends to spread from the structure outlet, as opposed to forming a jet as in the fully 2D solution which is conserving momentum. When using the SX link, a jet like effect can be created using “wing walls” in the 2D domain at the structure outlet by assigning flood free elevations to the ZU and ZV Zpts either side of where the 1D element discharges into the 2D domain.

HX Link: Momentum is not transferred into or out of the 1D element to/from the 2D domain, however the velocity field across the HX link is assumed to be undisturbed. Provided the HX link is appropriately located (i.e. perpendicular to the flow field) this produces the effect of preserving momentum, as illustrated by the dark blue arrows (Figure 7.31). Use of HX links at a structure may require a smaller 1D timestep than that required by a SX link.

2D Approach

2D representation of structures is preferred where the total structure width is greater than one or two 2D cells. The flow area must be adequately represented by the 2D Zpts and any adjustments to cell widths (see Section 7.2.9.3 and Section 7.2.9.2). The head drop across the structure during different flow regimes should be validated against other methods and/or literature. Some additional form losses are normally required to achieve correct head drop (see Syme (2001)). Momentum is transferred through the structure as shown in the top image (green arrows) in Figure 7.29, providing more realistic flow patterns than using a 1D element with a SX link, as illustrated by the middle image (red arrows) in Figure 7.30.

Figure 7.29: Flow patterns using 2D FC cells (i.e. a fully 2D solution)

Figure 7.30: Flow patterns using a 1D element connected to 2D SX links

Figure 7.31: Flow patterns using a 1D element connected to 2D HX links

7.2.9.2 2D Bridge Structures (2d_bg)

Modelling bridge structures in the 2D domain is possible using the BG Shape input layer (2d_bg). The 2d_bg layer facilitates the modelling of bridge structures based on research on energy losses (Collecutt et al., 2022). This method is similar to the Layered Flow Constrictions (described below in Section 7.2.9.3), however, some key updates have been made. The layers in a 2d_bg are called the “pier layer” and the “super structure layer”, which consists of the “deck layer” and the “rail layer”.

Bridges can be defined as either a line or polygon GIS feature using the .tgc command Read GIS BG Shape. The layer attribute description is provided in Figure 7.32 and Table 7.20.

Points can be used to spatially vary the bridge soffit, deck depth and rail depth, with the BG Shape points input layer (2d_bg_pts). The layer attribute description is provided in Table 7.21. This can be used to model an arch bridge or sloping deck. The point objects can be placed in a separate layer to the line or polygon GIS features.

The 2d_bg layer adjusts the FLC value in the vertical as follows:

\[ \zeta_{total} = \left( \zeta_{pier} + \zeta_{SuperS}\frac{y_{SuperS}}{D_{IP}} \right)\frac{(y_{pier} + y_{SuperS})}{y_{total}} \tag{7.17}\]

Where:

• \(\zeta_{pier}\) = FLC of the pier layer.
• \(\zeta_{SuperS}\) = FLC to the superstructure layer.
• \(y_{SuperS}\) = depth of water in the superstructure layer (cannot exceed \(D_{IP}\)).
• \(D_{IP}\) = depth to the inflection point from the bridge soffit.
• \(y_{pier}\) = depth of water in the pier layer (cannot exceed \(D_{pier}\)).
• \(D_{pier}\) = depth of the pier layer, i.e. from the bed elevation to the bridge soffit.
• \(\zeta_{total}\) = overall form loss coefficient applied.

The vertical distribution of the form loss coefficient has the following characteristics:

• Water level below the deck layer: The same result as the 2d_lfcsh approaches, i.e. a constant form loss based on that specified for the pier layer (\(\zeta_{pier}\)) is applied.
• Water level between the deck soffit and the inflection point: The FLC value is linearly increased from \(\zeta_{pier}\) to \(\zeta_{pier} + \zeta_{SuperS}\). The observations from the CFD and the field measurement indicate the inflection point is located around 1.6 times the bridge deck depth above the bridge soffit. The “Inflection Depth” (\(D_{IP}\)) is assumed as:

\[ D_{IP} = {IPf}_{SuperS}D_{deck}\phi_{Deck} \tag{7.18}\]

• Where:

  • \({IPf}_{SuperS}\) = a factor to set the elevation of the inflection point. (IP). Default is 1.6.
  • \(D_{deck}\) = depth of the bridge deck layer.
  • \(\phi_{Deck}\) = blockage of the bridge deck layer.

  Note that the effect of partial blockage at the rail layer can be considered by adding the rail layer depth (\(D_{Rail}\)) to the inflection depth proportionally based on the rail layer blockage (\(\phi_{Rail}\)):

\[ D_{IP} = {IPf}_{SuperS}\left( D_{deck}\phi_{Deck} + D_{Rail}\phi_{Rail} \right) \tag{7.19}\]

• Where:
  • \(D_{Rail}\) = depth of the rail layer.
  • \(\phi_{Rail}\) = blockage of the rail layer.
• Above the inflection point the FLC gradually reduces with increasing depth (in a similar manner to the 2D Layered Flow Constrictions METHOD B and METHOD C approaches). This is to simulate the transition to drowned flow and tendency to zero energy losses with increasing depth over the bridge deck.

Figure 7.32: 2d_bg Shape Attributes and Vertical Distribution of Form Loss Coefficient

The 2023-03-AD build implemented three extra approaches to adjust the FLC value between in the deck soffit and the inflection point:

BG FLC Default Approach == {LINEAR} | LINEAR-CONSTANT | PARABOLIC | INVERTED-PARABOLIC

• LINEAR: the default option introduced above. The depth averaged FLC value is linearly increased from the bottom of the bridge deck to the inflection point. After the inflection point it gradually reduces.
• LINEAR-CONSTANT: the depth averaged FLC is linearly increased from the bottom of the bridge deck to the top of the bridge deck. Between the top of the bridge deck and the inflection point, the FLC is kept constant. After the inflection point it gradually reduces.
• PARABOLIC: the depth averaged FLC is increased using a parabolic function from the bottom of the bridge deck to the inflection point. After the inflection point it gradually reduces.
• INVERTED-PARABOLIC: the depth averaged FLC is increased using an inverted parabolic function from the bottom of the bridge deck to the inflection point. After the inflection point it gradually reduces.

This approach can be changed for individual BG structure using the Options attribute (see Table 7.20). The vertical profiles of the depth averaged FLC values applied by the four methods are illustrated in the figure below:

Figure 7.33: Depth Averaged FLC Applied Using Four Adjustment Approaches

Figure 7.34: Depth Averaged FLC Applied Using Four Adjustment Approaches

Other key features of the BG layer that are different from the Layered Flow Constriction are:

• The BG layer does not apply geometry updates. This not only offers simplicity, but also is a desirable option considering that the size of TIN polygons used to modify the bridge bed elevations are often different from the size of polygon used to select BG cell faces.
• For a line layer a “Deck_Width” value of zero sets the layer as a thin line. It selects one row of faces in the direction of flow, and the FLC value is applied to the faces unchanged. A “Deck_Width” value larger than zero sets the layer as a thick line that selects a whole cell, i.e. two faces in the direction of flow. The FLC values at the faces are divided by two. The wide line feature is not supported by the BG layer. BG polygon shapes are recommended if more than 3 rows of faces must be selected.
• For a polygon layer, the “Deck_Width” defines the bridge width in the predominant direction of flow. This value is used to distribute the total FLC value to the selected faces, i.e.:

\[ \zeta_{face} = \left( \frac{\zeta_{total}}{Deck\_Width} \right) \times Cell\ Size \tag{7.20}\]

This means users no longer need to convert the FLC value to “form loss per metre”, which is required by the Layered Flow Constrictions polygons.

The overall blockage of the face is calculated using the same method as the Layered Flow Constriction:

\[ \phi_{total} = \frac{({y}_{pier}\phi_{pier} + y_{deck}\phi_{deck} + y_{rail}\phi_{rail})}{y_{total}} \tag{7.21}\]

Where:

• \(\phi_{n}\) = layer n blockage.
• \(y_{n}\) = layer n water depth (set to zero if dry and cannot exceed depth of layer).
• \(\phi_{total}\) = overall cell face blockage applied.

Note

Where a blockage of 100 percent is specified, a minimum cell face flow width factor of 0.001 (i.e. 0.1%) is applied. This prevents numerical instability that would occur if the width were zero and allows a very small discharge to pass through the constriction. The minimum width factor is applied to all faces where blockage is adjusted by 2d_bg, 2d_cwf, 2d_fcsh and 2d_lfcsh.

The FLC value applies an energy loss across 2D cell faces as:

\[ \mathrm{\Delta}h = \zeta_{total}\frac{V^{2}}{2g} \tag{7.22}\]

where \(V\) is the 2D cell face structure velocity (the velocity with area blocked out by the piers).

When setting the blockage values, considerations need to be taken into account whether the source of the FLC value is based on the approach velocity (the velocity in the absence of piers blockage) or structure velocity (the velocity with area blocked out by the piers).

When using FLC values from published guides, the Layer 1 blockage value should generally be set to zero (0), as the FLC may already account for blockage.

For example, Hydraulics of Bridge Waterways (Bradley, 1978) or Guide to Bridge Technology (Austroads, 2018) derives FLC values based on the cross-sectional averaged velocity in the absence of piers. For more information, please see the TUFLOW Wiki 2D Hydraulic Structures page.

Models using 2D bridge structures are provided in the 2D Structures Example Model Dataset on the TUFLOW Wiki.

Table 7.20: 2D Bridge Shape (2d_bg) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	ID	Unique identifier up to 32 characters in length. It may contain any character except for quotes and commas, and cannot be blank. As a general rule, spaces and special characters (e.g. \) should be avoided, although they are accepted.	Char(32)
2	Options	This attribute can be used to overwrite the global setting specified with the .tcf command BG FLC Default Approach to change FLC calulation method individually.	Char(32)
3	Pier_pBlockage	The percentage blockage of the pier Layer. For example, enter ‘5’ for a blockage of 5%. Note: there is no need to apply a blockage value if the FLC is estimated based on the approach velocity (the velocity in the absence of piers). For more discussion, please see the TUFLOW Wiki page on the 2D Hydraulic Structures.	Float
4	Pier_FLC	Pier layer total form loss coefficient.	Float
5	Deck_Soffit	The elevation of the bridge soffit (m or ft).	Float
6	Deck_Depth	The thickness of the bridge deck (m or ft).	Float
7	Deck_Width	Line: If equal to zero, sets the layer as a thin line and selects one row of faces in the direction of flow. The FLC value is applied to the faces unchanged. If larger than zero, sets the layer as a thick line that selects a whole cell, i.e. two faces in the direction of flow. The FLC value at the faces is the total FLC divided by two. Polygon: Defines the bridge width in the predominant direction of flow (m or ft). This value is used to distribute the total FLC value to the selected faces. The FLC applied at cell faces is calculated using Equation 7.20.	Float
8	Deck_pBlockage	The percentage blockage of the deck layer. Enter ‘100’ for a solid bridge deck obstruction.	Float
9	Rail_Depth	The depth of the rail layer (m or ft).	Float
10	Rail_pBlockage	The percentage blockage of the rail layer.	Float
11	SuperS_FLC	A postive SuperS_FLC specifes the combined form loss coefficient for the deck and the rail layers. Two layers are treated as a single “superstructure” layer in this new bridge method. If the SuperS_FLC attribute is left blank, and Blank BG FLC Approach is set to Method A or B, superstructure FLC is generated automatically. See Section 7.2.9.2.1. If a negative SuperS_FLC value is specified, the absolute value of SuperS_FLC is taken as a calibration factor, and is applied as a multiplier to the automatically generated superstructure FLC. If the SuperS_FLC attribute is blank, a multiplier of one is used.	Float
12	SuperS_Ipf	A factor to set the elevation of the inflection point (IP) at which the transition from pressure flow to drowned flow commences. The default value is 1.6.	Char(64)

Table 7.21: 2D Bridge Shape Points (2d_bg_pts) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Deck_Soffit	The elevation of the bridge soffit (m or ft).	Float
2	Deck_Depth	The thickness of the bridge deck (m or ft).	Float
3	Rail_Depth	The depth of the rail layer (m or ft).	Float
4	R1	Reserved for future use.	Float
5	R2	Reserved for future use.	Float
6	R3	Reserved for future use.	Float

7.2.9.2.1 Automatic Generation of Superstructure FLC

The 2023-03-AD build implemented a feature to automatically generate the superstructure FLC value based on the CFD benchmarking study (Collecutt et al., 2022) if the SuperS_FLC attribute is left blank and Blank BG FLC Approach is set to Method A or B.

Blank BG FLC Approach == {NONE} | METHOD A | METHOD B

• NONE: The default approach, i.e. the superstructure FLC must be specified manually.
• METHOD A: The blank superstructure FLC values will be calculated based on the ratio of the depth of the pier layer (\(D_{Pier}\) in the equations above) and the thickness of the super structure layer:

\[ T_{SuperS} = D_{Deck}{Blockage}_{Deck} + D_{Rail}{Blockage}_{Rail} \tag{7.23}\]

• METHOD B: The blank superstructure FLC values will be calculated based on the ratio of the depth of the pier layer (\(D_{Pier}\)) and the thickness of the deck layer:

\[ T_{SuperS} = D_{Deck}{Blockage}_{Deck} \tag{7.24}\]

With Method A and B, the superstructure FLC value is linearly interpolated within the \(D_{Pier}\)/\(T_{SuperS}\) range of 1~10 and capped at 0.6 and 0.16 outside this range, as shown in the Figure 7.36.

If a negative SuperS_FLC value is specified, the absolute value of SuperS_FLC is taken as a calibration factor, and is applied as a multiplier to the automatically generated superstructure FLC. If the SuperS_FLC attribute is blank, a multiplier of one is used.

Figure 7.35: CFD benchmarking study of FLC vs \(D_{Pier}\)/\(T_{SuperS}\) ratio (\(h_{B}\)/\(T\) in the original paper)

Figure 7.36: CFD benchmarking study of FLC vs \(D_{Pier}\)/\(T_{SuperS}\) ratio (\(h_{B}\)/\(T\) in the original paper)

7.2.9.3 Layered Flow Constrictions (2d_lfcsh)

Layered flow constrictions allow losses and blockages to be varied with water depth. This provides the opportunity to model the flow under and over a bridge deck, or a pipeline crossing a waterway.

Important

The 2023-03-AA release introduced the 2d_bg layer, see 2D Bridge Structures. It is recommended that it be used for representation of bridges. Prior to this, bridges were commonly modelled using a 2d_lfcsh.

Four vertical layers within a layered flow constriction (not GIS layers) are represented. The lower three layers each have their own attributes. Each layer is assigned its own percentage blockage and form loss coefficient. The top (fourth) layer assumes the flow in Layer 4 is unimpeded. Within the same shape, the invert of the bed, and thickness of each layer can vary in 3D.

For example, the layers of a bridge structure could be defined as follows.

• Layer 1: Beneath the bridge deck. Might be 5% blocked due to the bridge piers and have a small form loss for the energy losses associated with the piers.
• Layer 2: The bridge deck. This would be 100% blocked and the form loss coefficient would increase due to the additional energy losses associated with flow surcharging the deck.
• Layer 3: The bridge rails. These might be anything from 100% blocked (solid concrete rails) to 10% blocked (very open rails). Some form losses would be specified depending on the type of rails.
• Layer 4: Flow over the top of the rails - flow is assumed to be unimpeded.

Layered flow constrictions function by adjusting the flow width of the 2D cell so as to represent the combination of blockages of the four layers. When the flow is only within Layer 1, only the attributes of Layer 1 are applied. As the water level rises into Layer 2, the influence of the Layer 2 attributes increase as the water continues to rise. Similarly, for Layer 3 and Layer 4.

The cell side flow width is calculated by summing the flow areas of each layer (including the effects of layer blockages) and dividing by the water depth.

Note

Where a blockage of 100 percent is specified, a minimum cell face flow width factor of 0.001 (i.e. 0.1%) is applied. This prevents numerical instability that would occur if the width were zero and allows a very small discharge to pass through the constriction. The minimum width factor is applied to all faces where blockage is adjusted by 2d_bg, 2d_cwf, 2d_fcsh and 2d_lfcsh layers.

As of the 2020-10 release three options are available to specify the method in which form losses are applied:

• METHOD A (CUMULATE in releases prior to 2020-10): Accumulates losses through each of the layers in the 2d_lfcsh as the depth of water increases.
• METHOD B (PORTION in releases prior to 2020-10) – the default: Proportions the losses through each of the layers in the 2d_lfcsh based on the depth of water.
• METHOD C: Combines the METHOD A and METHOD B approaches by utilising METHOD A through to the top of Layer 3 and METHOD B above Layer 3.
• METHOD D: Allows the modeller to control the depth at which the losses start to reduce when the flow transitions between pressure flow and drowned flow. This approach is the same used by the 2d_bg layer (introduced in the 2023-03 release), and is discussed in Section 7.2.9.2. It is recommended to use the 2d_bg layer as it has the benefit of a simplified attribute table, for easier user input.

The following .tcf command can be used to set the default method to be applied to all structures in the model:

Layered FLC Default Approach == [ METHOD A | {METHOD B} | METHOD C | METHOD D ]

To specify the method on a structure by structure basis, populate the Shape_Options attribute (refer to Table 7.22) with either “Method A”, “Method B” (the default), “Method C” or “Method D”.

Whilst the default approach of using METHOD B is unchanged for backward compatibility purposes, it is recommended the new approaches of METHOD C or the TUFLOW HPC’s 2d_bg layer approach be considered as they are shown to emulate behaviour from CFD modelling more closely for pressurised flow conditions (Collecutt et al., 2022).

Detailed explanation outlining how the losses are applied is provided in the following sections.

METHOD A

If the form loss method has been set to METHOD A (CUMULATE), the losses are accumulated as the water level rises through the layers according to the following equation. This approach was replaced as the default setting for the 2016-03 release due to it producing inconsistent results where the bridge is substantially overtopped (drowned out) with a large percentage of the flow occurring through Layer 4, and the overall energy loss reducing with increasing water depth once the structure is submerged.

\[ \zeta_{total} = \zeta_{1} + \zeta_{2}\frac{y_{2}}{D_{2}} + \zeta_{3}\frac{y_{3}}{D_{3}} \tag{7.25}\]

Where:

• \(\zeta_{n}\) = layer n FLC
• \(D_{n}\) = depth of layer n
• \(y_{n}\) = layer n water depth (set to zero if dry and cannot exceed depth of layer)
• \(\zeta_{total}\) = overall form loss coefficient applied

METHOD B

If the form loss method has been set to METHOD B (PORTION), the losses are applied pro-rata according to the depth of water in each layer using the equation below. Note that Layer 4 (e.g. above the bridge deck rails) is always assumed to contribute a zero FLC. If a layer is not flooded the depth for that layer, y\(_n\), is set to zero.

\[ \zeta_{total} = \frac{{(y}_{1}\zeta_{1} + y_{2}\zeta_{2} + y_{3}\zeta_{3})}{y_{total}} \tag{7.26}\]

\[ y_{total} = y_{1} + y_{2} + y_{3} + y_{4} \tag{7.27}\]

Where:

• \(\zeta_{n}\) = layer n FLC
• \(y_{n}\) = layer n water depth (set to zero if dry and cannot exceed depth of layer)
• \(\zeta_{total}\) = overall form loss coefficient applied

METHOD C

METHOD C effectively combines the METHOD A (CUMULATE) and METHOD B (PORTION) approaches by utilising METHOD A through to the top of Layer 3 and METHOD B above Layer 3. METHOD C adjusts the FLC value in the vertical as follows.

• Water level below Layer 2: The same result as for all other approaches (i.e. a form loss based on that specified for Layer 1 is applied). This is typically used for the energy losses associated with bridge piers.
• Water level below top of Layer 3: Same approach as for METHOD A.
• Water level above Layer 3: Gradually reduces the energy loss with increasing depth by proportioning with depth. The energy loss is calculated as that accumulated from Layers 1 to 3 pro-rated by the depth of Layers 1 to 3, plus no losses pro-rated by the depth above the Layer 3.

\[ \zeta_{total} = \left( \zeta_{1} + \zeta_{2}\frac{y_{2}}{D_{2}} + \zeta_{3}\frac{y_{3}}{D_{3}} \right)\frac{{(y}_{1} + y_{2} + y_{3})}{y_{total}} \tag{7.28}\]

\[ y_{total} = y_{1} + y_{2} + y_{3} + y_{4} \tag{7.29}\]

Where:

• \(\zeta_{n}\) = layer n FLC
• \(y_{n}\) = layer n water depth (set to zero if dry and cannot exceed depth of layer)
• \(\zeta_{total}\) = overall form loss coefficient applied
• \(D_{n}\) = layer n thickness

Regardless of the FLC approach, the FLC value applies an energy loss across 2D cell faces as:

\[ \mathrm{\Delta}h = \zeta_{total}\frac{V^{2}}{2g} \tag{7.30}\]

where \(V\) is the 2D cell face structure velocity (the velocity with area blocked out by the piers).

When setting the blockage values, considerations need to be taken into account whether the source of the FLC value is based on the approach velocity (the velocity in the absence of piers blockage) or structure velocity (the velocity with area blocked out by the piers).

When using FLC values from published guides, the Layer 1 blockage value should generally be set to zero (0), as the FLC may already account for blockage.

For example, Hydraulics of Bridge Waterways (Bradley, 1978) or Guide to Bridge Technology (Austroads, 2018) derives FLC values based on the cross-sectional averaged velocity in the absence of piers. For more information, please see the TUFLOW Wiki 2D Hydraulic Structures page.

Layered FC shapes are defined as either a line or polygon GIS objects using the .tgc command Read GIS Layered FC Shape. The file attribute description is provided in Table 7.22.

Points can be used to vary the invert of the bed, and thickness of each layer in 3D. In this case, four attributes are required as outlined in Table 7.23. This can for example be used to model an arched bridge or sloping deck. The point objects can be placed in a separate layer to the line or polygon GIS features.

Important

Overlapping flow constriction inputs is not supported.

Models using layered FC shapes are provided in the 2D Structures Example Model Dataset on the TUFLOW Wiki.

Table 7.22: Layered Flow Constriction Shape (2d_lfcsh) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Invert	Performs same function as described for the Z attribute in Table 7.5 (point), Table 7.6 (line) and Table 7.7 (polygon) and is applied to the Invert elevation values. To leave the Zpt levels unchanged (i.e. use the existing Zpt elevations), enter a value of 99999. Point: Not used. Line: Invert of constriction (m above datum). If points are snapped to the line, set to 99999 and the invert specified in the points layer is used. Polygon: Used to set the elevation of the polygon perimeter only if NO MERGE option is specified and there are no points snapped to the perimeter. If points are snapped to the perimeter, set to 99999 and the invert specified in the points layer is used. Otherwise ignored and the Zpt elevations within the polygon will be adjusted by triangulated around the perimeter at every half-cell size.	Float
2	dZ	Performs same function as described for dZ in Table 7.5 and is applied to the Invert elevation values.	Float
3	Shape_Width_or_dMax	Point: Same as for Table 7.5. Line: Same as for Table 7.6. Polygon: Same as for Table 7.7.	Float
4	Shape_Options	Point: Same as for Table 7.5. Line or Polygon: Same as for Table 7.6 (line) and Table 7.7 (polygon). Additional options for Layered FLC Default Approach, this will overwrite the global setting specified with the .tcf command. METHOD A (CUMULATE in releases prior to 2020-10): Accumulates losses through each of the layers in the 2d_lfcsh as the depth of water increases. METHOD B (PORTION in releases prior to 2020-10): Proportions the losses through each of the layers in the 2d_lfcsh based on the depth of water. METHOD C: Combines the METHOD A and METHOD B approaches by utilising METHOD A through to the top of Layer 3 and METHOD B above Layer 3. METHOD D: Allows the modeller to control the depth at which the losses start to reduce when the flow transitions between pressure flow and drowned flow. From the 2023-03-AA release and onwards it is recommended to use the 2d_bg layer, see Section 7.2.9.2.	Char(20)
5	L1_Obvert	The obvert (soffit) of Layer 1. Set as an elevation (not depth). Units are in metres, or feet if using Units == US Customary.	Float
6	L1_pBlockage	The percentage blockage of Layer 1. For example, enter 5 for a blockage of 5%.	Float
7	L1_FLC	Layer 1 form loss coefficient (FLC). Used for modelling fine-scale “micro” contraction/expansion losses not picked up by the change in the 2D domain’s velocity patterns (e.g. bridge pier losses, vena-contracta losses, 3\(^{rd}\) (vertical) dimension etc.). This parameter should be used as a calibration parameter. So that this attribute is independent of 2D cell size, it has different treatment depending on the object it is attached to as follows: Line: The FLC value is the total form loss. The FLC value is automatically adjusted and applied to cell faces based on the width of the line: For thin lines, the FLC value is applied to the cell faces unchanged. For thick (whole cell) lines, the FLC value is divided by two (two cell faces in the direction of flow). For wide lines, the FLC value applied to all cell faces is calculated by the total FLC value divided by the line’s width and multiplied by the cell size. Polygon: The FLC value is the form loss per unit length. The FLC applied at cell faces is automatically adjusted to be the input FLC (form loss per unit length) multiplied by the cell size. To calculate the form loss per unit length, divide the total form loss for the structure by the length of travel in the direction of flow (i.e. bridge width). The unit of length is metres, or feet if using Units == US Customary. If a negative FLC value is specified, the absolute value is taken and applied unadjusted to all cell faces affected by the line or polygon. Note, this is not cell size independent, and if the 2D cell size is changed, this attribute also needs to be changed. The form loss coefficient is applied as an energy loss based on the dynamic head Equation 7.30 where \(\zeta_{a}\) is the FLC value.	Float
8	L2_Depth	The depth of Layer 2. Depth units are in metres, or feet if using Units == US Customary.	Float
9	L2_pBlockage	The percentage blockage of Layer 2. For example, enter 100 for a solid obstruction such as a bridge deck or pipe.	Float
10	L2_FLC	Layer 2 form loss coefficient. See notes for L1_FLC	Float
11	L3_Depth	The depth of Layer 3. Depth units are in metres, or feet if using Units == US Customary.	Float
12	L3_pBlockage	The percentage blockage of Layer 3.	Float
13	L3_FLC	Layer 3 form loss coefficient. See notes for L1_FLC	Float
14	Notes	Optional field for entering comments. Not used.	Char(40)

Table 7.23: Layered Flow Constriction Shape Point (2d_lfcsh…_pts) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Invert	Invert of constriction (metres above datum, or feet above datum if using Units == US Customary). To leave the Zpt levels unchanged (i.e. use the existing Zpt elevations), enter a value of 99999.	Float
2	L1_Obvert	The obvert (soffit) of Layer 1. Set as an elevation (not depth). Units are in metres, or feet if using Units == US Customary.	Float
3	L2_Depth	The depth of Layer 2. Depth units are in metres, or feet if using Units == US Customary.	Float
4	L3_Depth	The depth of Layer 3. Depth units are in metres, or feet if using Units == US Customary.	Float

7.2.10 Cell Modification

7.2.10.1 Storage Reduction (2d_srf)

The storage of 2D cells may be reduced (e.g. to model hypothetical filling, or reduced storage from buildings), or increased (e.g. for stability). For example, if a cell has a Storage Reduction Factor (SRF) value of 0.1, then its storage (surface area) is reduced by 10%. If the SRF value is less than zero, the storage is increased. The default SRF value is zero (i.e. no change in storage).

SRF values are assigned to cells in one or both of the following ways:

1. Using the Set SRF and Read GIS SRF (point, line or polygon) .tgc commands. The 2d_srf layer has only one attribute, a float or real value nominally called SRF (see Table 7.24).
2. Assigned to materials .csv/.tmf files as the 5\(^{th}\) and 10\(^{th}\) column values respectively (refer to Table 7.14 and Table 7.15).
3. Using the .tgc command Read Grid SRF to assign SRF values from a raster input grid.
4. Via SX connections to 1D. The default approach increases the storage in the 2D cells by the average associated 1D node storage. Refer to Section 10.2.2 and .tcf command SX Storage Approach. The applied SRF value is reported in the SRF attribute of the _grd_check layer.

Important

A user can apply a combination of material SRF values and Read GIS SRF layers. The item that is used by the simulation depends on the order of commands in the .tgc file. The last read input takes precedence. For example, if a cell’s storage is adjusted by layers read in using both the Read GIS SRF and Read GIS Mat commands, though Read GIS Mat is lower in the .tgc file, it will prevail.

Any cells that fall within the 2d_srf polygon will have the storage reduction factor applied. This reduction can be reviewed in the _grd_check file.

If a large reduction in cell storage is applied, for stability reasons it may be necessary to also reduce the cell conveyance using a higher roughness or cell flow width (see Section 7.2.10.2). For example, if a 10m cell size model has a SRF value of 0.99, without any SRF the cell area is 100m\(^2\) (10m x 10m), with a SRF of 0.99 this would be reduced to 1m\(^2\). If no reduction in conveyance for flow entering or leaving the cell is made, a significant volume of flow may enter the cell between timesteps, this would cause the water level in the cell to jump a large amount, potentially leading to oscillations in water levels.

Models using the storage reduction functionality are provided in the Boundary Condition Options Example Model Dataset on the TUFLOW Wiki.

Table 7.24: Storage Reduction Factor (2d_srf) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	SRF	The factor used to adjust the storage of 2D cells.	Float

7.2.10.2 Cell Width Factor (2d_cwf)

The .tgc commands Set CWF, Read GIS CWF (polygon) and Read Grid CWF can be used to adjust the 2D cell flow widths, in the same manner as for 2D flow constrictions (refer to Section 7.2.9.3). The CWF is a factor, for example a value of 0.1 will limit the flow width to 10 percent. The adjusted flow width applies at all depths. When reading from a polygon object in GIS format, any cell sides that fall within the polygon will have the width factor applied. The factor is entered using a 2d_cwf GIS layer (see Table 7.25). These factors can be reviewed in the uvpt_check file.

Any cell sides that fall within the 2d_cwf polygon will have the adjusted cell width applied. This can be reviewed in the uvpt_check file.

Note

Where a blockage of 100 percent is specified, a minimum cell face flow width factor of 0.001 (i.e. 0.1%) is applied. This prevents numerical instability that would occur if the width were zero and allows a very small discharge to pass through the constriction. The minimum width factor is applied to all faces where blockage is adjusted by 2d_bg, 2d_cwf, 2d_fcsh and 2d_lfcsh.

Table 7.25: Cell Width Factor (2d_cwf) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	CWF	The factor used to adjust 2D cell flow width at all depths.	Float

7.2.10.3 Form Loss Coefficient (2d_flc)

Form loss coefficients can be applied in the .tgc to add an additional energy loss at the 2D cell side (in the same manner as for 2D flow constrictions, see Section 7.2.9.3 and Section 7.2.9.2).

Unlike flow constrictions, the FLC applied using this approach does not vary with depth. The form or energy loss can be applied as fixed values with the .tgc commands Set FLC, Read GIS FLC and Read Grid FLC or on a form loss per unit length basis with the .tgc commands Set FLC/L, Read GIS FLC/L and Read Grid FLC/L. For Read GIS FLC and Read GIS FLC/L, the form loss coefficient is entered using a 2d_flc GIS layer (see Table 7.26).

The advantage to applying the FLC on a per unit length basis is that it makes these inputs independent of the 2D cell size when using regions (polygons). As such, if the 2D cell size is changed, the same energy loss will be applied to both models over the area of the region.

Any cell sides that fall within the 2d_flc polygon will have the additional form loss applied. These losses can be reviewed in the uvpt_check file.

Section 7.2.9 provides a detailed discussion about the application of additional form losses.

Table 7.26: Form Loss Coefficient (2d_flc) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	FLC	The form loss coefficient used to apply additional energy loss at 2D cell sides.	Float

7.2.10.4 Modify Conveyance

The. tgc command Read GIS Zpts Modify Conveyance can be used to modify the elevation of a series of cells based on an increase or decrease in conveyance. Three inputs are required for this operation. These are:

1. A GIS layer containing a region / polygon object within which the modification will apply;
2. A conveyance multiplication factor; and
3. A water level grid.

The command syntax is:

Read GIS Zpts Modify Conveyance == <gis layer> | <K_factor> | <Water_Level_Grid>

The terrain levels are adjusted based on the depth of water relative to the water surface provided. The elevations are modified by \(depth\left( 1 - f^{\frac{3}{5}} \right)\), which is effectively a change in conveyance of \(K\times f\) where \(K\) is conveyance and \(f\) is the multiplication factor. For example, f = 1.2 would be a 20% increase in conveyance resulting in a deepening of the waterway. The specified base water level grid is a grid surface in a supported input raster format. The grid does not need to be in the same resolution as the 2D domains, but is typically the _h_Max grid output by TUFLOW from a previous simulation (see Section 11.2.2.2).

This feature is useful for sensitivity testing the changes in flood behaviour due to a deepening or accretion of a river’s bed based on a change in conveyance.

This command would normally be applied in the .tgc, after all elevation (Zpt) commands as the Zpt adjustment is based on the Zpt elevations processed up until the location of this command in the .tgc file. An example is provided below:

Read GIS Zpts Modify Conveyance == gis\2d_mod_river.shp | 0.9 | base_h_Max.tif

Where:

• 2d_mod_river.shp contains polygons of the areas to be modified, in this case along the river. This can be in any supported vector format (.gpkg, .shp or .mif).
• 0.9 is the factor change in conveyance compared to a base case (i.e. 10% decrease in conveyance).
• base_h_Max.tif is a raster grid of the peak water level grid from a base case simulation. It is recommended to copy this grid layer from the results folder and place it somewhere under model\grid\. If this grid layer varies with AEP, it is recommended to use the one .tcf with a Variable to define the AEP, or alternatively a Scenario variable.

7.2.11 Modelling Urban Areas

The modelling of urbanised floodplains presents many more challenges to hydraulic modellers than those of rural floodplains. It requires careful consideration of the representation of buildings, walls and fences to accurately replicate the overland flow routes. With the advancement of modern computing hardware and software, 2D solutions are increasingly being used over 1D solutions to represent urban areas.

This section of the manual discusses the various methods available in TUFLOW to represent the features within an urban area. It does not aim to identify the most appropriate method, however, discusses the pros and cons to provide the modeller with guidance on selecting the most suitable method for their study.

Syme (2008) explores modelling approaches for buildings and fences, it can be found along with other publications in the TUFLOW Library.

7.2.11.1 Buildings

A number of methods can be used for representing buildings in a TUFLOW 2D model. These are:

• Utilising a higher bed resistance (roughness);
• Raising the elevation of 2D cells where buildings exist;
• Deactivating the 2D cells (code = 0) where buildings exist;
• Reducing the cell storage;
• Restricting the flow using a cell width reduction or additional form loss; or
• Applying a breakline along the upstream faces of buildings.

Utilising a high bed resistance parameter for buildings is a commonly used method to encourage the development of preferential flow paths around buildings during urban flooding scenarios. This method is commonly used and preferred over others as it accounts for the storage that the building provides once it becomes inundated (Syme, 2008). Depth varying roughness is recommended if using this approach in combination with direct rainfall modelling. A low roughness value is recommended at shallow depths, representing the rapid run-off response associated with rainfall on building roofs. High bed roughness is recommended for deeper flows when the building structure impedes overland flow. For TUFLOW HPC, the default Viscosity Formulation is Wu, which applies a default limiting Manning’s value of 0.1 in the viscosity calculation to avoid timesteps being driven down (see Section 7.3.2.4). The application of depth varying roughness is described in Section 7.2.7.

Buildings may also be “removed” from the 2D domain by categorising the buildings as “Land” cells (see Section 7.2.2) or raising the Zpt elevations above the predicted flood level (see Section 7.2.6.2). This method may be appropriate when the building/s have been designed to not flood. In this case, care should be taken to ensure the chosen cell size is able to appropriately represent the area of the building and the storage ‘removed’ from the floodplain. Representing buildings in this fashion will not take into account any below ground storage, such as that provided by underground car parks or basement properties. If using a direct rainfall approach (rain on grid), this method is most likely not suitable. Any land (inactive) cells will have no rainfall applied, meaning the flow volume will be incorrect, artificially high cells (e.g. 10m above ground levels) may cause stability issues as flows exit the building.

Flood studies often require the simulation of extreme flood events, where a great number of buildings are expected to become inundated. It may not be appropriate in these cases to “remove” the building entirely from the 2D domain. A more suitable approach may be to raise the Zpt elevations (see Section 7.2.6.2) of the building polygons to match the finished floor level, thereby specifying the elevation at which the building is able to flood. This may allow for direct rainfall simulations using a modification of the cell elevations.

The use of a thin breakline around the edges of the building can be used to prevent the building becoming a flow path, whilst not removing the storage volume (see Section 7.2.6.2). For example, three sides of the building could be blocked, this will still allow water to enter the building but not flow through it.

Layered flow constrictions (see Section 7.2.9.3) or Cell Width Factor (see Section 7.2.10.2) may be used to impede the flow of water as it passes through the building. The sides of 2D cells may be partially blocked to represent obstructions such as internal or external walls. Additionally, a form loss can be applied using form loss coefficients (see Section 7.2.10.3).

A storage reduction factor can also be used to reduce the area in the cell available for water using the SRF options (see Section 7.2.10.1). This can be used in conjunction with the other methods, for example, a building may have the storage reduced by 20% by applying a SRF of 0.2, this could be combined with a Cell Width Factor of 0.8 (20% reduction in flow area through the building).

In the absence of guidelines in the literature it is recommend that the approaches for modelling buildings be sensitivity tested or calibrated, if possible.

7.2.11.2 Roads

Roads typically represent the main flow routes through urban areas. The chosen cell size of the 2D domain is a key factor in ensuring an appropriate number of cells have been used to represent the width of the roads, as previously discussed in Section 3.2.2. Choosing too coarse a resolution may result in the roads not being accurately represented and minor flow routes not being represented at all. An unnecessarily fine resolution may translate to excessively long simulation times and stability problems if the Courant condition is exceeded. Ideally, it is recommended that the main flow routes be represented with a width of approximately 3-4 cells across. Note, when using the TUFLOW HPC SGS functionality (see Section 7.3.3) identification of locations where a breakline should be used to represent a flow control crest may be important. The SGS Breakline Detection Delta tool (see Section 7.3.3.1.3) can assist modellers in identifying where breaklines may need to be added.

For overland flood studies where a high proportion of sheet flow is predicted, differentiation between the road and pavement elevations may influence the predicted areas of inundation. A DTM based on LiDAR or ALS data, may not necessarily pick up this detail, and requiring manual modification of the topography using breaklines.

It is possible to apply the Log Law rule for very shallow flow through automatically varying Manning’s n. This is of particularly interest for modelling flow over smooth surfaces such as roads and concrete (see Section 7.2.7.2).

7.2.11.3 Fences and Walls

Fences and walls can cause significant blockages during a flood event, influencing the direction that flood waters take. Walls act to deflect the path of flood waters and if overtopped, may act as weirs. Fences may partially impede the flow of water and can be prone to becoming blocked with debris. There is also the potential for both walls and fences to collapse during a flood event.

Fences and walls are typically included within a TUFLOW model by digitising a series of breaklines. The 2d_zsh GIS layer permits the width of the breakline to be specified (see Section 7.2.6.2). A “thin” line modifying only the ZU and ZV (cell side) elevations may be preferred to represent fences, in situations where the width of the cell size is notably greater than the width of the obstruction. A “thin” line does not modify the ZC (cell centre) elevations hence has no impact on the cell storage. A “thick” line on the other hand, modifies the ZU, ZV and ZC Zpt elevations and may be more appropriate to represent wider obstructions such as railway embankments.

The 2d_vzsh GIS layer (see Section 7.2.6.3) may be used to represent an embankment that collapses during the simulation. The breach may be triggered at a specified time, by the water level at a defined location or the water level difference between two locations exceeding a specified amount.

7.3 TUFLOW HPC Specific

This section describes 2D domain functionality that is unique to TUFLOW HPC.

7.3.1 Quadtree

Quadtree grid refinement functionality allows users to vary the 2D cell resolution within a TUFLOW HPC model. It works by recursively dividing square TUFLOW cells into four smaller square cells, enabling higher resolution only where needed. Up to 9 levels of cell size refinement are permitted. All cells, for all levels of refinement, share a common orientation. A conceptual quadtree grid demonstrating recursive cell refinement is presented in Figure 7.37.

The Quadtree solver uses a modified version of the HPC solver. It uses an explicit finite volume solution that is 2\(^{nd}\) order in space and 4\(^{th}\) order in time. When applied to a single resolution grid (i.e. no refinement), the Quadtree solver produces near identical results to the HPC solver for the same resolution.

Note that the Quadtree solver uses a modified cell selection methodology for sampling datasets from line geometries, as described in Section 7.2.5.2.

For further background and implementation guidance, refer to:

• Future of 2D Hydraulic Modelling AWS Webinar - discusses the benefits and implementation of Quadtree
• How to Implement Quadtree video
• TUFLOW Tutorial Module 7
• Multiple Domain Model Design Example Model Dataset - contains numerous example models that use Quadtree

Figure 7.37: Conceptual TUFLOW HPC Quadtree Grid

7.3.1.1 Quadtree .tcf Commands

To use Quadtree, include the following in the .tcf (see Section 4.1.5):

• Set the solution scheme to HPC:

  Solution Scheme == HPC

• Set the hardware to GPU (recommended, but not mandatory):

  Hardware == GPU

• Specify a Quadtree Control File (.qcf):

  Quadtree Control File == <path_to_qcf>

For example:

Solution Scheme == HPC
Hardware == GPU
Quadtree Control File == ..\model\example_001.qcf

Note

Quadtree can be run without refinement (i.e. with a fixed cell size):

Quadtree Control File == Single Level

A simulation run in this manner will likely run slower than HPC, however it may use less GPU memory. This is because Quadtree only stores data for active cells, whereas fixed grid HPC uses the model domain (including inactive cells) to allocate GPU memory, see Section 16.4.3 for further discussion. CPU memory usage for HPC and Quadtree is likely to be similar if running without refinement.

7.3.1.2 Quadtree Control File (.qcf)

The .qcf defines grid refinement areas and optionally the model location and extent. Appendix H lists and describes the available .qcf commands.

7.3.1.2.1 Mandatory Commands

The following commands are mandatory in a .qcf file.

---

Read GIS Nesting == <gis_file>

Defines the GIS polygons used for grid refinement (i.e. different nesting levels). The GIS file must be in the 2d_qnl format (Table 7.27). See Section 7.3.1.3 for details.

---

Model Origin and Extent == AUTO | TGC

Controls how the model origin and extents are defined:

• AUTO: Uses the extents of the Level 1 GIS polygon. For this option, a Level 1 nesting polygon must be defined in the 2d_qnl GIS layer (Table 7.27), otherwise ERROR 2827 is produced.
• TGC: Uses the model origin and extent defined in the .tgc file (see Section 4.1.8).

Note:

• Default is AUTO when Quadtree Control File == <path_to_qcf>
• Default is TGC when Quadtree Control File == Single Level

This command does not need to be specified if the default behaviour is appropriate.

---

7.3.1.2.2 Optional Commands

The following commands are optional in a .qcf.

---

Base Cell Size == <cell size in m/ft> | {TGC}

Sets the Level 1 (parent) cell size, which is the largest cell size in the Quadtree grid:

• A numerical value (<cell size in m/ft>), overrides the cell size defined in the .tgc file
• If not specified, or set to TGC, the cell size defined in the .tgc is used

Note

Metric (SI) units are TUFLOW’s default. To set Base Cell Size in ft, ensure Units == US Customary has been specified in the .tcf.

---

Orientation Angle == <angle in degrees> | OPTIMISE | {TGC}

Defines the model orientation angle:

• <angle in degrees>: Defines the grid orientation angle and overrides any angle or location set in the .tgc.
• OPTIMISE: Automatically determines the orientation angle using the parent Level 1 polygon to minimise the area of the bounding rectangle. This can potentially reduce temporary memory usage during startup (no effect on runtime memory). Requires a Level 1 nesting polygon in the 2d_qnl layer (Table 7.27).
• If not specified, or if set to “TGC”, the orientation angle in the .tgc file is used.

---

Quadtree Mesh Processing Method == {FAST} | Memory Efficient

Controls how the Quadtree grid geometry is processed during pre-processing.

During this stage, a temporary hidden 2D domain is used for areas of refinement to allow fast processing of geometry on a regular grid. By default, each nesting level is treated as a separate 2D rectangular domain. For example, three nesting levels result in the .tgc being processed three times.

This command allows switching to a more memory efficient approach, where each separate polygon in each nesting level of the 2d_qnl layer (see Section 7.3.1.3) is processed as a separate domain.

An example of when this command could significantly improve memory efficiency is if the Quadtree model’s finest resolution nesting level includes several polygons areas that are very far apart. The default FAST method would require a temporary hidden 2D domain that bounds the entire nesting level and is processed at that finest resolution. Whereas the Memory Efficient method will split the nesting level to process each fine resolution polygon separately, requiring smaller temporary hidden 2D domains, but several more iterations of processing the .tgc.

• FAST (default): Processes geometry by nesting level. Faster startup, but higher temporary memory usage.
• Memory Efficient: Processes each polygon separately. Lower memory usage, but may increase model startup time.

The selected processing method does not affect the speed of the hydraulic computations or the memory demand during the hydraulic calculations.

---

7.3.1.3 Defining Grid Refinement Polygons

A 2d_qnl (Quadtree Nesting Level) GIS layer (see Table 7.27) is used to define the location and levels of grid refinement. The 2d_qnl layers can only contain polygon/region objects. All other GIS object types (lines, points, etc.) are ignored.

The nesting level attribute must be in the range 1 to 9. Each level represents a refinement of the base cell size:

• 1: Base cell size (Level 1) - see Base Cell Size.
• 2: Half of the base cell size (Level 2)
• 3: A quarter (1/4 \(^{th}\)) of the base cell size (Level 3)
• 4: An eighth (1/8 \(^{th}\)) of the base cell size (Level 4)
• …
• 9: 1/256 \(^{th}\) of the base cell size (Level 9)

Note that there should only be one Level 1 polygon defined in the 2d_qnl layer, but for all other levels there is no limit on the number of polygons.

Table 7.27: 2D Quadtree (2d_qnl) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Nest_Level	Defines the Quadtree nesting level used for grid refinement. The value must be in the range of 1-9. See Section 7.3.1.3 above.	Integer

A level 1 polygon is not required if the model origin and extent are defined in the .tgc file (Model Origin and Extent == TGC). In this case, the Level 1 temporary hidden 2D domain matches the .tgc computational domain, with the active area defined in the .tgc acting as the Level 1 polygon.

For example, if the .qcf contains the following commands and the only 2d_qnl polygon is Level 3 (red polygon in Figure 7.38), the grid created is based on the rectangular computational domain in the .tgc file (as shown by the thick dashed black line) with inactive cells removed from the grid to reduce memory (see Section 7.2.2).

Base Cell Size == TGC
Model Origin and Extent == TGC
Orientation Angle == TGC
Read GIS Nesting == gis\2d_qnl_example_R.shp

Figure 7.38: Example Showing Removal of Inactive Cells

7.3.1.4 Grid Refinement Behaviour

Nested Polygons (no automatic refinement required):

If a nesting polygon is fully contained within a polygon of the next higher level (e.g. a Level 3 polygon is defined within a Level 2), then no automatic grid refinement is required. An example of this is shown in Figure 7.39.

Figure 7.39: Example of Fully Nested Quadtree Nesting Level Polygons

Non-Nested Polygons (automatic refinement):

If a nesting polygon is not contained within a polygon of the next higher level (e.g. a Level 4 polygon is defined within a Level 1 or Level 2 polygon), TUFLOW automatically generates intermediate refinement areas. For example, Figure 7.40 shows the grid generated when transitioning from a Level 1 to Level 3 and a Level 1 to Level 5.

Figure 7.40: Examples of Automatic Quadtree Refinement

7.3.1.5 Quadtree Topographic Updates

The following summarises Quadtree specific considerations when implementing 2d_zsh topographic modifications (see Section 7.2.6.2):

• For 2d_zsh polygons with merged elevations at undefined vertices spanning multiple Quadtree nesting levels:
  • Review the _zsh_zpt check layer to ensure appropriate representation
  • CHECK/WARNING 2934 will be written in this instance
• For 2d_zsh line geometries:
  • A modified cell cross-hair selection method, (“spider web” approach), is used (see Section 7.2.5.2)
  • Instead of cross-hairs, imaginary lines are drawn from the cell centre to each connected cell face, ensuring continuous cell selection.

7.3.1.6 Quadtree Output Options

7.3.1.6.1 Check Files

For Quadtree models, additional attributes are included in the _zpt and _grd check layers:

• “qpc_id”: the Quadtree cell ID
• “dx”: the cell size

The dx value depends on the base cell size and the Quadtree nesting level. For example, with a base cell size of 5m and nesting level 2, the dx attribute would be 2.5m.

7.3.1.6.2 Map Outputs

In Quadtree models with varying cell sizes, care is required when interpreting the grid map outputs.

By default, the smallest Quadtree cell size is used to define the Grid Output Cell Size. This can result in many output grid cells within the largest Quadtree cells. For example, with 4 levels of refinement, the largest 2D cell will contain \(8^2 = 64\) output grid cells.

TUFLOW uses a TIN interpolation based on the output values at the 2D cell centres and corners. This may produce a “saw tooth” appearance where there are large differences between cell centre and corner values.

The following command can be used to control the interpolation method used for generating grid outputs.

Grid Output Interpolation == CC EXCLUDE H | CC ALL | {TIN}

7.3.2 HPC Turbulence

This section contains information relating to sub-grid-scale turbulence options for TUFLOW HPC.

7.3.2.1 Overview

For TUFLOW HPC, three options exist (Constant, Smagorinsky and Wu) for specifying eddy viscosity for the 2D domains to approximate the effect of sub-grid-scale turbulence. Collecutt et al. (2020) discusses the turbulence models in respect to model cell size sensitivity, resulting in the Wu approach (Section 7.3.2.4) becoming the default for TUFLOW HPC from the 2020-01 release and onwards. Discussion on the change in defaults is available in the TUFLOW Wiki HPC Turbulence FAQ page.

The Viscosity Formulation and Viscosity Coefficient commands are used to set the formulation and coefficients. Note that these commands may also be used to select the non-Newtonian module, which is discussed in Section 7.3.6 since it is not a sub-grid turbulence model.

7.3.2.2 Constant Eddy Viscosity

The Constant approach (Viscosity Formulation == CONSTANT) applies a constant value throughout the model, irrespective of velocity gradients and variations. This is generally satisfactory when the bed resistance is dominant or when modelling large ocean basins (in which case the recommended coefficient for the constant formulation is 1m\(^2\)/s), but is not recommended for riverine flow paths.

7.3.2.3 Smagorinsky Approach

TUFLOW’s Smagorinsky approach (Viscosity Formulation == SMAGORINSKY) applies a hybrid method that includes elements of both the Constant and Smagorinsky approaches, as given by the equation below. This formulation is better than the CONSTANT option, as the eddy viscosity is recalculated every timestep at every cell mid-side according to the change in velocity magnitude and direction. This results in higher coefficients being applied where there is greater turbulence.

However, this approach has two known limitations:

1. It assumes that all turbulent kinetic energy with length scales greater than the cell size is resolved within the velocity field. This assumption is reasonable when the 2D cell size is larger than the physical depth of the flow, but becomes incorrect when the 2D cell size is smaller than the depth of the flow – turbulence in the vertical direction contributes significantly to viscosity but is not resolved in the 2D velocity field. As a result, once the 2D cell size becomes smaller than the flow depth this approach can underpredict the viscosity, and can underpredict momentum coupling from river to flood plain.
2. It assumes that the physical length scales of the sub-grid turbulence scale with the 2D cell size. This is reasonable when the 2D cell size is similar to the depth of the flow, but becomes incorrect when the 2D cell size becomes significantly larger than the depth of the flow – the bed resistance and the finite depth of the flow prevent large scale turbulence. As a result, once the 2D cell size becomes much larger than the depth of the flow, the Smagorinsky turbulence model can overpredict the viscosity, and can overpredict the momentum coupling from river to flood plain.

The result of these two limitations, is that model results display some sensitivity to the 2D cell size. The Wu formulation (Section 7.3.2.4) has proven to be more capable in its ability to provide “out of the box” calibration and reduced sensitivity of model results to 2D cell size.

Note

If the formulation is changed, the user must also reset the coefficient using the command Viscosity Coefficient.

The Smagorinsky option remains the default when using TUFLOW Classic. The default viscosity coefficients are a Smagorinsky value of 0.5 plus a constant value of 0.05 m\(^2\)/s. The 0.5 Smagorinsky coefficient is dimensionless. The viscosity coefficient can be output using the Map Output Data Types command.

The hybrid Smagorinsky formulation used by TUFLOW is:

\[ 𝜐 = C_{c} + C_{s}A_{c}\sqrt{\left( \frac{\partial u}{\partial x} \right)^{2} + \left( \frac{\partial v}{\partial y} \right)^{2} + {\frac{1}{2}\left( \frac{|\partial u|}{\partial y} + \frac{|\partial v|}{\partial x} \right)}^{2}} \tag{7.31}\]

Where:

• \(u\ and\ v = \ \) Depth averaged velocity components in the X and Y directions
• \(x\ and\ y = \ \) Distance in the X and Y direction
• \(\nu = \ \) Horizontal diffusion of momentum (viscosity) coefficient, m\(^2\)/s or ft\(^2\)/s (if using Units == US Customary)
• \(A_{c} = \ \) Area of Cell
• \(C_{c} = \ \) Constant Coefficient (default = 0.05 m\(^2\)/s)
• \(C_{s} = \ \) Smagorinsky Coefficient (dimensionless, default = 0.5)

7.3.2.4 Wu Approach

The Wu approach (Viscosity Formulation == WU) is the default method for TUFLOW HPC from release 2020-01 and onwards. It is not available for TUFLOW Classic. This eddy viscosity (turbulence) model combines both 2D and 3D turbulence effects, and is a slightly adapted version of that described in Wu et al. (2005).

Like the Smagorinsky eddy viscosity model, it is a zero-equation model whereby the eddy viscosity coefficient can be diagnostically computed from the mean depth and velocity fields. However, unlike the Smagorinsky model, where the turbulent length scale is related to cell size, the length scales used in the Wu model are related to water depth, and hence the computed eddy viscosity is not related to or dependent on cell size. This has been shown to significantly improve the cell size convergence of the results compared to the Smagorinsky model (i.e. the results are nearly independent of the cell size provided there are enough cells across the waterways to adequately define the flows).

The computed eddy viscosity is the Pythagorean sum of 3D and 2D contributions:

\[ \begin{matrix} \nu_{T} = \sqrt{\nu_{3D}^{2} + \nu_{2D}^{2}} \\ \upsilon_{3D} = C_{3D}U^{*}L_{m} \\ U^{*} = |U|n\frac{\sqrt{g}}{d^{\frac{1}{6}}} \\ \upsilon_{2D} = C_{2D}L_{m}^{2}|\nabla U| \\ |\nabla U| = \sqrt{\left( \frac{\partial u}{\partial x} \right)^{2} + \left( \frac{\partial v}{\partial y} \right)^{2} + \frac{1}{2}\left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)^{2}} \end{matrix} \tag{7.32}\]

Where:

• \(\nu_{T}\ \) is the final total 2D viscosity
• \(\nu_{3D}\ \)is the vertical 3D viscosity contribution
• \(\nu_{2D}\ \)is the horizontal 2D viscosity contribution
• \(C_{3D}\) coefficient for 3D viscosity contribution
• \(C_{2D}\) coefficient for 2D viscosity contribution
• \(U^{*}\) is the friction velocity (is \(\sqrt{\frac{\tau_{bed}}{\rho}}\))
• \(L_{m}\) is the turbulent length scale (is lesser of flow depth or distance to wet/dry boundary)
• \(n\) is Manning’s bed friction (see note below on treatment of high Manning’s n materials)
• \(g\) is gravity
• \(d\) is water depth
• \(|\nabla U|\) is the magnitude of the 2D velocity strain tensor

The computed viscosity \(\nu_{T}\) can be output using the Map Output Data Types command. In our testing to date, we have found \(C_{3D} = 7\) and \(C_{2D} = 0\) to yield results that agree well with benchmark tests and that are not significantly dissimilar from those of the previous Smagorinsky method with its default coefficients for situations where the depth is not significantly greater than the cell-size. The values of \(C_{3D} = 7\) and \(C_{2D} = 0\) are the default values but can be changed using the Viscosity Coefficient command. This effectively ignores the 2D component in the Wu model. Alternatively, to use only the 2D component of the model (and ignore the 3D component), we have found \(C_{3D} = 0\) and \(C_{2D} = 4\) to be a suitable starting point. As always, calibration remains an essential step. If quality calibration data is available and all Manning’s values are within industry standard ranges, then we encourage the user to adjust the coefficients as necessary to achieve calibration. We anticipate that calibration should be possible with either \(C_{2D}\) or \(C_{3D}\) in the range 1 to 7. If values outside this range are required then this may indicate other errors (eg. boundary values or schematisation, poor input data, etc). As always, sensitivity testing of changes in parameters on the model results should also be performed.

It is a common practice to represent buildings in a model as roughness (i.e. cells with high Manning’s bed friction), see discussion on modelling buildings in Section 7.2.11. As the friction velocity includes the Manning’s value, this results in a very high viscosity value at those cells. This is not unreasonable, but for TUFLOW HPC, which handles the viscosity explicitly, this can drive the timestep down to small values. We have found that limiting the Manning’s values used in the Wu 3D model to a sensible maximum (for example 0.1) makes very little difference to the model results but can vastly improve solution time. This upper limit is included as an optional third value (default value of 0.1) to the Viscosity Coefficient when using the Wu approach. Note that the original high Manning’s value, as defined in the materials database, is used in the momentum equation where its primary influence is effected.

7.3.3 Sub-Grid Sampling (SGS)

Sub-grid sampling (SGS) stores elevation points sampled at a finer resolution than the 2D cell to more accurately represent the variations in terrain inside the cell, as shown in Figure 7.42. SGS functionality is introduced in Section 3.2.3. See Section 7.1.3 for a description of 2D cell control points and Section 7.2.4 for a description of the traditional approach to sampling datasets. Benchmark testing has shown the benefits of using SGS to be substantial, as discussed in Section 3.2.3.1.

SGS translates the high-resolution sampled elevation data into cell curve functions for its hydraulic calculations.

• The cell volume is stored as a non-linear function of elevation (i.e. a curve of cell volume versus elevation based on the wetted surface area versus depth measured from the lowest sampled elevation within the 2D cell).
• Across each cell side the flow area and flow width are calculated based on sampled elevations to produce a non-linear function of conveyance versus elevation (i.e. a curve of conveyance versus elevation. For determining the depth and width averaged velocity across each cell face the flow area is also retained as a curve of area versus elevation).

The elevation dataset sample resolution can be defined by the user. For example, using a 8m cell size and a 2m SGS Sample Target Distance the DEM is sampled at 5 points across each face. Shown in Figure 7.42 this translates to 25 elevation points being used to define the volume vs elevation relationship within the 2D cell, and 5 points along each cell side for defining the conveyance-elevation relationship of each face.

Figure 7.41: Traditional (non-SGS) Approach

Figure 7.42: SGS Approach

The flow or conveyance of water through 2D cells is determined through a full 2D solution of the momentum equation (see Section 7.1). The bed resistance term (e.g. Manning’s equation) and the determination of the 2D velocities are controlled by the conveyance properties of each cell face. If there are considerable variations in bed topography across the cell faces, SGS will provide a significantly more accurate representation of these variations, compared to a traditional 2D cell with a single elevation value per cell face.

If the resolution of the DEM data layer is coarser or similar to the 2D cell size, the use of SGS will have little or no benefit. For example, if the DEM resolution is 2m and the cell size is 2m there is little benefit in using SGS, although a slight improvement would result as the slope of the terrain across the cell would be picked up if using SGS. In contrast, if the cell size is 20m there would be considerable benefit in using SGS as the variations in terrain across the 2D cell will be represented. TUFLOW requires an odd number of sample locations across a cell face to ensure the cell mid-face and cell corners are sampled. The number of sample locations across a cell face is referred to as the SGS Sample Frequency.

Note

If the elevation input data has a reasonable representation of the channels, Read GIS Z Shape or Read GIS Z Line “MIN, GULLY or LOWER” lines are not required, or recommended, when using SGS. Including gully lines can overestimate conveyance in these channels.

The following section will focus on implementing SGS and options regarding SGS outputs.

7.3.3.1 SGS Methodology and Commands

To implement SGS the only additional command required in the .tcf is “SGS == ON”.

By default, using this command will set:

• the SGS Approach to Method C. Prior to 2023 release the default was Method B. The default Method C is the recommended approach and older versions are considered legacy though are still provided for backward compatibility.
• the SGS Sample Frequency to sample at the minimum grid elevation resolution that can be found from grids read in to the .tgc.

7.3.3.1.1 SGS Sample Distance

To manually set the SGS distance, the recommended approach is to use the .tcf SGS Sample Target Distance command. TUFLOW requires an odd number of sample locations across a cell face, the number of sample points is calculated from Equation 7.33

\[ \text{Sample Frequency} = \frac{\text{TUFLOW Cell Size}}{\text{Target Distance}} + 1 \tag{7.33}\]

For example, using a cell size of 10m and a sample target distance of 1m, TUFLOW will use 11 sample points along the cell face. If the target distance specified gives an even number, TUFLOW will round up to the next odd number. For example, using a cell size of 5m and a sample target distance of 1m, TUFLOW will use 7 sample points instead of 6 points along the cell face. The adopted sampling frequency is reported in the .tlf.

It is also possible to specify the SGS Sample Frequency directly instead of the target distance. When using this command with a quadtree model (see Section 7.3.1), frequencies can be set for each nested level. For example, a model with three quadtree nested levels could use SGS Sample Frequency == 9, 5, 3 and TUFLOW will use a sample frequency of 9 for the first nested level, 5 for the second, and 3 for the third. It is important to note that changing the base cell size (with or without quadtree) of the model would require a changing to the sample frequency to maintain the same sample distance. When using a target sample distance the same value can be used across all cell sizes.

To obtain the maximum resolution, the sampling frequency can be calculated as the TUFLOW cell size divided by the DEM resolution. If this yields an even number, then it should be rounded up as per Equation 7.33. Note, there is no benefit to specify an SGS target distance smaller than the DEM resolution.

If neither the SGS Sample Target Distance or SGS Sample Frequency commands are specified, a scan of the Geometry Control File (.tgc) is performed to find the minimum raster grid resolution used in any Read Grid Zpt commands, and this is used as the sampling target distance to compute the sampling frequency. If there are no gridded elevation datasets, a default sampling frequency of 11 is used. In summary, there are four ways for the SGS Sample Frequency to be defined; in decreasing order of priority they are:

1. “SGS Sample Frequency ==” command.
2. “SGS Sample Target Distance ==” command.
3. Target distance automatically based on the minimum grid elevation resolution in the .tgc.
4. Default Sampling Frequency of 11 per face / 121 per cell.

Note

The sampling frequency per face is capped at 31 by default to avoid long pre-processing times and high memory usage. In general, a sampling frequency of between 11 and 31 is sufficient for most natural water ways or artificial structures, and users may not benefit from applying a super fine sample distance against the model cell size. However, if required, high sampling resolution is possible and the set upper limit can be increased by using the following tcf command:

SGS Max Sample Frequency == <maximum_frequency>

Important

A hard limit of 127 per face (16,129 samples per cell) still applies.

7.3.3.1.2 SGS Topographic Updates

The approach by which topographic updates modify SGS elevations depends on the specific topographic feature. There are three types of SGS topographic update approaches:

1. The topographic update layer directly modifies the SGS elevations.
2. The topographic update layer modifies the cells ZC, ZU, ZV and ZH Zpts first. Then all SGS elevations within the cell are sampled from the sloping cell faces and sloping cell areas created by the modified control Zpts.
3. The entire cell/face is updated based on one flat elevation value.

Not all topography update features are currently supported for SGS. Table 7.28 summarises the update approach adopted by available features and lists any unsupported features.

Table 7.28: SGS Topography Update Approaches

Full SGS Sampling using Sample Distance	Uses ZC/ZU/ZV/ZH values	Sets Flat Cell/Face	Unsupported Features
Create TIN Zpts	Read GIS Z Shape (Breaklines)	Set Zpt	Read GIS Z HX Line
Read Grid Zpts	Read GIS Layered FC Shape (Breaklines)	Read GIS Variable Z Shape	Read GIS FC Shape
Read TIN Zpts	Read GIS Z Line	Read GIS Z Shape Route	Read GIS Zpts Modify Conveyance
Read GIS Z Shape (Regions)	Read GIS Zpts (Regions)	1D Nodes with SXZ Flag (SGS SX Z Flag Approach == Method A)	Read RowCol Zpts
Read GIS Layered FC Shape (Regions)	1D Nodes with SXZ Flag (SGS SX Z Flag Approach == Method B)	1D Pits with SXL Flag (SGS SX Z Flag Approach == Method A)	Interpolate ZC / ZHC / ZUV / ZUVC / ZUVH
	1D Pits with SXL Flag (SGS SX Z Flag Approach == Method B)		Set Code with Clean Zpt
			ZC == MIN(ZU,ZV)

The below provides a summary of SGS specific considerations for implementing 2d_zsh modifications (see Section 7.2.6.2):

1. 2d_zsh polygons and TINs directly modify all SGS resolution sample points within the polygon based on the TIN triangulation. However, 2d_zsh lines modify elevations based on the cell’s control Zpts, as described in Section 7.3.3.1.3 – SGS Z Shape Line Approach. See Table 7.28.
2. By default, a gradient is applied along the cell face (using the side and corner ZH, ZU and ZV Zpts) for all 2d_zsh line types. For thick and wide lines, a sloping cell area for the cell volumes (using ZC, ZU, ZV and ZH Zpts) is also applied. SGS points are sampled from these. See Section 7.3.3.1.3 – SGS Z Shape Line Approach for details.
   • For partially selected cells at a 2d_zsh wide line buffer’s boundary, a sloping cell area is determined based on the cell’s modified Zpts within the buffer and the unmodified Zpts outside the buffer. All SGS sample points will be updated within the partially selected cell.
   Note: It is recommended to consider alternative TIN 2d_zsh modification methods via Read GIS Z Shape 2d_zsh polygons or Create TIN Zpts to improve representation.
3. Use of the MIN or GULLY option is not recommended, as it can overestimate conveyance – see Section 7.3.4.
4. It is strongly recommended to add MAX or RIDGE lines (see Table 7.6) for every hydraulic control (e.g. levees (artificial or natural), road/rail embankments, etc.) to reduce “leakiness” caused by SGS sampling. See Section 7.3.3.1.3 – SGS Breakline Detection Delta. This includes enforcing any crests of hydraulic controls created through a polygon or TIN 2d_zsh modification.

7.3.3.1.3 Advanced SGS Commands

As noted above, the only command required to implement SGS is “SGS == ON”. One of the above commands to set the SGS sample distance or frequency is also recommended. The following section lists advanced commands related to SGS.

---

SGS Breakline Detection Delta == <maximum_delta_value>

When the SGS elevations are being processed, every cell face will store a curve of elevation versus flow width / flow area. Likewise, each cell has a curve defining the relationship of elevation versus cell volume. To make this memory efficient TUFLOW does not store information defining the specific locations within the cell or along the cell face where the low or high points in the topography occur. As a result, without breaklines, SGS will be more likely to allow “leaks” through a levee or embankment compared with running the model without SGS (see TUFLOW Wiki Quadtree and SGS FAQ for details). For example, if a 2D cell straddles a levee, SGS will sample the low elevations on either side of the levee and will flow at a lower level than if SGS was not applied and the 2D cell centre was located on top of the levee, causing the 2D cell to block flow until the levee is overtopped. Whilst the need to have breaklines is paramount, regardless of whether SGS is used, it is even more important for SGS models to accurately represent hydraulic controls such as levees (artificial or natural), road/rail embankments, fences, etc.

SGS Breakline Detection Delta is designed to identify locations where ridge lines may be required. As noted in Section 7.3.3, Read GIS Z Shape or Read GIS Z Line “MIN, GULLY or LOWER” lines are not recommended for use. See the TUFLOW Wiki Quadtree and SGS FAQ page for more details.

The SGS breakline detection check processes each 2D cell to identify the minimum water surface elevation needed for a continuous wet connection through the cell traversing between left to right and top to bottom. The maximum invert elevation of the left and right faces is subtracted from the left-right minimum water level and reports it as “uDelta”. Similarly, the top to bottom value is reported as “vDelta”. uDelta and vDelta represent the depth of water over the cell face inverts by which a natural ridge (break) line would block any flow. The maximum of uDelta and vDelta is output to the _grd_check file in the SGS_Delta_Z attribute. If either uDelta or vDelta exceeds the a spatial (GIS) CHECK 3543 message is issued at the cell location.

This command can add significant time to the initialisation process, it is recommended to perform this check during model set up only, and to turn it off once the creation of breaklines is complete.

Figure 7.43 shows an example of where CHECK 3543 indicates where the uDelta or vDelta has exceeded the (in this example, a value of 1m).

Figure 7.43: SGS Breakline Detection Delta Output

---

SGS Z Shape Line Approach == Method A | {Method B}

If set to Method A, cell faces modified by Read GIS Z Shape line are assumed flat (i.e. SGS is not applied and a rectangular section / flat cell is used). The default, Method B, applies a gradient along the face based on the cell corners and cell side Zpt values and for thick and wide lines uses the ZC, ZU, ZV and ZH values to apply a sloping cell area for the cell volume.

Figure 7.44: Diagram of SGS Z Shape Line Approach Options

---

SGS SX Z Flag Approach == Method A | {Method B}

If set to Method A, cells that are lowered by the “Z” flag on 2d_bc SX connections are assumed flat (i.e. as per the approach for no SGS). The default Method B retains the SGS information, but shifts it all to match the lowered elevation, as per Figure 7.45.

Figure 7.45: Diagram of SGS SX Z Flag Approach Options

---

SGS Infiltration Approach == {AUTO} | TOTAL AREA | WETTED AREA

It is possible to control the cell area used when infiltration is applied with SGS, using the SGS Infiltration Approach command. The three available options are:

• AUTO (the default): uses total area if there is a rainfall boundary and wetted area if no rainfall boundaries are present.
• TOTAL AREA: uses the whole cell area, regardless of the portion of the cell that is wet.
• WETTED AREA: uses only the wetted portion of the cell, for direct rainfall boundaries this may under predict the infiltration.

---

SGS Negative Rainfall Approach == TOTAL AREA | {WETTED AREA}

It is also possible to control the cell area used when negative rainfall is applied in conjunction with SGS, using the SGS Negative Rainfall Approach command. The two options available are:

• TOTAL AREA: uses the whole cell area, regardless of the portion of the cell that is wet.
• WETTED AREA (the default): uses only the wetted portion of the cell.

With SGS enabled, for cells that are partially wet the treatment of infiltration and direct rainfall is:

• For positive rainfall (i.e. rainfall on to the 2D cell), the volume source for each cell is the total cell area times the rainfall irrespective of whether the cell is partially wet or not.
• For negative rainfall (evaporation), the volume of evaporation is factored by the wet area fraction of the cell. That is, if the cell is only 10% wet, only one tenth of the cell’s total area contributes to the negative source term.
• For models with soil infiltration the infiltration rate is proportional to the wet area fraction of the cell. However, initial infiltration losses are based on the total area of the cell (i.e. infiltration will proceed at the maximum possible rate until the cumulative infiltration – also based on total cell area - equals the initial loss value) even if the infiltration occurs with the cell partially wet. This approach is adopted to conform with that required for direct rainfall, which assumes the rainfall is applied over the entire cell irrespective of whether the cell is partially wet or not. Likewise, soil capacity is based on the total cell area (i.e. infiltration will cease once the cumulative infiltration equals soil capacity), and the cumulative wet time for the Horton model will also increment for cells that are partially wet.

---

SGS Sheet Flow Approach == {Method A} | Method B

Method A (the default) assumes water is filling up from the bottom of the cell, and uses the depth above minimum cell elevation in the Manning’s bed friction calculation (see Equation 7.2 and Equation 7.3).

However, sheet flow conditions can occur in direct rainfall models on hill slopes, where water forms a thin layer and flows downstream. Using the depth above minimum cell elevation in such locations could underestimate bed friction. Method B checks whether SGS cells are subjected to sheet flow conditions or not, and uses cell averaged depth at the sheet flow cells to adjust the depth for the bed friction calculation.

Method B also fixes overestimation of water level gradient limiter for SGS models in cells experiencing sheet flow. Where flow in the 2D domain becomes upstream controlled, TUFLOW HPC triggers a weir flow approximation by applying a limiter to the water level gradient. This limiter can be wrongly triggered for SGS models in cells when water flows down a hill slope under sheet flow conditions. Method A (the default) calculates the gradient limiter using the same method as non-SGS models, while Method B applies an SGS specific calculation to address the issue above.

Figure 7.46: Diagram of SGS Sheet Flow Approach Options

7.3.3.2 SGS Output Options

7.3.3.2.1 Overview

When using SGS, care must be taken over the interpretation of the modelling results, particularly the map output results. Unlike non-SGS models where the cell bottom is assumed to be flat and only a single depth is stored per cell, SGS cells contain multiple elevation samples distributed across the cell area. The hydraulic solver calculates the water volume for each cell; however, the interpolation of output depth is no longer based on a uniform depth value across a cell (see Figure 7.48).

By sampling multiple elevations, SGS also allows cells to become wet once the water level exceeds the lowest elevation inside the cell. As a result, SGS will produce partially wet cells along wet/dry boundaries and sometimes open new flow paths that would otherwise be blocked. Compared to non-SGS models, more cells may become wet in the model output (the first cells from the left in Figure 7.48). This effect is more noticeable when cell sizes are large relative to the resolution of the terrain data.

Figure 7.47: Illustrations of cross-sections and depth in non-SGS (left) and SGS (right) cells

Figure 7.48: Illustrations of cross-sections and depth in non-SGS (left) and SGS (right) cells

Additionally, the distribution of water within a cell may vary depending on the hydraulic conditions. For most riparian flood models, water should be interpreted as filling from the bottom of SGS cells. However, for flow overtopping floodplain or direct rainfall models, water may form a thin layer and become evenly distributed across the entire cell (i.e. sheet flow). This variation in water distribution is illustrated in Figure 7.50.

To assist the interpolation of model outputs under these different conditions, a number of output options have been implemented. They are outlined in the following sections.

Figure 7.49: Illustrations of water distribution in SGS cells for different hydraulic conditions

Figure 7.50: Illustrations of water distribution in SGS cells for different hydraulic conditions

7.3.3.2.2 SGS Depth Output and Map Output Extent

An SGS model calculates the water volume at SGS cells that contain multiple elevations sampled across the cell area. The depth used for the model output can be selected using the following command:

SGS Depth Output == EXACT | {CELL AVERAGE} | MEAN | MEDIAN | MINIMUM | PERCENTILE

• EXACT: The cell volume versus elevation curve (Section 7.3.3) is used to convert the calculated water volume to the cell water level and the output depth is the water level minus the exact cell centre (ZC) and cell corner (ZH) elevations. This option is straightforward and would use the same ZC and ZH elevations as a non-SGS model to interpolate the depth. However, it may show partially wet cells as dry where the water level is below the exact ZC elevation.

• CELL AVERAGE (default): Calculates the depth by dividing the cell volume by the cell area. This option produces results similar to the EXACT option for fully wet cells, provided the bed slope is reasonably linear. In addition, it shows all partially wet cells in the model output and is ideal for displaying all inundated cells.

• MEAN / MEDIAN / MINIMUM: The cell volume versus elevation curve is used to convert the calculated water volume to the cell water level and the output depth is the water level minus the mean, median, or minimum elevations of the cells and the cell corners from the SGS sampling, with the MINIMUM option showing the maximum depth within a cell.

• PERCENTILE: The cell volume versus elevation curve is used to convert the calculated water volume to the cell water level and the output depth is the water level minus an elevation based on the specified percentile. This option requires a second argument, separated by vertical bar |, which is the percentile to use. For example, SGS Depth Output == Percentile | 25 will use the 25\(^{th}\) percentile of the SGS elevations sampled within the cell/around the cell corner.

Notes on the above options:

• For the CELL AVERAGE and MINIMUM approaches, positive depth values are output for all partially and fully wet cells. Both partially and fully inundated cells are therefore ‘mapped’ as wet in the output.
• For the EXACT, MEAN, MEDIAN and PERCENTILE approaches, output depth may become negative depending on the modelled water level and the elevation selected to interpolate the depth. Partially wet cells with negative output depths are trimmed by default from the mapped extent. Water level and all other Map Output Data Types are trimmed to this same depth (d) map extent. This default trimming behaviour can be controlled using the SGS Map Extent command.
• Map Cutoff Depth is also applied to the depth interpolated based on SGS Depth Output command.

Tip

Given that SGS models are geared towards finding all the possible flow paths and potential wet cells based on the lowest sampled elevations, it is recommended to initially apply the default CELL AVERAGE option to analyse the full inundation extent, before adjusting the Map Cutoff Depth and SGS Depth Output options.

7.3.3.2.3 SGS Water Level Output

Beta Functionality

SGS Water Level Output == {METHOD A} | METHOD B

• METHOD A (default): The cell volume versus elevation curve is used to convert the calculated water volume to the cell water level. This option assumes that water fills from the bottom of the SGS cell regardless of the hydraulic conditions. This option offers a robust water level output for riparian flood model applications where water levels within confined channels are of interest. However, this could result in the water level output ‘appearing’ below the ground in model outputs when the depth is small. Note that SGS Sheet Flow Approach == Method A always applies this output option.

• METHOD B: When SGS Sheet Flow Approach == Method B is turned on, cell averaged depth is used at sheet flow cells for hydraulic calculation (see Section 7.3.3.1.3). The water level output can be switched to cell average depth plus the cell centre elevation at sheet flow cells to reflect this hydraulic condition. At non sheet flow cells, the output water level is the same as METHOD A. METHOD B provides a more accurate representation of water level at sheet flow cells, which is particularly useful for presenting water level at steep slopes in direct rainfall models, or models with inundation over floodplains.

Warning

METHOD B relies on the elevation value that occurs exactly at cell centre. The value at the centre could contain an outlier elevation that causes the output to appear artificially high or low. For example, when the model is using a cell size much larger than features in the terrain (such as a road embankment or creek), the elevation in the centre of the cell may happen to fall on this feature. Sheet flow could be correctly triggered for the cell, but due to the elevation at the cell centre being anomalous, the water level output could appear incorrect. As this is currently tagged as beta functionality, it’s important to check water level results for anomalies which could be due to this feature, particularly if the model is using a coarse resolution.

7.3.3.2.4 SGS Velocity Based Outputs

An SGS model utilises the sub-grid elevation information to improve the velocity calculation at cell faces. However, it does not increase the model resolution by nature and only one face averaged velocity is stored per face. Therefore, the cell average depth is used to calculate the depth/velocity based map outputs, such as unit flow (q), bed shear stress (BSS), stream power (SP) and hazards, regardless of the SGS Depth Output setting (Section 7.3.3.2.2). If desired, the depth used to calculate these map outputs can be changed using the following command:

SGS Velocity Based Outputs =={CELL AVERAGE DEPTH} | SGS DEPTH OUTPUT

If set to the SGS DEPTH OUTPUT option, the output depth specified by the SGS Depth Output command is used. Note that the velocity output (v) is not affected by these commands.

7.3.3.2.5 High-Resolution Outputs

When using the default sub-grid sampling (SGS) approach (SGS Approach == Method C), the sampled elevations are retained, including topography modifiers such as breaklines, for result post processing. At the end of the geometry processing, this allows a high-resolution terrain grid to be written and used for high-resolution depth (d) or water level (h) map outputs and also for check grids such as the _DEM_Z_HR.

High-resolution outputs and the following commands are discussed further in Section 11.2.2.3 and in Appendix A.

• HR Grid Output Cell Size: sets the map output resolution for high-resolution grid outputs.
• HR Grid Output Use Face Elevations: controls the inclusion of thin breakline cell face elevations for high-resolution outputs.
• HR Interpolation Approach: can be used to control how the cell centre / corner data is interpolated across the cell.
• HR Thin Z Line Output Adjustment: for models with thin breaklines, the presence of breaklines can be used to alter the high-resolution output interpolation.

7.3.3.2.6 Check Files

When running a model with SGS enabled, some changes are made to the following check files and an additional check file is produced. Details of the changes are provided below:

• When the _zpt check layer is output with SGS turned on, the following additional attribute information is provided:

  • The “Zmin” attribute for ZC points, previously called “Elevation”, now represents the minimum elevation within the cell. The “Zmin” for ZU/ZV points is the minimum along the cell face. Note, these points are still located at the centre of the cell or cell face, but the minimum value is not necessarily at this location.
  • “ZExact” is the elevation at the exact location of cell centre, face mid-point, or cell corner.
  • “ZAvg” is the average (mean) elevation of the sampled values.
  • “ZMed” is the median elevation of the sampled values.
  • The “ZMax” attribute is the maximum elevation (i.e. the elevation at which the cell area or cell face flow width is fully wet).
  • “ZOut” is the elevation used for the SGS Depth Output (Section 7.3.3.2.2).

• When the _grd_check layer is output with SGS turned on, the following additional attribute information is provided:

  • “ZMin” is the minimum elevation within the cell.
  • “Z_Map_Cutoff” is the map output cutoff elevation for each cell. This is equal to the cell elevation basis used for the SGS Depth Output (i.e. “Zout”, as shown in the _zpt check file) plus the Map Cutoff Depth.
  • “SGS_Delta” is only output if the SGS Breakline Detection Delta command has been specified. See the TUFLOW Wiki _grd check page for details.

• When the _DEM_Z check layer is output with SGS turned on, it contains a raster based on the SGS elevations used for the depth output interpolation (based on the SGS Depth Output command)

• When SGS is turned on, an additional raster output is produced, _DEM_Zmin. It contains a raster based on the minimum SGS elevations.

7.3.4 2D Upstream Controlled Flow (Weirs and Supercritical Flow)

The original weir flow approximation in TUFLOW HPC was to adjust the water level gradient where flow in the 2D domain becomes upstream controlled. The 2023-03 release (and onwards) defaults to applying the full weir equation in the 2D domain along thin and thick breaklines based on the advanced weir equation (Section 5.7.4.3). However, users can revert to the original approach using the command:

HPC Weir Approach == Method A | Method B | {Method B Energy}

Where:

• “Method A” is the pre 2023-03 method that applies the water level gradient limiter only.
• “Method B” applies the weir equation and uses the upstream water level above the weir invert as H\(_u\), and downstream water level above the weir invert as H\(_d\).
• “Method B Energy” is the default approach. It applies the weir equation and uses the upstream energy level above the weir invert as H\(_u\), and downstream water level above the weir invert as H\(_d\).

The default weir parameters (C\(_d\), Ex, a and b) are based on the parameters used for a broad crest weir, see Table 5.16 (i.e. 0.577, 1.5, 8.55 and 0.556, respectively). The following .tcf commands can be used to change the coefficients globally.

For a thin breakline:

HPC Thin Weir Parameters == <Cd>, <Ex>, <a>, <b>

For a thick breakline:

HPC Thick Weir Parameters == <Cd>, <Ex>, <a>, <b>

A weir calibration factor can be specified globally using the Set WrF command. It can also be adjusted locally with the Read GIS WrF command, which uses a 2d_wrf GIS layer (see Table 7.29). Note WrF is not equivalent to C\(_f\) in the advanced weir equation. To be consistent with TUFLOW Classic, the weir equation is divided by WrF. This means a WrF of 1.0 represents no adjustment, while a factor greater than one will decrease the flow efficiency. The default value of WrF is 1.0.

Models using the Set WrF and Read GIS WrF functionality are provided in the 2D Structures Example Model Dataset on the TUFLOW Wiki.

Table 7.29: Weir Adjustment Factor (2d_wrf) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	WrF	The calibration factor used to adjust the weir equation.	Float

7.3.5 Infiltration and Groundwater Flow

This section contains information relating to soil infiltration and groundwater flow options for TUFLOW HPC.

7.3.5.1 Infiltration

7.3.5.1.1 SCS Curve Number (SCS)

The SCS method, added in the 2025 release, uses a “curve number” (USDA, 1986) (ranging from 100 for impervious rock down to approximately 30 for deep porous soils) to define the maximum infiltration that can soak into the ground. An initial loss for the soil can be defined, typically as a fraction (termed the “initial abstraction ratio”) of the maximum infiltration. The initial loss is not considered an extra loss - the remaining soil capacity after the initial loss is filled is equal to total soil capacity (as defined in Equation 7.34) less the initial loss. The SCS method is only available in the HPC solver - it is not available in the Classic solver. For the soil infiltration methods supported by both Classic and HPC, see Section 7.2.8.1.

SCS infiltration differs from the alternative methods (e.g. Green-Ampt, Horton and ILCL - see Section 7.2.8.1) in the following ways:

1. The maximum infiltration and the rate of infiltration are both defined by the curve number, and the initial loss by the abstraction ratio. As a result the soil capacity cannot be defined independently of the soil type, and similarly, initial water content of the soil cannot be defined independently of the soil type. Reading soil depth, porosity, and initial moisture as spatial layers are all disallowed with this method, and SCS soil types cannot be used in conjunction with any of the other soil types within the same model.
2. The infiltration is defined as a fraction of the applied rainfall, which is dependent only on the curve number and the cumulative rainfall at the location in question. Infiltration does not occur if there is no rain, and the amount that occurs is defined as a fraction of the rainfall, regardless of whether or not the soil has surface water over it. The SCS method is really a rainfall loss method defined according to soil type, rather rather than physics-based infiltration.
3. The SCS method is a rainfall loss method, and it can only be used with single layer soil files with groundwater depth unspecified, i.e. interflow (Section 7.3.5.2.2) is not allowed. The SCS method is not recommended for continuous simulations especially when sub surface flows may be of concern.

The parameters used in the tsoilf by the SCS method are shown in Table 7.30.

The maximum infiltration is defined (in inches) as

\[ S = \frac{1000 - 10C_N}{C_N} \tag{7.34}\]

Where:

• \(S\) is the maximum infiltration in inches (which is converted for models in SI units)
• \(C_N\) is the curve number

The initial loss is defined as a fraction of the maximum infiltration:

\[ IL = I_a S \tag{7.35}\]

Where:

• \(IL\) is the initial loss
• \(I_a\) is the “initial abstraction ratio”

For the duration where the cumulative rainfall is less than the initial loss, the cumulative infiltration must match the cumulative rainfall. After this time, the cumulative infiltration is defined as a fraction of the cumulative rainfall according to:

\[ CI = RFC - \frac{\left(RFC - IL\right)^2}{RFC + S - IL} \tag{7.36}\]

Where:

• \(CI\) is the cumulative infiltration
• \(RFC\) is the cumulative rainfall

Once the cumulative infiltration matches the maximum, \(S\), no further infiltration occurs. This occurs at a cumulative rainfall of

\[ RFC_s = S\left(\frac{1 - I_a + {I_a}^2}{I_a}\right) \tag{7.37}\]

This process is shown in Figure 7.51 for SCS curve number 50 with initial abstraction ratio 0.2. In this example the maximum infiltration, \(S\), is 10 inches, and is reached at 42 inches cumulation rainfall (marked as \(RFC_S\)).

Figure 7.51: SCS Infiltration for Curve Number 50, Ia = 0.2

The instantaneous infiltration rate, during the period \(IL < RFC < RFC_s\), can then be found to be:

\[ IR = \left(\frac{S}{S-IL+RFC}\right)^2 RF \tag{7.38}\]

Where:

• \(IR\) is the current infiltration rate
• \(RF\) is the current rainfall rate

For the SCS method, the infiltration rate is based on the rainfall as applied at ground level, i.e. after the material based initial loss and continuing losses have been applied. It is worth considering that the SCS method allows for an instantaneous absorption of rain up to the initial loss (initial abstraction ratio times maximum infiltration). In this case, any materials based initial loss and continuing loss should be thought of as interception losses only. With forested areas these may be large enough that they should be included, but in open grasslands they are probably insignificant. The infiltration losses are then also reduced according to the material fraction impervious, which effectively means that:

1. The soils based initial loss as defined by the initial abstraction ratio is reduced by the material fraction impervious, and
2. The total loss (i.e. equivalent to soil capacity) is also reduced by the material fraction impervious. When using sub-grid sampling, the wet area fraction method selected is not used as the losses apply to rainfall, not ponded water.

7.3.5.2 Groundwater Flow

From the 2023-03 release and onwards, horizontal flow (advection) of water (Section 7.3.5.2.1) and multiple vertical groundwater layers (Section 7.3.5.2.2) when using TUFLOW HPC are supported. This groundwater flow feature is designed to model unconfined subsurface flow movement, including the saturated groundwater flow within unconfined aquifer and the variably saturated flow in vadose zone.

Figure 7.52: Illustration of groundwater flow movement

Up to 10 vertical layers are supported, each with horizontal advection within the layer possible, and vertical fluid movement between groundwater soil layers. Movement of water within the ground in the real world is complex. The functionality introduced in the 2023-03 release is not intended to replace detailed groundwater modelling, but provides a functional mechanism by which:

• Catchment response for 2D direct rainfall (“rain on grid”) models can be better calibrated to available river flow data.
• Long-term catchment models can be constructed that offer some prediction of stream-flows long after rainfall events.
• The accuracy of real-time flood forecasts can be improved by taking into account the drying and draining of catchments that has occurred since the last rainfall event.
• Groundwater seepage can cause inundation, such as seepage under a levee bank through a porous layer (e.g., old gravel river course), resulting in flooding on the protected side of the levee.

The basic groundwater modelling approach is conceptualised by dividing the ground into vertical layers, each with a defined thickness (alternatively a z elevation for the bottom of each layer), soil type, and cumulative infiltration of water within the layer. The cumulative infiltration, combined with the porosity defined by the soil type, allows a theoretical “groundwater elevation” (also known as a “water table”) to be computed within the layer. Figure 7.53 illustrates the soil layers in adjoining cells with convective and advective fluxes.

Figure 7.53: Illustration of soil layers and flows in adjoining cells

The horizontal movement of water, \(Q^{h}\), is governed by the horizontal hydraulic conductivity value set in the tsoilf file (see Table 7.30). All groundwater layers allow for horizontal advection of the water within the layer. Horizontal advection is disabled if all horizontal hydraulic conductivities are set to zero (which is the default if they are not specified). An example of a model using horizontal hydraulic conductivity, is provided in the TUFLOW Example Model Dataset.

The vertical flow of water, \(Q^{z}\), from the surface water layer into the first (or top) groundwater layer is governed by the soil infiltration model. The infiltration into the top soil can use any of the Green-Ampt (GA), Horton (HO) and Initial/Continuing Loss (IlCL) infiltration methods, as listed in Table 7.30. For layers below the top vertical soil layer then the Convective (CO) type must be used, see Section 7.3.5.2.2. An example of a model using two vertical soil layers (i.e. the “CO” type), is provided in the TUFLOW Example Model Dataset.

The cumulative infiltration data and the vertical movement of water are stored and computed at cell centres. As such, all soil data and groundwater layer geometry is sampled from the GIS / Grid layers at cell centres. Likewise the material “fraction impervious” data is sampled at cell centres (unlike the mannings values which are sampled at face centres). For the horizontal advection calculation, which is performed at cell faces, the horizontal hydraulic conductivity used is the average of the cell centre values for the adjoining cells.

7.3.5.2.1 Horizontal Advection

For all layers, the horizontal advection of water from one cell to the next, \(Q^{h}\), is computed using the horizontal hydraulic conductivity of the soil and the horizontal gradient of a variable called the “groundwater pressure level”:

\[ Q^{h} = - k_{h}h\frac{de}{dx}dy \tag{7.39}\]

Where:

• \(k_{h}\) is the horizontal hydraulic conductivity of the soil (average of the adjoining cells used)
• \(h\) is the groundwater elevation less the layer bottom elevation
• \(\frac{de}{dx}\) is the gradient of the “groundwater pressure level”
• \(dy\) is the width of the face

Prior to the 2026.0.0 TUFLOW release, the soil porosity was included in the flux calculation. For the original formulation, please see the 2025.2 TUFLOW Classic/HPC User Manual.

Note

If advection of groundwater is not required, then the advection calculation may be disabled by setting \(k_{h}\) to zero for all soils. This is the default behaviour if \(k_{h}\) is not specified in the soil parameters.

For cells that are “unsaturated” the “groundwater pressure level” is exactly the groundwater elevation within that layer. For cells that are fully saturated, the groundwater pressure level is that of the cell in the next layer above (or the water surface elevation when considering the topmost groundwater layer). For cells that are “nearly saturated” the groundwater pressure level is transitioned between these two options. The threshold at which the transition begins is called the “groundwater blending threshold”, \(\phi\).

\[ \begin{matrix} \xi = \frac{\frac{h_{i}}{dz} - \phi}{(1 - \phi)} \\ e_{i} = (1 - \xi)\left( z_{i} + h_{i} \right) + \xi e_{i - 1} \\ e_{0} = {WSE}_{surface} \end{matrix} \tag{7.40}\]

Where:

• \(h_{i}\) is the depth of groundwater in the cell
• \(dz\) is the vertical thickness of the layer
• \(z_{i}\) is the layer bottom elevation
• \(e_{i - 1}\) is the groundwater pressure level for the layer above (or surface layer)
• \(e_{i}\) is the resulting groundwater pressure level for this layer

The layer ordering is such that the surface water is layer “0”, the highest groundwater layer is “1”, and any subsequent groundwater layers range from 2 to N (i.e. from top to bottom).

The groundwater blending threshold, \(\phi\), defaults to a value of 0.9 but can be overridden by the user in the .tcf using the command:

Groundwater Blend Threshold == <float> | {0.9}

If the sum of flows into a layer cell causes it to exceed its soil capacity, the excess is pushed upwards to the next layer. If this happens for the top-most groundwater layer, the excess is pushed into the surface water as “return flow”. This is to be expected in the creek beds but may also happen at the bottom of steep hills where the slope transitions from steep to shallow.

7.3.5.2.2 Vertical Flux For Convective Flow Layers (CO)

For layers below the top groundwater layer, the vertical flow of water, \(Q^{z}\), is computed using a “convective” hydraulic conductivity (effectively the same concept as continuing loss):

\[ Q^{z} = \ - k_{v}dxdyS_{e}^{Ex_{h}} \tag{7.41}\]

Where:

• \(Q^{z}\) is the convective vertical flux
• \(k_{v}\) is the vertical hydraulic conductivity
• \(dxdy\) is the cell area
• \(S_{e}\) is the relative saturation
• \(Ex_{v}\) is the exponents of the power functions applied to the vertical flux

Since the 2026.0.0 TUFLOW release, the adjustment based on the relative saturation \(S_{e}\) is applied by default to improve the ‘timing’ of vertical infiltration with multiple soil layers. \(S_{e}\) is calculated as:

\[ S_{e} = \frac{\theta(t)-\theta_{r}}{\theta_{s}-\theta_{r}} \tag{7.42}\]

Where:

• \(\theta(t)\) is the simulated water content (fraction)
• \(\theta_{r}\) is the residual water content (fraction)
• \(\theta_{s}\) is the saturated water content (fraction)

The residual and saturated water contents (\(\theta_{r}\) and \(\theta_{s}\)) are set to zero and the soil porosity by default. \(\theta_{r}\) and \(\theta_{s}\) can also be specified for each soil type in .tsoilf file (Table 7.30). The default exponents \(Ex_{v}\) is 2 and can be adjusted globally by including a second argument in the following command:

Groundwater Vertical Flux Model == Darcy’s Law | <Exv>

Note that \(Q^{z}\) is limited to the available cumulative infiltration stored in the groundwater layer above, and also limited to the remaining capacity of the groundwater layer below.

7.3.5.2.3 Unsaturated Groundwater Flow

The 2025.1.0 release implemented new features to adjust groundwater flux based on the saturation of soil layers. This enhancement is designed to simulate variably saturated flow (VSF) through vadose zones.

In addition to the adjustment to the vertical flux above, the same adjustment based on the power function of the relative saturation can be applied to the horizontal flux:

\[ Q^{h} = - k_{h}h\frac{de}{dx}dyS_{e}^{Ex_{h}} \tag{7.43}\]

Where:

• \(S_{e}\) is the relative saturation calculated from Equation 7.42
• \(Ex_{h}\) is the exponents of the power functions applied to the horizontal flux

This adjustment can be activated by setting \(Ex_{h}\) by including a second argument in the following command:

Groundwater Horizontal Flux Model == Darcy’s Law | <Exh>

The default \(Ex_{h}\) is zero (0), i.e. the horizontal groundwater flux is not adjusted based on \(S_{e}\).

The 2025.1 release also implemented a soil moisture dependent function for the horizontal flux based on (Van Genuchten, 1980).

Groundwater Horizontal Flux Model == Soil Moisture Dependent

With this option, the horizontal hydraulic conductivity is adjusted based on the following equation:

\[ k_{h}(S_{e}) = k_{h}S_{e}^{L}\{1-[1-S_{e}^{n/(n-1)}]^{1-1/n}\}^{2} \tag{7.44}\]

Where:

• \(S_{e}\) is the relative saturation (see Equation 7.42)
• \(n\) and \(L\) are the model parameters

Note, strictly speaking, groundwater pressure level \(e\) used to calculate the pressure gradient (\(\frac{de}{dx}\)) needs to be adjusted based on the water retention curve:

\[ S_{e} = [1+(\alpha e)^{n}]^{1/n-1} \tag{7.45}\]

Where:

• \(\alpha\) is the model parameter (\(m^{-1}\))

The model parameters \(\alpha\), \(n\) and \(L\) can be specified for each soil type in .tsoilf file (Table 7.30). The recommended model parameters for different soil classes can be found from Rosetta’s technical documentation (Schaap, 2002). However, the units used in this documentation differ from those required by TUFLOW, and some parameters require logarithmic transformation.

Note, when Groundwater Horizontal Flux Model is set to “Soil Moisture Dependent”, the adjustment based on Equation 7.44 is applied to all soil types. If users want to disable this feature for selected soil types, i.e. use standard Darcy’s Law equation, set all three parameters (\(\alpha\), \(n\) and \(L\)) to zero (0).

7.3.5.2.4 Implementation

Although the theory appears complex, the implementation of groundwater within the TUFLOW model files is straight forward. The essential steps are:

• Read a soils (.tsoilf) data file from within the .tcf (see Section 7.3.5.3).
• Define the soil types in the .tgc (see Section 7.3.5.2.5).
• Define the soil layer thickness in the .tgc (see Section 7.3.5.2.5).
• Specify locations of groundwater boundaries in the .tbc (see Section 8.4.1.2).
• Specify initial water levels within the soil layers (optional) (see Section 8.8.2).
• Add any additional output data types in the .tcf. For example:
  • Map Output Data Types (see Section 11.2.3): dGW, GWd, GWh, GWm, GWq, GWv in Table 11.4, or
  • Time series plot output data types (see Section 11.3): GWd, GWh, GWm, GWq, GWqa, GWqu, GWqv, GWVol in Table 11.9.

7.3.5.2.5 Soil Type

Setting the soil type for the groundwater layers activates the given layer. The Green-Ampt, Horton and Initial/Continuing Loss infiltration methods can be used for Layer 1. Any layers beneath Layer 1 (e.g. Layer 2, Layer 3) must be the Convective Type. For example, setting the soil type for layers 1 and 2 will activate 2 vertical groundwater layers:

Set Soil Layer 1 == 1 ! soil ID of infiltration layer e.g. GA type
Set Soil Layer 2 == 5 ! soil ID of a convective layer e.g. CO type

Multiple layers can be set simultaneously (the below will activate 3 interflow layers):

Set Soil Layer 1 == 1
Set Soil Layer 2,3 == 5 ! sets both layers 2 and 3 to Soil ID 5

This method also works for read GIS and Grid commands:

Read GIS Soil Layer 1 == <gis_layer>
Read Grid Soil Layer 2 == <gis_layer>

A maximum of 10 vertical layers is permitted and a soil ID should be assigned for each active layer. For example, if groundwater layer 2 is active, then soil IDs for layer 1 should be assigned. Note, if the layer number is not specified in a .tgc command (e.g. Set Soil), it is assumed to be layer 1.

7.3.5.2.6 Soil Thickness

The depth (thickness) of each vertical interflow layer can be set using the following commands:

[Set | Read GIS | Read Grid ] Soil Thickness Layer N == value | <gis_layer> | <grid_layer>
[Set | Read GIS | Read Grid ] Soil Base Elevation Layer N == value | <gis_layer> | <grid_layer>

If the soil layer thickness or base elevation is not set for a given layer, it is assumed to be infinite. The soil thickness sets the layer depth from the layer above and soil base elevation sets the absolute elevation of the bottom of the layer. If both methods are specified for a given grid cell, the highest of the two will be adopted. The input units should be in metres or feet (if using Units == US Customary).

7.3.5.2.7 Drying of Top Groundwater Layer

Evaporation is typically implemented with a global negative rainfall. Prior to the 2023-03 release this would only remove water in wet cells until they are dry. It would not draw from the cumulative infiltration layer. The 2023-03 release and onwards allows for negative rainfall to draw from the topmost groundwater layer under surface layer cells that are dry, thus mimicking evapotranspiration. Currently only one approach is implemented, called “factor”, which applies the defined negative rainfall, through a (globally) factored rate. To select this approach use the .tcf command:

Soil Negative Rainfall Approach == Factor

Important

A negative rainfall boundary must be applied to use this approach. This is typically applied using Global Rainfall BC.

The default approach is to not apply any negative rainfall to the topmost groundwater layer. If the “factor” method is selected, then the factor can be set using the Soil Negative Rainfall Factor .tcf command. The default factor is 1, meaning that the full value of the negative rainfall is applied. If the topmost groundwater layer is dry, then no negative rainfall is applied, and the approach does not look any further than the topmost layer. No water is drawn from subsequent ground water layers.

7.3.5.3 HPC Soils File (.tsoilf)

The .tsoilf is the same file mentioned in Section 7.2.8.2. However, when using TUFLOW HPC there are additional parameters available which can be used to enable horizontal advection, as well as two additional soil types: convective (CO) and SCS. These additional parameters are bold in Table 7.30.

Note the 2026.0.0 TUFLOW release implemented the following updates to soil type ‘NONE’ when horizontal groundwater advection is enabled:

1. surface infiltration is disabled;
2. horizontal groundwater advection is disabled; and
3. groundwater is not allowed to discharge to the surface.

Table 7.30: HPC .tsoilf Parameters

Column No.	No Infiltration	Green-Ampt	Green-Ampt	Horton	Initial/Continuing Loss	SCS Curve Number	Convective
1	Soil ID	Soil ID	Soil ID	Soil ID	Soil ID	Soil ID	Soil ID
2	NONE	GA	GA	HO	ILCL	SCS	CO
3		USDA Soil Type (see Table 7.16) (mandatory)	Suction (mm or inches) (mandatory)	Initial Loss (mm or inches) (mandatory)	Initial Loss (mm or inches) (mandatory)	Curve Number (see USDA (1986)) (mandatory)	Vertical Hydraulic Conductivity (mm/h or in/h) (mandatory)
4		Initial Moisture (Fraction) (optional)	Hydraulic Conductivity (mm/h or in/h) (mandatory)	Initial Loss Rate (\(f_0\)) (mm/h or in/h) (mandatory)	Continuing Loss (mm/h or in/h) (mandatory)	Initial Abstraction Ratio (Fraction) (mandatory)	Porosity (Fraction) (mandatory)
5		Max Ponding Depth (m or ft) (optional)	Porosity (Fraction) (mandatory)	Final Loss Rate (\(f_c\)) (mm/h or in/h) (mandatory)	Porosity (Fraction) (optional)	Residual Water Content (-) (optional)	Initial Moisture (Fraction) (optional)
6		Horizontal Hydraulic Conductivity (mm/h or in/h) (optional)	Initial Moisture (Fraction) (optional)	Exponential Decay Rate (k) (\(h^{-1}\)) (mandatory)	Initial Moisture (Fraction) (optional)	Saturated Water Content (-) (optional)	Horizontal Hydraulic Conductivity (mm/h or in/h) (optional)
7		Residual Water Content (-) (optional)	Max Ponding Depth (m or ft) (optional)	Porosity (Fraction) (optional)	Horizontal Hydraulic Conductivity (mm/h or in/h) (optional)	Parameter \(\alpha\) (\(m^{-1}\)) (optional)	Residual Water Content (-) (optional)
8		Saturated Water Content (-) (optional)	Horizontal Hydraulic Conductivity (mm/h or in/h) (optional)	Initial Moisture (Fraction) (optional)	Residual Water Content (-) (optional)	Parameter n (-) (optional)	Saturated Water Content (-) (optional)
9		Parameter \(\alpha\) (\(m^{-1}\)) (optional)	Residual Water Content (-) (optional)	Horizontal Hydraulic Conductivity (mm/h or in/h) (optional)	Saturated Water Content (-) (optional)	Parameter L (-) (optional)	Parameter \(\alpha\) (\(m^{-1}\)) (optional)
		Parameter n (-) (optional)	Saturated Water Content (-) (optional)	Residual Water Content (-) (optional)	Parameter \(\alpha\) (\(m^{-1}\)) (optional)		Parameter n (-) (optional)
		Parameter L (-) (optional)	Parameter \(\alpha\) (\(m^{-1}\)) (optional)	Saturated Water Content (-) (optional)	Parameter n (-) (optional)		Parameter L (-) (optional)
			Parameter n (-) (optional)	Parameter \(\alpha\) (\(m^{-1}\)) (optional)	Parameter L (-) (optional)
			Parameter L (-) (optional)	Parameter n (-) (optional)
				Parameter L (-) (optional)

7.3.6 Non-Newtonian Flow

TUFLOW HPC supports modelling of non-Newtonian fluids. The non-Newtonian flow theory, and an example of its use in TUFLOW is discussed in the Tsunami and Dam Failure and non-Newtonian Modelling AWS Webinar. This feature was introduced in the 2020-01 release. It is not supported in TUFLOW Classic. High-fidelity modelling of non-Newtonian fluids is complex and 3D in nature. However, with some assumptions, it is possible to model non-Newtonian fluids reasonably well in 2D. The assumptions are:

1. The 2D (horizontal) viscosity can be represented as a summation of laminar (derivative of non-Newtonian shear stress model) and turbulent (from Wu eddy viscosity model) contributions.
2. Acceleration effects are small and the laminar fluid shear stress is linear with depth.
3. The bed shear stress can be represented as a summation of laminar (from non-Newtonian shear stress model) and turbulent (from Manning’s equation) contributions.

The laminar fluid shear stress (for flow that is shearing), is assumed to follow the Hershel-Bulkley power law model:

\[ \tau = \tau_{0} + k{\dot{\gamma}}^{n} \tag{7.46}\]

Where:

• \(\tau_{0}\) = Yield shear stress of the fluid.
• \(k\) = Viscosity parameter (units Pa.s\(^n\))
• \(\dot{\gamma}\) = Shear strain rate of the fluid
• \(n\) = Power law exponent (1 = Newtonian, <1 is shear thinning, >1 is shear thickening)

For flows where the bed shear stress exceeds the yield stress, a ‘plug flow’ velocity profile is computed, as shown in Figure 7.54. For flows where the bed shear stress does not exceed the yield stress, the fluid is considered locked to the bed and does not flow. Within TUFLOW HPC, at the start of each timestep, a bed stress is calculated for each cell such that the predicted velocity profile (for the given bed shear stress) has a depth averaged mean velocity that matches the current model velocity.

Figure 7.54: Non-Newtonian Shear Stress and Velocity Profile

Note that for fluids with a low viscosity parameter, or for deep and/or faster flows, the bed friction that would arise through the use of Manning’s equation (which is based on turbulent water flow) can be equal to or greater than that from the Hershel-Bulkley equation. In this situation the mathematics is stating that the non-Newtonian flow is likely to be turbulent and behaving more like water, and the bed shear stress as computed by the Hershel-Bulkley powerlaw model is likely to be lower than what is occurring in reality. To account for this transition in a simple manner, the actual bed friction used is the sum of both non-Newtonian and Manning’s contributions (note the n in the final term on the righthand side of the following equation is Manning’s n, not the Hershel-Bulkley power law exponent):

\[ \tau = \tau_{0} + k{\dot{\gamma}}^{n} + \frac{\rho gn^{2}|U|^{2}}{d^{\frac{1}{3}}} \tag{7.47}\]

This provides a very stable and smooth transition from non-Newtonian bed shear stress to turbulent water-like bed shear stress as and when needed. The total shear stress is then utilised in the momentum equation to evolve the flow velocity. The momentum equation also considers the possibility that the fluid yield stress is sufficient to hold the fluid against the net momentum flux, in which case the fluid velocity is zeroed.

For the computation of 2D shear stresses, the yield shear stress is not used and the fluid is allowed to continue shearing down to low shear strain rates and correspondingly low shear stresses. The effective Non-Newtonian 2D viscosity is computed using the derivative of the Hershel-Bulkley equation:

\[ \begin{matrix} \mu_{2D} = kn\left| \dot{\gamma} \right|^{n - 1}\\ \left| \dot{\gamma} \right| = \max\left( \frac{2|U|}{d},\sqrt{\left( \frac{du}{dx} \right)^{2} + \left( \frac{dv}{dy} \right)^{2} + \frac{1}{2}\left( \frac{du}{dy} + \frac{dv}{dx} \right)^{2}} \right) \end{matrix} \tag{7.48}\]

Where:

• \(\left| \dot{\gamma} \right|\) = Magnitude of shear strain rate

Since the \(n - 1\) exponent can cause the viscosity to diverge to infinity at low strain rates for shear thinning fluids, it necessary to bound the resulting viscosity to a user defined upper limit, \(\mu_{high}\). For completeness a user defined lower limit, \(\mu_{low}\), is also applied.

Finally, for deeper and high velocity flows it is possible that the effective 2D viscosity of turbulent water, as computed by the Wu turbulence model, is greater than the value computed from the Hershel-Bulkey powerlaw model. Again, in this situation the mathematics is stating the non-Newtonian flow is likely to be turbulent and behaving more like water. To account for this transition the actual 2D viscosity used is the sum of both non-Newtonian and Wu turbulence model contributions:

\[ \nu_{T} = \frac{\mu_{2D}}{\rho} + 7U^{*}d \tag{7.49}\]

Where:

• \(U^{*}\) = Friction velocity = \(\frac{\sqrt{g}\ n|U|}{d^{\frac{1}{6}}}\)

This provides a very stable and smooth transition from non-Newtonian 2D viscosity to turbulent water-2D viscosity as and when needed. The value of 7 is the default Wu viscosity parameter and is not user-definable for the non-Newtonian module.

7.3.6.1 Implementation

The non-Newtonian flow model is invoked with the Viscosity Formulation command in the .tcf. If invoked there are no default values for the required coefficients (\(k\), \(n\), \(\mu_{low}\), \(\mu_{high}\), \(\tau_{0}\)). These must be supplied wit the Viscosity Coefficient command, where:

• \(k\): viscosity coefficient (Pa.s^n)
• \(n\): shear thickening exponent, must be non zero and positive (shear thinning fluids exhibit \(n\) < 1, shear thickening \(n\) > 1, and Newtonian fluids \(n\) = 1)
• \(\mu_{low}\): lower viscosity limit (Pa.s)
• \(\mu_{high}\): higher viscosity limit (Pa.s)
• \(\tau_{0}\): shear yield stress (Pa)

A step-by-step example of how to set up a non-Newtonian model is provided in TUFLOW Tutorial Module 5 on the TUFLOW Wiki. Additionally, an example model is provided in the Non-Newtonian Example Model Dataset on the TUFLOW Wiki.

Note

In TUFLOW HPC the 2D momentum diffusion is handled explicitly and can drive the model timestep to small values when the viscosity becomes large. For shear thinning models, it is important to define an upper limit, \(\mu_{high}\), that is only as large as necessary. We suggest starting with 1,000 Pa.s and adjusting lower if the model is being strongly controlled by the diffusion control number (Nd). Sensitivity testing should be done to test whether the upper viscosity limit is influencing results in the region of interest. The lower viscosity limit, \(\mu_{low}\), may be set to zero if desired.

7.3.6.2 Non-Newtonian Mixing

It is possible to consider mixing of the non-Newtonian fluid with water, for example a tailings dam failure into a flowing river. To implement this, a passive AD tracer must be included in the model, note that this requires a licence for the Advection-Dispersion (AD) Module. The AD setup is the same as detailed in Chapter 9. More than one tracer field is permissible, but the non-Newtonian mixing calculation will assume that the first tracer represents the fraction of non-Newtonian fluid, in the range 0 to 1. Care must be taken to ensure that boundary data and initialisation commands strictly adhere to this range. Where the tracer value is 1, the properties of the fluid will be exactly the non-Newtonian properties, as specified with Viscosity Coefficient command, where the tracer value is 0 the fluid will behave as water. For tracer values between 0 and 1, the non-Newtonian properties are blended with those of water using a power-law equation.

\[ \begin{matrix} \xi = a^{m} \\ \rho' = (1 - a)\rho_{w} + a\rho_{NN} \\ K' = (1 - a^{m})\nu_{w} + a^{m} K_{NN} \\ n' = (1 - a^{o})1.0 + a^{o} n_{NN} \\ \tau_{0}' = a^{p}\tau_{0,NN} \\ \end{matrix} \tag{7.50}\]

Where:

• \(a\) = the tracer value (0~1)
• \(m\) = the non-Newtonian mixing exponent for the viscosity parameter
• \(o\) = the non-Newtonian mixing exponent for the power law parameter
• \(p\) = the non-Newtonian mixing exponent for the yield stress
• \(\rho\) = density of fluid (kg/m³)
• \(\nu\) = kinematic viscosity (m²/s)
• \(K\) = the Hershel-Bulkley viscosity parameter (Pa sⁿ)
• \(n\) = the Hershel-Bulkley power law exponent
• \(\tau_{0}\) = the material yield stress (Pa)
• subscripts \(w\) denotes water, \(NN\) denotes undiluted non-Newtonian fluid
• superscript \(‘\) denotes blended value to use in the Hershel-Bulkley powerlaw model

The bed friction and 2D viscosity are computed using the non-Newtonian method with the blended parameters. To turn on the mixing calculation, the mixing exponents, \(m\), \(o\) and \(p\) can be specified with the HPC Non-Newtonian Mixing Exponent command. The values must be within the range 1 to 10. It is the user’s responsibility to source the non-Newtonian parameters that represent the undiluted mud, and also to source an appropriate values for the non-Newtonian mixing exponents.

To utilise non-Newtonian mixing in a model, the following is required:

1. Select the non-Newtonian viscosity formulation and specify the coefficients for pure non-Newtonian fluid:

Viscosity Formulation == Non-Newtonian
Viscosity Coefficient == k, n, muLow, muHigh, tau0

1. Include a passive tracer as per standard AD setup (Chapter 9), making sure boundary data and initialisation are scaled in range 0 to 1.
2. Specify the mixing exponents, m, o, and p to be used:

HPC Non-Newtonian Mixing Exponent == <float> | {0} ! used for all three of m, o, p exponents

or

HPC Non-Newtonian Mixing Exponent == <float>, <float>, <float> | {0, 0, 0} ! m, o, p exponents

Note that if the non-Newtonian mixing exponent is not specified it defaults to zero which turns the mixing formulation off, and all fluid in the model will be treated as non-Newtonian regardless of tracer value.

An example non-Newtonian mixing model is provided in the non-Newtonian section of the Solver Options Example Model Dataset on the TUFLOW Wiki.

7.3.7 Unsupported Features in TUFLOW HPC

The following features are not supported in TUFLOW HPC:

• Multiple 2D Domain (Section 7.4.4): This feature has been superseded by the TUFLOW HPC Quadtree functionality (Section 7.3.1).
• 2D Flow Constrictions (2d_fcsh and 2d_fc) (Section 7.4.6.1): This feature has been superseded by the 2d Bridge (2d_bg) (Section 7.2.9.2) and Layered Flow Constriction (2d_lfcsh) (Section 7.2.9.3) features.
• Coriolis Force (Section 7.4.5): The inclusion of the Coriolis term in the shallow water equations is possible with TUFLOW Classic only.
• Chézy and Fric options: TUFLOW HPC only supports Manning’s n bed resistance values (Section 7.4.3).

7.4 TUFLOW Classic Specific

This section describes functionality unique to TUFLOW Classic.

7.4.1 Classic Turbulence

For TUFLOW Classic, two options exist (Constant and Smagorinsky) for specifying eddy viscosity for the 2D domains to approximate the effect of sub-grid-scale turbulence. For TUFLOW Classic, the Smagorinsky approach (Section Section 7.4.1.2) is the default approach.

7.4.1.1 Constant

The Constant approach (Viscosity Formulation == CONSTANT) applies a constant value throughout the model, irrespective of velocity gradients and variations. This is generally satisfactory when the bed resistance is dominant or when modelling large ocean basins (in which case the recommended coefficient for the constant formulation is 1m\(^2\)/s), but is not recommended for riverine flow paths.

7.4.1.2 Smagorinsky

TUFLOW’s Smagorinsky approach (Viscosity Formulation == SMAGORINSKY) applies a hybrid method that includes elements of both the Constant and Smagorinsky approaches, as given by the equation below. This formulation is better than the CONSTANT option, as the eddy viscosity is recalculated every timestep at every cell mid-side according to the change in velocity magnitude and direction. This results in higher coefficients being applied where there is greater turbulence.

However, this approach has two known limitations:

1. It assumes that all turbulent kinetic energy with length scales greater than the cell size is resolved within the velocity field. This assumption is reasonable when the 2D cell size is larger than the physical depth of the flow, but becomes incorrect when the 2D cell size is smaller than the depth of the flow – turbulence in the vertical direction contributes significantly to viscosity but is not resolved in the 2D velocity field. As a result, once the 2D cell size becomes smaller than the flow depth this approach can underpredict the viscosity, and can underpredict momentum coupling from river to flood plain.
2. It assumes that the physical length scales of the sub-grid turbulence scale with the 2D cell size. This is reasonable when the 2D cell size is similar to the depth of the flow, but becomes incorrect when the 2D cell size becomes significantly larger than the depth of the flow – the bed resistance and the finite depth of the flow prevent large scale turbulence. As a result, once the 2D cell size becomes much larger than the depth of the flow, the Smagorinsky turbulence model can overpredict the viscosity, and can overpredict the momentum coupling from river to flood plain.

The result of these two limitations, is that model results display some sensitivity to the 2D cell size. The Wu formulation (Section 7.3.2.4) has proven to be more capable in its ability to provide “out of the box” calibration and reduced sensitivity of model results to 2D cell size.

Note

If the formulation is changed, the user must also reset the coefficient using the command Viscosity Coefficient.

The Smagorinsky option remains the default when using TUFLOW Classic. The default viscosity coefficients are a Smagorinsky value of 0.5 plus a constant value of 0.05 m\(^2\)/s. The 0.5 Smagorinsky coefficient is dimensionless. The viscosity coefficient can be output using the Map Output Data Types command.

The hybrid Smagorinsky formulation used by TUFLOW is:

\[ 𝜐 = C_{c} + C_{s}A_{c}\sqrt{\left( \frac{\partial u}{\partial x} \right)^{2} + \left( \frac{\partial v}{\partial y} \right)^{2} + {\frac{1}{2}\left( \frac{|\partial u|}{\partial y} + \frac{|\partial v|}{\partial x} \right)}^{2}} \tag{7.51}\]

Where:

• \(u\ and\ v = \ \) Depth averaged velocity components in the X and Y directions
• \(x\ and\ y = \ \) Distance in the X and Y direction
• \(\nu = \ \) Horizontal diffusion of momentum (viscosity) coefficient, m\(^2\)/s or ft\(^2\)/s (if using Units == US Customary)
• \(A_{c} = \ \) Area of Cell
• \(C_{\ c} = \ \) Constant Coefficient (default = 0.05 m\(^2\)/s)
• \(C_{s} = \ \) Smagorinsky Coefficient (dimensionless, default = 0.5)

7.4.2 2D Upstream Controlled Flow (Weirs and Supercritical Flow)

Where flow in the 2D domain becomes upstream controlled, TUFLOW Classic automatically switches between either weir flow or upstream controlled friction flow.

If Supercritical is set to ON (the default), the following rules apply.

• Where the bed slope at a ZU or ZV point is in the same direction as the water surface slope, tests are carried out to determine whether the flow is upstream controlled or downstream controlled. The adopted flow regime is determined by comparing the upstream and downstream controlled regime flows (preference to the lower flow) and whether the Froude number exceeds 1 (unless changed by Froude Check). The equation used for upstream controlled flow is the Manning equation, with the water surface slope set to the bed slope. This check can be disabled for backward compatibility using Froude Depth Adjustment.
• Weir flow only occurs if the bed slope is adverse (different direction) to the water surface slope. Weir flow across 2D cell sides is modelled by first testing whether the flow is upstream or downstream controlled. If upstream controlled, the broad-crested weir flow equation is used to replace the calculations for downstream controlled (sub-critical) flow conditions. Weir flow can be switched off using the Free Overfall options.

Note

The bed slope at ZU and ZV points is determined as the slope from the upstream ZC point to the ZU or ZV point in the direction of positive flow.

TUFLOW produces an increase in water level at transitions from supercritical flow to subcritical flow, as occurs with a hydraulic jump. It does not, however, model the complex 3D flow patterns that occur at a hydraulic jump, as it uses a 2D horizontal plane solution. Results in areas of transition should be interpreted with caution. It is also important to be careful presenting results in areas of supercritical flow, as complex localised flow interactions may occur that would yield higher localised water levels (such as supercritical flow surcharging against a building) – it is good practice to also view the energy levels when providing advice on flood planning levels.

If Supercritical is set to OFF, and Free Overfall is set to ON (the default), weir flow may occur on both adverse and normal bed slopes.

The weir flow switch may be adjusted globally using the .tcf command Global Weir Factor. It can also be varied spatially over the grid by setting a weir factor of zero where there is to be no automatic weir calculations using Read GIS WrF and a 2d_wrf GIS layer (see Table 7.29). The weir factor also allows calibration or adjustment where the broad-crested weir equation is applied. The weir factor is not the broad-crested weir coefficient. The broad-crested weir equation is divided by the weir factor. A factor of 1.0 represents no adjustment, while a factor greater than one will decrease the flow efficiency. Refer to Syme (2001) for further information.

Important

The global weir factor and the spatially varying value are multiplied together (i.e. one does not replace the other).

The .tcf command HX Force Weir Equation can be used to force the weir flow equation to be applied across all active HX cell sides when the flow is upstream controlled. Note, the default approach uses either weir flow or super-critical flow when the flow is upstream controlled, depending on whether the ground surface gradient from HX cell centre to cell side is adverse (weir flow) or not adverse (super-critical flow). When the flow is downstream controlled, regardless of the ground surface slope, the full 2D equations are applied, including allowance for momentum across the HX 1D/2D link.

7.4.3 Land Use (Materials)

For TUFLOW Classic, bed resistance can be set to Manning’s n, Manning’s M (1/n) or Chezy using the Bed Resistance Values command in the .tcf file. Chezy can be specified as a direct value or by bed ripple heights (using the Set FRIC, Read GIS FRIC or Read Grid FRIC commands). The default and recommended bed resistance formulation is Manning’s n. Manning’s n values can also be varied with depth, the VxD product or varied with depth using the Log Law formula (see Section 7.2.7).

If using the Chezy formula, a number of commands have been setup to provide backward compatibility. These are Depth/Ripple Height Factor Limit and Recalculate Chezy Interval.

7.4.4 Multiple 2D Domain

When using TUFLOW Classic with the Multiple 2D Domain Module (see Section 1.1.5.4), any number of 2D domains of different cell size and orientation can be combined to form one model. This TUFLOW Classic feature is discussed in Section 10.7.2.

7.4.5 Coriolis Term

When using TUFLOW Classic, it is possible to include the Coriolis term in the shallow water equations (see Section 6.3 of 2018 TUFLOW Manual). To enable this functionality, use the Latitude command.

7.4.6 Legacy Structures

The following features are considered legacy structures and are supported by TUFLOW Classic only. It is recommended, and preferred, to use the layered flow constriction shape (Section 7.2.9.3) instead, which is supported by both TUFLOW Classic and TUFLOW HPC. To convert from the flow constriction (2d_fc or 2d_fcsh) to the layered flow constriction (2d_lfcsh), refer to the discussion published in the TUFLOW Wiki 2D Hydraulic Structures.

7.4.6.1 2D Flow Constrictions (2d_fcsh and 2d_fc)

Important

These features are considered legacy structures and are supported by TUFLOW Classic only. It is recommended to use the layered flow constriction shape (Section 7.2.9.3) instead.

Flow Constriction (FC) attributes allow the user to constrict the flow across a 2D cell side as a way to define hydraulic structures, such as bridges and banks of box culverts, within a model. 2D cell sides can be modified in the following ways:

• Placing a lid (obvert or soffit) on the cell side.
• Changing the flow width of the cell side.
• Adding additional form (energy) losses.
• Including side wall friction (“BC” FC_Type only).

Note

Where a blockage of 100 percent is specified, a minimum cell face flow width factor of 0.001 (i.e. 0.1%) is applied. This prevents numerical instability that would occur if the width were zero and allows a very small discharge to pass through the constriction. The minimum width factor is applied to all faces where blockage is adjusted by 2d_bg, 2d_cwf, 2d_fcsh and 2d_lfcsh.

Flow constrictions can be digitised within GIS layers and read into the .tgc file using the command Read GIS FC Shape. This command is recommended over the previous (and now outdated) Read GIS FC.

Table 7.31 describes the different 2d_fcsh layer attributes as used by Read GIS FC Shape. If Write Check Files is set in the .tcf then the _fcsh_uvpt_check file will be created to allow you to cross-check that the changes to the cell-sides are as intended. The _fcsh_uvpt_check file contains information on the cell sides modified by the Read GIS FC Shape command. Details on the properties at ALL cell sides can be found in the _uvpt_check file.

Point objects associated with the 2d_fcsh layer can be placed in a separate layer. In this case, only the first two attributes are required. This is discussed in Section 7.2.6.7.1.

The original Read GIS FC feature continues to be supported and its attributes are documented in Table 7.32. If any 2d_fc inputs are included in the model and Write Check Files is set, then the _fc_check file will be created.

Important

Overlapping flow constriction inputs are not supported.

Table 7.31: Flow Constriction Shape (2d_fcsh) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Invert	Invert of constriction. Using metric units (default) the value is metres above the elevation datum . Using Units == US Customary, the value is feet above the elevation datum. To leave the Zpt levels unchanged (i.e. use the existing Zpt elevations), enter a value of 99999.	Float
2	Obvert_or_BC_Height	FC_Type: Blank or “BD”: Obvert (soffit) of constriction in m above datum. “BC”: Height of box culvert (m). “FD”: Floating depth (m) of the deck (i.e. depth below the water line). Enter a sufficiently high value (e.g. 99999) if there is no obvert constriction.	Float
3	Shape_Width_or_dMax	Point: Not used. Line: See same attribute in Table 7.6 for 2d_zsh layers. Polygon: Not used.	Float
4	Shape_Options	Point: Not used. Line or Polygon: MAX, RIDGE or RAISE: Only changes a Zpt (Invert) elevation if the Z Shape elevation at the Zpt is higher. MIN, GULLY or LOWER: Only changes a Zpt (Invert) elevation if the Z Shape elevation at the Zpt is lower. Line: TIN: Indicates the line is to only be used for generation of TINs within polygons (only sections of TIN lines that fall within polygons are used).	Char(20)
5	FC_Type	The FC_Type. Where: Blank for general (does not include allowances for any vertical walls or friction from underside of deck). “BC” for Box Culverts (BC is only available if Manning’s n option is specified in Bed Resistance Values). “BD” for Bridge Deck. “FD” for Floating Bridge Deck.	Char(2)
6	pBlockage	The percentage blockage of the cells. For example, if 40 is entered (i.e. 40%), the cell sides are reduced in flow width by 40% (i.e. is set to 0.6 times the full flow width).	Float
7	FLC_or_FLCpm_below_Obv	Form Loss Coefficient (FLC) to be applied below the FC obvert. Used for modelling fine-scale “micro” contraction/expansion losses not picked up by the change in the 2D domain’s velocity patterns (e.g. bridge pier losses, vena-contracta losses, 3rd (vertical) dimension etc.). This parameter should be used as a calibration parameter. Note: So that this attribute is independent of 2D cell size it has different treatment depending on the object it is attached to as follows: Line: For thin lines, the FLC value is applied to the cell sides unchanged. For thick (whole cell) lines, the FLC value is divided by two (two cell sides in the direction of flow). For wide lines the FLC value is divided by the number of cells across the line (i.e. the line’s width divided by the cell size) and applied to all cell-sides. Polygon: FLC is the form loss per metre length (in the predominant direction of flow). FLC values are not dependent on the flow width, but are on the length of travel in the direction of flow. However, if a negative FLC value is specified, the absolute value is taken and applied unadjusted to all cell-sides affected by the shape. Note, this is not cell size independent, and if the 2D cell size is changed, this attribute also needs to be changed. The form loss coefficient is applied as an energy loss based on the dynamic head Equation 7.30 where \(\zeta_{a}\) is the FLC value.	Float
8	FLC_or_FLCpm_above_Obv	Form Loss Coefficient (FLC) to be applied above the FC obvert. See FLC_below_Obvert attribute above for more information.	Float
9	Mannings_n	For box culverts (BC), the Manning’s n of the culverts (typically 0.011 to 0.015) should be specified. This value prevails over any other bed resistance values irrespective of where in the .tgc file they occur (the exception is if another FC BC object overrides this one). If set to less than 0.001, a default value of 0.013 is used. For bridge decks (BD), can be used to introduce additional flow resistance once the upstream water level reaches the bridge deck obvert (or soffit). For floating decks (FD) this is always the case as the deck soffit is permanently submerged. The additional flow resistance is modelled as an increase in bed resistance by increasing the wetted perimeter at the cell mid-side by a factor equal to (2.*Bed_n)/FC_n. For example, if the FC Manning’s n and the bed Manning’s n values are the same, the wetted perimeter is doubled, thereby reducing the conveyance and increasing the resistance to flow. This parameter can be used as a calibration parameter to fine-tune the energy losses across a bridge or floating structure. Ignored for blank FC_Type.	Float
10	BC_Width	The width of one BC culvert barrel in metres. For example, if there are 10 x 1.8metre wide culverts in a bank, enter a value of 1.8. Using metric units (default) the value is in metres. Using Units == US Customary, the value is in feet. Applicable to BC FC_Type only. Not used by other types of FCs.	Float

Table 7.32: Flow Constriction (2d_fc) Attribute Descriptions

No.	Default GIS Attribute Name	Description	Type
1	Type	Secondary flag identifier where: Blank for general (does not include allowances for any vertical walls or friction from underside of deck). “BC” for Box Culverts (BC is only available if Manning’s n option is specified in Bed Resistance Values). “BD” for Bridge Deck. “FD” for Floating Bridge Deck.	Char(2)
2	Invert	Invert of constriction (m above datum). Mandatory for box culverts (type = “BC”). If not a box culvert, and you wish to leave the Zpt levels unchanged (i.e. use the existing Zpt elevations), enter a value greater than the obvert level (see next attribute).	Float
3	Obvert_or_BC_height	Type: Blank or “BD”: Obvert of constriction (m above datum) “BC”: Height of box culvert (m). Values less than 0.01 are set to 0.01. “FD”: Floating depth (m) of the deck (i.e. depth below the water line). Enter a sufficiently high value (e.g. 99999) if there is no obvert constriction. If using Units == US Customary, the height unit is in feet.”	Float
4	u_width_factor	Flow width constriction factor in the X-direction (i.e. the flow width perpendicular to the X-direction). For example, a value of 0.6 sets the flow width at the left hand and right-hand sides of the cell to 60% of the cell width. Values less than 0.001 are set to 1. Use a value of 1.0 to leave the flow width unchanged. Values greater than 1 can be specified.	Float
5	v_width_factor	Width constriction factor in the Y-direction. See description above for u_width_factor.	Float
6	Add_form_loss	Form loss coefficient. Used for modelling fine-scale “micro” contraction/expansion losses not picked up by the change in the 2D domain’s velocity patterns (e.g. bridge pier losses, vena-contracta losses, 3\(^{rd}\) (vertical) dimension etc.). Can be used as a calibration parameter. The form loss coefficient is applied as an energy loss based on the dynamic head Equation 7.30 where \(\zeta_{a}\) is the add_form_loss value. The form loss coefficient is applied 50/50 to the right and left sides (u-points) of the cell, and similarly to the v-points.	Float
7	Mannings_n	For box culverts (BC), the Manning’s n of the culverts (typically 0.011 to 0.015) should be specified. This overwrites any previously specified Manning’s n values at the cell’s mid-sides. If set to less than 0.001, a default value of 0.013 is used. For bridge decks (BD), can be used to introduce additional flow resistance once the upstream water level reaches the bridge deck obvert or soffit. For floating decks (FD) this is always the case as the deck soffit is permanently submerged. The additional flow resistance is modelled as an increase in bed resistance by increasing the wetted perimeter at the cell’s mid-sides by a factor equal to (2.*Bed_n)/FC_n. For example, if the FC Manning’s n and the bed Manning’s n values are the same, the wetted perimeter is doubled, thereby reducing the conveyance and increasing the resistance to flow. To be used as a calibration parameter to fine-tune the energy losses across a bridge or floating structure. Ignored for “Blank” type FC.	Float
8	No_walls_or_neg_width	Number of vertical walls per grid cell. If set to zero (between 0.001 and 0.001) one vertical wall per cell is used. A non-integer value can be specified. Alternatively, and more easily, specify the width of one culvert in metres by using a negative value. For example, if the culverts are 1.8m wide, enter a value of 1.8 and the number of vertical walls per cell is automatically calculated. Applicable to Box Culverts only. Not used by other types of FCs.	Float
9	Blocked_sides	Indicates whether any of the walls of the constricted cell(s) are blocked off (i.e. no flow across/through the side wall). Specify one or more of the following letters in any order with in the field to indicated which wall(s) are blocked: “R” - block right hand side wall “L” - block left hand side wall “T” - block top side wall “B” - block bottom side wall Note: the quotes should be omitted.	Char(10)
10	Invert_2	leave blank (not used as yet)	Char(10)
11	Obvert_2	leave blank (not used as yet)	Char(10)
12	Comment	Optional field for entering comments. Not used.	Char(250)

7.4.6.2 Applying Form (Energy) Loss in 2d_fc and 2d_fcsh Layers

Important

2d_fcsh and 2d_fc are legacy structures and are supported in TUFLOW Classic only. It is recommended to use the layered flow constriction shape (Section 7.2.9.3) instead.

The following should be noted when adapting structure loss coefficients from a 1D model or from coefficients that apply across the entire waterway, for example, from “Hydraulics of Bridge Waterways” (Bradley, 1978) or “Guide to Bridge Technology Part 8, Hydraulic Design of Waterway Structures” (Austroads, 2018):

• The TUFLOW 2D solution automatically predicts the majority of “macro” losses due to the expansion and contraction of water through a constriction, or around a bend, provided the resolution of the grid is sufficiently fine (Barton, 2001; Syme, 2001).
• Where the 2D model is not of fine enough resolution to simulate the “micro” losses (e.g. from bridge piers, vena contracta, losses in the vertical (3\(^{rd}\)) dimension), additional form loss coefficients and/or modifications to the cells widths and flow height need to be added. This can be done by using flow constrictions (FC cells). Additional form loss can also be added using Read GIS FLC or the command Read Grid FLC.
• The additional or “micro” losses, which may be derived from information in publications, such as Hydraulics of Bridge Waterways, need to be either:
  • Distributed evenly over the FC cells across the waterway by dividing the overall additional loss coefficient by the number of cells (in the direction of flow); or
  • Assigned unevenly (e.g. more at the cells with the bridge piers), however, the total of the loss coefficients should be equivalent to the required overall additional loss coefficient.
• The head loss across key structures should be reviewed, and if necessary, benchmarked against other methods. Note that a well-designed 2D model will be more accurate than a 1D model, provided that any “micro” losses are incorporated.
• Ultimately the best approach is to calibrate the structure through adjustment of the additional “micro” losses – but this, of course, requires good calibration data!

An example of how to apply 2D FCs and a 2D FCSH to a bridge structure is shown in Figure 7.55 and Figure 7.56. The loss coefficient quoted in the figure is an example. Every structure is invariably different, and as such require structure specific parameters!

When applying FCs, the best approach is to view the structure as a collection of 2D cells representing the whole structure, rather than being too specific about the representation of each individual cell. A good approach is to use a 2d_po layer to extract time histories of the water levels upstream and downstream of the structure and of the flow and flow area upstream, downstream and through the structure (see Section 11.3.2).

Of particular importance is to check that the correct flow area through the structure is being modelled. Digitise a 2d_po QA line through the structure from bank to bank and use this output to cross-check the flow area of the 2D FC cells is appropriate (the QA line will take into account any adjustments to the 2D cells due to FC obverts and changes to the cell side flow widths). If the overall structure flow area is not correct, then the velocities within the structure will not be correct and the energy losses due to the changes in velocity direction and magnitude and additional form losses will not be well modelled.

Figure 7.55: Setting FC Parameters for a Bridge Structure

Figure 7.56: Setting FCSH Parameters for a Bridge Structure

7.4.7 Unsupported Features in TUFLOW Classic

The following features are not supported in TUFLOW Classic. TUFLOW HPC is required to access these features:

• Parallelisation / utilising GPU Hardware (Section 1.1.5.3)
• Wu turbulence approach (Section 7.3.2.4)
• Topography Sub-grid Sampling (Section 7.3.3)
• Quadtree grid refinement (Section 7.3.1)
• Groundwater with multiple sub-surface layers and horizontal hydraulic conductivity (Section 7.3.5.2)
• Non-Newtonian Flow (Section 7.3.6)

Austroads. (2018). Guide to Bridge Technology Part 8, Hydraulic Design of Waterway Structures. https://austroads.gov.au/

Barton, C. (2001). Flow Through an Abrupt Constriction – 2D Hydrodynamic Model Performance and Influence of Spatial Resolution (C. Barton, Ed.) [M.Eng.Sc Master Thesis]. https://tuflow.com/media/4983/2001-flow-through-an-abrupt-constriction-2d-hydrodynamic-performance-and-influence-of-spatial-resolution-barton-thesis.pdf

Boyte, C. (2014). The Application of Direct Rainfall Models as Hydrologic Models Incorporating Hydraulic Resistance at Shallow Depths [B.Eng Thesis]. https://tuflow.com/media/5009/2014-the-application-of-direct-rainfall-models-as-hydrologic-models-incorporating-hydraulic-resistance-at-shallow-depths-boyte-thesis.pdf

Bradley, J. (1978). Hydraulics of Bridge Waterways (HDS 1, FIWA). Bridge Division. https://www.fhwa.dot.gov/engineering/hydraulics/pubs/hds1.pdf

Collecutt, G., Baeumer, U., Gao, S., & Syme, W. (2022). Bridge Deck Afflux Modelling – Benchmarking of CFD and SWE Codes to Real-World Data. Hydrology & Water Resources Symposium. https://tuflow.com/media/7554/2022-bridge-deck-afflux-modelling-benchmarking-of-cfd-and-swe-codes-to-real-world-data-collecutt-et-al-hwrs.pdf

Collecutt, G., Gao, S., & Syme, W. (2020). Mesh-Size Insensitive Turbulence Modelling for the 2D Shallow Water Equations. River Flow 2020, pp. 1–9. https://tuflow.com/media/5023/2020-mesh-size-insensitive-turbulence-modelling-for-the-2d-shallow-water-equations-collecutt-et-al-iahr-river-flow-delft.pdf

Collecutt, G., & Syme, B. (2017). Experimental Benchmarking of Mesh Size and Time-Step Convergence for a 1st and 2nd Order SWE Finite Volume Scheme. Proceedings of the 37th IAHR World Congress. https://downloads.tuflow.com/_archive/Publications/07-31-1369.pdf

Schaap, M. G. (2002). Rosetta: A Computer Program for Estimating Soil Hydraulic Parameters with Hierarchical Pedotransfer Functions. University of Arizona. https://cales.arizona.edu/research/rosetta/download/rosetta.pdf

Stelling, G. (1984). On the Construction of Computational Methods for Shallow Water Flow Problems. Rijkswaterstaat Communications.

Syme, W. (2001). TUFLOW – Two & One-Dimensional Unsteady FLOW Software for Rivers, Estuaries and Coastal Waters. The Institution of Engineers, Australia 2D Seminar. https://tuflow.com/media/4985/2001-tuflow-two-and-one-dimensional-unsteady-flow-software-for-rivers-estuaries-and-coastal-waters-syme.pdf

Syme, W. (2008). Flooding in Urban Areas - 2D Modelling Approaches for Buildings and Fences. Engineers Australia, 9th National COnference on Hydraulics in Water Enginnering. https://tuflow.com/media/4997/2008-flooding-in-urban-areas-2d-modelling-approaches-for-buildings-and-fences-syme-hwe-aus.pdf

USDA. (1986). Urban Hydrology for Small Watersheds Technical Release 55 (TR-55). https://nationalstormwater.com/wp-content/uploads/2020/07/Urban-Hydrology-for-Small-Watersheds-TR-55.pdf

Van Genuchten, M. Th. (1980). A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal, 44(5), 892–898.

Wu, W., Shields Jr., F. D., Bennett, S. J., & Wang, S. S. Y. (2005). A Depth-Averaged Two-Dimensional Model for Flow, Sediment Transport, and Bed Topography in Curved Channels with Riparian Vegetation. Water Resources Research, 41(3). https://doi.org/10.1029/2004WR003730

---

1. Figure courtesy of University of Texas http://caee.webhost.utexas.edu/prof/mckinney/ce374k/Overheads/10-Infiltration.pdf#page=20