The complex orography of the Canary Islands favors the creation of microclimates, which cannot be studied using global climate models or regional models with moderate resolution. In this work, WRF is used to perform a dynamic climate regionalization in the archipelago, using the pseudo–global warming technique to compute the initial and boundary conditions from a reanalysis dataset and from results of 14 global climate models. The simulations were performed for three decades, one at present (1995–2004) and two in the future (2045–54 and 2090–99), and for two different greenhouse gas scenarios (RCP4.5 and RCP8.5), defined in phase 5 of the Coupled Model Intercomparison Project. The obtained results, at a 5-km horizontal resolution, show a clear dependence of temperature increase with height and a positive change in diurnal temperature range, which is mainly due to a reduction in soil moisture and a slight decrease in cloud cover. This negative change in soil moisture is mainly a consequence of a decrease in precipitation, although the evaluation of simulated reduction in precipitation does not show statistical significance in most of the Canary Islands for the analyzed periods and scenarios.
In recent decades, impacts on natural and human systems caused by changes in climate, independent of their causes, have been observed on all continents and across the oceans (IPCC 2014). To analyze the changes that could affect these systems in the future, global climate models (GCMs) have proven to be an invaluable tool, mainly for continental and hemispherical climate studies. However, their applicability to regional climate impact studies is limited because their typical spatial resolutions, on the order of hundreds of kilometers, are too coarse to provide useful climate information for applications at regional-scale regimes (Leung et al. 2003). This fact is especially relevant in climate studies on islands with a complex orography, where regional models should have a resolution of a few kilometers (Zhang et al. 2012). Pérez et al. (2014) showed that for the Canary Islands the model resolution should be of at least 5 km to reproduce the observed geographical distribution of temperature and, especially, precipitation.
The climatological differences between the islands that comprise the Canarian archipelago, even the differences between some areas of the same island, have been also studied from observational data (Mestre and Felipe 2012). Martín et al. (2012) studied the mean temperature trends during recent decades for Tenerife, the largest (2034 km2) and highest (3717 m MSL) island of the archipelago, finding that warming is larger for high mountain areas than for lower sites and that windward slopes show a more pronounced increase of temperature than leeward slopes. This larger trend for high altitude sites in Tenerife has been also discussed by other authors (Sanroma et al. 2010). García-Herrera et al. (2003) found an important decreasing trend in precipitation for Tenerife and northern Gran Canaria during the whole second half of the twentieth century, mainly due to a drop in the frequency of the most intense events. However, in the rest of the archipelago, lower nonsignificant trends, or even positive for Lanzarote, were found.
In this work the pseudo–global warming (PGW) method (Kimura and Kitoh 2007; Sato et al. 2007; Kawase et al. 2009) has been used to perform climate regionalization for the Canary Islands, taking the Weather Research and Forecasting (WRF) Model as the regional climate model (RCM). In the PGW approximation, initial and boundary conditions for a recent period, used to compute present climatology, are directly taken from reanalysis data. This is one of the main advantages of this method, because by using reanalysis data those errors in simulating observed climate caused by biases in the boundary conditions from global models are largely diminished (Kimura and Kitoh 2007). For future periods, initial and boundary conditions for the regional model integrations are given by the sum of a climate perturbation signal to the same reanalysis data used for the present simulation. This perturbation, or global warming increment, is estimated from GCM projections for some variables of interest based on their monthly mean values. The PGW method circumvents the uncertainty caused by the interannual variability, even when short periods, around 10 years, are used (Kawase et al. 2008). Another advantage is that this method allows the boundary conditions of the future climate to be based on an ensemble mean of GCM results, that is, the perturbation signal can be estimated not only from the results of one GCM but as the mean change computed from simulations of several GCMs. Kawase et al. (2009) showed that the results for rainfall in East Asia obtained from a PGW simulation that used an ensemble of the data provided by seven models from phase 3 of the Coupled Model Intercomparison Project (CMIP3) as the perturbation signal has good similarity to the ensemble of the results obtained by seven different PWG simulations, each of which takes the climate warming increment from a particular CMIP3 model. Lauer et al. (2013) performed a similar experiment, obtaining parallel results for the Hawaiian Islands, using 10 models from phase 5 of CMIP, and evaluating not only the rainfall but also changes in 2-m temperatures, 10-m wind speed, water vapor path, and trade wind inversion. The possibility of using shorter periods and the use of multi-GCM mean warming increments for high-resolution downscaled climate projections alleviate the computational costs, because in these cases the computational effort to downscale many separate global model projections might be prohibitive (Lauer et al. 2013). Despite the advantages of PGW, as indicated above, this method also has limitations that must be taken into account. It cannot adequately capture potential changes in the variability from daily to interannual time scales, because this method assumes unchanged variability in the future climate (Kawase et al. 2009). Furthermore, it assumes that frequency and intensity of weather perturbations that enter the regional simulation domain also remain unchanged, because they depend on the reanalysis data, not considering possible changes in storm tracks and baroclinicity as a result of warming (Rasmussen et al. 2011). However, previous studies based in CMIP5 results (Lehmann et al. 2014), focused on the study of the expected evolution of these synoptic processes under the different forcing scenarios, suggest that those possible changes will not significantly affect the Canary Islands region.
The outline of this article is as follows. In the next section, the configuration of WRF to perform the experiments and the details of the applied PGW method will be explained. The results for temperature and precipitation changes will be presented in section 3. Finally, in section 4, the main results are summarized.
2. Climate regionalization using pseudo global warming
WRF, version 3.4.1 (Skamarock et al. 2008), was used to perform the simulations of dynamical downscaling, using three domains in a double-nested configuration, with spatial resolutions of 45, 15, and 5 km (Fig. 1). The innermost domain, whose results have been used in this work to evaluate changes in temperature and precipitation, covers the seven islands of the archipelago, is located within 26.75°–30.29°N, 19.43°–12.62°W, and has 133 × 79 grid cells. The general model configuration is the same as in Pérez et al. (2014), where a study of the skills of different WRF simulations to replicate the observed maximum and minimum temperatures, together with the daily rainfall, was conducted. From all the parameterizations offered in WRF, those that provided the best results in that work were selected. Thus, radiation schemes were set to the Community Atmosphere Model, version 3 (CAM3), for computing both longwave and shortwave radiation fluxes (Collins et al. 2004). In the two outer domains, Kain–Fritsch cumulus parameterization (Kain and Fritsch 1990) was used, and no cumulus parameterization was applied in the innermost region because the fluxes can be explicitly resolved at resolutions under 10 km. The planetary boundary layer was characterized using the Yonsei University scheme (Hong et al. 2006) and the land surface using the Noah model (Chen and Dudhia 2001). Finally, the WRF double-moment 6-class (WDM6; Lim and Hong 2010) was used as a cloud microphysics scheme. To limit the divergence of the internal solution computed by WRF from the driven data, which supply the initial and boundary conditions, a grid analysis nudging technique was applied to the outer domain (D1). The higher-resolution inner domains are not nudged, in order to let the regional model create its own structures. For the same reason, in domain D1, nudging is only applied on vertical levels above the boundary layer.
WRF simulations have been conducted for three decades: the present (1995–2004) and the middle (2045–54) and end (2090–99) of the century. The selected present decade includes years with both positive and negative North Atlantic Oscillation (NAO) index (Osborn 2006), which is one of the main climatic patterns affecting the studied region. Four future simulations were performed, two for each decade, using GCM projections based on two different greenhouse gas concentration pathways, the CMIP5 representative concentration pathway 4.5 and 8.5 (RCP4.5 and RCP8.5) scenarios (Taylor et al. 2012). These scenarios represent middle and high emission assumptions, using emission pathways that lead to radiative forcings of 4.5 and 8.5 W m−2, which correspond to greenhouse gas concentrations of approximately 650 and 1370 ppm CO2 equivalent, at the end of this century, respectively (Van Vuuren et al. 2011). For each experiment the model was initialized on 1 January of the year before the first of the decade of interest and integrated continuously until 31 December of the last year. The first year of the 11-yr simulation period is considered as the spinup for the model physics and is excluded from any further analysis. The initial and boundary conditions for the present simulation were obtained from data of the ERA-Interim reanalysis (Dee et al. 2011), specifically data gridded to 0.75° × 0.75° spatial resolution, with 32 vertical pressure levels, and with a frequency of 6 h, were used.
To obtain the climate perturbation signal, present monthly 10-yr climatology was subtracted from the future monthly 10-yr climatologies for each period and emission scenario. The monthly mean climatologies were computed from the ensemble of 14 CMIP5 models (the detailed list of models is summarized in Table 1). All publicly available CMIP5 models, which provided all the needed variables at the time of this study, were selected. Only the first ensemble member of each model, r1i1p1, is used to compute the multimodel ensemble mean, considering that differences in the perturbation signals among the individual models are much larger than those between individual ensemble members from the same model (Lauer et al. 2013). The climate perturbation fields were computed for mean temperature, vapor mixing ratio, geopotential height, and wind at each vertical level including 2-m height, and for sea level pressure, pressure at the surface, and sea surface temperature. The monthly perturbation fields were then added to the present reanalysis fields, for the selected years, by interpolating linearly in time from the monthly climatologies to the 6-hourly time steps, assuming that the monthly mean is valid on the fifteenth of the month, creating the lateral boundary conditions for future WRF simulations.
The main results for simulated changes in temperature and precipitation in the Canary Islands at the middle and at the end of this century for both scenarios (RCP4.5 and RCP8.5) are presented in this section.
Annual and seasonal changes in the daily minimum (Tmin) and maximum (Tmax) temperatures for the two selected future decades (2045–54 and 2090–99) with respect to the present decade (1995–2004) have been analyzed. These two variables were chosen, instead of mean temperature, because they are affected by different climate conditions. To establish if the differences between future and present decades are statistically significant, in this case at a 0.05 significance level, a nonparametric technique has been used to avoid the restrictions of traditional parametric methods that assume that data distributions follow a particular model, as the normal one, and that these data are composed of independent samples from their parent populations. Particularly, a moving-block bootstrap algorithm has been implemented following Wilks (1997), which takes into account the effects of data autocorrelation, very common in atmospheric data, on the hypothesis tests. In this work, an autoregressive-moving average (ARMA) process of order 1 in both contributions, that is, ARMA(1, 1), has been assumed as a time series model for the studied variables. The simulated daily time series of temperature, for most of the simulation grid points, can be actually well represented by a simpler first-order autoregressive process, AR(1), but the ARMA(1, 1) model tends to AR(1) when the moving average contribution is negligible. So, the use of this model has been preferred to also take into account those points whose series are better adjusted by an ARMA process. For each of the grid points, the corresponding ARMA(1, 1) model is computed using daily time series and, based on its characteristics, the block length for the bootstrap test and the adjusted data variance for the test statistic are calculated (Wilks 1997).
Figure 2 shows the present mean values and the mean annual differences in daily minimum temperature for the four future simulations with respect to the present decade. The black dots shade those areas where the computed changes in minimum temperature are not statistically significant. As expected, the simulated changes are larger for the last decade of the century and for the simulations where the scenario with a larger concentration of greenhouse gases (RCP8.5) is used. These increments range from 0.5° to 2°C, showing a clear dependency of temperature change with height. This height dependence of temperature is also observed for Tmax (Fig. 3). That is, in general, an altitude dependence of temperature changes is observed in all the performed simulations. This behavior has already been observed in the Canary Islands, as mentioned in section 1, and has also been observed and modeled in other parts of the world (Giorgi et al. 1997; Pepin and Lundquist 2008; Rangwala and Miller 2012; Rangwala et al. 2013). In this case, the general elevation-dependent change in temperature is not only due to the orography or to the properties of the land surface but to a change in the atmospheric lapse rate, as observed, and simulated by CMIP5 models in tropical and subtropical regions (IPCC 2013; Lott et al. 2013). This change in the temperature profile is present in the climate perturbation signal obtained from CMIP5 models and therefore in the boundary conditions used for simulations of future decades. Temperature vertical profiles computed in an ocean region, north of the islands and inside the simulation innermost domain (Fig. 4), which is not affected by the orography, also show an increase in future simulated temperature with height.
From Figs. 2 and 3, the larger change in Tmax can also be noticed. That is, WRF simulations provide a positive change in diurnal temperature range (DTR) for the four future simulations with respect to the present decade. To facilitate the interpretation of this increase, the mean DTR changes for the different future experiments are shown in Fig. 5. DTR variations can be mainly caused by changes in clouds, soil moisture, precipitation, water vapor, snow/ice albedo, and/or aerosols (Dai et al. 1999; Rangwala and Miller 2012). All these possible modulators, except the aerosol concentration that remains unchanged in future simulations, have been evaluated, finding that, in this case, the DTR positive change is mostly produced by a decrease in the soil moisture and also, but to a lesser extent, by a lower cloud cover in future simulations. Figures 6 and 7 show the simulated changes in soil moisture of the uppermost layer (10 cm) and total cloud cover (TCC). TCC has been computed using the cloud fraction for each output level and assuming maximum random overlap, which corresponds to the same assumption used by the CAM3 radiation schemes (Collins et al. 2004). A decrease in the total cloud cover leads to an increase of the downward shortwave radiation at the surface, which contributes to an elevation of the diurnal maximum temperature. In addition, lower soil moisture reduces the evaporation, which is one of the responsible of limiting the Tmax. A clear correspondence between the areas more affected by changes in soil moisture and TCC, and those where the increase in DTR is larger, can be observed. To study this dependency, the monthly mean changes for each variable have been computed, and then the spatiotemporal correlations have been calculated using only land areas. The results are summarized in Table 2. The simple correlation coefficients are represented by r and the partial correlation coefficients by pr. Thus, for example, for the experiment of the decade 2045–54 and the RCP4.5 scenario, changes in TCC can explain the 20% variance of DTR future increase, assuming a linear relationship. The addition of the other variable, the change in soil moisture, implies that 32% of the variation that was left unexplained by the above-mentioned simple regression can be also explained, showing a larger contribution of soil moisture changes to the DTR increase. Moreover, the soil moisture in decadal periods, as used in this study, is closely related to precipitation and will be presented in the next subsection.
Both the dependence of the temperature change with height and the increment of the DTR were also found in the seasonal study (not shown). Moreover, from this seasonality study, a larger change in Tmin and Tmax can be observed during summer.
The mean change in annual precipitation, for the four future experiments, is shown in Fig. 8. Unlike what happens in the case of temperature, the changes are not statistically significant, except in very small areas and only at the end of the century. For the estimation of the statistical significance the same method used for temperature has been applied, assuming also an ARMA(1, 1) model. In any case, a decrease in annual precipitation is predicted in the whole archipelago, being more noticeable in the mountainous areas, which correspond to zones with larger annual rainfall in the present decade.
The seasonal behavior (not shown) is very similar to the annual results, where the winter season is the most important for the total precipitation and, therefore, for future precipitation decrease. To analyze the causes of this rainfall decrease, the mixing ratio vertical profiles of the condensates and water vapor provided by the different simulations corresponding to the winter season, when the rainfall events are predominant in the Canary Islands, have been computed. Figure 9 shows these vertical profiles averaged over an area that covers Tenerife (Fig. 4a), which is the island with the most complex orography and the largest climatic contrasts in the archipelago. Despite the increase in water vapor mixing ratio, as expected in a warmer climate, a global decrease in all condensate species can be observed. This reduction is more noticeable for higher altitudes, that is, for graupel, ice, and snow, because of the dependence of temperature increase with height. The reduction of cloud water content leads, in turn, to a decrease in rainwater.
Future climate high-resolution projections have been performed for the Canary Islands using WRF to implement a downscaling dynamical method based on the pseudo–global warming technique. Simulated changes in temperature, the daily maximum and minimum values, and precipitation have been analyzed. These changes, with respect to the present decade, have been studied for two periods, at the middle and at the end of the century, and for two different greenhouse gas scenarios. In all of these experiments a general rise of temperature was detected, finding a clear dependence of temperature with elevation. Thus, the largest temperature changes are expected in the high mountains. Furthermore, an increase in the DTR has been also detected, which is produced, in the performed simulations, by a decrease of soil moisture and, to a lesser extent, by a reduction in the total cloud cover.
A general decrease of annual precipitation is also expected, even though the obtained results are statistically significant only in small areas of the islands and at the end of the century. However, despite this lack of significance, this reduction of precipitation is the main cause of the decrease in soil moisture, which, in turn, causes the increase of DTR.
As indicated in section 1, despite the advantages of the PGW method, it has some drawbacks, from which the unchanged frequency and tracks of storms constitute the largest source of uncertainty, although in previous studies, based on CMIP5 data, no changes are predicted in the Canary Islands region. Thus, this study can be considered as an imposed warming experiment (Rasmussen et al. 2011) and not as a robust regionalization of future climate scenarios, which requires multiple simulations, using a different GCM to compute the boundary conditions for each of them. However, the obtained results, at a 5-km horizontal resolution, are the most comprehensive obtained so far for the archipelago, constituting an important dataset for the study of possible climate effects on ecological and social systems. In future works, it would also be desirable to increase the spatial resolution of the simulations to improve the representation of the complex orography of some of the islands and their effect on atmospheric processes.
The authors acknowledge the MEC (Ministry of Education and Science, Spain) for support of projects CGL2010-21366-C04-01 and CGL2010-22158-C02-01 (CORWES). The authors also thank the CLIMATIQUE Project cofinanced by the Programa de Cooperacin Transfronteriza España-Fronteras Exteriores 2008-2013 (POCTEFEX) and the European Union funds. We also acknowledge Program FEDER for economic support. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling institutions (listed in Table 1) for producing and making available their model output. For CMIP the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. The WRF simulations performed in this study were managed by WRF4G, which is an open-source tool funded by the Spanish government and cofunded by the European Regional Development Fund under Grant CGL2011-28864.