A Catalogue of Tropical Cyclone Induced Instantaneous Peak Flows Recorded in Puerto Rico and a Comparison with the World’s Maxima

: Peak streamﬂow rates from the Insular Caribbean have received limited attention in worldwide catalogues in spite of their potential for exceptionality given many of the islands’ steep topographic relief and proneness to high rainfall rates associated with tropical cyclones. This study compiled 1922 area-normalized peak streamﬂow rates recorded during tropical cyclones in Puerto Rico from 1899 to 2020. The results show that the highest peak ﬂow values recorded on the island were within the range of the world’s maxima for watersheds with drainage areas from 10 to 619 km 2 . Although higher tropical cyclone rainfall and streamﬂow rates were observed on average for the central–eastern half of Puerto Rico, the highest of all cyclone-related peaks occurred throughout the entire island and were caused by tropical depressions, tropical storms, or hurricanes. Improving our understanding of instantaneous peak ﬂow rates in Puerto Rico and other islands of the Caribbean is locally important due to their signiﬁcance in terms of ﬂooding extent and its associated impacts, but also because these could serve as indicators of the implications of a changing climate on tropical cyclone intensity and the associated hydrologic response.


Background
River flooding constitutes about a third of all geophysical hazards that annually occur throughout the world and it affects more people than any other type of natural hazard [1,2]. Between 1995 and 2004, floods were responsible for 20% of all worldwide natural hazard induced fatalities and for one-third of all economic damages [3]. Throughout the 1990s and the early 2000s, the annual number of flood disasters in the world doubled in relation to the period between the 1950s and 1980s [4]. Most flood-related deaths in the world are related to flash floods induced by tropical cyclones (TC) impacting landmasses along the western Pacific and Atlantic Oceans [5]. Although TCs also cause much harm and destruction due to high winds, waves, storm surges, and landslides [6,7], river flooding is the leading cause of TC-related deaths and economic loss [8,9]. TC-related river flooding occurs as a result of severe hourly rainfall intensities that can reach up to~150 mm h −1 [10][11][12], causing rivers to overflow [8].
River flooding impacts are expected to increase further due to human encroachment on flood-prone areas and to climate change-driven increases in rain intensities projected for many parts of the world [13,14]. This is particularly true for areas affected by TCs where largely unplanned development to accommodate rural to urban migration has been predominant [15] and where climate change-associated ocean temperature anomalies have already been linked to increased TC rain intensities [16,17]. In addition to increased rainfall rates, the weakening of summertime tropical circulation associated to climate change is potentially capable of slowing down TC translational speeds which can also lead to a greater likelihood of flooding [18]. These already noticeable manifestations of climate change reinforce the need for a better understanding of previously observed TC flooding magnitudes to aid in prescribing adequate adaptive measures [19].
This study concentrates on instantaneous peak flow rates (Q p ) which refer to the absolute maximum streamflow rates reached during rainstorms. Q p values considered in this study exclusively represent meteorological river flows as defined by O'Connor, et al. [20] in that they are " . . . notable flows of channelized water . . . induced by rainfall". Several Q p catalogues have been developed for different regions including the United States (U.S.), India, and China [21][22][23], as well as for the entire world [24,25]. Catalogues are used to define the extent of previously observed Q p maxima and encapsulate these in what is known as the flood envelope [26]. These datasets demonstrate an inverse relationship between area-normalized Q p (hereafter referred to as Q pa in units of m 3 s −1 km −2 ) and drainage area, as captured by the following equation: where a and b are fitted linear regression parameters, and A is area in km 2 [27]. Fitted values for b are always negative and greater than −0.66 for most regions [26,28] with a value of −0.4 found fitting to Mediterranean and inland continental settings in Europe, for example [28]. The decline in Q pa as a function of increasing drainage area is due to the area-scaling effect of precipitation and the attenuation of flow waves as they route through the fluvial network [20]. Flow wave attenuation is dependent on the physiographic factors related to stream channel hydraulic efficiency, channel slope, and the geometric configuration of the stream network [29,30]. Figure 1 shows a world envelope curve representing the maximum observed values based on the following equations developed by Herschy [24]: ln Q pa = 6.21 − 0.57 * ln(A) f or A < 100 km 2 4.61 − 0.20 * ln(A) f or A ≥ 100 km 2 (2) Figure 1. Relationship between area-normalized peak flow (Q pa ) and watershed drainage area (Area) according to the world envelope defined by Equation (2) [24] and for a selection of values from the Insular Caribbean. Values for Cuba, the Dominican Republic, Jamaica, and Puerto Rico (PR) are included in Herschy's (2003) catalogue. Values for Dominica are reported in Ogden [31] and those for the U.S. Virgin Islands (USVI) are from the U.S. Geological Survey's online database [32].
River flooding in the Insular Caribbean (IC) can be induced by a variety of rainstorm types including localized convective cells enhanced by orographic effects [33,34], cold fronts, troughs [35,36], and TCs. In fact, IC is one of the world's natural disaster hotspots [37,38] and river flooding is particularly important in terms of mortality, infrastructure damage, increased incidence of water-borne diseases, and for their lasting microand macro-economic impacts [39][40][41][42]. Climate change projections suggest a greater frequency and magnitude of extreme events for most of the IC [43,44]. In fact, climate records suggests a rise in the frequency of heavy rainfall events throughout the 20th century [45] and climate change is projected to induce a further increase in the number and intensity of major hurricanes for the region [46]. However, only a few Q p observations from the IC are included in worldwide catalogues of Q p maxima [25]. In the latest attempt to capture worldwide Q p maxima, out of hundreds of values, the IC was represented only by 10 entries from Cuba (only five TC related); six from the Dominican Republic (all related to Hurricane David in 1979), one non-TC related value for Jamaica, and four TC-related values for Puerto Rico (PR) (for a total of 21 values from the IC) [24]. Additional IC Q p values associated with TCs that merit recognition in worldwide catalogues include those provoked by TC Erika (2015) in Dominica [31] and several recorded in the U.S. Virgin Islands during TCs Hugo (1989) and Marilyn (1995) (Figure 1) [32], in addition to several more from PR.
The island of PR is no exception to the high propensity of TC rains, high peak flows, and impacts that typifies the IC [47,48]. In fact, many of the highest Q p values reported by the U.S. Geological Survey throughout the entire U.S. for watersheds smaller than~2000 km 2 have been recorded in PR [21]. The reasons why PR figures as a high peak flow magnitude hotspot include the relatively high density of stream gaging stations (e.g.,~55 active stations throughout the 1970s; 62 stations in 2017), the abundance of short and steep watersheds that can quickly transfer rainfall as runoff from hillslopes to streams [49][50][51], the high density and flow routing efficiency of its fluvial network [52,53], and the propensity of the island's precipitation regime to be heavily controlled by intense TC rainfall [34].
Even though non-TC storms can trigger significant peak flows [54][55][56], the individual intense rain cell lifetime of non-TC related extreme storms in PR is only about 20 min [57] and this is shorter than the time to concentration of most of the island's watersheds (i.e., 1 h for 10 2 km 2 watersheds, respectively) [58]. Therefore, TCs are presumed to provide the sufficient rainfall intensity-duration characteristics to generate some of the island's highest peak flows. However, other than describing the streamflow response of individual storms [58] and assessing peak flow rates from an average daily flow rate perspective [59], no study has catalogued the historical TC-related Q p values recorded on the island, nor compared them to worldwide maxima. Not all TCs generate extreme floods as rain intensities and totals are associated with the availability of abundant atmospheric moisture, proximity to the low-pressure center, TC translational speeds, and complex interactions with island topography [48,58]. Q p values are essential in understanding flood impacts in the type of flashy response that characterizes high-standing tropical islands like PR [60]. This became evident in 2017 during Hurricane María as instantaneous peak flow rates demarcated the extent of damaging flooding [61] and exerted a control on the stream channel geometric adjustments [62].

Objectives
This study intends to provide some context on the magnitude of TC-related Q pa values in PR. The specific objectives were to: (1) Develop a catalogue of instantaneous Q pa values and daily average rainfall rates for a selected group of TCs that affected PR from 1899 to 2020; (2) Compare the island's TC-related Q pa values to the world's flood envelope; and (3) Compare TC Q pa values and the associated 24 h rainfall rates for different regions of PR.
TC-related peak flows in PR are some of the highest recorded on the island (Ramos-Ginés, 1999), but also worldwide for watersheds with drainage areas ranging between 10 −1 and 10 2 km 2 [21,24]. Improving our knowledge of streamflow responses to TCs for the island is essential given that the projected increase in the frequency and magnitude of extreme events will undoubtedly have consequences on PR's future peak flow rates and associated flood damage [57]. However, it is not our intention here to develop an equation predicting Q pa as a function of explanatory variables.

Study Area
With a landmass area of about 8690 km 2 and located at roughly 18.25 • N and 66.50 • W, PR is the smallest and easternmost island of the Greater Antilles ( Figure 2). The island's physiography has been classified into three main provinces (i.e., upland, northern karst, and coastal plains) [63], which combined with the island's six subtropical life zones (i.e., from tropical dry forests to montane wet and rain forests) [64] and its wide variety of geologic terranes [65] has led to the recognition of a dozen distinct landscape units [66]. The Cordillera Central is the main topographic feature of the island's upland province traversing about two-thirds of its length before yielding into two northeast-and southeasttrending ranges (Sierra de Luquillo and Sierra de Cayey, respectively). Alluvial valleys typify the approach of most rivers towards all four coastlines, but the Caguas Valley is the only significant inland valley [66]. Even though forests currently cover about >40% of the island [67,68], most of these are secondary due to the widespread abandonment of sugar cane, coffee, tobacco, and cattle grazing lands that started during the second half of the 20th century [69,70]. Agricultural land presently covers only~11% of the island [71]. PR has a spatially varied yet consistently maritime tropical climate that is similar to other parts of the IC [72]. The island-wide mean annual temperature is~30 • C. PR-wide mean annual rainfall is 1690 mm year −1 and ranges from~700 mm year −1 in the southwest to~4600 mm year −1 in the northeast with~45% occurring during the peak of the TC season from August to October (Daly et al., 2003) (Figure 3a). Although not all of the largest rainstorms in PR are TC-related [73,74], TCs yielding an excess of 50 mm in mean island-wide rainfall occur on average 5-6 times per decade [48]. Major TC landfalls (>3 in the Saffir-Sampson TC scale) have occurred every five to six decades [50]. However, TC rainfall does not depend only on proximity to the low-pressure center, but also on the amount of precipitable water in the atmosphere, ground elevation, translational speed, and event duration [48]. Even though the largest amounts of TC rainfall have typically occurred in the eastern, southeastern, and central interior portions of the island [34,75], rainfall intensity and totals are largely dependent on the specific interactions of internal TC moisture, and TC wind direction with topographical features [58]. Runoff generation in PR is preferentially controlled by subsurface stormflow and saturation overland flow mechanisms for forested areas due to the high infiltration rates [76] and abundance of macropores [77] however, it is prone to the occurrence of excess precipitation overland flow on disturbed areas [78]. However, studies conducted elsewhere show that dramatic changes in land cover are typically undetectable by empirical flow data even in small watersheds (e.g., 35 ha watershed under controlled experiment conditions described by Birkinshaw, et al. [79]) and its effects dramatically diminish with both increasing drainage area and storm size [80]. Although PR underwent dramatic land cover changes during the 20th century [70], this study did not attempt to consider the effects of land cover changes because its focus is on Q p s induced by TCs for which land cover likely plays a minor role. Evapotranspiration occurring during the individual TCs is expected to play a minimal role in controlling Q p s given that the average daily evapotranspiration rates for the portion of the year when TCs predominantly affect PR (i.e., from August to November) ranges between 0.10 and 0.25 mm h −1 [81] and this equals only 3-7% and 0.8-2% of the average and maximum daily TC rainfall rates reported for the island, respectively, (3.5 and 12.7 mm h −1 , respectively [82]).

Development of TC Q p and 24 h Rainfall Rate Database
The USGS and other local government agencies have been operating stage-recording stream gauging stations in PR since 1958 with 10 and eventually 63 continuous recording stations in 1960 and 1970, respectively [83]. Automated gages have measured flow rates at 5-15 min intervals at stream reaches with drainage areas ranging from 10 −1 to 1.8 × 10 3 km 2 [84]. This study relied on data stored in the USGS' repository [32] and in USGS published reports (e.g., [85,86]) to construct a TC-associated Q p catalogue for PR. In this study, we included Q p values for 40 TCs that occurred between 1899 (i.e., TC San Ciriaco) and 2020 (i.e., TCs Isaías and Laura). For periods when recording stage-level equipment was active, Q p values rely on rating curves established at each station for flows that do not exceed stream channel capacities [87]. However, for flows that exceed the stream channel capacity and for those that occur when instruments were either not available or not operational, Q p values rely on field-checked high water marks and an application of the slope-area method [88] with these values prone to the highest errors [89,90].
We relied on the GIS-ready watershed polygons available through the USGS database to define the size and extent of the source area for each stream gauging station. However, these polygons are only available for the currently active stations (n = 67). For inactive stations, we developed a point-based GIS layer based on each station's coordinates and used that to define their drainage area using ArcGIS 10.7 s Flow Direction and Watershed tools on a sinkless, 90 m resolution DEM. Drainage areas were used to convert Q p (in m 3 s −1 ) into Q pa units (m 3 s −1 km −2 ).
The new catalogue allowed us to compare the various TCs with regard to their Q pa values. Since Q pa is strongly dependent on the drainage area (Equation (1) and Figure 1) and because the database represents observations from a variable collection of stream gaging stations with different drainage areas, a simple difference of mean Q pa values is not an appropriate comparison metric. Therefore, our comparisons rely on the criteria used by O'Connor et al. [20] to identify significant flows: where Q p sig is the cutoff value to declare a flow as significant in units of m 3 s −1 and A is in km 2 .
The area-normalized version of this equation is: where Q pa sig is in m 3 s −1 km −2 . Two metrics were used to assess the overall Q pa magnitude of every TC in our record. The first consisted of calculating the proportion of observations recorded during each storm that exceeded their respective Q pa sig values. The second assessed each TC based on its average Q pa /Q qpa sig ratio. It is important to emphasize that this is not meant to be a formal comparison of TCs but simply an aid to identify TCs with particularly high Q pa responses. Rainfall associated with each peak flow rate was also captured to further characterize every Q pa value and TC. Rainfall data relied on the interpolated rasters of maximum daily rainfalls for each TC calculated as the average 24 h rainfall intensities in mm h −1 (hereafter referred to as R 24h ) and developed by Ramos-Scharrón and Arima [82]. We relied on R 24h as rasterized hourly rain data were not available and because hourly rain data did not cover the entire period of interest since recording rain gauges were first installed on the island in the 1970s. The collection of storms with rasterized data included all TCs with direct hits on the island since 1899, TCs for which historical accounts and reports described significant rainfall and flooding, and TCs from 1970 onwards that produced more than 50 mm of rainfall [82].
We relied on the upstream drainage area polygon of each stream gauging station to calculate the average daily rainfall (R 24h ) corresponding to each Q pa value. This was calculated using the mean option of the Zonal Statistics as Table tool on the rasterized daily rainfall dataset [82]. Matching R 24h estimates were not calculated for all Q pa values because we lacked the coordinates for some inactive stream gaging stations or because of an inability of the GIS tools to define a source drainage area for some stations. R 24h values were used to identify TCs with exceptional daily rainfall. A value of 2.8 mm h −1 was chosen as a R 24h threshold to define exceptional storms as this represented the median island-wide daily rain intensity of 61 TCs for which this information was readily available [82]. We informally compared TCs based on the proportion of the R 24h values that exceeded 2.8 mm h −1 and on the average of all R 24h /2.8 mm h −1 ratios for each storm.

Comparison of PR's TC Related Instantaneous Peak Flows with the World's Envelope
We assembled a world instantaneous peak flow maxima database combining four existing catalogues. The database included values compiled by Costa [22] for the U.S., Li, et al. [91] for China, Rakhecha [23] for India, and Herschy [24] for the world (Supplementary Materials). The database considers only rain-associated flows and also contains drainage area information for every measurement location. Only values for drainage areas within 6.5 × 10 −3 and 2450 km 2 were considered for our analyses as these values span the drainage areas represented by flow rates values recorded in PR. Q pa and drainage area values were transformed into their natural logarithms before conducting linear regression analyses to calculate the best fit coefficients for both parameters a and b following Equation (1).
The PR peak flow envelope was estimated through a quantile regression model. Whereas ordinary least squares (i.e., linear regression) applied to the world envelope estimates the average relationship between peak flow conditional on area, quantile regressions can identify such relationships for different quantiles in the distribution. In our case, we estimated this relationship at the 95% quantile. This is similar to finding a linear model that splits the top 5% peak flows, conditional on upstream drainage area, from the bottom 95%. Following Cameron and Trivedi [92], the quantile regression estimator Q(β q ) finds the β coefficients that minimize the sum of the absolute deviation of the error according to the objective function: where q is the quantile (0.95), which pre-multiplies the absolute deviations in the formulation, placing more weight on predictions for portions where ln(Qpa) ≥ ln(A)b q + a q ; and a q is the intercept coefficient. The coefficients are subscripted with q to emphasize that each quantile has its own set of coefficients. This objective function is not differentiable, meaning that gradient-based algorithms cannot be used, and the solution (i.e., coefficients b q , a q ) that minimizes the function is found via linear programing. The standard errors of this regression were estimated using 1000 bootstrap resamples of the original set of observations with replacement. For each resample, the method runs a quantile regression and stores the coefficients of interest. Once all 1000 regressions are calculated, the final coefficients estimates are the average of all estimates and the standard errors are the square root of the variance of the bootstrapped estimators [93]. For this particular analysis, we only used those peak flows that were classified as significant according to O'Connor's Equation (3b), yielding 569 peak flow observations. Quantile regressions (QR) offer several advantages over ordinary least squares for this particular application. First, as described above, we can find the relationship between peak flow and watershed area for only the top 5% of our dataset. Second, the method is more robust to outliers (both positive or negative) because it minimizes the sum of absolute deviations instead of the sum of squared deviations. Third, QR is found to be more robust to non-normal errors, which is an important factor when the sample size is small. QR is also less parsimonious in the selection of what observations of peak flow are included in the regression. For instance, in our subsample, we included all peak flows that were "significant", including multiple entries for certain watersheds. Had we used standard OLS to model Q pa as a function of area, we would have to make an arbitrary decision on what observations from this subset to include in the regression. One option would be to only keep the highest values for each watershed and use OLS, however, questions remain as to whether, for example, two watersheds of similar size exhibit maximum values that are disparate, i.e., one with a high maximum value and one with a much lower maximum value. An arbitrary decision here would have to be made to either include the lower maximum value, which would lower the predicted fitted line, or bin the two watersheds into one single watershed class and only use the highest value of the two. This arbitrary selection problem was avoided with the use of QR because the method is less sensitive to outliers (low values in this case). An additional benefit of including more observations is the higher degrees of freedom. Finally, QR is invariant to monotonic transformations such as our case where the peak and area were log-transformed [94]. This means that all statistical results can be translated back to Qpa, which is not the case of OLS because the expectation operator does not "pass through" non-linear transformations (e.g., E[ Ln(Qpa)|Ln(area)] = Ln(E[Qpa|area]) ).

Regional Q pa and R 24h Tendencies within PR
Previous work has identified a propensity for some regions of the island to receive greater TC rainfall than others and is therefore logical to assume that these could lead to higher Q pa values. In fact, peak flow magnitude frequency analyses conducted for the island have established a positive relationship between Q p and the mean annual rainfall [84,95]. Annual rainfall on the island is highest in the vicinity of the Sierra de Luquillo in the Northeast and near the southeast hills marking the eastern end of the Sierra de Cayey, with some localized areas close to the highest peaks of the Cordillera Central [96] ( Figure 2). These regions mostly coincide with those where TC rainfall tends to be most abundant [34,48]. Here, we evaluated differences in both Q pa and R 24h for watersheds within different regions in PR. Some studies have zoned PR's climate into sub-areas based on convective storm hourly precipitation intensities [97,98] while others have divided the island based on both daily temperature and rainfall [99]. However, here we rely on the National Weather Service climate forecast zones (NWS zones; Figure 3b) [100], as these account for regional differences in rain patterns while also adhering to topographic boundaries to fulfill its streamflow forecasting purpose [101]. The expectation was that the northeastern, southeastern, eastern interior, and central interior zones (NWS zones 2, 3, 4, 5, and 6) should display higher Q pa and R 24h values than others. Source areas represented by each peak flow and daily rain observation were assigned to a single NWS zone using the majority option of ArcGIS' Zonal Statistics tool. Q pa and R 24h observations were grouped into their corresponding zones. Q pa observations for the different NWS zones were compared based on the overall average Q pa /Q pa sig ratio described above. After testing for normality (Shapiro Test) and equal variance (Levene's Test) we relied on one-way Kruskal-Wallis and Dunn's non-parametric tests [102] to compare the peak flow responses within each of the different NWS zones represented in our dataset. Given that average precipitation falling over a zone is not as strongly dependent on the drainage area as it is for Q pa (Average rain over an area = f (A −0/05 ) [20], we simply relied on the raw R 24h values to compare the differences among the zones also based on one-way Kruskal-Wallis and Dunn's tests after testing for normality and equal variance.  A total of seven out of the 11 NWS zones are represented by Q pa and R 24h values ( Table 2) (i.e., Zones 1,2,3,4,6,7,and 9). Each zone has between 40 and just over 600 Q pa and R 24h entries, with the lowest number of observations being for Zone 1 (San Juan) and the highest for Zone 4 (eastern interior). Historically, Zone 4 has contained the densest network of stream gaging stations due to the importance of its largest watershed (i.e., Río Grande de Loíza) in supplying much of the water needs of the San Juan Metropolitan area [103,104]. With the exception of Zone 1, all of these seven zones have observations from watersheds with drainage areas spanning over two to three orders of magnitude.  [20] (Figure 4a) and about 61% of the R 24h observations exceeded the 2.8 mm h −1 daily rainfall criteria (Figure 4b). The overall average Q pa /Q qpa sig and R 24h /2.8 mm h −1 ratios were 0.96 and 1.74, respectively. Notable TCs based on normalized instantaneous peak flows and daily rainfall include San Ciriaco, San Ciprián, and San Felipe II (Figure 5a-d). However, it is important to note that these three TCs are represented by only 3-9 observations (Table 1) and that these values most likely only exemplify the most extreme conditions that occurred on the island during each storm. Out of those three, Hurricane San Ciriaco is definitely a notable TC that still holds record rainfall records for portions of PR [82].

Worldwide and PR Q pa Envelopes
Linear regression analyses of the world's Q pa maxima (Figure 4a) led to the following result ( Figure 6): The 95% confidence intervals for coefficients a and b were 4.659 to 4.988 and −0.367 to −0.304, respectively. The b coefficient is only slightly greater than both the −0.4 calculated for Europe [28] and the values developed by combining values from the continental U.S., Hawaii, and PR [21] (−0.47 to −0.43). Trends in the data did not justify having separate regression curves for watersheds with drainage areas less and greater than 100 km 2 as performed by Herschy [24].
Application of the 95% quantile regression method to the PR database yielded this equation (Figure 6): Resulting a and b coefficients had 95% confidence intervals of 3.908 to 4.589 and from −0.358 to −0.219, respectively. A comparison of the world and PR peak flow envelopes suggest that recorded Q pa values in PR are below the low end of the world envelope for watersheds smaller than~10 km 2 (ln A =~2.3). This could be in part because only~5% of PR's database comes from watersheds smaller than 10 km 2 and the likelihood of recording the most extreme conditions is therefore restricted. The PR envelope is closer to the range of the 95% confidence interval of the world envelope for watersheds greater than 10 km 2 . A total of 20 Q pa observations for PR are within the 95% confidence interval of the world envelope ( Table 3)

Q pa and R 24h Values by Region
When grouped into the seven NWS forecast regions represented by the PR catalogue, Zones 2 and 4 displayed the highest average Q pa /Q qpa sig ratios among all zones. These were followed by Zones 1, 6, and 3, successively (Figure 7a,b). Zones 7 and 9 had the overall lowest values. However, the only statistically significant differences in average Q pa /Q qpa sig ratios were between Zone 9 and Zones 2, 3, 4, 6, and 7.
When ranked according to the average R 24h , Zones 3 and 4 displayed the highest values followed by Zones 2, 6, and 7 (Figure 7e,f). Statistically significant differences in average R 24h occurred only between Zone 9 and Zones 2, 3, 4, 6, and 7. These findings are in general agreement with previous studies that have suggested that higher TC rainfalls are typically experienced in the north-, south-, and central-eastern regions of the island due to the general westward approach of TCs towards the island and the depletion of precipitation potential as the storms move over PR due to the orographic effects of the Sierra de Luquillo, Sierra de Cayey, and the eastern end of the Cordillera Central [34,75]. However, most differences were not statistically significant.
Our findings also confirm the tendency for higher Q pa values to occur in areas with higher average annual precipitation (Figure 7c) [84,95]. However, Zone 9 displayed relatively low average Q pa /Q qpa sig values relative to its high annual rainfall. Annual rainfall for this zone is less controlled by TC rainfall than other parts of the island due to the greater relevance of localized convective storms in controlling rainfall patterns in the western parts of PR [73]. In contrast to the spatial pattern of Q pa observations, the highest average R 24h values were not consistently found in areas with a higher annual rainfall (Figure 7e-g). The highest average R 24h values occurred in Zones 3 and 4. The average R 24h in Zones 1 and 9 was relatively low in relation to all other zones. Figure 6. (a) The world envelope regression results (Equation (5)) and its 95% confidence interval. (b) The Puerto Rico instantaneous peak flow regression results based on 95% quantile regression (Equation (6)) and its 95% confidence interval. (c) The world and Puerto Rico peak flow envelopes (Equations (5) and (6)). Observed Q pa /Q qpa sig and R 24h maximum values in the PR catalogue showed a different spatial pattern than that based on average values, with Zone 6 (i.e., central interior) having the highest overall value. From this perspective, the central western, interior, and eastern zones as well as the southeast regions are the ones that recorded the most intense R 24h and the highest Q pa values on the island. The overall highest Q pa /Q qpa sig value in our database was registered in the central interior region (Zone 6) at the Río Grande de Manatí at Ciales station during the passing of Hurricane María in 2017 (Table 3). Another notable flow was registered in the southeastern region (Zone 3) during Hurricane Georges at the Río Grande de Patillas near Patillas station. The Río Valenciano, located in the eastern interior region (Zone 4) registered an extreme flow event during the passing of Hurricane Donna. Finally, TC San Ciriaco provoked a very high Q pa at the Río Grande de Arecibo-Dos Bocas dam station (i.e., prior to the establishment of the dam) located in the western interior region (Zone 9). Even though there is a clear propensity for TCs to provoke higher flows in the eastern half of the island than elsewhere, maximum Q pa values that are comparable to the most extreme values reported in the world have occurred within various regions of the island.

Implications on Flood Management for Puerto Rico
Flood management decisions somewhat rely on the expected frequency of instantaneous peak flows of a given magnitude [105]. This is true in PR where instantaneous peak flows with an expected recurrence interval of 100 years and 500 years are used to define the flood risk areas for the purposes of planning and insurance [106]. Calculating the magnitude of these events depends on the accuracy and completeness of past records [107]. Peak flow recurrence interval equations used to develop advisory information to the Federal Emergency Management Administration and the PR Government after the passing of Hurricanes Irma and María in 2017 [108] were those developed by Ramos-Ginés [95]. These equations presume a higher propensity for larger floods in watershed areas with a higher average annual rainfall and for those with relatively shallower soils. Since the depth-to-rock estimates are unavailable for most watersheds in PR [108], we simply relied on a depth-to-rock value of 111 cm (43.6 inches) as this was the average value of watersheds included in the original report [95]. Converted to Q pa units (m 3 s −1 km −2 ), these equations equal: where Q pa 100 and Q pa 500 refer to the 100-year and 500-year predicted area-normalized peak flow rate and R is the watershed-scale average annual rainfall with the low and high values being representative of relatively dry and wet portions of the island, respectively. Equations (7) and (8) relied on streamflow data collected only until 1994 (Ramos-Ginés, 1999). Our TC Q pa PR catalogue shows that between 1899 and 1994, there were 24 and 20 observations within the 100-year and 500-year peak flow range, respectively (Figure 8a,b). From 1994 to 2020, there were a total of 45 and 15 TC Q pa observations within and exceeding the 100-year and 500-year range predicted by Equations (7) and (8). These include observations during Hurricanes Hortense (n = 13 and 5 for Q pa 100 and Q pa 500 , respectively), Georges (n = 10 and 3), María (n = 10 and 5), and five other TCs (Kyle, Olga, Jeanne, Dean, and Grace; n = 5 and 0). In fact, the highest Q pa values in the entire record for watersheds larger than 100 km 2 have all been recorded since 1996. Therefore, the rather recent occurrence of flows with magnitudes matching these extreme recurrence interval values suggest the need for an updated version of these equations as it has been done in parts of the U.S. (e.g., Maurer, et al. [109]). This is particularly needed in PR because localized climate change projections suggest an increase in the frequency of the most extreme rainstorms [43]. In fact, studies claim that these projections already have become a reality as Hurricane María's intensity while in the vicinity of PR appears to have been enhanced by climate change [110]. The call for an update on PR's flood frequency curves echoes those that have already been made with regard to rainfall [111].

Conclusions
The highest area-normalized instantaneous peak streamflow maxima (Q pa ) recorded during tropical cyclones in Puerto Rico from 1899 to 2020 was relatively lower than the world's maximum envelope for watersheds smaller than 10 km 2 , yet within its span for those in the 10 to 619 km 2 range. The lack of much representation of watersheds < 10 km 2 in our catalogue likely led to limited opportunities to measure their most extreme conditions. A total of 20 Q pa values from watersheds in Puerto Rico ranging in drainage areas from 14 to 540 km 2 and induced by rainfall during eight different tropical cyclones were within the 95% confidence interval of the world's peak streamflow envelope. Some tropical cyclones with exceptionally high peak values in Puerto Rico include tropical depressions TD- 15 (1970) and the meteorological system that later developed into Hurricane Kyle (2008), tropical storm Eloísa (1975), and Hurricanes San Ciriaco (1899), Donna (1960), Hortense (1996), Georges (1998), and María (2017). Findings presented here support previous studies indicating that higher rainfall and peak streamflow values tend to be recorded near the central-eastern end of Puerto Rico. This is likely due to the predominant westward translation of most cyclones which tends to induce the most copious rainfall on the eastern end of the island as the storm interacts with the island's topography. However, the overall highest Q pa values in the catalogue were recorded throughout the entire island.
In terms of methods, we showed that quantile regressions hold promise in determining peak flow envelopes. This method is less parsimonious than ordinary least squares because the results of the latter method heavily rely on what observations are included in the regression. Should only the maximum peak flow for every arbitrarily selected watershed area bin be used or should more be included? These concerns are not as relevant in quantile regressions because the method is less sensitive to extreme values. Further work should compare both methods in terms of selection bias (i.e., what observations are included in the analysis) and overall performance in bounding the upper values of peak flow events.
The findings presented here highlight the importance of improving our quantitative understanding of peak flow rates observed in Puerto Rico, as well as throughout the Insular Caribbean. Peak rates observed in islands of the Caribbean have proven to be some of the most extreme values in the world for the size of the islands' watersheds and are likely to continue to display exceptional values due to projected regional manifestations of global climate change. Some of the most extreme events to have ever been recorded in Puerto Rico have occurred since the last update of the currently used flood frequency curves. Our findings therefore highlight the need for a revision given that these curves represent an essential component of flood hazard management for the island.

Conflicts of Interest:
The authors declare no conflict of interest.