The 2D ADI solver

The Alternating Direction Implicit (ADI) solver has a long track record in modelling shallow water hydraulics. Developed in the 1980s, when computers were less powerful, the ADI method is still popular because of its speed and robustness.

The method works by representing the model domain as a grid of square cells. Water levels are calculated at each cell centre, and the two components of velocity at cell edges. This allows the model to use the velocity components to calculate flow across cell edges and between cells.

Topography is represented as:

• The ground height at the centre of each cell. This height is used to determine when a cell becomes wet, for example.
• The ground height at each cell edge. This is used to calculate flows across the edge, between cells.

These heights are sampled from the DTM.

Since the model stores water depths at cell centres, and flows between cells are calculated at cell edges, the model must interpolate to determine water depths at these edges. This is actually harder than it sounds: we need to use either the average depth, or the average water level, from the two neighbouring cells, depending on the topography between them. Several methods can be used to determine these cell edge depths (see Figure 8.2). Method C is selected when water depths are below a threshold level (known as calculation depth), with a default value of 0.5m. This allows water to flow over obstacles correctly when the other interpolation methods may produce unrealistic results.

The spatial derivative terms in the shallow water equations are represented using finite differences , for example:

Due to the staggered grid, the solver is second order accurate in space, giving accurate solutions for flows where spatial variations are smooth, but causing possible oscillations where sudden changes in water elevation and velocity occur. The ADI solver may, therefore, not be suitable for trans-critical flows with hydraulic jumps. Even if it does produce a stable solution, this may not represent properly conservation of mass and momentum in the region of hydraulic jumps. For calculation of derivatives used in the advective acceleration term, upwinding is used to bias the estimated derivative in an upstream direction. This improves model stability.

The ADI solver solves the discretised equations by sub-dividing the computation at each time-step into x- and y-directions. On the first half time-step, the water depth and the flows between cells in the x-direction are solved implicitly, whilst the other variables are represented explicitly. Similarly, for the second half time-step, the water depth and the flows in the y-direction are solved implicitly, with other variables represented explicitly. The 2D solver problem is thus represented as a series of 1D calculations, which can be solved efficiently, both in terms of processor time and memory.

The 1D calculations in each direction are non-linear, and so the equations are linearised before being solved. For example, a square term may be approximated by:

This allows a linear solution method (Gaussian elimination and back substitution in 2D solver) to be used. Accuracy is improved by performing a number of sub-iterations (typically two, controlled by the Nadvit parameter) at each time-step, where the value produced by the last iteration is fed into the new iteration instead of the value from the previous time-step:

Linearisation is optionally applied to the friction term in the momentum equations:

The rationale behind this treatment is that as depths approach zero, the friction term becomes large and can destabilise the equations. By linearising, we ensure that velocity goes to zero as depth decreases. This linearisation is controlled by the w ater depth (friction) parameter, and is only applied below that depth. By default, water depth (friction ) is large (100m), meaning the linearisation is always applied.

Time-stepping is treated implicitly and, because of the two half time-steps used in each calculation direction, is second order in time. This allows Courant numbers up to 8 to be used.

Modifications to the advection term

As noted in the previous section, supercritical flows are harder to solve than subcritical. This is because the advection term in the shallow water equations causes instabilities in the model solution. Two modifications to the ADI solution method are used to improve stability:

1. For Froude numbers greater than 0.75, the advection term is reduced as the Froude number increases, becoming zero when the Froude number reaches one. This improves stability as flows become supercritical, but as Froude numbers approach 1, the solution is less accurate because of the reduction of the advection term.
2. An advection capping algorithm, designed to detect and fix problems caused by rapid changes in velocity. The advection term represents a spatial gradient of the kinetic energy of the flow:

The following limit is placed on the advection term:

For decelerating water, this has the effect of ensuring the advection term is sufficient to reduce velocity to zero (ie to develop a stagnation point), but no larger, reducing the likelihood of large oscillations. For accelerating flows, the advection term can double over each cell but no more. The advection cap is applied only to flows where the velocity head (v/2g ) is greater than the specified threshold (in the 2D model file). A default value for this parameter is applied in the 2D Simulation user interface, which can then be edited by the user if required.

Choose a very high value (>=5) to switch off advection capping; choose a lower value where accelerating flows are causing instabilities. However, note that a lower value may reduce the accuracy of the model output in the areas of your model with accelerating flows. If this area is an important part of your model then you could consider using the TVD solver instead (which tends to perform better for accelerating flows).

Wetting and Drying

Wetting and drying of cells is treated by keeping track of which cells are wet (in which case they are used in the ADI shallow water computation) and dry (cells are ignored in the computation). The wet/dry cells are updated as shown in Figure below. If a wet cell is next to a dry cell, and the water level is more than a depth threshold (known as w ater depth (dry)) above ground level in the dry cell, then flow is switched on between the two cells, allowing the dry cell to become wet at the next time-step. For drying, a cell is switched off if its depth drops below w ater depth (dry).

This approach can lead to a negative depth in some cells, especially for drying cells where there is flow down a steep slope. These are detected by 2D solver and corrected, where possible, by moving small amounts of water from neighbouring wet cells. This reduces the impact of these negative depths on mass balance.

A further modification to the solution method is made in very shallow water, to help prevent negative depths and improve model stability. For depths below a threshold (denoted d epth threshold ), the friction term in the shallow water equations is increased, to a maximum of twice the normal Manning friction, when water depth is equal to w ater depth (dry). This has the effect of reducing velocities as a cell dries, encouraging the cell to dry gradually, reducing the likelihood of negative depths.

Direct Rainfall Modelling

The solution method includes special methods for modelling direct rainfall.

Direct rainfall modelling works by adding an amount of water to each cell at each time-step, which varies with the time-step and rainfall rate. For small time-steps and low rainfall intensities, the amount added can be smaller than the floating point precision (typically 1 part in 10 million for single precision arithmetic) with which the water level is stored. The problem is worse for models with higher ground elevations, because the water level is stored for each cell, rather than the depth. Low rainfall rates may, therefore, cause mass balance errors or even no flooding at all.

This is unlikely to be a problem for double precision models. However if running in single precision, it can be an issue if the rainfall rate R, in mm/h, is (approximately) less than 0.5z/Δt, where z is a typical ground elevation. For example, if z≈100 m and the time-step is 1s, any rainfall with an intensity of less than 50mm/h may be subject to significant mass balance errors.

2D solver modifies the method to fix this issue by adding the rainfall in 'pulses', large enough to be represented at single precision. Water is not added at every time-step, but is saved up until enough can be added and not cause precision errors. Typically this pulse depth is 0.1 mm, and can be controlled by the user through the rainfall accumulation depth parameter. Adding rainfall in pulses, rather than at each time-step, makes very little or no difference to the model outputs, apart from fixing the mass balance errors.

A second modification addresses a problem that occurs because of the depth threshold required to activate a cell as wet. At the start of the simulation, when enough rainfall has accumulated to activate cells (typically 10 mm), all cells then become active. This can cause model instabilities and an inaccurate representation of catchment response.To fix this, a kinematic routing method moves water between cells at low depths, even when the depth is below the water depth (dry) threshold. This allows water to accumulate in low areas and potentially activate cells as wet, earlier in the storm. It also improves model stability.