Different ways to predict a flood

Last year, Bart de Beus was nominated for the 2025 CWI Best BSc Thesis in Applied Math Award.

Below are two photos of Itteren, a town in the south of the Netherlands. The first one shows the town as it looks normally. The second picture shows, well, quite an extreme: a flood. In 1995, this flood led to the evacuation of 250 000 people and a million animals. This was one of a handful of heavy floods in the 90's, so the Dutch regional water authorities (water boards, or "waterschappen" in Dutch) were tasked with calculating how often an area overflows. Based on these calculations, the Dutch government took better protective measures against floods. How do these calculations work?

Figure 1 and 2: The town of Itteren, Limburg, in 2016 (left; Van de Biezen, 2016) and 1995 (right; Rijkswaterstaat, n.d.)

These calculations are called flood frequency analyses. The objective: calculate if an area does not overflow "too often". Although floods are rarely desirable, they may be acceptable on occasion. For example, preventing a flood might be more expensive than dealing with the consequences if water reaches your house only once every hundred years. In this case, you may choose to let the flood happen, whenever that may be, and repair the damage when it does. Flood frequency analyses are done using one of two statistical methods; the "timeseries method" or the "stochastic method". Different water boards use different methods, the division being roughly 50/50. We don’t actually know, however, if these methods lead to different results. Spoiler alert: they do!

I had the opportunity to research both methods for Waternet, the company that works on behalf of water board Amstel, Gooi en Vecht. They use the timeseries method for the flood frequency analyses, but were also interested in implementing the stochastic method as there may be advantages to using one method over the other.

The objective

First, we must ask ourselves: how often is 'too often'? The Dutch provinces have decided on an acceptable return period $\tau$ , measured in years, for each type of area. According to the provinces, the type of area is allowed to flood once every $\tau$ years. For instance, the provincial bylaws note that a house is allowed to flood once every hundred years ( $\tau=100 \text{ years}$ ). By law, the water boards have to take measures making sure that water does not reach houses more often than that. Another example is farmland, which is allowed to flood only once every twenty-five years.

The water authorities use statistics and hydrological software to calculate the water level $w_\tau$ that corresponds to return period $\tau$ . What does this mean exactly? Say, we want to do a flood frequency analysis for farmland. That means that, given the provincial bylaws, we look for the water level corresponding to $\tau=25 \text{ years}$ . If the water level corresponding to $\tau =25 \text{ years}$ is higher than the level of the farmland, Waternet has to take measures to ensure this area overflows less often than once every tenty-five years. The number $\tau$ can take every value above zero, but common values (extracted from the bylaws) are ten, twenty-five, fifty and a hundred years.

A digital tool for mapping floods

The hydrological software D-Hydro can be used for both the stochastic and timeseries method. Using D-Hydro it is possible to map different areas into systems of canals, rivers, pumps, bridges and more. Figure 3 shows how this looks.

Figure 3: The Noordse Buurt and Westveense polder as simulated in D-Hydro.

D-Hydro can be given various input parameters and calculates water levels based on those inputs. Both the stochastic and timeseries method use the water levels calculated by D-Hydro. The difference between the methods is the input and the way its output is used to calculate how often an area floods.

The timeseries method

Let's start by looking at the timeseries method. In De Bilt the Royal Dutch Meteorological Institute (KNMI) records the amount of precipitation. These records are available as a dataset with hourly measurements, going back over a hundred years. These "rainfall events" are the input for D-Hydro, which then calculates the water levels over time.

After entering the rainfall events into D-Hydro and calculating the water levels, we select the highest water level per year. So, over the $n$ total years, we obtain with $n$ different maximum water levels. We denote these by $w_1, w_2, \dots, w_n$ . We call the maximum water level in a year $W$ , which is a random variable. The water levels $w_1, \dots, w_n$ are realizations of $W$ .

Ideally, we want to know the distribution of $W$ as it holds all kinds of relevant information. Extreme Value Statistics (EVS) can help with this. EVS has important results on the the maximum values of random variables, namely that the distribution of the maximum over many independent draws always converges to a specific probability distribution. If we assume that necessary conditions for EVS are met, we can use it to gain insight into the yearly water peak, in this case $W$ . EVS is a commonly used method for estimating flood events and it has many other applications.

Some extra clarification for the interested reader: As described, the first step is to take the $n$ different water levels, which are yearly maxima. We sort these water levels $w_1,\dots,w_n$ in a non-rising order. Denote by $m_k$ the rank of $w_k$ in the sorting. For example, if $w_3$ would be the highest water level recorded over all the years, $m_3=1$ . The second highest water level would get rank $2$ and so on. Next, we define

\begin{equation}
P_k = \frac{m_k}{n+1}, \quad, \text{for $k=1,\dots,n$}.
\end{equation}
Here, $P_k$ is an approximation of the probability of exceedance of $W$ . So,
\begin{equation}
P_k \approx \mathbb{P}(W > w_k),
\end{equation}
meaning the chance that a random maximum water level in an arbitrary year is higher than $w_k$ . Extreme value statistics tells us that, under certain assumptions, it holds that
\begin{equation}
1 - P_k = \mathbb{P}(W \leq w_k) \approx -e^{-e^{-(w_k-b)/a}},
\end{equation}
for some constants $a$ and $b$ . This distribution is known as the Gumble distribution. From this we can deduce that
\begin{equation}
w_k \approx a\big(-\ln(-\ln(1-P_k))\big) + b.
\end{equation}
If we define
\begin{equation}
Y_k = -\ln(-\ln(1-P_k)),
\end{equation}
we see that $w_k \approx aY_k + b$ , so it follows a linear trend.

In practice, we want to compute the water level $w_\tau$ for a given return period $\tau$ . We define $P_\tau = 1-e^{-1/\tau}$ and calculate $Y_\tau$ analogously to before, so

\begin{equation} Y_\tau = -\ln(1/\tau), \end{equation} and compute $w_\tau \approx aY_\tau + b.$

To use this result, however, we need to know $a$ and $b$ . A common approach is to use linear regression. In consultation with Waternet, a different approach was chosen: to do linear interpolation between the two surrounding datapoints.
Essentially, we plot the datapoints with $Y$ on the x-axis and $w$ on the y-axis. We connect the datapoints by straight line segments.
Precisely, to get $w_\tau$ , we compute
\begin{equation}
\hat{w}_\tau = w_L+\frac{Y_\tau-Y_L}{Y_R-Y_L}(w_R,w_L), \quad \tau>0,
\end{equation}
where $Y_L,Y_R$ are are values closest to $Y_\tau$ ( $Y_L$ on the left, $Y_R$ on the right hand side) and $w_L$ and $w_R$ the corresponding water levels ( $w_L$ is left, $w_R$ on the right hand side). Our $\hat{w}_\tau$ is an estimator of $w_\tau$ .

To give a more concrete, theoretical example: Say, we have a five-year series of rainfall events. From this five-year series, we get five different yearly maximums. See the table below for an example of these maximums, with the calculated $m_k$ , $P_k$ and $Y_k$ .

Input	Date	Water level (m+NAP)	Rank	Probability	Z-Score
1	2006	$-3.50$	2	$0.333$	$0.903$
2	2007	$-3.41$	1	$0.167$	$1.702$
3	2008	$-3.62$	5	$0.833$	$-0.583$
4	2009	$-3.58$	4	$0.667$	$-0.094$
5	2010	$-3.55$	3	$0.5$	$0.367$

We want to estimate the water level for a return period of $\tau=5 \text{ years}$ . This gives us a $P_\tau$ of $P_5 \approx 0.181$ , which gives $Y_\tau$ of $Y_5=1.609$ .

The closest $Y_k$ ’s are $Y_2$ on the right hand side and $Y_1$ on the left hand side. Filling this in, for our estimator $\hat{w}_\tau$ , gives: \begin{equation} \begin{split}\hat{w}_{\tau=5}&=-3.50+\frac{(1.609-0.903)}{(1.702-0.903)}(-3.41+3.50)\\&=-3.42.\end{split}\end{equation}

To conclude, using this theoretical dataset, we estimate the water level to reach $-3.42 \text{ m+NAP}$ every five years.

A note on assumptions: Here we assume that $W$ follows a Gumble distribution. Extreme Value Statistics tells us this is the case if certain assumptions are met. Let us look at some of those assumptions. The Fisher-Tippett-Gnedenko theorem or its more general counterparts gives us an asymptotic distribution of $W$ if the observations that induce $W$ (in this case the measurements of water levels) are independent. This is not necessarily the case, as the water level on a given day may be correlated with the water level on surrounding days. More recent work shows that the asymptotic behaviour still holds if the dependence is "weak enough". The theorem also gives other possible asymptotic distributions, other than Gumble. If water levels are heavy-tailed, for instance, another asymptotic distribution comes into play. The validity of the assumptions underlying this method are an important and interesting topic of discussion, but outside the scope of this article.

The stochastic method

With the stochastic method, we say that a whole plethora of explanatory variables can cause extremely high water levels. For example, if the water level in a river is extremely high, just a bit of rain can cause a flood. The same goes for if a pump fails: if water cannot be pumped out of a polder, canals might overflow, even without rain.

These (continuous) explanatory variables are discretized in different classes, after which every class gets a representative value. The goal being that every possible situation that might happen, gets simulated in D-Hydro. Every class gets allocated a probability, based on historical data from a case study area.

The advantage of this method is that we don’t just look at rainfall events to predict floods, but also look at other contributing factors. With enough historical data, we can add almost any contributing variable, like how high the groundwater level is, the water levels in neighboring rivers and lakes, and lastly also if certain pumping stations fail.

Let

$a_1, \dots, a_l$ be the explanatory variables. The probability that

$a_i$ takes the value

$a_{ij}$ , meaning it is in the

$j$ -th class, is denoted by by

$\mathbb{P}(a_i = a_{ij}).$ These are approximated using historical data.

The probability that an input vector $v$ occurs is approximated by multiplying the probabilities that variable $a_i$ falls in class $v_i$ .

\begin{equation} P_v = \mathbb{P}(a_1 = v_1) \cdots \mathbb{P}(a_l = v_l). \end{equation}

In other words, we assume that the random variables $a_1,\dots,a_l$ are independent.

For every possible input vector $v$ , we perform a simulation in D-Hydro, which gives the corresponding water level $w_v$ . Using these, we can compute the probability of exceeding a given water level $w$ as follows:

\begin{equation} \mathbb{P}(W \geq w) = \sum_{v:w_v \geq w} P_v. \end{equation}

Here $\mathbb{P}(W \leq w)$ is the approximation of the chance that the water level exceeds level $w$ at some point within one year. The return period can be calculated:

\begin{equation} \tau_v = \frac{1}{\mathbb{P}(W \geq w_v)} \end{equation}

Within this method we also use interpolation to estimate $w_\tau$ . Here we don’t interpolate based on a reduced variable (like we did with the timeseries method), but on our previously calculated $T_L$ and $T_R$ . So

\begin{equation} \hat{w}_\tau = w_L + \frac{\tau - T_L}{T_R-T_L}(w_R - w_L).\end{equation}

To give a more concrete, theoretical example: Say we have two theoretical variables $a_1 \in \{\text{summer}, \text{winter}\}$ and $a_2 \in \{\leq 20 \text{ millimeters}, > 20 \text{ millimeters}\}.$ The probability of it being summer or winter is equal, 0.5, for either both. However, the chance of less than 20 millimeters of precipitation is 0.7, and more then 20 millimeters is 0.3. This gives us the following table:

Input vector $v$	$P_v$	$w_v$
(summer,≤20)	0.35	−3.62
(summer,>20)	0.15	−3.55
(winter,≤20)	0.35	−3.65
(winter,>20)	0.15	−3.59

To calculate our $P(W\geq w)$ , we simply add the calculated probabilities of the water levels that are higher then the $w$ we seek. So, for $v = (\text{summer}, \leq 20)$ , the water levels that are higher or equal, are $w_{(\text{summer}, \leq 20)}$ itself, $w_{(\text{summer}, > 20)}$ and $w_{(\text{winter}, > 20)}$ . Which gives $P_v =0.35+0.15+0.15=0.65.$

Usually we have a couple of hundreds inputs, as opposed to only a few in this theoretical example. So we cannot really estimate our $w_\tau$ 's in this example. However, it works in the same way our timeseries method works.

What does the result look like?

Both methods ideally lead to the same result, which is a graph with on the logarithmical x-axis our return periods (the $\tau$ we are after), and on the y-axis the water level (the $w_\tau$ we seek). In Figure 4, the result is shown for one area (compartment OP-1-A) in the Noordse Buurt and Westveense polder, our case study region. In every of the thirty compartments in our case study, the timeseries method lead to higher predicted water level compared to the stochastic method. The average difference over all compartments was 4 to 9 centimeters, with $\tau$ ranging from 25 to a 100 years. The maximum difference was over 28 centimeters.

Figure 4: Results of the stochastic method (red) and the time series method (blue) for compartment OP-1-A in the Noordse Buurt and Westveense polder.

Now knowing that these two methods lead to varying results, it is interesting to know why we get these varying results. In turn, this could help the water boards make more informed decisions in what measures to take and when. During my research at Waternet, I investigated multiple explanations for these differences, about which I write in my thesis Schatten van wateroverlast, een vergelijking tussen statistiek vóóraf en statistiek áchteraf in de Noordse Buurt en Westveense polder (De Beus, 2025). As you can tell, it is written in Dutch.

My experience at Waternet

As an applied mathematician, I never imagined myself suddenly doing my thesis at a company which regulates and takes care of our water supply and systems. It is a fun example of how mathematics is almost everywhere! I found myself in a team of hydrologists, which is a whole other area of expertise, I thought. However, all models they use (including D-Hydro) were based on, as you might guess, mathematical models and formulas. Furthermore, these two methods I researched are largely mathematical problems. The hydrological part is mapping the water system in D-Hydro, and knowing how to take correct measures.

The models I implemented were immediately used for other case study areas. Through my research I hoped to give better insight on how and when to implement which method. Most importantly, I wanted to shed light on some of their drawbacks and when these drawbacks influence the results.

Cover photo taken from Orange County Register, June 19th, 2008

Would you like to stay up to date when a new post appears on the Network Pages?

Subscribe to our mailing list

The objective

A digital tool for mapping floods

The timeseries method

The stochastic method

What does the result look like?

My experience at Waternet

Related articles

Catching plastic waste from canals•

Gut flow: key to studying the microbiome••

Tales of the Tails: (Not so) mysterious Heavy Tail worldbeyond the Bell Curve••

Prizes for outstanding mathematics (Part 2)•