Flood Impacts on Critical Infrastructure in a Coastal Floodplain in Western Puerto Rico during Hurricane Mar í a

: Flooding during extreme weather events damages critical infrastructure, property, and threatens lives. Hurricane Mar í a devastated Puerto Rico (PR) on 20 September 2017. Sixty-four deaths were directly attributable to the ﬂooding. This paper describes the development of a hydrologic model using the Gridded Surface Subsurface Hydrologic Analysis (GSSHA), capable of simulating ﬂood depth and extent for the Añasco coastal ﬂood plain in Western PR. The purpose of the study was to develop a numerical model to simulate ﬂooding from extreme weather events and to evaluate the impacts on critical infrastructure and communities; Hurricane Mar í a is used as a case study. GSSHA was calibrated for Irma, a Category 3 hurricane, which struck the northeastern corner of the island on 7 September 2017, two weeks before Hurricane Mar í a. The upper Añasco watershed was calibrated using United States Geological Survey (USGS) stream discharge data. The model was validated using a storm of similar magnitude on 11–13 December 2007. Owing to the damage sustained by PR’s WSR-88D weather radar during Hurricane Mar í a, rainfall was estimated in this study using the Weather Research Forecast (WRF) model. Flooding in the coastal ﬂoodplain during Hurricane Mar í a was simulated using three methods: (1) Use of observed discharge hydrograph from the upper watershed as an inﬂow boundary condition for the coastal ﬂoodplain area, along with the WRF rainfall in the coastal ﬂood plain; (2) Use of WRF rainfall to simulate runoff in the upper watershed and coastal ﬂood plain; and (3) Similar to approach (2), except the use of bias-corrected WRF rainfall. Flooding results were compared with forty-two values of ﬂood depth obtained during face-to-face interviews with residents of the affected communities. Impacts on critical infrastructure (water, electric, and public schools) were evaluated, assuming any structure exposed to 20 cm or more of ﬂooding would sustain damage. Calibration equations were also used to improve ﬂood depth estimates. Our model included the inﬂuence of storm surge, which we found to have a minimal effect on ﬂood depths within the study area. Water infrastructure was more severely impacted by ﬂooding than electrical infrastructure. From these ﬁndings, we conclude that the model developed in this study can be used with sufﬁcient accuracy to identify infrastructure affected by future ﬂooding events. interconnected during extreme the is into the lower Añasco A Geological Survey stream discharge gauge is available within the watershed, which is used to calibrate and validate the hydrologic parameters in this study. The stream discharge gauge is located in (USGS-50144000); no gauges are available in the watershed’s coastal ﬂood plain The annual average rainfall in is 2160


Introduction
Worldwide, floods affect more people than any other natural disaster, costing $104 billion annually on average [1]. Harmsen and Goyal [2] cited several historically significant floods that resulted in extreme loss of life, for example, 12,000 lives lost in China in floods in 1980, 1996, and 1998, 30, When Hurricane María hit Puerto Rico (PR) in September 2017, most large and small powerlines were knocked down. Eighty percent of the island's electrical and communication systems were destroyed [17,18]. Downed powerlines and landslides made roads impassable for weeks or even months in specific locations. The water utility was not prepared for the hurricane and lacked power generators for running pumps. Water intakes on rivers were damaged due to extensive river meandering and sedimentation. Wastewater treatment facilities were damaged, resulting in the release of raw sewage into the floodwaters. The dam at the Guajataca Reservoir in northwest PR was severely damaged, threatening flash flooding and leaving the population without water for many months. As a consequence of the devastation of Hurricane María, the Joint Operational Catastrophic Incident Plan of PR [19] was developed. The plan represented the input from thirty State agencies, four Federal agencies, and representatives of the private critical infrastructure sectors.
The threat of flooding has increased with the relentless advance of climate change. The number of record-breaking precipitation events has increased by 12 percent between 1981 and 2010 [14]. Georgescu et al. [20] used WRF to dynamically downscale future midand end-of-Century rainfall patterns across the United States, focusing on the effect of weather and rainfall on metropolitan areas.
The study concluded that urban areas could expect to receive more extreme weather events and consequent increased flooding risk in the future. In general, the Caribbean region has a well-documented warming trend [21], and it is expected to experience a drying climate in the future [22], combined with more extreme rainfall [4] and drought-related events. Recently, however, Jury [23], using high-resolution satellite and reanalysis products, has hypothesized that mountainous island environments may become wetter in the future, which could lead to greater flood risk in some Caribbean island watersheds. Flooding on Caribbean islands will become aggravated with increasing sea-level rise from a warming climate. Díaz et al. [22] reported that sea levels have been rising along the coast of PR at a rate of 2 mm year −1 during the 20th Century. However, during the 2000s, the rate has

Previous Hydrologic Studies Conducted in Puerto Rico
Several hydrologic studies have been conducted in the Añasco watershed, including [42][43][44][45][46]. Prieto [42] developed a surface water/groundwater MIKESHE model for the Añasco, Yagüez, and Guanajibo watersheds. The surface and subsurface conceptual and numerical models and GIS database resulted in a greater understanding of the complexities of these systems in western PR. Rojas Gonzalez [43] (also [2]) performed an upscaling experiment with the Vflo model [47] in PR to determine how the model grid spacing could be made coarser using the same parameter settings without loss of accuracy. She found, based on a stochastic analysis using the Generalized Likelihood Uncertainty Estimation (GLUE, [48]) methodology, that the 10 m hydrologic model resolution can be upscaled to 100 m resolution without loss of accuracy and that the 100 m rainfall resolution can be upscaled to 1000 m without loss of accuracy. Using the coarser grid resolution, improved model computation efficiency can be achieved. Torres Molina [44] (also [2]) evaluated shortterm rainfall estimates in western PR, derived from a high-resolution radar (TropiNet, [49]) and stochastic analysis within the Mayaguez Bay drainage basin. The model served as a potential prototype for predicting real-time flooding using the Vflo model. Ramos and Gilbes [45] developed a Soil Water Assessment Tool (SWAT) model for the Añasco watershed. The study's main objective was to simulate runoff and sediment yields within the watershed and sediment transport to the coastal waters from 1998 to 2012. Neither the hydrologic model nor the sediment transport model were calibrated. Nevertheless, the hydrologic model was able to capture 50 percent of the variability of the observed discharge. The study showed a clear connection between coffee farms and their practices and urbanization with large sediment plumes in the coastal waters near the mouth of the Añasco River. Nilawar et al. [46] estimated soil moisture in the Añasco and Guanajibo watersheds using the SWAT model, comparing the AMSR2 soil moisture product results. The SWAT model was calibrated for 2010-2012 using USGS stream gauge data from the two watersheds. The authors concluded that SWAT fusion with AMSR2 soil moisture data can simulate streamflow in the two watersheds.
Several other hydrologic modelling studies have been conducted in PR for conditions similar to those considered in this study, including [28,[50][51][52][53][54][55][56]. The authors of [50] developed a real-time runoff model to evaluate runoff for the Loíza Basin in northeast PR. Water from this basin is captured in the Carraizo reservoir and constitutes the primary water source for the San Juan metropolitan area. The model was based on estimating excess rainfall using the Green-Ampt infiltration equation. Hydraulic routing was accomplished using the HYDRAUX model [57], which solves the unsteady one-dimensional discharge in a network of open channels. The authors of [51] simulated streamflow for the Loíza Basin to characterize forecast uncertainties and statistical measures of model performance as part of a National Weather Service experiment in the Distributed Model Intercomparison Project [58]. The authors of [53] applied the Water Supply Stress Index (WaSSI) model [54] at the Luquillo National Rainforest in northeast PR. WaSSI estimates the monthly runoff volume and routed water exiting a watershed outlet. The authors of [52] applied the Hydrologic Simulation Program-FORTRAN (HSPF) to a Río Caonillas watershed in north-central PR. The three-year analysis goal was to better understand the hydrology, soil erosion, and sediment transport of the watershed. Another purpose of the study was to evaluate the ability of HSPF to perform under the PR's extreme conditions (e.g., rainfall intensities greater than 25 mm h −1 and average slopes of 38%). The model successfully captured 85% of the monthly streamflow variability and 75% of the monthly suspended sediment concentration variability.
The authors of [55] used the Global Land Data Assimilation System (GLDAS) Noah Land Surface model to evaluate differences in pre-and post-Hurricane María hydrology over PR. GLDAS simulates overland and subsurface runoff based on 3 h and monthly increments. The study revealed that the widespread removal of leaf cover by Hurricane María resulted in elevated sediment losses for months after the passing of the hurricane. The study also showed that Hurricane María caused a change in the atmosphere's thermodynamic structure and consequently affected cloud development and rainfall. As noted in the previous section, [28] developed a GSSHA model to simulate compound flooding in eastern PR during Hurricane George. The authors of [56] describe the Geostationary Operational Environmental Satellite Puerto Rico Water and Energy Balance (GOES-PRWEB) model. The operational model provides twenty-five hydro-agro-climate variables at a 1 km spatial resolution daily. Monthly and annual averages and totals are also available online (https://pragwater.com). The authors of [56] describe several applications of the model, including monitoring of the 2014-2016 drought in PR [59], evaluation of irrigation requirements for a 7-day period for a sweet pepper crop in southern PR, and an analysis of aquifer recharge in PR's southern coastal aquifer from 2009-2020. Soil saturation data from GOES-PRWEB was used to initialize soil saturation values in this study.

Study Objective
The objective of this study is to develop a methodology to quantify flood conditions from extreme weather events and their impacts on critical infrastructure. As a case study, a flood inundation model is applied to the Añasco coastal flood plain of western PR. The model simulates fluvial flooding and the effect of sea rise due to storm surge. Geographic Information System (GIS) methods are used to identify potentially damaged critical infrastructure due to flooding. This study is part of a more extensive NSF-funded study to establish a framework to evaluate and enhance critical infrastructure resiliency in islanded communities. Output from the model will be incorporated into a water/power network model to assess recovery scenarios under extreme weather events.

Case Study Description
This section presents a short description of the effects of Hurricane María on PR and some characteristics of the study area such as location, land use land cover, and soil type distribution, among others.

Hurricane María
Puerto Rico is located in a susceptible zone affected by hurricanes and tropical storms ( Figure 1). On 7 September 2017, Hurricane Irma hit northeast PR as a Category 3 hurricane. The eyewall passed just above the northeast corner of PR. Just 14 days later, Hurricane María passed over PR, entering the southeast corner and exiting from the northwest. Hurricane María was a Category 4 when it passed over PR on 20 September 2017 [60]. The maximum sustained wind speed was 69.3 m s −1 (249.5 km h −1 ) and the heavy precipitation produced serious flooding and mudslides. One of the locations in PR had a storm total of almost 38 inches (965 mm) [61]. The neutral El Niño Southern Oscillation (ENSO) and unusually high sea surface temperature (SSTs) in 2017 [62] contributed to the strength of Hurricane María. Glenn et al. [21] have observed a 0.015 • C year −1 rise in SSTs per decade in the Caribbean region between 1982 and 2012, suggesting that the conditions that enhanced Hurricane María will increase in the future. model simulates fluvial flooding and the effect of sea rise due to storm surge. Geographic Information System (GIS) methods are used to identify potentially damaged critical infrastructure due to flooding. This study is part of a more extensive NSF-funded study to establish a framework to evaluate and enhance critical infrastructure resiliency in islanded communities. Output from the model will be incorporated into a water/power network model to assess recovery scenarios under extreme weather events.

Case Study Description
This section presents a short description of the effects of Hurricane María on PR and some characteristics of the study area such as location, land use land cover, and soil type distribution, among others.

Hurricane María
Puerto Rico is located in a susceptible zone affected by hurricanes and tropical storms ( Figure 1). On 7 September 2017, Hurricane Irma hit northeast PR as a Category 3 hurricane. The eyewall passed just above the northeast corner of PR. Just 14 days later, Hurricane María passed over PR, entering the southeast corner and exiting from the northwest. Hurricane María was a Category 4 when it passed over PR on 20 September 2017 [60]. The maximum sustained wind speed was 69.3 m s −1 (249.5 km h −1 ) and the heavy precipitation produced serious flooding and mudslides. One of the locations in PR had a storm total of almost 38 inches (965 mm) [61]. The neutral El Niño Southern Oscillation (ENSO) and unusually high sea surface temperature (SSTs) in 2017 [62] contributed to the strength of Hurricane María. Glenn et al. [21] have observed a 0.015 °C yr −1 rise in SSTs per decade in the Caribbean region between 1982 and 2012, suggesting that the conditions that enhanced Hurricane María will increase in the future.
According to the National Weather Service, in 2017, there were 64 flood-related fatalities in PR, while only 2.3 on average from 2010 to 2020 (not including 2017). For comparison, Texas reported 70 flood-related fatalities during 2017 related to Hurricane Harvey. The total estimated Hurricane María damages in PR were $90 billion [60], making it the third most costly hurricane in U.S. history, after Hurricanes Katrina and Harvey.

Study Area
The study area is the Añasco watershed, which is located in Western PR within the municipalities of Añasco, Las Marías, Maricao, Lares, and Adjuntas. The watershed (Figure 2) is composed of a wide diversity of conditions, including urban, agriculture, and Harvey. The total estimated Hurricane María damages in PR were $90 billion [60], making it the third most costly hurricane in U.S. history, after Hurricanes Katrina and Harvey.

Study Area
The study area is the Añasco watershed, which is located in Western PR within the municipalities of Añasco, Las Marías, Maricao, Lares, and Adjuntas. The watershed ( Figure 2) is composed of a wide diversity of conditions, including urban, agriculture, and forest landscapes, topographic relief of over 1000 m, a river course with extreme meander, and a sizeable coastal flood plain, part of which includes marshland. There are four interconnected lakes located in the upper part of the watershed (Lagos Toro, Prieto, Guayo, and Yahuecas), which transfer water to southwest PR, except during extreme rainfall, when the water is released into the lower Añasco Watershed. A U.S. Geological Survey stream discharge gauge is available within the watershed, which is used to calibrate and validate the hydrologic parameters in this study. The stream discharge gauge is located in San Sabastian, PR (USGS-50144000); however, no gauges are available in the watershed's coastal flood plain portion. The annual average rainfall in this watershed is 2160 mm. Average, minimum, and maximum yearly air temperatures are 26.1, 18.9, and 37.2 • C, respectively.

°C, respectively.
In order to reproduce the flooding during Hurricane María in this study, it was necessary to calibrate and validate some of the hydrologic parameters of the watershed; therefore, the watershed was divided into two areas, as shown in Figure 2. Sub-watershed 2 represents the region above the USGS stream discharge gauge, which was used to calibrate and validate the hydrologic parameters (essentially roughness and infiltration parameters). Subsequently, Sub-watershed 1 was used to evaluate the flooding using the calibrated parameters obtained in Sub-watershed 2. Figure 3 shows the variation of the elevations within the Añasco watershed, ranging from 0 to 1203 m above sea level. The rainfall event selected to calibrate the model was Hurricane Irma which happened days before Hurricane María, and a similar rainfall event which occurred in 2007 was selected to validate the results (this will be explained in more detail in the methodology section). It is worth mentioning that the raster precipitation data for those two events were obtained from the Next Generation Radar (NEXRAD), known as Weather Surveillance Radar 1988 Doppler (WSR-88D), which is located at 18°06'56.16'' N−66°04'41.03" W. The star symbol in Figure 2 represents the location of the radar in PR.  In order to reproduce the flooding during Hurricane María in this study, it was necessary to calibrate and validate some of the hydrologic parameters of the watershed; therefore, the watershed was divided into two areas, as shown in Figure 2. Sub-watershed 2 represents the region above the USGS stream discharge gauge, which was used to calibrate and validate the hydrologic parameters (essentially roughness and infiltration parameters). Subsequently, Sub-watershed 1 was used to evaluate the flooding using the calibrated parameters obtained in Sub-watershed 2. Figure 3 shows the variation of the elevations within the Añasco watershed, ranging from 0 to 1203 m above sea level. The rainfall event selected to calibrate the model was Hurricane Irma which happened days before Hurricane María, and a similar rainfall event which occurred in 2007 was selected to validate the results (this will be explained in more detail in the methodology section). It is worth mentioning that the raster precipitation data for those two events were obtained from the Next Generation Radar (NEXRAD), known as Weather Surveillance Radar 1988 Doppler (WSR-88D), which is located at 18 • 06'56.16 N−66 • 04'41.03 W. The star symbol in Figure 2 represents the location of the radar in PR.

Land Use Land Cover and Soil Distribution
Land use land cover and soil dataset maps are essential for hydrologic simulations, allowing the association of land and soil properties into the model, such as roughness values, hydraulic conductivity, porosity, capillary head, and field capacity. Figure 4a shows the 30 m resolution land use land cover map for 2010 [63]. It was the most recent high-resolution land use land cover classification map developed for PR. The predominant land cover categories for the study area are Mixed Forest, Cultivated, and Low Intensity Developed Areas. On the other hand, Figure 4b shows the Soil Survey Geographic Dataset (SSURGO) over the Añasco watershed, which was acquired from the United States Department of Agriculture (USDA [64]). For this dataset, Clay and Clay Loam are the most representative soil categories in the Añasco watershed.

Methodology
The following chart ( Figure 5) summarizes the main steps used in this study:

Land Use Land Cover and Soil Distribution
Land use land cover and soil dataset maps are essential for hydrologic simulations, allowing the association of land and soil properties into the model, such as roughness values, hydraulic conductivity, porosity, capillary head, and field capacity. Figure 4a shows the 30 m resolution land use land cover map for 2010 [63]. It was the most recent highresolution land use land cover classification map developed for PR. The predominant land cover categories for the study area are Mixed Forest, Cultivated, and Low Intensity Developed Areas. On the other hand, Figure 4b shows the Soil Survey Geographic Dataset (SSURGO) over the Añasco watershed, which was acquired from the United States Department of Agriculture (USDA [64]). For this dataset, Clay and Clay Loam are the most representative soil categories in the Añasco watershed.

Land Use Land Cover and Soil Distribution
Land use land cover and soil dataset maps are essential for hydrologic simulations, allowing the association of land and soil properties into the model, such as roughness values, hydraulic conductivity, porosity, capillary head, and field capacity. Figure 4a shows the 30 m resolution land use land cover map for 2010 [63]. It was the most recent high-resolution land use land cover classification map developed for PR. The predominant land cover categories for the study area are Mixed Forest, Cultivated, and Low Intensity Developed Areas. On the other hand, Figure 4b shows the Soil Survey Geographic Dataset (SSURGO) over the Añasco watershed, which was acquired from the United States Department of Agriculture (USDA [64]). For this dataset, Clay and Clay Loam are the most representative soil categories in the Añasco watershed.

Methodology
The following chart ( Figure 5) summarizes the main steps used in this study:

Methodology
The following chart ( Figure 5) summarizes the main steps used in this study:

Hydrologic Model
The GSSHA is a process-based model, developed by the U.S. Army Engineer Research and Development Center (ERDC) Coastal and Hydraulics Laboratory (CHL), and is a multidimensional model that couples overland, surface, and subsurface flow. GSSHA simulates two-dimensional (2D) overland and groundwater flow, one dimension (1D) streamflow, and soil moisture [65]. A feature of GSSHA is the inclusion of spatial and temporal precipitation, allowing the use of NEXRAD precipitation products. It should be mentioned that GSSHA is approved by the Federal Emergency Management Agency (FEMA) for hydrologic modelling [66,67]. Some of the hydrological processes and approximations simulated in GSSHA are precipitation distribution, snowfall accumulation and melting, precipitation interception, overland water retention, infiltration, overland flow routing, evapotranspiration, soil moisture in the vadose zone, lateral groundwater flow, stream/ground water interaction, and exfiltration [68]. One of the key characteristics of GSSHA is the detailed modelling of the soil water profile in the unsaturated zone; this dynamic area partitions rainfall into infiltration, runoff, groundwater recharge and ET [68], and the GSSHA solves the one-dimensional (vertical direction) head-based form of Richards' equation [69], as follows: where C is the specific moisture capacity, is the soil capillary head (cm), z is the vertical coordinate (downward positive) (cm), t is the time (hr), ( ) is the effective hydraulic conductivity (cm hr −1 ), and W is the source/sink term (cm hr −1 ). It should be noted that the Richards' equation (RE) is highly nonlinear because C and K both depend on the soil water content; therefore, in GSSHA, the RE is linearized by making C and the intercell K constant during each time step [68]. In addition, GSSHA uses similar two-step, finite-volume schemes to route water for both 1D and 2D overland flow, where flows are computed based on heads and volumes and are updated based on the computed flows [68]. The intercell flows / and / (m 3 s −1 ) in the longitudinal x direction are computed from the depths d (m) at the n time level using the Manning equation for the head discharge relationship [68]: where n is the manning roughness coefficient, A is the cross-sectional flow area (m 2 ), R is the hydraulic radius (m), Sf is the friction slope, and So is the land surface slope in the longitudinal x direction. GSSHA uses the same method described for the 1D channel routing, but the calculation is made in two dimensions.

Calibration
The following steps describe the calibration process in GSSHA:

Hydrologic Model
The GSSHA is a process-based model, developed by the U.S. Army Engineer Research and Development Center (ERDC) Coastal and Hydraulics Laboratory (CHL), and is a multidimensional model that couples overland, surface, and subsurface flow. GSSHA simulates two-dimensional (2D) overland and groundwater flow, one dimension (1D) streamflow, and soil moisture [65]. A feature of GSSHA is the inclusion of spatial and temporal precipitation, allowing the use of NEXRAD precipitation products. It should be mentioned that GSSHA is approved by the Federal Emergency Management Agency (FEMA) for hydrologic modelling [66]. Some of the hydrological processes and approximations simulated in GSSHA are precipitation distribution, snowfall accumulation and melting, precipitation interception, overland water retention, infiltration, overland flow routing, evapotranspiration, soil moisture in the vadose zone, lateral groundwater flow, stream/ground water interaction, and exfiltration [67]. One of the key characteristics of GSSHA is the detailed modelling of the soil water profile in the unsaturated zone; this dynamic area partitions rainfall into infiltration, runoff, groundwater recharge and ET [67], and the GSSHA solves the one-dimensional (vertical direction) head-based form of Richards' equation [68], as follows: where C is the specific moisture capacity, ψ is the soil capillary head (cm), z is the vertical coordinate (downward positive) (cm), t is the time (hr), K(ψ) is the effective hydraulic conductivity (cm h −1 ), and W is the source/sink term (cm h −1 ). It should be noted that the Richards' equation (RE) is highly nonlinear because C and K both depend on the soil water content; therefore, in GSSHA, the RE is linearized by making C and the intercell K constant during each time step [67]. In addition, GSSHA uses similar two-step, finite-volume schemes to route water for both 1D and 2D overland flow, where flows are computed based on heads and volumes and are updated based on the computed flows [67]. The intercell flows Q n i−1/2 and Q n i+1/2 (m 3 s −1 ) in the longitudinal x direction are computed from the depths d (m) at the n time level using the Manning equation for the head discharge relationship [67]: where n is the manning roughness coefficient, A is the cross-sectional flow area (m 2 ), R is the hydraulic radius (m), S f is the friction slope, and S o is the land surface slope in the longitudinal x direction. GSSHA uses the same method described for the 1D channel routing, but the calculation is made in two dimensions.

Calibration
The following steps describe the calibration process in GSSHA: (1) Delineation of the watershed area above the USGS stream discharge station (Subwatershed 2): A Digital Elevation Model (DEM) of 10 m resolution was required to delineate the watershed and obtain the flow path accumulation. Subsequently, an outlet point was placed in the exact location of the USGS stream gauge (Latitude 18 • 17 03.00 , Longitude 67 • 03 03.46 NAD27) to obtain the delineated watershed above the USGS station. (2) Generation of the stream channel cross-section: The 3DEP Lidar Explorer tool was utilized from the USGS website to download the Light Detection and Ranging (LiDAR) data for the area of interest. The data acquired was from 2015 with a horizontal resolution of 100 cm and a vertical resolution of 10 cm. Subsequently, the data were processed using QGIS, and the cross-sections were created using QGIS and a python script. It should be mentioned that Lake Guayo collects water that is diverted to southwest PR through tunnels; however, water is released into the watershed when there are events with high rainfall intensity and duration, such as hurricanes or tropical storms. Thus, a downstream cross-section channel was used as the crosssection of the lake; this hypothesis assumes that the precipitation that falls into the lake area will be counted as a volume in the outlet point previously defined. Additionally, the smooth stream option was used to prevent adverse slopes. (3) Grid Selection: The grid spacing is essential to define the model's resolution; small grid cells lead to a high-resolution model, which requires high computational resources and increases computational time. A 100 m grid spacing was chosen, based on watershed size and the time required to run the simulation. (4) Land Use and Soil Data: The soil type and land use maps were added to the model using the datasets previously described. It should be noted from Tables A1 and A2 that the land and soil categories for Sub-watershed 1 ( Figure 2) and Sub-watershed 2 ( Figure 2) were predominantly Mixed Forest and clay, respectively. The percentages were calculated using a python script. On the other hand, the initial values of the roughness for each land class were assumed using the values in Table A1 [69]. In addition, Table A3 shows the values of the soil parameters required to use the Green Ampt method, such as hydraulic conductivity, capillary head, and porosity, among others [70]. (5) Soil Moisture Saturation: The soil moisture saturation product from GOES-PRWEB (pragwater.com, [56]) was employed to estimate the average initial soil moisture for the study area, which was determined to be 98.9%. Therefore, initial soil moisture saturation was assumed to be 98.9% of the soils' porosity values. (6) Precipitation Data: Hurricane Irma precipitation data was chosen to calibrate the model, representing an extreme weather event for PR. Therefore, NTP (Storm Total Precipitation) Level 3 product from NEXRAD was used to obtain the rainfall distribution from 6 September at 12:00 to 8 September at 11:00 of 2017 (AST). According to NOAA, the NTP product, which uses the Precipitation Processing System (PPS) algorithm, is useful for locating flood potential over urban or rural areas and estimating the total basin runoff [71]. In addition, a python code was implemented to convert the cumulative rainfall into hourly rainfall so that GSSHA can easily read it. (7) Calibration process: The Secant Levenberg-Marquardt (SLM) method is an efficient enhancement to the Levenberg-Marquardt (LM) method [72]. The Multilevel Single Linkage (MLSL) method was created to reduce the probability of not finding a local minimum. The SLM and the MLSL methods are implemented in GSSHA for conducting model calibrations. These two methods use the observed and the simulated hydrograph to find the objective function minimum, varying the parameters to be calibrated, where the objective function consists of the sum of weighted squared differences between modeled and transformed flow values, with all weights assigned a value of 1 [73]. The Box-Cox transformation with λ = 0.3 [74,75] was employed to transform the observed and modeled flows. Additionally, a time step of 10 s was used for the computational process, and overbank flow was allowed. The following parameters were selected to calibrate the model: • Roughness parameter for land use land cover: The Mixed Forest class was selected because, as shown in Table A1, this class is the most representative of the land classes for Sub-watershed 1 and Sub-watershed 2 of the Añasco watershed. To define the ranges, it was assumed a range of values between 0.1 and 0.16 which corresponds to areas dominated by trees generally greater than 5 m tall and greater than 20% of total vegetation cover. Neither deciduous nor evergreen species are greater than 75 percent of total tree cover [76]. It should be mentioned that the urban and suburban drainage systems were not explicitly modelled in this study. Overland runoff in these areas was modelled using roughness values associated with high and medium developed land use land covers, respectively. • Soil parameters: The hydraulic conductivity was identified in previous studies [28] as one of the most sensitive soil parameters. Therefore, it was decided to calibrate the hydraulic conductivity for the most representative soil type (clay) in the watershed.

•
Streambed roughness values: Two-channel bed roughness values were proposed for the Añasco watershed; channel roughness 1 corresponds to the streams of the basin and subbasins surrounding the mainstream. Channel 1 values were assumed to be between 0.02 and 0.05, based on an initial value used in previous studies [25]. On the other hand, channel roughness 2 was proposed for the upper area of Sub-watershed 2 with values ranging between 0.051 and 0.075, because this area is composed of shallow water with some gravel and rock and irregular cross-sectional channels, which generate significant friction [77].
The following table (Table 1) shows a summary of the parameters that were selected to calibrate the model.

Validation
For the validation process, a rainfall event was chosen with the criteria of having a similar or higher peak discharge in the magnitude of Hurricane Irma, similar initial soil saturation as Hurricane Irma, and similar volume of rainfall. After an exhaustive search using the USGS stream discharge gauge historical data, a rainfall event was found which occurred between 11-13 December 2007. The soil moisture could not be verified using the PRAGWATER website, which has a starting date of 2009. Therefore, the NOAA Advanced Hydrologic Prediction Service (AHPS) website was used to confirm that the previous days experienced significant rainfall, and therefore, we assumed that the soil was saturated. Table 2 shows a comparison of the spatial resolution, temporal resolution, maximum discharge recorded in the USGS stream discharge station, the volume from the observed hydrograph of the USGS stream discharge station, and the volume of the rainfall that fell in both of the sub-watersheds using the NEXRAD data. The performance of the models was evaluated by measuring the coefficient of determination (R 2 ), Nash-Sutcliffe efficient coefficient (NSE), and percentage of bias (PBIAS) between the simulated discharge at the outlet of the Sub-watershed 2 and measured discharge by the USGS stream gauge station. According to previous studies [78], a model can be evaluated as satisfactory for flow simulations if daily, monthly, or annual R 2 > 0.60, NSE > 0.5, and PBIAS ≤ ±15%, where: • The coefficient of determination (R 2 ) indicates the degree of correlation between the simulated and the observed data, describing the portion of the variance in the observed data explained by the simulation; the ranges are between 0 and 1, where 1 indicates less error in variance [79].
• The Nash-Sutcliffe efficiency coefficient (NSE) measures the perfect match between the predicted values and the observed values; the ranges are between −∞ to 1, where a value of 1 means a perfect match and values less than zero indicate unacceptable performance [79].
• Percent Bias (PBIAS) measures the average tendency of the simulated data to be larger or smaller than their observed counterparts, with the optimal value being 0; negative values indicate overestimation bias and positive values underestimation bias [80].
where, n is the number of observations, Y obs i is each of the observed data, Y sim i is each of the simulated data, Y obs is the mean observed values, and Y sim is the mean of the simulated values. Each of the metrics was used to evaluate the performance of the calibrated and validated hydrographs.

Flooding Event
The general methodology implemented to reproduce the flooding event during Hurricane María is described in the following steps: • Delineation of the watershed area below the USGS stream discharge station (Subwatershed 1): Using the same methodology as for Sub-watershed 2, a DEM of 10 m resolution was utilized to delineate the area below the USGS stream discharge station (Sub-watershed 1 of Figure 2). During this process, an outlet point was placed at the end of the Añasco watershed in the coastal zone to obtain a first approximation of the delineated area. Subsequently, the Arc tool from GSSHA was used to extend the delineated area, ensuring the total cover of Sub-watershed 1. This process is required for watershed delineation when the results from the automated delineation option from GSSHA do not match the actual delineated area. • The same procedure as described for Sub-watershed 2 was used to generate the crosssections of the stream channels for Sub-watershed 1 and define land use, soil type, and grid selection. It should be mentioned that the best-calibrated parameters obtained from Sub-watershed 2 were also used for Sub-watershed 1.

•
Definition of the inflow hydrograph and rainfall during Hurricane María: An inflow node was placed at the location of the USGS stream gauge in Sub-watershed 1 to simulate the water coming from Sub-watershed 2. Additionally, the hourly time series of the WRF simulated rainfall from Hurricane María was used as the precipitation data.

•
Definition of the storm surge: In the grid cells located along the coast, a variable stage (water surface elevation) boundary condition was created to simulate the storm surge during Hurricane María. For the simulations without using storm surge, a constant stage boundary condition was used.

Rainfall and Inflow Hydrograph Definition during Hurricane María
During Hurricane María, the WSR-88D radar located in PR was damaged due to extreme winds [60]; thus, level 3 NEXRAD data were unavailable. Therefore, the rainfall was defined using data from a concurrent study [61]; a Weather Research and Forecasting Model (WRF) was used to reproduce atmospheric variables during Hurricane María, such as precipitation and wind speed. The WRF precipitation output has a 1 km spatial and 1 h temporal resolution in this study, covering the period 09/19/2017, 12 pm to 09/21/2017, 11 am, a total of 48 h. The criteria used for choosing 48 h was to be consistent with the same number of hours used in the calibration and validation process. Figure 6a shows the cumulative rainfall from the simulated Hurricane María using the WRF model. In addition, 63 different rain gauges from the National Weather Service (NWS), which reported the cumulative rainfall for 48 h, were used to generate a rainfall distribution map using the inverse distance weighted (IDW) interpolation method with a spatial resolution of 700 m. Subsequently, this map was cropped to show the cumulative rainfall for the Añasco watershed, as shown in Figure 6b. A comparison made between these two cumulative maps indicates a difference in the volume of 65 million cubic meters, which translates to an error of approximately 35 % (Table 3).  Delineation of the watershed area below the USGS stream discharge station (Subwatershed 1): Using the same methodology as for Sub-watershed 2, a DEM of 10 m resolution was utilized to delineate the area below the USGS stream discharge station (Sub-watershed 1 of Figure 2). During this process, an outlet point was placed at the end of the Añasco watershed in the coastal zone to obtain a first approximation of the delineated area. Subsequently, the Arc tool from GSSHA was used to extend the delineated area, ensuring the total cover of Sub-watershed 1. This process is required for watershed delineation when the results from the automated delineation option from GSSHA do not match the actual delineated area.  The same procedure as described for Sub-watershed 2 was used to generate the crosssections of the stream channels for Sub-watershed 1 and define land use, soil type, and grid selection. It should be mentioned that the best-calibrated parameters obtained from Sub-watershed 2 were also used for Sub-watershed 1.  Definition of the inflow hydrograph and rainfall during Hurricane María: An inflow node was placed at the location of the USGS stream gauge in Sub-watershed 1 to simulate the water coming from Sub-watershed 2. Additionally, the hourly time series of the WRF simulated rainfall from Hurricane María was used as the precipitation data.  Definition of the storm surge: In the grid cells located along the coast, a variable stage (water surface elevation) boundary condition was created to simulate the storm surge during Hurricane María. For the simulations without using storm surge, a constant stage boundary condition was used.

Rainfall and Inflow Hydrograph Definition during Hurricane María
During Hurricane María, the WSR-88D radar located in PR was damaged due to extreme winds [60]; thus, level 3 NEXRAD data were unavailable. Therefore, the rainfall was defined using data from a concurrent study [61]; a Weather Research and Forecasting Model (WRF) was used to reproduce atmospheric variables during Hurricane María, such as precipitation and wind speed. The WRF precipitation output has a 1 km spatial and 1 hr temporal resolution in this study, covering the period 09/19/2017, 12 pm to 09/21/2017, 11 am, a total of 48 h. The criteria used for choosing 48 h was to be consistent with the same number of hours used in the calibration and validation process. Figure 6a shows the cumulative rainfall from the simulated Hurricane María using the WRF model. In addition, 63 different rain gauges from the National Weather Service (NWS), which reported the cumulative rainfall for 48 h, were used to generate a rainfall distribution map using the inverse distance weighted (IDW) interpolation method with a spatial resolution of 700 m. Subsequently, this map was cropped to show the cumulative rainfall for the Añasco watershed, as shown in Figure 6b. A comparison made between these two cumulative maps indicates a difference in the volume of 65 million cubic meters, which translates to an error of approximately 35 % (Table 3).    Figure 7 shows WRF precipitation rainfall, NWS cumulative rainfall, and the error in percentages between the WRF and NWS cumulative rainfall. A correction of the WRF time series was performed using the cumulative WRF rainfall and the NWS interpolated rainfall for the whole of PR; for this correction, it was assumed that the temporal distribution of the rainfall in the WRF model was correct. Therefore, the only variable to modify was the spatial distribution of the rainfall. To correct rainfall, Equation (7) was used to compute the distributed error between the two cumulative maps with respect to the WRF cumulative rainfall. Subsequently, Equation (8) was implemented to adjust the original hourly WRF precipitation values using the distributed errors previously described to obtain the corrected hourly WRF precipitation.
where, ε i,j is the distributed error between the two cumulative maps in position ith, jth, (t) is the hour WRF rainfall at t time, and R correct i,j (t) is the hour corrected rainfall at t time.   Figure 7 shows WRF precipitation rainfall, NWS cumulative rainfall, and the error in percentages between the WRF and NWS cumulative rainfall. A correction of the WRF time series was performed using the cumulative WRF rainfall and the NWS interpolated rainfall for the whole of PR; for this correction, it was assumed that the temporal distribution of the rainfall in the WRF model was correct. Therefore, the only variable to modify was the spatial distribution of the rainfall. To correct rainfall, Equation (7) was used to compute the distributed error between the two cumulative maps with respect to the WRF cumulative rainfall. Subsequently, Equation (8) was implemented to adjust the original hourly WRF precipitation values using the distributed errors previously described to obtain the corrected hourly WRF precipitation.
, ( ) = , ( ) + , ( ) * , where, , is the distributed error between the two cumulative maps in position ith, jth, , is the WRF cumulative rainfall, , is the NWS cumulative rainfall, , ( ) is the hour WRF rainfall at t time, and , ( ) is the hour corrected rainfall at t time. Three optional hydrographs were considered to recreate the flooding event during Hurricane María:


During Hurricane María, the USGS stream discharge gauge stopped working when it reported a maximum discharge measurement of almost 5200 m 3 s −1 , assumed to be the peak discharge for the whole event. Consequently, there is a gap of nearly 24 h where there is no data available, and after that, the USGS reported an estimated discharge. A simple linear regression was employed to interpolate the data to fill the Three optional hydrographs were considered to recreate the flooding event during Hurricane María:

•
During Hurricane María, the USGS stream discharge gauge stopped working when it reported a maximum discharge measurement of almost 5200 m 3 s −1 , assumed to be the peak discharge for the whole event. Consequently, there is a gap of nearly 24 h where there is no data available, and after that, the USGS reported an estimated discharge. A simple linear regression was employed to interpolate the data to fill the gap and obtain time series values every 15 min. Figure 8a shows the final built hydrograph superimposed on the USGS discharge data. The time series, limited to 48 h (Figure 8b), was used to simulate the inflow boundary condition at the upstream location in Sub-watershed 1 for two scenarios (defined below).

•
The calibrated parameters obtained in the calibration process were used for Hurricane María to obtain the output hydrograph for Sub-watershed 2 using the WRF rainfall and the corrected rainfall separately. Following this, each output hydrograph obtained was used as an inflow hydrograph at the upstream location for Sub-watershed 1 to reproduce the flooding event during Hurricane María.


The calibrated parameters obtained in the calibration process were used for Hurricane María to obtain the output hydrograph for Sub-watershed 2 using the WRF rainfall and the corrected rainfall separately. Following this, each output hydrograph obtained was used as an inflow hydrograph at the upstream location for Sub-watershed 1 to reproduce the flooding event during Hurricane María.

Storm Surge
According to the Hurricane María report developed by NOAA, the storm surge in Western PR reached values between 1 and 3 feet (0.3 and 0.9 m) [60]. A NOAA tide station (9759394) was identified close to the coast of the Añasco watershed, which contained data from predicted and partially observed water height levels during Hurricane María. It was noticed that the maximum height of the water level (using the Mean Lower Low Water, MLLW) was approximately 0.85 m. Subsequently, it was decided to merge the predicted and the observed data to obtain a height water level variation over time. The last point of the observed data was merged with the closest predicted water level height to create a smooth curve, as shown in Figure 9.

Storm Surge
According to the Hurricane María report developed by NOAA, the storm surge in Western PR reached values between 1 and 3 feet (0.3 and 0.9 m) [60]. A NOAA tide station (9759394) was identified close to the coast of the Añasco watershed, which contained data from predicted and partially observed water height levels during Hurricane María. It was noticed that the maximum height of the water level (using the Mean Lower Low Water, MLLW) was approximately 0.85 m. Subsequently, it was decided to merge the predicted and the observed data to obtain a height water level variation over time. The last point of the observed data was merged with the closest predicted water level height to create a smooth curve, as shown in Figure 9.
(9759394) was identified close to the coast of the Añasco watershed, which contained data from predicted and partially observed water height levels during Hurricane María. It was noticed that the maximum height of the water level (using the Mean Lower Low Water, MLLW) was approximately 0.85 m. Subsequently, it was decided to merge the predicted and the observed data to obtain a height water level variation over time. The last point of the observed data was merged with the closest predicted water level height to create a smooth curve, as shown in Figure 9.

Hydrologic Modelling Approaches (HMAs)
Six different approaches were created to evaluate and analyze the flooding results in the lower area of the Añasco watershed (Table 4).

Flood Depth Validation
To validate the results from the flooding model, a survey was carried out to obtain an estimate of the flooding depth from the residents around the area affected during the Hurricane María. It should be mentioned that the flooding extent could not be obtained due to the lack of satellite images during the event. Several trips were made to the study area to interview people to obtain information about their experience during Hurricane María. People were asked, based on their recollection, how high the floodwater rose. Typically, they would point to a position on a wall or indicate the depth relative to a part of their body, for example, their knee or their waist. Table A4 gives the survey results, including the latitude and longitude of the interview location and the depth of flooding. In the results section of this paper, the simulated flooding will be compared to the depths obtained from the face-to-face interviews provided in Table A4; the geolocation and the distribution of the data collected can be observed in Figure 10. A python script was created to transform the coordinate reference system from WGS 84 (ESPS:4326) to the same coordinate reference system used in the GSSHA output flooding map, WGS 84/UTM zone 19N. Subsequently, the same python script was used to obtain the flooding depth at each of the survey locations in the output flooding map from GSSHA, in order to compare the results. their body, for example, their knee or their waist. Table A4 gives the survey results, including the latitude and longitude of the interview location and the depth of flooding. In the results section of this paper, the simulated flooding will be compared to the depths obtained from the face-to-face interviews provided in Table A4; the geolocation and the distribution of the data collected can be observed in Figure 10. A python script was created to transform the coordinate reference system from WGS 84 (ESPS:4326) to the same coordinate reference system used in the GSSHA output flooding map, WGS 84/UTM zone 19N. Subsequently, the same python script was used to obtain the flooding depth at each of the survey locations in the output flooding map from GSSHA, in order to compare the results.

Methodology for Assessing Impacts to Critical Infrastructure
A python script was written to identify the infrastructure affected by the flooding. Input to the script included shapefiles with the geolocated infrastructure and each of the flooding maps. In this study, we assumed that a flood depth greater than 20 cm was defined as the criteria of failure for each of the water, electrical and public infrastructures. Table 5 summarizes all the infrastructures that were considered and their quantity within the study area.

Calibration Results
As was mentioned previously, two methods available in GSSHA were used to calibrate the model. The SLM calibration method took almost 72 h to finish compared to the MLSL method, which took more than 147 h. The final calibrated parameters are summarized in Table 6. Having output hydrographs for each of the calibrations and the observed hydrograph, GSSHA computed the performance evaluation criteria to indicate which models had the better results. As shown in Table 7, the SLM calibration method achieved all the evaluation criteria to be satisfactory. Unfortunately, the MLSL calibration method did not achieve two of the evaluation criteria. Therefore, the parameters obtained in the SLM calibration method were utilized to validate the 2007 rainfall event.  Figure 11a,b shows the comparison of the observed and modeled hydrograph during Hurricane Irma using the calibrated parameters from the SLM calibration method and the MLSL calibration method, respectively. Additionally, it can be observed from Table 7 that the PBIAS values are positive, which mean that the modeled hydrograph underestimates the amount of volume below the hydrograph.

Validation Results
In this section, the calibrated parameters obtained using the SLM calibration method were validated using a rainfall event from 11-13 December 2007. In this case, two simulations were made (Figure 12), using a time step of 2 s and a time step of 10 s. The 10 s time step was used first; however, there was a region in the hydrograph where an oscillation in the discharge occurred. For this reason, it was decided to decrease the time step to 2 s, and as a result, the oscillatory behavior was eliminated. In Table 8, it can be seen that all the evaluation criteria are satisfied, and the evaluation metrics improved when a time step of 2 s was employed. Moreover, it was evidenced how the computational time increased when a time step of 2 s was used compared to a time step of 10 s, from a 3.6 h runtime (time step of 10 s) to a 16 h runtime (time step of 2 s). In addition, it should be noted from the PBIAS value in Table 8 that the validation results produced an overestimation in com-

Validation Results
In this section, the calibrated parameters obtained using the SLM calibration method were validated using a rainfall event from 11-13 December 2007. In this case, two simulations were made (Figure 12), using a time step of 2 s and a time step of 10 s. The 10 s time step was used first; however, there was a region in the hydrograph where an oscillation in the discharge occurred. For this reason, it was decided to decrease the time step to 2 s, and as a result, the oscillatory behavior was eliminated. In Table 8, it can be seen that all the evaluation criteria are satisfied, and the evaluation metrics improved when a time step of 2 s was employed. Moreover, it was evidenced how the computational time increased when a time step of 2 s was used compared to a time step of 10 s, from a 3.6 h runtime (time step of 10 s) to a 16 h runtime (time step of 2 s). In addition, it should be noted from the PBIAS value in Table 8 that the validation results produced an overestimation in comparison to the measured discharge. In addition to the validation results, the calibrated model was used with the Hurricane María WRF rainfall (with and without correction) to obtain the hydrographs at the USGS stream discharge gauge location. These hydrographs were used in the next section to reproduce the flooding in the coastal flood plain. The hydrograph 1 and hydrograph 2 in Figure 13 show the output hydrographs obtained using the WRF rainfall and the corrected rainfall over Sub-watershed 2. Figure 13. Hydrograph 1 is the output hydrograph at USGS stream discharge location using the WRF rainfall without any correction (blue line); Hydrograph 2 is the output hydrograph at USGS stream discharge location using the WRF corrected rainfall (red line). Figure 14 show how the significant amount of discharge reported by the USGS stream discharge station affected the lower area of the Añasco watershed, especially along the coastal zone. The model produced a large extent and depth compared to the other HMAs (Figures 15 and 16). It should be noted that the incorporation of the storm surge in HMAs 2, 4, and 6 produced less than 1 m of flooding in the areas close to Añasco Bay. However, incorporating the storm surge did not produce considerable changes in the central part of Sub-watershed 1, where much of the flooding occurred. On the other hand, it can be noted in Figure 16 how the bias corrections are noticeable in the flooding map, producing a more significant extent of flooding than the HMAs of Figure  15, which represent the original result of WRF without any correction.  In addition to the validation results, the calibrated model was used with the Hurricane María WRF rainfall (with and without correction) to obtain the hydrographs at the USGS stream discharge gauge location. These hydrographs were used in the next section to reproduce the flooding in the coastal flood plain. The hydrograph 1 and hydrograph 2 in Figure 13 show the output hydrographs obtained using the WRF rainfall and the corrected rainfall over Sub-watershed 2. In addition to the validation results, the calibrated model was used with the Hurricane María WRF rainfall (with and without correction) to obtain the hydrographs at the USGS stream discharge gauge location. These hydrographs were used in the next section to reproduce the flooding in the coastal flood plain. The hydrograph 1 and hydrograph 2 in Figure 13 show the output hydrographs obtained using the WRF rainfall and the corrected rainfall over Sub-watershed 2. Figure 13. Hydrograph 1 is the output hydrograph at USGS stream discharge location using the WRF rainfall without any correction (blue line); Hydrograph 2 is the output hydrograph at USGS stream discharge location using the WRF corrected rainfall (red line). Figure 14 show how the significant amount of discharge reported by the USGS stream discharge station affected the lower area of the Añasco watershed, especially along the coastal zone. The model produced a large extent and depth compared to the other HMAs (Figures 15 and 16). It should be noted that the incorporation of the Figure 13. Hydrograph 1 is the output hydrograph at USGS stream discharge location using the WRF rainfall without any correction (blue line); Hydrograph 2 is the output hydrograph at USGS stream discharge location using the WRF corrected rainfall (red line). Figure 14 show how the significant amount of discharge reported by the USGS stream discharge station affected the lower area of the Añasco watershed, especially along the coastal zone. The model produced a large extent and depth compared to the other HMAs (Figures 15 and 16). It should be noted that the incorporation of the storm surge in HMAs 2, 4, and 6 produced less than 1 m of flooding in the areas close to Añasco Bay. However, incorporating the storm surge did not produce considerable changes in the central part of Sub-watershed 1, where much of the flooding occurred. On the other hand, it can be noted in Figure 16 how the bias corrections are noticeable in the flooding map, producing a more significant extent of flooding than the HMAs of Figure 15, which represent the original result of WRF without any correction.      This section presents the comparison between the flood depths obtained from the face-to-face interviews with the model-simulated flood depths for the six HMAs ( Figures  17-19). It is noted that there were no noticeable differences between the simulations with   This section presents the comparison between the flood depths obtained from the face-to-face interviews with the model-simulated flood depths for the six HMAs ( Figures  17-19). It is noted that there were no noticeable differences between the simulations with

Flood Validation Using Survey Data
This section presents the comparison between the flood depths obtained from the faceto-face interviews with the model-simulated flood depths for the six HMAs (Figures 17-19). It is noted that there were no noticeable differences between the simulations with and without storm surge. For example, HMAs 1 and 2 (Figure 17a,b) simulated identical conditions, the only difference being the storm surge in the HMA 2. In general, the simulated flood depths in HMAs 1 and 2 were larger than for the other HMAs. The linear regression equations and the R 2 and p values are provided in each of the graphs with the 95% confidence interval. It should be noted in Figures 17-19 how the p-values are less than the significance level (p-value ≤ 0.05), meaning that the data provide enough evidence to reject the null hypothesis, and that no correlation exists between the two datasets. The flood depth estimates can be improved by adjusting the flood depths using the linear regression equations. Figures 20-22 show the flood depth maps for HMAs 2, 4, and 6, corrected using the flooding results from Figures 14b, 15b and 16b Figure 18 the maximum flooding depth in Sub-watershed 1 using WRF rainfall for Hurricane María; t flow hydrograph for Sub-watershed 1 was obtained by simulating runoff using WRF rainfall in watershed 2. Additionally, the storm surge was incorporated (HMA 4).

Extent of Flooded Area
The lack of satellite images, such as Synthetic Aperture Radar (SAR), made it im sible to obtain details of the flood extent during Hurricane María. Therefore, the ar the 500-year annual chance of flood hazard map from FEMA was used to compar extent of flooding. Figure 23 shows the 500-year extent of flooding, which has an ar 34,789,370 m 2 .  Figure 19b) and the maximum flooding depth in Sub-watershed 1 using bias-corrected WRF rainfall for Hurricane María; the inflow hydrograph for Sub-watershed 1 was obtained by simulating runoff using biascorrected WRF rainfall in Sub-watershed 2. Additionally, the storm surge was incorporated (HMA 6).

Extent of Flooded Area
The lack of satellite images, such as Synthetic Aperture Radar (SAR), made it impossible to obtain details of the flood extent during Hurricane María. Therefore, the area of the 500-year annual chance of flood hazard map from FEMA was used to compare the extent of flooding. Figure 23 shows the 500-year extent of flooding, which has an area of 34,789,370 m 2 .  Table 9 shows the flooded area for each HMA and the percentage differences from the FEMA 500-year flood area ( Figure 23). The predicted flooded areas from all HMAs were less than FEMA's 500-year flood area. The lowest percentage difference was for HMA 2 corrected (−13.51%). In addition, the HMAs 5 and 6, where rainfall bias corrections were applied, had smaller differences compared to the results using the original WRF rainfall data (HMAs 3 and 4). We compared the areas from the different HMAs with the 500-year annual chance of flood hazard map, assuming that the character of Hurricane María corresponded to this kind of event and taking into account that it is the maximum return period available from FEMA.   Table 9 shows the flooded area for each HMA and the percentage differences from the FEMA 500-year flood area ( Figure 23). The predicted flooded areas from all HMAs were less than FEMA's 500-year flood area. The lowest percentage difference was for HMA 2 corrected (−13.51%). In addition, the HMAs 5 and 6, where rainfall bias corrections were applied, had smaller differences compared to the results using the original WRF rainfall data (HMAs 3 and 4). We compared the areas from the different HMAs with the 500-year annual chance of flood hazard map, assuming that the character of Hurricane María corresponded to this kind of event and taking into account that it is the maximum return period available from FEMA.

Critical Infrastructure
In this section, different water, electrical, and public infrastructure locations were used to evaluate the flooding impacts, assuming infrastructure failure when the flood depth is greater than 20 cm. In this analysis the simulated flood depths were corrected using the regression equations presented in Figures 17b, 18b and 19b. Although the regression analyses did not result in R 2 values close to 1, the application of the regression equations will significantly improve the flood depth estimation. Tables 10 and 11 show the number of water and public/electrical infrastructures impacted by flooding, respectively. In Table 10, it is evident how some of the potable and sewer water components are potentially affected by the flooding caused by Hurricane María; specifically, sewer pump stations, potable water pump stations, and some of the wells owned and managed by the water authority (Autoridad de Acueductos y Alcantarillado de PR [AAA]). Additionally, it can be shown in Table 10 that most of the sewer pump stations are affected by flooding in HMAs 2 and 6. On the other hand, four potable water pump locations were impacted by the flooding event during Hurricane María in all the HMAs. Furthermore, some wells not located in the flood-prone areas were affected by HMAs 5 and 6. The flooding depth increased in some other areas due to the increment of the bias-corrected WRF rainfall. These failures are critical, especially the potable water components because they cause the disruption of potable water in the system. For example, pump stations are needed to fill water in the system tanks, deliver water to higher location customers, or increase the flow rate to achieve customer requirements.  A flood analysis of some other public and electrical infrastructure was conducted (Table 11). From the evaluation, we determined that there were no electrical substations or transmission centers affected by the flooding event during Hurricane María for any of the HMAs. On the other hand, Table 11 shows a public school was affected by Hurricane María flooding, occurring in four of the six HMAs. Tables 12 and 13 show the number of water and public/electrical infrastructures impacted, respectively, using the linear regression approximation for flooding. With the corrected flooding, it is noted (Table 12) how the quantity of affected infrastructure increases, especially the sewer and potable pump stations in HMA 2 (relative to Table 10); additionally, a number of wells were impacted. On the other hand, electrical infrastructure shows no damage in the transmission center and substations (Table 13), while the number of public schools affected by flooding increased.

Discussion
Based on the evaluation criteria, the validation results can be considered very good, with R 2 = 0.946, NSE = 0.903, PBIAS (%) = −1.87, and volume error = 1.9%. The model was able to successfully capture 95% of the variation in the streamflow variability. On the other hand, the calibration results were not as good, indicating that the model could only capture 61% of the variation in the streamflow variability and have a 15.2% error in the volume estimate. Only one other study calibrated hourly stream discharge for the Añasco River. Using the Vflo model, [27] adjusted model parameters for eight different storms. The storms that were evaluated were significantly smaller than the storms we considered in this study. Although they did not use the same evaluation metrics that we used, their model performance appears to be comparable. All other hydrologic model studies for the Añasco watershed simulated daily or monthly average discharge. The authors of [26,28,29] reported coefficients of determination (R 2 ) of 0.24, 0.50, and 0.57, respectively. The authors of [28] did not perform a formal model calibration but rather compared their model, configured with best-estimated model parameters, with available stream gauge data. Although reporting on a model calibration, [25] did not provide calibration statistics.
There are various sources of uncertainty that may explain the lower metrics for Hurricane Irma calibration, including:

•
The rating curve used by the USGS to convert river stage to discharge may not have been accurate during the Hurricane Irma event.

•
Although an effort was made to accurately account for the water released from the lakes in the upper Añasco watershed, the PR Energy and Power Authority (PREPA), the agency that manages the lakes, could not provide any information on discharges from the lakes. • NEXRAD radar rainfall data may also introduce errors into the hydrologic simulation. Rojas Gonzalez [26] reported a simulation of streamflow for the Añasco River in which NEXRAD significantly overestimated the rainfall. A known problem associated with radar rainfall in Western PR is related to the radar's inability to see low clouds (less than 3000 m) at large distances from the radar facility due to the earth's curvature [27]. This problem would be expected to produce underestimates in rainfall. However, one of the authors of this paper has considerable experience using radar rainfall in Western PR and has observed both under and overestimation compared with NWS rain gauges.
Time series rainfall from rain gauges was not available for Hurricane María. Consequently, it was necessary to use a WRF model rainfall dataset. In general, the WRF model of [61] did a good job of simulating rainfall over PR, with a Normalized Root Mean Square Error (RMSE) of 0.2. However, we observed inconsistencies in the total rainfall depths within our study area compared with the NWS interpolated rain gauge data. A correction of the WRF time series rainfall was made using the NWS rain gauge data. There was little observable difference in flooding between the corrected and noncorrected rainfall (Figures 15 and 16).
Neither the WRF nor the bias-corrected WRF rainfall datasets could produce the large volume measured by the stream gauge during Hurricane María (Figure 8), and the use of the USGS stream gauge data for the inflow boundary condition to Sub-watershed 1 resulted in significantly greater flooding ( Figure 14) than resulted from the WRF noncorrected rainfall (Figure 15) or the corrected rainfall ( Figure 16).
The authors of [28] simulated the combined effect of storm surge and overland runoff in Eastern PR for Hurricane George. Although the peak of the storm surge, estimated to be 3 m, did not correspond with the peak runoff, the compound effect did aggravate flooding. In this study, we also used GSSHA to simulate the combined effects of storm surge and runoff. The authors of [60] estimated the storm surge in the study area to be 0.3 to 1 m. Figures 14-16 show flooding within and without storm surge included. The figures indicate that the impact from the storm surge was minimal. However, it should be noted that individuals interviewed in the study living near the coast observed the storm surge. Figures 17-19 show the Hurricane María flood depths obtained from face-to-face surveys versus the simulated flood depths for the six simulated HMAs. The data pairs associated with data points Nos. 17 and 32 were considered outliers and removed from the dataset. Datapoint No. 17 was removed because there was another data point within the same model grid pixel (No. 18) with an observed depth three times that of No. 17.
There was considerable small-scale topographic variation, which the model could not capture with the 100 m grid spacing. Datapoint No. 32 was located at a car dealership. The grid pixel contained two large buildings with roof heights of approximately 10 m and a raised highway approximately 5 m above the car dealership parking lot. The complexity of the topography at the location may have introduced an error in the grid elevation. It should also be noted that the survey data was acquired from resident interviews and not through post-flooding analysis which introduces some uncertainties in the surveyed depth; additionally, the flooding extent could not be validated because of the lack of satellite images during the Hurricane María event.
In general, the GSSHA model did not provide highly accurate flood depths, as evidenced from the low coefficients of determination (R 2 ) values of 0.547, 0.399, and 0.441 for HMAs 1 and 2, 3 and 4, and 5 and 6, respectively (note that the r 2 values were the same with or without storm surge due to the negligible effect storm surge had on flooding at the locations of the comparisons). To improve the model's ability to estimate flood depths, the linear regression equations were applied to enhance the depth estimates. Finally, a GIS intersection analysis was performed to evaluate infrastructure impacts, assuming to occur for flood depths greater than 20 cm. The analysis revealed numerous water components that were potentially impacted. No electrical infrastructure was identified as being impacted, while several public schools were potentially impacted. Electric infrastructure is known to have been damaged in the study area during Hurricane María; however, the cause is likely attributable to wind damage.

Model Limitations
Some of the limitations of the model are described:

•
The model used in this study was developed for large storm events proceeded by saturated soil conditions. Consequently, the model has not been calibrated for conditions when soil saturation is low and infiltration becomes a significant factor. Future efforts should consider calibrating the average matric suction at the wetting front during periods of low soil saturation. • The hydrologic model of Sub-watershed 1 was not calibrated in this study. However, the depths and extent of the flood were compared with data obtained from a homeowner survey. Although the analysis did not yield a high correlation between the observed and simulated flood depths, we applied the resulting regression equations, which significantly improved the performance of the flood depth estimation. The flood depths obtained from the face-to-face interviews were limited to flood plain areas where people were present during Hurricane María.

•
The model was based on a simplified conceptualization of the land surface parameters. Single values of roughness and hydraulic conductivity were used for the dominant land use/soil in Sub-watershed 2. More accurate distribution of these parameters in the watershed may help to improve the model calibrations.

Future Work
Recommendations for future research include: • Couple the flood model with a WRF rainfall forecast model to obtain 2-to 3-day advanced flooding warnings.

•
Couple the flood modelling results with a Power/Water network model. • Expand the model domain to include the Culebrinas, Yagüez, and Guanajibo watersheds. This will expand the area of coverage from 350 km 2 to over 900 km 2 . • Continue to work with the water and electric authorities to refine our approach for identifying impacted infrastructure. • Use the model to investigate ways to mitigate flooding in Sub-watershed 1. For example, evaluate changes in agricultural land-use practices or incorporate natural infrastructure.

Conclusions
A GSSHA model was developed in this study to simulate overland runoff, streamflow, and coastal storm surge in Western PR during Hurricane María. The upper gauged watershed (Sub-watershed 2) was calibrated with stream discharge data for Hurricane Irma and validated using a comparable storm from 11-13 December 2007. The Secant Levenberg-Marquardt (SLM) method was used for model calibration. Because the NEXRAD radar was damaged during the hurricane, we used the rainfall generated from the WRF model. Additionally, a bias-corrected WRF rainfall data set was used. The bias corrections were derived from the NWS total rainfall dataset. Flooding was simulated in the coastal flood plain (Sub-watershed 1) for six different HMAs. Storm surge did not appear to aggravate the flooding, at least in the areas considered in this study. Although the simulated flooding depth did not correspond exactly with the observed data, a set of calibration equations was employed to improve the model's performance. An analysis was conducted to determine the impact of flooding on water, electricity, and public infrastructure. Numerous water infrastructure components were impacted, including sewer pump stations, potable water pump stations, and wells. No electrical infrastructure sustained damage from flooding (per our criteria: water depth > 20 cm), while several public schools were damaged by flooding. Future research will focus on expanding the model domain and coupling the flood model with a rainfall prediction model to produce a flood forecast alarm tool.   Table A1. Percentage distribution of the land classes within Sub-watershed 1 and Sub-watershed 2 from Figure 2, and the Manning's values for each of the classes [67].