Economic damages due to extreme precipitation during tropical storms: evidence from Jamaica

This study investigates economic damage risk due to extreme rainfall during tropical storms in Jamaica. To this end, remote sensing precipitation data are linked to regional damage data for five storms. Extreme value modelling of precipitation is combined with an estimated damage function and satellite-derived nightlight intensity to estimate local risk in monetary terms. The results show that variation in maximum rainfall during a storm significantly contributes to parish level damages even after controlling for local wind speed. For instance, the damage risk for a 20 year rainfall event in Jamaica is estimated to be at least 238 million USD, i.e. about 1.5% of Jamaica’s yearly GDP.


Introduction
Tropical cyclones (TCs) are among the most destructive natural disasters, estimated to have caused over 800 USD billion in damages globally over the last 20 years. 1 Damages attributable to these storms are mainly due to three factors, namely extreme wind, storm surge, and torrential rainfall (Bakkensen et al. 2018;Park et al. 2013;Lin et al. 2010). The literature modelling the economic impact has, however, mainly focused on damages as a result of wind, and to a lesser extent storm surge (Emanuel 2005;Nordhaus 2006;Hu et al. 2016; Baradaranshoraka et al. 2017;Masoomi et al. 2019;Hatzikyriakou and Lin 2018;Do et al. 2020), due to the difficulties associated with large-scale flood modelling as a result of extreme precipitation (Murnane and Elsner 2012;Zhai and Jiang 2014). This is often justified on the grounds that wind is strongly correlated with rainfall during a TC and thus that wind exposure will capture damages due to rainfall as well. However, recent evidence suggests that rainfall is not absolutely dependant on TC intensity, "...suggesting that stronger TCs do not have systematically higher maximum rain rates than weaker storms." (Yu et al. 2017). It has also been shown that the positive correlation between wind speed and rainfall may only be true over the ocean and not on land (Jiang et al. 2008).
The failure to take account of extreme precipitation in damage estimation arguably neglects an important driver of the economic costs due to TCs (Czajkowski et al. 2011;Rezapour and Baldock 2014;Rappaport 2014;Park et al. 2015;Bakkensen et al. 2018). For instance, available data for the Caribbean suggest that rainfall is either the primary or secondary cause of damages in 75% of TC events. 2 While the influence of anthropogenic climate change on TC frequency, general intensity, and affected areas is still a matter of debate, there is a general consensus that rainfall-heavy TCs will likely become more frequent in the future (Emanuel 2013;Knutson et al. 2019;Van Oldenborgh et al. 2017;Villarini et al. 2014b;Grossmann and Morgan 2011;Walsh et al. 2016). Nevertheless, the link between extreme rainfall and economic damages is not well understood yet (Villarini et al. 2014a;Rosenzweig et al. 2018;Rözer et al. 2019). While non-hazard measures such as risk awareness, building type, and topography are important aspects, the most influential determinants for whether a building is damaged by flooding are local water depth and accumulated rainfall (Spekkers et al. 2013;Van Ootegem et al. 2015, 2018Rözer et al. 2019). At the same time, short duration and high-intensity rainfall has been shown to be a major factor for the occurrence of landslides (Dou et al. 2019). Case studies have demonstrated that hazard scales that include rainfall in addition to wind speed are able to better predict the cost of TCs (Rezapour and Baldock 2014) and that wind only damage functions systematically underestimate damages by rainfall heavy typhoons in the Philippines (Eberenz et al. 2021).
The current study explicitly focuses on estimating the damage and risk associated with extreme precipitation during TCs, using Jamaica as a case study. Jamaica is arguably a particularly interesting setting for this purpose since the island is frequently afflicted by TCs. Excess rainfall is a major cause of destruction by inducing flash floods and landslides (Laing 2004). Moreover, due to the absence of large rivers, the volcanic origin of many hill-slopes, such floods and landslides, is common throughout most of the island (Miller et al. 2009). This allows the exploitation of spatial variation in rainfall intensity during TCs for the estimation of a rainfall-based damage function.
Importantly, the Planning Institute of Jamaica (PIOJ) has compiled consistent and detailed damage reports for most storms since the turn of the century. The lack of consistently reported damages often introduces additional uncertainty gitepbakkensen2018impact. Here, reports for five major TCs over the period of 2001-2012 are used to construct parish level information on economic damages. 3 Given that excess rainfall damages during TCs are often due to short, high intensity events (Larsen and Simon 1993). These damage data are coupled with high resolution, high frequency remote sensing precipitation data from the Global Precipitation Measurement Mission (Huffman et al. 2015a;Hou et al. 2014). Remote sensing information like the Tropical Rainfall Measuring Mission (TRMM) has been used extensively for the monitoring of TCs (Lonfat et al. 2004;Chen et al. 2006;Lau et al. 2008;Hendricks et al. 2010;Jiang et al. 2011;Villarini et al. 2011;Hence and Houze Jr 2011;Matyas and Silva 2013).
Using our data, we estimate a precipitation-dependent damage function via regression analysis, while controlling for the maximum wind speed during the storm. We then employ the precipitation data to generate return level maps using an extreme value model in the spirit of Demirdjian et al. (2018). Such a statistical approach has been shown to model extreme rainfall risk well in several case studies (Sugahara et al. 2009;Beguería et al. 2011;Tramblay et al. 2013). In conjunction with the estimated damage function and a proxy of the economic activity distribution derived from remote sensing nightlight intensity, information on the severity of a rainfall event occurring with a certain probability yields risk maps of economic damages due to extreme rainfall for Jamaica. The outlined method does not rely on large-scale hydrological models and should easily extend to other regions where gauge-based precipitation data are scarce.
The remainder of the paper is organized as follows: Sect. 2 presents the study region and describes the data. Section 3 details the methodology. Section 4 presents results while Sect. 5 discusses the findings. Finally, Sect. 6 concludes.

Study region
Jamaica is an island country situated in the Caribbean and is the third-largest island of the Greater Antilles after Cuba and Hispaniola. With a population of 2.9 Mio. (World Bank Group 2020), Jamaica is one of the larger states in the region and is ranked as a high development country by the UN human development index (Conceição 2019). Jamaica consists of 14 parishes, which is the highest regional administrative unit for the island. The country's economy is mixed but increasingly based on tourism and finance, while the export of agricultural commodities is declining (Johnston and Montecino 2012). Jamaica has two major cities, the capital Kingston in the south-east and Montego Bay in the north-west, known for its tourism. The islands geography is dominated by its central high plains and mountains. The Blue Mountains in the east, famous for their coffee plantations, constitutes Jamaica's highest point at 2.256 m. Jamaica is highly exposed to natural disasters such as TCs, earthquakes and droughts (Carby 2018). Between 1988 and 2012, 11 named storms made landfall in Jamaica and caused severe destruction (The World Bank Group 2018). The climate in Jamaica is tropical with little seasonality in temperature. However, due to the north-east trade winds, the period from November to March is dry and the rainy season lasts from April to October. Highest rainfall on the island occurs in the east, the north-eastern coast averaging more than 400 mm a year. Near the peak of the Blue Mountains, the average is more than 625 mm a year (Nkemdirim 1979). Summary of the data of Jamaica subsequently presented is given in Table 1.

Damages
The Planning Institute of Jamaica (PIOJ) has produced consistent socio-economic damage reports after damaging tropical storms since Hurricane Michelle in October 2001. In this regard, the PIOJ follows the United Nations Economic Commission for Latin America and the Caribbean (UNECLAC) Damage and Loss Assessment (DaLA) methodology (Bradshaw 2003). Each report contains a precise description of the event, preliminary rainfall reports, description of social trauma, emergency actions, and an assessment of economic damages. The assessment of economic damages considers damages in terms of three 'sectors,' namely agriculture (crop and livestock), public infrastructure (roads, schools and others), and housing.
Since not all reports provide this information at the parish level, the analysis in this study is restricted to events for which there is such a regional breakdown, i.e. five storms: Hurricane Michelle 2001, Hurricane Dean 2007, Hurricane Gustav 2008, Tropical Storm Nicole 2010 and Hurricane Sandy 2012. Additionally, even when there is regional information, in some reports, not all three of the main sectors are necessarily covered sub-nationally. Table 2 shows which sectoral damages have been reported by parish. When not all three sectors are covered sub-nationally, the damage measure was calculated as the proportion of damages with data at the parish level and scaled to total reported damages (including countrywide aggregates). For instance, during Hurricane Sandy in 2012, 25% ($3.7 M) of the damage to agriculture occurred in the parish of Portland and 30% (5,190) of all damaged houses were located there. Thus, it was assumed that Portland accounted for 27.5% or 29.4 M USD of the estimated total damage of 106.6 M USD across all three sectors. One should note that in most reports, the parish of Kingston is counted together with the parish of St. Andrew due to its small size and their proximity, leaving 13 distinct regions for the analysis. The per parish damages originally reported in USD are inflation adjusted with the Federal Reserve Economic Data (FRED) consumer price index (U.S. Bureau of Labor Statistics 2019, CPI), normalized to February 1st 2020 values. The data from the PIOJ reports for 13 Jamaican parishes and 5 TCs provide a total of 65 regional level damage observations.

Precipitation
The source for precipitation data is version 06B of the Global Precipitation Measurement (GPM) Integrated Multi-satellitE Retrievals (IMERG), a satellite based estimate (Huffman et al. 2015b). The satellite precipitation algorithm combines various microwave and infrared precipitation measurements to produce precipitation estimates, adjusted with surface gauge data. The resulting product is a half-hourly data set with near global coverage at a 0.1 • × 0.1 • resolution for the sample period 01. June 2000-31. May 2019. The data have been pre-processed according to the accompanying tech-report (Huffman et al. 2020). One should note that the GPM's predecessor, the Tropical Rainfall Measuring Mission (Huffman et al. 2007, TRMM), has been regularly used in the context of extreme value modelling in meteorology (Furrer and Katz 2008;Demirdjian et al. 2018) and hydrology (Collischonn et al. 2008;Li et al. 2012).

Storm tracks
Information on the TCs in our analysis comes from the HURDAT Best Track Data (Landsea and Franklin 2013). The data provide six hourly observations on all tropical cyclones in the North Atlantic Basin, including the position of the eye and the maximum wind speed of the storm. Additional information on the spatial extent of the TCs is taken from the Extended Best Track Dataset (Demuth et al. 2006).

Nightlights
The economic exposure to extreme precipitation during a TC is unlikely to be evenly distributed even within a parish. Ideally, this should be taken into account when estimating a damage function using parish level damages. Given the lack of localized economic data from official statistical sources, nightlight intensity is often instead used as an alternative proxy (Chen and Nordhaus 2011;Elliott et al. 2015). To this end, gridded nightlight activity data from the Defense Meteorological Satellite Program (DMSP) Operational Linescan System (OLS) are used. The DMSP-OLS is available monthly at a high resolution of 30 �� × 30 �� (approx. 1 km × 1 km). Since the nightlight data are at a higher resolution than the GPM and damage data, the former is employed to aggregate the latter. 4

Framework
The risk of natural disasters consists of a combination of multiple factors (Peduzzi et al. 2012): the hazard itself (frequency and intensity), exposure (economic value, number of people) and vulnerability (the degree of loss given a hazard). Few attempts haven been made to incorporate extreme rainfall into the TC damage function. Li et al. (2019) define a precipitation intensity index analogue to the common used wind speed-based power dissipation index. However, they do not derive or validate this functional form. Bakkensen et al. (2018) utilizes the natural logarithm of maximum TC lifetime-rainfall of any weather station in a province, analogue to their use of the natural logarithm of maximum wind speed anywhere in the province. While both approaches highlight the importance of rainfall as a potential source of TC damages, they mimic their functional form of max. wind speed. Our approach is to first construct parish level wind and rainfall measures by TC, taking local exposure into account via the use of nightlight intensity measures. These measures are then combined with the parish level damage reports to estimate a precipitation specific damage function. The coefficient estimates and an extreme value analysis of precipitation are then combined to create spatial risk maps for Jamaica.

Windfield model
Following Strobl (2011) in calculating the local wind exposure during a storm, the Boose et al. (2004) version of the well-known Holland (1980) wind field model is utilized. The model estimates the location specific wind speed by taking into account the maximum sustained wind velocity anywhere in the storm, the forward path of the storm, the transition speed of the storm, the radius of maximum winds and the radial distance to the storm's eye. The model further adjusts for gust factor, surface friction, asymmetry due to the forward motion of the storm, and the shape of the wind profile curve. The source of storm data used is the HURDAT Best Track Data (Landsea and Franklin 2013). These 6-hourly track data are linearly interpolated to hourly observations. WIND c,s,t , the wind experienced at any point c, during storm s at time t is given by: where V m,s,t is the maximum sustained wind velocity anywhere in the storm, T c,s,t is the clockwise angle between the forward path of the storm and a radial line from the storm centre to the c-th pixel of interest, V h,s,t is the forward velocity of the TC, R m,s,t is the radius of maximum winds, and R c,s,t is the radial distance from the centre of the storm to point c.
The remaining ingredients in Eq.
(1) consist of the gust factor G and the scaling parameters D for surface friction, S for the asymmetry due to the forward motion of the storm, and B for the shape of the wind profile curve. Points c are set equal to the centroid-coordinates of the GPM rainfall data and the storm-wise maximum wind speed for point c is given by MWIND c,s . Appendix A.2 provides additional information on the model parameters. (1)

Rainfall
Rainfall measures for the five storms are obtained for every GPM cell whose centroid is within the storm's reach at a certain observational time. 5 This is judged by comparing the radius of the outermost closed isobar provided by the Extended Best Track Dataset (Demuth et al. 2006) with the distance to the storm's eye. If the distance of a GPM cell's centroid to the storm's eye is smaller than the radius of the outermost closed isobar, that cell is currently affected by the storm. The Extended Best Track Dataset provides 6-hourly observations which we linearly interpolate to match the half-hourly GPM precipitation measurements. Then, all rainfall observations RAIN c,s,t can be summarized as the total rainfall SRAIN c,s or maximum rainfall during that storm, MRAIN c,s . These two statistics are the extremes for describing the full distribution of RAIN c,s,t .

Parish aggregation
The remote sensing information of rainfall, and thus, the calculated wind speed does not match the parish borders. 6 The monetary damages are, however, reported on the parish level. An intuitive way to aggregate rain and wind is to weight the cells by the share of nightlight activity in that cell relative to total parish nightlights. The grid of the nightlight activity is approx. 100× finer compared to the rainfall and wind speed grid. This allows one to down-weight cells that lie not entirely within a parish if weighted only by the within parish nightlight activity. It further controls for the value at risk of areas with more nightlight activity compared to areas that are unlit at night, e.g. large population centres. Cells are only added to the sum when their centroid lies within the parish. Specifically, j denotes a 30 �� × 30 �� nightlight activity cell and w j,s is the associated 3-month mean nightlight intensity prior to storm s and c is a 0.1 • × 0.1 • rainfall or wind speed cell. 7 The nightlight weight is the fraction of the sum of nightlight activity in cell c that lies in parish i, divided by the total nightlight activity in parish i, With these weights, the gridded data can be aggregated to the parish level: Even though many studies focus on landfalling TCs only, it has been shown that non-landfalling TCs can be equally destructive (López-Marrero and Castro-Rivera 2019). 6 The data of rainfall and wind speed are on a 0.1 • × 0.1 • grid. 7 Note that, many nightlight cells j are in one rainfall (wind) cell c and many such cells c are in one parish i.
where SRAIN i,s and MRAIN i,s are the parish i, storm s total and maximum rainfall and MWIND i,s the respective maximum wind measure.

Damage function regression
The functional form which is assumed for the damage function has to be imposed in the regression. In this regard, it is well established that the destruction of TCs relates roughly to the cube of the maximum wind speed (Emanuel 2005). It is used as control for a tropical storms' power dissipation (Bertinelli and Strobl 2013). A priori, it is unclear whether maximum rainfall during a TC or the total rainfall adequately represents the effect of heavy rainfall during a TC. It follows that both measures should be included. Additionally, differences in parish economic endowment and thus potential damage is accounted for by including parish indicators. The linear regression is then given by where DAMAGES i,s are the reported monetary damages, MRAIN i,s is the maximum hourly rainfall in a parish, SRAIN i,s is the sum of total rainfall, MWIND 3 i,s is the cube of the maximum wind speed and PARISH i are indicator variables for parish i. 8 Potentially, other functional forms like polynomials or nonlinear models might provide a better fit to the data. Given the small sample size (n = 65), more evolved models could overfit. 9 Estimates obtained with an ordinary least squares regression of Eq. 6 are likely to generalize well and are interpretable.
We run several alternative specifications of Eq. 6. In Eq. 7, we use parish population POP i measured by Jamaica's 2011 census to control for local economic endowment (Statistical Institute of Jamaica 2011) instead of parish indicators. This has the advantage of leaving more degrees of freedom for the model estimation. Robustness to outliers driving our results by omitting extreme observations as determined by Cook's distance (Cook 1977). 10 The data without outliers are then used to estimate Eq. 7.

Extreme value modelling
The objective of a statistical extreme value analysis (EVA) is to quantify the tail behaviour of a process, for instance extreme rainfall occurring at a certain location. Such an EVA can be categorized into either a block maxima or a peaks-over-threshold approach. Both allow to fit a distribution to the extreme values of any sequence of i.i.d. random variables but differ in how they make use of the available information. The half-hourly time series X i 10 We define a observation as outlier if it is 6 × more influential than the average observation. of rainfall in mm/h of GPM cell i with 333'072 observations from June 2000 to May 2019 are the subject of this EVA. Extreme rainfall in Jamaica happens often during TCs which occur 10-17 times in the Atlantic basin (NOAA 2019), not all affecting Jamaica. Due to this irregular pattern, a peak-over-threshold approach appropriately represents the physical phenomenon while making better use of the available information. The peaks-over-threshold approach is based on the Pickands-Balkema-de Haan theorem stating that threshold excesses y = (X − u|X > u) have a corresponding generalized Pareto distribution (GPD) if the threshold u is sufficiently high (Coles et al. 2001, p. 75): where The parameters of the GPD distribution are location , scale and shape . Here, a Poisson point-process (PPP) representation of threshold excesses is used, which allows one to a), directly estimate the likelihood function in terms of location , scale and shape parameter, b) model non-stationarity in these parameters, and c) include the Poisson distributed threshold exceedance rate together with the threshold excesses in the inference. The EVA is carried out on the cell level by separately estimating GPD(y i ) as in Eq. 8 where y i = (X i − u i |X i > u i ) is the threshold excess of GPM cell i. The vector X i contains the GPM rainfall measurements of cell i and u i is the cell-specific threshold derived in Sect. 3.3.1. Estimation is carried out via the extRemes library in R (Gilleland and Katz 2016).

Threshold selection
The selection of an appropriate threshold above which an observation is considered extreme is a classical case of the bias-variance trade-off. Too low a threshold violates the underlying asymptotic basis and leads to bias, while too high a threshold discharges valid observations resulting in an unnecessarily high variance. A plethora of competing threshold selection techniques exist (Scarrott and MacDonald 2012). The traditional procedure is graphical via mean residual life plot (Davison and Smith 1990, MRL) which is subjective and becomes quickly infeasible for a large number of time series. An alternative by Northrop et al. (2017) proposes the use of Bayesian cross-validation, comparing thresholds based on their predictive ability at extreme levels. This method scales well while directly addressing the desired property for a threshold-if the threshold is too low, threshold excesses will not follow a GPD and predictions will be off. For this study, the size of the posterior sample simulated at each threshold is set to be 50'000 and 100 different thresholds are considered for the estimation of the quality of predictive inference. The training thresholds correspond to the 0.95-0.9995 quantiles at each GPM cell, in increments of 0.0005.

Dependence
Rainfall tends to occur in temporal clusters, violating the independence assumption necessary for extreme value modelling. To overcome the issue of temporal dependence, the common method of declustering is implemented (Demirdjian et al. 2018;Gilleland and Katz 2006). The time series are declustered cell-wise such that the resulting series are nearindependent if the observations are sufficiently distant in time. Such a series will have a dependence structure with no effect on the limit laws for extremes. This "runs" declustering algorithm requires first some threshold u where the values below are not considered extreme (Coles et al. 2001). Second, it starts a cluster at every first entry v with v > u which runs until r consecutive observations are below u. Third, only the cluster maxima are retained and all other observations are rendered to zero. The same declustering scheme as in Demirdjian et al. (2018) is adopted, with a threshold equal to the 99-th percentile and r = 5.

Non-stationarity
The Pickands-Balkema-de Haan theorem requires the extreme values to be i.i.d. random variables which implies stationarity. Rainfall observations constitute by their very nature a non-stationary process. The PPP representation enables one to incorporate non-stationarity in the GPD parameters. In this regard, a first-order sinusoidal-function of the day of the year in the location and scale parameters is allowed, analogous to Demirdjian et al. (2018). 11 Additionally, a binary variable STORM c,t for cell c at time t is included to allow for a different distribution of rainfall observations that are within the spatial extent of a TC: Note that, for calculating return levels of a non-stationary EV model, one has to assume values for the non-stationary parameters. Here, this is done by fixing the day of the year t and storm dummy STORM c,t . The day of the year t is assumed to be t = 0 and to not influence multi-year return levels. Jamaica is on average for around 3% of the year within the outermost closed isobar of a TC. Thus, we assume that the frequency of TCs stays the same and set the storm dummy STORM c,t parameter equal to the sample probability Pr[STORM c = 1] ≈ 0.03 for the calculation of the return levels.

Monetary return levels
The estimated return levels can be used in the damage function to obtain the n year extreme rainfall return level in monetary terms. Location of economic endowment plays also a role in the monetary risk. Average nightlight activity W c,i serves as weight for economic endowment, analogously to Eq. 3, where it was employed to obtain rainfall by parish per storm. 12 Projected parish i damage for a return period of r is then given by where RAIN c,r is the estimated rainfall return level in cell c and is the coefficient of rainfall from the damage function estimation. PDAMAGES i,r is the projected extreme rainfall damage during TCs in parish i that is expected to occur every r years. These parish level estimates can then be aggregated to the country level and compared to Jamaica's GDP. Table 3 shows the mean damages, rainfall and wind speed per storm across parishes, whereas Appendix A.1 shows the damages by storm and parish. Hurricane Dean 2007 is associated with the strongest winds, while storm Nicole 2010 brought the heaviest rainfall to Jamaica. These two storms were also the most damaging of the five storms. Figure 1 displays the co-occurrence of rainfall, wind speed, nightlight activity and damages of Hurricane Sandy. Both rainfall and wind are the strongest in the north-east. Heaviest rainfall has been observed slightly more inland compared to wind speed, which peaked at the coast. While there is a strong relationship between the two, it is evident that these are two distinct physical phenomena with distinct spatial distribution. The map depicting nightlights clearly outlines the capital city of Kingston as the largest convolution of nightlights. A comparison with the map of damages suggests that parishes which undergo severe economic damages are also those with more nightlight activity, high winds and much rainfall.

Regression results
Estimates of the regression specified in Sect. 3.2.4 are displayed in Table 4. Baseline regression (1) (given by Eq. 6) shows max. wind speed to be a significant predictor of parish level damages while both maximum and total rainfall are insignificant. Since the latter two are highly collinear, 13 we consider them separately by estimating Eq. 6 without either MRAIN i,s or SRAIN i,s in columns (2) and (3), respectively. Accordingly, total rainfall during a storm is an imprecise predictor of damages. Maximum rainfall in regression (3), in contrast, is a statistically significant (10%-level) determinant of damages. Switching from parish fixed effects, which reduce the degrees of freedom considerably, to alternatively including parish population as in Eq. 7 for regression (4) yields virtually the same result. Furthermore, results do not change significantly in regression (5), based on the same Eq. 7 but without outliers as defined by Cook's Distance. The R 2 and adjusted R 2 are higher in models (1)-(3) compared to (4) and (5). 14 Since there are twelve parish indicators compared to one variable for population, the R 2 is expected to fall. The adjusted R 2 takes this The linear correlation coefficient is 0.9 and Spearman's rank correlation is 0.89. Note that, the other variables are near independent. Pearson's correlation coefficient is − 0.07 between max. wind speed and max. rainfall, 0.05 between max. wind speed and population and 0.11 between max. rainfall and population. Thus, we do not expect to have issues with collinearity besides the one already mentioned. 14 The R 2 is the proportion of variation in the dependent variable that is explained by the independent variables. The adjusted R 2 further takes into account the reduction in degrees of freedom. change in degrees of freedom into account-one would thus expect that the adjusted R 2 in (4) and (5) should not be lower than in (1)-(3). This is not the case. Likely, the parish indicators contain more information than just the population number and thus a model that contains these indicators will have a higher R 2 . From a model selection perspective, the Bayesian Information Criterion (BIC) favours model (4) and (5). 15

Threshold and parameter estimates
The average threshold chosen by the Bayesian cross-validation is 3.24 mm/h. However, the selected thresholds are spatially heterogeneous as depicted in Fig. 2. This highlights the extent to which meteorological conditions vary across the island. Two examples of the model fitting are given in Fig. 3. Panel (a) shows the diagnostic plot for the cell that covers the capital city of Kingston close to the Blue Mountains, while panel (b) shows the diagnostic plot for the cell with the lowest 40 year return level (which is a candidate for potential threshold-misspecification). The plots show estimated measures of predictive performance, normalized to sum to 1, against training threshold. In panel (a), the highest threshold weight and thus the selected threshold is at 0.56 mm/h, while panel (b) peaks at 0.49 mm/h. 16 We see that both plots give most mass of threshold weight to values in the range of 0.5 − 3 mm/h. As such, the specific choice of threshold appears not to have a large impact on the results as long as the threshold is within a certain range. Parameter estimates from the extreme value GPD models are shown as boxplots in Fig. 4. The estimated GPD's are shifted to the right (large ), smoothly decreasing ( ≫ 0 ), are not degenerate 17 and have a mean ( ≤ 1 ). Estimates of the non-stationary parameters can be found in Appendix A.3.

Return levels
Return levels are summarized in the maps in Fig. 5. The spatial pattern over different return periods stays the same with an increase in the average level as we extrapolate to less frequent events. The highest return levels can be found in the north-eastern parish of Portland, around the city of Port Antonio.
We use the coefficient of 0.22 from the fourth column damage function regression in Table 4 as the damage function parameter to calculate predicted monetary damages in 15 Note that, a smaller BIC means a smaller information loss relative to the true (data generating) model. Thus, one typically selects the model with smallest BIC. 16 The x-axis is scaled on the quantiles, thus the relative distance between the two plots is hard to compare. 17 If ≪ 0, then the GPD's support is 0 ≤ x ≤ −1∕ and is, in the case of extreme rainfall, degenerate.

Discussion
This paper makes three contributions. It adds to the recent literature that acknowledges the importance of rainfall in a TC damage function, especially with regard to climate change interactions. Employing the recent, high-resolution GPM rainfall data for an extreme value analysis in Jamaica identify the spatial distribution of extreme rainfall. Lastly, we propose  a simple approach to combine the damage function with return levels and transform them into monetary risk. The estimated linear coefficient suggests that if maximum rainfall during a TC increases by 1 mm/h in a parish, the inflicted damage increases by 0.22 Mio. USD. Depending on the TC, this can constitute a sizeable part of damages. On average for these five storms, ex post projected max. rainfall damages are 25.5% of total damages. 18 However, there can be considerable differences, depending also on the wind exposure during the storm: the average contribution of max. wind speed is 34%, while residual damages are 40.5%. 19 For example, during tropical storm Nicole in 2010, the average parish max. rainfall was 40.5 mm/h, which translates into 8.9 Mio. USD, while average parish damages were 21.9 Mio. USD. In contrast, when Hurricane Dean brought strong winds over the island of Jamaica MWIND 3 i,s ∕10 3 0.01 * * * (0.00) 0.01 * * * (0.00) 0.01 * * * (0.00) 0.01 * * * (0.00) 0.01 * * * (0.00)  18 That is calculated as the fraction of projected damages against total reported damages. 19 Damages not explained in the model come from factors such as storm surge and residual variation in the data.
in 2007, the average parish max. rainfall was only 18.3 mm/h resulting 4 Mio. USD damages (compared to the average 32 Mio. USD per parish), while 67% (21.4 Mio. USD) can be attributed to max. wind speed. This supports the view expressed in Eberenz et al. (2021) that wind-only damage functions fail to predict the destruction of rainfall heavy TCs. When comparing the observed events to the risk estimates, the extreme value analysis suggests that a rainfall heavy event like Nicole in 2010 has a lower return period than 10 years. If For details about the procedure to produce the plots, see Eqs. (7) and (14)  year return level to be 110 and 128 mm/h, respectively. In comparison, if we consider the location of these two airports (Montego Bay and Kingston), the 20 year return level estimated in this analysis is around 80 to 90 mm/h for both. 21 Unsurprisingly, we find that the area north of the Blue Mountains is most susceptible to extreme rainfall during a TC. This is likely because TCs often first make landfall in the north-east of Jamaica. Orographic lift of the Blue Mountains in the north-eastern parishes increases the prevalence of heavy rainfall in that area further (Laing 2004  Note that, the TRMM data is 3-hourly compared to the half-hourly GPM. 21 The cities are in more than one grid cell. Our risk estimates could arguably be useful for policy. More precisely, from a damage mitigation perspective, modelling the spatial distribution of risk is necessary for the planning and execution of ex-ante mitigation interventions (Holcombe and Anderson 2010;Anderson et al. 2011). For example, we find that the north-eastern parish of Portland is most susceptible to extreme rainfall and should therefore be a priority area for preventive measures such as proper surface water management to decrease landslide risk. Nevertheless, while our damage estimates serve to give an indication of the loss of assets from extreme precipitation during a TC, one should note that they only constitute part of the economic impact. More specifically, the direct damages are likely to result also in indirect damages through disruption of economic activity (Hallegatte and Przyluski 2010). Strobl (2012) for example estimates the macroeconomic reduction in output for the average hurricane strike in the Central American and Caribbean region to be at least 0.83 percentage points, although this estimate was solely based on damages due to wind exposure. The estimates here suggest that projected direct damages for a 5 and 40 year return period event are in value 0.9% and 1.9% of Jamaica's GDP. However, these figures cannot a priori tell us how the damages will translate into changes in aggregate output as this will depend on the indirect consequences, disaster relief and whether the damaged assets are replaced, among others. For instance, Lenzen et al. (2019)  disruption through the supply chain to regions and industries in Australia that were not themselves hit by the storm. These spillover effects could be an avenue for future research.
There are a number of weaknesses inherent in our methodological approach that future research could address. Firstly, the magnitude of the projected damages crucially depends on the damage function estimate for extreme rainfall. While the coefficient of interest in Table 4 does not change statistically significant across specifications (3)-(5), 22 there are still a number of caveats to be made in this approach. First, both the dependent variable (reported damages) and the independent variables are subject to measurement error. 23 Measurement error in the dependent variable is less problematic (if it is not systematic) since it "only" increases the unexplainable part of the regression and will simply result in less precise estimates. Measurement error in the extreme rainfall, our independent variable of interest, is more of a concern since a sufficient amount of such will induce attenuation bias of the estimates towards zero (Frost and Thompson 2000). A second caveat is the focus on statistical modelling and the omission of a hydrological perspective. Extreme rainfall during TCs is itself not a hazard but the cause for hazards such as floods and landslides (Yumul et al. 2012;Nolasco-Javier et al. 2015;Nolasco-Javier and Kumar 2018). Hydrological models would be able to decipher this relationship more precisely. The challenge of employing a hydrological model instead of a statistical approach is, however, in terms of feasibility and generalizability. Complete multi-hazard large-scale hydrological models are still an active research area (Lung et al. 2013;Koks et al. 2019). With the additional complexity of a model that makes use of the extent of river catchments, local soil type and hill slopes, any result from the analysis could only be in terms of these specific conditions. The precision gained by a realistic model is thus paid with a narrower transferability of the results. The last restriction concerns the valuation of non-monetary damages. More specifically, the monetary damages figures used here do not account for other nonmonetized impacts like the erosion and deterioration of soil which can be linked to extreme rainfall events (Rawlins et al. 1998), as well as psychological and mental health impacts to the exposed population (Lindell and Prater 2003;Bourque et al. 2006).

Conclusion
While the overall effect of climate change on the frequency and severity of TCs is a matter of debate and may be ocean basin specific, there is a general consensus that rainfall heavy TCs will likely become more common with global warming (Grossmann and Morgan 2011;Walsh et al. 2016;Knutson et al. 2019). In considering what greater precipitation in future TCs will mean in terms of economic impact, much of the current literature has focused on estimating the impact in terms of wind-induced damages. Even if one assumes that wind is a sufficient proxy for TC damages, a wind-only damage function likely underestimates future damage in climate change scenarios. The analysis presented in this case study of Jamaica shows that extreme rainfall during TCs is also potentially an important driver. Depending on the specification, a conservative estimate suggests that rainfall during a TC causes direct damage of 1.5% of GDP for a 20 year event. Thus, failing to take account of extreme precipitation in risk assessments may not only underestimate future damage but substantially bias current risk assessment.

Wind field model
In terms of implementing Eq.
(1), one should note that the maximum sustained wind velocity anywhere in the storm V m,s,t is given by the storm track data, the forward velocity of the storm V h,s,t can be directly calculated by following the storm's movements between successive locations along its track, and the radial distance R c,s,t and the clockwise angle T c,s,t are calculated relative to the point of interest c. All other parameters have to be estimated or values assumed. For instance, we have no information on the gust wind factor G, but a number of studies (see, e.g. Paulsen and Schroeder 2005) have measured G to be around 1.5, and we also use this value. For S, we follow Boose et al. (2004) and assume it to be 1. While we also do not know the surface friction to directly determine D, Vickery et al. (2009) note that in open water, the reduction factor is about 0.7 and reduces by 14% on the coast and 28% further 50 km inland. We thus adopt a reduction factor that decreases linearly within this range as we consider points c further inland from the coast. Finally, to determine the shape of the wind profile curve B, we employ the approximation method of Holland (1980) where B is assumed to be in the range of 1.5-2.5 and negatively correlated with central pressure. We use the parametric model estimated by Xiao et al. (2009) to estimate the radius of maximum winds R m,s,t .

PPP parameter estimates
Figures 7 and 8.  on the free R library extRemes. Data from remote sensing sources have been prepared in version 1.9 of the free Climate Data Operators (CDO) program by the Max-Planck-Institute on Linux. On reasonable request, R scripts that were used to perform the analysis can be obtained from the corresponding author.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.