# Quantizing rare random maps: application to flooding visualization

Charlie Sire
IRSN, BRGM, Mines Saint-Etienne, Univ. Clermont Auvergne, CNRS, UMR 6158 LIMOS
and
Rodolphe Le Riche
Mines Saint-Etienne, Univ. Clermont Auvergne, CNRS, UMR 6158 LIMOS
and
Didier Rullière
Mines Saint-Etienne, Univ. Clermont Auvergne, CNRS, UMR 6158 LIMOS
and
Jérémy Rohmer
BRGM
and
Lucie Pheulpin
IRSN
and
Yann Richet
IRSN

###### Abstract

Visualization is an essential operation when assessing the risk of rare events
such as coastal or river floodings. The goal is to display a few prototype
events that best represent the probability law of the observed phenomenon, a
task known as quantization. It becomes a challenge when data is expensive to
generate and critical events are scarce, like extreme natural hazard. In the
case of floodings, each event relies on an expensive-to-evaluate hydraulic
simulator which takes as inputs offshore meteo-oceanic conditions and dyke
breach parameters to compute the water level map. In this article, Lloyd’s
algorithm, which classically serves to quantize data, is adapted to the
context of rare and costly-to-observe events. Low probability is treated
through importance sampling, while Functional Principal Component Analysis
combined with a Gaussian process deal with the costly hydraulic simulations.
The calculated prototype maps represent the probability distribution of the
flooding events in a minimal expected distance sense, and each is associated
to a probability mass. The method is first validated using a 2D analytical
model and then applied to a real coastal flooding scenario. The two sources of
error, the metamodel and the importance sampling, are evaluated to quantify
the precision of the method.

Keywords: Visualization, Quantization, Flood maps, Rare events, Metamodel

##  1 Introduction

Political authorities typically support their decisions regarding flood risk
management based of flood maps calculated with hydraulic simulations. Beven et
al. 2015 provide a realistic application, and a review of the existing methods
for graphical representations of flooding uncertainty is given by Seipel and
Lim 2017.

The design and analysis of numerical experiments is a generic procedure to
capture and explore the uncertainties related to such heavy numerical models
(Santner et al. 2003). It consists in performing a series of runs of the
hydrodynamic flow simulator where the forcing conditions (the river flow rate
and the offshore wave characteristics), the breach characteristics (location
and dimension of the breach) or the embankments characteristics (geometry,
failures conditions, etc.) are one realization of a random law. Rohmer et al.
2018 provide an example of the setting up of the forcing scenarios in coastal
simulations and specify a method to determine the probabilistic laws
describing the natural phenomena.

Visualizing simulations for risk assessment is a real challenge because of
their time-consuming aspect, the low probability of the flooding event, and
the need to communicate the results to authorities who do not necessarily have
a scientific background. To this end, a simple-but-efficient approach consists
in providing a small size set of representative flood maps. Such prototype
maps (Hastie et al. 2009, chap. 13) are chosen to best represent the
probability law associated to the flooding event. Calculating the prototype
maps is a quantization task (Gray and Neuhoff 1998), implemented here through
the Lloyd’s algorithm (Pagès 2014). Lloyd’s algorithm searches for Voronoi
cells centroids such that the mean distance (with respect to the uncertainty
distribution) between an observation and its nearest centroid is minimized.

This procedure yields results that have both a probabilistic meaning and can
be visually interpreted by a wide audience. For instance, in the flooding case
tackled in this article, the four non-empty prototype maps and their
probability masses expressed as frequencies are given in Figure 1. The empty
prototype map occurs with $97\%$ probability.

(a) $2.0\times 10^{-2}$, 1 in 50

(b) $5.2\times 10^{-2}$, 1 in 192

(c) $6.2\times 10^{-3}$, 1 in 161

(d) $1.4\times 10^{-3}$, 1 in 714

Figure 1: The four non-empty prototype maps, computed to minimize the expected
distance between a random map and its nearest prototype. The coast is at Les
Boucholeurs, France, and the maps are $1.6\times 1.6\leavevmode\nobreak\
\mathrm{km}^{2}$. The probability mass associated to each prototype map is the
probability that a random map is closer to this map than to the other
prototypes, for a given distance. The frequence, i.e., the inverse of the
probability is added for an easier interpretation. More details about this
case study are given in Section 6.

With this visualization, one can picture the most typical flooding scenarios,
like Figure 1(c) where water invades a large area, or more intense and
concentrated flooding like Figure 1(d).

Lloyd’s algorithm is not a new quantization technique (MacQueen 1967), and has
evolved since its introduction (Bock 2008). This paper presents new extensions
to quantization that address the following bottlenecks. First, the
quantization is done here in the context of costly numerical simulations,
which makes the computation of the centroids more difficult. Indeed, the flood
maps are computed with hydraulic codes that can take several hours to simulate
the water levels associated to each set of input parameters. Plus, flooding is
a rare event. Classical Monte Carlo sampling is inefficient for extreme events
with a very large return period (typically 1000-10000 years). Finally, the
quantization is performed in a space of pixelated maps ($64\times 64$) which
raises metamodeling and storage questions when a very large number $(>10^{6})$
of realizations has to be handled.

The paper is structured as follows. Section 2 defines the quantization problem
and the importance sampling estimator to compute the prototype maps. Section 3
introduces the metamodeling scheme built to overcome the duration of the
hydraulic simulations. Section 4 presents the algorithm and the performance
metrics that will be evaluated on the two applications described in Section 5
(analytical test case) and Section 6 (real coastal flooding). Finally, Section
7 discusses the main results and gives possible extensions to the method. A
list of symbols is also provided at the end.

##  2 Finding prototype maps

###  2.1 Quantization

Let us consider a random vector $X\in\mathcal{X}$, with $\mathcal{X}$ a subset
of $\mathbb{R}^{d},d\in\mathbb{N}$. In our context, $X$ is a parametric
description of the offshore hydrodynamic (like the tide and surge parameters)
and the dyke breach conditions. We define a random function $Y$ from
$\mathcal{X}$ the space of offshore conditions and breach parameters to
$\mathcal{Y}$ the space of possible flood maps, with an inner product
$\langle.,.\rangle_{\mathcal{Y}}$, and a metric $\lVert.\lVert_{\mathcal{Y}}$
induced by the inner product:

$\displaystyle Y\colon\mathcal{X}$ $\displaystyle\to\mathcal{Y}$
$\displaystyle x$ $\displaystyle\mapsto Y(x).$

In other terms, with this general formulation, $Y$ associates a distribution
of flood maps to every vector of input conditions $x$. In practice, $Y$ can be
deterministic and thus associate a single flood map to every $x$ which we
write $y(x)$.

The quantization problem is to find a set
$\Gamma_{\ell}=\\{\gamma_{1},\gamma_{2},\dots,\gamma_{\ell}\\}\in\mathcal{Y}^{\ell}$
of $\ell\in\mathbb{N}$ representatives of $Y(X)$ in a mathematically relevant
way. The notation $Y(X)$ should be understood as the map $Y$ conditioned by
the random vector $X$ as in regression, $Y(X)\equiv Y\\!\\!\mid\\!\\!X$. The
prototype maps will be these representatives in our study. More precisely, we
first define for $\Gamma_{\ell}\in\mathcal{Y}^{\ell}$ the closest candidate
representative as the function $q_{\Gamma_{\ell}}$:

$\displaystyle q_{\Gamma_{\ell}}\colon\mathcal{Y}$
$\displaystyle\to\Gamma_{\ell}$ $\displaystyle y$ $\displaystyle\mapsto
q_{\Gamma_{\ell}}(y)=\underset{\gamma_{i}\in\Gamma_{\ell}}{\operatorname*{arg\,min}}\lVert
y-\gamma_{i}\lVert_{\mathcal{Y}}\leavevmode\nobreak\ .$

A measure of the quantization error is the average distance between a randomly
chosen map and its closest representative (Pagès 2014),

$e(\Gamma_{\ell})=\left[\mathbb{E}\left[\lVert
Y(X)-q_{\Gamma_{\ell}}(Y(X))\lVert_{\mathcal{Y}}^{2}\right]\right]^{\frac{1}{2}}\leavevmode\nobreak\
.$

To each centroid $\gamma_{j}$ in $\Gamma_{\ell}$, we associate a Voronoi cell
$C^{\Gamma_{\ell}}_{j}$ which is the subset of points in $\mathcal{Y}$ closest
to $\gamma_{j}$,

$C^{\Gamma_{\ell}}_{j}=\\{y\in\mathcal{Y},q_{\Gamma_{\ell}}(y)=\gamma_{j}\\}\leavevmode\nobreak\
.$

The objective is to minimize this quantization error, i.e., to find

$\displaystyle\Gamma_{\ell}^{\star}=\\{\gamma^{\star}_{1},\dots,\gamma^{\star}_{\ell}\\}$
$\displaystyle\in\underset{\Gamma_{\ell}\in\mathcal{Y}^{\ell}}{\operatorname*{arg\,min}}\>(e(\Gamma_{\ell}))$
$\displaystyle\in\underset{\Gamma_{\ell}\in\mathcal{Y}^{\ell}}{\operatorname*{arg\,min}}\>\left[\mathbb{E}\left[\lVert
Y(X)-q_{\Gamma_{\ell}}(Y(X))\lVert_{\mathcal{Y}}^{2}\right]\right]^{\frac{1}{2}}$
$\displaystyle\in\underset{\Gamma_{\ell}\in\mathcal{Y}^{\ell}}{\operatorname*{arg\,min}}\>\left[\mathbb{E}\left[\underset{i\in\\{1\dots\ell\\}}{\min}\>\lVert
Y(X)-\gamma_{i}\lVert_{\mathcal{Y}}^{2}\right]\right]^{\frac{1}{2}}\leavevmode\nobreak\
.$

Finding the optimal centroids is in general a NP hard problem (Mahajan et al.
2012). However, the solution has a strong attribute (Pagès and Yu 2013):

######  Theorem 2.1 (Optimal quantization, Kieffer, Cuesta-Albertos).

If  $\bullet$ $\mathcal{Y}$ is of finite dimension m , $\bullet$ $\forall
y,y^{\prime}\in\mathcal{Y},\langle
y,y^{\prime}\rangle_{\mathcal{Y}}=\sum_{i=1}^{m}\lambda_{i}y_{i}y^{\prime}_{i}$
with $\forall i,\lambda_{i}>0$ , $\bullet$ $\mathbb{E}\left[\lVert
Y(X)\lVert_{\mathcal{Y}}^{2}\right]<+\infty$ ,
then $\forall i\in\\{1\dots\ell\\},\quad\mathbb{E}\left[Y(X)\mid Y(X)\in
C^{\Gamma_{\ell}^{\star}}_{i}\right]=\gamma^{\star}_{i}$ .

This means that the representatives of an optimal quantification coincide with
the Voronoi cell centroids. The proof is given in Appendix A. In our context,
the distance induced by the scalar product $\langle
y,y^{\prime}\rangle_{\mathcal{Y}}=\sum_{i=1}^{n}\lambda_{i}y_{i}y^{\prime}_{i}$
between two maps is the mean squared difference between the water depth at
each pixel of the maps, weighted by $\lambda_{i}$ which can represent the
importance of the zone (for example the housing density, the presence of
factories, etc.).

The popular Lloyd’s algorithm, also called K-means, is based on Theorem 2.1.

1:$\Gamma_{\ell}^{[0]}\leftarrow\\{\gamma^{[0]}_{1},\dots,\gamma^{[0]}_{\ell}\\}$
, $k\leftarrow 0$

2:while stopping criterion not met do

3:$\gamma^{[k+1]}_{j}\leftarrow\mathbb{E}\left[Y(X)\mid Y(X)\in
C^{\Gamma_{\ell}^{[k]}}_{j}\right],j\in\\{1,\dots,\ell\\}.$

4:$k\leftarrow k+1$

5:end while

Algorithm 1 Lloyd’s algorithm

The stopping criterion of the algorithm can be a number of iterations or a
lower bound on a distance between $\Gamma_{\ell}^{[k]}$ and
$\Gamma_{\ell}^{[k+1]}$. The principle of the iterations is illustrated in
Figure 2 in 2D: the dots are 5 possible starting points
$\\{\gamma_{1}^{[0]},\dots,\gamma_{5}^{[0]}\\}$, their associated Voronoi
cells are delimited by the black lines, the centroids of these cells are
computed for a uniform distribution and are drawn as the triangles. This
Figure shows only one iteration and does not represent the convergence of the
algorithm.

Figure 2: An example of Voronoi cells associated to initial points and their
centroids in 2D.

###  2.2 Dealing with low probability events

As stated before, Lloyd’s algorithm requires to compute at each iteration and
for each Voronoi cell $C^{\Gamma_{\ell}}_{j},j\in\\{1,\dots,\ell\\}$,
$\mathbb{E}\left[Y(X)\mid Y(X)\in C^{\Gamma_{\ell}^{[k]}}_{j}\right]$. The
most natural approach would be to sample $Y(X)$ from its density and then to
calculate the empirical mean of $Y(X)$ given $Y(X)\in
C^{\Gamma_{\ell}^{[k]}}_{j}$. However, flooding being a rare event, we expect
to have one Voronoi cell with almost only empty maps associated to a
probability close to 1, and other Voronoi cells with more critical floodings
associated to very low probabilities maps. Thus, the straightforward sampling
is not relevant to locate the Voronoi cells of low probabilities because their
member maps will almost surely be missed. Importance sampling is thus needed
to perform a sampling less concentrated on a single cluster.

We will need to estimate centroids and probabilities, which reduces to
calculating $\mathbb{E}\left[h(Y(X))\right]$ with
$h\colon\mathcal{Y}\mapsto\mathcal{H}$ for some set $\mathcal{H}$. To compute
the probability that $Y$ is in $A$, $P(Y(X)\in A)$, $h$ is the indicator
function $h\colon y\mapsto\mathds{1}_{y\in A}$ with $A\subset\mathcal{Y}$ and
$\mathcal{H}=\\{0,1\\}$. Alternatively $h$ is $h\colon y\mapsto
y\mathds{1}_{y\in A}$ with $\mathcal{H}=\mathcal{Y}$ to calculate the mean of
$Y$ when it belongs to $A$.

The importance sampling technique is based on the fact that the representation
of
$\mathbb{E}\left[h(Y(X))\right]$ as an expectation is not unique (L’Ecuyer et
al. 2009). Indeed, if we consider $\tilde{X}$ a random variable with density
function $g$ with a support such that
${\mathrm{supp}}(f_{X})\subset{\mathrm{supp}}(g)\subset\mathcal{X}$, then

$\mathbb{E}\left[h(Y(X))\right]=\mathbb{E}\left[h(Y(\tilde{X}))\frac{f_{X}(\tilde{X})}{g(\tilde{X})}\right]\leavevmode\nobreak\
.$

The importance sampling estimator of the mean comes from this last
representation:

$\hat{E}_{n}^{IS}(h)=\frac{1}{n}\sum^{n}_{k=1}h(Y(\tilde{X}^{k}))\frac{f_{X}(\tilde{X}^{k})}{g(\tilde{X}^{k})}\in\mathbb{R}^{p}\leavevmode\nobreak\
,$ (1)

with $(\tilde{X}^{k})^{n}_{k=1}$ i.i.d. random variables of density function
$g$. It is well known that such importance sampling estimator is unbiased and
consistent (Tabandeh et al. 2022) .

The covariance matrix of the estimator is

$\mathbb{V}\left(\hat{E}_{n}^{IS}(h)\right)=\frac{1}{n}\mathbb{V}\left(h(Y(\tilde{X}))\frac{f_{X}(\tilde{X})}{g(\tilde{X})}\right)\leavevmode\nobreak\
.$ (2)

The idea of the importance sampling method is to choose a relevant density
function $g$, i.e., a density function that minimizes the terms of the
diagonal of the above covariance matrix.

###  2.3 Importance sampling combined with quantization

In Lloyd’s algorithm, we need to repeatedly compute the expectations of $Y$
within the candidate Voronoi cells. At the $p^{th}$ iteration, these
expectations are,

$\gamma^{[p]}_{j}=\mathbb{E}\left[Y(X)\mid Y(X)\in
C^{\Gamma_{\ell}^{[p-1]}}_{j}\right]=\frac{\mathbb{E}\left[Y(X)\mathds{1}_{Y(X)\in
C^{\Gamma_{\ell}^{[p-1]}}_{j}}\right]}{\mathbb{P}(Y(X)\in
C^{\Gamma_{\ell}^{[p-1]}}_{j})}$ (3)

As for a given $\Gamma_{\ell}$, $\mathbb{P}(Y(X)\in
C^{\Gamma_{\ell}}_{j})=\mathbb{E}\left[\mathds{1}_{Y(X)\in
C^{\Gamma_{\ell}}_{j}}\right]$ is not known, we rely on the importance
sampling estimators for both the numerator and the denominator and we
introduce

$\left\\{\begin{array}[]{ll}\hat{P}_{n}(\Gamma_{\ell},j,Y)&=\frac{1}{n}\sum^{n}_{k=1}\mathds{1}_{Y(\tilde{X}^{k})\in
C^{\Gamma_{\ell}}_{j}}\frac{f_{X}(\tilde{X}^{k})}{g(\tilde{X}_{k})}\leavevmode\nobreak\
,\\\
\hat{E}_{n}(\Gamma_{\ell},j,Y)&=\frac{\frac{1}{n}\sum^{n}_{k=1}Y(\tilde{X}^{k})\mathds{1}_{Y(\tilde{X}^{k})\in
C^{\Gamma_{\ell}}_{j}}\frac{f_{X}(\tilde{X}^{k})}{g(\tilde{X}_{k})}}{\hat{P}_{n}(\Gamma_{\ell},j)}\leavevmode\nobreak\
,\end{array}\right.$ (4)

the estimators of $\mathbb{P}(Y(X)\in C^{\Gamma_{\ell}}_{j})$ and
$\mathbb{E}\left[Y(X)\mid Y(X)\in C^{\Gamma_{\ell}}_{j}\right]$ respectively.
The choice of the distribution $g$ is delicate and will be further discussed
in Section 6.

The overall approach was described so far in the very general context of a
random field $Y$. In the following, these estimators will be computed with $y$
the deterministic output of the hydraulic simulators, and $\hat{y}$ its
predictor.

##  3 Dealing with time-consuming experiments

In order to have an estimator of $\mathbb{E}\left[{y}(X)\mid{y}(X)\in
C^{\Gamma_{\ell}}_{j}\right]$ of sufficiently small variance for all cells
$j\in\\{1,\dots,\ell\\}$, we need to compute ${y}(\tilde{X}^{k})$ for a sample
$(\tilde{X}^{k})^{n}_{k=1}$ of large size $n$, with $y$ the deterministic
output of the hydraulic simulator. Realistic physical simulations would be too
costly for that purpose, and a metamodel is required to approximate ${y}(X)$.
Gaussian processes (GP) (Williams and Rasmussen 2006) are used here to replace
some of the simulations. GPs are readily applicable to low dimensional outputs
(most often scalars) but not to ${y}(X)$, the flood map, which is in
$\mathbb{R}^{4096}$ ($64\times 64$). To overcome that, Perrin et al. 2021
implemented a method based on functional principal component analysis (FPCA).
Their strategy, which we rely on, is to work with the projections of the maps
in a space induced by a basis of $n_{\mathrm{pc}}\ll 4096$ maps only. The
predictions are expressed as coordinates of a map in this basis, through a
Gaussian process approach. In the following, we call “metamodel“ the whole
procedure: the FPCA followed by a Gaussian process regression.

###  3.1 Projection on a basis of maps

The objective here is to find a basis of $n_{\mathrm{pc}}$ maps, with
$n_{\mathrm{pc}}$ small enough, and to compute the projection of the maps onto
this basis, thus representing every map $y(x)$ as a linear combination of the
$n_{\mathrm{pc}}$ basis maps:

$y^{\mathrm{proj}}(x):=t_{1}(x)y^{\mathrm{pca}}_{1}+t_{2}(x)y^{\mathrm{pca}}_{2}+\dots+t_{n_{\mathrm{pc}}}(x)y^{\mathrm{pca}}_{n_{\mathrm{pc}}}\leavevmode\nobreak\
.$

For this purpose, a functional principal components analysis (FPCA) is carried
out on a training database of maps
$({y}(x^{i}))_{i=1,\dots,n_{\mathrm{train}}}$. The steps are the following,
and further details are given in Appendix B.1. :

  * •

Find a decomposition of these maps, seen as realizations of a spatial
function, on a functional orthonormal basis
$\Phi=(\Phi_{1},\dots,\Phi_{K})^{\top}$ :
${y}(x^{i})=\Phi^{\top}\alpha(x^{i})$

  * •

Reduce the number of basis functions by keeping the $\tilde{K}$ most important
ones, $\tilde{\Phi}$, in terms of their contribution to the energy of the
coordinates in the database
$(\tilde{\alpha}(x^{i}))_{i=1,\dots,n_{\mathrm{train}}}$.

  * •

Perform a PCA on the selected coordinates to finally get the principal
components $(t(x^{i}))_{i=1,\dots,n_{\mathrm{train}}}$ with
$t(x^{i})=(t_{1}(x^{i}),\dots,t_{n_{\mathrm{pc}}}(x^{i}))^{\top}$ and the
$n_{\mathrm{pc}}\times\tilde{K}$ projection matrix $\Omega$ :
$\tilde{\alpha}(x^{i})\simeq\Omega^{\top}t(x^{i})$.

The obtained FPCA-basis maps are $\Omega\tilde{\Phi}$. The two basis maps
$y^{\mathrm{pca}}_{1}$ and $y^{\mathrm{pca}}_{2}$ of the coastal flooding
example are shown in Figure 3. Negative values can appear on these FPCA maps,
but positive depths are recovered with the linear combinations of the
principal components.

Figure 3: $y^{\mathrm{pca}}_{1}$ and $y^{\mathrm{pca}}_{2}$, the two FPCA maps
of the coastal flooding application without breach.

Figure 4 shows the map $y(x^{\mathrm{worst}})$ the most different from its
projection $y^{\mathrm{proj}}(x^{\mathrm{worst}})$, among the 1300 maps of the
training database $\mathcal{B}_{\mathrm{train}}$,

$y(x^{\mathrm{worst}})=\underset{y(x)\in\mathcal{B}_{\mathrm{train}}}{\operatorname*{arg\,max}}\lVert
y(x)-y^{\mathrm{proj}}(x)\lVert_{\mathcal{Y}}\leavevmode\nobreak\ .$

The norm $\lVert\leavevmode\nobreak\ \lVert_{\mathcal{Y}}$ selected here is
induced by the usual inner product $\langle
y,y^{\prime}\rangle_{\mathcal{Y}}=\sum_{i=1}^{m}y_{i}y^{\prime}_{i}$.

One can see that even in the worse case, the map $Y$ and its projection on the
FPCA basis $y^{\mathrm{proj}}$ are quite close. This highlights that FPCA
indeed captures most of the information. A map $y(x)$ is satisfactorily
compressed into a vector of 2 components, $\left(t_{1}(x),t_{2}(x)\right)$.

Figure 4: Most distant maps before and after FPCA: true
$y(x^{\mathrm{worst}})$ and its projection
$y^{\mathrm{proj}}(x^{\mathrm{worst}})$.

###  3.2 Gaussian process regression

We now aim at predicting the flood map for a new input $x^{\star}$ by
predicting its coordinates in the FPCA basis
$(t_{1}(x^{\star}),\dots,t_{n_{\mathrm{pc}}}(x^{\star}))$. For simplicity
these coordinates are assumed independent, but the study of their dependency
may be relevant in some cases (Gu and Shen 2020). Therefore, $n_{\mathrm{pc}}$
one-dimensional Gaussian processes define the metamodel for the flood maps.
Let $t\colon\mathcal{X}\subset\mathbb{R}^{d}\longrightarrow\mathbb{R}$ be a
deterministic function. $t$ is approximated as a realization of a Gaussian
process $Z(x)$ of prior mean $\mu(x)$ and prior covariance kernel
$k(x,x^{\prime})$. The kernel contains the spatial dependencies between $y(x)$
and $y(x^{\prime})$. In our study we use the stationary Matérn 5/2 kernel
(Palar et al. 2023).

We consider a design of experiments
$\mathbb{X}_{\mathrm{train}}=(x^{1},\dots,x^{n_{\mathrm{train}}})$ and the
observations $Z_{\mathrm{train}}=(z_{i})_{i\in
I_{\mathrm{train}}}=(t(x^{i}))_{i\in I_{\mathrm{train}}}$ with
$I_{\mathrm{train}}=\\{1,\dots,n_{\mathrm{train}}\\}$. $t$ can be approximated
using the distribution of $Z(x)\mid
Z(\mathbb{X}_{\mathrm{train}})=Z_{\mathrm{train}}$, which is a Gaussian
process with the following mean (Williams and Rasmussen 2006):

$\hat{t}(x)=\mu(x)+k(x,\mathbb{X}_{\mathrm{train}})k(\mathbb{X}_{\mathrm{train}},\mathbb{X}_{\mathrm{train}})^{-1}(Z_{\mathrm{train}}-\mu(\mathbb{X}_{\mathrm{train}}))$

where
$k(\mathbb{X}_{\mathrm{train}},\mathbb{X}_{\mathrm{train}})=[k(x^{i},x^{j})]_{i,j\in
I_{\mathrm{train}}}$ and $k(x,\mathbb{X}_{\mathrm{train}})=[k(x,x^{i})]_{i\in
I_{\mathrm{train}}}$. This method is applied to all the FPCA components
$i\in\\{1,\dots,n_{\mathrm{pc}}\\}$, for which the GP mean is a predictor
$\hat{t}_{i}(x^{\star})$ of $t_{i}(x^{\star})$. Then, $\hat{y}(x^{\star})$, a
predictor of $y(x^{\star})$ is obtained by inverse FPCA as explained in
Appendix B.2. The metamodel errors are investigated in the following sections.

##  4 Algorithm and performance metrics

Sections 2 and 3 described how to estimate the expectations in Lloyd’s
algorithm and how to predict a large number of maps that cannot be computed
with the time-consuming hydraulic simulators. This section details the overall
scheme for calculating the prototype maps and their probability masses, and
then explains how to evaluate its robustness.

###  4.1 Prototype Maps Algorithm

The Prototype Maps Algorithm incorporates, inside Lloyd’s algorithm, the
metamodeling of spatial output and the importance sampling scheme.

1:$({y}(x^{i}))_{i=1,\dots,n_{\mathrm{train}}}$, $f_{X}$, $g$, minDistance,
$\ell$

2:Sample $(\tilde{x}^{k})_{k\in\\{1,\dots,n_{\mathrm{maps}}\\}}$ i.i.d. of
density function $g$

3:Compute $(\hat{y}(\tilde{x}^{k}))_{1\leq k\leq n_{\mathrm{maps}}}$ from
$({y}(x^{i}))_{i=1,\dots,n_{\mathrm{train}}}$ (GP & FPCA)

4:Compute $(\frac{f_{X}(\tilde{x}^{k})}{g(\tilde{x}_{k})})_{1\leq k\leq
n_{\mathrm{maps}}}$

5:Initialize
$\Gamma_{\ell}^{[0]}\leftarrow\\{\gamma^{[0]}_{0},\dots,\gamma^{[0]}_{\ell}\\}\in\mathcal{Y}^{\ell}$

6:while
$\lVert\Gamma_{\ell}^{[k+1]}-\Gamma_{\ell}^{[k]}\rVert>\text{minDistance}$ do

7:
$\gamma^{[k+1]}_{j}\leftarrow\hat{E}_{n_{\mathrm{maps}}}(\Gamma_{\ell}^{[k]},j,\hat{y}),j=1,\dots,\ell$

8:$k\leftarrow k+1$

9:end while

10:$\hat{\Gamma}_{\ell}^{\star}=\Gamma_{\ell}^{[k]}$

11:Compute
$\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},j,\hat{y}),j=1,\dots,\ell$

Output: $\hat{\Gamma}_{\ell}^{\star}$ and
$\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},j,\hat{y})$

Algorithm 2 Prototype Maps Algorithm

The expression for
$\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},j,\hat{y})$ and
$\hat{E}_{n_{\mathrm{maps}}}(\Gamma_{\ell}^{[k]},j,\hat{y})$ stems from
applying Equations (4) to the Gaussian Process prediction of the maps. The
calculation of
$\hat{P}_{n_{\mathrm{maps}}}(\hat{\Gamma}_{\ell}^{\star},j,\hat{y})$ at step 8
could be extracted from the last iteration since it is the denominator of
$\hat{E}_{n_{\mathrm{maps}}}(\Gamma_{\ell}^{\star},j,\hat{y})$. However, these
quantities are the estimators of the probabilities that a flood map is part of
each Voronoi cell, for which we want to be as precise as possible, we thus re-
compute $\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},j,\hat{y})$ at the
end with $\tilde{n}>n_{\mathrm{maps}}$.

The overall complexity of the Prototype Maps Algorithm is
$\mathcal{O}((\tilde{n}+n_{\mathrm{maps}})n_{\mathrm{train}}^{2}n_{\mathrm{pc}}+\ell
s^{2}(n_{\mathrm{maps}}n_{\mathrm{it}}+\tilde{n}))$ with $n_{\mathrm{it}}$ the
number of iterations. The calculation of the complexity is detailed in Sire et
al. 2023a. The algorithm requires the storage of
$(\hat{y}(\tilde{X}^{k}))_{1\leq k\leq n_{\mathrm{maps}}}$ of size
$n_{\mathrm{maps}}s^{2}$ (with $s=64$ here), which is a limiting factor.

If the map simulations are not expensive, the algorithm can of course be
performed directly with the true maps,
$({y}(x^{i}))_{i=1,\dots,n_{\mathrm{maps}}}$, instead of the GP predictions.
minDistance $=10^{-16}$ is chosen here in the applications.

###  4.2 Performance metrics

We expect that Lloyd’s algorithm returns $\ell$ prototype maps associated to
the minimum quantization error, as well as $\ell$ probabilities of maps to
belong to each of the induced Voronoi cell. The errors will be evaluated both
on the centroids and on each cluster probability. To ensure the robustness of
the metrics, we consider
$(\Gamma_{\ell}^{r})_{r\in\\{1,\dots,n_{\Gamma_{\ell}}\\}}$ a family of
$n_{\Gamma_{\ell}}$ different random quantizations chosen as random
perturbations of a priori prototypes. Details about the random quantizations
are provided in Sire et al. 2023b. We also introduce $\Gamma_{\ell}^{\star}$
the optimal quantization obtained with the real maps coming from the
simulators, and $\hat{\Gamma}_{\ell}^{\star}$ the optimal quantization
obtained with the predicted maps.

#### Metamodel performance metrics.

The following performance metrics quantify the errors caused by the metamodel
(MM) on the quantization and on the cluster probability:

  * •

The excess in quantization error due to the metamodel,

$\epsilon_{\Gamma}^{MM}=\frac{\hat{e}(\hat{\Gamma}_{\ell}^{\star})-\hat{e}(\Gamma_{\ell}^{\star})}{\hat{e}(\Gamma_{\ell}^{\star})}\leavevmode\nobreak\
,$ (5)

where
$\hat{e}(\Gamma_{\ell})=\left(\frac{1}{n_{e}}\sum_{k=1}^{n_{e}}\lVert{y}(\tilde{X}^{k})-q_{\Gamma_{\ell}}({y}(\tilde{X}^{k}))\lVert_{\mathcal{Y}}^{2}\frac{f_{X}(\tilde{X}^{k})}{g(\tilde{X}^{k})}\right)^{\frac{1}{2}}$
is an empirical quantization error. We take $n_{e}=10^{6}$ here. This error
can only be computed in the analytical case (see Section 5.2) as
$\Gamma_{\ell}^{\star}$ can not be computed in the flooding case whose
computational cost is too high.

  * •

The relative probability error,

$\epsilon_{P}^{MM}(n,\Gamma_{\ell}^{r},j)=\frac{\mid\hat{P}_{n}(\Gamma_{\ell}^{r},j,{y})-\hat{P}_{n}(\Gamma_{\ell}^{r},j,\hat{y})\mid}{\hat{P}_{n}(\Gamma_{\ell}^{r},j,{y})}\leavevmode\nobreak\
.$ (6)

computed for all $r=1,\dots,n_{\Gamma_{\ell}}$ and $j=1,\dots,\ell$. Notice
that the prototypes are identical in both probabilities, and only the effect
of the metamodel on the cluster probability is assessed.

#### Importance sampling performance metrics.

The effect of the importance sampling (IS) on the quantization is quantified,
for a given Voronoi cell $C^{\Gamma_{\ell}^{r}}_{j}$, through two metrics:

  * •

The IS coefficient of variation of the membership probability,
$\epsilon_{P}^{IS}(\tilde{n},\Gamma_{\ell}^{r},j)$, which is an estimator of
$\frac{\sqrt{\mathbb{V}\left(\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{r},j,{y})\right)}}{\mathbb{E}\left(\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{r},j,{y})\right)}$.
It measures an error due to the importance sampling in the estimation of the
probability of the Voronoi cell number $j$. .

  * •

The IS centroid standard deviation,
$\epsilon_{\Gamma_{\ell}}^{IS}(n_{\mathrm{maps}},\Gamma_{\ell}^{r},j)$, which
is the $90\%$-quantile of the diagonal terms of the covariance matrix
$\mathbb{V}\left(\hat{E}_{n_{\mathrm{maps}}}(\Gamma_{\ell}^{r},j,{y})\right)$.
It characterizes the error coming from importance sampling in the estimation
of the centroid of the Voronoi cell number $j$. Among other things, it will
show potential variations due to the denominator being close to zero in the
centroid formula (3).

These quantities are calculated for $r=1,\dots,n_{\Gamma_{\ell}}$ and
$j=1,\dots,\ell$ and their distributions serve as the metrics for the
importance sampling performance. Details about the computation of these errors
are given in Sire et al. 2023b.

##  5 Application to an analytical case

###  5.1 Campbell test case

The performance of the method is first evaluated on an analytical test case,
called the Campbell2D function (Marrel et al. 2011). Contrarily to realistic
flooding simulators, a very high number of calls to the Campbell function are
possible, which makes it an appropriate case for testing the Prototype Maps
Algorithm. This function has eight inputs ($d=8$) and a spatial map as output,
i.e., its output is a function which depends on two inputs $(z=(z_{1},z_{2}))$
corresponding to spatial coordinates,

$\displaystyle h\colon\;\;\;\;\;\;\;\;\;\;[-1,5]^{8}$
$\displaystyle\longrightarrow\mathbb{L}^{2}([-90,90]^{2})$ $\displaystyle
x=(x_{1},\dots,x_{8})$ $\displaystyle\longmapsto h_{x}(z)\leavevmode\nobreak\
,$

where $z=(z_{1},z_{2})\in[-90,90]^{2}$, $x_{j}\in[-1,5]$ for $j=1,\dots,8$ and

$\displaystyle h_{x}(z_{1},z_{2})$
$\displaystyle=x_{1}\exp\left[-\frac{(0.8z_{1}+0.2z_{2}-10x_{2})^{2}}{60x_{1}^{2}}\right]$
$\displaystyle+(x_{2}+x_{4})\exp\left[\frac{(0.5z_{1}+0.5z_{2})x_{1}}{500}\right]+x_{5}(x_{3}-2)\exp\left[-\frac{(0.4z_{1}+0.6z_{2}-20x_{6})^{2}}{40x_{5}^{2}}\right]$
$\displaystyle+(x_{6}+x_{8})\exp\left[\frac{(0.3z_{1}+0.7z_{2})x_{7}}{250}\right]$

Examples of Campbell maps are shown in Figure 5.

(a)

(b)

Figure 5: Two examples of Campbell maps. 5(a) for
$x=(-0.8,4.7,1,-0.1,-0.4,3.2,0.7,-1)$ and 5(b) for
$x=(4.5,-0.3,3.5,1.4,0.3,2.6,-0.6,-1)$

As in the forthcoming flooding case, the output space is
$\mathcal{Y}=\mathcal{M}_{64,64}(\mathbb{R})$. $f_{X}$, the density function
of $X$, is the same as the density function of the input variables of the
flooding application. In order to have the same input domain as in the
forthcoming flooding example, we impose $x_{8}=-1$ and we introduce an affine
transformation, denoted $a$, from ${\mathrm{supp}}(f_{X})$ to $[-1,5]^{7}$, so
that the input space is
$\mathcal{X}={\mathrm{supp}}(f_{X})\subset\mathbb{R}^{7}$. Details are given
in Appendix C. In the end, the true map function is

$\displaystyle{y}\colon\mathcal{X}$ $\displaystyle\to\mathcal{Y}$
$\displaystyle x$
$\displaystyle\mapsto(h_{(a(x)_{1},\dots,a(x)_{7},-1)}(\mathbf{z}))_{z\in\mathcal{Z}}\leavevmode\nobreak\
,$

with
$\mathcal{Z}=\\{(-90+\frac{i}{64}180,-90+\frac{j}{64}180)\\}_{i,j\in\\{1,\dots,64\\}^{2}}$.

The sampling density, $g$, is the same as in the flooding case and is
described in Section 6.1. The number of prototypes $\ell=5$ is chosen a priori
here but its selection is discussed in Section 7.

###  5.2 Quantizing the Campbell function

#### Campbell metamodel errors.

As explained is Section 4.2, we gauge the performance of the empirical
quantization with the two following metrics.

First, the excess in quantization error due the metamodel,
$\epsilon_{\Gamma}^{MM}$ (Equation (5)), is equal to $3.6\%$ when evaluated
from $n_{\mathrm{maps}}=10^{6}$ maps.

Regarding the error on the probability estimators, we evaluate
$\epsilon_{P}^{MM}(10^{5},\Gamma_{\ell}^{r},j)$ for all the cells
$j=1,\dots,\ell$ of all the tested quantization $\Gamma_{\ell}^{r}$,
$r=1,\dots,n_{\Gamma_{\ell}}$. Figure 6(a) illustrates the distribution of
these errors. The results are really satisfying: the median and the 75th
percentile are below $1\%$. The tests on the Campbell function thus show that
the metamodel is sufficiently precise in terms of both the prototype maps and
the computed probability mass of the Voronoi cells.

(a) Relative probability error

(b) IS relative std. deviation

(c) IS centroid std. deviation

Figure 6:  Boxplots of the distributions of the various errors in the Campbell
test case. The distributions correspond to all $r=1,\ldots,100\text{ and
}j=1,\ldots,5.$ (a) $\epsilon_{P}^{MM}(10^{5},\Gamma_{\ell}^{r},j)$, the
relative probability error. (b)
$\epsilon_{P}^{IS}(10^{7},\Gamma_{\ell}^{r},j)$, the IS coefficient of
variation of the membership probability. (c)
$\epsilon_{\Gamma_{\ell}}^{IS}(10^{6},\Gamma_{\ell}^{r},j)$, the IS centroid
standard deviation.

#### Importance sampling errors.

As stated in Section 4.2, the errors related to the importance sampling scheme
are measured through the IS coefficient of variation of the membership
probability, $\epsilon_{P}^{IS}(\tilde{n},\Gamma_{\ell}^{r},j)$, and the IS
centroid standard deviation,
$\epsilon_{\Gamma_{\ell}}^{IS}(n_{\mathrm{maps}},\Gamma_{\ell}^{r},j)$, for
all choices of prototypes $r=1,\dots,n_{\Gamma_{\ell}}$ and all clusters
$j=1,\dots,\ell$.

The distributions of these metrics are presented in Figures 6(b) and 6(c) for
$\tilde{n}=10^{7}$ and $n_{\mathrm{maps}}=10^{6}$.
$\epsilon_{P}^{IS}(\tilde{n},\Gamma_{\ell}^{r},j)$ is around $1\%$, which is
of the same order of magnitude as the error from the metamodel, and which is
sufficiently accurate.
$\epsilon_{\Gamma_{\ell}}^{IS}(n_{\mathrm{maps}},\Gamma_{\ell}^{r},j)$ is
around 0.04, which is also very precise. These results validate the precision
of the importance sampling estimators.

#### Campbell prototype maps.

The outputs of the adapted Lloyd’s algorithm are the $\ell=5$ representative
maps $\Gamma_{\ell}^{\star}=\\{\gamma^{\star}_{1},\dots,\gamma^{\star}_{5}\\}$
drawn in Figure 7 which were obtained by using the true maps. The probability
mass of each Voronoi cell is given below each map with two estimates, one
calculated with the true maps and the other with the metamodeled maps. The
closeness of the results confirms the precision of the metamodel-based
approach. The same broad backslash-like pattern runs across the 5 maps, which
differ by their $h$ levels. It is explained by the density function $f_{X}$
illustrated in Appendix C, that concentrates $x_{2},x_{4}$ and $x_{5}$ in
specific zones that boost the occurence of this pattern.

(a) aa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,{y})=6.6\times 10^{-2}\\\ $
aaaaaa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,\hat{y})=6.2\times 10^{-2}$

(b) aa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,{y})=2.2\times 10^{-1}\\\ $
aaaaaa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,\hat{y})=2.2\times 10^{-1}$

(c) aa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,{y})=2.7\times 10^{-1}\\\ $
aaaaaa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,\hat{y})=2.7\times 10^{-1}$

(d) aa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,{y})=2.6\times 10^{-1}\\\ $
aaaaaa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,\hat{y})=2.6\times 10^{-1}$

(e) aa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,{y})=2.0\times 10^{-1}\\\ $
aaaaaa$\hat{P}_{\tilde{n}}(\Gamma_{\ell}^{\star},1,\hat{y})=2.0\times 10^{-1}$

Figure 7: The 5 prototype maps obtained from $n_{\mathrm{maps}}=10^{6}$ true
maps (without metamodel) in the Campbell case, and their probability masses
estimated with $\tilde{n}=10^{7}$ true maps (top of subcaptions) and estimated
with the predicted maps (bottom of the subcaptions).

##  6 Application to flooding

###  6.1 Description of the flooding test case

#### Parameterization.

The main application of this article is about the Boucholeurs district located
on the French Atlantic Coast (Figure 8(a)) near the city of La Rochelle. This
site was hit by the Xynthia storm in February 2010, which caused the
inundation of several areas and severe human and economic damage (Bertin et
al. 2012).

The main coastal flooding processes correspond to overflow and are simulated
with the numerical code MARS (Lazure and Dumas 2008), where adaptations were
made by the BRGM to take into account specificities of local water circulation
flows around connections (nozzles, spillways, etc. and breaching phenomena).
We focus on the interplay between tide and storm surge on the spatial
distribution of the maximum water depth after flooding in the urban area. As
shown in the Figure 8(b), the tide is simplified and assumed to be represented
by a sinusoidal signal, parameterized by the high-tide level $T$. The surge
signal is modeled as a triangular function defined by four parameters: the
peak amplitude $S$, the phase difference between the surge peak and the high
tide $t_{0}$, the time duration of the rising part $t_{-}$ and the time
duration of the falling part $t_{+}$. In total, 5 variables describe the
offshore conditions, $X=(T,S,t_{0},t_{-},t_{+})$.

Each of these input variables is assumed to be independent and described by an
empirically defined probability distribution based on the analysis performed
by Rohmer et al. 2018 of tide gauge measurements from Ile d’Aix 111data
available on data.shom.fr (the tide gauge closest to the study site) and La
Rochelle-La Pallice. A constant surge level offset (chosen at 65 cm i.e.,
compatible with future sea level rise in this region 222see e.g.
https://sealevel.nasa.gov/data_tools/17)) is added to its empirical
distribution to ensure that a sufficient number of flooding events are
induced. Appendix C shows the obtained probability distributions.

Regarding the variables representing a possible breach, two characteristics
are taken into account. First, the location of the breaches. Ten different
locations are selected. Six of them (B1-B6, Figure 8(a)) are related to
natural dunes and are selected based on observations made during the Xynthia
event (Muller et al. 2016). Four other locations (B7-B10) are related to
artificial dykes and are added because they are near vulnerable zones. The
other breach-related variable is the topographic level after failure,
parametrized as the fraction $p$ (erosion rate) of the initial crest level
left.

The distribution of the breach characteristics (location and erosion rate)
depends on the offshore conditions. If the maximum of the sum of the signals
(tide and surge) over time is higher than $70\%$ of the embankment’s height,
the probability of breaching is $0.5$. Vice versa, the probability of
breaching is assumed to be $10^{-4}$. Finally, we consider that if there is a
breach, then the location and the erosion rate follow a uniform distribution
$\mathcal{U}_{\\{1,\dots,10\\}}$ and $\mathcal{U}_{]0,1]}$, respectively. This
breach model is a compromise between simplicity and realism.

(a) Study site of Les Boucholeurs on the Atlantic coast in France, with the
locations of the breaches (triangle symbols) (Rohmer et al. 2018). The red
rectangle represents the investigated zone where quantization is performed.

(b) Parameterization of the tide and surge temporal signals (see Rohmer et al.
2018).

Figure 8:

#### Data set generation.

Our study relies on 1300 simulations of flood maps generated as follows:

  * •

The offshore conditions are sampled as the beginning of the Sobol Sequence in
5 dimensions (Cheng and Druzdzel 2013).

  * •

500 of the 1300 simulations are run without breach, and the breach parameters
of the 800 others are sampled uniformly according to
$\mathcal{U}_{\\{1,\dots,10\\}}$ for a single breach location and
$\mathcal{U}_{[0,1]}$ for the associated erosion rate.

These 1300 maps form the training database, and each map belongs to the space
$\mathcal{Y}=\mathcal{M}_{64,64}(\mathbb{R})$. The inner product
$\langle.,.\rangle_{\mathcal{Y}}$ of Theorem 2.1 is such that all
$\lambda_{i}$’s are set to 1 so that
$\forall(y,y^{\prime})\in\mathcal{Y}^{2},\langle
y,y^{\prime}\rangle_{\mathcal{Y}}=\sum_{i=1}^{64}\sum_{j=1}^{64}y_{i,j}y^{\prime}_{i,j}$.

#### IS sampling density.

The biased density function is $g=g_{1}\times\dots\times g_{7}$ with
$g_{i}=\frac{\mathds{1}_{I_{i}}}{\int_{I_{i}}dx},i=1,\dots,5$, where $I_{i}$
is the support of the marginal density function of the $i^{\mathrm{th}}$ input
variable (cf. Appendix C), $g_{6}=\frac{1}{10}\sum_{k=1}^{10}\delta_{k}$ and
$g_{7}=\frac{5}{13}\delta_{0}+(\frac{8}{13})\mathds{1}_{]0,1]}$.
The idea is to choose a single density for all the estimators related to all
the Voronoi cells.

The optimal biased density for which the IS estimators have zero variance is
theoretically known (L’Ecuyer et al. 2009) to be a normalized version of
$\mathds{1}_{{y}(\tilde{x})\in C^{\Gamma_{\ell}}_{j}}f_{X}(\tilde{x})$ or
$y(\tilde{x})\mathds{1}_{{y}(\tilde{x})\in
C^{\Gamma_{\ell}}_{j}}f_{X}(\tilde{x})$. It changes with each cell $j$,
affecting a non-zero probability mass only to the cell. It is not practical as
it requires the knowledge of the result of the estimation (probability or mean
map) for the normalization. The choice of a uniform distribution for the
variables related to the offshore conditions is an approximation to these
optimal densities that doesn’t depend on the cell. The uniform density is
relevant for quantization as it provides maps covering the whole domain
irrespective of their occurrence probability.

#### Dealing with the concentration of empty maps.

One of the main difference between the Campbell function and the flooding case
is a Dirac mass in the output distribution. Indeed, the probability that the
map is free of water ($0$ at every pixel of the matrix) is close to 1. As
Gaussian processes are not adapted to predict empty maps, their direct use in
the Prototype Maps Algorithm would overestimate the flooding probability. As a
counter-measure, a random forest classification step is added (Ayyadevara
2018), creating a hurdle model (McDowell 2003):

  * •

A first step predicts if the water level is below a threshold $t$ (it is a
parameter to tune) or not with the classification method.

  * •

If the water level is forecasted to be above $t$, then use the spatial
metamodel described in Section 3. Otherwise, the predicted map is empty.

A final post-processing step is to convert all negative water depth values
into null values.

Another difference with the Campbell case is the existence of two behaviours
with or without breach. Two metamodels are built independently for each of
these two situations.

#### Initializing Lloyd.

Another implementation aspect is the choice of the initial prototype maps at
the beginning of Lloyd’s algorithm (Step 4, Algorithm 2). Our a priori is that
the future Voronoi cells will be associated with diverse flooding levels. The
Prototype Maps Algorithm is first run with only the 1300 realistic maps
instead of the $n_{\mathrm{maps}}=10^{6}$ metamodeled maps. But this degraded
version of the algorithm, which provides the initial centroids, also needs to
be initialized, which is done with the maps associated to the $\ell$ equally
spaced quantiles of the training water volumes.

###  6.2 Quantizing floodings

#### Flooding metamodel error.

As presented in Section 4.2, to evaluate the impact of the metamodel on the
precision of the probability estimators, we compute the relative probability
error, $\epsilon_{P}^{MM}(1300,\Gamma_{\ell}^{r},j)$. Note that, unlike the
analytical Campbell function, the excess in quantization error,
$\epsilon_{\Gamma}^{MM}$, cannot be calculated since the true centroids are
not known.

The calculation of the relative probability error is composed of the following
steps:

  * •

The flood predictions, $(\hat{y}(\tilde{X}^{(j)}))_{1\leq j\leq 1300}$, are
computed from the $({y}(\tilde{X}^{(j)}))_{1\leq j\leq 1300}$ hydraulic
simulations within a 10-fold cross-validation procedure. Each of the fold
involves a classification and a metamodeling from 1170 maps followed by a
prediction of 130 maps.

  * •

The probabilities $\hat{P}_{1300}(\Gamma_{\ell},j,{y})$ and
$\hat{P}_{1300}(\Gamma_{\ell},j,\hat{y})$, $j\in\\{1,\dots,5\\}$, are computed
for the perturbed centroids (generated as explained in Section 4.2)
$\Gamma_{\ell}\in(\Gamma_{\ell}^{r})_{r\in\\{1,\dots,n_{\Gamma_{\ell}}\\}}$.

  * •

The relative probability error $\epsilon_{P}^{MM}(1300,\Gamma_{\ell},j)$ is
computed following Equation (6) for
$\Gamma_{\ell}\in(\Gamma_{\ell}^{r})_{r\in\\{1,\dots,n_{\Gamma_{\ell}}\\}}$
and $j=1,\dots,5$.

The distribution of $\epsilon_{P}^{MM}(1300,\Gamma_{\ell}^{r},j)$, the
relative errors, is shown in Figure 9(a) for $n_{\Gamma_{\ell}}=100$. The
median of the error is around $2.5\%$. It reveals that the FPCA and the
Gaussian process together are effective in predicting the Voronoi cell in
which a flood map is located, i.e., in predicting the type of flooding of a
new input vector. This performance is obtained with only two FPCA basis maps
that explain $98\%$ of the variance, see Figure 3, which highlights the low
intrinsic dimension of the phenomenon.

(a) Relative probability error

(b) IS relative std. deviation

(c) IS centroid std. deviation

Figure 9:  Boxplots of the distributions of the various errors in the flooding
case. The distributions correspond to all $r=1,\ldots,100\text{ and
}j=1,\ldots,5.$ (a) $\epsilon_{P}^{MM}(1300,\Gamma_{\ell}^{r},j)$, the
relative probability error. (b)
$\epsilon_{P}^{IS}(10^{7},\Gamma_{\ell}^{r},j)$, the IS coefficient of
variation of membership probability. (c)
$\epsilon_{\Gamma_{\ell}}^{IS}(10^{6},\Gamma_{\ell}^{r},j)$, the IS centroid
standard deviation.

#### Flooding importance sampling error.

We complete this study by an analysis of the importance sampling error as
explained in Sire et al. 2023b. A bootstrap approach provides, for each of the
$j=1,\ldots,\ell$ centroids of each $\Gamma_{\ell}^{r}$, the IS coefficient of
variation of the membership probability,
$\epsilon_{P}^{IS}(\tilde{n},\Gamma_{\ell}^{r},j)$, and the IS centroid
standard deviation,
$\epsilon_{\Gamma_{\ell}}^{IS}(n_{\mathrm{maps}},\Gamma_{\ell}^{r},j)$. The
distributions of these metrics are given in Figure 9(b) and 9(c) for
$\tilde{n}=10^{7}$ and $n_{\mathrm{maps}}=10^{6}$.
$\epsilon_{P}^{IS}(\tilde{n},\Gamma_{\ell}^{r},j)$ is lower than $1\%$,
showing that the error on the probability is mainly driven by the error of the
metamodel.
$\epsilon_{\Gamma_{\ell}}^{IS}(n_{\mathrm{maps}},\Gamma_{\ell}^{r},j)$ is less
than $0.02$ meters for $n_{\mathrm{maps}}=10^{6}$ which is quite satisfying in
a flooding context.

#### Flooding prototype maps.

The application of the Prototype Maps Algorithm to Les Boucholeurs with
$n_{\mathrm{maps}}=10^{6}$ gives the $5$ representative maps of Figure 10.
Each of them represents a type of flooding. The probability masses are
computed with $\tilde{n}=10^{7}$ to increase precision, which was not possible
in the Lloyd’s algorithm because of the computational cost of the multiple
iterations. The most serious the flooding, the lowest the probability, except
for the Voronoi cell number 3 (Figure 10(c)) which is a bit less likely than
Voronoi cell number 4 (Figure 10(d)), even if it is associated to more
moderate floodings. Although this seems counter-intuitive, it is possible that
a specific type of small flooding is less likely than another type of more
severe flooding. Among the Voronoi cells associated to the strongest
floodings, the prototype number 4 (Figure 10(d)) is a map with shallow water
covering a large area, while the flooded area is smaller for cell number 5
(Figure 10(e)) but with a greater water depth near the coastline. These 5
Voronoi cells can be visualized on the two first FPCA axes that summarize 98%
of the inertia of the 1300 flooding maps. Figure 11 illustrates this (and
Appendix D discusses these results more thoroughly), showing that the fourth
and fifth cells associated to the most extreme floodings contains very
different types of maps.

(a) $\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},1,\hat{y})=9.7\times
10^{-1}$, 1 in 1.03

(b) $\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},2,\hat{y})=2.0\times
10^{-2}$, 1 in 50

(c) $\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},3,\hat{y})=5.2\times
10^{-3}$, 1 in 192

(d) $\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},4,\hat{y})=6.2\times
10^{-3}$, 1 in 161

(e) $\hat{P}_{\tilde{n}}(\hat{\Gamma}_{\ell}^{\star},5,\hat{y})=1.4\times
10^{-3}$, 1 in 714

Figure 10: Centroids calculated by the Prototype Maps Algorithm with
$n_{\mathrm{maps}}=10^{6}$ predicted maps, and their probability mass
estimated with $\tilde{n}=10^{7}$. The frequency, i.e., the inverse of the
probability is provided for an easier interpretation. The water depths (in
meters) are drawn with colors ranging from light to dark blue. Blue contour
lines represent a water depth of 15 cm.

##  7 Summary and perspectives

This work proposes a new visualization method based on quantization adapted to
rare events and for which a limited database is available, as it occurs with
computationally costly simulations. Without loss of generality, the method,
called Prototype Maps Algorithm, is developed around an application to flood
maps. The visualization takes the form of a set of prototype maps and their
probability masses that best represent the continuous probability distribution
on the space of maps. Importance sampling and a metamodel scheme based on
Functional PCA combined with Gaussian process regression are used to address
the challenges of the scarcity of the (flooding) events and the small data set
size. The obtained prototype maps along with their probabilities are an
intuitive representation of the flooding risk while remaining meaningful
mathematically. Indeed, the prototype maps are a discrete approximation of the
continuous random flooding event. The distance between the discrete and the
continuous distribution is the quantization error, which is easily
interpretable as it is the expected squared distance between a random map and
its closest prototype map.

Performance metrics used in the article validate the precision of our method
on an analytical test problem and a real flooding case.

The following developments are worthwhile, the first being related to the
tuning of the algorithm, the others dealing with complementary approaches.

#### Tuning of the algorithm.

First, a more sophisticated biased density function of importance sampling
could be sought. We have integrated importance sampling inside Lloyd’s
algorithm to deal with the low flooding probability when estimating the
centroid and the probabilities of each cluster. The choice of a uniform
density $g$ proposed here for the offshore conditions is a compromise between
simplicity and the exploration of the whole domain and was empirically seen to
keep the variance of the estimators small. Yet a real optimization of this
density has not been performed. While it is a difficult problem with
functional unknowns, the search for a better biased density, defined
independently for each cell, is a relevant continuation to this work.

Second, the number of prototype maps, $\ell$, is chosen a priori to remain
usable in practice while capturing sufficiently diverse floodings. Strategies
to identify the number of Voronoi cells could be implemented based, for
instance, on the evolution of the quantization error with the number of
representatives.

#### Complementary approaches.

The method described here characterizes representative floodings. It would now
be valuable to know more about the forcing and breaching conditions leading to
these scenarios, i.e., to study how the inputs are related to each prototype.

An inherent feature of the obtained quantization is that none of the prototype
maps is the no-flood map. The centroid of the first Voronoi cell has a
negligible water depth ($3\times 10^{-5}$ meters on average), but it is not
really empty. However, fortunately, it is known that the probability that a
flood map is empty is close to 1. The formulation of quantization could then
be adapted so as to force the no-flood map (or more generally a Dirac of a
continuous probability law) to be one of the representatives.

In this article, the metamodel is the kriging mean. The uncertainties of the
Gaussian process were not considered. Another future development could account
for this metamodel uncertainty in the computation of the conditional
expectations and the probabilities.

## Acknowledgement

This research was conducted with the support of IRSN and BRGM, through the
consortium in Applied Mathematics CIROQUO
(https://doi.org/10.5281/zenodo.6581217), gathering partners in technological
and academia in the development of advanced methods for Computer Experiments.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

SUPPLEMENTARY MATERIAL

Supplementary material description:

Note describing the git repository (readme.md).

Codes related to the Campbell2D case:

Git repository containing R notebooks to reproduce all the experiments related
to the Campbell2D function that are described in the article.
(https://github.com/charliesire/quantization_Campbell2D.git)

R package FunQuant:

R package containing functions to perform the Prototype Maps Algorithm and
compute the associated performance metrics.
(https://github.com/charliesire/FunQuant.git)

Importance sampling performance metrics:

Note providing details about the computation of the importance sampling
performance metrics.(is_perf_metrics.pdf).

Complexity of the Prototype Maps Algorithm:

Note explaining the computation of the complexity of the Prototype Maps
Algorithm (complexity_pma.pdf).

## References

  * Ayyadevara (2018) Ayyadevara, V. K. (2018).  Random Forest, pp. 105–116.  Berkeley, CA: Apress.
  * Bertin et al. (2012) Bertin, X., N. Bruneau, J. Breilh, and K. M. Fortunato, A. (2012).  Importance of wave age and resonance in storm surges: The case xynthia, bay of biscay.  Ocean Model 42, 16–30.
  * Beven et al. (2015) Beven, K., R. Lamb, D. Leedal, and N. Hunter (2015).  Communicating uncertainty in flood inundation mapping: a case study.  International Journal of River Basin Management 13(3), 285–295.
  * Bock (2008) Bock, H.-H. (2008, 01).  Origins and extensions of the k-means algorithm in cluster analysis.  Journal Électronique d’Histoire des Probabilités et de la Statistique [electronic only] 4, Article 14.
  * Cheng and Druzdzel (2013) Cheng, J. and M. J. Druzdzel (2013).  Computational investigation of low-discrepancy sequences in simulation algorithms for bayesian networks.  CoRR abs/1301.3841, 72–91.
  * Daubechies (1988) Daubechies, I. (1988).  Orthonormal bases of compactly supported wavelets.  Communications on Pure and Applied Mathematics 41, 909–996.
  * Durrett (2010) Durrett, R. (2010).  Probability: Theory and Examples.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Gray and Neuhoff (1998) Gray, R. and D. Neuhoff (1998).  Quantization.  IEEE Transactions on Information Theory 44(6), 2325–2383.
  * Gu and Shen (2020) Gu, M. and W. Shen (2020).  Generalized probabilistic principal component analysis of correlated data.  Journal of Machine Learning Research 21(13), 1–41.
  * Hastie et al. (2009) Hastie, T., R. Tibshirani, and J. Friedman (2009).  The elements of statistical learning: data mining, inference and prediction (2 ed.).  Springer.
  * Lazure and Dumas (2008) Lazure, P. and F. Dumas (2008, 02).  An external-internal mode coupling for a 3d hydrodynamical model for applications at regional scale (mars).  Advances in Water Resources 31, 233–250.
  * L’Ecuyer et al. (2009) L’Ecuyer, P., M. Mandjes, and B. Tuffin (2009, 03).  Rare Event Simulation using Monte Carlo Methods, pp. 17 – 38\.  John Wiley and Sons.
  * MacQueen (1967) MacQueen, J. B. (1967).  Some methods for classification and analysis of multivariate observations.  In L. M. L. Cam and J. Neyman (Eds.), Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1, pp. 281–297. University of California Press.
  * Mahajan et al. (2012) Mahajan, M., P. Nimbhorkar, and K. Varadarajan (2012).  The planar k-means problem is np-hard.  Theoretical Computer Science 442, 13–21.  Special Issue on the Workshop on Algorithms and Computation (WALCOM 2009).
  * Marrel et al. (2011) Marrel, A., B. Iooss, M. Jullien, B. Laurent, and E. Volkova (2011).  Global sensitivity analysis for models with spatially dependent outputs.  Environmetrics 22(3), 383–397.
  * McDowell (2003) McDowell, A. (2003).  From the help desk: Hurdle models.  The Stata Journal 3(2), 178–184.
  * Muller et al. (2016) Muller, H., A. Van Rooijen, D. Idier, and R. J. Pedreros, R. (2016).  Assessing storm impact on a french coastal dune system using morphodynamic modeling.  J.Coast Res. 33(2), 254–272.
  * Pagès and Yu (2013) Pagès, G. and J. Yu (2013, 12).  Pointwise convergence of the Lloyd algorithm in higher dimension.  working paper or preprint.
  * Pagès (2014) Pagès, G. (2014, 07).  Introduction to optimal vector quantization and its applications for numerics.  Technical report, LPMA.  54 pages.
  * Palar et al. (2023) Palar, P. S., L. Parussini, L. Bregant, K. Shimoyama, and L. R. Zuhal (2023).  On kernel functions for bi-fidelity gaussian process regressions.  Structural and Multidisciplinary Optimization 66(2), 37\.
  * Perrin et al. (2021) Perrin, T., O. Roustant, J. Rohmer, O. Alata, J. Naulin, D. Idier, R. Pedreros, D. Moncoulon, and P. Tinard (2021, 07).  Functional principal component analysis for global sensitivity analysis of model with spatial output.  Reliability Engineering and System Safety 211, 107522\.
  * Rohmer et al. (2018) Rohmer, J., D. Idier, F. Paris, F. Pedreros, and L. J. (2018).  Casting light on forcing and breaching scenarios that lead to marine inundation: Combining numerical simulations with a random-forest classification approach.  Environmental Modelling & Software 11, 64–80.
  * Santner et al. (2003) Santner, T., B. Williams, and W. Notz (2003, 01).  The Design and Analysis Computer Experiments.  Springer.
  * Seipel and Lim (2017) Seipel, S. and N. J. Lim (2017).  Color map design for visualization in flood risk assessment.  International Journal of Geographical Information Science 31(11), 2286–2309.
  * Sire et al. (2023a) Sire, C., R. Le Riche, D. Rullière, J. Rohmer, L. Pheulpin, and Y. Richet (2023a).  Complexity of the prototype maps algorithm.  https://github.com/charliesire/quantization_Campbell2D.git.
  * Sire et al. (2023b) Sire, C., R. Le Riche, D. Rullière, J. Rohmer, L. Pheulpin, and Y. Richet (2023b).  Importance sampling performance metrics.  https://github.com/charliesire/quantization_Campbell2D.git.
  * Tabandeh et al. (2022) Tabandeh, A., G. Jia, and P. Gardoni (2022).  A review and assessment of importance sampling methods for reliability analysis.  Structural Safety 97, 102216.
  * Williams and Rasmussen (2006) Williams, C. K. and C. E. Rasmussen (2006).  Gaussian processes for machine learning, Volume 2.  MIT press Cambridge, MA.

##  Appendix A Proof of theorem 2.1

In this Appendix, we generalize the proof of optimal quantization of Kieffer
and Cuesta-Albertos (Pagès and Yu 2013 ) to random fields $Y$. Let us consider
a probability space $(\Omega_{0},\mathcal{F}_{0},\mathbb{P})$, a random
variable $Y\colon\Omega_{0}\to\mathcal{Y}$ and the optimal quantization of
$\ell$ prototype maps
$\Gamma_{\ell}^{\star}=\\{\gamma^{\star}_{1},\dots,\gamma^{\star}_{\ell}\\}\in\underset{\Gamma_{\ell}\in\mathcal{Y}^{\ell}}{\operatorname*{arg\,min}}\>(e(\Gamma_{\ell}))$.

We suppose that

  * •

$\mathcal{Y}$ is of finite dimension $m$ ,

  * •

$\forall z,y\in\mathcal{Y},\langle
z,y\rangle_{\mathcal{Y}}=\sum_{i=1}^{m}\lambda_{i}z_{i}y_{i}$ with $\forall
i,\lambda_{i}>0$ ,

  * •

and $Y\in
L^{2}_{\mathcal{Y}}(\mathcal{F}_{0})=\\{W\colon\Omega_{0}\to\mathcal{Y}\in\mathcal{F}_{0},\mathbb{E}\left[\lVert
W\lVert_{\mathcal{Y}}^{2}\right]<+\infty\\}$ .

Let us show that $\forall i\in\\{1\dots\ell\\},\quad\mathbb{E}\left[Y\mid Y\in
C^{\Gamma_{\ell}^{\star}}_{i}\right]=\gamma^{\star}_{i}$.

We introduce the inner product of the space
$L^{2}_{\mathcal{Y}}(\mathcal{F}_{0})$: $\langle
W,Y\rangle=\mathbb{E}\left[\langle W,Y\rangle_{\mathcal{Y}}\right],\>W,Y\in
L^{2}_{\mathcal{Y}}(\mathcal{F}_{0})$. We set
$\mathcal{F}=\sigma(q_{\Gamma_{\ell}}(Y))\subset\mathcal{F}_{0}$ the
$\sigma$-algebra induced by $q_{\Gamma_{\ell}}(Y)$.

######  Proposition A.1.

Let $W\in L^{2}_{\mathcal{Y}}(\mathcal{F}_{0})$,
$\pi_{L^{2}(\mathcal{F})}(W)=\mathbb{E}\left[W\mid\mathcal{F}\right]$, with
$\pi_{L^{2}(\mathcal{F})}(W)$ the orthogonal projection of $W$ on
$L^{2}(\mathcal{F})$ (Durrett 2010, pp. 229)

###### Proof.

Let $W\in L^{2}_{\mathcal{Y}}(\mathcal{F}_{0}).\leavevmode\nobreak\ \forall
Z\in L^{2}_{\mathcal{Y}}(\mathcal{F}),$

$\begin{split}\mathbb{E}\left[\langle
Z,W\rangle_{\mathcal{Y}}\right]&=\mathbb{E}\left[\mathbb{E}\left[\langle
Z,W\rangle_{\mathcal{Y}}\mid\mathcal{F}\right]\right]\\\
&=\mathbb{E}\left[\mathbb{E}\left[\sum_{i=1}^{n}\lambda_{i}z_{i}w_{i}\mid\mathcal{F}\right]\right]\\\
&=\mathbb{E}\left[\sum_{i=1}^{n}\lambda_{i}\mathbb{E}\left[z_{i}w_{i}\mid\mathcal{F}\right]\right]\\\
&=\mathbb{E}\left[\sum_{i=1}^{n}\lambda_{i}z_{i}\mathbb{E}\left[w_{i}\mid\mathcal{F}\right]\right]\>\mbox{because
}\forall i,\>z_{i}\in\mathcal{F}\\\ &=\mathbb{E}\left[\langle
Z,\mathbb{E}\left[W\mid\mathcal{F}\right]\rangle_{\mathcal{Y}}\right]\leavevmode\nobreak\
.\\\ \end{split}$

∎

We have $q_{{\Gamma_{\ell}}^{\star}}(Y)-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]\in\mathcal{F}=\sigma(q_{{\Gamma_{\ell}}^{\star}}(Y))$,
so $Y-\mathbb{E}\left[Y\mid q_{{\Gamma_{\ell}}^{\star}}(Y)\right]$ is
orthogonal to $q_{{\Gamma_{\ell}}^{\star}}(Y)-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]$ for the inner product of the space
$L^{2}_{\mathcal{Y}}(\mathcal{F}_{0})$ according to proposition A.1. Then, the
Pythagorean theorem gives

$\displaystyle\mathbb{E}\left[\lVert
Y-q_{{\Gamma_{\ell}}^{\star}}(Y)\lVert_{\mathcal{Y}}^{2}\right]=$
$\displaystyle\mathbb{E}\left[\lVert Y-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]\lVert_{\mathcal{Y}}^{2}\right]+$ (7)
$\displaystyle\mathbb{E}\left[\lVert
q_{{\Gamma_{\ell}}^{\star}}(Y)-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]\lVert_{\mathcal{Y}}^{2}\right]\leavevmode\nobreak\
.$

We set
$\tilde{\Gamma_{\ell}}=\\{\tilde{\gamma}_{1},\dots,\tilde{\gamma}_{\ell}\\}$
such that $\forall i\in\\{1\dots N\\},\tilde{\gamma}_{i}=\mathbb{E}\left[Y\mid
Y\in C^{\Gamma_{\ell}^{\star}}_{i}\right]$.

  * •

$\mathbb{E}\left[\lVert
Y-q_{\tilde{\Gamma_{\ell}}}(Y)\lVert_{\mathcal{Y}}^{2}\right]\leq\mathbb{E}\left[\lVert
Y-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]\lVert_{\mathcal{Y}}^{2}\right]$ :

$\displaystyle\forall y\in\mathcal{Y},\exists i\in\\{1,\dots,N\\},\lVert
y-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)=q_{{\Gamma_{\ell}}^{\star}}(y)\right]\lVert_{\mathcal{Y}}$
$\displaystyle=\lVert y-\mathbb{E}\left[Y\mid Y\in
C^{\Gamma_{\ell}^{\star}}_{i}\right]\lVert_{\mathcal{Y}}$
$\displaystyle=\lVert y-\tilde{\gamma}_{i}\lVert_{\mathcal{Y}}$
$\displaystyle\geq\lVert y-q_{\tilde{\Gamma_{\ell}}}(y)\lVert_{\mathcal{Y}}$
$\displaystyle\text{ by definition of }q_{\tilde{\Gamma_{\ell}}}.$

  * •

$\mathbb{E}\left[\lVert
Y-q_{\tilde{\Gamma_{\ell}}}(Y)\lVert_{\mathcal{Y}}^{2}\right]\geq\mathbb{E}\left[\lVert
Y-q_{{\Gamma_{\ell}}^{\star}}(Y)\lVert_{\mathcal{Y}}^{2}\right]$ by definition
of $\Gamma_{\ell}^{\star}$

Then, $\mathbb{E}\left[\lVert
Y-q_{{\Gamma_{\ell}}^{\star}}(Y)\lVert_{\mathcal{Y}}^{2}\right]\leq\mathbb{E}\left[\lVert
Y-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]\lVert_{\mathcal{Y}}^{2}\right]$.
We conclude from Equation (7): $\mathbb{E}\left[\lVert
q_{{\Gamma_{\ell}}^{\star}}(Y)-\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]\lVert_{\mathcal{Y}}^{2}\right]=0$ so
that $q_{{\Gamma_{\ell}}^{\star}}(Y)=\mathbb{E}\left[Y\mid
q_{{\Gamma_{\ell}}^{\star}}(Y)\right]$.
This is a known result (Pagès 2014) that, therefore, applies to our study with
the random field $Y(X)$.

##  Appendix B Projection on a basis of maps

###  B.1 Functional Principal Components Analysis

Every flood map $y(x)$ that we want to model can be considered as a matrix
indicating the water depth at each pixel, i.e., as a function

$\displaystyle y_{x}\colon\mathbb{R}^{2}$
$\displaystyle\longrightarrow\mathbb{R}$ $\displaystyle z=(z_{1},z_{2})$
$\displaystyle\longmapsto y_{x}(z).$

We work with a training database of $n_{\mathrm{train}}$ maps
$(y_{x^{i}}(z))_{i=1,\dots,n_{\mathrm{train}}}$ and the objective here is to
compute the projection of the flood maps in a FPCA basis of $n_{\mathrm{pc}}$
maps. We first consider a functional orthonormal basis
$\Phi(z)=(\Phi_{1}({z}),\dots,\Phi_{K}({z}))^{\top}$. D4 Daubechies wavelets
(Daubechies 1988), which are widely used in image processing, are chosen here,
as Perrin et al. 2021 shows their accuracy for a similar problem. The
procedure consists in the following steps:

  1. 1.

Find the coordinates $(\alpha(x^{i}))_{i=1,\dots,n_{\mathrm{train}}}$ of the
training maps $(y_{x^{i}}(z))_{i=1,\dots,n_{\mathrm{train}}}$ in the $\Phi$
basis:
$y_{x^{i}}(z)=\sum_{k=1}^{K}\alpha_{k}(x^{i})\Phi_{k}(z)=\Phi(z)^{\top}\alpha(x^{i})$

  2. 2.

Compute
$\lambda_{k}=\mathbb{E}\left[\frac{\alpha_{k}(X)^{2}}{\sum_{j=1}^{K}\alpha_{j}(X)^{2}}\right],\>k=1\dots
K$, order them by decreasing values and select the $\tilde{K}\ll K$ most
important vectors of the basis (the $\tilde{K}$ vectors associated with the
highest $\lambda_{k}$) such that $\sum_{k=1}^{\tilde{K}}\lambda_{k}\geq p$,
with $p$ a given proportion of energy. Then, we have
$y_{x^{i}}(z)\simeq\sum_{k=1}^{\tilde{K}}\alpha_{k}(x^{i})\Phi_{k}(z)=\tilde{\Phi}(z)^{\top}\tilde{\alpha}(x^{i})$

  3. 3.

Apply PCA in $\mathbb{R}^{\tilde{K}}$ on the dataset of coefficients evaluated
at the design points $(\tilde{\alpha}(x^{i}))_{i=1,\dots,n_{\mathrm{train}}}$.
Denote
$t_{1}(x^{i}),\dots,t_{n_{\mathrm{pc}}}(x^{i}),i=1,\dots,n_{\mathrm{train}}$,
the coordinates of the first $n_{\mathrm{pc}}$ principal components and
$\omega_{1},\dots,\omega_{n_{\mathrm{pc}}}$ the associated eigenvectors in
$\mathbb{R}^{\tilde{K}}.$ We have
$\tilde{\alpha}(x^{i})\simeq\sum_{j=1}^{n_{\mathrm{pc}}}t_{j}(x^{i})w_{j}=\Omega^{\top}t(x^{i})$,
with $\Omega$ the $n_{\mathrm{pc}}\times\tilde{K}$ projection matrix.

Finally we can write
$y_{x^{i}}(z)\simeq(\Omega\tilde{\Phi}(z))^{\top}t(x^{i})$, with
$\Omega\tilde{\Phi}$ a basis of $n_{\mathrm{pc}}$ maps. $\tilde{K}$ (or,
equivalently, $p$) and $n_{\mathrm{pc}}$ are tuned by cross-validation.

###  B.2 Inverse FPCA

As described in Section 3, Gaussian processes are used to predict
$\hat{t}(x^{\star})$ for a new input $x^{\star}$. The following steps provide
the map $\hat{y}(x^{\star})$ from the coordinates $\hat{t}(x^{\star})$:

  1. 1.

Compute the coefficients of the vector $\alpha(x^{\star})$

     * •

for
$k=\tilde{K}+1,\dots,K,\hat{\alpha}_{k}(x^{\star})=\frac{1}{n_{\mathrm{train}}}\sum_{i=1}^{n_{\mathrm{train}}}\alpha_{k}(x^{i})$

     * •

for
$k=1,\dots,\tilde{K},\hat{\alpha}_{k}(x^{\star})=\sum_{j=1}^{n_{\mathrm{pc}}}\hat{t}_{j}(x^{\star})\omega_{j}$

  2. 2.

Compute $(\hat{y}_{x^{\star}}(z))$ with
$\hat{y}_{x^{\star}}(z)=\sum_{k=1}^{K}\hat{\alpha}_{k}(x^{\star})\Phi_{k}(z)$

The computation of $y_{x^{\star}}(z)$ at the targeted family of pixels
provides the predicted map $\hat{y}(x^{\star})$.

##  Appendix C Probability Distributions

The density functions of the offshore variables are obtained through empirical
observations as described in Section 6.1. They are shown in Table 1, along
with their support. The location of the breach is an integer in
$\\{1,\dots,10\\}$, and the erosion rate is in $[0,1]$.

Variable $X_{i}$ |  Distribution $I_{i}$ |  Support |  Law type
---|---|---|---
$T$ a Tide |  a
|  [0.52,3.59] |  Empirical
$S$ a Surge |  a
|  [0.65,2.5] |  Combination of a genera- lised Pareto distribution above a threshold of 0.55m and the empirical distribu- tion derived from the tide gauge measurements
$t_{0}$ a Phase difference |  a
|  [-8.05,8.05] |  Empirical
$t_{-}$ a Duration rising part |  a
|  [-12,0] |  Empirical
$t_{+}$ a Duration falling part |  a
|  [0,12.2] |  Empirical
Table 1: Variables describing the offshore conditions and their probability
distribution.

The affine transformation $a$ from ${\mathrm{supp}}(f_{X})$ to $[-1,5]^{7}$,
introduced in the Campbell case, is then

$\displaystyle a\colon{\mathrm{supp}}(f_{X})$ $\displaystyle\to[-1,5]^{7}$
$\displaystyle x$
$\displaystyle\mapsto\begin{bmatrix}-1+6\frac{x-0.52}{3.59-0.52}\\\
-1+6\frac{x-0.65}{2.5-0.65}\\\ -1+6\frac{x-8.05}{8.05+8.05}\\\
-1+6\frac{x+12}{12}\\\ -1+6\frac{x-12.2}{12.2}\\\ -1+6\frac{x-1}{9}\\\
-1+6x\end{bmatrix}\leavevmode\nobreak\ .$

##  Appendix D Voronoi cells on the FPCA axes

In the flooding case, the two first FPCA axes summarize $98\%$ of the inertia
of the training 1300 maps. It is relevant to visualize the quantization cells
through their two first FPCA coefficients. Among other things, this allows to
have a better understanding of their size. Figure 11(a) shows the distribution
of $10^{5}$ maps sampled with the uniform density $g$ and the prototype maps
in the FPCA basis calculated for the maps without breach. As expected, cell 1
is almost a Dirac at the empty map whose coefficients are $(-27,8)$. Moreover,
the largest cells are the ones associated to the most extreme floodings. For
instance, cell 5 contains very diversified maps in terms of water depth
magnitude and spatial distribution. It is an expected property of the optimal
Voronoi cells: the cells with high probabilities are small and vice versa.
However, this characteristic may not be true is some specific contexts, as
shown in Appendix E. Figure 11(b) shows the probability density of the maps
expressed in the 2 first FPCA coefficients. It is a histogram where each bin
has a colour corresponding to the weight
$\log\left(\sum_{k=1}^{n_{\text{bin}}}\frac{f(\tilde{X}^{k})}{g(\tilde{X}^{k})}\right)$
associated to the IS density ratio of the maps, $n_{\text{bin}}$ being the
number of maps falling within each bin. It is essential to understand the
position of the prototype maps and the probability mass of the cells. It
illustrates the concentration of the mass in cell 1, with a high
$\log\left(\sum_{k=1}^{n_{\text{bin}}}\frac{f(X^{k})}{g(X^{k})}\right)$
compared to the other cells. It also explains the position of the centroid of
the cell 5, which is located near the left limit. This is due to the low
weight terms at the right of the cell.

(a)

(b)

Figure 11:  (a) Distribution of $10^{5}$ points sampled with density $g$ in
the different Voronoi cells in the two first FPCA axes of maps without breach.
(b) 2D histogram where for each bin, the quantity
$log\left(\sum_{k=1}^{n_{\text{bin}}}\frac{f(X^{k})}{g(X^{k})}\right)$ is
displayed with the color scale. The prototype maps are represented by the red
squares. The limits of cells 2, 3 and 4 are drawn with black lines. Grey
pixels are associated to a null weight
$\left(\sum_{k=1}^{n_{\text{bin}}}\frac{f(X^{k})}{g(X^{k})}\right)$.

##  Appendix E Relation between cells size and probability

As the Lloyd’s algorithm minimizes the quantization error
$e(\Gamma_{\ell})=\left[\mathbb{E}\left[\lVert
Y-q_{\Gamma_{\ell}}(Y)\lVert_{\mathcal{Y}}^{2}\right]\right]^{\frac{1}{2}}$,
we can expect that the largest Voronoi cells are associated to a small
probability mass. However, it may not be the case with specific random
variables $Y$.

For instance, let consider a random variable $Y\in[-20,20]$, with a density
function
$d=\frac{1}{10}\mathds{1}_{[-20,-10]}+\frac{9}{10}\mathds{1}_{[0,20]}$.

If we perform the Lloyd’s algorithm with two centroids $\gamma_{1}$ and
$\gamma_{2}$, it gives $\gamma_{1}=-15$ and $\gamma_{2}=10$ and two Voronoi
cells: $C_{1}=[-20,-2.5]$ with probability mass $p_{1}=\frac{1}{10}$ and
$C_{2}=[-2.5,20]$ with probability mass $p_{2}=\frac{9}{10}$.