Recent changes and drivers of the atmospheric evaporative demand in the Canary Islands

. We analysed recent evolution and meteorological drivers of the atmospheric evaporative demand (AED) in the Canary Islands for the period 1961–2013. We employed long and high-quality time series of meteorological variables to analyse current AED changes in this region and found that AED has increased during the investigated period. Overall, the annual ET o , which was estimated by means of the FAO-56 Penman–Monteith equation, increased signiﬁcantly by 18.2 mm decade  1 on average, with a stronger trend in summer (6.7 mm decade  1 ). In this study we analysed the con-tribution of (i) the aerodynamic (related to the water vapour that a parcel of air can store) and (ii) radiative (related to the available energy to evaporate a quantity of water) components to the decadal A relative homogeneity method was applied to identify possible inhomogeneities.


Introduction
The atmospheric evaporative demand (AED) is one of the key variables of the hydrological cycle , with multiple implications for agriculture, hydrology and the environment (Allen et al., 2015). Several studies have indicated that current global warming is increasing the intensity of the hydrological cycle, mainly as a consequence of an intensification of the AED (Huntington, 2006). Sherwood and Fu (2014) suggested that mechanisms driving the AED over land regions could be the main driver of increasing climate aridity in world semi-arid regions under a global warming scenario.
Warming may play an important role in increasing the AED via the aerodynamic component (McVicar et al., 2012a). Following the Clausius-Clapeyron relationship, the quantity of water vapour that a given mass of air can store increases exponentially with the air temperature. Nevertheless, there are other climate variables whose temporal evolution could compensate for the increased AED induced by increasing air temperature, such as wind speed and vapour pressure deficit (McVicar et al., 2012a). In addition, the radiative component of the AED, which is related to the available solar energy that transforms a unit of liquid water into vapour, may compensate for or accentuate the increase in AED associated with warming. Wild et al. (2015) noted that solar radiation had increased over large regions since the 1980s as a consequence of changes in cloud cover and/or atmospheric aerosol concentrations.
These large number of variables interact in a non-linear manner to determine the AED (McMahon et al., 2013), so assessing recent changes in the AED and defining their determinant factors is not an easy task. For this reason, while several studies have analysed the AED at the global scale using different datasets and methods, there is no general con-sensus on the recent AED evolution (Sheffield et al., 2012;Matsoukas et al., 2011;Dai, 2013). In this context, the few existing direct AED observations, based on evaporation pans, show a decrease since the 1950s at the global scale (Peterson et al. 1995;Farquhar, 2002, 2004), a finding that adds more uncertainty regarding the behaviour of the AED under current global warming. These issues stress the need for new studies that employ high-quality datasets to assess the time evolution of the AED at the regional scale.
There are a number of studies published in the last decade that analysed the AED evolution across different regions of the world. Some of them are based on AED estimated using empirical formulations, mostly based on air temperature data (e.g. Thornthwaite, 1948;Hargreaves and Samani, 1995). However, to adequately quantify the AED evolution it is necessary to use long time series of the meteorological variables that control its radiative and aerodynamic components (e.g. air temperature, vapour pressure deficit and wind speed). Although these variables are generally poorly measured and highly inhomogeneous over both space and time, numerous regional studies analysed the evolution of the AED by means of the robust Penman-Monteith (PM) equation using long time series of these variables. The available regional studies show quite contradictory results, where some studies showed AED negative trends, including those in China (Xu et al., 2006;Ma et al., 2012;Zhang et al., 2007;Liu et al., 2015) and northwestern India (Jhajharia et al., 2014). In contrast, other regional studies have found positive trends in AED, including those in central India (Darshana et al., 2012), Iran (Kousari and Ahani, 2012;Tabari et al., 2012), Florida (Abtew et al., 2011), continental Spain (Espadafor et al., 2011;Vicente-Serrano et al., 2014a;Azorin-Molina et al., 2015), France (Chaouche et al., 2010) and Moldova (Piticar et al., 2015).
The contrasting trends among world regions would be a consequence of the evolution of the different meteorological variables that control the AED. Specifically, some studies suggest that temporal variability and changes in the AED are related to changes in the relative humidity, mainly in semi-arid regions Espadafor et al., 2011;Vicente-Serrano et al., 2014b), whereas others stress the importance of solar radiation Ambas and Baltas, 2012;Fan and Thomas, 2013) or wind speed (McVicar et al., 2012b).
Among these studies, few analysed the AED variability and trends and their possible drivers in the eastern North Atlantic region (Chaouche et al., 2010;Vicente-Serrano et al., 2014a;Azorin-Molina et al., 2015). Nevertheless, there are no studies about this issue in the subtropical areas of the North Atlantic region. In this area, there are very few meteorological stations measuring long-term series of the variables necessary to make robust calculations of the AED. This uneven distribution of meteorological observatories constrains the high interest to know the evolution of atmospheric pro-cesses in this region, where climate variability is strongly controlled by changes in the Hadley circulation (Hansen et al., 2005) that affect the position and intensity of the subtropical anticyclone belt. Knowing the evolution of AED and its main drivers in this region is highly relevant given the general climate aridity of the region and the low availability of water resources (Custodio and Cabrera, 2002). In this work we analyse the recent evolution and meteorological drivers of the AED in the Canary Islands. The main hypothesis of the study is that, in opposition to other continental temperate regions of the North Hemisphere, the warm and humid climate of the subtropical Canary Islands provides the water supply to the atmosphere needed to maintain the AED constant under the current global warming scenarios; consequently, only wind speed and solar radiation could affect the observed decadal variability and trends of the AED. Thus, the availability of long and high-quality time series of meteorological variables in the Canary Islands provides an opportunity to analyse current AED changes in the subtropical northeastern Atlantic region and the role played by different meteorological variables.

Dataset
We used the complete meteorological records of the Spanish National Meteorological Agency (AEMET) in the Canary Islands for the following variables at the monthly scale: maximum and minimum air temperature (308 stations), wind speed (99), sunshine duration (42) and mean relative humidity (139). The majority of the stations cover short periods or are affected by large data gaps. As the number of meteorological stations before 1961 was very little for several variables, we restricted our analysis to the period between 1961 and 2013. Specifically, only eight meteorological stations had data gaps of less than 20 % of the months in all the necessary variables. As illustrated in Fig. 1, these stations are distributed between the islands of Tenerife (three stations), Gran Canaria (two), La Palma (one), Lanzarote (one) and Fuerteventura (one). Given that some series included records for a longer period (e.g. Izaña from 1933 and Santa Cruz de Tenerife from 1943), neighbouring stations with shorter temporal coverage were used to reconstruct the existing data gaps in the selected observatories, using a regression-based approach. Details of the site names, coordinates, relocations, data gaps and inhomogeneities of the selected meteorological stations can be found in Table 1.
Then, the time series were subject to quality control and homogenization procedures. The quality control procedure was based on comparison of the rank of each data record with the average rank of the data recorded at adjacent stations (Vicente-Serrano et al., 2010). A relative homogeneity method was applied to identify possible inhomogeneities.  For this purpose, we used HOMER (HOMogenization soft-warE in R), which compares each candidate series with a number of available series (Mestre et al., 2013). The method provides an estimation of break points in the time series relative to other stations, indicating high probabilities of the presence of inhomogeneities. This method was applied to the different variables and time series following Mestre et al. (2013). Finally, a single regional series for the different variables was obtained using a simple arithmetic average of data values at the available eight stations.

Calculation of ET o
The Penman-Monteith (PM) equation is the standard technique for calculation of ET o from climatic data (Allen et al., 1998), and it is the method officially adopted (with small variations) by the International Commission for Irrigation (ICID), the Food and Agriculture Organization (FAO) of the United Nations, and the American Society of Civil Engineers (ASCE). The PM method can be used globally, and has been widely verified based on lysimeter data from diverse climatic regions (Itenfisu et al., 2000;López-Urrea et al., 2006). Allen et al. (1998) where ET o is the reference evapotranspiration (mm day 1 ), R n is the net radiation at the crop surface (MJ m 2 day 1 ), G is the soil heat flux density (MJ m 2 day 1 ), T is the mean air temperature at 2 m height ( C), u 2 is the wind speed at 2 m height (m s 1 ), e s is the saturation vapour pressure (kPa), e a is the actual vapour pressure (kPa), e s e a is the saturation vapour pressure deficit (kPa), 1 is the slope of the vapour pressure curve (kPa C 1 ), and is the psychrometric constant (kPa C 1 ). The FAO-56 PM is an equation initially designed for crop monitoring and irrigation operation at daily and subdaily scales. This equation involves non-linear relationships among the variables used for calculation and averaging these variables for long-term intervals could affect the reliability of the ET o estimations. Nevertheless, Allen et al. (1998) indicated that the FAO-56 PM equation can be used for daily, weekly, 10-day or monthly calculations, and several previous studies have computed the PM ET o using monthly values for some variables (e.g. Sheffield et al., 2012;Dai, 2013). We have found that using monthly averages instead of daily records for the different variables has no relevant influence on the ET o estimations in the Canary Islands. Figure 2 shows an example using two of the available stations (Los Rodeos and Izaña) for the 1978-2010 period. The relationship between the monthly sum of the daily ET o calculations and the ET o calculation from the monthly averages justifies the equality of applying both procedures. This is observed for the ET o monthly values (including seasonality but also considering monthly standardized anomalies in which seasonality is removed). Moreover, there are other technical reasons that lead to recommending the use of monthly instead daily records to calculate ET o since testing and correcting the temporal homogeneity of the necessary variables on a daily basis is highly problematic, whereas testing and correcting homo-geneity using monthly records is reliable (e.g. Venema et al., 2012). Therefore, the monthly ET o was calculated from data of the monthly averages of five meteorological parameters: maximum and minimum air temperature, relative humidity (which allows calculating the vapour pressure deficit), wind speed at a height of 2 m, and daily sunshine duration (which allows estimating the net radiation). García et al. (2014) compared the capability of sunshine duration series to reconstruct long-term radiation in the observatory of Izaña (Tenerife), showing very good temporal agreement between sunshine duration and radiation, independently of the season of the year. Further details on the required equations to obtain the necessary parameters from meteorological data can be consulted in Allen et al. (1998).
We also calculated the evolution of the radiative (Eq. 2) and the aerodynamic components (Eq. 3) of the ET o , as follows:

Analysis
Using the time series of ET o , we determined the seasonal (winter: December-February; spring: March-May; summer: June-August; autumn: September-November) and annual ET o averages. To analyse changes in ET o we used the nonparametric Mann-Kendall statistic that measures the degree to which a trend is consistently increasing or decreasing. The Mann-Kendall statistic is advantageous compared to parametric tests as it is robust to outliers and it does not assume any underlying probability distribution of the data (Zhang et al., 2001). For these reasons, it has been widely used for trend detection in a wide range of hydrological and climatological studies (e.g. Zhang et al., 2001;El Kenawy and McCabe, 2015). Autocorrelation was considered in the trend analysis applied to the series of ET o , the series of the aerodynamic and radiative components of the ET o , and the series of the different climate variables (temperature, relative humidity, wind speed and sunshine duration). This was applied using the FUME R package, which performs the modified Mann-Kendall trend test, returning the corrected p values after accounting for temporal pseudoreplication (Hamed and Rao, 1998;Yue and Wang, 2004). To assess the magnitude of change in ET o , we used a linear regression analysis between the series of time (independent variable) and the ET o series (dependent variable). The slope of the regression indicated the amount of change (ET o change per year), with higher slope values indicating greater change. We also calculated the trend observed in the different meteorological variables (air temperature, relative humidity, sunshine duration and wind speed) at both the seasonal and annual scales.
To get insight into the influence of changes in the different meteorological variables on ET o , we related the evolution of ET o to relative humidity, maximum and minimum air temperature, wind speed, and sunshine duration by means of correlation analyses. To assess the importance of trends in the different meteorological variables on the observed trends in ET o between 1961 and 2013, we applied the PM equation while keeping one variable stationary (using the average from 1961 to 2013) each time. This approach provided five simulated series of ET o , one per input variable, which could be compared to the ET o series computed with all the data to determine the isolated influence of the five variables. Significant differences between each pair of ET o series (the original one and the alternative one in which one variable was kept constant) were assessed by comparing the slopes of the linear models, with time as the independent variable. A statistical test for the equality of regression coefficients was used (Paternoster et al., 1998). The significance of the difference was assessed at a confidence interval of 95 % (p < 0.05). Figure 3 shows a box plot with the seasonal and annual values of ET o at the different meteorological stations across the Canary Islands, which are also summarized in Table 2. There were strong seasonal differences in ET o , as all different meteorological stations show their maximum values in summer and minimum in winter, albeit with strong differences among them. In winter, the highest average values were recorded in the most arid islands (i.e. Fuerteventura and Lanzarote) and at the station of Los Rodeos (northern Tenerife). In summer, the stations of Izaña and Los Rodeos showed the highest av- erage values (663.8 and 612.9 mm, respectively). The lowest summer ET o averages were recorded at the stations of Gran Canaria (San Cristóbal and Gran Canaria/Airport). At the annual scale, there were very few differences in the average values between the stations of Los Rodeos, Izaña, Fuerteventura and Lanzarote, with very high ET o values ranging between 1693 and 1784 mm ( Table 2). The observatory with the lowest ET o values is located at Gran Canaria/Airport, although the observatory of San Cristóbal (also in Gran Canaria) records the minimum values in summer. The magnitude of the differences can be quite important (up to 34 %) between the highest ET o values recorded in Los Rodeos, Izaña, Fuerteventura and Lanzarote and the lowest ET o values (Gran Canaria and San Cristóbal). In general, variability, as revealed by the coefficient of variation, was higher at the meteorological stations that recorded the highest ET o values at the annual scale, but there was no clear spatial pattern at the seasonal scale as different stations showed few differences in terms of the coefficients of variation (Table 2).

Average ET o values
At the majority of weather stations the seasonal and annual ET o magnitude was mostly driven by the aerodynamic component. The average aerodynamic fraction was higher than the radiative fraction at the weather stations that record the highest ET o values (Los Rodeos and Izaña) in all seasons around the year (Fig. 4). At other weather stations (Santa Cruz de Tenerife and San Cristóbal), the ET o associated with the radiative component was much higher than that observed for the aerodynamic component (Table 3). The temporal variability in the aerodynamic component was much higher than that observed in the radiative one, regardless of the season of the year or the meteorological station.

Long-term evolution of ET o
The regional ET o series for the whole of the Canary Islands (Fig. 5) shows a significant increase at the annual scale (18.2 mm decade 1 ), which is stronger in summer (6.7 mm decade 1 ) (Table 4). Nevertheless, there was a strong variability between the different meteorological stations, since most meteorological stations experimented significant increases of ET o between 1961 and 2013. The largest annual increase was recorded in Los Rodeos (34.8 mm decade 1 ), La Palma (29.8 mm decade 1 ) and Lanzarote (29.7 mm decade 1 ). Considering a longer period (1933( -2013( for Izaña, and 1943( -2013 for Santa Cruz de Tenerife), the changes are not statistically significant, although it was not possible to check the homogeneity of the climate records prior to 1961 and thus the results for the longer period must be carefully considered. For the period 1961-2013, there is no general spatial pattern in the observed changes; thus, some differences can be observed. For example, in Gran Canaria, San Cristóbal station shows a statistically non-significant negative change in ET o on the order of 8.4 mm decade 1 , while there is a general significant increase of 28.4 mm decade 1 at Gran Canaria Airport.
Trends in the aerodynamic and radiative components showed clear differences among stations and for the average of the Canary Islands (Fig. 6). Main changes were recorded in the aerodynamic component. The regional series showed an increase of 16.2 mm decade 1 in the aerodynamic component, but it only showed an increase of 2 mm decade 1 in the radiative component (Table 5). This can be translated to an average increase in the ET o of 89 % over the whole period due to changes in the aerodynamic component, and of 11 % due to changes in the radiative component. However, there are spatial differences between the meteorological stations, since the aerodynamic component showed a decrease of 21 mm decade 1 in San Cristóbal, compared to an increase of 44.6 mm decade 1 in Los Rodeos. However, the radiative component showed lower differences among stations, with values ranging from 9.9 mm decade 1 in Los Rodeos to 12.7 mm decade 1 in San Cristóbal. Nevertheless, and regardless of the observed trends, the results indicate that the interannual variability in ET o between 1961 and 2013 was   Table 6). The temporal correlation between ET o and the aerodynamic component was statistically significant for the different meteorological stations in the seasonal and the annual series, with correlation coefficients higher than 0.95 in most cases. The correlation for the regional series was also strong and statistically significant. In contrast, the correlation coefficients calculated between ET o and the radiative component were much lower, and generally non-significant (p < 0.05). Los Rodeos is the unique weather station where the correlation between ET o and the radiative component was statistically significant at both the seasonal and annual scales, but showing a negative correlation. Overall, the results show that the correlation between the annual radiative component and the total annual regional series of ET o is statistically nonsignificant. Table 7 shows the correlation between the different meteorological variables and ET o at the seasonal and annual scales at the eight meteorological stations. Maximum and minimum air temperatures were positively correlated with ET o , and this relationship was statistically significant at some stations, and the correlation coefficients tended to be higher for maximum air temperature. In Los Rodeos and La Palma, the ET o variability could not be explained by the variability in air temperature, with correlation coefficients weaker than 0.3. Overall, the results indicate that the seasonal and annual series of ET o were significantly correlated with variations in sunshine duration and wind speed, suggesting that these two variables are the key drivers of ET o variability in the Canary Islands. The variable that showed the strongest correlation with the evolution of ET o in the seasonal and annual series of the different meteorological observatories was relative humidity, with negative coefficients. Only in the annual series of Santa Cruz de Tenerife was the correlation non-significant. Moreover, there were no significant differences in the magnitude of correlations among seasons. The regional series summarize the pattern observed at the individual meteorological stations (Fig. 7). In winter, relative humidity had the strongest correlation with ET o (r = 0.85), with a mostly linear relationship. Minimum air temperature and sunshine duration showed significant positive correlations with ET o (r = 0.40 and 0.36, respectively). Maximum air temperature and wind speed showed weaker cor-relation with the winter ET o . In spring, the magnitude of the correlations was similar among the different variables, and the highest correlation corresponded again to relative humidity (r = 0.72). A similar pattern was found in summer, when relative humidity showed the strongest correlation (r = 0.74) followed by maximum and minimum air temperature. In autumn, relative humidity also showed the strongest correlation and wind speed showed more importance than both maximum and minimum air temperature. As expected, relative humidity showed the strongest correlation with ET o (r = 0.83) at the annual scale, followed by wind speed (r = 0.62). However, the correlation with maximum air temperature was statistically non-significant. The general increase observed in ET o in the Canary Islands was largely determined by changes in the different meteorological variables (Table 8). The maximum air temperature does not show noticeable changes, with the exception of Gran Canaria/Airport, Lanzarote and San Cristóbal stations, where significant trends were found. The regional average did not show significant changes. However, the minimum air temperature showed an average increase of 0.12 C decade 1 in summer and 0.09 C decade 1 at the annual scale between 1961 and 2013. The significant increase recorded in summer was found at six meteorological stations, with a maximum of 0.25 C decade 1 in Izaña. Changes in relative humidity were also significant. There was a significant decrease in winter, summer and annually, which represent a decline of 0.47 % decade 1 , although there were differences among stations. Sunshine duration and wind speed did not show noticeable changes, and the unique remarkable pattern was the significant increase in the summer sunshine duration at the regional scale (0.12 h decade 1 ) and the significant increase in wind speed at the station of Los Rodeos in the four seasons and also annually.

Drivers of ET o variability and trends
With respect to the sensitivity of changes in ET o to its five driving meteorological drivers (Fig. 8), substantial differences were found between variables. The differences between observed ET o and simulated ET o with average maximum and minimum air temperature were small irrespective of the season, indicating a low sensitivity to these two variables. In contrast, ET o was more sensitive to setting sunshine duration and wind speed at their mean values. Thus, at the station of Los Rodeos, the predicted magnitude of change in winter, autumn and annually was different from the observed magnitude of change. The highest sensitivity was, however, to relative humidity. In general, the different meteorological stations showed an important increase in observed ET o with respect to predicted ET o keeping relative humidity constant. This was observed at the seasonal and annual scales. Thus, at three meteorological stations the observed magnitude of change on annual basis is between 2 and 3 times higher than that predicted considering relative humidity as stationary. This pattern was also found in the regional series ( Fig. 9). Considering air temperature, sunshine duration and wind speed as constant, there were no statistical differences between the observed and predicted magnitudes of change, both seasonally and annually. However, with relative humidity left constant, the magnitude of the trend was quite different to the observations, and temporal trends would not be statistically significant. Thus, the magnitude of change of ET o , considering relative humidity as constant, is significantly different from the observed magnitude of change in winter and annually.

Discussion
This work analyses the recent evolution   logical stations in which the necessary meteorological variables for calculation of the ET o were available. The results showed a general increase in ET o , although different magnitudes of change were found between the different meteorological stations. These differences did not follow any specific geographic pattern, so they must be considered to be due to either random effects and uncertainty at various levels or micro-geographic effects that were not considered in this study. There is no general pattern that may connect the observed trends in a certain forcing variable with the observed trend of ET o at each of the eight analysed stations, although those that showed a higher increase in ET o (i.e. Lanzarote, Los Rodeos and Gran Canaria) displayed a higher increase in the aerodynamic component, a process which is in agreement with the significant reductions observed in relative humidity. Nevertheless, with the exception of the observatory of San Cristóbal in the north of Gran Canaria, other meteorological observatories showed positive changes in ET o , with annual trends statistically significant at six stations. In any case, we must also stress that trends in ET o at the regional scale are mostly significant because of the low values in the beginning of the study period starting in the 1960s. Thus, the results of the two sites with longer temporal coverage (i.e. Izaña and Santa Cruz de Tenerife) do not show significant trends. This makes necessary to consider these trends with caution since they could be driven by variability processes at the decadal scale.
The few existing studies in northwestern Africa (Ouysse et al., 2010;Tekken and Kropp, 2012) are not comparable with our findings, since the variables required to apply the PM equation were not available. Instead, these studies relied on simplified methods that just employ air temperature records. Despite the difference in methods, these studies also found a general increase in the ET o . The closest region in which it is possible to make a direct comparison using the same method is the Iberian Peninsula, where a general increase of 24.5 mm decade 1 was found between 1961 and 2011 (Vicente-Serrano et al., 2014a). This study also found that the variability and trends in the aerodynamic component determined most of the observed variability and the magnitude of change of ET o at the majority of the meteorological stations in the Iberian Peninsula. The radiative component showed much lower temporal variability than the aerodynamic component did. Thus, more than 90 % of the observed ET o variability at the seasonal and annual scales can be associated with the variability in the aerodynamic component. This is in agreement with the results obtained in previous studies. For example,  showed that recent ET o variability at the global scale was mainly driven by the aerodynamic component. Equally, other studies in southern Europe indicated a higher importance of the aerodynamic component Azorin-Molina et al., 2015). It could be argued, however, that quantification of the radiative component in our study was based on a simplified assumption since it was calculated from sunshine duration, which is mostly determined by cloud coverage (Hoyt, 1978). Nevertheless, it is also worth noting that global radiation measurements and sunshine duration records contain a signal of the direct effects of aerosols (Sanroma et al., 2010;Sanchez-Romero et al., 2014;Wild, 2015) in the Canary Islands. Nevertheless, the Canary Islands is a region mostly free of anthropogenic aerosols given the large frequency and intensity of trade winds (Mazorra et al., 2007), and it is not expected that the frequency of Saharan dust events, which Figure 7. Relationship between the regional annual and seasonal ET o and the regional series of the different meteorological variables.
Pearson's coefficients are included in each plot. In bold are the coefficients that statistically significant at the 95 % confidence level. could affect incoming solar radiation, has noticeably changed over the last decades (Flentje et al., 2015;Laken et al., 2015). Consequently, in the Canary Islands we can consider high accuracy in determining the radiative component using sunshine duration series. In continental Spain, Azorin-Molina et al. (2015) also found strong positive correlations between interannual variations in solar radiation and sunshine duration at different meteorological stations. Overall, in the Canary Islands there is a positive and significant correlation between interannual variations in ET o and sunshine duration, although this correlation did not explain the observed trends of ET o in the region.
We showed that the temporal variability in ET o is strongly controlled by the temporal variability in relative humidity. Specifically, seasonal and annual series of ET o at the differ-ent stations showed very strong negative and significant correlations with those of the relative humidity. Thus, the magnitude of correlations was much higher than those obtained for other meteorological variables, and this finding was common to the whole set of meteorological stations. This strong control of relative humidity on the temporal variability in ET o has been already identified in some studies in the Iberian Peninsula (Vicente-Serrano et al., 2014b;Azorin-Molina et al., 2015;Espadafor et al., 2013).
Among the variables that control the aerodynamic component, wind speed and maximum air temperature did not show significant trends at the regional scale, and only a few stations recorded significant trends in these variables at either the seasonal or the annual scale. Significant trends were obtained for minimum air temperature, mainly in summer. Recently, Cropper and Hanna (2014) analysed long-term climate trends in the Macaronesia region, and for the Canary Islands they showed an increase in air temperature during summer for the period 1981-2010. Martín et al. (2012) analysed air temperature changes in Tenerife from 1944 to 2010, and they also showed that night-time air temperature increased rapidly compared to daytime temperature. Nevertheless, they found strong spatial contrasts between the high mountains, which showed a higher increase, and the coastal areas in which the air temperature regulation of the ocean could be reducing the general air temperature increase.
In any case, the variable that recorded more significant changes in the Canary Islands was relative humidity, and among the different meteorological variables used to calculate ET o , relative humidity was the main driver of the observed ET o trends. Significant negative humidity trends were recorded in winter, summer and autumn, but also annually.
Thus, simulation of ET o series considering the different meteorological variables as constant produced few differences in relation to the observed evolution of ET o , with the exception of the relative humidity. Leaving relative humidity as constant for the period 1961-2013 showed no significant ET o changes at seasonal and annual scales and also statistically significant differences with changes obtained from observations. In continental Spain, Vicente-Serrano et al. (2014b) showed a general decrease in relative humidity from the decade of 1960, mainly associated with a general decrease in the moisture transport to the Iberian Peninsula as well as a certain precipitation decrease. Similarly, Espadafor et al. (2011) andVicente-Serrano et al. (2014b) showed that the strong increase in ET o in the last decades is associated with the relative humidity decrease due to air temperature rise, which caused more severe drought events (Coll et al., 2016; Gallardo et al., 2016). In the Canary Islands, no precipitation changes have been identified during the analysed period (Sánchez-Benitez et al., 2016). Therefore a lower moisture supply from the humidity sources to the islands should explain the observed pattern toward a relative humidity decrease. Sherwood and Fu (2014) suggested that differences in the air temperature increase between oceanic and continental areas could increase land aridity, as a consequence of the subsaturation conditions of the oceanic air masses that come to the land areas, given higher warming rates in maritime regions in comparison to continental areas. The results of this study confirm this pattern in the Canary Islands, since this region should not be constrained by constant moisture supply from the surrounding warm Atlantic Ocean. Overall, Willett et al. (2014) recently found a general decrease in relative humidity at the global scale, including several islands and coastal regions in which the moisture supply was expected to be unlimited. This finding suggests that contrasted mean air temperature and trends between land and ocean areas could also play an important role in explaining this phenomenon, even at local scales.

Conclusions
We found that the reference evapotranspiration ET o increased by 18.2 mm decade 1 -on average -between 1961 and 2013 over the Canary Islands, with the highest increase recorded during summer.
Although there were noticeable spatial differences, this increase was mainly driven by changes in the aerodynamic component, caused by a statistically significant reduction of the relative humidity.
This study provides an outstanding example of how climate change and interactions between different meteorological variables drive an increase in the ET o event in a subtropical North Atlantic archipelago. Given the general aridity conditions in most of the Canary Islands and the scarcity of water resources, the observed trend could have negative consequences in a number of water-depending sectors if it continues in the future.
Data used here are licensed by the Spanish National Meteorological Agency (AEMet) and cannot be redistributed. Potential users of the raw data can directly contact AEMet (http://www.aemet.es/).