Estimating flood hazard curves with a hierarchical log-increment model

[jackelsey]

I’m working for a US government program that is developing a new process for generating maps of flooding hazards. I’d appreciate feedback on the Bayesian formulation below, particularly its separation of aleatoric and epistemic uncertainty.

The current process involves selecting a single historical storm as a representative of a 100-year rainfall event, running that rainfall data through a mechanistic flood simulation model, then designating all areas that experienced simulated flooding as within the 100-year floodplain. Apart from ignoring many sources of uncertainty, the main practical issue with this approach is that characterization of risk is binary. You’re either within the 100-year floodplain or not; there is no gradation in estimated hazards.

The new process we’re developing will estimate a flood hazard curve with an uncertainty band for every location in the map. This is a curve where flooding depth at that location is on the y-axis, and annual exceedance probability (AEP) is on the x-axis. The depth resulting from a 100-year flood is given by the curve’s value at an AEP of 1/100. (This information could also be expressed with a CDF, but the engineers are used to flood hazard curves.)

We first tried brute-force Monte Carlo simulation, but this strategy is too computationally expensive for wide-scale implementation. The brute-force approach would generate R=50 realizations of synthetic rainfall data, where each realization contains Y=2000 years, and each year contains N\approx10 large rainfall events. Aleatoric sources of uncertainty (e.g. rainfall event magnitude and location, antecedent moisture conditions) would vary from event to event, and a plausible set of mechanistic flood simulation model parameters \theta (derived from fitting the model to historical flooding event data) would be assigned to each realization to account for epistemic uncertainty. This process would yield R empirical flood hazard curves at each location and allow for simple estimation of uncertainty bands, but the cloud computing costs of running RYN 10-15 hour mechanistic flood simulations for every watershed in the US was prohibitive.

Our current objective is to find ways of approximating the output of brute-force Monte Carlo at a lower computational cost. Other people on the project team are testing strategies involving importance sampling and surrogate models; I’m developing a Bayesian model with Stan which can leverage output from a much cheaper mechanistic model for predicting river flowrates. Here’s an overview of my Bayesian approach:

First, any flood hazard curve can be expressed as a cumulative sum of depth increases.

Assume that the log-transformed depth increment at AEP index i and location j is normally distributed.

y_{ij} \sim N(\eta_{ij}, \sigma _i )

The mean log-increment is a sum of fixed and random effects.

\eta_{ij} = \mu + c_i + a_j + d_{ij}

• \mu is an intercept
• c_i is a random effect capturing how some locations (cells in a simulation mesh) have consistently larger or smaller log-increments
• a_j is a fixed effect describing how some AEPs have consistently larger or smaller log-increments
• d_{ij} is a random effect that describes remaining variability in hazard curve shapes

I found that when you take the flowrate predictions at several junctures on the river, perform PCA, sort the rainfall events by the first principal component score, and then sub-sample events centered around a PC1 score percentile corresponding with an AEP of interest, the flooding depths are biased from the empirical hazard curve depths in a predictable way.

Two types of observations would be used to fit the model:

• Depths from empirical hazard curves derived from running full realizations
• Approximated depths from the PCA sub-sampling approach

I still need to figure out how to account for epistemic uncertainty. One possible approach would select 3-5 target AEPs, sub-sample 30–50 rainfall events around each corresponding PC1 percentile, and rerun those events using R plausible parameter vectors \theta_r. Holding the rainfall event fixed while varying \theta_r would reveal the additional variance attributable to epistemic uncertainty. A Bayesian model could estimate how this variance changes with AEP and use it to widen the hazard-curve uncertainty band.

Here are some questions. I’m open to any input.

1. Am I overcomplicating this in some way? I originally suggested using 2D splines to model a GEV distribution with spatially variable parameters, but was told that complex patterns in flooding depth (e.g. flat expanses in hazard curves followed by a sudden jump when a levee is overtopped) necessitate a more nonparametric approach. Is there maybe some clever way to use bootstrapping?
2. Are there any other ways one could use to nonparametrically model upper tails of a distribution? I’ve never seen this log-increment approach used anywhere.
3. Are there any nested or block designs that would allow me to more effectively separate out epistemic uncertainty and Monte Carlo noise in the empirical hazard curves?

[kyle.hofmann]

Take everything I say with a grain of salt; I’m not a statistician, and I don’t know anything about your problem domain. That said, here are some thoughts.

The plots you showed of flood depth versus annual exceedance probability look, by eye, close to logarithmic. Since the horizontal axis is logarithmic, that suggests that the true behavior is roughly linear. It’s surely not exactly linear, but maybe it’s close enough that that’s somehow a useful set of coordinates to be thinking in: Perhaps the upper tail can be modeled as something linear (which could be fit using data not in the tail) plus a residual (which, if all goes well, is somehow easier to model). Perhaps it’s convenient to fit a linear model prior to taking logs and look at logarithmically spaced residuals.

A depth increment, if I’m understanding you right, is a (discrete approximation to a) derivative. The derivative of a curve tends to be more difficult to estimate than the curve itself because any little wiggle in the original curve gets magnified. On the other hand, a good physical model for predicting depth increments is likely to perform better than a bad physical model for predicting total depth. You haven’t specified how your fixed and random effects are computed (and I’m guessing it’s complicated), but if they naturally produce total depths, then it might be better to stick with those rather than work with increments.

It often happens in PCA that there are a number of components with non-trivial coefficients and then a lot of noise. Possibly that is what you are seeing here: The first PCA component is the most important, but perhaps there are another five, ten, etc. which are not noise. Including those should make your empirical and approximate depths match better.

The fact that the first PCA component does so well suggests to me that the cells in the simulation are correlated: When some have a large flow, then others also have a large flow. This could be extremely useful! It suggests that instead of modeling cell-by-cell, you could instead model “flow-by-flow”, i.e., by PCA components (or some other method of dimension reduction). For this to be effective, you would need to have a way to find the flows without a huge amount of simulation. I have not thought about this; it might require domain-specific knowledge or some kind of two-part strategy (where you run a little bit of an expensive simulation to figure out what components you need to keep and then only simulate those for the rest of the computation). This could be interpreted as a kind of answer to your block design question.

[jackelsey]

Kyle, thank you very much for sharing your thoughts, and sorry that my response was delayed!

@kyle.hofmann

The plots you showed of flood depth versus annual exceedance probability look, by eye, close to logarithmic.

@kyle.hofmann

Perhaps the upper tail can be modeled as something linear (which could be fit using data not in the tail) plus a residual (which, if all goes well, is somehow easier to model).

The example flood hazard curves I shared were for a relatively undeveloped watershed without flood control structures (e.g., levees). For other watersheds, the flood hazard curves will have more irregular shapes. For example, they may exhibit almost flat regions where the flood depth doesn’t increase until some nearby levee is overtopped. I’ll definitely keep your suggestions in mind, but I’m wary of introducing parametric assumptions regarding the tail shape.

@kyle.hofmann

A depth increment, if I’m understanding you right, is a (discrete approximation to a) derivative.

You make an interesting point. Modeling the flood hazard curve as the cumulative sum of exponentiated normal variables was just my (possibly crude) way of enforcing monotonicity. Maybe there’s a better way of doing that with splines or Gaussian processes…

@kyle.hofmann

It often happens in PCA that there are a number of components with non-trivial coefficients and then a lot of noise. Possibly that is what you are seeing here: The first PCA component is the most important, but perhaps there are another five, ten, etc. which are not noise. Including those should make your empirical and approximate depths match better.

I’m only using the first component, because I needed a single metric to rank all the simulations from the cheaper mechanistic model for predicting river flowrates. I’m not sure how that would work if I used more than one component?

@kyle.hofmann

The fact that the first PCA component does so well suggests to me that the cells in the simulation are correlated: When some have a large flow, then others also have a large flow. This could be extremely useful! It suggests that instead of modeling cell-by-cell, you could instead model “flow-by-flow”, i.e., by PCA components (or some other method of dimension reduction).

This idea is similar to what another project team did to develop a surrogate model. They implemented the approach in this paper for “upskilling” a low-resolution version of the mechanistic flood simulation model with PCA and GPR.

[kyle.hofmann]

sorry that my response was delayed!

No worries! Actually the delay might have worked out. In the meantime, I went to a seminar talk that could be relevant to what you’re doing.

The speaker was Dimitrios Psaltis, an astrophysicist from Georgia Tech. One of the things he described in his talk was automatically discovering equations for physical laws. He is interested in the magnetohydrodynamics of black holes. The equations governing them are, as normally written down, impossible to solve numerically in some conditions of interest: For his questions, there are important contributions from both very small and very large scales, and if you discretize enough to model the small scales, then there are not enough computational resources anywhere on Earth to reach the large scales. What one wants is some kind of new approach; that could be a different thing to model, or an additional conserved quantity, or something that can be integrated analytically, or maybe something else. Either way, the challenge is basically to discover new physics in an area that’s already been picked over.

Physical laws tend to be relatively simple; the number of terms they include is vastly smaller than the number of terms that one could write down. If one is sufficiently clever, it is possible to discover equations automatically. Several groups have worked on this, and his group made some advances which helped them model black holes. The web page for his work in this area is here. Perhaps this work would be useful to you.

For other watersheds, the flood hazard curves will have more irregular shapes. For example, they may exhibit almost flat regions where the flood depth doesn’t increase until some nearby levee is overtopped.

In that case, I agree that a parametric model is unlikely to perform well. The situation where a levee is overtopped looks difficult to model with any standard parametric method because the derivative is discontinuous. I suppose that one try to find where the derivative is discontinuous and then use a piecewise model with knots at the discontinuities; but that requires detecting the discontinuities. That’s basically a change-point detection problem. Plenty has been written about these, but I know someone who wrote a thesis in this area and he once described the literature as “chaos”.

Modeling the flood hazard curve as the cumulative sum of exponentiated normal variables was just my (possibly crude) way of enforcing monotonicity.

There is a thing called “monotone regression” or “monotonic regression”. I’m not overly familiar with it, but a colleague once pointed me to the paper: J. O. Ramsay, “Monotone Regression Splines in Action”, Statist. Sci. 3(4): 425-441 (November, 1988). DOI: 10.1214/ss/1177012761. To the extent to which I understand it, the idea is to model the logarithm of the derivative of a function using splines; as long as the function is strictly increasing (or, after flipping a sign, strictly decreasing), you get monotonicity for free.

If I recall correctly, there are situations where the MLE for a monotically increasing sequence of measurements is just the cumulative maximum. But regardless of whatever frequentist properties that estimate might have, it violates my prior that a monotonic function which I want to model is gonna be strictly increasing and have several derivatives almost everywhere.

I needed a single metric to rank all the simulations from the cheaper mechanistic model for predicting river flowrates. I’m not sure how that would work if I used more than one component?

You could take the length of the projection onto the first however many PCA components. But how many components to take is something I don’t know how to answer.

They implemented the approach in this paper for “upskilling” a low-resolution version of the mechanistic flood simulation model with PCA and GPR.

Oh, that sounds interesting. I’ll take a look!

[Bob_Carpenter]

@jackelsey

particularly its separation of aleatoric and epistemic uncertainty.

I wouldn’t restrict myself to these. One of the harder things is taking into account model uncertainty. Usually the true data generating process isn’t one of the models under consideration.

@jackelsey

You’re either within the 100-year floodplain or not; there is no gradation in estimated hazards.

Is that because of some downstream binary requirement for decision making?

@jackelsey

The new process we’re developing will estimate a flood hazard curve with an uncertainty band for every location in the map.

This seems like the right approach even if there is a downstream binary decision requirement. Why isn’t depth 0 at 100%?

@jackelsey

each realization contains Y=2000 𝑌 =2000 years, and each year contains N\approx10 𝑁 ≈10 large rainfall events

If you’re interested in 100-year floods, why are there 10 per year?

@jackelsey

too computationally expensive for wide-scale implementation.

How efficiently are you doing the compute and on what kind of hardware?

@jackelsey

importance sampling and surrogate models;

I would give the latter a much higher chance of success than the former. Importance sampling only works when the proposals are very close, which they rarely are in situations where people want to try importance sampling.

@jackelsey

Assume that the log-transformed depth increment at AEP index i 𝑖 and location j 𝑗 is normally distributed.

This is the same as taking the depth increment to be lognormally distributed.

In your model, there is a lot of non-identifiability to deal with, as you typically get in a mixed effect model like this. In particular, \mu, c, a, d are not identified, only heir sum is. That is, if I add a value to all of the a, I can subtract it from \mu and get the same answer. Mitigating this non-identifiability is the main thing you have to do for efficiency. One way to do this is to have the random effects sum-to-zero. That’s better and more symmetric for priors than using the traditional pinning of one value to zero.

@jackelsey

the flooding depths are biased from the empirical hazard curve depths in a predictable way.

Do you know where the bias is coming from? I’m not sure what you mean by bias being predictable. By definition, bias is expected error, so it always comes with a sign when it’s non-zero.

@jackelsey

2. I’ve never seen this log-increment approach used anywhere.

This is how Stan implements an ordered parameter under the hood. See:

Constraint Transforms

I used the underlying ordered vector type in my Covid sensitivity case study for ISBA last year to model the sensitivity of Covid tests as a function of time since symptom onset. Log is not the only function that maps (0, \infty) to (-\infty, \infty), so just like in a GLM, you can change the link function (e.g., softplus).

@jackelsey

3. nested or block designs that would allow me to more effectively separate out epistemic uncertainty

I’m not sure why you want to try to separate these things. You can take a given model, though, and look at what happens as you simulate more and more data. The “aleatoric uncertainty” (from only having a finite sample) will go from some value with no data at all to a value after you’ve seen all the data. In models where the parameter size doesn’t grow with the data and a few other conditions hold, the Bernstein-von Mises theorem ensures things are normal in the limit concentrating on the true parameter values. This is characterized in a single run by the Monte Carlo standard error (MCSE), which tells you how much uncertainty there is in an expectation estimate, like for the posterior mean of a parameter.