Research ArticleCLIMATOLOGY

Climate models predict increasing temperature variability in poor countries

See allHide authors and affiliations

Science Advances  02 May 2018:
Vol. 4, no. 5, eaar5809
DOI: 10.1126/sciadv.aar5809


Extreme events such as heat waves are among the most challenging aspects of climate change for societies. We show that climate models consistently project increases in temperature variability in tropical countries over the coming decades, with the Amazon as a particular hotspot of concern. During the season with maximum insolation, temperature variability increases by ~15% per degree of global warming in Amazonia and Southern Africa and by up to 10%°C−1 in the Sahel, India, and Southeast Asia. Mechanisms include drying soils and shifts in atmospheric structure. Outside the tropics, temperature variability is projected to decrease on average because of a reduced meridional temperature gradient and sea-ice loss. The countries that have contributed least to climate change, and are most vulnerable to extreme events, are projected to experience the strongest increase in variability. These changes would therefore amplify the inequality associated with the impacts of a changing climate.


Anticipating how anthropogenic climate change will affect natural and human systems requires an understanding not only of the changes in the mean climate but also whether and how climate variability will change. For instance, present temperature fluctuations in the tropics are strongly correlated with negative agricultural, economic, and political impacts (13). The increasing frequency and persistence of extreme weather events (4) like heat waves or droughts can have severe impacts on many biological systems showing nonlinear responses to environmental change (5). These effects involve changes in the physiological responses of individual species (6), the risk of species extinction (7), the functioning of global terrestrial ecosystems (8), and human food security (9). Furthermore, climate variability can also be an indicator of the sensitivity of climate to perturbations, and variability changes have been proposed to be a forewarning of upcoming tipping points (10, 11). It is therefore important to understand the patterns of variability change and the mechanisms behind them.

Previous studies analyzing variability changes predicted by climate models have focused on changes in variance (12), mostly in high and mid-latitudes (13) and the global mean temperature (14, 15). Here, we analyze changes in the standard deviation (SD) of monthly temperature anomalies on a local scale, with the aims of (i) identifying regional hotspots where models agree on the largest magnitude changes and (ii) investigating the physical mechanisms behind the identified changes. Because of the variability of climate on multiple time scales, observational time series are usually too short to detect a clear human influence on climate variability (1618). Here, instead, we analyze all available state-of-the art climate models from the Coupled Model Intercomparison Project 5 (CMIP5) (19), which were used to inform the last report of the Intergovernmental Panel on Climate Change. Our time series consist of the historical simulation starting in 1850, followed by the future scenario with the largest greenhouse gas emissions, RCP8.5 (20). The 37 models in total are used to calculate changes until the year 2100, and simulations with nine of the models have been extended until the year 2300. In the following, we examine monthly temperature anomalies, where each anomaly value represents the deviation from the 30-year running mean of the same month at a particular grid cell in a model (see Materials and Methods). We then analyze the changes in SD over time that are associated with global warming.


Unequal regional changes

We loosely define hotspots as regions of relatively large variability changes that also show a large model agreement on the sign of the change. Figure 1 shows the projected relative changes in SD of monthly temperature anomalies until the end of the 21st century, averaged over all 37 models (see fig. S1 for model agreement). Several hotspots of increases in SD are revealed over (sub)tropical land, namely, Amazonia, Southern Africa, Australia, the Sahel, the Arabian Peninsula, India, and Southeast Asia. In boreal summer, large increases occur over Europe, North America, and near the Arctic coast (Fig. 1A). In all other seasons, SD decreases over northern mid-latitudes (fig. S2). Several of the identified hotspots are consistent with previous modeling studies (2123).

Fig. 1 Relative changes of SD of monthly temperature anomalies until the end of the 21st century.

(A) Boreal summer [June, July, and August (JJA)], (B) austral summer [December, January, and February (DJF)], and (C) the whole year, averaged over 37 climate models. In hatched areas, at least 30 of the 37 models agree on the sign of change from 1850 to 2100 expressed in the Kendall τ value [for example, see fig. S1A for the model agreement in (C)].

Absolute changes in variability are greatest in mid- and high latitudes in boreal autumn and winter (fig. S2), but these are also regions with the largest background levels of temperature variability. In relative terms, the largest changes in SD are increases of up to 40% in Amazonia and Southern Africa in austral spring and summer. Using temperature anomalies from all months of the year, Amazonia stands out as a region with a particularly large relative SD increase (Fig. 1C) because increases occur in all seasons. Although the models agree on the sign of change, the magnitude differs substantially between the models, with some models showing a doubling of SD in the tropical hotspots, especially in central and northern Amazonia. The changes continue in the nine models extended beyond 2100 and scale approximately with global warming. In the model average, each degree of global warming leads to an ~15% increase in SD in Amazonia, Southern Africa, and the Arctic coast during the local summer season, and a 10% increase in the subtropical hotspots of the Northern Hemisphere (fig. S3).

Previous studies have isolated the role of selected mechanisms in individual models by suppressing certain sources of variability (2426). This is not feasible across all of the CMIP5 models; hence, we take a different approach and analyze the changes in the terms contributing to the surface energy balance in the model mean, and identify important regional differences. In high latitudes, decreases in temperature SD are particularly pronounced over regions with sea-ice loss because the exposed ocean introduces a much larger heat capacity. This leads to more persistent but damped anomalies (12, 15, 27), in other words, SD decreases. Reductions in SD in mid-latitudes are not directly related to Arctic sea-ice loss. The largest SD decrease in the Arctic occurs in autumn [September, October, and November (SON)], the season with the largest sea-ice loss until 2100, whereas the SD decrease in mid-latitudes is most pronounced during Northern Hemisphere winter (DJF), the season where Arctic amplification is largest (Fig. 1 and fig. S2). This amplification of mean warming in high northern latitudes leads to a reduced meridional temperature gradient (fig. S4). If one assumes that typical wind patterns do not change, then incoming air from the North or South is less different in temperature than it was when the temperature gradient weakens, leading to smaller local temperature anomalies. Calculating the monthly anomalies of advection (the local temperature tendencies due to horizontal wind) indeed reveals a decrease of the advection-induced temperature variability (fig. S4). We thereby confirm previous studies showing that the decreased temperature variability in mid-latitudes is associated with a decreased meridional temperature gradient due to the polar amplification of global warming (13, 22, 27, 28). The CMIP5 models do not support the alternative hypothesis that Arctic amplification causes increased meandering of the mid-latitude westerlies and, hence, greater variance and persistence of mid-latitude weather conditions (29).

The increases of temperature variability that we find are largest in the season of maximum insolation (fig. S2). In the following, we therefore focus on the summer season. In the Southern Hemisphere, the results are qualitatively similar in all seasons. In the Northern Hemisphere, the summer season is the only season with an increased SD of temperature anomalies, whereas SD decreases in all other seasons (for the reasons described above). The large increases in SD over tropical land are particularly striking and have not been explained in previous analyses.

Changes in the Southern Hemisphere

Although sea surface temperature (SST) anomalies are known to affect temperature over land (24), SST variability appears not to be responsible for the increased variability on land, because of the low model agreement over the oceans in contrast to the large agreement in the identified land hotspots (Fig. 1 and fig. S1). Furthermore, in the Southern Hemisphere summer, advection only plays a minor role (fig. S4C). We therefore focus on the local energy balance over land.

From a local perspective, fluctuations in the downwelling radiation (short-wave plus long-wave fluxes) at the surface, Rd, can be seen as a direct atmospheric driver of temperature changes, which determines how much energy is available at the surface. We do not consider changes in surface albedo here because its contribution to monthly variability is likely small. We recognize that local conditions can feed back to the atmosphere (for example, via changes in cloud cover) to affect Rd, but note that Rd is also determined by processes on larger scales. A second important parameter is the evaporative fraction (EF) (the part of the available energy at the surface used to evaporate water), here calculated as EF = LH/(Ru + Rd), with upwelling radiation Ru and latent heat flux (LH). Because an increase in sensible heat flux is associated with temperature increase, whereas increased LH is not, EF determines by how much a change in atmospheric forcing translates into temperature change. The smaller the EF becomes, the greater the effect of variations in available net radiation on temperature, because more radiation is converted into sensible instead of latent heat.

Previous studies have shown that changes in soil moisture affect the partitioning of surface fluxes mostly in intermediate moisture regimes (30). Where evapotranspiration is limited by moisture but still has a substantial contribution to the energy balance, soil moisture has the largest effect on LH or temperature (and hence their variability). In this context, it is striking that regions with substantial increases in temperature variability coincide with regions identified as hotspots of strong soil moisture temperature coupling in previous studies (30, 31). We therefore analyze how increases in temperature variability are related to changes in soil moisture and EF, in comparison to variations in downwelling radiation Rd. Figure 2A shows that the grid cells where temperature variability increases (red on the color scale) tend to have decreasing soil moisture (horizontal axis) as well as increasing variability in Rd (vertical axis). The same relation emerges when plotting EF instead of soil moisture (Fig. 2B). Therefore, the changes in the Southern Hemisphere hotspots are consistent with the explanation that soil moisture loss leads to smaller EF and thus larger temperature variability. Even in the humid climate of Amazonia where evapotranspiration is limited by radiation in many places, the large future drying brings the region closer to a soil moisture–limited regime, which is associated with a large reduction in the positive correlation between evapotranspiration and incoming radiation (fig. S5).

Fig. 2 Changes in variability of the surface energy balance until the end of the 21st century.

The horizontal axes in (A), (B), (D), and (E) show changes in the mean, whereas all other axes and the color bar show changes in SD. Each dot represents the difference between the periods 1875–1904 and 2055–2084 at a particular land grid cell in the Southern Hemisphere without Antarctica (A to C) and the Northern Hemisphere without Greenland (D to F) for the multimodel average. All figures show summer conditions (DJF in the Southern Hemisphere, JJA in the Northern Hemisphere). The blue coordinate system in (C) and (F) defines the index shown in Fig. 3, with negative values indicating larger increases in LH variability and positive values indicating larger radiative variability.

Such a transition from a wet to a moisture-limited regime is associated with an increase in the sensitivity of LH to soil moisture fluctuations (even in the absence of any changes in available radiation) (30, 32). The variability changes in LH and Rd can therefore potentially reveal information on the mechanisms behind changes in variability. In agreement with the soil moisture hypothesis, increased temperature variability is associated with substantial increases in the variability of LH (Fig. 2C). The SD of Rd and LH both increase to a similar extent, and thus, both contribute to increased temperature variability. However, as the changes in EF and soil moisture show, Rd cannot be the sole driver of the increased temperature variability because there is no clear reason why the increased SD of Rd would systematically deplete the mean soil moisture. Instead, we suggest that the SD increase in downwelling radiation in these areas is at least partly associated with land-atmosphere feedbacks. For example, surface drying can affect cloud distribution and, hence, the variability of Rd. To separate these effects quantitatively, specific additional model simulations would be needed.

Changes in Northern Hemisphere summer

The Northern Hemisphere shows a very different pattern (Fig. 2, D to F). There, soil moisture is not a good explanatory variable in general (Fig. 2D). This occurs because soil moisture does not limit EF in many places, so that decreases in soil moisture do not necessarily cause increased temperature variability. Moreover, changes in EF alone cannot explain the increased temperature variability in the Northern Hemisphere (Fig. 2E), and increases in Rd variability are needed as an explanatory variable at many places. The few exceptions with decreased temperature SD in the range of +4 to 6 W/m2 of downwelling radiation SD (Fig. 2E) occur at grid boxes in very high northern latitudes (mainly Northeast Canada) where the fluctuations in available energy are often compensated by phase changes (melting snow and water, or freezing water) instead of temperature fluctuations.

This pattern suggests that changes in the atmosphere on a regional scale are an important driver of the changes in variability at the surface. When we compare changes in the variability of Rd and LH in the Northern Hemisphere summer, we find two distinct regimes (Fig. 2F). The horizontal (bottom-right) branch with substantial increases in LH variability is consistent with the soil-drying mechanism explained above. Because Rd variability does not increase or even decreases, the soil moisture loss must be solely responsible for the enhanced temperature variability. In contrast, the vertical (upper) branch in Fig. 2F consists of dry regions where other mechanisms dominate, whereas changes in LH variability are small.

To illustrate the geographical distribution of these differences, we quantify the relation between SD changes in LH and net downwelling surface radiation (Rd) as an index (Fig. 3). This index (I) is proportional to the angle in the diagram spanned by the two properties (blue coordinate system in Fig. 2, C and F)Embedded Imagewhere ΔSD(Rd) and ΔSD(LH) stand for the changes in the SD of Rd and LH. The factor 2/π and the shift by 0.5 transform this angle in such a way to obtain a value of −1 if the increase in LH SD equals the decrease in Rd SD, 0 if both increase by the same amount (most points in Fig. 2C lie in this direction), and +1 if Rd SD increases as much as LH SD decreases. In Fig. 3, all clearly negative values (blue areas) are associated with variability increases in the latent heat flux, which cannot be explained with radiation alone. All clearly positive values (red areas) are associated with an increased variability of the downwelling radiation Rd, which cannot be explained with local hydrological changes. For example, heat waves in European summer (33) show up as dominated by soil moisture loss (the horizontal branch in the scatter cloud of Fig. 2F), whereas, in the subtropics of the Northern Hemisphere (mostly the Sahel, Arabia, and India), radiation often dominates over LH changes (vertical branch in Fig. 2F). The latter indicates that local changes in hydrology are too small to explain the radiative changes and that processes other than soil moisture loss are causing increased temperature fluctuations.

Fig. 3 Index related to the changes in the SD of LH versus downwelling net radiation.

See main text and Fig. 2 (C and F). Negative values indicate that changes are driven by soil drying, whereas positive values indicate a larger contribution of atmospheric variability. Ice sheets and locations where changes in temperature SD are smaller than 0.1K have been masked out (white areas). Each hemisphere shows summer conditions (JJA in the North, DJF in the South). All changes are calculated from the multimodel average and the period 1875–1904 versus 2055–2084.

Because model agreement tends to be lower in India and the Sahel than in other hotspots (Fig. 1 and fig. S1), we also calculate model averages for only the models showing an increase in temperature SD in a given region. In India, the SD of rainfall, soil moisture, cloud cover, and LH all increase despite little changes in mean soil moisture. This indicates that the Indian summer monsoon in the models is becoming more variable under greenhouse forcing, a result that is confirmed by other studies (34, 35). A similar effect occurs in Southeast Asia (36), in combination with a mean soil moisture decrease, causing an additional SD increase in temperature. In the Sahel, summer rainfall protrudes further north in those models showing a clear increase in temperature SD. Where the preindustrial climate has been almost permanently dry, an increased occurrence of cloudier and wetter months promotes increased variability in both radiation and LH. Moreover, because the Sahara warms more than the vegetated regions to the South, the meridional temperature gradient in the Sahel region is enhanced (fig. S4B). This enhancement of variability due to the increased spatial temperature contrast is hence the reversed mechanism of the decreased variability in mid-latitudes. An increased temperature gradient is also the reason for the zone of large SD increase surrounding the Arctic Ocean (fig. S4A), which warms much less than the land.


Model evaluation and time of emergence

The mechanisms that we identify above are in line with previous studies that have usually focused on temperature extremes (extreme compared to current climate) and on shorter time scales (37, 38). When comparing the occurrence of such extreme events between models and observation-based data sets, previous studies have usually found good agreement regarding their frequency and trend patterns (39, 40). Here, we briefly compare the modeled variability of monthly mean temperatures with four reanalysis data sets (see Materials and Methods for more detail): ERA-Interim (41), NOAA-CIRES 20th Century Reanalysis V2c (42), NCEP/DOE Reanalysis 2 (43), and MERRA2 (44). Although locally constrained by observations, these reanalysis data sets allow a global and temporally complete coverage. We calculate the temperature SD in the period 1980–2010 for all models and the four reanalyses by linear detrending for each month and computing the SD of the residuals. Overall, our comparison confirms that the current temperature SD of the CMIP5 models is within the observational range (fig. S6). Regional biases occur mostly in high latitudes over areas of substantial variability in sea-ice cover and over some tropical land regions (fig. S6A). However, these model mean biases are comparable in magnitude with the typical differences between models (fig. S6B) and, over some regions, are also comparable to the differences between reanalysis data sets (fig. S6C).

CMIP5 models are particularly uncertain and biased in the Amazon. Although most models show a pronounced hotspot of future temperature SD increase, this is also the region where the magnitude of this increase is most uncertain. Overall, the regions that the models identify as hotspots of future temperature SD increase tend to have a positive bias, that is, the model average shows larger present-day variability than reanalyses (fig. S6A). However, selecting groups of models based on their overall sign of biases hardly affects the general pattern of changing temperature variability. Although there is a tendency for models with positive bias to also show larger increases in temperature SD in many places, this bias is relatively small compared to the predicted changes in most hotspot areas, and the position of hotspots and the relative magnitude of temperature variability changes in these regions are not substantially affected.

Donat et al. (37) recently raised the question whether the predictions of the models are reliable or whether they share common biases in their predictions of temperature extremes, despite the overall agreement regarding the magnitude of observed variability. In the case that the model predictions are correct, the question arises whether the anthropogenic fingerprint is already apparent in observations of temperature variability. We therefore recalculate the changes in temperature SD for the period from 1958 to 2017 and compare the result to two observational data sets based on in situ measurements from meteorological stations: HadCRUT4 (45) and the GISS (NASA Goddard Institute for Space Studies) Surface Temperature Analysis (GISTEMP) (see Materials and Methods) (46). These data sets report temperatures whose mean and variance have been standardized with a reference period. This procedure is known to introduce biases in variance changes (47), although the effect is likely small here because the trends that we find in the observational data sets tend to also show up in the reanalyses that have no such bias (17).

At first sight, the observations do not agree well with the signal in the model predictions (Fig. 4). The most apparent features in observations are an increase in temperature SD in the Mediterranean region and a decrease in some parts of North America and Northern Europe (Fig. 4, A and B). In contrast, the tendencies of changes in temperature SD in the models already show the overall pattern that we analyzed above when all models are considered (Fig. 4C). However, the model agreement is still relatively low, and the characteristic contrast between high and low latitudes is not apparent in individual models with very few exceptions. If we assume that each model is one realization of climate variability in the real world, then this means that we cannot expect to see the anthropogenic signal in recent observations. It is only during the 21st century that the signal emerges in the large majority of the models, with more than 80% of the models agreeing on the sign until 2100 over subtropical land and the mid- and high latitudes (fig. S1). Although the sign of variability change in observations can be very clear at places (indicated by a large absolute Kendall τ value), these observations still provide only one realization.

Fig. 4 Trends in temperature variability in the period 1958–2017 in observations and models.

(A) Observations from HadCRUT4, (B) observations from GISTEMP, and (C) CMIP5 models. (A) and (B) show Kendall τ values, whereas (C) shows the number of models agreeing on the sign of the Kendall τ value. The window size for removing the annual cycle was 15 years, and the window to compute the Kendall τ value for the change in SD was 20 years.

In a recent companion study (17), we have shown that the sign of the observed temperature trends substantially depends on the choice of the time period, which confirms that the observed changes are largely related to internal climate variability and not a forced trend. The relatively robust decrease of observed variability over high northern latitudes in the North Atlantic region, however, might already be part of the anthropogenic signal (associated with sea-ice loss). Unfortunately, the regions with the largest model agreement in Fig. 4C (tropical land areas and high latitudes)—where the anthropogenic signal can be expected to emerge the earliest—are also among the regions with the fewest observations. It is therefore hard to conclude whether current observations confirm or refute the predictions of models regarding trends in monthly temperature variability.

A matter of climate justice

The predicted increases in temperature variability in the identified hotspots suggest the potential for substantial impacts on top of the impacts caused by changes in the mean climate. Poor countries are located in the identified hotspot regions of increasing temperature variability. A growing body of climate-economy literature shows that poor countries are particularly vulnerable to weather and climate shocks, whereas rich countries tend to show no significant vulnerability (2, 3, 48). Tropical temperature fluctuations are not only strongly correlated with a temporal loss of agricultural output, as might be expected from anomalous weather affecting crop productivity and yields (3, 9). They also affect the economic growth rate of countries, leading to long-term economic disadvantages compared to the development of other nations (3). Moreover, it has been argued that high temperatures favor political instability and conflict in tropical countries (3, 49), although this claim remains controversial (50). Like our own analysis, these studies are often based on monthly climate anomalies over many years and, thus, integrate many individual weather and climate events. The unequal projected regional changes in climate variability therefore constitute a matter of climate justice.

To illustrate this, we calculate the relative changes in temperature SD for individual countries, omitting small countries that are usually not resolved in the models (see Materials and Methods). It is particularly striking that the geographical pattern of the changes in SD coincides closely with the current distribution of wealth on the national level (Fig. 5). Although temperature fluctuations are projected to become smaller in countries having a large per capita gross domestic product (GDP), they tend to increase in countries with low GDP. Hence, countries facing the largest increases in temperature variability also have the least economic potential to cope with the impacts. The countries affected by this dual challenge of poverty and increasing variability already share half of the world’s population (Fig. 5), and population growth rates are particularly large in these countries (see Materials and Methods for data sources).

Fig. 5 Relative change in SD of monthly temperature anomalies until the end of the 21st century versus per capita GDP and greenhouse gas emissions.

(A) Per capita GDP and (B) per capita greenhouse gas emissions (without land use and forestry) between 1990 and 2013 in different countries. The red line marks zero change in temperature variability. The blue lines mark half of the world population (A) and half of the total worldwide emissions (B).

The unequal changes in projected climate variability are also a challenge for climate justice regarding the contributions of different nations to global climate change. Considering the average per capita emissions in the period 1990–2015 (excluding land use change and forestry), countries with the smallest emissions tend to have the largest increases in variability (fig. S7A), whereas those with the largest emissions tend to have smaller increases or even decreases in variability (Australia being a notable exception). The picture changes somewhat when including recent emissions from land use change and forestry (fig. S7B), bearing in mind that land use–related emissions are much more uncertain than emission estimates from fossil fuel burning. In the selected time period, tropical countries contributed a large fraction to land use–related emissions, associated with tropical deforestation especially in sub-Saharan Africa, Amazonia, and Southeast Asia.

On the basis of the econometric studies mentioned above, it is hard to conclude what the exact impacts of changing climate variability will be. For example, it is important to assess to what extent the impacts caused by warm climate anomalies are compensated by cold anomalies. In general, an increase in temperature variability is arguably more difficult to cope with than a decrease in variability. Societies and ecosystems are typically adapted to a certain long-term climate at any particular location—the more an anomaly deviates from the local mean, the more damage it can thus cause. This is particularly important in the presence of critical thresholds, for example, in species performance (6). Several studies find nonlinearity in the effect of temperature anomalies, for example, on agriculture, mortality, and energy demand (2), indicating that the effects of warm and cold anomalies will not simply cancel out. Moreover, the distribution of daily temperature is projected to become more skewed to warm anomalies in many of the identified hotspot locations (37). It therefore seems plausible that, overall, the projected increases in temperature variability would have negative consequences for nature and society.


Our analysis shows marked predicted increases in temperature variability in tropical land regions including many of the world’s poorer countries, with the Amazon being a particular hotspot of concern. These hotspots result mainly from soil drying in the Southern Hemisphere and from increased atmospheric variability in the subtropics of the Northern Hemisphere. Temperature variability is projected to decrease in mid-latitudes because of a decreasing meridional temperature gradient and in high latitudes because of sea-ice loss. Overall, the pattern of predicted changes agrees with previous studies (2123, 51, 52) but differs from observed changes in temperature variability. Models and observations indicate that the anthropogenic signal may not yet have emerged from the “noise” of natural variability but can be expected to do so during the 21st century. The identified hotspots are robust across models and mechanistically plausible. The predicted increase in the amplitude of temperature variability in much of the developing world, and the corresponding decrease in temperature variability in much of the developed world, could have substantial social, economic, and ecological consequences that add to the inequality in the impacts of a changing climate.


Models and scenarios

We analyzed all 37 models for which surface air temperature was available from the historical and RCP8.5 simulations: ACCESS1-0, ACCESS1-3, bcc-csm1-1, bcc-csm1-1-m, BNU-ESM, CanESM2, CCSM4, CESM1-BGC, CESM1-CAM5, CMCC-CESM, CMCC-CM, CMCC-CMS, CNRM-CM5, CSIRO-Mk3-6-0, FGOALS-g2, GFDL-CM3, GFDL-ESM2G, GFDL-ESM2M, FIO-ESM, EC-EARTH, GISS-E2-H, GISS-E2-R, HadGEM2-AO, HadGEM2-CC, HadGEM2-ES, inmcm4, IPSL-CM5A-LR, IPSL-CM5A-MR, IPSL-CM5B-LR, MIROC5, MIROC-ESM, MIROC-ESM-CHEM, MPI-ESM-LR, MPI-ESM-MR, MRI-CGCM3, NorESM1-M, and NorESM1-ME. The historical simulation starts in the year 1850, with the exception of the three versions of HadGEM2, where it starts in 1860. The RCP8.5 simulation ends in 2100 and has been extended to 2300 in the models bcc-csm1-1, CCSM4, CNRM-CM5, CSIRO-Mk3-6-0, GISS-E2-G, GISS-E2-H, IPSL-CM5A-LR, and MPI-ESM-LR, and to 2299 in HadGEM2-ES. For GFDL-ESM2M, years 2101 to 2200 were also available but were not used for our analysis.

In the RCP8.5 scenario, atmospheric CO2 shows an accelerated increase until the year 2100, when a radiative forcing of approximately 8.5 Wm−2 is reached. Thereafter, the CO2 concentration stabilizes at almost 2000 parts per million (20), yielding the largest warming of all CMIP5 simulations. For each model, we combined the historical simulation, the RCP8.5 simulation, and (if available) the extended RCP8.5 simulation, obtaining one long time series per grid cell and model.

Although temperature output was available from all selected models, this is not true for all variables that we analyzed. Specifically, LH was not available for EC-EARTH, downwelling net radiation was not available for FIO-ESM and EC-EARTH, and surface soil moisture (shown in Fig. 2, A and D) was not available for the CMCC models, FIO-ESM, EC-EARTH, GISS-E2-R, HadGEM2-AO, HadGEM2-ES, inmcm4, the IPSL models, MIROC-ESM, and the MPI-ESM models.

Observation and reanalysis data sets

We used two observational temperature data sets, based on in situ measurements from meteorological stations: HadCRUT4 (45) and the GISTEMP (46). The latter data set has been interpolated over a scale of 1200 km to fill spatial gaps between stations, which is the reason for the very uniform changes in the Arctic in fig. S8 (46). Because these data sets report temperature as standardized anomalies, we did not calculate absolute changes in variability, but only the direction of change as Kendall τ values. The GISTEMP data set has been obtained as version 5 from the GISS, The HadCRUT data set (version has been obtained from Both data sets have been downloaded on 31 January 2018. The four reanalysis data sets and their sources are the following:

(1) NOAA-CIRES (National Oceanic and Atmospheric Administration–Cooperative Institute for Research in Environmental Sciences) 20th Century Reanalysis V2c (42),

(2) NCEP/DOE (National Centers for Environmental Prediction/Department of Energy) 2 Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, CO, USA (43),

(3) MERRA2 (Modern-Era Retrospective Analysis for Research and Applications, version 2) (44),

(4) ERA (European reanalysis)–Interim (41),

All reanalysis data sets have been downloaded on 19 January 2018.

Calculation of anomalies and trends in SD

We used the monthly mean surface air temperature and removed the annual cycle at each grid cell within a running window of 30 years length, centered at the time point under consideration (fig. S8A). Choosing a baseline in such a local time window and considering each month separately are crucial to our analysis because it circumvents two problems: (i) the typical annual cycle changes over time (53, 54), for example, because the minimum and maximum sea-ice area occurs later in the year in a warmer world (54), and (ii) the choice of any specific reference period would lead to biases due to the limited sample size (47). Within the window and for each month, we subtracted the mean value of all 30 temperature values to obtain a temperature anomaly. Because of the necessity of 15 years of data to both sides of a data point, the procedure shortens each time series by 30 years.

Thereafter, we calculated the SD within a sliding window of 100 years that we moved along each anomaly time series at each grid cell in each model (fig. S8B). Hence, the window covers almost half the length of the whole time series like in previous studies (55). We calculated the trend in the obtained SD values based on the Kendal τ rank correlation coefficient (fig. S8C). A positive Kendal τ signals increasing SD, whereas a negative τ indicates decreasing SD. For example, a monotonous increase over time would yield τ = 1, whereas a monotonous decrease would yield τ = −1, and no trend would yield τ = 0. The Kendall τ values were used to assess the model agreement in fig. S1 and the hatching in Fig. 1. These sign agreements hence reflect the changes over the whole time series from 1850 to 2100.

For comparing trends in models and observations in the period 1958–2017, we used the same procedure but with a sliding window size of 15 years to remove the annual cycle, and 20 years to calculate the SD of the residuals. We took 1958 as the start year because we expected the anthropogenic fingerprint to be small before the mid–20th century and because only few observations are available from earlier years.

For comparing models and reanalysis data sets in the time period 1980–2010, we used a slightly different procedure: Instead of using a running window, we removed a linear trend from each month’s yearly time series. This linear trend removal is adequate because of the short time period and has the advantage that no data points are lost in the beginning and end of the time series. As we removed trends for each month separately, we still took into account that the typical annual cycle can change within the time period.

Time periods and multimodel averages

To calculate geographic maps (Fig. 1 and figs. S2 and S3) and scatter plots (Fig. 2) with changes in variability, we selected two time periods (1875–1904 and 2055–2084 to assess changes until 2100) from the anomaly time series mentioned above in each model.

When calculating seasonal changes, we only selected the anomalies from the specific months belonging to a season (hence removing all other months from the time series). For each model, season, and time period, we then calculated the SD of the anomaly time series for the variable under consideration and took the absolute and relative difference between the two periods for each season and model. Because the historical simulation from 1850–2005 includes the radiative forcing from volcanic eruptions (in particular the strong eruption of Mt. Krakatoa in 1883), whereas the future scenario does not, we verified the robustness of our results. For example, when using a base period from 1888 to 1917 (that is, excluding the Krakatoa eruption) or the preindustrial control simulations (which do not contain volcanic eruptions), our results remain essentially the same.

Thereafter, we averaged the resulting maps over all models. To do so, we brought the data on the same grid by interpolating each map to the grid of the model with the highest resolution (CCSM4), using a bilinear interpolation (CMCC-CM has even higher resolution but not all variables that we analyzed were available in that model, see above). The same method and grid were also used for interpolating the four reanalysis data sets. All multimodel averages are based on the same time period in each model.

Economic and geographic data

The GDP per person (based on purchasing power parity) for Fig. 5 has been obtained from for the year 2015. Countries smaller than 60,000 km2 were omitted. For four countries, data from the year 2015 were not available, in which case we chose the most recent year for which data were available: Eritrea (2011), Libya (2011), Papua New Guinea (2014), and Venezuela (2013). To calculate changes in SD in individual countries (Fig. 5), we interpolated the multimodel average of the relative annual SD changes (Fig. 1C) to a global grid of 0.1 × 0.1 degrees using a bilinear interpolation. We then picked all grid cells in a particular country and calculated the area average. The country boundaries for creating the selection mask have been obtained from Population estimates and projections have been obtained from Estimates for the year 2015 have been used for Fig. 5. Greenhouse gas emissions have been obtained from the World Resources Institute ( For Sudan, emissions from 2011 were used, whereas South Sudan was omitted. Country acronyms in Fig. 5 and fig. S7 are International Organization for Standardization alpha-3 codes.


Supplementary material for this article is available at

fig. S1. Model agreement on changes in SD in monthly anomalies for all months of the year.

fig. S2. Absolute and relative changes in SD of monthly temperature anomalies until the end of the 21st century.

fig. S3. Relative change in SD of monthly temperature anomalies per global mean warming.

fig. S4. Changes in temperature gradients and SD due to advection.

fig. S5. Temporal correlation between monthly anomalies in latent heat flux and downwelling net radiation.

fig. S6. Model evaluation of temperature SD between 1980 and 2010.

fig. S7. Greenhouse gas emissions per year and person between 1990 and 2013 versus relative change in temperature SD in different countries (using all seasons).

fig. S8. Documentation of our time series analysis method.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We are grateful to two anonymous reviewers for constructive comments; M. Flach, E. van Nes, and B. van der Bolt for valuable feedback; and P. Cox for his advice and inspiring discussions about the surface energy balance, climate variability, and the smell of rain. Funding: S.B. and M.S. acknowledge support by the Netherlands Earth System Science Centre, financed by the Ministry of Education, Culture, and Science. T.M.L. was supported by a Royal Society Wolfson Research Merit Award, the Natural Environment Research Council (NE/P007880/1), and the European Union Seventh Framework Programme FP7/2007–2013 under grant agreement no. 603864 [High-End cLimate Impacts and eXtremes (HELIX)]. V.D. was supported by an Adaptation to Changing Environments ETH (Eidgenössische Technische Hochschule Zürich) fellowship. Author contributions: All authors designed the research. S.B. analyzed the data and performed the analysis with input from V.D., M.S., and T.M.L. S.B. wrote the article with contributions from T.M.L. All authors discussed and proofread the article. Competing interests: The authors declare that they have no conflicts of interest. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups for producing and making available their model output. Observational and reanalysis data sets are freely available online. Support for the NOAA-CIRES Twentieth Century Reanalysis Project version 2c data set is provided by the U.S. DOE, Office of Science Biological and Environmental Research and by the NOAA Climate Program Office.

Stay Connected to Science Advances

Navigate This Article