Estimating spatial flooding risk in Ellerslie

This repo contains code to simulate rainfall related flooding across a small open source GIS dataset sourced from Auckland Council. The purpose of this project was to learn how to parse, model and visualise GIS data and derived measures. I have no prior experience with GIS data or models, but have a reasonable grasp of foundational mathematical modelling so thought it would be interesting to learn about specific GIS applications.

I implemented a flux based discretisation of the Transport equation across a regular grid, using the Forward Euler method to solve forwards in time. The velocity field which determines how water volume moved across the land area was computed used Manning flow on a sloped surface, while the digital elevation map was interpolated from open source contour data from the Ellerslie area in Auckland. Rainfall intensity for a 6H storm event was calculated using a cumulative density function based on open source data from NIWA, from the Alexander Park Race Course. Spatial risk was calculated with a flooding threshold of 0.4 m above ground level, the Monta Carlo method was used to model the variation in possible rainfall intensity, with 150 6-hour storm events simulated to compute the spatial risk map.

I chose a small section of land 600 m in the E-W by 370 m in the N-S direction. The land is between Peach Parade and Green Lane East, with the Ellerslie Race Course to the south. The approximate area is shown courtsey of Google Maps below. Low resolution satellite imagery of the area of interest

I was interested in exploring flooding risk in the area as I've driven through many a puddle on Peach Parade in times of rain. Auckland Council provided a rich variety of data (https://geomapspublic.aucklandcouncil.govt.nz/) that was present across the area, but for my purposes I just wanted the contours, the selected shape, and the building outlines. I wasn't setting out to build a professional mathematical model, so was not planning on incorporating soil dynamics or drainage.

Fortunately, the data is provided in the New Zealand Transverse Mercator 2000 coordinate reference system, so axes are N/S and E/W aligned and have a spacing of metres. Parsing the contours into Python the metadata indicated they were sourced from a 2024 Lidar survey, which explains the density of the data. Given my own computational limitations I sampled 5000 points across the line sequences that were defined in the contours.shp file, then used radial basis functions to interpolate across a regular grid so that I had a regular digital elevation map with spacing of 10 m in both axes. To make up for the relatively sparse sampling of the very dense contour data, the field was smoothed with a Gaussian kernel. Digital elevation map interpolated from the contour data, smoothed with a Gaussian Kernel.

I was interested in examining how rainfall moves and accumulates across a surface geometry. My initial feeling was that an appropriate model would be something like a 2D Poisson equation, where rainfall acts as a source, and drainage would behave like a sink term in the domain. However, lacking advective terms or a time derivative makes the equation ill-suited for my needs. I settled on the transport/continuity equation which is a simple PDE that describes local conservation of mass for some quantity ($C$). Assuming water is incompressible, the spatial velocity field and the source term ($S$) are a function of space and time:

$\frac{\partial C}{\partial t} + \frac{\partial C}{\partial x}u + \frac{\partial C}{\partial y}v = S(x,y,t)$

Where $u$ and $v$ are the x, and y components of the velocity fields respectively.

I initially chose to model this with an upwind finite difference scheme with forward Euler time stepping and first order upwards difference for spatial derivatives.

$\frac{\partial C}{\partial t} \approx \frac{C_{i,j}^{n+1} - C_{i,j}^n}{\Delta t}$

$\frac{\partial (u C)}{\partial x} \approx\dfrac{u_{i,j}^n C_{i,j}^n - u_{i-1,j}^n C_{i-1,j}^n}{\Delta x}$ Where $u_{i,j}^n \geq 0$ Alternatively when $u_{i,j}^n < 0$ $\frac{\partial (u C)}{\partial x} \approx\dfrac{u_{i+1,j}^n C_{i+1,j}^n - u_{i,j}^n C_{i,j}^n}{\Delta x}$

Discretising at the nodes in this way does not strictly enforce conservation of mass and was unsurprisingly very numerically unstable. Re-discretising the equations in terms of interface fluxes between control volumes across cells in the grid was very helpful for model stability. For a source term I used homogenous rainfall across the domain, effectively multiplied by a step function in time that set the source everywhere to zero after the rainfall period was complete. Fluxes were set to zero at the boundary interfaces.

My initial attempts to calculate a velocity field across the domain used Poiseuille flow, a simplification of the Navier-Stokes equation under the assumptions of incompressibility, laminar flow, and fully developed flow across the domain. For a flat plane with a water film of thickness $h$, this flow resolves to: $V=g*sin(θ)*h^2 /3ν$, where $g$ is acceleration due to gravity at $9.81~ms^-2$, $θ$ is the slope angle,$h$ is the height of the water film across the cell, and $v$ is the kinematic viscosity of water. As my model only models the advection of water, water moves from high to low elevations, but low points with high volumes of water do not spill or diffuse into adjacent cells. This means that low cells have VERY high water volumes, which unsurprisingly break the stability of the model as when using Poiseuille flow, velocity scales quadratically with the fluid height.

A quick google lead me to Manning flow, an empirical relationship (presented in SI units): $V = \frac{1}{n}R_{h}^{2/3}S^{1/2}$ Where $n$ is the Manning roughness coefficient, it is approximately 0.012 for concrete, $R_{h}$ is the hydraulic radius which is approximately $h$ for uniform sheet flow, $S$ is the channel bed slope which is approximately $\text{tan}(\theta)$. This reduced my velocity enough to be physically reasonable and within the CFL stability condition.

Fluid flow across the surface of synethic valley geometry, initialised with 1 mm of fluid across the surface

Using publicly available Intensity Duration data from NIWA (https://data.niwa.co.nz/hirds) sourced from the nearby Alexandra Park Raceway for a 6-hour storm event. I computed a simple linearly interpolated cumulative density function of rainfall intensity and sampled 150 times from it to parameterise a Monte-Carlo simulation of rainfall across my land area.

As I only modelled advection of water across the space, I needed to redistribute the water across cells to calculate the water level for any given basin or basins in the domain. This was done using a connected spatial graph with elevation and volumes as nodal properties and conservation of mass strictly enforced. I used a variation of a flood fill algorithm that recursively solved for the water level across connected inundated cells.

Using a water level of 0.4 m above the surface as a threshold for 'flooded' - I looked at the probability that a grid cell of 10 by 10 m would flood in a given year. 0.4 metres is relatively arbitrary, but it seems from casual reading to be about the point where vehicles would begin to float and structures begin to weaken. Flood risk across the simulation grid for a 6-hour storm event, overlaid with the regions contour and building outline data

We can see that for a given year the model predicts that the central block of houses will flood due to a 6-hour storm event 65% of the time. There are a few reasons that this shouldn't be alarming. The model does not take into account any soil dynamics, so all water that lands within my land area will contribute to the final surface water calculations. I also assume that all boundary cells have zero flux, which is essentially saying that no matter how flooded my boundaries are, the neighbouring cells are also that flooded - this deprives the model of what is effectively the reservoir of the Ellerslie Race Course, although it should be noted the building outlines do not show the newly built apartment along Peach Parade. Without modelling any sink terms in the model, i.e. no stormwater clearance system, or a stormwater clearance system where the local network capacity is overloaded, the model predicts significantly elevated flood risk.

The limitations of the model, that of impervious surfaces, no functioning stormwater clearance and local reservoirs at capacity closely align with worst practices for stormwater management compared to the reccomendations given in this Auckland Council planning document (see here). Skimming the document I was not able to gauge a simple estimate of local stormwater capacity for the Ellerslie zone, and looking at the influence of stormwater clearance through the drain network would be a natural extension of the model and simple to implement the model is numerically stable for sink terms already.