Articles | Volume 16, issue 9
Development and technical paper
08 May 2023
Development and technical paper |  | 08 May 2023

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

Matthew D. Wilson and Thomas J. Coulthard

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 the RGB colour space, while continuing to allow depth to be resolved. The number of water sources that may be traced is limited only by the 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 method to address environmental issues relating to water sources and mixing, such as water quality and contamination problems.

1 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 (2D) hydrodynamic models have been developed (Hunter et al.2007; Neal et al.2018). These model codes have enabled the assessment of flood inundation at high spatial resolutions (e.g. Yu and Coulthard2015) 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.2022; Wing et al.2017), and global scales (Dottori et al.2016; Sampson et al.2015). The LISFLOOD-FP model (Bates and De Roo2000; 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, in addition to 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 by, for example, determining the relative contribution of tributaries or where waterborne contamination is an issue. This 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 modelling in a lake (Kopmann and Markofsky2000) and solute and viral dispersal in estuaries (Robins et al.2014, 2019). For example, in TELEMAC, non-buoyant tracers can be added and their course and concentration followed (Sect. 9; Ata et al.2014). Similarly, in PCSWMM, Qi et al. (2021) and Qi et al. (2022) developed a module to assess the relative contributions to downstream flood waters from upstream source catchments, with tracer sources generated by the PCSWMM water quality routing module. In their approach, each upstream source is assigned a constant tracer concentration, which is then routed downstream using PCSWMM methods for the transportation of pollutants. The source proportion for a downstream catchment is then determined from the relative mass of the total tracer amount which is in that catchment and multiplied by the total flood volume. However, approaches using such 2D 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 without the addition of tracers 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 (available from, last access: 25 April 2023) 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.

In earlier work with CAESAR, Coulthard and Macklin (2003) demonstrated how sediment eroded from mining waste deposits could be traced downstream. The model worked by including different types of sediments which were used to represent contaminated and uncontaminated sediments as separate arrays for each sediment diameter included; during the erosion and deposition of different sediment sizes, equal proportions of contaminated and uncontaminated sediment were transported. This enabled the prediction of patterns and levels of floodplain contamination over a period of ∼400 years. However, the approach was limited to sediment tracing only and did not account for different sources of water, which may carry contaminants with it. As noted by Coulthard et al. (2013), the inclusion of an unsteady flow formula within CAESAR-Lisflood has enabled the possibility, for example, to simulate water balances and solute fluxes using the model code. Here, we present the formulation of a simple methodology to enable this functionality by accounting for the source of water within model cells throughout the simulation.

We use only the hydraulic and hydrological functionality of CAESAR-Lisflood to demonstrate the tracing method and visualisation of water sources independently of sediment routing. We have used CAESAR-Lisflood for this purpose, as the software is fully open source, integrating a graphical user interface (GUI) in addition to the hydraulic methods from LISFLOOD-FP. This enables our visualisation methods to be incorporated but does not limit our tracing method to CAESAR-Lisflood. The equations and pseudo-code examples provided make it a straightforward task for users and researchers to add this functionality to LISFLOOD-FP (and other variants using similar methods).

2 Methods

2.1 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 and for cell i,j, is determined using the following (Bates and De Roo2000):

(1) h i , j t = h i , j t - Δ t + Δ t Q x , i - 1 , j t - Q x , i , j t + Q y , i , j - 1 t - Q y , i , j t Δ x 2 ,

where hi,jt-Δt is the cell depth at the end of the previous time step (Δ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 the following (de Almeida et al.2012; Bates et al.2010):

(2) Q x t = q x t - Δ t - g h [ flow ] t Δ t Δ ( h t + z ) Δ x ( 1 + g h [ flow ] t Δ t n 2 | q x t - Δ t | / ( h [ flow ] t ) 10 / 3 ) Δ y ,

with flow in the y direction obtained analogously. In Eq. (2), qxt-Δt represents the flux between the cells from the previous time step (Qxt-Δt/Δ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 Bates2013; Neal et al.2012).

2.2 Water source tracing

Given that the sum of depths in a cell from each source, w, make up its total flood depth, h,

(3) w = 1 W h w = h ,

the volume fraction of each source, ϕw, may be obtained from

(4) ϕ w = h w h ,


(5) w = 1 W ϕ w = 1 .

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, our water tracing method 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 time step. (3) Contributing inflow depths from each source for all neighbours is added to each source depth. (4) Updated water source fractions are calculated, as per Eq. (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 as follows:

(6) ϕ w t = h Q [ out ] t ϕ w t - Δ t + D ( Q [ in ] t , D ϕ w t - Δ t , D ) Δ t Δ x 2 h t ,

where hQ[out]t represents the depth remaining in the cell after the removal of outflows at the current time, Q[out]t (and before the addition of any inflows), which is scaled according to the fraction from this source from the previous time step, ϕwt-Δt. To this depth the total amount of water from this source is added, which is contributed from neighbouring cells and given by D(Q[in]t,Dϕwt-Δt,D), where Q[in]t,D is the inflow from a particular direction, D, and ϕwt-Δt,D 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, ht. Finally, hQ[out]t is given by the following:

(7) h Q [ out ] t = h t - Δ t - D Q [ out ] t , D Δ t Δ x 2 ,

where Q[out]t,D are the outflows in each direction. To improve the computational efficiency, the method only requires knowledge of the source volume fractions for the previous time step in each neighbouring, upstream cell; the actual depths or volumes from each source in each cell are not saved. The calculation cell depth after the removal of outflows, hQ[out]t, enables each source fraction to be updated based on the flows in (Q[in]t,D), the fraction from each source in direction, D (ϕwt-Δt,D), and the updated depth, ht. Using this approach, the calculation of water source fractions flowing out of the cell is not necessary. For completeness, using the notation of Eq. (1) and expanding D into cell indices, i,j, the updated depth from each source can be obtained using the following:

(8) ϕ w , i , j t h i , j t = ϕ w , i , j t - Δ t h i , j t - Δ t + Δ t ϕ w , i - 1 , j t - Δ t Q x , i - 1 , j t - ϕ w , i + 1 , j t - Δ t Q x , i , j t + ϕ w , i , j - 1 t - Δ t Q y , i , j - 1 t - ϕ w , i , j + 1 t - Δ t Q y , i , j t Δ x 2 .

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

(9) ϕ w = h ϕ w + Q w Δ t Δ x 2 h + w Q w Δ t Δ x 2

for input sources or

(10) ϕ w = h ϕ w h + w Q w Δ t Δ x 2

for other sources, where Qw is the inflow added to the cell from source w at the start of the time step. Thus, fractions from sources where water is added to the cell are adjusted upwards, while fractions for non-source volumes are adjusted downwards. This means that additional inflow from a certain source results in an increase in the fraction of that source in that cell, which is reflected in a higher value of the water source fraction variable, ϕw; the opposite occurs for fractions of water from other sources, as they are diluted.

It should be noted that, similar to the approach of Qi et al. (2021), for simplicity, this method treats cell water volumes as being 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 time step via ϕwt-Δt,D. As a result, caution should be given to the interpretation of very small water source fractions.

2.3 Water tracing implementation

Importantly, within the model code, the additional equations required for the water source tracing presented in Sect. 2.2 can be piggybacked on top of the existing finite-volume model code, requiring minimal modification of the original schemes. The water source tracing method was implemented in C#, building from the code base in CAESAR-Lisflood version 1.9j (Coulthard2019). 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 of water as point sources or via rainfall through a hydrological model based on a spatially distributed version of TOPMODEL (Beven and Kirkby1979; for implementation, see Coulthard and Skinner2016, and Coulthard and Van De Wiel2017). The zonal input also allows tidal sources to be traced, although it is worth noting 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 grid cells sized W×X×Y for the current and previous time step, 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 time step.

The main functional flow of CAESAR-Lisflood, as related to water source tracing, is shown in Fig. 1. For each time step, inflows from various sources are first added. The addition of traced water requires the 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 (Eqs. 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 time step, ht−1 and ϕt−1. For the water source tracing implementation, the tracer states are saved via the function save_tracer_states() (Fig. 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 the volumes of water moving between cells are stored to the δhδt in x and y arrays, which are calculated using

(11) δ h δ t t , D = Q t , D Δ t Δ x

for a single direction (C# code in Fig. A4). This is used as part of the solution to Eqs. (6) and (7) to update the water source tracers, enabling a simplified form to be implemented for the tracer update as follows:

(12) ϕ w t = h Q [ out ] t ϕ w t - Δ t + D ( δ h δ t [ in ] t , D ϕ w t - Δ t , D ) h t ,


(13) h Q [ out ] t = h t - Δ t - D δ h δ t [ out ] t , D .

Water source tracers are updated in the function update_tracer_states(), which is the main addition to the CAESAR-Lisflood code. Algorithm 2 provides the details of the procedure used. For every cell in the model which contains water (hi,jt>0.0), the depth remaining in the cell after the outflows (hQ[out]t) is first calculated by subtracting the outflows from each of the four flow directions. We sum δhδt from each direction (previously calculated in depth_update()), which represents the contributions to depth from the right (x) or up (y). Thus, outflows will be negative for the right and upward directions and positive in the left and downward directions.

The procedure then updates each source for the cell by first calculating the sum of depths added from this source (i.e. D(δhδt[in]t,Dϕwt-Δt,D)), which is again done by scanning the δhδt arrays in each direction. Finally, the fraction is calculated by dividing the total new depth from the source by the depth from all sources, ht.

Algorithm 1Updating water source fractions after adding the input sources. Inputs required are the depth in the cell for the previous iteration (ht−1), a vector of all water sources for the cell for the previous iteration (ϕt−1), input from this source (Qw), the index of this source (w), the number of sources (W), the cell size (Δx), and the time step (Δt). A C# implementation is shown in Fig. A2.

Input: ht−1, ϕt−1, Qw, w, W, Δx, Δt
Output: ϕt
δhδtwQwΔt/Δx2 {get change in depth resulting from this source}
htht-1+δhδtw {get new depth after adding input from this source}
if ht>0.0 then
if ht=δhδtw then
ϕwt1.0 {if depth was zero, then all water in cell is from this source}
else if ϕwt-1<1.0 {if cell contains multiple sources} then
hwtϕwt-1ht-1 {get depth after addition of this source}
ϕwt(hwt+δhδtw)/ht {update this water source fraction}
for src=1 to W do
if src≠w then
ϕsrct(ϕsrct-1ht-1)/ht {update fraction of other water sources}
end if
end for
end if
end if

Algorithm 2Updating water source fractions after the flow routing and depth update. Inputs required are the updated depth in the cell (ht), water sources for the previous iteration (ϕt−1), depth contributions horizontally (δhδtx) and vertically (δhδty), the number of sources (W), and the number of cells in the grid horizontally (X) and vertically (Y). A C# implementation is shown in Fig. A5.

Input: ht, ϕt−1, δhδtx, δhδty, W, X, Y
Output: ϕt
for all i, j cells in ranges 1iX and 1jY do
if hi,jt>0.0 then
δhδtsumQ[out]0.0 {obtain total change in depth after flows out (source not important)}
if δhδtxi+1,j<0.0, then δhδtsumQ[out]δhδtsumQ[out]+δhδtxi+1,j; end if {right, ve = outflow}
if δhδtxi,j>0.0, then δhδtsumQ[out]δhδtsumQ[out]-δhδtxi,j;     end if {left, +ve = outflow}
if δhδtyi,j+1<0.0, then δhδtsumQ[out]δhδtsumQ[out]+δhδtyi,j+1; end if {up, ve = outflow}
if δhδtyi,j>0.0, then δhδtsumQ[out]δhδtsumQ[out]-δhδtyi,j;     end if {down, +ve = outflow}
hQ[out]tht-1+δhδtsumQ[out] {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)}
if δhδtxi+1,j>0.0, then δhδtsumQ[in]δhδtsumQ[in]+δhδtxi+1,jϕi+1,j,srct-1; end if {right, +ve = inflow}
if δhδtxi,j<0.0, then δhδtsumQ[in]δhδtsumQ[in]-δhδtxi,jϕi-1,j,srct-1; end if {left, ve = inflow}
if δhδtyi,j+1>0.0, then δhδtsumQ[in]δhδtsumQ[in]+δhδtyi,j+1ϕi,j+1,srct-1; end if {up, +ve = inflow}
if δhδtyi,j<0.0, then δhδtsumQ[in]δhδtsumQ[in]-δhδtyi,jϕi,j-1,srct-1; end if {down, ve = inflow}
if δhδtsumQ[in]=0 and ϕsrct-1=0 then
ϕi,j,srct0.0 {if there is no contribution or existing water from this source}
ϕi,j,srct(hQ[out]tϕi,j,srct-1)+δhδtsumQ[in]/hi,jt {([src depth in cell] + [src depth added]) / [updated depth]}
end if
end for
end if
end for

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


2.4 Water source visualisation

Visualisations of the 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 the 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 the following:

(14) RGB i , j = 255 ϕ i , j , w β .

The power term, β, enables a visual enhancement of the 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, Eq. (14) may be modified to the following:

(15) RGB i , j = ( ( 1 - h i , j / h [ range ] ) 127 ) + 128 ϕ i , j , w β ,

where hi,j is the cell depth, and h[range] is the range of depths to which the visualisation is scaled. Thus, the RGB value for a cell which only has one source of water in it will range from 128 for deeper water, where h=h[range], towards 255 for shallow depths close to 0. Figure 2 provides a detailed illustration of the RGB values obtained using Eqs. (14) and (15) for different values of ϕ, β, and h and the resulting colour scales obtained when visualising multiple sources. The mixing of water sources with red and blue shading gives a pink colour scale, green and blue give cyan, red and green give yellow, and shades towards white would indicate the mixing of all three sources. Note that, while it is possible to trace an arbitrary number of water sources using the method described in Sect. 2.2, this visualisation method is limited to a maximum of three water sources. However, using Eq. (15), if a cell only contains water from sources different to those selected for visualisation, then the depths can still be resolved in greyscale in the RGB value range from 0 to 127. This also means that it is possible to visualise only one water source at a time in a selected primary colour.

Figure 2Illustration of the RGB colour scaling used for the visualisation of water source fractions. (a) RGB values for water fractions ϕ from source w, obtained using Eq. (14) (upper plot; no depth visualisation) and Eq. (15) (lower plot, with depth visualisation, where h[range]=10) and a 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 Eq. (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 Eq. (15), with and without enhancement of lower water fractions, showing reduced lightness as the water depth approaches h[range]. Note that the bottom colour scale is used for rendering in the absence of any water from the sources selected for visualisation (i.e. where more than three sources are present in the simulation).


2.5 Evaluating performance overhead

In order to test the computational performance of the water tracing code, a simple planar test case was developed, consisting of a constant 0.001 m m−1 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 (Fig. 3a). Manning's roughness, n, was set at 0.05. At 250 m intervals downslope, 1 m “walls” were added across the slope, each with between four and eight gaps of 5 m width added to allow water to flow through. This was done to ensure that water mixed as it flowed downslope. Simulations were conducted for a period of 24 h, with constant inflow of 10 m3 s−1 to eight locations at the top of the slope. Eight simulations were conducted in total, including one with no tracing and seven with between two and eight water sources being traced. To ensure that each of these simulations had identical total inflow, the 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 the model performance in a real-world scenario, the effect of the number of tracers on model performance was tested on the Carlisle example described below.

Figure 3Planar test case set-up and results. (a) The 2000×1000m planar slope (0.001 m m−1) elevation model (5 m grid spacing; 400×200 cells), with 1 m “walls” added across the slope at 250 m intervals (each wall had between four and eight evenly spaced gaps of 5 m, as indicated by the numbering along the top of the grid). Injection points for each of the eight water sources are indicated by the numbering on the left. (b) Water depth prediction at t=8400s (140 min), with a constant inflow of 10 m3 s−1 and n=0.05. (c) Visualisation of water sources numbered 2 (red), 4 (green), and 6 (blue) at t=8400s, using Eq. (14), with β=0.2. (d) Water source fractions throughout the 24 h simulation for each of 1–8 sources at locations A–H, as indicated in panel (a). For an animated version, please see (last access: 25 April 2023) or the video supplement provided.


Table 1Planar slope test case simulations. In each of the eight simulations, the total volume of inflow was identical, and each of the eight sources was a constant 10 m3 s−1 (where a source was used for multiple inflow points, its total volume was added to all). Inflow points are shown in Fig. 3a. Note that SD is the standard deviation.

Download Print Version | Download XLSX

3 Results

First, 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, namely (1) a major flood event in 2005 at Carlisle, United Kingdom, which entailed fluvial flooding in an urban area that is situated at the confluence of three rivers, (2) a shallow estuary in Christchurch, New Zealand, combining tides with inflow from two small rivers, and (3) flooding at the tributary junction of two large rivers on the Amazon floodplain, Brazil. Second, we present the computational performance of the method for the planar test case and Carlisle example.

3.1 Carlisle, United Kingdom

The city of Carlisle in northern England experienced significant flooding in 2005 (estimated annual exceedance probability from 0.5 % to 0.57 %) that resulted from heavy rains in the headwater catchments of rivers running through the city and affecting more than 2500 properties (Environment Agency2006). With the benefit of an extensive set of observational data obtained from field collection, the 2005 flood event at Carlisle 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 (Fig. 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 panel (b) show river flows for the January 2005 flood event simulated. The plot in panel (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) and indicates that the River Caldew reached the flood peak ∼5.5 h ahead of the River Petteril and ∼17 h ahead of the River Eden. This figure contains OS data © Crown copyright and database right (2020).

A simulation was run using a model grid of 5 m (domain size of 951×612 cells; 14.6 km2), topography from lidar (Environment Agency2020), and inflows from gauging station records (Environment Agency2021). The simulation began on 8 January 2005, 00:00 UTC for 120 h (5 d). Visualisations of the model output are shown in Fig. 5. Flood water mixing is shown in the RGB colour space, using Eq. (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.

Figure 5Model results for Carlisle, showing the flood extent on 10 January 2005, 06:00, close to the flood peak. Map colours indicate the mixing of water from different sources, with darker shades indicating increased depth. Colour scales show depths with selected source combinations (left to right), namely ϕ [Eden]=1.0 (blues), ϕ [Eden]=0.8, ϕ [Petteril]=0.2 (pinks/purples), ϕ [Petteril]=1.0 (reds), and ϕ [Caldew]=1.0 (greens). β = 0.2 is used to emphasise lower water fractions from the Caldew and Petteril rivers. Depths and fractions for four selected locations, marked A–D, are shown. This figure contains OS data © Crown copyright and database right (2020). For an animated version, please see (last access: 25 April 2023) or the video supplement provided.

The ability to trace inputs from discrete tributaries allows us to determine how differently sourced flood waters contributed to the flooding. This is especially pertinent to the Carlisle 2005 flood events, where the most flood damage was considered to have come from the two tributary inputs, namely the Caldew and Petteril (Environment Agency2006). For example, the time series plots in Fig. 5 illustrate the 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 h earlier on the Petteril than the Eden; ∼17 h earlier on the Caldew).

3.2 Avon–Heathcote estuary, Christchurch, New Zealand

On the eastern edge of the city of Christchurch, New Zealand, the estuary of the Avon and Heathcote rivers (Fig. 6) is a large (∼8km2), shallow (mean depth ∼1.4 m at spring high tide and a tidal range of ∼1.7–2.2 m) intertidal area, with about 85 % of its elevation above the spring low tide level (Findlay and Kirk1988). The estuary is semi-enclosed and separated from the Pacific Ocean by a ∼4.5km long spit (the Christchurch suburbs of New Brighton and/or Southshore. It has been designated as a site of ecological significance in the Christchurch District Plan (Christchurch City Council2016) and is internationally recognised for its importance as a habitat for migratory bird species (Woodley2018). 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 km2, 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 of 483×675 cells; 32.6 km2), derived from lidar (Land Information New Zealand2015). 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 southwest (Environment Canterbury, 2023) and a downstream tide boundary condition at the estuary outlet at the southeast (Fig. 6b). The simulated period of July 2017 included both the 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 m3 s−1 of flow, close to the mean annual flood level (Environment Canterbury2023), and Heathcote River peaked at 14:00 with 28.1 m3 s−1 of flow, which is an annual exceedance probability of less than 0.1 (Environment Canterbury2023).

As noted in Sect. 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 variations of any kind in the vertical salinity profile of estuarine water are 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 of 1.87 and 1.04 m3 s−1, respectively (LAWA2023), compared with tidal inflow rates of up to 500–800 m3 s−1 during an incoming tide (van der Peet and Measures2015), it is likely that it is well mixed, with little vertical salinity difference (Hansen and Rattray1966).

The results of the model application are shown in Fig. 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 the 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 storm water 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. McMurtrie2015).

During the high-flow conditions on 22 July, shortly after the high tide and inflow peak (image C) and low tide (image D), the results indicate a substantially greater riverine contribution to the estuary, such that, after the fall of the tide, the estuary remained substantially inundated but with almost all water being from riverine sources. After high tide, the Avon and Heathcote waters are largely confined by the tide to the western shore, with a clear boundary between the two river sources being 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.

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 Fig. 7. Tide data are from the National Institute for Water and Atmospheric Research (NIWA) Sumner Head gauge, which is around 1 km from the estuary outlet. Flow data are from the Environment Canterbury gauges at Gloucester Street (Avon) and Buxton Terrace (Heathcote). This figure contains data sourced from the Land Information New Zealand (LINZ) Data Service and is licensed for reuse under CC BY 4.0.

Figure 7Model results for the Avon–Heathcote estuary. Low-flow conditions are shown during (a) low tide on 8 July 2017 and (b) after the spring high tide on 16 July. High-flow conditions are shown on 22 July shortly after high tide (c) and low tide (d). The images were produced with β=0.2. Their timing and flow and/or tide levels are shown in Fig. 6. Map colours indicate the mixing of water from different sources, with darker shades reflecting increased depth. The colour scales show depths for individual sources (left to right), namely ϕ [Tide]=1.0 (blues), ϕ [Avon]=1.0 (reds), and ϕ [Heathcote]=1.0 (greens). The background images contain NationalMap Basemaps (CC-BY-ND license) and lidar data shaded in grey. For an animated version, please see (last access: 25 April 2023) or the video supplement provided.

3.3 Amazon River, Brazil

The flow of the central Amazon River is characterised 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 in around June or July (Mitsch and Gosselink2015), with the extensive surface water flow recognised as being 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). The inter-annual variability in the flood level and flood wave timings creates spatial and temporal complexities on the floodplain (Melack and Hess2011), with the source of the water being 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 (Amazon River mainstem) at its confluence with the Purus River (Fig. 8), using a model grid of ∼278m (domain size of 966×479 cells;  35 650 km2) for the period of October 2013 through December 2014. Operating at a far larger spatial scale compared 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 Purus (indicated by orange; Fig. 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 (equal mix), and to green (Purus). The cross-sectional profile plots in Fig. 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. Furthermore, it may enable a rapid assessment of issues of importance for human health, such as identifying areas which may be susceptible to mercury contamination from gold mining activities in the Amazon basin, for which the dynamics of fluvial and pluvial waters are key factors (Maurice-Bourgoin et al.2003; Maia et al.2009). However, for a robust analysis, it is likely that additional tracing of sediment erosion, transportation, and deposition is required (e.g. Haddadchi et al.2019).

Figure 8(a) Study site 3 at the confluence of the Solimões (Amazon River mainstem ) and Purus rivers in the central Amazon, Brazil, showing elevation derived from SRTM (Shuttle Radar Topography Mission; O'Loughlin et al.2016) and bathymetry (Wilson et al.2007). The plots in panels (b) and (c) show river flows for the period simulated from October 2013 through December 2014 (Agência Nacional de Águas e Saneamento Básico2023). The lower plot of panel (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 Purus, which has a later, steep rise from around mid-February onwards. Peak flow levels on Purus (indicated by letter B) occurred around 1 month earlier than the peak flow on Solimões. The dotted box in panel (a) and letters A–D in panels (b) and (c) indicate the location and timing of the model results shown in Fig. 9. This figure contains Natural Earth data in the public domain.

Figure 9Model results for the Amazon, with colours on the map at Purus flood peak (top) indicating the mixing of the waters along the channel and across the floodplain. Map colours indicate the mixing of water from each source, with darker shades reflecting increased depth. Colour scales show depths with selected source combinations (left to right), namely ϕ [Solimões]=1.0 (reds), ϕ [Purus]=1.0 (greens), and ϕ [Solimões]=0.75, ϕ [Purus]=0.25 (yellows). Channel cross-sectional 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 (last access: 25 April 2023) or the video supplement provided.


3.4 Computational performance

The results of the planar test case are shown in Fig. 3b–d. Flow depth at t=8400s is shown in Fig. 3b, and a water source visualisation at the same time for the simulation with all eight sources is shown in Fig. 3c. Here, flood water mixing is shown in the RGB colour space, as per Eq. (14) (i.e. in which the mixing of sources 2 (red) and 4 (green) gives yellow and the mixing of sources 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 and is helped by water being forced to flow through the gaps left in the walls placed across the slope. Once the 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 Fig. 3a). The fractions of water sources 1–8, at locations A–H, throughout the simulation are shown in Fig. 3d. As the domain has mirror symmetry along the west–east axis at y=500m, 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).

Code profiling, using the diagnostic tools of the Microsoft® Visual Studio 2019 compiler software (version 16.11.2), based on this test application with three sources, indicates that update_tracer_states() uses around 17.9 % of the CPU resources allocated to CAESAR-Lisflood during a simulation, with 3.0 % used for save_tracer_states(). This compares to 39.6 % for qroute() and 11.5 % for depth_update(). This profiling was completed on a Microsoft® Windows® 10 (build 18363) computer system with an Intel® Xeon® E-2278G CPU running at 3.4 GHz base (5.0 GHz maximum turbo) frequency with 8 cores or 16 threads, with 128 Gb of system memory (4×32Gb Samsung M391A4G43MB1-CTD ECC RAM UDIMM running at 2666 Mbps), writing outputs to a Samsung NVMe™ Solid State Disk (SSD; MZVLB1T0). All simulations presented in this paper were completed on this computer system.

Table 2Summary of the simulations completed for each of the three study sites presented and their computational costs.

* Number of computation hours needed per day of the period simulated.

Download Print Version | Download XLSX

Figure 10Computational efficiency of water source tracing for the Carlisle simulations with and without water source tracing. (a) Computation time required per second of the 5 d (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.


The total time of each of the eight simulations is presented in Table 1, ranging from 1362.5 s with two tracers up to 1925.2 s with eight tracers, as compared to 800.6 s without water source tracing. The ratio of the 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. The mean computational requirement time ranged from 0.016 s s−1 for two tracers to 0.022 s s−1 with eight tracers, compared to 0.009 s s−1 without water source tracing. In addition, the variability in the 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.

Table 2 summarises each of the domains and simulations completed for the Carlisle, Avon–Heathcote, and Amazon study sites. Although the computational cost for the Amazon simulation was highest (∼32 d of computation time), when normalising for the larger number of flooded cells and the longer simulation period, it had the lowest computational cost. This higher computational efficiency likely results from the longer model time step used, due to the larger cell size (Δx=277.6 m). For the Avon–Heathcote domain (Δx=10 m), the computational cost per flooded cell per day of period simulated was 2.4 times the cost for Amazon; for the Carlisle domain (Δx=5 m), this ratio increased to 4.2.

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.2× and 1.5× the simulation without tracing, which is less than the 1.8 ratio found for the planar test case with three water sources. As with the planar test case, the computational requirements of water source tracing were found to increase with the area of inundation. Early in the simulation, when the flood area was around 1 km2 (40 000 cells), the computational overhead of water source tracing was approximately 0.01 s s−1. This increased nearly linearly to ∼0.05s s−1, as the flood extent reached its peak extent of around 5 km2 (200 000 cells).

4 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 and applicable in urban, coastal, and fluvial environments over a wide range of spatial scales. An important advantage of the method is that the number of 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 the 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 or water from different sources to propagate quickly downstream, since fluxes into a cell would be included for all of the fractions assigned as an inflow to a downstream neighbouring cell. In effect, this represents a tiny amount of the 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. chap. 9.5 Ata et al.2014). Additionally, as each cell is treated as 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 the simplified simulation of water, solute, and waterborne fluxes across river catchments and in estuaries. In particular, with the CAESAR-Lisflood implementation combining a hydrological and hydraulic model, the contributing effects of subbasins 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) contaminants within a computationally fast model opens up many opportunities to simulate water quality and pathogen and/or contaminant issues. Assessments of water sources may be especially useful for other fluvial or estuarine sites with similar human health considerations to the examples presented here. For example, Robins et al. (2019) assessed water quality risks from viral dispersal from the Conwy catchment, north Wales, and highlighted the importance of river flow contributions to exposure risk. Therefore, this simple and efficient water tracing algorithm within the open-source CAESAR-Lisflood model provides a powerful tool for studying the effect of water quality and contaminant issues on environmental and human health.

5 Conclusions

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 an 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.

Appendix A: Code description

The tracing component requires a series of additional variables representing the water state prior to the routing (Fig. A1). The variable water_depth_prev is an array of X×Y grid cells containing the depth from the previous time step. dhdt_x and dhdt_y are arrays containing the change in cell depth in the last time step in for columns and rows, respectively. The fraction of each cell from each source (i.e. the water tracers) are recorded in a series of three-dimensional (W×X×Y grid cells) or stacked arrays containing the tracer fractions per cell for present and previous iterations (watertracer, watertracer_prev) the rain zonations (watertracerRainZone, watertracerRainZone_prev), and a series of flags to enable or disable tracing. Note that the rainfall source as a whole (i.e. from any zone) is contained within the main watertracer arrays, while the fraction of rainfall from each zone is recorded in the watertracerRainZone arrays.

The addition to the reach inputs is shown in Fig. A2, from the function reach_water_and_sediment_input() (rainfall inputs and tidal or stage inputs are added similarly, but the code is omitted for brevity). In this function, the updated water source fractions for the water source being added to the cell (Eq. 9) and other sources already in the cell (Eq. 10) are implemented on lines 33 and 40, respectively.

After the addition of boundary condition inputs from each source, the new tracer states (prior to surface water routing) are saved within the function save_tracer_states(), in which the waterdepth, watertracer, and watertracerRainZone arrays are copied into the waterdepth_prev, watertracer_prev, and watertracerRainZone_prev arrays, respectively (Fig. A3).

Flow rates between all cells are then calculated as normal in the function qroute() (i.e. solving Eq. 2 between all cells), with no modifications needed to the code for water source tracing. Once flows are calculated, depths are updated according to Eq. (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 (Fig. A4). Here, the dhdt_x and dhdt_y arrays are updated on lines 18–24, before the water_depth array is updated according to Eq. (1) on line 27. For brevity, the code presented in Fig. A4 omits some additional processing for updates to suspended sediment and error checking.

Finally, the main addition to the CAESAR Lisflood code for water source tracing is within the function update_tracer_states() that deals with the mixing component and updates the different source tracer fractions as water is moved between cells during each iteration (Fig. A5). In the update_tracer_states() code, the cell depth after the removal of outflows, hQ[out]t (Eq. 7), is first calculated (lines 31–43). Then, the updated volume fractions (Eq. 6) obtained for each source are calculated (lines 45–81), with the total volume contributed from each source first calculated (lines 47–58), and then the updated cell volume fractions are obtained (lines 77–79). Finally, the function will update the volume fractions for each rainfall zone source in a similar way (code omitted for brevity).

Figure A1Additional variables required for water source tracing.


Figure A2Updating cell water fractions according to Eqs. (9) and (10) from the function reach_water_and_sediment_input().


Figure A3The function save_tracer_states().


Figure A4Modified depth_update() function.


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


Code and data availability

The CAESAR-Lisflood v1.9j-WS source code is freely available for download from Zenodo under the GNU General Public License v3.0, along with the planar test case and a 15 m version of the Carlisle model (; Wilson and Coulthard2023).

Video supplement

Animations for each of the case studies are accessible by following the YouTube links provided in the paper and from Zenodo (; Wilson and Coulthard2021).

Author contributions

MDW designed the methods, developed the code, and performed the model simulations. TJC verified and validated the methods and code. MDW and TJC jointly prepared the paper.

Competing interests

The contact author has declared that neither of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


An earlier version of this work was presented at the GeoComputation 2019 conference in Queenstown, New Zealand (Wilson and Coulthard2019). Lidar elevation data (Environment Agency2020) and river flow data (Environment Agency2021) for the Carlisle site were provided by the Department for Environment, Food and Rural Affairs (Defra) and the Environment Agency under the Open Government Licence. For the Avon–Heathcote site, lidar data (Land Information New Zealand2015) and mapping data were provided by Land Information New Zealand (LINZ; Toitū Te Whenua), river flow data were provided by Environment Canterbury (2023), and tide data were provided by the National Institute for Water and Atmospheric Research (NIWA). For the Amazon site, Shuttle Radar Topography Mission (SRTM) data were used, corrected for vegetation by O'Loughlin et al. (2016), with additional bathymetry data from Wilson et al. (2007); river flow and stage data were provided by the Agência Nacional de Águas e Saneamento Básico (ANA), Brazil. We gratefully acknowledge Rob Thomas for his comments on an earlier version of the paper. We would also like to thank the two referees for their time and effort; we appreciate their constructive comments, which were valuable in helping us to improve the quality of the paper.

Financial support

This research has been supported by Research Councils UK (grant nos. NE/R009007/1, NE/V004247/1, and NE/K008668/1).

Review statement

This paper was edited by Min-Hui Lo and reviewed by Saswata Nandi and one anonymous referee.


Adams, J. M., Gasparini, N. M., Hobley, D. E. J., Tucker, G. E., Hutton, E. W. H., Nudurupati, S. S., and Istanbulluoglu, E.: The Landlab v1.0 OverlandFlow component: a Python tool for computing shallow-water flow across watersheds, Geosci. Model Dev., 10, 1645–1663,, 2017. a

Agência Nacional de Águas e Saneamento Básico: Hidroweb v3.2.7, Agência Nacional de Águas e Saneamento Básico [data set], (last access: 5 May 2023), 2023. a

Ata, R., Goeury, C., and Hervouet, J.: TELEMAC-2D Software Release 7.0 User Manual, Tech. rep., EDF-R&D, (last access: 5 May 2023), 2014. a, b

Baronas, J. J., Torres, M. A., Clark, K. E., and West, A. J.: Mixing as a driver of temporal variations in river hydrochemistry: 2. Major and trace element concentration dynamics in the Andes‐Amazon transition, Water Resour. Res., 53, 3120–3145,, 2017. a

Bates, P. and De Roo, A.: A simple raster-based model for flood inundation simulation, J. Hydrol., 236, 54–77,, 2000. a, b

Bates, P. D., Horritt, M. S., and Fewtrell, T. J.: A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling, J. Hydrol., 387, 33–45,, 2010. a, b, c, d, e

Beven, K. J. and Kirkby, M. J.: A physically based, variable contributing area model of basin hydrology/Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant, Hydrol. Sci. B., 24, 43–69,, 1979. a

Christchurch City Council: Christchurch District Plan, Site of Ecological Significance: Site Significance Statement – Avon Heathcote Estuary/Ihutai & Environs, Tech. rep., Christchurch City Council, of Ecological Significance/SES LP 14.pdf (last access: 5 May 2023), 2016. a

Coulthard, T. J.: CAESAR-Lisflood landscape evolution model version 1.9j, (last access: 8 April 2023), 2019. a

Coulthard, T. J. and Macklin, M. G.: Modeling long-term contamination in river systems from historical metal mining, Geology, 31, 451,<0451:MLCIRS>2.0.CO;2, 2003. a

Coulthard, T. J. and Skinner, C. J.: The sensitivity of landscape evolution models to spatial and temporal rainfall resolution, Earth Surf. Dynam., 4, 757–771,, 2016. a

Coulthard, T. J. and Van De Wiel, M. J.: Modelling long term basin scale sediment connectivity, driven by spatial land use changes, Geomorphology, 277, 265–281,, 2017. a

Coulthard, T. J., Macklin, M. G., and Kirkby, M. J.: A cellular model of Holocene upland river basin and alluvial fan evolution, Earth Surf. Proc. Land., 27, 269–288,, 2002. a

Coulthard, T. J., Neal, J. C., Bates, P. D., Ramirez, J., de Almeida, G. A. M., and Hancock, G. R.: Integrating the LISFLOOD-FP 2D hydrodynamic model with the CAESAR model: implications for modelling landscape evolution, Earth Surf. Proc. Land., 38, 1897–1906,, 2013. a, b, c

de Almeida, G. A. M. and Bates, P.: Applicability of the local inertial approximation of the shallow water equations to flood modeling, Water Resour. Res., 49, 4833–4844,, 2013. a

de Almeida, G. A. M., Bates, P., Freer, J. E., and Souvignet, M.: Improving the stability of a simple formulation of the shallow water equations for 2-D flood modeling, Water Resour. Res., 48, W05528,, 2012. a

Dottori, F., Salamon, P., Bianchi, A., Alfieri, L., Hirpa, F. A., and Feyen, L.: Development and evaluation of a framework for global flood hazard mapping, Adv. Water Resour., 94, 87–102,, 2016. a

Dottori, F., Alfieri, L., Bianchi, A., Skoien, J., and Salamon, P.: A new dataset of river flood hazard maps for Europe and the Mediterranean Basin, Earth Syst. Sci. Data, 14, 1549–1569,, 2022. a

Environment Agency: Cumbria Floods Technical Report – Factual report on meteorology, hydrology and impacts of January 2005 flooding in Cumbria., Tech. rep., Environment Agency Bristol, UK, 2006. a, b

Environment Agency: National LiDAR Programme DSM, Environment Agency [data set], (last access: 5 May 2023), 2020. a, b

Environment Agency: Hydrology Data Explorer, Environment Agency [data set], (last access: 5 May 2023), 2021. a, b

Environment Canterbury: River flow data, Environment Canterbury [data set], (last access: 8 April 2023), 2023. a, b, c

Findlay, R. H. and Kirk, R. M.: Post‐1847 changes in the Avon‐Heathcote Estuary, Christchurch: A study of the effect of urban development around a tidal estuary, New Zeal. J. Mar. Fresh., 22, 101–127,, 1988. a

Galland, J.-C., Goutal, N., and Hervouet, J.-M.: TELEMAC: A new numerical model for solving shallow water equations, Adv. Water Resour., 14, 138–148,, 1991. a

Haddadchi, A., Hicks, M., Olley, J. M., Singh, S., and Srinivasan, M.: Grid‐based sediment tracing approach to determine sediment sources, Land Degrad. Dev., 30, 2088–2106,, 2019. a

Hansen, D. and Rattray, M.: New dimensions in estuary classification, Limnol. Oceanogr., 11, 319–326,, 1966. a

Horritt, M. S., Bates, P. D., Fewtrell, T. J., Mason, D. C., and Wilson, M. D.: Modelling the hydraulics of the Carlisle 2005 flood event, P. I. Civil Eng.-Wat. M., 163, 273–281,, 2010. a, b

Hunter, N. M., Bates, P. D., Horritt, M. S., and Wilson, M. D.: Simple spatially-distributed models for predicting flood inundation: A review, Geomorphology, 90, 208–225,, 2007. a

Kopmann, R. and Markofsky, M.: Three-dimensional water quality modelling with TELEMAC-3D, Hydrol. Process., 14, 2279–2292,<2279::AID-HYP28>3.0.CO;2-7, 2000. a

Land Information New Zealand: Christchurch and Selwyn, Canterbury, New Zealand 2015, OpenTopography [data set], (last access: 5 May 2023), 2015. a, b

LAWA: Land, Air, Water Aotearoa,, last access: 8 April 2023. a

Maia, P. D., Maurice, L., Tessier, E., Amouroux, D., Cossa, D., Pérez, M., Moreira-Turcq, P., and Rhéault, I.: Mercury distribution and exchanges between the Amazon River and connected floodplain lakes., Sci. Total Environ., 407, 6073–6084,, 2009. a

Maurice-Bourgoin, L., Quemerais, B., Moreira-Turcq, P., and Seyler, P.: Transport, distribution and speciation of mercury in the Amazon River at the confluence of black and white waters of the Negro and Solimões Rivers, Hydrol. Process., 17, 1405–1417,, 2003. a

McMurtrie, S.: Heavy Metals in Christchurch Fish and Shellfish: 2014 Survey, Tech. Rep. 08002-ENV01-04, EOS Ecology, New Zealand, (last access: 5 May 2023), 2015. a

Melack, J. M. and Hess, L. L.: Remote sensing of the distribution and extent of wetlands in the amazon basin, in: Amazonian Floodplain Forests, edited by: Junk, W. J., Piedade, M. T. F., Wittmann, F., Schöngart, J., and Parolin, P., vol. 210, Ecological Studies, 43–59, Springer Netherlands, Dordrecht,, 2011. a

Mitsch, W. J. and Gosselink, J. G.: Wetlands, Wiley, Hoboken, NJ, 5th edn., 752 pp., ISBN 978-1-118-67682-0, 2015. a

Neal, J., Villanueva, I., Wright, N., Willis, T., Fewtrell, T., and Bates, P.: How much physical complexity is needed to model flood inundation?, Hydrol. Process., 26, 2264–2282,, 2012. a

Neal, J., Dunne, T., Sampson, C., Smith, A., and Bates, P.: Optimisation of the two-dimensional hydraulic model LISFLOOD-FP for CPU architecture, Environ. Modell. Softw., 107, 148–157,, 2018. a

Neal, J. C., Bates, P. D., Fewtrell, T. J., Hunter, N. M., Wilson, M. D., and Horritt, M. S.: Distributed whole city water level measurements from the Carlisle 2005 urban flood event and comparison with hydraulic model simulations, J. Hydrol., 368, 42–55,, 2009. a

O'Loughlin, F., Paiva, R., Durand, M., Alsdorf, D., and Bates, P.: A multi-sensor approach towards a global vegetation corrected SRTM DEM product, Remote Sens. Environ., 182, 49–59,, 2016. a, b

O'Loughlin, F., Neal, J., Schumann, G., Beighley, E., and Bates, P.: A LISFLOOD-FP hydraulic model of the middle reach of the Congo, J. Hydrol., 580, 124203,, 2020. a

Qi, W., Ma, C., Xu, H., Chen, Z., Zhao, K., and Han, H.: Low Impact Development Measures Spatial Arrangement for Urban Flood Mitigation: An Exploratory Optimal Framework based on Source Tracking, Water Resour. Manag., 35, 3755–3770,, 2021. a, b

Qi, W., Ma, C., Xu, H., and Zhao, K.: Urban flood response analysis for designed rainstorms with different characteristics based on a tracer-aided modeling simulation, J. Clean. Prod., 355, 131797,, 2022. a

Richey, J. E., Melack, J. M., Aufdenkampe, A. K., Ballester, V. M., and Hess, L. L.: Outgassing from Amazonian rivers and wetlands as a large tropical source of atmospheric CO2., Nature, 416, 617–620,, 2002. a

Robins, P. E., Lewis, M. J., Simpson, J. H., Howlett, E. R., and Malham, S. K.: Future variability of solute transport in a macrotidal estuary, Estuar. Coast. Shelf S., 151, 88–99,, 2014. a

Robins, P. E., Farkas, K., Cooper, D., Malham, S. K., and Jones, D. L.: Viral dispersal in the coastal zone: A method to quantify water quality risk., Environ. Int., 126, 430–442,, 2019. a, b

Sampson, C. C., Smith, A. M., Bates, P. D., Neal, J. C., Alfieri, L., and Freer, J. E.: A high-resolution global flood hazard model., Water Resour. Res., 51, 7358–7381,, 2015. a

Shaw, J., Kesserwani, G., Neal, J., Bates, P., and Sharifian, M. K.: LISFLOOD-FP 8.0: the new discontinuous Galerkin shallow-water solver for multi-core CPUs and GPUs, Geosci. Model Dev., 14, 3577–3602,, 2021. a

Trigg, M. A., Wilson, M. D., Bates, P. D., Horritt, M. S., Alsdorf, D. E., Forsberg, B. R., and Vega, M. C.: Amazon flood wave hydraulics, J. Hydrol., 374, 92–105,, 2009. a

Van De Wiel, M., Coulthard, T., Macklin, M., and Lewin, J.: Embedding reach-scale fluvial dynamics within the CAESAR cellular automaton landscape evolution model, Geomorphology, 90, 283–301, 2007. a

van der Peet, M. and Measures, R.: Avon-Heathcote Tidal Barrier Pre-Feasibility Study, Tech. rep., Christchurch City Council, Christchurch, New Zealand, (last access: 5 May 2023), 2015. a

Vauchel, P., Santini, W., Guyot, J. L., Moquet, J. S., Martinez, J. M., Espinoza, J. C., Baby, P., Fuertes, O., Noriega, L., Puita, O., Sondag, F., Fraizy, P., Armijos, E., Cochonneau, G., Timouk, F., de Oliveira, E., Filizola, N., Molina, J., and Ronchail, J.: A reassessment of the suspended sediment load in the Madeira River basin from the Andes of Peru and Bolivia to the Amazon River in Brazil, based on 10 years of data from the HYBAM monitoring programme, J. Hydrol., 553, 35–48,, 2017.  a

Vopel, K., Wilson, P. S., and Zeldis, J.: Sediment-seawater solute flux in a polluted New Zealand estuary, Mar. Pollut. Bull., 64, 2885–2891,, 2012. a, b

Vousdoukas, M. I., Bouziotas, D., Giardino, A., Bouwer, L. M., Mentaschi, L., Voukouvalas, E., and Feyen, L.: Understanding epistemic uncertainty in large-scale coastal flood risk assessment for present and future climates, Nat. Hazards Earth Syst. Sci., 18, 2127–2142,, 2018. a

Wilson, M. and Coulthard, T.: Tracing and visualisation of contributing water sources in a model of flood inundation. GeoComputation 2019, The University of Auckland,, 2019. a, b

Wilson, M., Bates, P., Alsdorf, D., Forsberg, B., Horritt, M., Melack, J., Frappart, F., and Famiglietti, J.: Modeling large-scale inundation of Amazonian seasonally flooded wetlands, Geophys. Res. Lett., 34, L15404,, 2007. a, b, c, d

Wilson, M. D. and Coulthard, T. J.: Tracing and visualisation of contributing water sources in a model of flood inundation: video supplement, Zenodo [video],, 2021. a

Wilson, M. D. and Coulthard, T. J.: CAESAR-Lisflood v1.9j-WS (water source tracing and visualisation), Zenodo [code and data set],, 2023. a

Wing, O. E. J., Bates, P. D., Sampson, C. C., Smith, A. M., Johnson, K. A., and Erickson, T. A.: Validation of a 30 m resolution flood hazard model of the conterminous United States, Water Resour. Res., 53, 7968–7986,, 2017. a

Wittmann, F., Junk, W. J., and Piedade, M. T.: The várzea forests in Amazonia: flooding and the highly dynamic geomorphology interact with natural forest succession, Forest Ecol. Manag., 196, 199–212,, 2004. a

Woodley, K.: Avon-Heathcote Estuary/Ihutai New Zealand: Site Information Sheet on East Asian-Australasian Flyway Network Sites (SIS) – 2017 version, Tech. rep., Partnership for the East Asian – Australasian Flyway (EAAFP), (last access: 5 May 2023), 2018. a

Yu, D. and Coulthard, T. J.: Evaluating the importance of catchment hydrological parameters for urban surface water flood modelling using a simple hydro-inundation model, J. Hydrol., 524, 385–400,, 2015. a

Short summary
During flooding, the sources of water that inundate a location can influence impacts such as pollution. However, methods to trace water sources in flood events are currently only available in complex, computationally expensive hydraulic models. We propose a simplified method which can be added to efficient, reduced-complexity model codes, enabling an improved understanding of flood dynamics and its impacts. We demonstrate its application for three sites at a range of spatial and temporal scales.