Tracing and visualisation of contributing water sources in the LISFLOOD-FP model of flood inundation (within CAESAR-Lisflood version 1.9j-WS)

. We describe the formulation of a simple method of water source tracing for computational models of flood inundation and demonstrate its implementation within CAESAR-Lisflood. Water source tracing can provide additional insight into flood dynamics by accounting for flow pathways of each model boundary condition. The method developed is independent of the hydraulic formulation used, allowing it to be implemented in other model codes without affecting flow routing. In addition, we developed a method which allows up to three water sources to be visualised in RGB colour-space, while continuing to allow 5 depth to be resolved. The number of water sources that may be traced is limited only by computational resources. We show the application of the methods developed for example applications of a major flood event, a shallow estuary, and Amazonian wetland inundation. A key advantage is that the method is independent of the hydraulic formulation, meaning that it is relatively straightforward to add to existing finite volume codes including those based on or developed around the LISFLOOD-FP method. This method enables water tracing with a minimal computational overhead, allowing users of the LISFLOOD-FP 10 method to address environmental issues relating to water sources and mixing, such as water quality and contamination problems.


Introduction
Flood inundation models relate upstream river inflow and other dynamic boundary conditions (e.g.downstream stage) to inundation extent, depth and velocity, and have become invaluable tools for the assessment and understanding of flood dynamics and risk.Recently a series of fast and effective two dimensional hydrodynamic models have been developed (Hunter et al., 2007;Neal et al., 2018).These model codes have enabled assessment of flood inundation at high spatial resolutions (e.g.Yu and Coulthard, 2015) and large spatial scales including applications on the Amazon (Wilson et al., 2007) and Congo (O'Loughlin et al., 2020), at continental scales (Dottori et al., 2021;Wing et al., 2017) and global (Dottori et al., 2016;Sampson et al., 2015).The LISFLOOD-FP model (Bates and De Roo, 2000;Bates et al., 2010) uses a fixed raster grid structure and includes several flow formulations of different physical complexity (Shaw et al., 2021), which have been widely used for flood hazard assessment in fluvial (Horritt et al., 2010) and coastal (Vousdoukas et al., 2018) applications, as well as within models of landscape evolution (Adams et al., 2017;Coulthard et al., 2013).However, the ability to trace the sources of water in the model domain is presently missing from reduced-complexity 2D flood models such as LISFLOOD-FP.
Understanding the contribution of different water sources to flooding and river flows is important when managing river basins, for example determining the relative contribution of tributaries or where water borne contamination is an issue.This 1 https://doi.org/10.5194/gmd-2021-340Preprint.Discussion started: 18 October 2021 c Author(s) 2021.CC BY 4.0 License. is a function found in more complex two-and three-dimensional models based on the full shallow water equations, such as TELEMAC (Galland et al., 1991) and Ansys Fluent, and has been used for a variety of applications including water quality modeling in a lake (Kopmann and Markofsky, 2000) and solute and viral dispersal in estuaries (Robins et al., 2014(Robins et al., , 2019)).
For example in TELEMAC, non-buoyant tracers can be added and their course and concentration followed (Ch.9, Ata et al., 2014).However, such 2 and 3D codes have a considerable computational overhead compared to simpler schemes such as LISFLOOD-FP.
In this paper, we build on the work of Wilson and Coulthard (2019) to propose a simple and efficient method which can be added to finite volume codes to track the contribution of multiple water sources to predicted flood depths, and demonstrate it within the CAESAR-Lisflood model (Coulthard et al., 2013) for several example flooding case studies which each represent the mixing of water from multiple sources.CAESAR-Lisflood 1 is a development of the CAESAR model (Coulthard et al., 2002;Van De Wiel et al., 2007) and the LISFLOOD-FP 2D hydrodynamic model (Bates et al., 2010) that simulates landscape evolution by coupling a hydrological model, a surface water flow model, fluvial erosion and deposition and slope processes.
Here, we focus only on surface water, and additionally develop a method for the visualisation of water sources.

Flow formulation
CAESAR-Lisflood implements the LISFLOOD-FP inertial formulation (Bates et al., 2010) to estimate depths across the domain.This is a reduced physical complexity mass-conservative hydraulic model structured on a Cartesian grid, in which the updated water depth h, at time t, for cell i, j is determined using (Bates and De Roo, 2000): where h t−∆t i,j is the cell depth at the end of the previous timestep (∆t), Q represents the flows into or out of the cell in the x or y directions, and ∆x is the cell size.Flows are decoupled in each direction; flow in the x direction is determined using (de Almeida et al., 2012;Bates et al., 2010): with flow in the y direction obtained analogously.In Eqn.
(2), q t−∆t x represents the flux between the cells from the previous timestep (Q t−∆t x /∆y), g is acceleration from gravity, h [flow] is the maximum depth of flow between two cells, n is Manning's roughness and z is the cell bed elevation.The LISFLOOD-FP inertial formulation has been benchmarked against other formulations and industry standard codes and showed that, in appropriate flow conditions (i.e.gradually varied flow, Froude number < ∼1), the model performed favourably and with high efficiency (de Almeida and Bates, 2013;Neal et al., 2012).

Water source tracing
Given that the sum of depths in a cell from each source, w, make up its total flood depth, h: the volume fraction of each source, ϕ w , may be obtained from: ensuring: Thus, the fraction of each water source in each cell is defined as the depth of water from that source in the cell, divided by the total cell depth.
Following this, once flows between cells are calculated and depths updated, water tracing proceeds as follows: (1) for each cell, the remaining depth following removal of water from any outflows is found; (2) the amount of the depth belonging to each water source is obtained by scaling the depths according to water source fractions from the previous timestep; (3) contributing inflow depths from each source for all neighbours is added to each source depth; and (4) updated water source fractions are calculated as per Eqn. 4 by dividing the fractions from each source by the updated total cell depth.More formally, for each cell, for each water source, w, at time, t, the volume fraction ϕ is: where represents the depth remaining in the cell after the removal of outflows at the current time, Q t [out] (and before the addition of any inflows), which is scaled according to the fraction from this source from the previous timestep, ϕ t−∆t w .
To this depth is added the total amount of water from this source which is contributed from neighbouring cells, given by D  [in] is the inflow from a particular direction, D, and ϕ t−∆t,D w is the fraction of flow for source w from that direction.The updated source volume fraction is obtained by dividing by the updated cell depth, h t .Finally, h t Q [out] is given by:  [out] are the outflows in each direction.To improve computational efficiency, the method only requires knowledge of the source volume fractions for the previous timestep in each neighbouring, upstream cell; the actual depths or volumes from each source in each cell are not saved.The calculation cell depth after removal of outflows, h t Q [out] , enables each source fraction to be updated based on the flows in (Q t,D [in] ), the fraction from each source in direction, D (ϕ t−∆t,D w ), and the updated depth, h t .
Using this approach, calculation of water source fractions flowing out of the cell is not necessary.For completeness, using the notation of (1) and expanding D into cell indices, i, j, the updated depth from each source can be obtained using: Boundary condition inputs are added prior to the routing of surface water.Adjusting the water volume fractions for the cells in which water is added is straightforward.The boundary modified water volume fractions, ϕ ′ w , are obtained using: for input sources, or: for other sources, where Q w is the inflow added to the cell from source w at the start of timestep.Thus, fractions from sources where water is added to the cell are adjusted upwards, while fractions for non-source volumes are adjusted downwards.
It should be noted that, for simplicity, this method treats cell water volumes as fully mixed.Consequently, it may be possible for small fractions to propagate quickly in a downstream direction, since fluxes into a cell would be included in the fractions assigned as inflow to a downstream neighbouring cell in the next timestep, via ϕ t−∆t,D w .As a result, caution should be given to the interpretation of very small water source fractions.

Water tracing implementation
Importantly, within the model code, the additional equations required for the water source tracing presented in Section 2.2 can be piggybacked on top of existing finite volume model code, requiring minimal modification of the original schemes.
The water source tracing method was implemented in C#, within CAESAR-Lisflood version 1.8f (Coulthard, 2015).However, the method is readily transferable to other programming languages where different versions of LISFLOOD-FP (Bates et al., 2010) are implemented, or to other similar flow models.The additional C# code added to CAESAR-Lisflood is included in Appendix A. Here, we summarise the algorithms developed.
Our method allows water from different sources to be tracked according to point sources (e.g. a reach input) or from spatial areas that can be user-defined (e.g. to represent different rain zones or subcatchments).CAESAR-Lisflood allows the addition https://doi.org/10.5194/gmd-2021-340Preprint.Discussion started: 18 October 2021 c Author(s) 2021.CC BY 4.0 License. of water as point sources, or via rainfall through a hydrological model based on a spatially distributed version of TOPMODEL (Beven and Kirkby, 1979), for implementation see Coulthard and Skinner (2016) and Coulthard and Van De Wiel (2017).The zonal input also allows tidal sources to be traced, although note that the model is depth-averaged and does not account for the higher density of saltwater.Water tracers are saved as two 3D or stacked arrays of size W * X * Y grid cells, for the current and previous timestep, which are initialised as zeros.Additionally, two arrays for δhδt in the x and y directions are used to save the volumes of water moving between cells in units of change in depth in the timestep.
The main functional flow of CAESAR-Lisflood as related to water source tracing is shown in Figure 1.For each time step, inflows from various sources are first added.The addition of traced water requires code to be modified such that, after the addition of water from each source, the water source fractions are updated to account for the increase in volume from a source (Eqn.9 and 10).In CAESAR-Lisflood, inflows to the modelled domain may reflect reach inputs, rainfall inputs and tidal or stage inputs, which are added in three separate functions (reach_water_and_sediment_input(), catchment_water_input_and_hydrology() and stage_tidal_input(), respectively).Each of these are updated in a similar way using Algorithm 1.
Following the addition of inputs from sources, the implementation in CAESAR-Lisflood copies the water depths and water tracers into arrays which represent the previous timestep, h t−1 and ϕ t−1 .For the water source tracing implementation, the tracer states are saved via the function save_tracer_states() (Figure A3).Next, flows are updated in the function qroute(), for which no modifications are required (meaning that the water source tracing does not affect flow).Depths are then updated via depth_update() with no modifications except that volumes of water moving between cells are stored to the δhδt in x and y arrays, calculated using: for a single direction (C# code in Figure A4).This is used as part of the solution to ( 6) and ( 7) to update the water source tracers, enabling a simplified form to be implemented for the tracer update: where: Water source tracers are updated in the function update_tracer_states(), which is the main addition to the CAESAR-Lisflood code.Algorithm 2 provides details of the procedure used.For every cell in the model which contains water (h t i,j > 0.0), the depth remaining in the cell after outflows (h ) is first calculated by subtracting outflows from each of the four flow directions.We sum δhδt from each direction (previously calculated in depth_update()), which represent contributions to https://doi.org/10.5194/gmd-2021-340Preprint.Discussion started: 18 October 2021 c Author(s) 2021.CC BY 4.0 License.
Algorithm 1 Updating water source fractions after adding input sources.Inputs required are the depth in the cell for the previous iteration (h t−1 ), a vector of all water sources for the cell for the previous iteration (ϕ t−1 ), input from this source (Q w ), the index of this source (w), the number of sources (W ), the cell size (∆x) and the timestep (∆t).A C# implementation is shown in Figure A2.depth from the right (x) or up (y): outflows will be negative for the right and up directions and positive in the left and down directions.
The procedure then updates each source for the cell, by first calculating the sum of depths added from this source (i.e. )), again by scanning the δhδt arrays in each direction.Finally, the fraction is then calculated by dividing 135 the total new depth from the source by the depth from all sources, h t .

Water source visualisation
Visualisations of flood model outputs are often produced for depths in the form of animated maps.Here, we derived a simple colour scaling which permits up to three water sources to be visualised in RGB colour-space.For each cell, the RGB colour index in the range 0 to 255 is obtained for the desired water source, w, using: Algorithm 2 Updating water source fractions after flow routing and depth update.Inputs required are the updated depth in the cell (h t ), water sources for the previous iteration (ϕ t−1 ), depth contributions horizontally (δhδtx) and vertically (δhδty), the number of sources (W ), and number of cells in the grid horizontally (X) and vertically (Y ).A C# implementation is shown in Figure A5.
{get depth remaining in cell after removal of outflows} for src = 1 to W do δhδtsumQ [in] ← 0.0 {obtain total change in depth after for inflows of this source (ignoring outflows)} i,j,src ← 0.0 {if there is no contribution or existing water from this source} else The power term, β, enables visual enhancement of lower water fractions, in the range ∼0.1 to 1.0, where 1 would represent no enhancement.In order to allow water depth to be resolved through the use of reduced lightness for deeper water, ( 14) may be modified to: 145 where h i,j is the cell depth and h [range] is the range of depths to scale the visualisation to.Thus, the RGB value for a cell which has only one source of water in it will range from 128 for deeper water where h = h [range] , towards 255 for shallow depths close to zero. Figure 2 provides a detailed illustration of RGB values obtained using ( 14) and ( 15) for different values of ϕ, β and h, and the resulting colour scales obtained when visualising multiple sources.Mixing of water sources with red and blue shading gives a pink colour scale; green and blue gives cyan; red and green give yellow; and shades towards white 150 would indicate mixing of all three sources.Note that, while it is possible to trace an arbitrary number of water sources using the method described in Section 2.2, this visualisation method is limited to a maximum of three water sources.However, using (15), if a cell only contains water from sources different to those selected for visualisation, depths can still be resolved in greyscale in the RGB value range 0 to 127.This also means that it is possible to visualise only one water source at a time in a selected primary colour.

Evaluating performance overhead
In order to test computational performance of the water tracing code, a simple planar test case was developed, consisting of a constant 0.001 m/m slope of length 2000 m and width 1000 m, on a grid with a spatial resolution of 5 m, giving a total of 80,000 cells (Figure 3a).Manning's roughness, n, was set at 0.05.At 250 m intervals down slope, 1 m "walls" were added across the slope, each with between 4 and 8 gaps of 5 m width added to allow water to flow through.This was done to 160 ensure that water mixed as it flowed downslope.Simulations were conducted for a period of 24-hours with constant inflow of 10 m 3 /s to 8 locations at the top of the slope.Eight simulations were conducted in total: one with no tracing and seven with between 2 and 8 water sources being traced.To ensure that each of these simulations had identical total inflow, sources were grouped for tracing purposes, allowing the variation in computation requirements to be assessed.The simulations conducted are summarised in Table 1.Furthermore, to benchmark model performance in a real world scenario, the effect of the number 165 of tracers on model performance was tested on the Carlisle example described below.

Results
Firstly, we demonstrate water source tracing and visualisation for three flood inundation case studies that are examples of flooding at three different spatial scales in three contrasting contexts: (1) a major flood event in 2005 at Carlisle, United Kingdom, fluvial flooding in an urban area which is situated at the confluence of three rivers, (2) a shallow estuary in Christchurch, 170 New Zealand, combining tides with inflow from two small rivers, and (3) flooding at the tributary junction of two large rivers Table 1.Planar slope test case simulations.In each of the eight simulations, the total volume of inflow was identical: each of the eight sources was a constant 10 m 3 /s (where a source was used for multiple inflow points, its total volume was added to all).Inflow points are shown in Figure 3a.on the Amazon floodplain, Brazil.Secondly, we present the computational performance of the method for the planar test case and Carlisle example.

Carlisle, United Kingdom
The city of Carlisle in northern England experienced significant flooding in 2005 (estimated annual exceedance probability 0.57% to 0.5%) that resulted from heavy rains in the headwater catchments of rivers running through the city and affecting more than 2500 properties (Environment Agency, 2006).With the benefit of an extensive set of observational data obtained from field collection, the 2005 flood event at Carlisle event has been used extensively as a test case for hydraulic model development (Horritt et al., 2010;Neal et al., 2009).Here, the site is of particular interest as it is at the confluence of three separate rivers: the main River Eden, which runs from east to west through the city, is joined by the Rivers Petteril and Caldew which flow into the city from the south (Figure 4).
A simulation was run using a model grid of 5 m (domain size: 951 x 612 cells, 14.6 km 2 ), topography from LiDAR and inflows from gauging station records.The simulation began on 08 January 2005, 00:00 AM, for 120 hours (5 days).
Visualisations of the model output are shown in Figure 5. Flood water mixing is shown in RGB colour space, using (15) to enable darker shades to indicate deeper water.The mixing of Eden (blue) and Petteril (red) gives pink shades, which as it moves downstream transitions towards purple due to the larger volume of flow on the Eden; water from the Caldew is shaded in green and contributes to the greyer shades in the floodplains close to the main channel as the three waters are mixed.On the eastern edge of the city of Christchurch, New Zealand, the estuary of the Avon and Heathcote Rivers (Figure 6) is a large (∼8 km 2 ), shallow (mean depth ∼1.4 m at spring high tide and a tidal range of ∼1.7-2.2 m) and intertidal area, with about 85% of its elevation above the spring low tide level (Findlay and Kirk, 1988).The estuary is semi-enclosed and separated  recognised for its importance as a habitat for migratory bird species (Woodley, 2018).The estuary acts as a sediment trap for inflows from the Avon and Heathcote Rivers, the catchments of which cover a combined area of 188 km 2 which is largely urbanised, meaning that urban pollution is a considerable issue (Vopel et al., 2012).
Here, we demonstrate water source tracing and visualisation for the Avon-Heathcote estuary for the month of July, 2017, using a model grid of 10 m (domain size: 483 x 675 cells, 32.6 km 2 ).The sources of water included were the two rivers, with inflow from the Avon River from the north and inflow from the Heathcote from the south-west, and a downstream tide boundary condition at the estuary outlet at the south-east (Figure 6b).The simulated period of July 2017 included both neap and spring tides and a high flow event on 22 July on both rivers, where the Avon River peaked at 13:00 with 14.7 m 3 /s of flow, close to the mean annual flood level (LAWA), and Heathcote River peaked at 14:00 with 28.1 m 3 /s of flow, an Annual Exceedance Probability of less than 0.1 (LAWA).
As noted in Section 2.3, as a depth-averaged model, CAESAR-Lisflood does not account for the higher density of the saline tidal water compared to fresh riverine water, meaning that any variations in the vertical salinity profile of estuarine water is not accounted for.However, given that the Avon-Heathcote estuary receives only small volumes of riverine inflow, with mean inflow from the Avon and Heathcote rivers 1.87 and 1.04 m 3 /s, respectively (LAWA), compared with tidal inflow rates of up to 500-800 m 3 /s during an incoming tide (van der Peet and Measures, 2015), it is likely that it is well-mixed with little vertical salinity difference (Hansen and Rattray, 1966).
Results of the model application are shown in Figure 7.During low flow conditions during at low-tide on 8 July 2017 (image A), the estuary is mostly drained, with the water remaining in the estuary containing a greater proportion of water from the Avon River as a result of its higher baseflow.Shortly after spring-high tide on 16 July (image B), it can be seen that water from both the Avon and Heathcote is pushed back by the tide, particularly along the eastern shoreline.If accurate, this situation may play a significant role in determining where pollutants, such as those from stormwater runoff which are then discharged from the Avon and Heathcote Rivers (e.g., see Vopel et al., 2012) are deposited.This would affect human health through various activities in the estuary, such as the gathering of food sources such as fish and shellfish which accumulate heavy metals contained in riverine discharge (e.g.McMurtrie, 2015).
During the high flow conditions on 22 July shortly after high tide and inflow peak (image C) and low tide (image D), results indicate substantially greater riverine contribution to the estuary, such that, after the fall of the tide, the estuary remained substantially inundated but with water almost all from riverine sources.After high tide, Avon and Heathcote waters are largely confined by the tide to the western shore, with a clear boundary between the two river sources visible.After low tide, this boundary is sharper due to the absence of tidal mixing, and gradually dissipates as water mixes towards the estuary outlet.

Amazon River, Brazil
The flow of the central Amazon River is characterized by a large annual flood wave of around 10 m amplitude, resulting primarily from the distinctive wet and dry seasons (Trigg et al., 2009).Each year during the wet season, the Amazon rises and inundates its low lying floodplain forests, creating an internationally significant wetland ecosystem which reaches its maximum extent around June or July (Mitsch and Gosselink, 2015), with the extensive surface water flow recognised a key factor in the functioning of the habitats created and exerting a strong control on biological and biogeochemical processes (Richey et al., 2002;Wittmann et al., 2004).Inter-annual variability in the flood level and flood wave timings creates spatial and temporal complexities on the floodplain (Melack and Hess, 2011), with the source of water recognised as a key determinant for sediment transport (Vauchel et al., 2017) and trace element concentration (Baronas et al., 2017) dynamics across the basin.
We applied the model of Wilson et al. (2007) to a ∼300 km reach of the Solimões River (mainstem Amazon River) at its confluence with the Purus River (Figure 8), using a model grid of ∼278 m (domain size: 966 x 479 cells, ∼35,650 km 2 ) for the period of October 2013 through December 2014.Operating at a far larger spatial scale to examples 1 and 2, here water source tracing results show that, downstream (west) of the river confluence, the southern floodplain is primarily inundated with water from the River Purus (indicated by green), while the floodplain to the north contains a mix of water from the Rivers Solimões and the Purus (indicated by orange) (Figure 9).The shared floodplain upstream of the confluence displays a clear gradient transition between the two water sources, passing from red (Solimões), through yellow (equally mixed) to green (Purus).The cross-section profile plots in Figure 9 show that, after flowing into the Solimões, water from the Purus remains largely on the right (southern) bank, but mixes with the water from the Solimões as the river progresses downstream.As the Purus peaks sooner than the Solimões, the proportion of water from the Purus close to the right bank decreases through the Solimões flood peak (C) and mid-falling water (D).These demonstrate the ability of the LISFLOOD-FP and the new tracing functions to simulate the convergence and mixing of water from these two major river systems.The proposed water source tracing method  therefore provides the means to assess the complex dynamics of variable contributions in river-floodplain systems such as those on the Amazon, enabling the estimation of important variables which are associated with each water source through time, such as nutrient availability.

Computational performance
Results of the planar test case are shown in Figure 3b-d.Flow depth at t = 8400 s is shown in Figure 3b, and a water source visualisation at the same time for the simulation with all 8 sources is shown in Figure 3c.Here, flood water mixing is shown in RGB colour space, as per ( 14) i.e. the mixing of sources 2 (red) and 4 (green) gives yellow; mixing source 4 and 6 (blue) gives cyan.Other water sources are shown in black.As the water flows down the slope, it becomes increasingly mixed within each section, helped by water being forced to flow through the gaps left in the walls placed across the slope.Once water flow reaches the bottom of the slope, a high proportion of the domain contains water (77,901 out of 80,000 cells).The flow is constrained to flow out of the domain through eight gaps in the downstream boundary wall, marked A-H in Figure 3a.The fractions of water sources 1-8, at locations A-H, throughout the simulation are shown in Figure 3d.As the domain has mirror symmetry along the west to east axis at y = 500 m, the fractions at each of the locations are also symmetrical with respect to water sources.For example, at location A, water from source 1 increases to a fraction of around 0.5 of the water flowing through this location at steady-state, while water from source 8 does not reach it at all; the reverse is true for these water sources at location H.For the two central locations, D and E, all water sources are present with steady-state fractions ranging from around 0.02 (for source 8 at location D or source 1 at location E), to 0.2 (for sources 4 and 5).The total time of each of the eight simulations is presented in Table 1, ranging from 1362.5 s with 2 tracers, up to 1925.2 s with 8 tracers, as compared to 800.6 s without water source tracing.The ratio of computation time for the full simulations with tracing, against computation time without tracing, ranges from 1.7 for two sources to 2.4 for eight.Mean computational requirements time ranged from 0.016 s/s for 2 tracers to 0.022 s/s with 8 tracers, compared to 0.009 s/s without water source tracing.In addition, the variability in computational costs increases with the number of sources being traced, likely due to the complexity of the mixing introduced by constraining flow through gaps within this test case.
Figure 10 illustrates the computational efficiency of our scheme for the Carlisle flood model, with and without water source tracing.Computational requirements for the simulation with water source tracing were found to be between 1.2x and 1.5x the simulation without tracing, which is less than the 1.8 ratio found for the planar test case with 3 water sources.As with the   planar test case, the computational requirements of water source tracing were found to increase with the area of inundation.

285
Early in the simulation, when the flood area was around 1 km 2 (40,000 cells), the computational overhead of water source tracing was approximately 0.01 s/s.This increased nearly linearly to ∼0.05 s/s as the flood extent reached its peak extent of around 5 km 2 (200,000 cells).

Discussion
Our method for water tracing is simple, effective and has a low additional computational overhead of 1.2-1.5 for our real world test cases.The three case studies above show that the approach is flexible, applicable in urban, coastal and fluvial environments over a wide range of spatial scales.An important advantage of the method is that the number or water sources which may be traced is limited only by computational constraints.Furthermore, as the formulation records the proportion of different tracers or water sources per cell, it is fully independent of the hydraulic calculations, meaning it is straightforward to add to existing finite volume codes and similar 2D flow models to LISFLOOD-FP.This also means that predicted water levels with or without the tracers are identical.
However, it is important to remember that the speed and simplicity of the method is at the expense of treating cell water volumes as being fully mixed.This means that it is possible for small fractions of a tracer/water from different sources to propagate quickly downstream, since fluxes into a cell would be included for all the fractions assigned as an inflow to a downstream neighbouring cell.In effect, this represents a tiny amount of numerical diffusion and the volumes moved are very small, but as a result caution should be applied to the interpretation of very small water source fractions.This is also an issue for tracing functions in other codes such as TELEMAC (e.g.Ch. 9.5 Ata et al., 2014).Additionally, as each cell is fully mixed we cannot yet incorporate the effects of different water densities, such as saline water.
These issues notwithstanding, this method provides excellent opportunities for simplified simulation of water, solute and water borne fluxes across river catchments and in estuaries.In particular, with the CAESAR-Lisflood implementation combining a hydrological and hydraulic model, the contributing effects of sub basins on flooding can easily be simulated.Furthermore, the ability to combine both point source (e.g. from combined sewer outflows) and diffuse source (e.g.fields) contaminant sources within a computationally fast model opens up many opportunities to simulate water quality and pathogen/contaminant issues.Flow rates between all cells are then calculated as normal in the function qroute() (i.e.solving Eqn. 2 between all cells), with no modifications needed to the code for water source tracing.Once flows are calculated, depths are updated according to (1) in the function depth_update().Here, the only modification required for water source tracing is that the volumes of water moving between cells is recorded in the dhdt_x and dhdt_y arrays (Figure A4).Here, the dhdt_x and dhdt_y arrays are updated on lines 18-24, before the water_depth array is updated according to (1) on line 27.For brevity, the code presented in Figure A4 omits some additional processing for updates to suspended sediment and error checking.

Figure 1 .
Figure 1.The main CAESAR-Lisflood control showing additions for water source tracing.Additional functions related to erosion and deposition are not shown.

Figure 2 .
Figure 2. Illustration of RGB colour scaling used for visualisation of water source fractions: (a) RGB values for water fractions ϕ from source w, obtained using (14) (upper plot, no depth visualisation) and (15) (lower plot, with depth visualisation, where h [range] = 10), and visual enhancement of lower water fractions using β; (b) RGB values with changing depth, h, for different values of ϕw and β; (c) colour rendering resulting from (14) for three water sources, R, G and B (note that ϕR + ϕG + ϕB = 1.0);(d) colour scales obtained for various combinations of water sources using (15), with and without enhancement of lower water fractions, showing reduced lightness as water depth approaches h [range] .Note that the bottom colour scale is used for rendering in the absence of any water from sources selected for visualisation (i.e.where more than three sources are present in the simulation). 155

Figure 3 .
Figure 3. Planar test case setup and results: (a) the 2000 by 1000 m planar slope (0.001 m/m) elevation model (5 m grid spacing, 400 x 200 cells), with 1 m "walls" added across the slope at 250 m intervals (each wall had between 4 and 8 evenly-spaced gaps of 5 m, as indicated by the numbering along the top of the grid); injection points for each of the 8 water sources are indicated by numbering on the left; (b) water depth prediction at t = 8400 s (140 minutes), with a constant inflow of 10 m 3 /s and n = 0.05; (c) visualisation of water sources numbered 2(red), 4 (green) and 6 (blue) at t = 8400 s, using (14) with β = 0.2; (d) water source fractions throughout the 24-hour simulation for each of 1-8 sources at locations A-H indicated in (a).For an animated version, please see: https://youtu.be/DTw8ysJtx8o.

Figure 4 .
Figure 4. (a) Study site 1 at Carlisle, UK, at the confluence of the Rivers Caldew, Petteril and Eden, showing elevation from LiDAR.The plots in (b) show river flows for the January 2005 flood event simulated.The plot in (c) shows the relative level of the three rivers (calculated as a proportion of the difference between the flow on 8 January and the flood peak), indicates that the River Caldew reached flood peak ∼5.5 hours ahead of the River Petteril, and ∼17 hours ahead of the River Eden.Contains OS data © Crown copyright and database right (2020).The ability to trace inputs from discrete tributaries allows us to determine how different sourced flood waters contributed to the flooding.This is especially pertinent to the Carlisle 2005 flood events, where most flood damage was considered to have come from the two tributary inputs, the Caldew and Petteril(Environment Agency, 2006).For example, the time-series plots in Figure5illustrate locations where water is mixed from multiple sources.The larger volume of water from the River Eden dominates, but at point A, a railway embankment prevents flooding from the Eden but the area is flooded by the Caldew.Points B and D are initially flooded by the Caldew and Petteril, respectively, until flooding from the Eden arrives.This is likely due to the timing of the flood peaks (∼11.5 hours earlier on the Petteril than the Eden; ∼17 hours earlier on the Caldew).

Figure 5 .
Figure 5. Model results for Carlisle, showing flood extent on 10 January 2005, 06:00 AM, close to the flood peak.Map colours indicate mixing of water, with β = 0.2 to emphasise lower water fractions from the Caldew and Petteril Rivers.Depths and fractions for four selected locations, marked A-D, are shown.Contains OS data © Crown copyright and database right (2020).For an animated version, please see: https://youtu.be/xOtOi06cXvA.

Figure 6 .
Figure 6.(a) Study site 2 at the Avon-Heathcote estuary in Christchurch, New Zealand showing elevation from LiDAR; (b) inflow and downstream tide levels for the period simulated in July 2017, which includes a high-flow event on 22 July.The four vertical lines labelled A-D represent the timing of the images shown in Figure 7. Tide data are from the NIWA Sumner Head gauge, around 1 km from the estuary outlet; flow data are from the Environment Canterbury gauges at Gloucester St. (Avon) and Buxton Terrace (Heathcote).Contains data sourced from the LINZ Data Service licensed for reuse under CC BY 4.0 https://doi.org/10.5194/gmd-2021-340Preprint.Discussion started: 18 October 2021 c Author(s) 2021.CC BY 4.0 License.

Figure 7 .
Figure 7. Model results for the Avon-Heathcote estuary: low flow conditions during (A) low-tide on 8 July 2017 and (B) after spring-high tide on 16 July; high flow conditions on 22 July shortly after high tide (C) and low tide (D).The images were produced with β = 0.2; their timing and flow/ tide levels are shown in Figure 6.The background images contain the NationalMap Basemap (CC-BY-ND license) and LiDAR data shaded in grey.For an animated version, please see: https://youtu.be/Fczr5tczzXU.

Figure 8 .
Figure 8.(a) Study site 3 at the confluence of the Solimões (mainstem Amazon) and Purus rivers in the central Amazon, Brazil, showing elevation derived from SRTM (O'Loughlin et al., 2016) and bathymetry (Wilson et al., 2007).The plots in (b) and (c) show river flows for the period simulated from October 2013 through December 2014.The lower plot of (b), which shows the relative level of Solimões and Purus flow (calculated as a proportion of the difference between the flow on 1 November 2013 and the 2014 peak flow), indicates that the Solimões rises more steadily and before the Purus, which has a later, steep rise from around mid-February onwards.Peak flow levels on the Purus (indicated by letter B) occurred around one month earlier than peak flow on the Solimões.The dotted box in (a) and letters A-D in (b) and (c) indicate the location and timing of the model results shown in Figure 9. Contains Natural Earth data in the public domain.

Figure 9 .
Figure 9. Model results for the Amazon, with colours on the map at the Purus flood peak (top) indicating the mixing of the waters along the channel and across the floodplain.Channel cross-section plots (bottom) at locations indicated by a-a ′ and b-b ′ show the water source fractions (upper plots) for each channel (A = mid-rising; B = Purus flood peak; C = Solimões flood peak; D = mid-falling), and depth/ source profiles (lower plots) for image B (Purus flood peak).For an animated version, please see: https://youtu.be/PknAL_8fd1I.

Figure 10 .
Figure 10.Computational efficiency of water source tracing for the Carlisle simulations with and without water source tracing: (a) computation time required per second of the 5-day (432,000 s) simulation, with colours indicating the area of inundation; the dashed line shows the ratio between the two (secondary y-axis); (b) computation time difference (i.e.cost of tracing) against flood area, with a linear model fitted for reference; colours indicate the simulation time.
https://doi.org/10.5194/gmd-2021-340Preprint.Discussion started: 18 October 2021 c Author(s) 2021.CC BY 4.0 License.We developed a simple method for tracing water sources through a flood inundation model and demonstrated its application for a simple test case and three example case studies at different spatial and temporal scales.A key advantage of the approach lies in its independence from the hydraulic methodology used, meaning that it is relatively straightforward to add to existing in finite volume codes.The method enables effective and informative visualisation of flood inundation, providing additional insight into flood dynamics.In addition, it may potentially provide the ability to simulate the movement of solute contaminants and pathogens within a computationally efficient model, although additional testing is required to verify the reliability of the method, given the assumption of full mixing within cells.https://doi.org/10.5194/gmd-2021-340Preprint.Discussion started: 18 October 2021 c Author(s) 2021.CC BY 4.0 License.The addition to the reach inputs is shown in FigureA2, from the function reach_water_and_sediment_input() (rainfall inputs and tidal or stage inputs are added similarly but code is omitted for brevity).In this function, the updated water source fractions for the water source being added to the cell (Eqn.9) and other sources already in the cell (Eqn.10) are implemented on lines 33 and 40, respectively.

Figure A5 .
Figure A5.The main water source tracing calculations within the update_tracer_states() function