The potential ecological and economic effects of climate change for tropical islands were studied using output from 12 statistically downscaled general circulation models (GCMs) taking Puerto Rico as a test case. Two model selection/model averaging strategies were used: the average of all available GCMs and the average of the models that are able to reproduce the observed large-scale dynamics that control precipitation over the Caribbean. Five island-wide and multidecadal averages of daily precipitation and temperature were estimated by way of a climatology-informed interpolation of the site-specific downscaled climate model output. Annual cooling degree-days (CDD) were calculated as a proxy index for air-conditioning energy demand, and two measures of annual no-rainfall days were used as drought indices. Holdridge life zone classification was used to map the possible ecological effects of climate change. Precipitation is predicted to decline in both model ensembles, but the decrease was more severe in the “regionally consistent” models. The precipitation declines cause gradual and linear increases in drought intensity and extremes. The warming from the 1960–90 period to the 2071–99 period was 4.6°–9°C depending on the global emission scenarios and location. This warming may cause increases in CDD, and consequently increasing energy demands. Life zones may shift from wetter to drier zones with the possibility of losing most, if not all, of the subtropical rain forests and extinction risks to rain forest specialists or obligates.
Anthropogenic increases in atmospheric CO2 are projected to drive significant changes in temperature and precipitation patterns, means, and extremes (IPCC 2014), which will in turn affect a wide range of social and ecological systems (Mideksa and Kallbekken 2010; Grimm et al. 2013). Tropical islands are among the most vulnerable areas in terms of climate change effects on biodiversity, vegetation, and economy (e.g., Fordham and Brook 2010; Jennings et al. 2014; Singh and Bainsla 2015). Assessing and communicating these potential effects in a way that allows land managers, planners, governments, and the public to make informed decisions requires a cascade of modeling efforts that can include several tiers: 1) physically based general circulation models (GCMs) tied to greenhouse gas (GHG) emission scenarios that reflect potential future human and societal choices, 2) downscaled climate model projections based on GCMs and GHG scenarios that take finer-scale information and physical processes into account, 3) modeling the conditional effects of projected climatic change onto a wide range of social and ecological drivers, and 4) model output that specifically addresses decision-maker questions—for example, timing, thresholds, cumulative effects, human and other species’ responses to potential future conditions, feedbacks from potential response scenarios to climate, and confounding effects of complex decision scenarios. Making global climate projections useful to decision-makers is a multifaceted endeavor that requires expertise in many fields, from physics to ecology to social science, and communication with decision-makers to develop knowledge relevant to decision needs (Gould et al. 2016). The scale and content of the management information needs are typically finer and more variable as information cascades to lower tiers and targeting users at all stages improves science delivery. While global-scale products are by nature globally applicable, there is considerable regional variation in the availability and applicability of appropriately scaled information for management decisions. Tropical islands often have both finescale spatial variability with steep gradients in climate and heterogeneous land use, and lack detailed downscaled climate projections. These islands harbor a large number of endemic species that are acutely vulnerable to global climate change (Fordham and Brook 2010), highlighting the importance of interpreting the geospatial aspects of the global-scale climate projections at the local scales within the islands. An example of such a heterogeneous tropical island is the Caribbean island of Puerto Rico (Fig. 1). Modeling finescale implications of global projections has been problematic in Puerto Rico as the resolution of global-scale predictions is too coarse for evaluating the effects on local dynamics. Harmsen et al. (2009) used the statistical downscaling approach of Miller et al. (2008) to downscale the twenty-first-century predictions from a global-scale GCM to three locations in Puerto Rico, while Hayhoe (2013) used an asynchronous regional quantile regression-based model (Stoner et al. 2013) for the statistical downscaling of global-scale GCMs to the Caribbean region and Puerto Rico. The development of these downscaled datasets has allowed for the first time a more decision-relevant analysis of the potential effects of anthropogenic climate change for Puerto Rico’s socioecological ecosystems.
In this study, we used the predictions from 12 GCMs under three global GHG emission scenarios of the United Nation’s Intergovernmental Panel on Climate Change (IPCC) (Nakicenovic and Swart 2000), downscaled to station locations in Puerto Rico by Hayhoe (2013). The objectives were 1) to create continuous gridded surfaces of the predicted climatic variables by interpolating the downscaled GCM predictions at station locations to the entire island; 2) to assess how climate scenarios and model selection strategies influence subsequent effect projections; 3) to estimate the potential societal consequences of changes in precipitation and temperature in Puerto Rico through the end of the current century by using climatic indices; and 4) to estimate the ecological consequences of precipitation and temperature changes by mapping the potential bioclimatic zones for certain vegetation types called life zones (Holdridge 1947) under the climate change scenarios. The Holdridge life zones system has been used to map life zones at a given time based on measurements of climatic variables [e.g., Ewel and Whitmore (1973) in Puerto Rico and the U.S. Virgin Islands; Sawyer and Lindsey (1963) in the eastern and central United States; and Lugo et al. (1999) in the conterminous United States]. It was also used as a criterion to map the ecological effects of climate change at a global scale (e.g., Emanuel et al. 1985; Sisneros et al. 2011) or at the scale of some countries (e.g., China; Yue et al. 2001; Chen et al. 2003). But the modeled climate change effects on life zone distribution in tropical islands have not been examined so far. The Holdridge life zone of a particular area implies the potential vegetation type that can dominate in terms of the basic climatic conditions if local conditions (e.g., edaphic, topographic, and socioeconomic conditions) also coexist with the climatic condition. Therefore, changes in the Holdridge life zones due to climate change imply eventual shifts in the potential climatic support for different vegetation types and species. The remainder of the paper is structured as follows. In section 2, we describe the statistically downscaled climate data used and model selection strategies, and then we describe the methods used for interpolating climatic variables and mapping Holdridge life zones. Section 3 presents the projected maps of precipitation, temperature, and life zones. In section 4, we discuss some further implications of the projected changes in the climatic variables, indices, and life zones, and describe potential sources of uncertainties. Section 5 presents our conclusions.
2. Materials and methods
a. Study area
The island of Puerto Rico is located at ~(17°45′–18°30′)N, ~(65°45′–67°15′)W in the Caribbean region (Fig. 1). Its size is about 8740 km2, and it has a predominately maritime climate (Daly et al. 2003). Island temperatures have a linear negative relationship with elevation (Goyal 2014). Puerto Rico is topographically diverse in terms of elevation and slope. Elevation ranges from sea level to 1338 m above the sea level in the central mountains. Therefore, climatic conditions across the island are highly variable. The latest ecological life zone map of Puerto Rico showed six life zones ranging from subtropical dry forests to subtropical rain forests (Ewel and Whitmore 1973) (Fig. 1).
b. Downscaled climate data
Downscaled climate model output for individual locations in Puerto Rico were taken from Hayhoe (2013) based on the Asynchronous Regional Regression Model (ARRM; Stoner et al. 2013). This is a statistically downscaled dataset of GCMs from phase 3 of the Coupled Model Intercomparison Project (CMIP3) using three IPCC global GHG emission scenarios: mid-high (A2), mid-low (A1B), and low (B1) (Nakicenovic and Swart 2000). The ARRM method uses a piecewise regression method to predict daily temperature and precipitation values based on the relationship between observations and coarse-scale GCM output. An attractive feature of this method is improved accuracy in simulating extremes, an important consideration given the occurrence of very large 24-h precipitation events (>100 mm) across the island. The GCM outputs for precipitation and temperature were downscaled to 72 and 29 station locations (Fig. 1) respectively by Hayhoe (2013).
For each GCM (Table 1) the downscaled daily outputs were available for 140 yr (1960–2099) and three variables: precipitation and maximum and minimum temperature. Because life zones represent successional climaxes that evolve over decades, we calculated multidecadal means for each variable for five time intervals: 1960–90, 1991–2010, 2011–40, 2041–70, and 2071–99. The downscaled projections were expected to match the reality better if used as 20–30-yr averages (Hayhoe 2013). The duration of the first period was set to match the climate layer means used for spatial variability (Daly et al. 2003). For the rest of the time span the 30-yr averages were set for later periods to allow better comparison between the first and last periods of the century. The difference among the climate scenarios does not apply to the first time interval (1960–90). The mean annual temperature was derived from maximum and minimum temperatures and the averages were calculated for the same intervals. Then, for each variable after calculating the multidecadal mean for a downscaled GCM, we calculated the multimodel mean of all selected climate models under each scenario and time period (Table 1).
c. Climate model selection
Differing GCM parameterizations lead to substantial structural uncertainty for many climate variables at local and regional scales (cf. Masson and Knutti 2011; Cook et al. 2010). The result is that in some instances, a subset of the full GCM ensemble may be more skillful at simulating important regional climate features, which may warrant the inclusion of model selection or model weighting strategies (Knutti 2010). Annual rainfall season in the Caribbean region has a bimodal pattern with the early precipitation mode from May to June, a “midsummer drought” (MSD) in July, and the late precipitation mode from August to November (Angeles et al. 2010). Ryu and Hayhoe (2014) examined CMIP3 and CMIP5 GCM performance in the context of the ability to reproduce the bimodal pattern in the region and found distinct differences between GCMs. Models with higher skill had features that included more accurate representation of the westward extension of the North Atlantic subtropical high into the Caribbean basin in early summer and higher sea surface temperatures that more closely match observed conditions. Accordingly, GCMs were segregated into three categories: “bimodal” models that simulated both precipitation modes and the MSD, “single with MSD” models that simulated the MSD but only the late precipitation mode, and “single” models without the MSD and only the late precipitation mode.
Given these results, we used two different types of GCM ensembles: skill-based ensemble and an ensemble of all available downscaled GCMs. The first model selection strategy uses the subset of CMIP3 bimodal GCMs available in the downscaled dataset from Hayhoe (2013). Figure 2 shows annual precipitation at two locations with different climatic conditions: Rio Blanco Lower station (18°14′35.88″N, 65°47′5.9994″W) near El Yunque National Forest in the Luquillo mountains and Ensenada station (17°58′22.0794″N, 66°56′44.88″W) in southern dry region. The downscaled precipitation by the three categories of the models for A2 scenario is compared with the available observations from 1960 to 2014. Bimodal models show a higher range of variability through time than do the other two categories [standard deviations (SD) are 684.6, 494.1, and 392.1 mm for average bimodal, single with MSD, and single models at Rio Blanco Lower station; SD: 287.5, 269.5, and 200.7 mm for the same model sets at Ensenada station]. Table 2 compares the projected precipitation changes from the period 1960–90 to 2070–99 for the same sets of models, the climate scenario, and stations used in Fig. 2. The multimodel precipitation projection for the bimodal GCMs showed larger declines than for the other two categories (Table 2). The second model selection strategy we employed was to use all 12 available downscaled CMIP3 GCMs for the three scenarios regardless of the model categorization to maximize the sample space of projected climate change.
d. Climatologically aided interpolation
Interpolation is required to transform the downscaled climatic predictions at station points to a continuous surface. Kriging is a widely used interpolation method that results in unbiased optimal estimates of climatic variables based on the spatial correlation between values, assuming a Gaussian distribution of the variable (Cressie 1993). However, kriging and other traditional interpolation methods may result in too much smoothing of surfaces (Brown and Comrie 2002) and in cases of inadequate station densities with distributions unrepresentative of topographic variability and geographic features the outputs may become even smoother and less representative (Daly et al. 2002). An alternative to the kriging-only method is climatically aided interpolation (CAI) (Willmott and Robeson 1995), which results in the interpolation errors not highly correlated with the number of stations. The CAI has been used by some studies as a multivariate method to reduce the between station uncertainty by incorporating the spatial variability related to elevation and topography (Willmott and Robeson 1995; Hunter and Meentemeyer 2005; Wilson and Silander 2014). In the CAI approach spatial variability can be incorporated from other spatially rich datasets. We used output from the Parameter–Elevation Regressions on Independent Slopes Model (PRISM; Daly et al. 2003) as the basis of the CAI in order to transform the original downscaled GCM outputs into interpolated surfaces of climate variables. PRISM is based on spatial dependence of climatic variables on elevation, topography, and coastal proximity. Using PRISM, Daly et al. (2003) mapped the monthly averages of maximum and minimum temperature and precipitation of the island from 1963 to 1995 at a spatial resolution of ~450 m. Our use of the PRISM output for CAI assumes that the relationship between climate and the physiographic covariates remains stationary through time. The final outputs are the interpolated downscaled climate variables, with most of the spatial variability being captured by the PRISM layers. First, the downscaled climate variable means for each time interval at each location were transformed into anomalies from the 1963–95 PRISM gridded means [Eqs. (1) and (2)]:
where Pa is the transformed anomaly for precipitation, Pd is the downscaled precipitation mean at each station, Pp is the mean precipitation from the PRISM grid at the station’s location, Ta is the temperature anomaly, Td is the downscaled temperature mean at each station, and Tp is the mean temperature from the PRISM grid at the station’s location.
We used ordinary kriging to interpolate the anomalies across the island to the same spatial resolution as the PRISM grid. The semivariogram was plotted for each variable and time interval and an exponential model was chosen for all of the variograms because their change patterns in autocorrelation by changes in distance lags were exponential. Two fitting methods were compared to fit the exponential model to the variograms: the residual maximum likelihood and ordinary least squares method (Stein 1999). We applied one-out cross-validation of the observations in both methods using the “krige.cv” function in the “gstat” package (Pebesma 2004) in R (R Development Core Team 2013). We selected the fitting method based on the goodness of fits of the cross-validation procedure—that is, the method which resulted in higher Pearson’s R2 (the full goodness-of-fit results are given in the supplemental material). The chosen kriging model was then applied to the means of the anomalies at station points. The interpolated anomalies were then converted back to the original units by inverting the relationships in the Eqs. (1) and (2):
where PCAI is the final precipitation projection for each scenario and time interval, Pai is the transformed anomaly for precipitation interpolated to the surface grid, Ppg is the PRISM surface for mean precipitation, TCAI is the final temperature projection for each scenario and time interval, Tai is the transformed anomaly for temperature interpolated to the surface grid, and Tpg is the PRISM surface for mean temperature.
e. Accuracy assessment
We compared the interpolated downscaled GCM means for the two historic time intervals (1960–90 and 1991–2010) with the observed means taken from NOAA’s Global Historical Climatology Network (GHCN) daily dataset (Menne et al. 2012). We calculated the total annual precipitation and mean annual maximum and minimum temperature from the daily observations at each station and calculated the average for each time period after removing the years with missing records for more than 10% of the days. For each variable we plotted the interpolated and downscaled GCM mean in each time interval against the observed mean in all measured years of the time interval.
f. Climatic indices
We derived three annual indices from the climatic variables (mean daily temperature and precipitation) to transform the changes in the variables to more ecologically and economically oriented measures. We calculated them for each available downscaled station, and also did a trend analysis on the changes in the indices for the average of all available stations versus time.
Cooling degree-days (CDD) were calculated annually using the mean daily temperature [Eq. (5)]. It is defined as the annual sum of differences between mean daily temperature and a base temperature Tb above which indoor cooling is assumed to be needed. Note that Tb is based on the geographic locations. We used the commonly used Tb of 18°C (Rosenthal et al. 1995; Isaac and van Vuuren 2009). CDD is widely used as a proxy for air-conditioning energy demand as it indicates building cooling requirements (e.g., Zhou et al. 2013; Zubler et al. 2014):
where Ti is mean daily temperature (°C), Tb is the base temperature (18°C), n is the total number of days each year, and di = 1 if Td,i > Tb or else di = 0.
Annual maximum number of consecutive dry days (MCDD) measures the maximum length of dry spell each year. It is the longest period of consecutive days with no precipitation or less than 1 mm of precipitation. MCDD was listed as a climate index for estimating drought in the IPCC special report on extreme events (Seneviratne et al. 2012), and it has been used as an index of extreme drought conditions (e.g., Nakaegawa et al. 2014b).
Total annual dry days (TDD) is a measure of drought intensity, unlike MCDD, which is a measure of drought extremes. We calculated MCDD and TDD values using the “run length encoding” (rle) function in R that counts the consecutive repetitions of the values for a variable, in this case days with no or less than 1 mm precipitation. We calculated the percentage of changes in the climate indices from 1960 to 2099. We applied the nonparametric Mann–Kendall (MK) test (Mann 1945; Kendall 1975) to see if there was a monotonic trend of increase in the climate indices. Tau (τ) in the MK test is analogous to the correlation coefficient in regression analysis. We also applied linear regression to assess for linear trends after checking for normality of the calculated indices (Shapiro test, p < 0.05).
g. Life zones mapping
The boundaries of life zones in the Holdridge system are based on three climatic measurements: annual precipitation, biotemperature, and the ratio of potential evapotranspiration (PET) to annual precipitation (Holdridge 1947). The all model selection averages were used for life zone mapping. Biotemperature is an annual temperature index with values below 0°C adjusted to 0°C and values above 30°C adjusted to 30°C. We calculated sea level temperature from the mean annual temperature and elevation layers using a lapse rate of 6.0°C km−1. Then we created a biotemperature layer for each scenario and time interval from the sea level temperature layer by adjusting values above 30°C to 30°C. To estimate the annual PET we applied the Hargreaves–Samani temperature difference method (Hargreaves and Samani 1985) for daily PET and multiplied the result by 365 to get the annual amount [Eq. (6)]:
where PET is annual potential evapotranspiration (mm), RA is extraterrestrial radiation (mm day−1), TD is (Tmax − Tmin) (°C), and T is mean annual temperature (°C).
Extraterrestrial radiation (RA) is the solar radiation at the top of the atmosphere on a horizontal surface. It is a function of latitude, sunset angle, solar declination, and the relative Earth–sun distance (Allen et al. 1998). It was calculated for the 15th day of each month at different latitudinal degrees in the Northern and Southern Hemispheres (Allen et al. 1998). For the purpose of calculating average annual PET in Puerto Rico we used the average of all monthly values of RA at 18° latitude in the Northern Hemisphere. The Hargreaves–Samani method is considered a simple and practical method (Lu et al. 2005; Ramírez et al. 2011). The ratio of PET to mean annual precipitation acts as an aridity index in mapping Holdridge life zones.
We used the three variables for classifying the life zones in eCognition Developer software v8 (www.ecognition.com), which provides hierarchical classification in an object-based classification approach. First, we applied a multiresolution segmentation process (Blaschke et al. 2004) to create image objects as the areas with homogeneous values of PET ratio, biotemperature, and precipitation. Then, we assigned the created objects to different life zone hexagons based on the thresholds in each of the three variables in the Holdridge framework. Since the three variables are not orthogonal, the life zones can be mapped by any two of the three variables (Lugo et al. 1999; Sisneros et al. 2011). We labeled the life zones based on the PET ratio or humidity index and the latitudinal regions defined by biotemperature. We applied a fuzzy adjustment defined in Lugo et al. (1999) to biotemperature thresholds before using in naming the life zones.
a. Accuracy assessment of downscaled climate means
Downscaled and interpolated variables for A2 scenario in the first two time intervals are plotted against their observed values (Fig. 3). The climate scenarios do not apply to the first time interval and the scenario differences were very low at station locations in the second time interval. For precipitation 59.85% and 43.90% of the stations in the first and second periods had modeled values with under 100-mm difference with the observations (Figs. 3a,b; 132 stations in the first and 82 stations in the second time interval). The temperature difference between the model projections and the observations in the first and second periods were under 1°C at 90.38% and 68.57% of the stations for maximum temperature, and at 86.54% and 68.57% of the stations for minimum temperature (Figs. 3c–f; 52 stations in the first and 35 stations in the second period).
b. Precipitation and drought
GCM selection made large differences in the calculations for precipitation. The mean pixel decline in rainfall was 510.67, 354.60, and 312.57 mm for A2, A1B, and B1 scenarios respectively from the first to the last time interval based on the multimodel average of all 12 models (Fig. 4a), whereas the corresponding declines were 916.30, 842.62, and 619.58 mm using the bimodal models (Fig. 4b). Precipitation was projected to decrease faster in the wetter regions of the island such as Luquillo and the central mountains. Figure 5 shows the projected annual precipitation in the five time intervals using the average of all model predictions (the results for the bimodal models are in the supplemental material). Precipitation is projected to change gradually in the entire area and the GHG emission scenarios showed more difference in later periods of the century (Fig. 5). The TDD index as a measure of drought intensity showed a monotone and linear increase for all three scenarios at all station locations (Fig. 6). There was a gradual increase in the mean of MCDD as the measure of drought extremes (Table 3, Fig. 7) with the highest total change for the A2 scenario (71.61%). The mean station trends for all three scenarios were significant and linear (Table 3). However, the trends of change were different for different geographic areas. The interannual variability increased and dominated from the wet to dry stations with the highest fluctuations in the south coast (Fig. 7). The climate scenarios showed similar patterns of changes in precipitation and drought intensity and extremes but total increases from 1960 to 2099 were higher for the A2 scenario (Figs. 6–8, Table 3).
The projected means from the all-model ensemble showed increases of 7.5°–9°, 6.4°–7.6°, and 4.6°–5.4°C under the A2, A1B, and B1 scenarios respectively. The bimodal models gave similar results except for the B1 scenario, which showed a 5.5°–7°C increase (Fig. 4). The scenarios differentiate the most from each other in the last time interval (Figs. 8e,j,o). The last two time intervals in the A2 scenario (Figs. 8d,e) showed higher rates of changes relative to the other period of the century (Figs. 8a–c) and also relative to the last two periods in the A1B and B1 scenarios. This higher temperature change rate in the A2 scenario in the late time intervals is also shown in the CDD index in the late decades of the century (Fig. 9). The total change in CDD from 1960 to 2099 was dependent on the GHG emission scenarios. There was more than 100% total increase in CDD in the A2 scenario (Fig. 9, Table 3). The CDD increased in all 28 stations located across the island. The south coast showed the highest annual CDD from 1960 to 2040 after which the San Juan metropolitan area was the highest location (Fig. 9).
d. Life zones under changing climate
Dramatic changes were projected in the life zone distributions in Puerto Rico in this century (Fig. 10). The general pattern is similar to the changes in the climatic variables since life zones were derived from those variables. The humidity shifting from rain, wet, and moist zones to drier zones was identified by the changes in annual precipitation and PET ratio. The changes from subtropical to tropical were identified by the changes in adjusted biotemperature. The subtropical rain forest, which was mapped in the life zone maps of 1973 (Fig. 1), disappears after the first time interval. Generally, decreasing trends were observed in the areas of wet and moist zones while increasing trends were observed in the areas of dry zones in all three scenarios.
Temperature projections showed good accuracy based on the comparison of projected model means with real observations from station locations. The precipitation accuracy assessment showed that 59.85% and 43.90% of the modeled stations fall under 100-mm precipitation difference with the measured observations for the same stations in the first and the second time intervals. This level of accuracy shows a high potential uncertainty in modeling precipitation in the region. Caribbean rainfall is influenced by complexities in large-scale atmosphere and ocean dynamics and the results are impacted by systematic biases in model structure (Ryu and Hayhoe 2014). The projected amount of precipitation decline was highly dependent on the model selection strategies. The average of bimodal models showed twice as much precipitation loss as the average of all models. The all-model selection showed a gradual effect of GHG emission increase with lower rainfall for higher emission scenarios through time (Fig. 5). In total the warming and drying were more prominent in the late periods of the century, which may be because most climate models predict increased warming of the atmosphere by shortwave radiation after midcentury (Donohoe et al. 2014). The immediate potential consequences of precipitation decline are shown by the measures of drought intensity and extremes. The annual total dry days increased almost monotonically at all station locations in the island (Fig. 6). Drier locations such as the south coast showed more interannual fluctuations in the annual maximum number of consecutive dry days when compared to wet regions such as the central mountains (Fig. 7). MCDD is an extreme drought index and changes in its severity are most notable in drier areas (e.g., Nastos and Zerefos 2009). Polade et al. (2014) examined the trends of TDD at the global scale using CMIP5 model predictions from 1960 to 2089 and suggested a future increase in the index in Central America and the Caribbean. The increases in predicted MCDD and TDD imply rising probability of annual wildfires. Changes in TDD also imply interannual precipitation variability as more dry days means fewer precipitation events (Polade et al. 2014). Warming will tend to accelerate the overall hydrological cycles, which means intensifying the wet extremes in addition to dry extremes (Allan and Soden 2008), that is, increasing the likelihood of floods.
In contrast to precipitation, the projected changes in temperature were not dependent on model selection because GCMs could not be differentiated based on their ability to reproduce the climatology of temperature in the region (Hayhoe 2013). The projected warming in Puerto Rico was greater than the global averages. For example, for A2 scenario there was 8.16°C mean pixel warming for the entire island from the first to the last time interval, whereas the predicted global warming is 3°–4°C for the same climate scenario in this century (IPCC 2007). This projected warming may further impact the changes in precipitation in both space and time (Hayhoe 2013). The projected increase in temperature and CDDs suggests energy demand will increase on this tropical island. The total increase in CDD was high for all three scenarios at all locations. Air-conditioning energy consumption is often assumed to be approximately proportional to CDD in a way that a one percent increase in degree-days corresponds to one percent increase in energy demand (e.g., Zubler et al. 2014). However, detailed calculation of building energy demand would also require predicting wind climatology, which is not possible with the statistically downscaled GCMs. The total climate change effect on building energy demand also depends on the scale of the observation. Zhou et al. (2013) found that in the conterminous United States and China the increased cooling requirements in hot areas are balanced with decreased heating demands in cold areas and the overall effect is modest (6% or less). Rosenthal et al. (1995) suggested that climate change could reduce rather than increase energy costs to the U.S. economy at large. But at the local scale for Puerto Rico with its hot and humid climate, the effect is focused on air-conditioning energy demand, which at present is met by burning fossil fuels to produce electricity and is significantly more expensive than the average heating energy (Energy Information Administration 1994).
The projected changes in the life zones of Puerto Rico based on the all-model selection showed a shift among nine life zones, switching from humid to drier life zones (Fig. 10). This includes changes in the relative area and distribution pattern of the life zones and disappearance of the humid life zones. The Holdridge life zone names of an area do not necessarily imply that a particular vegetation type dominates as they only depict the conditions that regulate ecosystem functions (Lugo et al. 1999). Other local conditions such as soil, proximity to water, and even other climatic conditions affect the actual vegetation types of a given area at a particular time. However, changes in the life zones show the shifts in climatic support for different vegetation types. For example, the areas projected to shift from moist forest to dry forest will not have sufficient climatic support for moist forests regardless of the local conditions. This may suggest future loss of large forested areas particularly rain and moist forests. The predicted forest loss may further contribute to global climate change by 1) decreasing local carbon sequestration capacity (e.g., Backéus et al. 2005) and 2) exacerbating effects of land-use and land-cover change on air temperature (e.g., Comarazamy et al. 2013).
The wide variety of potential local consequences would be particularly severe for the island’s natural environment and society. Puerto Rico has unique ecosystems and species that provide natural and cultural capital for the economic, social, and environmental well-being of the island’s communities. Decreased rainfall may turn the rain forest in northeastern Puerto Rico to wet forests, which will alter the environmental gradients and cause migration, distribution changes, and potential extirpation of many species that depend on the unique environmental constraints of the rain forest (Weaver and Gould 2013). Endemic species abound in this region, as exemplified in the El Yunque National Forest, which harbors 8 tree species, 12 avian species, 1 reptile species, and 11 amphibians (Weaver 2012). The existence of these species will depend on their potential to survive in wider climatic ranges.
Results from this study have implications for impacts on socioeconomic sectors of the region such as agriculture, health, energy, and natural resources. In the case of Puerto Rico and other Caribbean islands, impacts will be exacerbated because small island states are especially vulnerable, given limited resources, markets, and their relative isolation (IPCC 2014). The anticipated drier and hotter conditions can be expected to produce severe yield losses in agriculture production (e.g., Jones and Thornton 2003). Animal production in the future may require enclosed structures with climate control, which is currently not a common practice in Puerto Rico. Increasing tree mortality in forests can be expected due to physiological stress related to climate-related interactions such as insect pressure and forest fires (Allen et al. 2010). Complex changes can be expected to occur in the ecosystems of the island related to altered timing of growth stages in flora and fauna (McMichael and Haines 1997). Increased temperatures may increase illness and mortality among the young, elderly, and poor (McGeehin and Mirabelli 2001). Insufficient water supply is currently a significant problem in the island. Based on experience from droughts (previous and currently during 2015), a large-scale redesign of the island’s water supply system will be required. Coupling the threat of drier conditions from climate change with the need to significantly increase local food production (approximately 85% of food in Puerto Rico is currently imported), providing a reliable future water supply represents one of the greatest challenges to the island.
The first and largest source of uncertainty for end-of-century projections is associated with the choice of long-term GHG emission pathway. The results are also accompanied by systematic biases in the structure of the CMIP3 GCMs (Ryu and Hayhoe 2014). Furthermore, these uncertainties are difficult to quantify in GCM ensembles where traditional forecast verification is typically not possible and model dependence, bias, and tuning are all present (Tebaldi and Knutti 2007; Wilson and Silander 2014). The choice of downscaling method adds yet another source of uncertainty. For example, the maximum average temperature of the PCM for Lajas station using the ARRM method of Stoner et al. (2013) was more than 1°C higher than in Harmsen et al. (2009), which used the method of Miller et al. (2008) for statistical downscaling.
Another limitation is the limited number of locations where downscaled projections were available. For example, we had no model output in the El Yunque National Forest, which experiences the peak rainfall amounts on the island and is an important ecological site. The predicted values of the climatic variables at this location were derived from the values of the closest stations and spatial covariance dependence using PRISM datasets in the CAI approach. Other uncertainties are related to projecting the possible shifts in future life zones. We chose to apply the Hargreaves–Samani method for estimating PET ratios based on the limited number of downscaled climatic variables. The Hargreaves–Samani method has been shown to overestimate PET in humid regions and underestimate PET in arid regions when compared to the Penman–Monteith method (Jensen et al. 1990; Harmsen et al. 2002). Furthermore, life zones are a broad-scale overview of the environmental conditions. The Holdridge life zones are based on annual means of climatic variables (here the annual means of multidecadal means) and ignore extreme events and seasonality, which can have more influence in the distribution of woody plant species than annual means (Ohmann and Spies 1998). We stress that the Holdridge system is not used in this research to map the real life zones. Instead, we applied it as a criterion or scale to map the possible ecological effects of a changing climate. More work is needed to periodically update the life zone maps using the observed climatic variables. However, the question of how often the life zones should be updated remains open and is dependent upon estimating the time scale under which significant life zone changes can occur due to shifts in the underlying climatic conditions.
The original data and the projected results of this study are publicly accessible (data accessed 15 September 2015 from http://caribbeanlcc.org/data-center/). The methods can be repeated to reproduce new maps using other climate models and downscaling methods and to compare the new results with the maps presented in this study. The ability of the CMIP ensembles as a whole has improved over time. It is likely that CMIP3 models overestimate the projected precipitation decreases in wet seasons and in extreme rainfall events because they are driven by changes in carbon dioxide as compared with the CMIP5 model projections, which are driven by changes in both carbon dioxide and aerosols (Hayhoe 2013). However, the downscaled CMIP5 projections are not still available for Puerto Rico. Also, different projections would almost certainly result from the use of dynamic downscaling instead of statistical downscaling, as dynamic downscaling directly models the entire system instead of regressing spatial patterns. Although, Nakaegawa et al. (2014a) indicated that dynamic downscaling does not necessarily provide better information and improvements will depend on the type of applications and variables being analyzed. Future work will examine these differences upon completion of a suite of high-resolution (2 km) dynamically downscaled models. The dynamically downscaled results include more climatic variables, which can improve the potential evapotranspiration calculation and life zone projection.
We examined the potential effects of climate change on ecological life zones, increasing energy demands, and drought indices in Puerto Rico from 1960 to 2099 using a suite of statistically downscaled CMIP3 GCMs. The full model ensemble projected a 130- to 1397-mm decrease in precipitation (from 1960–90 mean to 2071–99 mean) depending on location and GHG emission scenario. The bimodal model selection projected almost twice this decline. Furthermore, there were projected gradual linear increases in drought intensity and extremes with their end of the century increases depending on the emission scenarios. The consequences could include water supply deficits in the future and attendant consideration of adaptation strategies to mitigate these impacts (e.g., water harvest practices). The projected warming in the same period was 4.6°–9°C depending on the geographic location and emission scenarios, and was insensitive to the model selection criteria. We illustrated potential implications using energy demands for air conditioning, and projections suggested that there could be a greater than a 100% increase in energy demand in the absence of intervening changes in technology and consumer behavior. Climate change may alter the life zones of the island with shifts from rain, wet, and moist zones to drier zones. This includes the loss of the subtropical rain, moist, and wet forests and appearance of tropical wet, moist, dry, and very dry forests. New ecological conditions may result in new ecosystems and new communities. For example, present trees that require soil moisture throughout the year may be replaced by other tree and shrub species. Our results point to severe vulnerabilities that exist for the rain forests in the area in the future, which could endanger the endemic species of the island and changes their distributions.
We thank our anonymous reviewers. Any use of trade, product, or firms names is for descriptive purposes only and does not imply endorsement by the U.S. government. All research at the USFS International Institute of Tropical Forestry is done in collaboration with the University of PR.
Publisher’s Note: This article was revised on 8 April 2016 to correct the institutional affiliation of the fifth author, which was misidentified when originally published.