Hurricane Maria in the U.S. Caribbean: Disturbance Forces, Variation of Effects, and Implications for Future Storms

Abstract

The impact of Hurricane Maria on the U.S. Caribbean was used to study the causes of remotely-sensed spatial variation in the effects of (1) vegetation index loss and (2) landslide occurrence. The vegetation index is a measure of canopy ‘greenness’, a combination of leaf chlorophyll, leaf area, canopy cover and structure. A generalized linear model was made for each kind of effect, using idealized maps of the hurricane forces, along with three landscape characteristics that were significantly associated. In each model, one of these characteristics was forest fragmentation, and another was a measure of disturbance-propensity. For the greenness loss model, the hurricane force was wind, the disturbance-propensity measure was initial greenness, and the third landscape characteristic was fraction forest cover. For the landslide occurrence model, the hurricane force was rain, the disturbance-propensity measure was amount of land slope, and the third landscape characteristic was soil clay content. The model of greenness loss had a pseudo R² of 0.73 and showed the U.S. Caribbean lost 31% of its initial greenness from the hurricane, with 51% lost from the initial in the Luquillo Experimental Forest (LEF) from Hurricane Maria along with Hurricane Irma. More greenness disturbance was seen in areas with less wind sheltering, higher elevation and topographic sides. The model of landslide occurrence had a pseudo R² of 0.53 and showed the U.S. Caribbean had 34% of its area and 52% of the LEF area with a landslide density of at least one in 1 km² from Hurricane Maria. Four experiments with parameters from previous storms of wind speed, storm duration, rainfall, and forest structure over the same storm path and topographic landscape were run as examples of possible future scenarios. While intensity of the storm makes by far the largest scenario difference, forest fragmentation makes a sizable difference especially in vulnerable areas of high clay content or high wind susceptibility. This study showed the utility of simple hurricane force calculations connected with landscape characteristics and remote-sensing data to determine forest susceptibility to hurricane effects.

1. Introduction

It has been estimated that when Hurricane Maria hit Puerto Rico on 20 September 2017, 23 to 31 million trees were severely damaged or killed [1]. Hurricane Maria was a category 4 storm on the Saffir-Simpson scale. Its path went diagonally across the island, maximizing the area of impact. Two weeks earlier the larger category 5 hurricane, Irma, passed 90 km to the northeast of the island. It is accepted that hurricanes cause massive vegetation disturbance and are the primary driver of forest transformation in the hurricane-belt of the tropics [2]. Hurricanes and their effects have been widely studied because they are a major cause of loss of life and property [3]. Studies of hurricane effects on forests and ecosystem function overwhelmingly look at specific ecosystem components or are syntheses of such studies. These run the gamut from short- and long-term effects of wind and rain to secondary effects of altered post-hurricane conditions [4,5]. The spatial variation in hurricane effects are widely observed due to the obvious correlation of effects with a hurricane’s path and intensity. However, there is variation beyond the distance from the hurricane eye, and few studies have looked at the causes of this variation. The determinants of forest susceptibility to hurricane effects still lack landscape-wide characterizations that capture within-landscape patterns of effects related to hurricane force intensity, topography, and forest structure.

Recently, availability of satellite data has made these landscape-wide studies effects more quantifiable, illuminating patterns to assess cause and effect of hurricane disturbance [6]. There are landscape-scale studies on forest susceptibility to hurricane (or large storm) wind forces causing vegetation loss. These studies used disturbance patterns from specific storms to draw conclusions. Both Boose et al. [2,7] and Negrón-Juárez et al. [8] assessed the role of topography and the resulting local increases in wind strength. Schwartz et al. [6] assessed the role of forest fragmentation. Wang and Xu [9] studied the role of topography, forest density, and soil stability. Hu and Smith [10] compared the role of topography and distance from the hurricane eye versus landform. Landscape-scale studies on forest susceptibility to hurricane rain forces causing landslides are more numerous but most are in the form of regional or global susceptibility maps for landslides, not studies of the effects of specific storms. One exception is the work of Philpott et al. [11], which compared the role of topography and landform (using proxy of river proximity) versus vegetation density and complexity. Landslide susceptibility maps, e.g., [12,13,14], often do not consider rain amounts, but instead compare locations of previous landslides for commonality in topography, forest density, soil, landform, and road proximity. Most of these studies could not directly compare how much of the hurricane disturbance variation could be attributed to the local hurricane force magnitude variation versus the land characteristic variation, as only three studies [2,7,8] had wind maps available and none had maps of rain. No study could be found that assessed both causes of vegetation loss and landslide occurrence; assessment on the same landscape with the same methods would allow a more direct comparison of the nature of forest vulnerabilities to the different hurricane forces of wind and rain.

Projections of increasingly warm atmosphere and oceans indicate increased likelihood of more intense wind and rain associated with stronger hurricanes [15,16,17], with an expected increased disturbance of forest systems. Thus, better understanding of within-landscape variation in hurricane effects can lead to better planning of activities such as infrastructure development, land management and reforestation efforts that are more resistant and resilient to hurricanes. The severity with which Hurricanes Irma and Maria struck the U.S. Caribbean provides an exceptional case study of the effects of large hurricane energy expended on a tropical ecosystem. Assessing these modeled effects with parameters from previous storms of wind speed, storm duration, rainfall, and forest structure over the same storm path and topographic landscape provides examples for predicted future scenarios.

This study, using broad-scale maps of vegetation index loss and landslides and idealized maps of wind and rain, seeks to estimate through generalized linear models how much of the effect in the U.S. Caribbean from Hurricane Maria can be attributed to innate landscape characteristics, and how much can be attributed to forest fragmentation and the distribution of trees across the landscape. Based on these attributions, and using forcings of potentially stronger, slower, and wetter future storms, this study makes estimations on the immutable vulnerabilities of future forests, as well as vulnerabilities that can be managed by with land resource decisions. Novel to this study is a demonstration of the reproducibility of the remotely-sensed spatial pattern of hurricane effects with fairly simple models, specifically in recreating the hurricane wind force. It is demonstrated here that remotely-sensed data used along with physics-based calculations and field data can characterize within-landscape variations of vegetation susceptibility to hurricane effects.

2. Materials and Methods

2.1. Study Area

The U.S. Caribbean (Puerto Rico and the U.S. Virgin Islands) has a land area of 9238 km², and is located in the Northern Caribbean (18°N, 66°W). Four mountain ranges cover the main island of Puerto Rico, with elevations up to 1338 m. The Luquillo Mountain range is in the northeast, and the three ranges making up the west-east spread of the Cordillera Central dominate the landscape of the middle two-thirds of the island. Erosion of these mountains led to coastal plains on the north and south, with valleys running perpendicular to the coasts, especially on the north side of the island. The bottoms of the valleys and the coastal 10–20 km are very flat. The lesser islands of Puerto Rico and the U.S. Virgin Islands also have mountainous terrain. Soil in the mountains is a result of the weathering of the volcaniclastic rock, creating clay, whereas soil on the coastal plains is mostly sand from sedimentary sources [18]. The steep elevational gradients localize weather and climate that forest ecosystems subsist in. These ecosystems have been challenged in their existence by a history of considerable land cover change [18,19]. Prior to European settlement in the 16th century, the region was likely 95% forested with very dense vegetation [20]. Then, agricultural development on the main island of Puerto Rico led to a low point on 6% forested land in 1940. A switch to an industrial economy meant land was unmanaged or protected, so the forested land was back up to 42% by the year 2000. However, most of this reforestation happened around the forest edges, with more deforestation still happening on the forest interiors [21].

Natural disturbance further complicates the forest-climate interaction: the frequent hurricanes in the Atlantic are often at their maximum intensity in this location [7]. This region is ideal to study landscape-wide characterizations of forest susceptibility to hurricane effects because of the diversity of landscapes on the island, with landforms (i.e., plain, valley, wetland, mountain top or side) and the characteristic forest structure and topography that go with these landforms changing over short distances. This means that similar hurricane forces, expended as the storm moves across the region, are enacted on a diversity of landscapes.

2.2. Effects Maps

Vegetation indices are a transformations of remotely sensed spectral bands used to measure canopy ‘greenness’, a combination of leaf chlorophyll, leaf area, canopy cover and structure [22]. The Enhanced Vegetation Index (EVI) has been shown to be more responsive to canopy structural variations than other indices, and work well for measuring density of green vegetation in high biomass areas such as the tropics [22]. Values are generally in the range of 0.2 to 0.8 for healthy vegetation, with higher values indicating more greenness [22]. An effect map (Figure 1a) of the greenness loss was computed from Moderate Resolution Imaging Spectroradiometer (MODIS) satellite imagery [23].

EVI data at 250 m resolution with 16-day compositing period 8 days apart from the Aqua and Terra satellites were used here to calculate greenness lost from Hurricane Maria, with contribution from Hurricane Irma. Low-quality and cloud-influenced pixels were removed. Due to the high frequency of cloud cover in the season leading up to the hurricane and the relative stability of the vegetation during this time, four months of data (14 MODIS scenes, 1 May through 21 August) were averaged to compute the pre-hurricane EVI (Figure 2e). The 14 scenes were all before Hurricane Irma passed near the island (6 September 2017). Spatial gaps in the data were filled in with the ‘quiet year’ average EVI value for May 1 through August 21 for the particular 250 m cell location. A quiet year is defined as a year without significant disturbance, so no tropical storms or droughts in the U.S. Caribbean; there are seven of these quiet years in the MODIS collection years of 2002 to 2017. The post-hurricane EVI was calculated from the average of one month of data after the passage of Hurricane Maria (four MODIS scenes, 22 September through 31 October). Re-greening of vegetation started happened within days of the hurricane passing, so only a small number of scenes were used. The EVI lost from the hurricane was computed as the post-hurricane EVI subtracted from the pre-hurricane EVI, with the average quiet year (spatially computed) difference from 1 May–21 August to 22 September–31 October removed. Note, this means if there is usually vegetation greening from May to October in a location, that location could lose EVI from the hurricane and from the lost growing season, so it could lose more EVI than it started with pre-hurricane.

An effect map of landslides (Figure 1b) was taken directly from previous research [24]. Landslides were estimated from aerial imagery, defined on a 2 km grid as no landslides, <25 landslides per 1 km², or >25 landslides per 1 km². Landslides were identified by areas of exposed rock or soil, which could be easily seen in imagery because so much of the vegetation had lost all leaves [24]. Landslides were assumed to have occurred as a result of Hurricane Maria, although the authors concede it is possible they occurred from Hurricane Irma or other rainstorms. This study uses the data in a binary form, as landslide occurrence or no occurrence, and assumes all landslides can be attributed to heavy rain from Hurricane Maria.

2.3. Hurricane Forces

In this study, the hurricane forces of wind and rain are treated as independent variables that interact with forest canopy and soil, respectively, to cause disturbance [25]. Wind causes disturbance through in situ tree falls, branch breakage, and defoliation, while rain causes disturbance through landslide mass transfer (including vegetation). The hurricane force of wind (Figure 2a) was computed on each MODIS 250 m cell as total gale-force wind kinetic energy of Hurricane Maria and Irma over Puerto Rico and the U.S. Virgin Islands [2]. The kinetic energy was modeled from topographically sheltered continuous wind fields, with wind from forecasts by the National Weather Service (NWS) in 3-h time steps of the hurricane progression (10 forecasts for Maria and five for Irma across the region of interest). These forecasts were updated as the hurricanes progressed and are the closest available data to observational data (personal communication with Ernesto Rodriguez, NWS; 14 March 2018).

First, the forecasts were interpolated using theoretical wind profiles to 10-min wind fields, which is approximately continuous on the MODIS 250 m grid. The wind profiles (wind velocity with distance from the eye center) area approximated assuming air in cyclostrophic balance [26] (Equation (5)) following the HURRECON model equations [7]. The HURRECON model [2,7] further assumes the hurricane eye and surrounding isobars are circular with asymmetrical wind speeds (higher to the right of the eye and lower to the left) due to forward velocity of the storm (assumptions of the model are discussed in [2]). Friction forces and cross-isobar inflow angle for wind over land, versus over ocean, were accounted for with parameters taken from HURRECON [7] and other parameters are fit to the NOAA forecasts. A double eye wall was modeled for Hurricane Irma, following the forecasts.

Next, the continuous wind fields were modified for topographic sheltering [27] on a 125 m sheltering distance (half of a MODIS cell). Average elevation, slope, and aspect of each MODIS cell were summarized from the U.S. Geological Survey National Elevation Database (30 m resolution). Using these values, wind direction was first deflected from the HURRECON values [7] using the slope of the downwind terrain, and then sheltered with upwind slope (need to consider aspect and slope of the MODIS cell) [27] (Equations (15) and (16)). Higher elevations have linearly less shelter due to air mass characteristics [27] (Equation (17)). These winds were then used to compute kinetic energy of gale winds per m³ of air over the land, or the sum of half of velocity squared for the time of hurricane passage using all winds over 17 m/s.

The hurricane force of rain (Figure 2b) was computed on each 2 km landslide grid cell using station-observed rain totals from the NWS for Hurricane Maria only. There were 42 station observations at a set of wide-ranging elevations; these data were made continuous spatial field by the NWS with natural neighbor interpolation method (personal communication with José Alamo, NWS; 2 May 2018). Natural neighbor methods gave a more conservative spatial field that was more resistant to parameter changes, than did co-kriging with elevation as a covariate. This would be expected in an area with many different local topography orographic influences [28]. This rainfall map is a coarse resolution map without the resolution needed to characterize the intricacies of the hurricane impact on the ecological system, but instead with a resolution that can be used to model the effects of the rain force on landslide density.

2.4. Landscape Characteristics

Landscape characteristics were summarized to the resolution of the effects maps. Characteristics used in the models are shown in Figure 2c–g. Percent clay content in the depth 0–15 cm (Figure 2d) was provided by the Natural Resources Conservation Service, (NRCS) (personal communication with Manuel Matos NRCS; 2 May 2018). Fraction forest cover in an area (Figure 2f) was computed from the percent of woody vegetation from the PRGAP land-cover map, a 15 m raster [29].

Increased vegetation loss and landslide occurrence from hurricanes in fragmented forests has long been suggested and contested with localized studies [11,30,31]. Fragmented forest are thought to experience more extreme hurricane wind forces due to isolation [30] and stronger effects of rain because of more roads and agriculture diminishing vegetation complexity [11,32]. Fraction forest edge (Figure 2g) was used as a measure of forest fragmentation, and was computed with morphological spatial-pattern analysis on the PRGAP land-cover map [33]. Other fragmentation measures were investigated (fraction holes in forest and fraction core forest), but these were found to be highly correlated with fraction forest and slope, thus not useable in a generalized linear model along in conjunction with the use of forest cover or slope (slope is shown in Figure 2c). Landform type (Figure 2h; wetland, valley, top, side, or plain) did not influence model results and was excluded from all models.

Note, characteristics of aspect, slope, and elevation are contained in the wind energy model (as discussed above) and elevation is contained (implicitly) in the rain natural-neighbor surface. Thus, models using a particular force will not use those topographical characteristics separately from the force.

2.5. Generalized Linear Models

Two generalized linear models (GLMs) were made with independent (not significantly correlated with each other) variables to model the dependent variable of effect (greenness loss and landslide occurrence). In a GLM, the dependent variable is assumed to come from a distribution in the exponential family with a link function of the expected value equaling a linear combination of the dependent variables. All variables were scaled to be of order one. Models were built by adding and removing variables to find the smallest ‘best’ model. ‘Best’ here is used in the sense that the addition of a variable must have led to substantial increase in the pseudo adjusted R² and decrease in Akaike information criteria (AIC; [34]) for it to be included in the model. Here, the nonlinear regression model pseudo adjusted R² [35] is used, a measure of the improvement from the null (intercept only) model to the fitted model. The AIC measures how much (relative) information is lost from using one regression model or another regression model to represent the ‘true model’, and works well for model selection when the true model is infinite-dimensional (the finite-parameterized regression model is then far from being accurate) [36]. Both the AIC and the pseudo adjusted R² include a penalty for added parameters to discourage model over-fitting. After model selection, models were validated with a re-calculation of the regression parameters using 75% of the data, with the other 25% of the data held out for model validation. Model validation was based on metrics of normalized root mean square (NRMSE; root mean square divided by observed data range as a percentage) for the greenness loss model and misclassification percentage for the landslide occurrence model.

All independent variables (hurricane forces and landscape characteristics) were made from data summarized at the resolution of the dependent variable, since they were finer resolution. If a variable was used in both GLMs, it was summarized once at the 250 m resolution of the EVI data, and once at the 2 km resolution of the landslide data. The dependent variable of greenness loss was treated as positive from 0 to infinity, as it may be larger than the initial EVI (due to greenness gain from May to October during quiet years in some locations). It is right skewed, so it was modeled as having a Gaussian distribution with a log link. The dependent variable of landslide occurrence has a Bernoulli distribution with values 0 or 1. Thus, the probability of landslide occurrence equaling 1 is modeled as having a binomial distribution with a logit link.

2.6. Stronger, Slower, Wetter Storm Events

These GLMs were used to model the observed hurricane effects and estimate the importance of the hurricane forces and the landscape characteristics on the magnitude of the effects, as well as to see the result of changing these forces and characteristics. While it is not advisable to use a regression models to predict responses from extremes of independent variables, they can be used to explore smaller changes that are representative of past scenarios and could be representative of future scenarios. Four experiments were run for conceivable future scenarios. First, the wind profile of Hurricane Maria (maximum speed 62–75 m/s near Puerto Rico with an eye radius of ~24 km) was changed into the larger Hurricane Irma (maximum speed 83 m/s near Puerto Rico with an inner eye radius of ~24 km and an outer radius at ~43 km). The average forward speed of Maria was 4.5 m/s, an average forward speed for Atlantic hurricanes at this latitude, whereas Irma traveled quickly, at 7.1 m/s. Second, the forward speed of Maria was slowed to 3.1 m/s, to match the duration with which Hurricane Mitch would have passed (a 1998 category 5 Atlantic hurricane that did not hit Puerto Rico but was notable for the effects of its slow and stalling speed [37]). Third, the rainfall of Maria was multiplied by a factor of 1.62 to match the maximum rainfall totals of Hurricane Harvey (a 2017 category 4 Atlantic hurricane that did not hit Puerto Rico but was notable for the effects of its high rainfall [38]). Fourth, the fraction forest edge was increased. Historical patterns of forest change in Puerto Rico show that while reforestation in the last 20 years was approximately equal to deforestation, deforested sites tend to be located in the forest interiors and reforestation sites on the edges of forests, with deforestation more strongly impacting forest fragmentation than reforestation [21]. So, in the fourth experiment, holes were added to areas with >90% forest to remove 10% of the forest and add 80% edges, i.e., adding maximum edge raster cells for every forest raster cell removed. Although this should alter canopy greenness, the magnitude of EVI change was not clear from the available data and so was not considered.

Effects were calculated spatially over the entire region, but also specifically examined in the Luquillo Experimental Forest (LEF, centered at 18.3°N and −65.8°E), a protected research area of 113 km² that can be seen as the ~100% forested area in northeast Puerto Rico in Figure 2f.

3. Results

The gale strength wind kinetic energy in one meter of air over the U.S. Caribbean from the passage of Hurricane Maria and Irma was calculated as 232 petajoules (PJ, 1 PJ = 10¹⁵ J), 87 PJ of which came from Irma alone. The total rainfall over the land region from Maria alone was calculated as 3.13 km³, a mean value of 34 cm at every location.

The best model for greenness loss (mathematically symbolized as G_L) with the logarithm transformed out for readability, was:

E (G_L) = 0.6 (3.0^W) (1.5^F) (1.3^D) G₀^3.5

(1)

where E is the expected value, W is the kinetic wind energy 100 MJ/m³, F is the fraction forest, D is the fraction forest edge, and G₀ is the initial EVI. All independent variables (right side of Equation (1)) are shown in Figure 2 in the same units, and all had high significance levels (p-values < 10⁻¹⁵). The result of the model validation was good; the NRMSE of the model fit on the training data set (75% of the data) was 10.74%, and on the test set (25% of the data held out in fitting) it increased very slightly to 10.78%. The modeled values (left side of Equation (1)) are shown in Figure 1c. This model has a pseudo R² of 0.73. The model indicates that an average of 0.17 EVI was lost over the U.S. Caribbean from Hurricane Maria and Irma. This is 31% of the initial EVI of the region and 51% of the LEF; 6% of the region and 41% of the LEF lost more than half of its initial (pre-hurricane) EVI.

Let the probability of landslide occurrence (L) being one be mathematically symbolized as P(L). Then, P(L)/(1 − P(L)) is the odds of a landslide occurrence. The best model for the odds with the logarithm transformed out for readability, was:

E (P(L)/(1 − P(L))) = 0.4 (26.0^R) (298.9^C) (80.9^D)S^2.6

(2)

where R is the total rain in meters, C is the fraction clay, D is again the fraction forest edge, and S is the slope in radians. All independent variables are shown in Figure 2 in the same units, and all had high significance levels (p-values < 10⁻⁴). The modeled values of P(L) are shown in Figure 1d. A probability cutoff for landslide occurrence was set to match the observed ratio of landslide occurrence to no occurrence; this model said if P(L) > 0.53 then L = 1. The result of the model validation was good; the misclassification percentage of the model fit on the training dataset was 19.26%, and on the test set it increased slightly to 19.58%, indicating this is a valid model for making predictions of effects of future hurricanes. This model of L is shown in Figure 1e, and has that an average of 34% of the U.S. Caribbean, and 52% of the LEF, had a density of at least 1 landslide per 1 km² as a result of Hurricane Maria.

The results of experiments of conceivable future scenarios are shown in Figure 3 and Figure 4. The figures are split into experiments on force changes (Figure 3) and on landscape changes (Figure 4), which can possibly be managed by with land resource decisions. For experiment one, changing Maria to the size of Irma, with the larger Maria having the duration of Irma (with Irma still passing earlier as observed to the northeast of Puerto Rico), the energy total increased to 308 PJ and 3% more of the initial EVI over the region would have been lost, than in modeled reality (as in Figure 1c) (note: these, and all presented regional average numbers, did not change (to >3 significant digits) if the EVI lost was capped at the initial EVI). This model is not shown, as it would be altering the size and duration of Maria at the same, complicating the result. Instead, Figure 3a shows the wind energy from the larger Irma-sized Maria made to have the observed duration of Maria. Here, 22% more of the region and 34% more of the LEF would have lost over half of its initial EVI than in modeled reality. The larger-sized hurricane did increase the overall amount of disturbance versus the actual hurricane, but it also increased the footprint of the disturbance with substantial greenness loss on the western half of the island (Figure 1c to Figure 3d). For experiment two, changing Maria to have the duration of Hurricane Mitch (Figure 3b), the model overall result of EVI loss (Figure 3e) was comparable to an Irma-sized Maria with Irma-duration (but the spatial distribution was different); 4% more of the region and 11% more of the LEF would have lost over half of its initial EVI, than in modeled reality. The disturbance pattern was largely the same as the actual hurricane, with the greenness loss increasing in the areas that had losses during the actual hurricane. For experiment three, increasing Maria’s rain by a factor to match the maximum rain of Hurricane Harvey (Figure 3c), there was a large increase in the spatial spread of the landslide occurrences from that of the actual hurricane (10% was added from Figure 1e to Figure 3g). For experiment four, increasing fragmentation in heavily forested areas (Figure 4a,b), smaller average changes of 2–3% were seen in the modeled effects of EVI lost and landslide occurrence (Figure 4c–e). Since the areas that were the most forested were the areas that had the most greenness loss from the actual hurricane and also the areas where fragmentation was increased, the greenness loss pattern appears as though the wind force of the hurricane was increased (similar to experiment 2). The landslide occurrences pattern changed the most from the actual hurricane pattern in the areas nearer the coasts; areas that were fragmented and previously below the threshold to slide.

4. Discussion

This assessment allows qualitative insight on the immutable vulnerabilities of future forests, over vulnerabilities that can possibly be managed by with land resource decisions. Overall, there was a very good ability to model the spatial variation of the hurricane effects with the simple models. Furthermore, the percentage of greenness loss in the region (31%) and the much larger amount in the LEF (51%), agrees with the small number of field studies that have been conducted so far [39,40]. The studies are not quite correspondent, since these other studies were focused on litterfall and tree damage (not equivalent to canopy greenness loss) and calculated greenness loss here was up to a month after the hurricane (because of satellite image quality limitations discussed in the methods). Additionally, EVI as measured from MODIS pixels, is a subpixel mix of grass and tree green vegetation. As grasses were unaffected by defoliation, field, or other of tree defoliation measures may be expected to be greater than the loss of EVI modeled here. The good model results in this study indicate that at this scale, the storm was intense enough to not require high-resolution climate simulations of the forces, or the physical actions of the forces, to explain the majority of the variation of the effects. The topographic position and elevation of the land were shown to alter the distribution of hurricane wind energy to substantially influence the magnitude of hurricane effects on forest canopies. Hurricane effects were greater in unsheltered areas of the landscape and at higher elevations, as the wind model has less sheltering at higher elevation because of air-mass characteristics [27]. Additionally, the patterns of forest cover in a landscape also were shown to influence hurricane effects. Taken together, these results indicate it is possible to identify spatial patterns of more or less vulnerability to hurricane effects based on topography and land use decisions. Improved resolution of hurricane effects, i.e., higher resolution imagery or extrapolated field measurements would likely capture finer scale relationships to topography and forest structure. Additionally, more computationally-expensive climate simulations and detailed field studies can provide additional quantitative information about the physical relation of forest-vulnerabilities to the disturbance from Hurricane Maria. Specific cause and effect studies based on Hurricanes Irma and Maria will be useful in modeling effects of other (perhaps weaker) storms, other forests, and future climate projections [8,41,42,43,44]. The conclusions from this study can be used to inform the research hypotheses of these more intensive studies.

The most important variable in this study for modeling hurricane effects was the propensity for disturbance, or the initial greenness amount (EVI) for greenness (EVI) loss, and the land slope for slope-failure (landslide) occurrence. Fundamentally, this disturbance-propensity variable was a measure of hazard, with Equations (1) and (2) showing the traditional power-law function of a stronger than linear growth of effects in relation to compounding influences of physical processes in natural systems [45,46]. The landslide occurrence model explained less of the variability than the EVI loss model, possibly because the landslide data were not as accurate in spatial resolution or magnitude (it was binary) as the EVI data. Furthermore, the relationship between compounding storms (e.g., Irma rain) was not considered in the landslide occurrence model but was considered in the EVI loss model, because the relationship between successive storms and landslides [42] was beyond the scope of the linear models used here.

Landform type (wetland, valley, top, side, or plain; Figure 2h) might be expected to encompass a wide variety of land characteristics that would influence hurricane force interaction with the forest interfaces of canopy and soil. Previous research found land cover (with overlap to the landform divisions here) made a difference in the amount of vegetation disturbance caused Hurricane Maria, with possibly some influence from fragmentation [10]. This study agrees, but found that once wind force with sheltering considerations was included in the model, the significance of landform type was negligible, and a few other variables could explain the differences in disturbance seen among the landforms. The topographic differences in landforms are reflected in the wind values, highlighting the utility of the simple wind model. Sides had the most greenness loss; this was attributed by the GLM to be caused by higher amounts of forest cover and initial greenness here, as sides, valleys, and tops all experienced a similar amount of wind on average. Wetlands experienced twice as much wind on average versus any other landform but did not have much forest to cause greenness loss. These results agree somewhat with field studies from previous hurricanes in the Caribbean [47,48,49]; areas with more trees had more disturbance, but with satellite data it was found that these areas were on sides more than valleys and tops. For landslide occurrence, areas of landslides were far more likely to be seen on sides due to slope, as would be expected.

Both models showed landscape-wide significant positive correlation (which does not mean causation) of fragmentation and hurricane effects, a correlation that was much higher than the correlation of landform and hurricane effects. A two-times increase of fragmentation in an area had 87% of the effect on the model of EVI loss there and a 27% increase on the model of odds of landslide occurrence (in the 1 km²) as did an equivalent increase in forest cover or clay content, respectively. Of note, in the four experiments for possible future scenarios, increasing fragmentation even with a somewhat extreme scenario did not change the average effect nearly to the extent of increasing hurricane force. The reasons for this are trifold. Firstly, when building the GLMs, the fragmentation index that was used did not take into account a threshold on the size of the fragment, which surely must be exceeded for fragmentation to matter. This likely diminished the importance of fragmentation in the GLM and is a model design limitation. Secondly, the reasonable range of force magnitude is much larger than the reasonable range of fragmentation magnitude. This is logical; it is expected to hold in the future. Thirdly, the models’ behaviors in the areas of experimentally increased fragmentation were dominated by the influence of other variables. In the EVI loss model Equation (1), wind force has a much larger impact than fragmentation, and fragmentation is inversely related to the other predictor of forest cover. In the landslide occurrence model Equation (2), the areas of experimentally-increased fragmentation (chosen to be in heavily forested areas) did not correlate well with locations of high clay content, the second best predictor of landslide occurrence after slope. Where clay content was high and the area was heavily forested, the probability of landslide occurrence increased dramatically (see isolated areas at the edge of the mountain-coast transition, seen in Figure 2c,d as areas of high slope next to zero slope with high clay content).

The third set of issues illustrates the complexity of the vulnerability of forests to fragmentation: where it happens matters, in relation to other characteristics and hazards of that area. In light of future climate scenarios calling for more intense or frequent hurricanes, the implications of this study are that land resources should reduce forest edges in areas of high hurricane force vulnerability. For rain causing landslides, these are areas of high soil clay content. For winds causing greenness loss, these are the north northeast areas of Puerto Rico with high wind susceptibility as hurricanes arrive here from the east, southeast, and winds are strongest to the right of the forward velocity, and are strongest immediately as they hit land.

5. Conclusions

A landscape-wide study of hurricane effect was made possible through remotely-sensed data. The spatial variation of this data was used, along with simple hurricane force models connected with landscape characteristics through generalized linear models, to determine forest susceptibility to hurricane effects. First, the effect of the specific hurricane was quantified: it was found for the U.S. Caribbean, Hurricane Maria and caused 31% of the greenness to be lost (including the earlier September disturbance of Irma) and 34% of the region to have a landslide density of at least one in 1 km², with higher percentages for both types of effects in the LEF. This loss was up to a month after the hurricane; immediate loss was likely higher as the forest started recovering right after disturbance [10]. Next, the forest susceptibility to hurricane effects described in the generalized linear models was assessed for landscape-wide characterizations: it was found that while intensity of the storm makes by far the largest scenario difference, forest fragmentation makes a sizable difference, especially in vulnerable areas of high clay content or high wind predisposition. Although more disturbance is seen at higher elevations and on topographic sides, differences in disturbance for different landforms could be explained by the amount of hurricane force sheltering on these areas, and the characteristics the landforms for greenness, forest cover, fragmentation, slope, and soil clay content. Furthermore, the reproducibility of the remotely-sensed spatial pattern of hurricane effects with the fairly simple models suggests that such models could be useful in quickly predicting the effects of approaching hurricanes, especially with large hurricanes that are likely to push landscapes past their abilities to absorb the forces without disturbance.