Flood hazard model calibration using multiresolution model output

While our focus is calibrating an expensive computer model using model runs at multiple resolutions, we also address the computational challenges associated with calibrating computer models with high-dimensional spatial outputs. At the 50 m and 10 m resolutions, each model run contains observations at 14,214 locations and 126,791 locations, respectively. Although Gaussian process emulation-calibration can be useful for dealing with computationally prohibitive model run times (hall_etal_2011), high-dimensional model outputs pose computational challenges of their own (katzfuss2020vecchia; chang2014), particularly in training the emulator and generating emulator output. For example, fitting an emulator calls for a Cholesky decomposition of a dense positive definite matrix with $\mathcal{O}(\frac{1}{3}n^{3})$ floating-point operations (FLOPs); hence, emulation can be computationally expensive in high-dimensional settings. Many researchers have developed methods to address these challenges which are prevalent in environmental sciences (cite:10; cite:26). Methods to deal with emulation-calibration for high dimensional model runs include basis expansion approaches (cite:26; higdon2008), dimension reduction approaches (chang_binary; chang2014; cite:27), and composite likelihood approaches (cite:25). In this study, we extend the dimension-reduction approach of chang2014.

Consider ${\bf Y}$, a matrix of concatenated expensive and cheap model runs such that

with $C$ indicating cheap and $E$ indicating expensive. Here $p^{E}$ and $p^{C}$ are the numbers of expensive and cheap design points. We require cheap model runs at all of the same design points as the expensive model runs, i.e. $\{\mbox{\boldmath${\theta}$}_{1}^{E},...,\mbox{\boldmath${\theta}$}^{E}_{p^{E}}\}\subset\{\mbox{\boldmath${\theta}$}_{1}^{C},...,\mbox{\boldmath${\theta}$}^{C}_{p^{C}}\}$. Each design point contains a setting for channel roughness and river width error. So that ${\bf Y}^{E}$ and ${\bf Y}^{C}$ are observed at the same locations, we interpolate the cheap model runs onto the same locations as the expensive model runs using bilinear interpolation.

Following chang2014, we first center the columns of ${\bf Y}$ then apply singular value decomposition to obtain the scaled eigenvectors, ${\bf k}_{1}=\sqrt{\lambda_{1}}{\bf e}_{1},...,{\bf k}_{p}=\sqrt{\lambda_{1}}{\bf e}_{p}$, where $\lambda_{1}>...>\lambda_{p}$ are the ordered eigenvalues and ${\bf e}_{1},...{\bf e}_{p}$ are the corresponding eigenvectors of the sample covariance matrix of ${\bf Y}$. We select the leading $J_{y}$ principal components that explain 95% of the variation, then we form the basis matrix ${\bf K}_{y}=({\bf k}_{1},...,{\bf k}_{J_{y}})$. Next, we compute the $p\times J_{y}$ dimension-reduced matrix of model runs as:

where $\{{\bf Y}^{R}\}_{ij}$ corresponds to the $i$th design point and the $j$th principal component.

Next, we fit a modified version of the multiresolution emulator of jck_ml to each column of ${\bf Y}^{R}$. The emulator for all principal components $\mbox{\boldmath${\eta}$}({\bf Y}^{R},\mbox{\boldmath${\theta}$})$ is the collection of predictive processes for the $J_{y}$ principal components at ${\theta}$. These predictive processes are defined by the mean (Equation (10)), covariance function (Equation (11)), and emulator parameter estimates detailed in Section 3.2.2. We project the value of an expensive computer model run at new design point $\mbox{\boldmath${\theta}$}_{0}$ in the original space by computing $Y_{0}={\bf K}_{y}\mbox{\boldmath${\eta}$}({\bf Y}^{R},\mbox{\boldmath${\theta}$}_{0})$.

At the observation-level scale, the calibration (statistical) model can be represented as ${\bf Z}={\bf K}_{y}\mbox{\boldmath${\eta}$}(\mbox{\boldmath${\theta}$},{\bf Y}^{R})+{\bf K}_{d}\mbox{\boldmath${\nu}$}+\mbox{\boldmath${\epsilon}$}$, where the discrepancy ${\delta}$ is replaced by its kernel convolution representation ${\bf K}_{d}\mbox{\boldmath${\nu}$}$ (higdon1998) and ${\epsilon}$ is the usual observational error term. Here $\mbox{\boldmath${\nu}$}\sim N(0,\kappa_{d}I_{J_{d}})$ with $J_{d}<n$, where $\kappa_{d}$ describes the magnitude of the discrepancy. More information on selection of the kernel basis ${\bf K}_{d}$ can be found in chang2014.

In this study, we chose to leverage the low-dimensional space of the principal components as well as its computational benefits. Here, we model the dimension-reduced observations ${\bf Z}^{R}=({\bf K}^{T}{\bf K})^{-1}{\bf K}^{T}({\bf Z}-\mbox{\boldmath${\mu}$}_{\bf Y})$, where $\mbox{\boldmath${\mu}$}_{\bf Y}$ is the vector of column means of ${\bf Y}$ and ${\bf K}=({\bf K}_{y},{\bf K}_{d})$. Then the dimension-reduced probability model becomes

where $\mbox{\boldmath${\mu}$}_{\mbox{\boldmath${\eta}$}}$ and $\mbox{\boldmath${\Sigma}$}_{\mbox{\boldmath${\eta}$}}$ are the mean and covariance of $\mbox{\boldmath${\eta}$}(\mbox{\boldmath${\theta}$},{\bf Y}^{R})$. In our case, ${\bf Z}$ is an expensive computer model run contaminated by independent and identically distributed Gaussian noise. Then there is no need to model a discrepancy term, and the probability model for ${\bf Z}^{R}=({\bf K}_{y}^{T}{\bf K}_{y})^{-1}{\bf K}_{y}^{T}(Z-\mbox{\boldmath${\mu}$}_{\bf Y})$ becomes

We extend the multiresolution emulation approach from jck_ml which was developed for a temporal stochastic computer model to the context of a spatial deterministic computer model. Our adaptation of their approach models the expensive computer model runs as homoscedastic rather than heteroscedastic. jck_ml’s approach links the cheap and expensive computer model runs through a similar autoregressive structure to that of cite:13. cite:13 employs a multi-stage emulator training procedure that first estimates the emulator parameters for the cheap model then plugs in those estimates to the equation to estimate the emulator parameters for the expensive model. On the other hand, jck_ml provides a way to infer all emulation parameters simultaneously by modeling the concatenated cheap and expensive model runs in one multivariate normal distribution as shown in Equation (7). When training the cheap and expensive model emulators, simultaneous inference allows there to be feedback between the cheap and expensive model emulators. Moreover, jck_ml analytically integrates out the regression parameters for the emulator mean function. The resulting predictive processes using cite:13 and our adaptation of jck_ml’s approaches are similar in structure, yet differ in how they estimate their emulation parameters.

As in jck_ml, we aim to construct a high-fidelity emulator that accurately represents the true computer model. We follow their general emulation approach and notation. Here we fit an emulator separately to each column of ${\bf Y}^{R}=({\bf Y}^{R}_{1},...,{\bf Y}^{R}_{J_{y}})$, the principal components of the computer model runs. To simplify notation we denote column $j$ of ${\bf Y}^{R}$ as ${\bf t}={\bf Y}^{R}_{j}$. We also denote ${\bf t}^{E}={\bf Y}^{RE}_{j}=(Y^{R}_{1,j},...,Y^{R}_{p^{E},j})^{T}$ as the entries corresponding to the expensive model runs at design points $\{\mbox{\boldmath${\theta}$}_{1}^{E},...,\mbox{\boldmath${\theta}$}^{E}_{p^{E}}\}$, and ${\bf t}^{C}={\bf Y}^{RC}_{j}=(Y^{R}_{p^{E}+1,j},...,Y^{R}_{p^{E}+p^{C},j})^{T}$ as the entries corresponding to the cheap model runs at design points $\{\mbox{\boldmath${\theta}$}_{1}^{C},...,\mbox{\boldmath${\theta}$}^{C}_{p^{C}}\}$.

Following jck_ml, we define the covariance of ${\bf t}^{C}$ and ${\bf t}^{E}$ and the cross-covariance between ${\bf t}^{C}$ and ${\bf t}^{E}$ using the squared exponential covariance function.

By assigning the priors $(\mbox{\boldmath${\beta}$}^{C},\mbox{\boldmath${\beta}$}^{E})^{T}\sim{\mathcal{N}}({\bf b},{\bf B})$, where ${\bf B}$ is a block diagonal matrix with ${\bf B}_{C}$ at the top left and ${\bf B}_{E}$ at the bottom right, we can integrate out $\mbox{\boldmath${\beta}$}^{C}$ and $\mbox{\boldmath${\beta}$}^{E}$ (jck_ml). We define $\mbox{\boldmath${\Omega}$}_{-\beta}$ to be the set of emulator parameters without $\mbox{\boldmath${\beta}$}^{C}$ and $\mbox{\boldmath${\beta}$}^{E}$. The distribution of ${\bf t}$ conditional on $\mbox{\boldmath${\Omega}$}_{-\beta}$ is:

where ${\bf H}$ is the block matrix:

We compare these emulators’ performances both in the context of a 10-fold cross-validation study and in the context of an edge case study. In the 10-fold cross validation study 10% of the expensive model runs are held out at random, the emulators are trained on the rest of the model runs, and their predictive abilities are compared on the held out expensive model runs. In the the edge case study, we hold out the expensive model runs at design points where $RWE\in(1.045,1.05)$ or $n_{ch}\in(0.02,0.028)$ so that the only model runs at edge design points are cheap. We then train the emulator on the remaining model runs, and compare the predictive abilities of the MR and HR emulators at the held out design points. We perform both the cross-validation study and edge case study for the combinations of expensive and cheap model runs previously listed.

At the original scale of the model runs, we compare the performances of the emulators in terms of the root mean-squared error (RMSE). We define $D_{MR-HR}=RMSE_{MR}-RMSE_{HR}$, i.e. the difference in RMSE between the MR emulator and the HR emulator at the same design point. A negative $D_{MR-HR}$ indicates that the MR approach had smaller RMSE, and a positive $D_{MR-HR}$ indicates that the HR approach had smaller RMSE. In the cross-validation study (Table 1), the summary statistics for $D_{MR-HR}$ indicate that the MR emulator performs better in terms of RMSE in the 50e200c combination: here the median and mean $D_{MR-HR}$ are -0.03 and -0.033, while for the 100e400c combination the median and mean of $D_{MR-HR}$ are both 0.015. These summary statistics reveal a small difference in performance of the two emulators in terms of RMSE for either combination, as the unit of measurement is in m and flood heights projected by the computer model range from 0 to over 12.5 m.

We also consider the the values of $D_{MR-HR}$ for the edge case study (Table 2). For 50e200c but not 100e400c the MR emulator outperforms the HR emulator in terms of RMSE at the edge design points. For 50e200c, the difference in performances (median $D_{MR-HR}$ = -0.533, mean $D_{MR-HR}$ = -0.736) is much larger than for 100e400c (median $D_{MR-HR}$ = 0.035, mean $D_{MR-HR}$ = 0.054). In both the cross-validation and edge case studies, the values of the raw RMSEs reveal that the 100e400c HR and MR emulator have much greater predictive accuracy than the the 50e200c HR and MR emulators.

To sidestep the computational costs of evaluating other aspects of the emulator’s performance at the original scale, we consider the uncorrelated standardized prediction errors (USPEs) (bastos_ohagan2009) for the leading principal components for each emulation approach, as done by chang2014. We define leading principal components as defining at least 75% of the variation. We plot the uncorrelated standardized prediction errors against the theoretical quantiles of the appropriate T distribution (degrees of freedom = number of points in test set - number of parameters estimated in the mean) (bastos_ohagan2009), the emulator prediction, channel roughness, and river width error to compare MR versus HR emulator fit. We provide these plots for the cross-validation and edge case study in the supplement.

First we consider these plots in the cross validation study. For 50e200c, the $|USPE|$s tend to be smaller for the MR approach (Q1= -0.148, Q3=0.219, |Max|= 8.29) than for the HR approach (Q1= -0.614, Q3=0.462, |Max|= 15.2). Considering the QQ plots and confidence bounds, both approaches exhibit good fit in the middle and heavy tails. However the MR approach’s lower tail is closest to appropriate thickness. In the plots of the USPEs versus $n_{ch}$, we see that for both the MR and HR approaches there is no general pattern but slightly increased spread in some areas. In the plots of the USPEs versus $RWE$, there are ranges of $RWE$ with increased variability in the USPEs for both approaches: for the MR approach this is near the upper boundary and for the HR approach this is slightly above the lower boundary. For the MR approach, the plot of the predicted values versus the USPEs gives no indication of heteroscedasticity. For the HR appraoach, this plot shows slight increased spread for medium-high predictions. Overall, the USPE plots indicate better performance of the MR emulator.

For 100e400c, the $|USPE|$s again tend to be smaller for the MR approach (Q1= -0.075, Q3= 0.111, |Max|= 1.44) than for the HR approach (Q1= -0.570, Q3= 0.575, |Max|= 9.95). Considering the QQ plots and confidence bounds, for the MR emulator, the fit is good in the middle but unsatisfactory in the tails. For the HR emulator, the same is true except the lack of fit in the upper tail is more pronounced. Considering all plots of the USPEs versus $RWE$, we see that in both approaches there is no general pattern. In the plots of the USPEs versus $n_{ch}$, we see that for the neither approach there is a strong pattern, but for the HR approach there may be slightly decreased spread at the highest $n_{ch}$ values. Concerning the plots of USPEs versus the predicted values, for the MR approach there may be slightly increased spread at the highest and lowest predicted values, while for the HR approach there may be slightly increased spread at the middle predicted values. However, neither plot shows a strong pattern. Based on the USPEs, we show preference to the MR emulator. In the supplement, we provide a detailed examination of the USPEs for the edge case model runs for the 50e200c and 100e400c settings. Similar to the non-edge cases, the $|USPE|$s tend to be smaller for the MR approach than the HR approach; hence, the MR emulator is preferred over the HR counterpart.

Next, we compare the calibration results using the MR versus the HR approach. For both approaches, we ran the MCMC algorithm for 300,000 iterations to ensure satisfactory effective sample size ($\geq 3500$). We display the calibration results using the 100e400c combination of model runs in the main text (Figure 2) as these had greatest predictive accuracy in terms of RMSE for the cross validation and edge case studies. We provide the calibration results using other combinations of model runs in the supplement.

For the 100e400c combination when $\mbox{\boldmath${\theta}$}^{*1}=(n_{ch}^{*1}=0.0305,RWE^{*1}=1)$, $n_{ch}^{*1}$ is inferred nearly accurately by both approaches, but slightly overestimated by both. The uncertainty of the MR approach is slightly greater, and the HR approach places slightly more posterior density at $n_{ch}^{*1}=0.0305$. Both approaches yield a posterior mean of 0.0308. For both the MR and HR approaches, the posterior for $RWE^{*1}$ is fairly deconcentrated and contains the truth, but the MR approach is slightly more accurate, having its center closer to and more posterior density at $RWE^{*1}=1$. This is reflected in the posterior means from the MR and HR approaches: 1.011 and 1.015, respectively.

For $\mbox{\boldmath${\theta}$}^{*2}=(n_{ch}^{*2}=0.0249,RWE^{*2}=1.0452)$ with 100e400c, the HR approach estimates $n_{ch}^{*2}$ very accurately with posterior mean 0.0252. The MR slightly overestimates $n_{ch}^{*2}$, giving a posterior mean of 0.0283, with too sharp of a posterior such that almost no posterior density is allocated to $n_{ch}^{*2}=0.0249$. The MR approach estimates $RWE^{*2}=1.0452$ fairly accurately with much of the posterior density being allocated there, but with the posterior being skewed left, the posterior mean is 1.031. The HR approach yields a wider posterior density that contains $RWE^{*2}=1.0452$ but is not as informative, giving a posterior mean of 1.019.

We also obtain projections from the expensive model using thinned posterior parameter samples for each $\mbox{\boldmath${\theta}$}^{*}$ from each calibration approach. These projections come from using 100e400c since this combination generally produced the most reasonable emulation and calibration results across approaches. We thin the MCMC samples because obtaining a model projection at all sampled parameter values would be computationally infeasible. We selected 100 thinned samples to be consistent with the number of expensive model runs afforded for training the emulators. We run LISFLOOD-FP at the parameter values of each selected step. For each calibration method, we calculate the mean of the projections from the thinned posterior and call this the calibrated projection. We compare the mean projections to the expensive projections for each $\mbox{\boldmath${\theta}$}^{*}$.

For $\mbox{\boldmath${\theta}$}^{*1}$ and $\mbox{\boldmath${\theta}$}^{*2}$ we compare the calibrated projections for each method to the simulated observation using root mean-squared error, fit, percent bias, and correctness. We calculate fit and correctness by $\textnormal{fit}=(A_{rm})/(A_{r}+A_{m}-A_{rm})$ and $\textnormal{correctness}=(A_{rm})/(A_{r})$, where $A_{r}$ is the observed flooded area, $A_{m}$ is the projected flooded area, and $A_{rm}$ is the area both observed and projected to be flooded (Rajib_etal_2020). $\textnormal{Percent bias}=100\times(\sum_{j=1}^{N_{L}}(P_{j}-Z_{j}))/(\sum_{j=1}^{N_{L}}Z_{j})$, where $Z_{j}$ is the simulated observation at location $j$ and $N_{L}$ is the number of locations (PBiasSource). $P_{j}$ is the mean projection at location $j$ resulting from the given calibration method.

The plotted results (Figure 3) and the metrics (Table 3) show that for $\mbox{\boldmath${\theta}$}^{*1}=(n_{ch}=0.0305,RWE=1)$, both calibration approaches yield projections very similar to the simulated observation. The RMSE for the MR approach is slightly higher, but by less than 5 cm. Both approaches have nearly identical percent bias but the MR approach’s bias (-3.73%) is negative and the HR approach’s (3.72%) is positive, as reflected by the plots of the residuals in space. However, the fit indicates that the the MR calibrated projection was notably closer to the simulated observation in terms of predicting the same areas as flooded versus not flooded (98.5% versus 78.6% overlap). Correctness was equal to one for both approaches, so the calibrated projections from both approaches predicted all flooded areas to be flooded. In summary, if the primary interest is in projecting flooding in the right areas, the MR approach is preferable, but if the primary interest is in predicting high enough flood heights, the HR approach is preferable.

For $\mbox{\boldmath${\theta}$}^{*2}=(n_{ch}=0.0249,RWE=1.0452)$, the plotted results (Figure 3) and the metrics (Table 3) indicate that the MR approach yields a calibrated projection most similar to the simulated observation considering the plotted results and metrics. Both approaches yield calibrated projections with more flooding than in the observation, as reflected by the percent bias and in the plots of the residuals. Specifically, the RMSE is smaller, the absolute value of the percent bias is smaller, and the overlap in projected flooded areas was much larger. Based on calibrated projections for this simulated truth, we would recommend the MR approach.

With regards to computational cost, our approach is sensible for cases where the difference in model run times between the expensive and cheap model runs is larger. We compute the amount of time each calibration approach needed to reach 300,000 steps using the variable-at-a-time random walk Metropolis-Hastings algorithm for the 100e400c combination (Table 5). Generally, the multiresolution calibration methods are 5.79-6.31 times slower, likely due to the larger matrix operations associated with the additional cheap model runs and a more complicated emulator structure.

To obtain the LISFLOOD-FP model runs, obtaining 100 expensive model runs sequentially takes about 10.2 hours. Obtaining 400 cheap model runs sequentially takes 1.78 hours. Some of these model runs may be obtained in parallel depending on the computing system available, however parallelizing the expensive model runs makes model failure more likely. We also computed the amount of time needed for interpolating the cheap model runs and computing principal components and emulation for each approach. Interpolation took 7.81 minutes for the 100e400c combinations. Computing principal components took no more than 30 seconds for any approach. Fitting the emulators to all principal components sequentially takes 1.25 to 6.05 minutes for the MR approaches, and 5.86 to 9.36 seconds for the HR approaches depending on the combination of model runs.