Potential evaporation trends over land between 1983 – 2008 : driven by radiative fluxes or vapour-pressure deficit ?

We model the Penman potential evaporation (PE) over all land areas of the globe for the 25-yr period 1983– 2008, relying on radiation transfer models (RTMs) for the shortwave and longwave fluxes. Penman’s PE is determined by two factors: available energy for evaporation and ground to atmosphere vapour transfer. Input to the PE model and RTMs comprises satellite cloud and aerosol data, as well as data from reanalyses. PE is closely linked to pan evaporation, whose trends have sparked controversy in the community, since the factors responsible for the observed pan evaporation trends are not determined with consensus. Our particular interest is the temporal evolution of PE, and the provided insight to the observed trends of pan evaporation. We examine the decadal trends of PE and various related physical quantities, such as net solar flux, net longwave flux, water vapour saturation deficit and wind speed. Our findings are the following: Global warming has led to a larger water vapour saturation deficit. The periods 1983–1989, 1990– 1999, and 2000–2008 were characterised by decreasing, increasing, and slightly decreasing PE, respectively. In these last 25 yr, global dimming/brightening cycles generally increased the available energy for evaporation. PE trends seem to follow more closely the trends of energy availability than the trends of the atmospheric capability for vapour transfer, at most locations on the globe, with trends in the Northern hemisphere significantly larger than in the Southern. These results support the hypothesis that global potential evaporation trends are attributed primarily to secular changes in the radiation fluxes, and secondarily to vapour transfer considerations. Correspondence to: C. Matsoukas (matsoukas@aegean.gr)


Introduction
There have been many reports on significant changes in pan evaporation for widespread locations across the globe, e.g.European Russia, Siberia, and the USA (Peterson et al., 1995;Golubev et al., 2001), India (Chattopadhyay and Hulme, 1997), the USA (Lawrimore and Peterson, 2000), Israel (Cohen et al., 2002), China (Liu et al., 2004), Australia (Roderick and Farquhar, 2004), the Tibetan plateau (Zhang et al., 2007), and many others (references in Roderick et al. (2007), Roderick et al. (2009), and Fu et al. (2009)).However, a study for Australia (Jovanovic et al., 2008) cast some doubt on the veracity of the Australian trends.It reported that unaccounted discontinuities (e.g. the installation of bird meshes) produced a spurious declining trend, which disappears if the time series are homogenised.There is a need for similar studies in other regions, but until they are carried out, the emerging general picture is one of worldwide decreasing pan evaporation, with some exceptions, e.g.east USA, a single pan in Israel, central Australia, etc.However, Global Circulation Model (GCM) runs on the one hand (Wetherald and Manabe, 2002), and reanalyses (Trenberth et al., 2005) and empirical hydrological evidence (Huntington, 2006) on the other, dictate that in a warming world the hydrological cycle should be enhanced and evaporation should increase.There is also a GCM study by Roeckner et al. (1999), predicting instead a reduction of evaporation in the changing climate of the next decades.The responsible physical process is the dampening of turbulent heat fluxes due to aerosol-caused reconfiguration of energy exchange.
C. Matsoukas et al.: Potential evaporation trends when Brutsaert and Parlange (1998) brought forward the complementary hypothesis, already proposed by Bouchet (1963) and applied by Morton (1975).It consists of the following reasoning: in humid environments, with ample supply of moisture to the surface, the actual evaporation takes values close to the "true" potential evaporation.Also, the pan evaporation corresponds directly to the potential evaporation, after multiplication with the pan coefficient.However, in arid environments actual evaporation cannot reach the values of potential evaporation, and therefore a large portion of the available energy takes the form of sensible heat flux, thus warming the atmosphere.The warmer atmosphere now has a larger water vapour deficit and is characterized by the "apparent" potential evaporation, which is larger than the "true" potential evaporation.In other words, the same location has a "true" potential evaporation if there is ample water on the surface, and a larger "apparent" potential evaporation in water-limited conditions.In arid conditions, the pan evaporation is related to the apparent potential evaporation.In these lines, if the available energy for evaporation is constant, an increase in actual evaporation decreases the sensible heat flux, decreases the apparent potential evaporation, and finally decreases the pan evaporation.Lawrimore and Peterson (2000) after examining precipitation and pan evaporation trends in the USA, and assuming that precipitation and actual evaporation have to be related over large areas, lend support to the complementary hypothesis.Golubev et al. (2001) after deriving actual evaporation rates over the former Soviet Union and the USA, also find the complementary hypothesis reasonable.Zhang et al. (2007) calculated potential and actual evaporation from 16 watersheds in the Tibetan plateau and analysed them for the validity of the complementary hypothesis.They found indications for the existence of a complementary relationship, although weaker than the one originally proposed by Bouchet (1963).
However, during the same period an alternate theory appeared.Stanhill and Cohen (2001) showed that the solar irradiance had been declining the past decades (global dimming), with potential influences on the evaporation.Cohen et al. (2002) in an Israel case study, argued that the known sensitivity of pan evaporation to net radiation at the surface is enough to explain its decrease in a globally dimming world.They proposed that aerosol and cloud-induced global dimming is the main reason for the general downward trend of the pan evaporation, because the vapour pressure deficit (VPD) was displaying increasing trends, contrary to the expectations of the complementary hypothesis.Roderick and Farquhar (2002) supported the global dimming solution to the paradox, showing that the recent solar flux decrease is enough to account for the pan evaporation trend in a former Soviet Union area.Linacre ( 2004) in a study with simplified global average changes of temperature, dew point temperature, solar radiation, agrees that global dimming is the major factor in the decreases of pan evaporation.Wild et al. (2004) showed Global Energy Balance Archive (GEBA) data, indi-cating that the net shortwave and longwave radiation available for evaporation has been declining.In order to reconcile the decreasing radiative energy with the increasing temperature observations, they conclude that the hydrological cycle (and thus the evaporative cooling) has to weaken.Therefore, they disagree with the complementary hypothesis, which requires enhanced evaporation.The global dimming has reversed to brightening in the 1990s, and according to Wild et al. (2008), Wild (2009) and references therein, the hydrological cycle shows signs of transition from weakening to strengthening at the same time.
To further complicate the picture, case studies in Australia (Roderick et al., 2007) and the Tibetan plateau (Zhang et al., 2007) showed that neither radiation trends, nor humidity issues were the major factor in pan evaporation trends.Instead the authors attribute the change to wind speed decreases.Johnson and Sharma (2010) calculated the pan evaporation trend in Australia from station, reanalysis and GCM data.Although the station data identified the wind change as the factor with the strongest contribution to the trend, the reanalysis data attributed the trend mostly to water vapour deficit change.However, it still remains to be seen if the effect of the wind has a regional or even global character.Brutsaert (2006) has proposed that the reported decreases in pan evaporation can be attributed partly to solar dimming and partly to actual evaporation increases via the complementary relationship.In other words the two hypotheses do not have to be mutually exclusive.Teuling et al. (2009) arrive to similar conclusions.They have found that evaporation depends on different drivers, in regions of Europe and North America.In drier water-limited regions, such as South Europe and the US Southwest, evaporation follows the interannual fluctuations of precipitation, while in wetter energylimited regions, such as central and North Europe and the American Northeast, evaporation follows global dimming and brightening.Therefore, evaporation can be increasing or decreasing with decreasing pan evaporation, depending on the location.For example, global dimming in energy-limited areas, will cause both pan evaporation and actual evaporation to decrease.On the other hand, a positive precipitation trend in water-limited regions, will increase actual evaporation, but decrease pan evaporation.
A recent review by Fu et al. (2009) on the subject lists the significant problems associated with each of the two explanations of the paradox, i.e. solar dimming and complementary hypothesis.In order to address these problems, they propose five priorities.First, homogenization of the data sets or quantification of the errors.Second, investigation on which of the three major factors (net radiation, vapour pressure deficit, and wind) is most significant in a regional context.Third, extension of the analysis to post 1990 yr, when a reversal from global dimming to global brightening is noted.Fourth, the relationship between pan evaporation and reference potential evaporation needs to be further clarified.Fifth, the separate study of land and ocean evaporation trends.Finally, they highlight that the crux of the problem is the trend of the actual evaporation and we need to find ways to advance our knowledge there.
In this work, we focus on the potential (open-water) evaporation rate over all land areas of the planet.In order to estimate it, we assume the existence of a small shallow water body in each land location.The hypothetical water body has to be small, so that the regional climate can be considered undisturbed, and shallow, therefore heat storage considerations can be ignored.The open water evaporation from this small, shallow water body can be assumed to be the (apparent) potential evaporation at the same location.Our objectives here are the accurate calculation of potential evaporation over land areas, derivation of its regional trends, and quantification of the contribution of solar decadal fluctuations and water vapour transfer changes to the potential evaporation trends.The link of potential evaporation to the pan evaporation will provide some insight to the pan evaporation decadal changes.
In more detail, we will analyse the monthly potential evaporation for years 1983-2008, in a spatial 2.5 • ×2.5 • resolution, over all land areas of the globe.We use our radiation transfer models for solar and thermal longwave radiation, already employed in a variety of climatic studies (Hatzianastassiou andVardavas, 1999, 2001a,b;Hatzianastassiou et al., 1999Hatzianastassiou et al., , 2004aHatzianastassiou et al., ,b, 2005Hatzianastassiou et al., , 2007a,b;,b;Pavlakis et al., 2004Pavlakis et al., , 2007Pavlakis et al., , 2008;;Vardavas and Taylor, 2007;Matsoukas et al., 2010) for the calculation of the necessary radiation fluxes.We also employ reanalysis humidity, wind and temperature data, in order to calculate the evaporation component due to mass transfer processes.We are thus in a position to include both parameters in our analysis: net radiation fluxes and water vapour transfer factors.Moreover, out of the five priorities proposed by Fu et al. (2009), our approach addresses priorities one, two, and three.Priority one, because we partly homogenise the data, using the same global data sources and same methodologies.Priority two, because all proposed quantities are included in our trend analysis.Priority three, because our period of study extends well beyond 1990.In Matsoukas et al. (2005Matsoukas et al. ( , 2007)), we touched upon priority five, but due to lack of inter-annual, spatially distributed heat storage data in the oceans, the trend of general ocean evaporation is beyond our reach.

Methodology and data
Our objective is to estimate the monthly potential evaporation globally over all land areas in 2.5 • ×2.5 • resolution.In each 2.5 • ×2.5 • cell the aerodynamic evaporation rate E a (in m s −1 units) is estimated from where U is the scalar wind, ρ the water density, e s the saturation water vapour pressure at the 2 m temperature, and e the observed water vapour pressure at 2 m.The difference e s − e is the vapour pressure deficit (VPD).C w is a turbulent exchange coefficient, estimated by a variety of methods, e.g.Brutsaert (1982, chap. 4), Winter et al. (1995).The derivation of C w is relatively simple for neutral atmospheric stability conditions, but once atmospheric instability is included in the analysis, it becomes considerably more complex.Neglecting air stability issues in the estimation of potential evaporation can lead to serious errors when the time resolution is finer than 24 h.However, for resolutions coarser than daily, the errors tend to cancel out and the assumption of neutral stability works well (Mahrt and Ek, 1993;Brutsaert, 1982).Even though our approach uses monthly values, we take into account atmospheric instability using Monin-Obukhov similarity theory to derive C w .Our approach is based on Beljaars and Holtslag (1991) for the calculation of Monin-Obukhov theory functions and on Brutsaert (1982) for the surface roughness parametrisation.The aerodynamically derived evaporation E a takes into account only the drying power of air.It requires as input few and easily available physical quantities, i.e. air temperature, humidity, wind speed, air pressure.All these data were taken in a monthly, 0.5 • ×0.5 • resolution from the European Centre for Medium Range Forecast (ECMWF) Re-Analyses, namely ERA-40 (Uppala et al., 2005) up to August 2002 and the latest ERA Interim for the period January 1989-June 2008.The main objectives of ERA Interim were "to improve on certain key aspects of ERA-40, such as the representation of the hydrological cycle, the quality of the stratospheric circulation, and the handling of biases and changes in the observing system" (Berrisford et al., 2009), including wind measurements, which is of specific interest in this study.ERA data were regridded to 2.5 • ×2.5 • resolution, in order to match the radiation transfer model resolution.
The energy balance method estimates the available energy R n for the turbulent fluxes (evaporation and sensible heat) using the principle of energy conservation: where Q s is the net solar energy flux, Q l the net terrestrial flux corresponding to radiative cooling, G is the energy flux advected away from the surface, and H is the stored heat.Since our model is applied to a shallow water body, H can be neglected (Brutsaert, 1982).Also, for monthly resolutions G is relatively small and can be considered negligible (Shuttleworth, 1993).
The energy balance evaporation rate E r (in m s −1 units), which corresponds to R n is where L its latent heat of evaporation, estimated by with the air temperature T a in • C. The energy balance evaporation rate E r , is a biased estimator of evaporation rate, because it neglects the sensible heat flux.However, E r is a necessary step and is correctly used in the framework of Penman's method, which will be presented below.Chow et al. (1988) state that "the evaporation may be computed by the aerodynamic method when energy supply is not limiting and by the energy balance method when vapour transport is not limiting.But, normally, both of these factors are limiting so a combination of the two methods is needed".Doing just that, the method developed by Penman (1948) gives the evaporation rate E p as a weighted average of E r and E a , using the formula where is the slope of the saturation water vapour pressure curve at the air temperature T a , and γ is the psychrometric constant with c p the moist air heat capacity and p the atmospheric pressure.
Penman's method has consistently ranked among the best methods for the calculation of potential evaporation over water bodies.Chow et al. (1988) classified it as the best evaporation method when all relevant data are available and the necessary assumptions are justified.Winter et al. (1995), Rosenberry et al. (2004Rosenberry et al. ( , 2007) ) found it one of three best out of at least eleven methods in the cases of Lake Williams in Minnesota, a prairie wetland in North Dakota, and a small mountain lake in the northeastern USA, respectively.Their comparison baseline was the energy balance method.Tanny et al. (2008) in a case study of a small reservoir in northern Israel ranked the Penman method between the two best out of five evaporation models.Their comparison was performed against eddy covariance evaporation measurements.However, its application requires many physical quantities, some of them not readily available.The ones that are not measured at every meteorological station are the radiation fluxes Q s and Q l .Therefore, using radiation transfer models to calculate them everywhere on the globe, is our only recourse.

Radiative transfer model description
The deterministic 1-D spectral radiative transfer model used here was developed from a radiative-convective model (Vardavas and Carver, 1984;Vardavas and Taylor, 2007).The sky is divided into clear and cloudy fractions.The cloudy fraction includes three non-overlapping layers of low, mid and high-level clouds.The model input data include cloud amounts (for low, mid, high-level clouds), cloud scattering and absorption optical depths, cloud-top pressure and temperature (for each cloud type), cloud geometrical thickness and vertical temperature and specific humidity profiles.For the total amount of ozone, carbon dioxide, methane, and nitrous oxide in the atmosphere, we used the same values as in Hatzianastassiou and Vardavas (2001a).
All of the cloud climatological data for our radiation transfer model were taken from the International Satellite Cloud Climatology Project (ISCCP-D2) data set (Rossow and Schiffer, 1999), which provides monthly means for 72 climatological variables in 2.5 • ×2.5 • , monthly resolution, for 15 cloud types and for the 25-yr period July 1983-June 2008.ISCCP converts 30 km × 30 km cloud data every 3 h to an equal-area map grid with 280 km resolution.The Stage D2 data product is produced by further averaging over each month, first at each of the eight 3 h time slots and then over all time slots.The cloud-top temperature is derived from the infrared radiances, while the cloud-top pressure from the vertical temperature profile of the atmosphere.We decided to use ISCCP and radiation transfer models for the calculation of radiative fluxes and not reanalyses, because the latter depend strongly on model generated clouds.Using models to create clouds introduces large uncertainties in GCMs and reanalyses.On the other hand, the use of satellite-observed cloud properties, such as the ISCCP dataset, saves us from the errors introduced by the cloud microphysical models.
The water vapour and temperature vertical atmospheric profiles, used in the radiative transfer model, come from the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) global reanalysis project (Kistler et al., 2001), corrected for topography as in Hatzianastassiou and Vardavas (2001a).These data are also on a 2.5 • resolution, monthly averaged and cover the same 25-yr period as the ISCCP-D2 data.
Our radiation transfer model is accordingly customised and applied in two different radiation domains, the shortwave (SW, solar) and the longwave (LW, terrestrial).Below, we present briefly the two separate models.

Shortwave radiation transfer model
The incoming solar irradiance conforms to the spectral profile of Thekaekara and Drummond (1971) and corresponds to a solar constant S 0 of 1367 Wm −2 (Willson, 1997;Hartmann, 1994).The model makes adjustments for the elliptical Earth orbit and apportions 69.48 % of the incoming spectral irradiance to the ultra violet-visible-near infrared (UV-Vis-NIR) part (0.20-1 µm) and 30.52 % to the near infraredinfrared (NIR-IR) part (1-10 µm).Then, the radiative transfer equations are solved for 118 separate wavelengths for the UV-Vis-NIR part and for 10 bands for the NIR-IR part, using the Delta-Eddington method of Joseph et al. (1976).For a more detailed model description the reader is referred to Hatzianastassiou et al. (2004aHatzianastassiou et al. ( ,b, 2007a,b),b).
The model takes into account Rayleigh scattering due to atmospheric gas molecules, as well as absorption from O 3 , CO 2 , H 2 O, and CH 4 .The O 3 column amount is taken from the Television Infrared Observational Satellite (TIROS) Operational Vertical Sounder (TOVS).Complete aerosol data are provided by the Global Aerosol Data Set (GADS) (Köpke et al., 1997).
The model output can include downwelling and upwelling fluxes at the top of atmosphere, at the surface and at any atmospheric height.The focus of this study is the solar energy absorbed at the surface, or net (downwelling minus reflected) flux at the surface.As mentioned above, in this study we use the fluxes absorbed by a hypothetical shallow water body on the ground.Therefore, in the estimation of the net shortwave flux we use the water surface albedo, modelled using Fresnel reflection as a function of the solar zenith angle and corrected for a non-smooth surface.

Longwave radiation transfer model
The detailed radiative-convective model developed for climate change studies of Vardavas and Carver (1984) is modified for the radiation transfer of terrestrial infrared radiation, in order to compute the downwelling longwave radiation (DLR) and upwelling fluxes at the surface of the Earth.The model has monthly, 2.5 • ×2.5 • resolution (dictated from the resolution of ISCCP-D2) and a vertical resolution of 5 mb, from the surface up to 50 mb, to ensure that the atmospheric layers are optically thin with respect to the Planck mean longwave opacity.The skin temperature, as well as the humidity and temperature vertical atmospheric profiles come from the NCEP/NCAR global reanalysis project.We have assumed that the skin temperature of the theoretical shallow body of water is the same as the skin temperature of the soil, as given by the NCEP/NCAR reanalysis.This is probably a good approximation, if indeed the water body is small and shallow, and the evaporative cooling is weak.In any case, in a sensitivity analysis we reduced drastically the skin temperature by 10 • C over the Sahara and the Arabian peninsula.The weaker radiative cooling resulted in an E p increase about 15 % over these regions.The global land E p increased by 3 %.However, the trends analysis presented later, is unaffected.
The atmospheric molecules considered are H 2 O, CO 2 , CH 4 , O 3 , and N 2 O.The sky is divided into clear and cloudy fractions.The cloudy fraction includes three nonoverlapping layers of low, mid and high-level clouds.The model input data include cloud amounts (for low, mid, highlevel clouds), cloud scattering and absorption optical depths, cloud-top pressure and temperature (for each cloud type), cloud geometrical thickness and vertical temperature and specific humidity profiles.
A full presentation and discussion of an earlier version of the model can be found in Pavlakis et al. (2004).There, a series of sensitivity tests were performed to investigate how much uncertainty is introduced in the model DLR by uncertainties in the input parameters, such as air temperature, skin temperature, low, middle or high cloud amount as well as the cloud physical thickness, cloud overlap schemes, and the use of daily-mean instead of monthly-mean input data.The model DLR was also validated against BSRN station measurements for the entire globe (Pavlakis et al., 2004;Matsoukas et al., 2005).Here, we use a newer version with an improved spectral resolution of 28 bands for the longwave spectrum.

Long-term average
We start by using Eq. ( 1) to compute the bulk aerodynamic evaporation E a for the globe, with input data (air temperature, humidity, wind speed, air pressure) originating from ERA-40 and ERA Interim.This quantity corresponds only to mass transfer procedures and assumes unlimited energy availability.A long-term average (July 1983-June 2002) of the annual E a , calculated from ERA-40 data is presented in Fig. 1 (top left), showing maxima over generally dry areas, such as deserts and minima over wetter and colder areas.A direct comparison with the long-term annual average of the VPD in Fig. 1 (top right), shows the qualitative resemblance with E a and highlights the dominance of the VPD in the geographical distribution of aerodynamic evaporation.The long-term annual average of wind speed is shown in Fig. 1 (bottom).The regional patterns of U and E a do not correlate very well, indicating that the wind speed U plays only a secondary role in the bulk aerodynamic evaporation regional distribution.This is also true for the the exchange coefficient C w .
We proceed by calculating the available energy R n for turbulent processes, i.e. evaporation and sensible heat flux.As mentioned before, this energy is derived from our radiative models, run over fictitious small and shallow water bodies at each location.If we assume that all this energy flux is used up in evaporation, we obtain the evaporation rate E r , whose global distribution is shown in Fig. 2 (top left).A comparison of E a in Fig. 1 (top left) with E r in Fig. 2 (top left) shows that there are two very distinct processes that govern potential evaporation.In some dry locations E a is larger than E r , meaning that the available energy does not suffice to maintain the potential evaporation rate dictated by mass transfer and potential evaporation is energy limited there.This is highlighted in Fig. 2 (top right), where we present the ratio E a /E r .In areas where this ratio increases above one, e.g. the Sahara and Australian deserts, potential evaporation is limited by energy.In other locations, e.g.North Eurasia and America, South America, E r is larger than E a and the ratio drops below one, meaning that the potential evaporation there is vapour transfer limited.
Penman's method (Eq.5) takes into account the two processes and derives an accurate estimate E p of the potential evaporation.The long-term average E p is shown in Fig. 2    (bottom), with its expected latitudinal gradient, its maximum values close to the Equator, its poleward decrease and its minimal values in Antarctica and Greenland.E p follows closely E r , because is usually larger than γ in Eq. ( 5) and therefore the E r contribution to E p dominates the E a contribution in all but the coldest regions.The global average of the potential evaporation of small shallow water bodies over land areas for July 1983 to June 2002 is 3.4 mm day −1 .

Global trends
We now focus on the interannual behaviour of global potential evaporation, as well as the quantities that affect it, namely wind speed, VPD, net solar flux and net terrestrial flux.In this study "global" means the area-weighted average of land-only values.We first present in Fig. 3a the evolution of global mean bulk aerodynamic evaporation E a , calculated by Eq. (1), using data from ERA-40 (blue line) and ERA Interim (red line).There is a consistent difference between the two time series during the temporal overlap between January 1989 and August 2002, when data were available from both ERA-40 and Interim.The reason for the larger Interim values is the consistently increased wind speeds compared to ERA-40, as can be seen in Fig. 3b.Except from the wind speed, the global means of the other relevant quantities (2 m air temperature, dew point temperature, VPD) do not show significant differences between ERA-40 and Interim.For example, in Fig. 3c we show the time series of the VPD global mean, where the two lines are very close.
Fitting linear trends to raw time series, such as the ones in Fig. 3, is not a good practice, because the large seasonal fluc-tuations hide possible trends and keep them from attaining statistical significance.Moreover, the switch from ERA-40 to ERA Interim after 2002 would introduce a spurious trend, especially if the same quantity has different long-term mean in the two data bases, as for example wind speed.In order to perform a more robust trend analysis, we take the deseasonalised time series, normalised by the standard deviation of the interannual variability for each specific month.For example, the normalised anomaly of June 2001 is the difference between the June 2001 value and the mean of all June values, divided by the standard deviation of all June values.In this fashion we have derived normalised plots for wind speed U , VPD, and bulk aerodynamic evaporation E a (Fig. 4), net solar flux Q s , net terrestrial flux Q l , and energy balance evaporation E r (Fig. 5), and Penman potential evaporation (Fig. 6).
As expected, the seasonality of the quantities is not apparent in these plots and we can proceed to check if decadal trends are present globally over land.Also, the differences between normalised anomalies of ERA-40 and Interim for quantities U , VPD, E a , and E p (Figs. 4 and 6), seem small and non-systematic.If we compare the trends calculated separately by ERA-40 and Interim in their common time period (January 1989-August 2002), we see that they are very similar, a finding also visually supported by the proximity of the blue and red lines of Figs. 4 and 6.Therefore, there is evidence that using ERA-40 even before 1989 to calculate trends does not produce significant errors.For each one of the aforementioned quantities, we generated a "blended" time series with only ERA-40 normalised anomalies for July 1983-December 1988, the average of ERA-40 and Interim  normalised anomalies for the overlapping period of January 1989-August 2002, and ERA Interim only for September 2002-June 2008.These "blended" time series will be used for the rest of the paper in order to calculate linear trends for the relevant variables U , e s − e, E a , and E p .
There seems to be a wind speed decrease in the 90's, which has stopped (if not reversed) around 2000.The VPD seems to be increasing since the 90's, as does E a (but less steeply).Trends in Q s , Q l , E r , and E p are not very obvious and de-pend on the choice of the period of interest and its start and end points.Q l in particular displays an increasing trend in the last two decades.
Fitting linear trends to the generated normalised anomaly time series produces the trend values shown in Table 1.The numbers in italics mean that the 95 % confidence interval of the slope does not contain zero.The normalised anomalies have also been averaged over North and South Hemispheres and their trends are presented separately.Restricting ourselves to the statistically significant trends at the 95 % confidence level, we can draw some conclusions.
The VPD has the steepest increasing trend, already apparent in Fig. 4b, driving the bulk aerodynamic evaporation E a trend to positive statistically significant values.Although Q s shows a positive trend in the late eighties and throughout the nineties (solar brightening), after 2000 this trend has levelled off (Wild et al., 2009).In the full period (July 1983-June 2008) there remains a weak but statistically significant increase in the net solar heating of the surface.The radiative cooling of the surface Q l is increasing more steeply than Q s .
The same sign of trends in solar heating and thermal radiation cooling cause a statistically non significant change in the net energy flux at the surface, therefore E r has a weak, non significant trend.However, our result for the apparent potential evaporation E p over land areas, estimated by Penman's method, is a statistically significant increase over the last decades.The ratio E a /E r is also increasing with statistical significance, indicating a potential evaporation that is less limited by vapour transfer, but more limited by energy fluxes.In other words, the VPD has increased, resulting in a more vapour-hungry atmosphere, but the net energy at the surface has not increased as much as to satisfy this demand.
All above trends are stronger in the North Hemisphere than either globally or in the South Hemisphere.In SH the only statistically significant trends are found for VPD and E a .On the other hand, in NH all examined quantities seem to be significantly changing, with the wind speed decreasing and everything else increasing.
It is physically more meaningful to present the trends of the actual physical quantities than the trends of the normalised anomalies, even though the latter can be derived with smaller statistical errors.A first set of quantities (Q s , Q l and E r ) originate from a single, consistent data set for the whole period July 1983 -June 2008, in contrast with a second set (U , e s − e, E a , E p ) coming from a blend of ERA-40 and Interim data.For the first set with quantities from only one data source, the trend of the original timeseries can be produced safely (Table 2).We also show the trends of the original timeseries of the second set separately for ERA-40 and ERA Interim, because merging the two sets is problematic when we deal with the actual (and not the normalised) timeseries.Our trends seem rather small compared to the ones derived by Wild et al. (2008).However, they correspond to a different period and for comparison purposes we would need to shorten our temporal window from 1983-2008 to 1986-2000.If we do that, the trends can change significantly, e.g. the Q s trend increases by an order of magnitude.
There is a general reluctance to derive trends from reanalysis data.The main reason is that the observational data assimilated by the reanalysis scheme may at some times be less dense than other times, or that the assimilation scheme may be revised.These changes could produce spurious, statistically significant trends.However, we decided to use reanalysis data for the calculation of E p trends for two reasons.First, ERA Interim starts in 1989, when satellite data were mature, global and abundant, leading us to expect small changes in data coverage and density.Second, reanalysis data trends affect mainly E a and not E r .As we note above, E p depends weakly on E a , except in the very cold regions.Therefore, spurious reanalysis trends, will affect E p only www.atmos-chem-phys.net/11/7601/2011/Atmos.Chem.Phys., 11, 7601-7616, 2011  slightly.We should mention that there are many reanalyses available, as well as many surface radiation flux datasets, which sometimes show different trends in the physical quantities related to potential evaporation.It is encouraging to see good agreement between the U , VPD, E a , and E p normalised anomaly trends produced with ERA-40 and Interim data during the common period (January 1989-August 2002) in Figs. 4 and 6).A thorough comparison between many more datasets and a sensitivity study with respect to the estimation of potential evaporation trends would be very useful, but is out of the scope of this work.Also, we cannot preclude the existence of spurious trends in our radiation model inputs, such as the ones reported by Evan et al. (2007) for ISCCP.However, radiation flux trends derived by our earlier works compare well with observations from ground stations and satellites (Hatzidimitriou et al., 2004;Hatzianastassiou et al., 2005;Fotiadi et al., 2005).With these caveats in mind, we continue our analysis.

Geographically resolved trends
We compile monthly timeseries at every 2.5 • × 2.5 • cell of the globe for the period July 1983-June 2008.Each individual timeseries is transformed to a normalised anomaly timeseries and fitted with a least-squares line.For quantities originating from both ERA-40 and Interim, a blended normalised anomaly time series has been calculated and all plots and results reported here are derived from it.Each slope is tested for difference from zero with statistical significance at 95 % confidence level.Below, we present the regional distribution of these trends for potential evaporation and all other relevant physical quantities.Only statistically significant slopes are presented.The normalised anomalies are unitless.An interpretation of the values of quantities in the following trend plots (Figs. 7 and 8) is the change per decade of the difference of the quantity from its interannual monthly mean over its interannual monthly standard deviation.
Examining the regional wind trends, it is apparent that in the last 25 yr wind speeds are generally decreasing everywhere except Antarctica, scattered parts of North America and East Asia, extended parts of South America and North Africa, central Europe, the Kalahari, Indochina and Indonesia.These trends are in agreement (at least qualitatively) with Roderick et al. (2007), Zhang et al. (2007), and McVicar et al. (2008).McVicar et al. (2008) point out that Australian stations and ERA-40 reanalysis data agree well in wind climatologies, but the reanalysis trends are weaker than the station ones.Our ERA analysis is also in agreement with Pryor et al. (2009), who however find differences in reanalysis and station trends for the contiguous US and show scepticism in using reanalysis data for trend detection.There is also some similarity with the the investigation of station data by Vautard et al. (2010), because renalyses such as ERA Interim capture the large-scale circulation changes.Differences with Vautard et al. (2010) can be attributed to processes such as land-use changes, which are not taken into account by reanalyses.We understand the difficulties in substituting local wind measurements with reanalysed data, but we should also note that there is no consensus on the magnitude of wind trends in Australia, even between station-only analyses (Jovanovic et al., 2008).
The bulk aerodynamic evaporation E a and the VPD (i.e. e s − e) are globally on the rise except (not necessarily with statistical significance) in India, West Australia, South Africa and various scattered small regions throughout the globe.The VPD trend is caused by the fact that in ERA both air temperature T and dew point temperature T d are rising, but T is doing so faster than T d (results not shown).The positive trend of T d shows that the water cycle is intensifying, but not enough to decrease the VPD in a warming world.This increase in the "drying capacity" of air runs contrary to the complementary hypothesis.
The globe seems to be divided in half with respect to the prevalence of global dimming or global brightening in the period July 1983-June 2008.Decadal changes in net solar heating Q s led to less radiation (dimming) in parts of Canada, Greenland, North Eurasia, East Asia, most of Australia, and Antarctica (not shown, but very similar to Fig. 8 (top)).All other areas witnessed brightening.The surface seems to be radiating increasingly net longwave flux Q l during the examined period, contributing to less available energy for evaporation.Exceptions are West South America, West Africa, West Australia, Central Eurasia, Indonesia and scattered spots in various locations (not shown).The combined trends of Q s and Q l produced the trend of E r , shown in Fig. 8 (top).The general image is similar to the trends in Q s , indicating that the major player in radiative flux changes is the solar energy and not the terrestrial longwave.Also, Penman potential evaporation E p regional trends in Fig. 8 (bottom) agree quite well with E r in Fig. 8 (top).In some regions, opposite signs of trends in E r and E a have removed statistical significance in E p changes, e.g.parts of Greenland and North-east Asia, but the differences between the two figures are small.
In the Introduction we highlighted that generally the pan evaporation in observations has been decreasing, while our results so far show a general increase in potential evaporation.This disagreement is based on the fact that the majority of our referenced observations correspond to data before the early 1990s, when dimming was worldwide still prevalent.Things are quite different in the late 1990s, when potential evaporation increased globally, and in the 2000s, when it slightly decreased.We compare here our results with recent pan evaporation observations.Roderick et al. (2009)   and for different time periods in their Table 1.We select from there the trends that extend to 2000 and later, and compare them with the (not necessarily statistically significant) normalized trends of our E p , in Table 3.This comparison aims to provide a first outlook of our model behaviour and is not meant to be a point-by-point quantitative comparison.
We did not go into details such as the exact location of the pans.There are countries with opposing regional trends, but we assigned one trend per country taking into account the geographically prevalent trend.With respect to the observational time series, we should note that one starts as early as the 1950's, while almost all start in the 1970's or before.Our modelled potential evaporation time series begins in 1983, so the observational time series contain a dimming period of at least one decade, which is not there in our E p .We do not have access to the observational data and we cannot isolate the post-1983 observational data.Unfortunately, a more meaningful comparison of contemporaneous modelled and observational time series cannot be performed.With the above in mind, only two out of seventeen studies do not agree with our results, both of them coming from analysis of only one site: (1) in Turkey, where Roderick et al. (2009) state "This pan was located in an expanding irrigation area", and (2) in Ireland, where we have two studies agreeing and one disagreeing with us.These recent observations extending at least to 2000, tend to correspond geographically with areas where we find decreasing potential evaporation.Our comparison would be more complete if we had a wider coverage including more regions with increasing potential evaporation.
However, the large majority of comparisons in Table 3 shows generally good agreement of observational studies with our modelling approach.
4 Discussion: is potential evaporation driven by energy or mass transfer issues?
Let us revisit the two competing hypotheses on the decreasing pan evaporation trends.On the one hand the complementary hypothesis proposes that changes in the water vapour mass transfer are responsible for the observed pan evaporation trends.On the other hand, the secular global dimming/brightening is supposedly enough to account for the pan evaporation trends.In order to provide some insight into this problem, it is important to elucidate the relationships between the deseasonalised, normalised trends of U , e s −e, E a , Q s , Q l , E r and E p .To this end, we compile Table 4 with the cross-correlation coefficients R 2 between all possible pairs of the above quantities.We will examine the values greater than 0.5, namely the ones corresponding to the pairs E and E p -E r .The large correlation of E a and VPD has been inferred previously in this study, by the similarities in both their regional distributions and their global normalised anomalies.Q s and Q l are correlated because of their common relationship to cloud cover.For example, when clouds are present, they obstruct both net solar heating and terrestrial cooling, leading to smaller values of  2) and ( 3), the close relationship of E r with Q s is not surprising.Finally, it is expected that E p would be correlated with E r and E a through Penman equation (Eq.5).We have already pointed out the stronger relationship of E p with E r , rather than with E a , due to the usually larger coefficient of E r in Eq. ( 5).Notably, with a 0.90 correlation coefficient, we can see that Penman's potential evaporation E p trend is strongly dictated by the trend of the available energy E r .The E p -E a correlation is important, but weaker than this of the pair E p -E r .This result is examined further, by presenting the geographically distributed cross-correlation between deseasonalised and normalised trends of E p on the one hand and E r or E a on the other, in Fig. 9.The very strong crosscorrelation found for the global normalised, deseasonalised trends of E p and E r (top of Fig. 9) in the previous paragraph seems very robust.In most areas this correlation coefficient is above 0.7, with exceptions at most deserts, parts of Canada and Siberia, Greenland, and Antarctica.E p trends are dominated by E a trends only in West Australia, Sahara, Kalahari, Greenland, Antarctica, parts of Russia and North America, while everywhere else in the globe E r is more important.E a trends are dominant in the deserts, because E a is larger than E r there, and in the coldest areas, because the coefficient of E a is larger there than the coefficient of E r in Eq. ( 5).This result can be extended to pan evaporation, since potential and pan evaporations are generally considered proportional, with "pan coefficient" being the coefficient of proportionality for a specific pan.
This finding comes to the support of the hypothesis that secular dimming/brightening controls the global trends in potential and pan evaporation.The bulk aerodynamic evaporation E a is also positively correlated with the potential evaporation E p , but this effect is of principal importance only in deserts and the coldest regions.The energy balance evaporation E r is not totally independent of E a , so all three quantities are linked.However, the degree of E r and E p trend correlation is prominent and provides a clear answer.In this study, at most locations the available energy (mainly solar) is the driver for potential (and consecutively pan) evaporation trends in our dataset and not water vapour transfer.

Conclusions
Pan and potential evaporation are directly proportional and therefore a study on potential evaporation can shed some light on the factors that affect the much discussed observed trends of pan evaporation.We use Penman's method to calculate the potential evaporation for all land areas of the globe in a monthly 2.5 • × 2.5 • resolution.Penman's method is widely recognised as one of the most accurate calculations of potential evaporation.It takes into account two relevant processes: the drying power of the air (the water vapour transfer potential) and the energy available to the evaporation process.Penman's method then compromises the two sometimes conflicting processes and produces an estimate for potential evaporation.Since it requires modelling of both radiative and turbulent fluxes, it is data intensive and needs quantities (e.g.net radiation fluxes), which are not readily measured in all world regions.We employ radiation transfer models as a way to circumvent the problem of radiative flux www.atmos-chem-phys.net/11/7601/2011/Atmos.Chem.Phys., 11, 7601-7616, 2011 availability.In order to run the radiative transfer models we need data for clouds (ISCCP-D2), aerosol (GADS), and atmospheric temperature and humidity profiles (NCEP/NCAR reanalysis).ECMWF reanalysis data (ERA-40 and ERA Interim) are used for the turbulent flux model.The temporal evolution of potential evaporation for our examined period (July 1983-June 2008) is of particular interest because it provides insight to the observed trends of pan evaporation.We examine the decadal trends of potential evaporation and various related physical quantities, such as net solar flux, net longwave flux, vapour pressure deficit (VPD) and wind speed.Trends of reanalysis quantities, which are related to turbulent fluxes, such as the VPD, wind, and bulk aerodynamic evaporation E a , should be reported carefully because of possible changes in data coverage and assimilation techniques during the reanalysis period and the consequent generation of spurious trends.In our case, we think spurious effects are minimal because data coverage changes are small after 1989, with satellites routinely providing uninterrupted global coverage.Moreover, Penman potential evaporation depends weakly on turbulent fluxes but more strongly on radiative fluxes, which are less sensitive to reanalyses.
Atmospheric temperature, dew point temperature, and vapour pressure deficit (VPD) appear to be globally increasing.The increase of dew point temperature is an indication of an enhanced water cycle, in line with the complementary hypothesis.However, the agreement stops here, because the temperature is increasing even faster than the dew point temperature, resulting in decreasing relative humidity and increasing VPD and E a .
Global dimming/brightening cycles in these 25 yr increased the available solar energy, driving the total energy available for evaporation to larger values.Radiative longwave cooling also increased, having the opposite to solar effect on total available energy.The energy balance evaporation E r is proportional to this available energy for evaporation and therefore it has a positive trend.The ratio E a /E r is also increasing with statistical significance, indicating a potential evaporation that is less limited by vapour transfer, but more limited by energy fluxes.All above trends are stronger in the North Hemisphere than either globally or in the South Hemisphere.In SH the only statistically significant trends are found for VPD and E a .On the other hand, in NH all examined quantities seem to be significantly changing, with the wind speed decreasing and everything else increasing.
Potential evaporation trends seem to follow more closely the trends of energy availability and less so the trends of the atmospheric capability for vapour transfer.This finding is geographically rather robust, with the exceptions of deserts and the coldest regions, such as Antarctica, Greenland, most of Australia, the Sahara, the Kalahari, and parts of N. America and Russia.The results above tend to support the hypothesis that secular changes in the radiation fluxes are responsible for global potential evaporation trends, and not vapour trans-fer issues, such as the ones proposed by the complementary hypothesis.

Fig. 1 .
Fig. 1.Long-term averages (July 1983-June 2002).Top left: the bulk aerodynamic evaporation rate E a (mm day −1 ) computed from wind speed, air temperature and humidity data from ERA-40.The E a values reach as high as 16 mm day −1 over the Sahara.Top right: the VPD (mb) computed from air temperature and humidity data from ERA-40.Bottom: the wind speed U (m s −1 ) from ERA-40.

Fig. 2 .
Fig. 2. Long-term averages (July 1983-June 2002).Top left: energy balance evaporation rate E r (mm day −1 ) computed by radiation transfer models.Top right: ratio E a /E r .Bottom: Penman evaporation rate E p (mm day −1 ), as a weighted average of E r and E a .

Fig. 3 .
Fig. 3. Global means of: (a) bulk aerodynamic evaporation E a (in mm day −1 ) calculated from input data by ERA-40 (blue line) and ERA Interim (red line).(b) 10 m wind speed U (in m s −1 ) by ERA-40 (blue line) and ERA Interim (red line).(c) VPD (in mb) by ERA-40 (blue line) and ERA Interim (red line).

Fig. 4 .
Fig. 4. Normalised anomalies of (a) 10 m wind speed U by ERA-40 (blue line) and ERA Interim (red line) (b) VPD by ERA-40 (blue line) and ERA Interim (red line) and (c) bulk aerodynamic evaporation E a , computed by data from ERA-40 (blue line) and ERA Interim (red line).

Fig. 5 .
Fig. 5. Normalised anomalies of (a) net solar flux Q s , (b) net terrestrial flux Q l , and (c) energy balance evaporation E r , all calculated by our radiation transfer model.

Fig. 6 .
Fig. 6.Normalised anomaly of Penman potential evaporation, computed by data from our radiation transfer model, and ERA-40 (blue line) or ERA Interim (red line).

Table 1 .
Global, North Hemisphere (NH) and South Hemisphere (SH) normalised anomaly trends for potential evaporation over land areas and other relevant physical quantities, for the period July 1983-June 2008.The numbers are the slopes of Figs.4-6 (except the E a /E r slope) in units of decade −1 .Numbers in italics correspond to trend slopes significantly different from zero at the 95 % confidence level.

Fig. 7 .
Fig. 7. Statistically significant (95 % confidence interval) trends of the normalised anomaly of (top left) wind speed U , (top right) VPD e s − e, and (bottom) bulk aerodynamic evaporation E a .All three time series are blended from ERA-40 and Interim, at each 2.5 • ×2.5 • cell.Units are decade −1 .

Fig. 9 .
Fig. 9. Geographically distributed correlation coefficient between trends of E p and E r (top) and between trends of E p and E a (bottom).

Table 3 .
Roderick et al. (2009) of pan evaporation trends from Table1ofRoderick et al. (2009)with E p trends from our study.The first column has units of mma −1 a −1 ."√" in the last column means that the model reproduces the sign of the observed trend, while "×" means the opposite.

Table 4 .
Cross-correlation coefficients between potential evaporation E p and relevant quantities.The analysed timeseries are the normalised anomalies for the period July 1983-June 2008.l and Q s .The long-term global average values of Q s and Q l are respectively 150 W m −2 and 48 W m −2 , so Q s dominates Q l .If we also take into account Eqs. ( Q