25/07/2017 11:46:41

*by Konrad Adams, Senior Developer for Flood Modeller Pro at Jacobs*

### Why should we care?

One of the fundamental equations solved in the Shallow Water Equations is the conservation of mass (the other(s) being conservation of momentum in however many dimensions you are solving it in). It also stands to reason – the volume you put in the model should equal the volume that leaves the model plus the additional volume stored in the model.

Put simply, if your volume is not conserving mass adequately, how can you possibly believe any of the results you are getting?

New for Flood Modeller version 4.3 is the inclusion of the mass balance display in both the 1D run-time information and as a summary in the diagnostics file.

### Why is it not zero?

In an ideal world, your mass error (essentially, the difference between your net inflow volume and your model storage volume increase) would be zero. Unfortunately, this is not an ideal world; or let me turn this around and be thankful that there are jobs for modellers, reviewers, etc. to lend our expertise and knowledge to such issues. Since the shallow water equations are solved numerically, and even that is solved approximately, there is no guarantee that the mass (or indeed momentum) equation is solved identically; what we need to do is to solve it as close to exactly as possible.

There are also some approximations used in the discretization of the volume calculations and the flux calculations, which could also lead to minor differences.

There are also further modelling and numerical tricks we can use to ensure we do have a solution, which could compromise the integrity to an extent. However, a bad result is about as much use as no result at all; therefore, one of our jobs is to ensure that it is solved “well enough”, i.e. the integrity is not compromised detrimentally and the results are representative of reality. One of those criteria is to ensure that mass is conserved to within our accepted tolerance.

### What do we (not) show?

There are several ways in which we could display the mass error. We could show it as an absolute value, either cumulatively, or an a per timestep basis. However, this needs to be set against the context and size of the model; a large model with a small relative mass error would show a larger absolute error against a small model with a high relative mass error.

In order to be consistent across different models, we show a percentage mass error… but even this is not straightforward – as a percentage of what? Current volume? Maximum volume? Change in volume? The volume passing into the model at boundaries?

We have chosen to show the relative error as the difference in mass relative to the maximum system volume, both as the simulation progresses and as a final summary.

We also show an alternative summary final mass error, as relative to the total inflow volume (throughput).

### How is it calculated

**Volume**

There are two type of units which store water in Flood Modeller 1D – channel sections (RIVERs and CONDUITs) and RESERVOIRs / PONDs. The structures, junctions, etc. are assumed to have zero length, and hence contain no volume.

The storage volume between two adjacent sections is calculated as the mean wetted cross-sectional area multiplied by the distance between them; hence the total system *channel* volume, *V _{CHAN} * is given by:

Where *A _{u}* and

*A*are the upstream and downstream wetted areas of each consecutive cross-section pair, and Δ

_{d}*x*is the distance between them.

Storage volume in a reservoir or pond is calculated from the area-stage table – similarly, this is the sum of mean surface area multiplied by the vertical distance between each data pair, over all pairs up to and including the current water level (the current surface area being calculated by interpolation of surface areas in the table); hence the total system *reservoir* volume, *V _{RES}* is given by:

Where *A*_{i }and *h _{i}* are surface area-elevation pairs in the reservoir data table, and the inner summation is over all elevations up to the reservoir water level.

The total system volume at any time is therefore the sum of the volumes of all channel and reservoir units (except for routing sections – see note below); i.e. *V _{CHAN }*+

*V*

_{RES}.**Net inflow**

The boundary flow volumes, *V _{IN}* and

*V*, are calculated as:

_{OUT}- For inflow, the cumulative sum of inflow rate (in m
^{3}/s [or ft^{3}/s]), multiplied by the timestep, over all inflow boundaries (e.g. QTBDY units, etc.) and timesteps in the simulation

Where Δ*t* is the model timestep, and the inner summation is over all model boundaries flowing into the model.

- Conversely for the outflow, the cumulative sum of outflow rate (in m
^{3}/s [or ft^{3}/s]), multiplied by the timestep, over all outflow boundaries (e.g. HTBDY units, etc.) and timesteps in the simulation

Where Δ*t* is the model timestep, and the inner summation is over all model boundaries flowing out of the model.

Note: a volume entering the system is *always* counted as an inflow, and *vice versa*, e.g. in a typical model schematisation, a QTBDY with negative flow is decreed an outflow, and an inflowing tidal HTBDY (with *Q<0)*, is decreed an inflow.

### 2D / linked flow

Link volume between 1D and 2D model components is treated in a similar way to the boundary inflows/outflows, i.e. the cumulative sum of exchange flow rate (between 1D and 2D) multiplied by the timestep.

Again, irrespective of whether it is an HX/H-link or SX/Q-link, if the flow is into the 1D model, it contributes to the link *inflow* volume, and vice versa for the outflow. This is supplied separately from the **boundary** flow volumes and is included within the final net inflow.

### Routing

Muskingum routing (reaches of) models calculate discharge, but do not always explicitly contain cross-sectional information in order to calculate storage volume; or put another way, the mass balance / storage is an integral component of the routing model. Therefore, these are implicitly assumed to be mass conservative and excluded from the output calculations. Thus, the inflow into routing reaches is removed from the system inflow, and the outflow from routing reaches is added to the system inflow. This means effectively the routing sections are removed from any mass balance calculations, which is thus confined to the fully hydrodynamic component of the model.

Therefore, the total net inflow volume, *V _{NET} *is given by the total inflow – total outflow, excluding routing sections, i.e.

Where *V _{LINKIN }* and

*V*are the link inflow and outflow volumes, respectively, and

_{LINKOUT}*V*

_{MUSKIN}and

*V*are the respective routing inflow and outflow volumes.

_{MUSKOUT}### Benefits and drawbacks of each measure

- Consider a typical event-based simulation, where we start from and return to baseflow conditions, having run a large hydrograph through the system. The difference between initial and final system volume will therefore be close to zero. However, due to the reasons stated above, there may be a discrepancy, however slight, between net inflow volume and system volume. If this were expressed as a percentage of
*net volume increase*, then the percentage mass error would be misleadingly high. - A constant boundary (steady-state) model – in this case we are simulating constant boundary conditions, using an unsteady simulation, in order to arrive at a steady-state solution. In this case, the maximum increase in volume should be negligible; any volume discrepancy expressed as a proportion of
*maximum*volume increase again would be exaggerated.

There is no one-size-fits-all solution, but for the best all-round picture, we express the % mass error as a proportion of *maximum system volume*. Again, this is not perfect; for a start, the maximum system volume itself varies until the peak is reached (so the same absolute error would give a larger % error in the rising limb than on the falling limb). Also, this may give over-optimistic representations of percentage error, especially if the initial volume is large compared with the increase in volume. However, in typical cases of flood modelling simulations (low initial conditions; large simulated event) this should be sufficient, although the numbers presented should always be taken in the correct context.

One issue with the above can occur with cyclical models. This includes tidal models, and I will also include multiple-event models under this umbrella. For instance, if one performs a single cycle and the reported mass error is 0.1% (volume discrepancy = maximum volume/1000), then this looks good. If we repeat the cycle to a total of ten cycles, then the maximum system volume does not increase, but assuming the same mass error accumulates, then 1% is not looking as good… then after 100 cycles a *reported* 10% error is starting to look quite bad, even though the model is returning to its initial state after each cycle (and the actual error is nowhere near as bad as it looks). In this case, the percentage error is better expressed as a proportion of the volume passing through the boundaries (which does increase each cycle).

The above was from a real situation during testing, with a constant upstream inflow and a repeated tidal boundary (see output in Figure 1 and Figure 2). In this case, the error expressed as a % of throughput volume was deemed more reasonable, and the model results did show almost exact repetition during each cycle (as opposed to any cumulative change).

### Things that affect and improve mass conservation

The following have been identified as potentially having a significant effect on mass balance:

**Parameters**

*Alpha / Minitr*

The purist in me is not a great fan of the *concept* of the parameter alpha (α), but if you try removing it, you tend to appreciate its worth a bit more. This is the parameter which weights the calculated solution a certain proportion (1-α) towards the previous solution. Effectively what it says is, you (or at least the solver does) go to the trouble of calculating a solution, but you don’t use it 100%, and you need to take a proportion of the previously calculated solution. The practical effect of this is that it smooths the solution between iterations, reduces the chances/effect of sudden spikes, and gives more chance of obtaining a valid solution. Now because we are replacing a chunk of our calculated solution with a chunk calculated at some other time, this has the potential to be non-mass-conserving. The more solution iterations performed, the more the effect of α is diluted; therefore if you reduce α, it is always wise to increase the minimum number of iterations, *minitr*, to say 3 or 4 (from its default of 2).

*Timestep*

A smaller timestep *usually* increases the accuracy of calculations, both in the overall simulation and also specifically the mass balance calculations (although see later notes on *Numerical precision*). This stands to reason – we’re increasing the resolution of our solution, so would expect a higher accuracy. However, it also has an impact on the model simulation time, so the trick here is to achieve an appropriate compromise. We may do this by altering the timestep and checking the solution is still the same. The mass error calculation will give a very quick idea whether this is so, although other key indicators (e.g. the results themselves) will also need checking.

*Tolerance*

The tolerance parameters **htol** (head) and **qtol** (flow) determine how “good” the numerical solution needs to be before being considered as solved (the smaller the number, the tighter the tolerance). If these parameters are increased (from their respective defaults of 0.01), then generally the less accurate the solution is – which may manifest itself as a mass error. However, decrease them and the solution may take more time to converge (or not converge at all) and *may* not improve the solution significantly. Usually the default values are sufficient, but increasing them is not generally recommended.

*Numerical precision*

The default numerical engine uses Single Precision arithmetic in its calculations. This means that the numbers used in the computations are accurate to about 7 significant figures. Usually this is more than enough, although there are at least two typical instances where this is not the case, and Double Precision arithmetic (approximately 15 significant figures) gives a significant improvement. These are:

- Rainfall modelling, where the typical resultant inflows, expressed in m
^{3}/s are relatively small - Modelling reservoirs, particularly if [most of] the following apply: high elevation [~100m AD], large surface area, small timestep. Again, the relatively small increase in water level per timestep (in m) can be lost alongside the large volumes using single precision arithmetic.

Because of the size of numbers being stored, double precision simulations will use more memory, and consequently may take more time. However, the mass error may be significantly improved, especially in the above cases.

*Spatially varying density*

So far, I have used the terms mass and volume interchangeably – and one can safely do this in most simulations, where we are simulating water flow under normal atmospheric conditions. However, it is possible in Flood Modeller 1D to apply a density gradient (spatially varying density), in order to simulate, primarily, the higher density (due to salinity) within tidal models towards the tidal end of the model. In the hydraulic model alone (i.e. not the WQ model), this is a simplification, and does not take into account transport of the salt. However, the “mass” error is actually quoted as a volume, whereas the mass equation *does* solve conservation of mass. Hence, if this feature is used, a larger mass error may be reported.

### Summary

In summary, the mass balance output is a useful tool to help us assess model performance and identify when it may need improving.

One can view this as a single final summary value, or as a time series during the progress of the simulation (see Figure 3), which were previously difficult to obtain without performing manual calculations.

The final summary volume figures are output to the diagnostics (*.zzd) file, as well as an alternative calculation of the % error, which give the numbers more context and further assist the modeller/reviewer in ascertaining whether this is a cause of concern.

The obvious question is *What is an acceptable mass balance error*? And it’s a very difficult one to answer generally. As we have seen above, the number needs to be viewed in context, and it very much depends on what our goal is and what the mechanisms for the error are. It’s probably fair to say that numbers approaching and over 10% would definitely warrant further investigation; conversely, since we’re generally happy with a relative tolerance (*qtol*) of 0.01, which normally equates to a flow tolerance of 1%, then maybe we shouldn’t be too expectant of a mass error any less than that (although of course, this can often be achieved). In between those two values, there are many shades of grey, and almost certainly an area for further discussion. However, it should be viewed as one of many possible indicators of model quality, and others, such as convergence statistics and the results themselves should not be ignored.