Stronger warming amplification over drier ecoregions observed since 1979

Observations show that the global mean surface temperature has increased steadily since the 1950s and this warming trend is particularly strong and linear over land after 1979. This paper analyzes the relationship between surface temperature trends observed over land for the period 1979–2012 and enhanced vegetation index (EVI), a satellite measured vegetation greenness index, by large-scale ecoregion. The land areas between 50°S and 50°N are classified into various large-scale ecoregions based on the climatological EVI values. The regional mean temperature trends exhibit significant spatial dependence on the regional mean EVI. In general, the warming rate increases dramatically with decreasing EVI, with the strongest warming rate seen over the driest ecoregions. When anthropogenic and natural forcings are included, climate models are generally able to reproduce observed major features of the spatial dependence. When only natural forcings are used, none of the observed features are simulated. Furthermore, the simulated temperature changes in the latter are mostly far outside the range of those in the former. These results suggest stronger warming amplification over drier ecoregions in the context of global warming, pointing mainly to human influence.


Introduction
The global mean surface temperature has increased by about 0.85°C since 1880 and much of the warming has occurred in the last several decades (IPCC 2013). Both observations and model simulations have attributed the global warming since 1950 mostly to the increase in anthropogenic greenhouse gases (GHGs) in the atmosphere (IPCC 2013). The GHGs-induced positive radiative forcing is global-scale and roughly spatially uniform but can be amplified or dampened by various feedbacks at regional scales (Hansen et al 2010, Thorne et al 2011. For example, the warming is much greater over land than oceans, over higher latitudes (Lu and Cai 2010) and over higher elevations (e.g., Naud et al 2013.
Over non-polar land surfaces, the warming rate differs substantially among ecosystems due to their differences in thermal and hydrological properties.
The GHGs-induced positive radiative forcing at the surface can be largely weakened by evapotranspiration (ET) with increasing amounts of vegetation and soil moisture (Koster et al 2004, Zhou et al 2007, Jeong et al 2009, Georgescu et al 2011. Lim et al (2005) showed that the warming rate is larger for barren lands. Zhou et al (2007Zhou et al ( , 2009Zhou et al ( , 2010 found that the warming trend of minimum temperature and the decreasing trend of diurnal temperature range are stronger over drier regions. ET is a primary process driving energy and water exchanges between the hydrosphere, atmosphere and biosphere (Wang and Dickinson 2012) in regions where soil moisture is the main controlling factor for ET. With decreasing soil moisture content, particularly over very dry ecosystems, the soil moisture-ET coupling weakens dramatically. Thus, such spatial variations of ET may determine surface warming rates among different ecosystems. However, few studies have quantified such spatial variations in the warming rates over land.
Quantifying and attributing this spatial dependence is important for understanding climate sensitivity to anthropogenic forcings, uncovering mechanisms of spatiotemporal patterns of climate change, and assessing climatic impacts. Here we use the enhanced vegetation index (EVI), an optimized vegetation greenness index measured by the MODerate resolution Imaging Spectroradiometer (MODIS) satellite sensor (Huete et al 2002), to examine the spatial dependence of observed surface warming rates on large-scale ecoregions. Similar analyses are also performed to historical simulations of global coupled atmosphere-ocean general circulation models (AOGCMs) developed for the Coupled Model Intercomparison Project phase 5 (CMIP5) (Taylor et al 2012), which are generally able to reproduce the observed warming, especially after the 1950s (Simmons et al 2010, IPCC 2013. EVI does not saturate, even over dense forests, and correlates highly with ET, particularly at large scales (Suzuki andMasuda 2004, Nagler et al 2005). This study focuses only on the modern satellite data era for the period 1979-2012 to maximize spatial coverage of measurements. Furthermore, after adjusting observations to remove the estimated impacts of known factors on short-term temperature variations, the global warming signal from 1979 through 2010 is particularly strong and steadily linear (Foster and Rahmstorf 2011).

Data and methods
This study uses the most recent versions of global gridded monthly surface air temperature (T) datasets: CRU (Harris and Jones 2014), GISS (Hansen et al 2010) and NCDC (Vose et al 2012) for the period 1979-2012, which have been widely used for long-term T variability and trend analysis. Despite sharing some similarities in input data sources, these three datasets differ substantially in their data processing approaches. For example, satellite data is used extensively in GISS, used very limited in NCDC, and not used at all in CRU (Vose et al 2012). In addition, monthly data of GPCP precipitation (Adler et al 2003) for the period 1979-2012 and MODIS EVI (Huete et al 2002) for the period 2000-2012 are also used. More specifics of each dataset are detailed in supplementary table S1.
We analyze monthly output of historical simulations from the CMIP5 archives, which include timeevolving changes in anthropogenic (GHGs and sulfate aerosols) and/or natural (solar and volcanic) forcing agents for the period 1860-2006, extended for the years 2006-2012 with the RCP4.5 scenario runs. These simulations are divided into two groups: one with anthropogenic and natural forcings (referred to as ALL) and the other with only natural forcings (referred to as NAT). Although most of CMIP5 models have performed multi-member ensemble runs, we choose only one ensemble 'r1i1p1' from 12 AOGCMs (supplementary table S2) that have required variables in ALL, NAT and RCP45 simulations. As averaging over multiple members enhances the forcing signal and reduces noise from internal variability and errors from individual models, we simply average the 12 simulations to obtain the multi-model ensemble mean for the period 1979-2012 in ALL and the period 1979-2005 in NAT.
All variables are spatially re-projected into grid boxes of 2.5°latitude by 2.5°longitude. The monthly EVI and precipitation are aggregated to create the annual climatology, and the monthly temperatures are temporally averaged to generate annual mean anomalies. For a time series at a given grid box, a linear trend is estimated using least squares fitting and a two-tailed student's t test is used to test whether the trend differs significantly from zero. For simplicity, we refer to the temperature trend as T trend .
We focus only on the land areas that have at least 10 months of data for each year and at least 26 years of data (i.e., 75%) during the period 1979-2012 in all three observational T datasets following the selection criteria of Vose et al (2012). The land beyond 50°N and 50°S is excluded to minimize snow/ice-albedo feedbacks that dominate the high-latitude warming. Annual mean T anomalies are used to maximize large-scale long-term T trends as the GHGs-induced radiative forcing is likely a dominant driver on these scales (e.g., Zhou et al 2007Zhou et al , 2009Zhou et al , 2010. Seasonal or monthly T anomalies are also affected by other factors such as sea surface temperature (Cohen et al 2012) and cloud cover (Tang and Leng 2013). The three datasets show similar T changes and thus their ensemble mean is mostly used to reduce redundancy. In addition, only those grids with annual T trend that is statistically significant (p < 0.1) are considered to exclude regions with no trends or unrealistic trends due to data uncertainties or gaps. In total, among the 1537 land grid boxes of 2.5°by 2.5°between 50°S and 50°N, 1338 are considered as the study region. CMIP5 simulations are sampled so that coverage corresponds to that of the observations. This study evaluates T trend at the aggregated largescale level to minimize grid-level data noise and variability. We classify the 1338 grids into 7, 14, 21, 28, and 35 large-scale ecoregions from barren to dense vegetation based on the climatological EVI values (referred to as EVI ), and then analyze how T trend varies with EVI by ecoregion via least squares fitting. The goodness of fit (R 2 ) is used to measure how successful the fit is in approximating the fraction of the data variations. Different classifications are used to test whether the fitted T trend −EVI relationship is robust. For each classification, every ecoregion contains about the same number of grid boxes. The regional mean time series is calculated using area-weighted averaging over land grids within each ecoregion, and its trend is estimated as done at the grid level. For brevity, only results of 7 and 35 ecoregions, which represent the least and most ecoregions classified, are shown in figures, while those of other classifications are listed in tables. Figure 1(a) shows the observed T trend over all the 1537 land grids between 50°S and 50°N for the period 1979-2012, together with the climatology of EVI (figure 1(d)) and precipitation (figure 1(e)). T increased almost everywhere except for parts of Northern Australia. Among the 1537 land grids, 87 and 84% exhibit a linear trend that is statistically significant at the 10 and 5% level, respectively. The three T datasets of CRU, GISS and NCDC display very similar warming features in both spatial patterns and magnitude (not shown), except for several particular regions in the CRU data (e.g., Amazon and Central Africa) where limited observations were used to interpolate the gridded data (Harris and Jones 2014). Significant warming occurs mostly in arid and semiarid regions such as Northern Africa, Middle East, Northern China, and Western US. EVI (figure 1(d)) and precipitation (figure 1(e)) resemble each other geographically (with a spatial correlation of 0.90; n = 1338)-the region with the least amount of precipitation has the least amount of vegetation or vice visa. Overall the warming trend is generally strongest over the driest or least vegetated regions such as the Sahara desert and the Arabian Peninsula. Figure 2 illustrates the 7 and 35 large-scale ecoregions classified based on the climatology of EVI, representing the range of the least and most ecoregions considered in this study. As expected, most ecoregions are spatially coherent and resemble the geographic distribution of the climatology of precipitation (figure 1(e)). Figure 3 shows regional mean T anomalies from 1979 to 2012 for the three individual T datasets. Four ecoregions were chosen to represent the ones with the least and most amount of vegetation or the driest and wettest climate. In every ecoregion, the three T datasets exhibit similar interannual variability and comparable warming trends that are statistically significant (p < 0.001), with stronger warming rates over drier ecosystems. For example, the warming trend in NCDC is 0.39 ± 0.05 and 0.19 ± 0.03°C/10 yr for the driest and wettest regions, respectively, in the case of 7 ecoregions. The corresponding values are 0.39 ± 0.06 and 0.21 ± 0.03°C/10 yr for the 35 ecoregion case. The GISS and CRU data show very similar warming trends as those in NCDC. Figure 4 shows T trend as a function of EVI in terms of 7 and 35 ecoregions. Evidently, the warming rate depends strongly on ecoregions. T trend increases dramatically with decreasing EVI, particularly over very dry areas, indicating the lower the vegetation greenness, the stronger the warming trend. Five different regression lines (exponential, linear, logarithmic, polynomial and power) are fit between T trend and EVI but the negative logarithmic function works best, with R 2 = 96 and 89% for the 7 and 35 ecoregions, respectively (figures 4(a) and (b)). Other clarifications of ecoregions (table 1) agree that the negative logarithmic function best describes the T trend −EVI relationship by ecoregion. Furthermore, the CRU, GISS and NCDC datasets display very similar relationships individually (figures 4(c) and (d)). Note that the R 2 values in CRU are smaller than those in GISS and NCDC because CRU has missing or low-quality data over several data-scarcity regions as mentioned previously. To examine whether the negative logarithmic T trend −EVI relationship varies with observational datasets and depends on how the ecoregions are classified, we perform the same fitting by ecoregion for other ecoregion classifications (table 2) and for each dataset (supplementary table S3). The negative logarithmic fit performs better than the linear fit and this is consistent across all classifications and in all the datasets. The fitted coefficient (A 0 ) decreases slightly with the increasing number of ecoregions, and so does R 2 . When more ecoregions (or EVI bins) are considered, R 2 decreases because more small-scale factors affect the spatial variations of T trend . As the estimated trends may be sensitive to start or end points of data record, we also perform the same analyses as done above but consider other study periods (e.g., 1979-2010 and 1981-2012) by changing the start and/or end dates. Again, the negative logarithmic relationship remains robust.

Observational temperature trends
One way to validate the robustness of the above T trend −EVI relationship by ecoregion is to simply reconstruct the spatial patterns of T trend from the ecoregion map of EVI clusters (figures 2(a) and (b)) multiplied by the average temperature trend for each of these clusters (figures 4(a) and (b)). If the reconstructed T trend resembles (the first order) the observed trend ( figure 1(a)), this supports the idea that the warming rate depends mostly on the ecoregions and the classification of warming rates based on ecoregions is useful. Figure 5 shows the spatial patterns of reconstructed T trend , which exhibit warming patterns similar to the observed at large scales. The strongest warming rates are consistently seen over the driest ecoregions. The spatial correlation at the grid level is 0.55 (between figures 1(a) and 5(a)) and 0.58 (between figures 1(a) and 5(b)), respectively. These correlations are statistically significant at p < 0.001 (n = 1338). Figure 1(b) shows the spatial patterns of ALL T trend over land for the period 1979-2012. The warming is ubiquitous and the strongest warming occurs primarily in arid and semi-arid regions such as Northern Africa, Middle East, and Western US. All the land grids (50°S-50°N) exhibit a statistically significant linear trend (p < 0.05), which is much higher than the observations. The observed T trend −EVI relationship remains true also in ALL (figure 6). The R 2 value for the logarithmic (linear) fit is 99 and 92% (92 and 85%) for the 7 and 35 ecoregions, respectively. Other ecoregion clarifications (table 2) support consistently  that the negative logarithmic fit better describes the T trend −EVI relationship by large-scale ecoregion than the linear fit. Overall ALL largely reproduces the observed warming spatial patterns and the negative logarithmic T trend −EVI relationship by ecoregion. At the grid level, the spatial correlation is 0.75 (n = 1338) between the observed and ALL T trend , indicating that some regional biases in ALL simulations due to models' difficulties in simulating T change at small scales. Figure 1(c) displays the spatial patterns of NAT T trend over land for the period 1979-2005. The warming trend is strongest in several arid and semi-arid zones. For the 1537 land grids (50°S-50°N), only 27 and 13% exhibit a linear trend that is statistically significant at the 10 and 5% level, respectively.  The corresponding large-scale T trend −EVI relationship (figure 6) displays larger warming rates over drier (or lower EVI) ecoregions as shown in ALL, but has a much smaller warming magnitude than ALL. The spatial pattern of T trend (figure 1(c)) and the T trend −EVI relationships by ecoregion in NAT differ substantially from the observations and ALL (figure 6; table 2), indicating that anthropogenic forcings should play a key role in determining the observed T trend −EVI relationships.

Temperature trends in CMIP5 simulations
The histogram of T trend for the CMIP5 ALL and NAT simulations, together with that from the three observational T datasets, is shown in figure 7. In general, ALL T trend is within the observed T trend range but has a narrower distribution skewed to the right (warming) side, indicating higher warming rates and stronger spatial coherence in ALL than observed. This is expected as the multi-model ensemble mean in ALL represents primarily the forced signal (Dai 2013), and ALL slightly overestimates the warming trends for the period 1979-2012 because the simulated T anomalies are higher than observed during the global warming hiatus (IPCC 2013). In contrast, NAT shows small and equal possibilities of cooling and warming trends, which are mostly skewed to the left side (cooling and small warming) and outside the range of ALL simulations and observations, indicating again that the role of anthropogenic forcing in reproducing the observed warming patterns.

Discussion and conclusions
This paper explores spatial relationships between the observed land surface warming rates and the climatological EVI in terms of large-scale ecoregions for the last three decades. The regional mean temperature trends exhibit significant spatial dependence on the regional mean EVI. In general, the warming rate increases dramatically with decreasing EVI, with the strongest warming rate over the driest ecoregions. When anthropogenic and natural forcings are included, climate models are generally able to reproduce observed major features of the spatial dependence. When only natural forcings are used, none of the observed features are simulated. Furthermore, the temperature changes simulated in NAT are mostly far outside the range of those simulated in ALL and observed. These results indicate stronger warming amplification over drier ecoregions, pointing mainly to human influence.
We focus our analysis mostly on the negative logarithmic T trend −EVI fit because it generally has the highest R 2 than other fitting functions. However, all fit functions indicate consistently the strongest warming rates over the driest ecosystems for the period 1979-2012. The warming rate of T trend is determined by the magnitude of the surface radiative forcing and any feedbacks involved related to land surface and atmospheric boundary layer processes (Hasselmann 1976, Betts et al 1996, Zhou et al 2007, McNider et al 2012, Guo and Dirmeyer 2013, Davy and Esau 2014. Evidently, the finding of the strongest warming in very dry regions needs a new mechanism to explain it. Our further analyses of multiple observational variables, combined with reanalysis data and CMIP5 simulations, could attribute the spatial dependence of warming on ecosystem to much stronger water vapor and ecosystem feedbacks in response to the positive large-scale Table 2. The fitted coefficients and goodness of fit (R 2 ) of the linear and logarithmic fits between EVI and T trend (°C/10 yr) from observations and CMIP5 simulations by large-scale ecoregion.  GHGs enhanced longwave radiative forcing over drier areas, which will be detailed in future work. Our results suggest a fundamental pattern of global warming over land that depend on the dryness of ecosystems in mid-and low-latitudes. The ecoregiondependent warming may reflect primarily the largescale thermodynamic component of global warming linked to changes in the water and energy cycles over different ecosystems. It is very likely that the reconstructed T trend (figure 5) provides the first order estimate of this pattern and its differences from the observed T trend are related to the dynamical component linked to changes in mean circulation. So this analysis provides a new way to break down the total global warming patterns into the dynamical and thermodynamic components and thus help us to understand whether they reinforce or counteract each other over different ecosystems. For example, the reconstructed T trend from observations (figure 5) resembles the CMIP5 ALL T trend (figure 1(b)) because they share similar large-scale thermodynamic responses and feedbacks to global warming; the stronger warming trend over Central USA and Amazonia in the observations than in the reconstructions have been linked to increases in local drought conditions associated with sea surface temperature patterns (e.g. McCabe et al 2004, which in turn have been linked to aerosol changes (e.g. Cox et al 2008, Booth et al 2012. This warming amplification over very dry surfaces might help interpret and attribute global warming and assess climate impacts. The presence of the much stronger warming over arid regions in ALL than in NAT indicates that increasing GHGs may play a crucial role in the enhanced warming in the driest ecosystems. Hence our work points to new potential fingerprints of anthropogenic warming that could be used to further distinguish observed temperature trends from natural variations.