2.1 Top-down catastrophe modeling

According to the United Nations Office for Disaster Risk Reduction (UNDRR) Sendai Framework Terminology on Disaster Risk Reduction (UNDRR, 2017), hazard, vulnerability, and exposure are defined as (IPCC, Reference Semenov, Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger and Caud2021a uses a similar terminology):

• “Hazard: a process, phenomenon or human activity that may cause loss of life, injury or other health impacts, property damage, social and economic disruption or environmental degradation;

• Vulnerability: the conditions determined by physical, social, economic and environmental factors or processes which increase the susceptibility of an individual, a community, assets or systems to the impacts of hazards;

• Exposure: the situation of people, infrastructure, housing, production capacities and other tangible human assets located in hazard-prone areas.”

Risk is at the intersection of hazard, vulnerability, and exposure. This is well illustrated in Figure 1 of UN (2023) or Figure TS.4 of IPCC (2021b). For example, a property is exposed to flooding if it is located in an area of flood hazard, whereas it is vulnerable to flooding if there are possible entries where water can enter into a house (basement windows or doors). Catastrophe modeling is based on such decomposition of risk and aims to model each of these three components, providing in the end what is known as the ground-up loss; that is, losses before the application of any insurance or reinsurance. This is illustrated at the bottom of Figure 1.

Figure 1. Top-down catastrophe modeling approach with climate on top.

The hazard component represents the frequency, intensity, duration, and footprint of an event. The exposure includes the geographical location of the property, its size (e.g., square footage, number of floors) and value (e.g., market value, reconstruction costs). Vulnerability represents the characteristics of a house that magnifies or attenuates the impacts of the hazard. In the context of flooding, this includes for example whether the basement is finished/unfinished or there is a crawlspace, the first floor elevation, and the height of basement windows, etc. Damage curves typically link the intensity of an hazard with the vulnerability of a home to yield dollars of losses or percentage of damage.

2.2 Climate models

Hazard modeling of climatic events such as flooding or tropical cyclones is founded on an understanding of the interactions between the climate and, for example, the frequency and intensity of a climatic hazard (top of Figure 1). CCRA adds another layer of modeling as we need to relate greenhouse gas (GHG) emissions and concentration to impacts on key climate variables such as temperature and precipitation. A natural approach is therefore the integration of climate models (general circulation models and regional climate models) into a PRA.

General circulation models (GCMs) are numerical models that simulate the evolution and interactions between the components of the climate system (atmosphere, land, ocean, ice, etc.) using physical equations and empirical relationships (Chen et al., Reference Chen, Rojas, Samset, Cobb, Diongue Niang, Edwards, Emori, Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger and Caud2021). They are at the core of climate and climate change studies and are thus widely used to study global temperature and precipitation (among other variables) patterns over the present and future (Chen et al., Reference Chen, Rojas, Samset, Cobb, Diongue Niang, Edwards, Emori, Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger and Caud2021). In many respects, regional climate models (RCMs) are similar to GCMs yet they model atmospheric phenomena at continental and regional scales allowing simulations at higher spatial and temporal resolutions, which can resolve smaller-scale physical processes (Chen et al., Reference Chen, Rojas, Samset, Cobb, Diongue Niang, Edwards, Emori, Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger and Caud2021). Climate models are computationally intensive and typically run on large supercomputing clusters.

Climate models are forced with GHG emissions scenarios that are designed to capture the impacts of future socioeconomic growth and energy alternatives. Those emissions scenarios are in turn converted into radiative forcings (see top of Figure 1). Changes in the radiative forcing represent the extra warming in the atmosphere due to GHGs and are measured in watts per square meter. For example, the IPCC AR6 SSP2-4.5 scenario, which represents a “Middle of the Road” future (see Fricko et al., Reference Fricko, Havlik, Rogelj, Klimont, Gusti, Johnson, Kolp, Strubegger, Valin and Amann2017 for the detailed storyline), corresponds to an extra energy flux of 4.5 $\mathrm{Wm}^{-2}$ to the atmosphere by 2100. Over available climate models, this emissions scenario typically leads to an approximate global warming of +2.5 degrees in 2100 compared to preindustrial.

Typical outputs of climate models include (surface) temperature, (liquid, snow, and convective) precipitation, (surface) relative/specific humidity, eastward and northward (surface) winds, (surface) air pressure, etc. Outputs typically are stored as grids that may take up to terabytes and petabytes of storage depending on the vertical (in the atmosphere), horizontal (over the surface), and temporal (hourly, daily) resolution of the data and the number of variables stored.

Integrating climate models into hazard modeling should also take into consideration the biases in the outputs that can affect the results at the end of the modeling chain. This is typically approached using what is called preprocessing or postprocessing. Preprocessing implies bias correction prior to its application in a hazard module. This is done by comparing, for example, simulated precipitation outputs with past observations and applying an appropriate correction. Postprocessing implies comparing a hazard feature simulated from climate model outputs (e.g., hazard frequency) with what was observed in the past. If preprocessing does not succeed in eliminating all biases in the hazard component, then postprocessing could also be applied.

Reporting and regulatory requirements are often based on assessing the overall impacts to an organization of a given temperature increase (e.g., an increase of 2 degrees compared to preindustrial temperature). Temperature increase and global warming are endogenous in a climate model and result from radiative forcings. Integrating climate models into PRA and CCRA therefore requires using available runs of climate models along with corresponding climate scenarios.

2.3 Proposed approach

We take a data science approach to CCRA for pluvial flooding in an insurance portfolio. The core of the work lies in the hazard modeling of pluvial flooding occurrence (Y) as a function of a set of atmospheric and socioeconomic variables (X). We first fit and validate the $Y|X$ relationship using statistical and machine learning approaches for binary responses using past observed data for X and Y. This is the fitting and validation step. Climate (change) risk assessment is then performed by computing flood occurrence probabilities using outputs from climate models for X over different time intervals. This is the projection and simulation step. When computing predictions, one may hold socioeconomic variables fixed to isolate the effects of climate change from socioeconomic growth, or one may use projections of the socioeconomic variables.

This study is not meant to provide a detailed account of the impacts of pluvial flooding at the street level. There is inevitably a trade-off between the financial and/or computing resources necessary and the resolution of the information needed. Here, we focus on the large-scale impacts of climate change on pluvial flooding in portfolios covering Canada or the US. As such, we do not explicitly model vulnerability, yet we proxy exposure as the number of people living or insured in each grid cell. More details about portfolio modeling are provided in Section 5.

3.1 Data

This section characterizes the historical flood occurrence data (Section 3.1.1) used as the response variable in the statistical and machine learning models. Then we examine the predictors derived from atmospheric and socioeconomic variables (Sections 3.1.2 and 3.1.3).

3.2 Methods

Flood occurrence is a classification problem and as such, we have used the generalized linear model (GLM) (namely, the logistic regression), the GAM, and random forests (RF). The GLM extends multiple regression models for, for example, binary or count responses. It does so by modeling a transformation of the expected response as a linear function of the predictors. The GAM expands on the GLM by allowing nonlinear transformations of the predictors. Finally, RF combines decision trees by sampling observations and predictors to derive empirical relationships. We focused on the latter three methods because the resulting models are flexible and interpretable while allowing for nonlinearities (GAM and RF) and interactions (RF). More details about GLM, GAM, and RF methods can be found in Chapters 4, 7, and 8 of James et al. (Reference James, Witten, Hastie and Tibshirani2021) and in the SM Section 1.1.

The historical occurrence data described in Section 3.1.1 has more than 2.4 million observations, of which 0.27% are ones (pluvial flood observed in a given month and grid cell) and the rest are zeroes (no pluvial flood observed in a given month and grid cell). As such, the dataset is imbalanced and we focus on avoiding problems related to the overestimation of the probability of no flood (false negatives).

There are few solutions to deal with data imbalance (Ganganwar, Reference Ganganwar2012). It is possible to oversample ones or undersample zeroes. Given the size of the dataset, undersampling zeroes was a more appropriate approach than oversampling ones, since it reduces the database size and accelerates computations. In other words, we randomly (over months and grid cells) eliminated zeroes from the dataset so that the resulting proportion of zeroes was either 90% or 50%. We tested two proportions to determine whether the outcomes are sensitive to this choice. Whenever undersampling was used, predicted flood occurrence probabilities were adjusted following Saerens et al. (Reference Saerens, Latinne and Decaestecker2002) to match observed probabilities.

3.3 Models

The response variable is flood occurrence in the US, measured over grid cells and months. We assume that given a set of covariates, flood occurrence is independent (over grid cells of $0.1^{\circ} \times 0.1^{\circ}$ and months) and as such the problem becomes a typical classification problem (a discussion on this assumption is included in Section 6). The set of covariates are (as described in Section 3.1) five precipitation variables (6-hourly, 9-hourly, 12-hourly, 24-hourly, and 48-hourly), temperature, snow cover, climate classification, eight land use proportions, and population density.

To eliminate the adverse effects of multicollinearity, we constructed a smaller set of covariates. We found high correlation in the precipitation variables (by construction, as one precipitation variable was often included in another) and only kept the 24-h precipitation since higher frequency precipitation is not available in the climate model projections described in Section 4.1. Moreover, we found high correlation between monthly temperature and snow cover, which led us to exclude snow cover since temperature is readily available in climate models whereas snow cover is not. We also found high correlation between the proportion of urban extent (one of the land use covariates) and population density and decided to keep the latter since it is readily available in future population projections. We also combined forest and grassland proportions into a single land use category for similar reasons. Overall, the smaller set of covariates comprises of 24-h precipitation, average daily temperature, climate classification, five proportions of land use (combining forest and grassland, dropping water, and urban extent), and population density.

Given the imbalanced dataset problem and correlation between many of the covariates (especially between the precipitation variables), we analyze each of the following using the GLM, GAM, and RF:

• All covariates, no undersampling;

• All covariates, undersampling with 90% zeroes;

• All covariates, undersampling with 50% zeroes;

• Smaller set of covariates, undersampling with 90% zeroes;

• Smaller set of covariates, undersampling with 90% zeroes, with logged population density.

Consequently, there are 15 models overall.

3.4 Validation

We first fit the 15 flood occurrence models (described in Section 3.3) over the US with 70% of data and assess the predictive capability using a test set made of the remaining 30% (Section 3.4.1). The test set has been generated by randomly selecting 30% of the 2.4 million observations, where one observation is either a 1 or 0 for each grid cell and month. In other words, the random sampling is achieved over both the temporal and spatial dimensions simultaneously. The training set is therefore made of the remaining observations. More details about the implementation of the GLM, GAM, and RF can be found in the SM.

A model fitted over the US (with Canadian KG climates available in the US) is then used to predict pluvial flooding in Canada. The model quality over Canada is investigated in two ways (Section 3.4.2). First, we used flood claims data from a major Canadian insurance company, which yielded a purely out-of-sample predictive analysis. Second, we perform a qualitative assessment of the model over major historical flood events in Canada.

3.4.1 United States

We can evaluate the predictive power of occurrence models by determining whether a model can predict an occurrence when one was observed, or vice versa, predict a non-occurrence when nothing occurred. This is done by computing true/false positive/negative rates. Popular approaches involve analyzing the receiver operating characteristic (ROC) and precision-recall (PR) curves. The ROC curve plots the true positive rate as a function of the false positive rate, whereas the PR curve focuses on the relationship between precision and recall, which are two different ways to express the true positive rate. There are however limitations about the ROC curve (Cook, Reference Cook2007; Saito and Rehmsmeier, Reference Saito and Rehmsmeier2015; Muschelli, Reference Muschelli2020) notably when the dataset is imbalanced, whereas the ROC is only a valid model selection criterion under auto-calibration (see Wüthrich, Reference Wüthrich2023 and references therein). We will therefore use the PR curve as an additional metric to assess models’ performance. For more details about the ROC and PR metrics, the reader should refer to Chapter 4 of James et al. (Reference James, Witten, Hastie and Tibshirani2021).

To evaluate the predictive power of the 15 occurrence models, we computed the areas under the ROC and PR curves with the test set in the US. The results are shown in Table 1.

Table 1. Area under the ROC and PR curves with the test set over the US for all 15 models considered. Note that “u/s” stands for undersampling.

We find that models perform very well in the test set with an area under the ROC curve in the range of 0.89–0.93, with a slight advantage to RFs. As for the area under the PR curve, values range from 0.05 to 0.11, which is above the baseline for a non-informative model (0.005, computed with historical occurrences). The GAM model performs better than the GLM, while RF again shows the best predictive capability in the test set. Across all models considered, RF yields the largest area under both the ROC and PR curves.

We find that there is no clear advantage from a predictive standpoint to undersample zeroes (Table 1). When using all covariates, undersampling with a target of 50% or 90% zeroes provided a very similar area under the ROC curve for the GLM, GAM, and RF. The latter result is different though for the PR curve, where undersampling slightly worsens predictive capability for the GAM, and more significantly for the RF. That said, the more parsimonious models with undersampling still yielded strongly comparable performance to cases without undersampling. When computation times matter, undersampling the dataset and using less covariates therefore yield very similar predictive performance (in a shorter amount of time and with less computational cost).

We compare in Figure 2 the predicted probabilities from one RF model with historical flood probabilities. We see that pluvial flood occurrence is concentrated in urban areas and that the model does very well in characterizing the spatial characteristics of pluvial flood, which is important to distinguish where climate change might have a greater impact. In the case illustrated here, the RF appropriately captures areas of low risk (white on top vs. yellow in the bottom) and pluvial flood probabilities in urban areas are at very similar levels and locations in both panels.

Figure 2. Flood probabilities over the US: empirical (Panel A, top) versus predicted (Panel B, bottom) using the RF model (undersampling with 90% of zeroes, smallest set of covariates, and logged population). Empirical flood probability is calculated as the number of months with flood occurrence over the total number of months. A white cell means no occurrence has been observed. Predicted flood probabilities are computed as an average over months and years between 2007 and 2020.

The SM includes 15 plots, one per model (GLM, GAM, and RF) and one per set of covariates (5), plus the empirical probability. They show that adding the logged population density was important for GLM and GAM, since it appears that predictions were too sensitive to slight changes in population otherwise. We also added to the SM a summary of model outputs for the GLM and GAM trained with the smaller set of covariates, with undersampling and logged population density. For both methods, precipitation is by far the most significant predictor. With the GLM, the coefficient for precipitation is positive and with the GAM the relationship is increasing nonlinearly. Other statistically significant predictors include population density, the proportion of wetlands and temperature, among others.

3.4.2 Canada

We would like to perform a validation exercise of the 15 models but over Canada. There are however no formal datasets in Canada that record flood events (or other weather events) at a level of granularity that we can find in the NOAA SED (with latitude and longitude of location). The Canadian Disaster Database (Canada, Reference Canada2022a), maintained by Public Safety Canada, has approximately the same level of information as the EM-DAT dataset (Guha-Sapir et al., Reference Guha-Sapir, Below and Hoyois2022). The Flood List website (Davies et al., Reference Davies, Behrend and Hill2021) also provides information about flooding globally, but in all three cases, location information is much too vague to be able to formally validate the flood models.

We thus perform a quantitative and a qualitative validation of the models over Canada. The quantitative assessment was based on a sample (non-random and non-divulged in this article) of data from Co-operators, a Canadian insurance company part of the top six P&C insurers in Canada. Data specific to clients were not used for this analysis, and only aggregate information about a flood event was compiled. Such assessment is feasible since the predictors used (Section 3.1) are available globally or over North America, with the exception of population which covers the US only. In this case, we used the Gridded Population of the World (GPW) v4.11, which is also available from the SEDAC (CIESIN 2018).

We have also recalibrated the 15 models over the US, but only over regions whose Canadian KG climates are available in the US (therefore excluding grid cells whose climate would not contribute significantly to predicting flood dynamics in Canada such as Southern US states). As such, we show ROC and PR metrics with the claims data available from 2012 to 2020. The qualitative assessment compares time series of predicted flood probabilities between 2007 and 2020 knowing that major historical flood events took place in Toronto (2013, 2018) and Calgary (2013, 2019).

Table 2 shows the area under the ROC and PR curves for the 15 models applied in Canada. In bold face, we highlight the method with the highest metric for each line. We now see a different picture, as is usually the case when performing out-of-sample prediction exercises. The areas under the ROC curves now range within 0.77–0.90, which is lower than what we obtained over the US. That said, the performance is very good, particularly for the models with the smallest set of covariates and logged population, with metrics in the range of 0.83–0.90. It is particularly surprising to observe a value of 0.9 with the fifth model under the GLM; it appears that simpler specifications are performing well out-of-sample in Canada and that logged population captures the reality that urbanization beyond some level should not have the same impact on pluvial flood probability. We therefore have a solid case for the fifth set of models (smaller set of covariates, undersampling with 90% zeroes, with logged population density) which has the highest scores while being the fastest to fit (because of fewer covariates and smaller sample size due to undersampling). As for the areas under the PR curves, RF have a slight advantage over the GAM.

Table 2. Area under the ROC and PR curves with flood claims from a Canadian insurer (2012-2020) for all 15 models considered. Note that “u/s” stands for undersampling.

We continue this section with a qualitative assessment of the models using selected flood events in Canada. Figure 3 shows major peaks in probabilities in July and August 2013 as well as in August 2018 which coincide with major flooding events in downtown Toronto. The 2013 floods in Toronto were among the most expensive for the insurance industry in Canada.

Figure 3. Validation of pluvial flood models with predicted flood probabilities in Toronto over July and August (top row), and Calgary over June (bottom row) between 2012 and 2020. Models with the smallest set of covariates, 90% of zeroes, and logged population density were used.

Looking at the bottom row for Calgary, we can distinguish significant peaks over June 2013 and June 2019. Although heavy rain triggered flooding in Calgary in June 2013, heavy snow accumulation upstream in the prior months magnified the intensity of the event, which might explain why the maximum observed in 2019 is higher than in 2013, because of the longer-term snow melt processes. Calgary also saw floods in June 2019 due to heavy rain from thunderstorms.

Moreover, in the three subplots of Figure 3, we observe that the GAM typically generates the largest range of flood probabilities, indicating that the model is the most responsive to changes in precipitation patterns. Finally, on the basis of the quantitative and qualitative validations, we find that the pluvial flood model fitted in the US provides strong predictions in Canada.

We conclude this section by providing a map of pluvial flood probabilities over Canada. Figure 4 shows that pluvial flooding is also concentrated in urban areas as was the case for the US. This is especially true in the Greater Vancouver area, Southern Quebec, and Ontario (including, Montreal and Toronto), as well as many urban areas of New Brunswick and Nova Scotia. Note that we cannot show historical claims patterns in Canada to protect the confidentiality of the clients. Moreover, we have more blank cells in Canada because more cells have no population (or too few) and because Northern Canadian climates do not exist in the US climate groups.

Figure 4. Predicted flood probabilities over Canada for the RF model (Panel A, top) and GLM (Panel B, bottom) using undersampling with 90% of zeroes, the smallest set of covariates, and logged population. Note that we cannot show historical flood probabilities to protect the confidentiality of the data. Similar plots for GAM are available in the SM.

4.1 Data

4.2 Methods

Because hourly precipitation is not available in the CRCM5 runs we analyzed, not all 15 models from Section 3.3 will be used for CCRA. Due to the unavailability of some variables, the predictive capability of the models, and computation times, we used the smaller set of covariates along with a target of 90% zeroes for each of the GLM, GAM, and RF. Covariates thus comprise 24-h precipitation, temperature, climate classification, five proportions of land use, and population density.

The first step consists in updating the covariates with output from the climate model. For each of the six runs of the CRCM5 over 2006 until 2060, we extract 24-h precipitation and compute the average daily temperature. Afterward, we record the maximum daily precipitation over the month and compute the average monthly temperature. The static variables such as Köppen–Geiger climate classifications and proportions of land use were held fixed until 2060 because no future projections were available. As for population density, we used two different assumptions: (1) SSP2 population projections until 2060 or (2) fixed from 2020 levels (SSP5 was also considered but results were not materially different over the time horizon considered). The latter thus fixes land use and urban extent and allows us to focus strictly on changes in future atmospheric conditions, whereas the former allows for interactions between increased urbanization and possibly more heavy rain.

The second step consists in computing flood probabilities with the updated covariates. For each month (12), year (54), and run (6) of the CRCM5, we computed flood probabilities using outputs of the CRCM5 as simulated covariates. We call these simulated flood probabilities, and they are available over the present and future climates. We interpret climate simulations over the present climate as alternate and plausible trajectories of the climate.

To mitigate the need to apply some kind of postprocessing (see Section 2.2) on flood probabilities, our analyses focuses on differences between two time periods (rather than looking at raw probabilities); therefore, assuming that any bias found in the CRCM5 over the historical period is likely to be of a similar order in future projections. The SM provides an analysis of the CRCM5 over 2007–2020, and we find that such bias is very small in most areas.

4.3 Maps

Our analysis first compares flood probabilities between two time periods: 2010–2030 (centered on 2020) and 2040–2060 (centered on 2050), which are 30 years apart. We therefore average simulated probabilities across months, years, and runs of each time period for both the US (Figure 5) and Canada (Figure 6). Both figures therefore highlight the combined impacts of climate change and future urbanization on pluvial flooding hazard.

Figure 5. Difference in simulated pluvial flood probability between 2040–2060 and 2010–2030 computed with the GLM (Panel A, top), GAM (Panel B, middle), and RF (Panel C, bottom) models over the US. Blank cells represent either too small population (in the past observations or future projections) or missing data.

Figure 6. Difference in simulated pluvial flood probability between 2040–2060 and 2010–2030 computed with the GLM (Panel A, top), GAM (Panel B, middle), and RF (Panel C, bottom) models over Canada. Blank cells represent either too small population (in the past observations or future projections) or missing data.

Whereas all three models agree that the West Coast of the US and Canada will be the most affected by changes in flood probabilities, there are however large discrepancies between predictions of the GAM/GLM and the RFs. The GAM/GLM families of models yield increases of pluvial flooding elsewhere in the US and Canada, concentrated in urban areas of Eastern US, Southern Quebec, and Ontario, whereas the RF shows close to no changes elsewhere. In fact, RF shows the smallest increases over the West Coast.

One should caveat the results from the RF models since RF and other tree-based methods are unable to extrapolate beyond the range of the training set (Hengl et al., Reference Hengl, Nussbaum, Wright, Heuvelink and Gräler2018). This becomes a major issue for CCRA because atmospheric variables such as precipitation and temperature patterns in the future may very well be different from their historical values.

4.4 Selected cities

We plot in this section the time series of the average (taken over months, grid cells of the city, and the six runs of the CRCM) annual simulated pluvial flooding probability for selected cities in the US (Figure 7) and in Canada (Figure 8). We observe increasing trends with different slopes across the different cities. Even though we averaged results over the six climate model runs, we still find substantial interannual variability. We also included in the SM similar plots for the GLM and RF. We find notably that there is nearly no trend with the RF and weak interannual variability, thus confirming the inability of the RF to extrapolate beyond the original training set.

Figure 7. Annual simulated pluvial flood probability from 2006 to 2060 over New York, Houston, Chicago and Denver with the GAM model. Similar plots for GLM and RF are available in the SM.

Figure 8. Annual simulated pluvial flooding probability from 2006 to 2060 over Montreal, Toronto and Vancouver with the GAM model. Similar plots for GLM and RF are available in the SM.

In both Figures 7 and 8, we also isolated the effects of climate change from increased urbanization with the continuous and dotted lines. That is, the dotted line represents a scenario where population remains fixed after 2020, whereas the continuous line represents a scenario where population increases according to the SSP2 scenario. In the latter scenario, the population of New York, Chicago and Denver will continuously increase in the future, whereas Houston should see a population decrease from 2020 to 2030 and an increase thereafter. Although the increasing trend in flood probability seems primarily driven by changing patterns in temperature and precipitation, urbanization also plays an important role on flood hazard.

5.1 Methodology

The occurrence model applied to the CRCM5 and SSP2 population projection yields simulated flood probabilities for each grid cell, month and year between 2006 and 2060. We can therefore directly use these probabilities to simulate monthly flood occurrences over the future. Yet an important piece remains, linking flood occurrence to impact in terms of losses.

The typical flood risk modeling chain (see e.g., Boudreault et al., Reference Boudreault, Grenier, Pigeon, Potvin and Turcotte2020 for fluvial flooding or Figure 1 for a general setup) entails using the location (latitude and longitude) and characteristics of each building in combination with a damage model (e.g., a damage curve linking water depth and dollars of losses) to represent losses from flooding. But with hazard information available on a 10–25 km grid in the form of flood occurrence instead of water depth, we do not aim to analyze impacts at the street level, and as such, the exact location of each building is not necessary for this exercise. For similar reasons, we will also ignore the vulnerability of each building and rather assume that each flooded property suffers a fixed or random loss amount. Aggregating exposure value, or the number of households insured per grid cell at a resolution similar than that of the climate model is straightforward for an insurer. But for this paper, we will rather build generic insurance portfolios based on the population data described in Sections 3.1 and 4.1.

It remains to determine the number of homeowners that are flooded when there is flood occurrence in a given grid cell. One can fix that number as a given percentage, but we instead modeled it as a beta-distributed random variable with fixed mean and a fixed upper percentile. As such, this adds randomness as to how extreme precipitation may locally affect a community, replicating the effects of spatial contagion within each grid cell.

The specific methodology is as follows. We have split the time horizon until 2060 into two time periods: present climate (2010–2030), centered around 2020, and future climate, centered around 2050 (2040–2060). A 30-year time horizon is reasonable for an insurer for strategic decision-making and solvency analyses, while avoiding the considerable uncertainty tied to considering climate up to 2100, which is heavily dependent on current climate policies. Each model run (6) and year (20) within each time period is assumed to be independent and identically distributed. This yields 120 climate simulations under the present climate and another 120 climate simulations under the future climate.

For the present climate, we draw 10,000 random numbers in order to randomly select a climate from the 120 available. For the selected climate, we compute simulated flood probabilities for each grid cell and month over Canada and the US. We then draw Bernoulli random variates according to these probabilities over each grid cell assuming that flood occurrences conditional upon the climate is spatially independent. If there is a flood in a given grid cell, we then randomly draw from a beta distribution with mean 2% and 99-th percentile equal to 20% to represent the percentage of homeowners affected. For each household affected by a flood, we assume a loss of $25,000. The value of $25,000 is based on the average damage given a pluvial flood per property, whereas 2% is selected to replicate industry losses. These choices do not make a material difference to understanding the relative impacts of climate change.

The previous steps were then repeated for the future climate as well. In both cases, we worked with the GAM and the smaller set of covariates fitted with a targeted 90% of zeroes. GLM yields similar results, whereas the RF has been excluded for the reasons explained in Section 4.

5.2 Results

To meet the objectives presented at the beginning of this section, we construct three scenarios for future change in the hazard and exposure. The baseline scenario represents our best estimate of the current loss distribution. It is based on present-day hazard (2020) and exposure (2020). The second scenario assesses the sensitivity of the insurer’s current exposure (2020) to changes in hazard (2050), including future projections of temperature, precipitation, and urbanization. It represents what would be typically asked for reporting and regulating purposes to assess the impacts of future pluvial flood hazard. In this case, the insurer’s portfolio is held fixed in the future, as if the insurer would not underwrite additional risks. Finally, the third scenario includes changes in both the hazard (2050) and exposure (2050), depicting a situation where a company underwrites in a similar manner and fixes its future market share instead. Such a scenario also highlights possible interactions between hazard and exposure where population could increase or decrease in riskier or safer areas. In all cases, we fixed the market share to 100% of the corresponding geographic region, which therefore proxies industry losses. Note that it is possible to create a fourth scenario in which hazard is fixed as of 2020 but exposure is that of 2050. Of lesser interest for insurance applications, this fourth scenario would support understanding the contribution of projected exposure to future flood risk. This is left for future research.

Table 3 shows the results of the latter three scenarios for four portfolios, fully underwritten in Quebec, Ontario, Canada, or over the US. Note that $ amounts have not been adjusted for inflation and reflect losses as of 2020. All risk measures were computed with 10,000 simulations, whereas the mean and standard deviations were validated with closed-form expressions that are straightforward to derive. The SM shows the equivalent of Table 3 but for each of the 10 Canadian provinces and the 10 most populous US states.

Table 3. Portfolio loss statistics for four portfolios and three scenarios for changes in hazard and exposure (in millions of 2020 dollars). Relative difference in % shown between parentheses (compared to the baseline scenario).

With the four portfolios illustrated in Table 3, we see that even in the aggregate, changes in hazard can be profoundly different across regions. Under the second scenario (hazard of 2050 but exposure of 2020), losses are expected to increase by nearly 50% in both Quebec and Ontario, whereas the increase is lower Canada-wide (about +40%) or in the US (+30%). Expressed differently, such increases represent 0.88% to 1.36% per year when compounded annually. Across states and provinces, SM Section 3 shows more homogeneity across Canadian provinces (increases of about 40–50% with the exception of BC and PEI) than in the US, where increases range from 15 to 50%.

Although we find benefits to diversification, country-wide effects of climate change on pluvial flood are expected to be relatively more significant in Canada than in the US. According to ECCC (Bush and Lemmen, Reference Bush and Lemmen2019), the average temperature increase in Canada is expected to be greater than in the US. This indicates that the air over Canada, being warmer, could hold more moisture that would in turn drive more intense rainfall, all other things being equal.

In the third scenario, both the hazard and exposure change in the future. There are however nontrivial interactions between changes in hazard and exposure depending on where population will live. Indeed, if for example current and future population move to areas with increasing hazard, then portfolio losses will increase at a faster pace than population growth. Table 3 shows that the third scenario yields losses much greater than the second scenario with significant heterogeneity. For example, the Quebec portfolio losses nearly double, while the US-wide portfolio losses increase by about 60%. Expressed on an annual basis, the compounding effects of increasing hazard and exposure mean that losses should increase by a rate of 1.6–2.2% annually. Across states and provinces, SM Section 3 shows variations between 40 and 95% over Canadian provinces and the top 10 US states, which is significant.

Significant trends in future losses should not be a surprise. Adding inflation of about 3% (which is slightly above the historical inflation over the last 40 years, but still below the inflation observed in 2022–2023) could yield a compound annual rate of increase in losses of over 5% (all else being equal, namely adaptation). Note that according to the Parliamentary Budget Officer in Canada, claims to the Disaster Financial Assistance Arrangements due to flooding have quadrupled over the last 40 years (Office of the Parliamentary Budget Officer, 2016), and therefore, the figures we present are clearly not unrealistic.

We conclude this section by analyzing the (kernel-smoothed) loss distributions in each of the three scenarios for the four portfolios. We clearly see rightward shifts in Figure 9 (as well as in the log-log version of Figure 9 provided in the SM) as we move from the baseline scenario (hazard of 2020, exposure of 2020) to the third scenario (hazard of 2050, exposure of 2050). The Quebec and Ontario portfolios are right-skewed and heavy-tailed, even more so than the Canadian and US portfolios. Judging by the upper percentiles, there does not appear to be a significant thickening of the right tail under the third scenario, but this is based on 10,000 simulations founded on 120 different climates, which could limit the potential to capture extreme losses. The conditional independence assumption may have played a role in limiting the potential for extreme portfolio losses; a longer discussion is included in Section 6. Therefore, we caveat interpreting and extrapolating future impacts of climate change on extreme portfolio losses based on Table 3 and Figure 9.

Figure 9. Probability density functions of portfolio losses for each portfolio and scenario.

This portfolio application shows the value for (re)insurance companies to invest in better underwriting practices and/or work with communities to attenuate the financial impacts of climate change on flooding. It also highlights the importance of quantitative CCRA to support such strategic decision-making at the organization level.