Is there a relationship between the magnitude and spatial extent of surface temperature anomalies and global warming? Why has climate change intensified in recent years?

Global mean surface temperature (GMST) is calculated by taking a spatial average of surface temperature across the globe. Traditionally reported in anomalies, GMST is widely used to describe the magnitude of climate change and often serves as the basis for defining policy targets, such as keeping anthropogenic warming below 1.5 or 2°C relative to preindustrial levels (Gulev et al., 2021). While GMST provides a useful benchmark for scientists and policy makers, as a global average, its meaning is far removed from the local climate effects experienced by individuals around the world (Hansen et al., 2012). Differences in the levels of warming experienced at different locations create disparities in the level of risk from various climate hazards (Arnell et al., 2019). Additionally, spatial heterogeneity in temperature variability drives differences in community vulnerability to climate change (Thornton et al., 2014) as well as differences in the physical manifestation of the warming trend—for instance, low intrinsic variability in the tropics means the frequency of hot extremes is at risk to increase faster (Simolo & Corti, 2022). In light of this, as anthropogenic climate change has continued to worsen into the 21st century, a significant amount of attention has been paid to understanding the evolving spatial effects of increasing global mean temperatures.

A number of metrics have been used to understand the spatial extent of modern warming, exploring concepts ranging from the distribution of extreme events (Diffenbaugh et al., 2017) to irreversibility in the response of surface temperature to CO2 forcing (Kim et al., 2022). One of the most straightforward ways to understand the spatial effects of climate change is to assess the percentage of global area exceeding a certain temperature threshold. Earlier uses of this framework (Coumou et al., 2013; Hansen et al., 2012; Rahmstorf & Coumou, 2011) focused on extremes in seasonal land temperatures, assessing the percentage of land area experiencing temperatures three or five standard deviations above a common reference period. Using the statistical framework proposed by Rahmstorf and Coumou (2011), further studies were able to describe the crossing of temperature thresholds in historical and future model simulations using models based upon a normal distribution .

Preparing for climate change requires an understanding of the degree to which global warming has regional implications. Here we document a strong relationship between the magnitude and extent of warming and explain its origin using a simple model based on binomial statistics. Applied to HadCRUT5 instrumental observations, the model shows that 96% of interannual variability in the proportion of regions experiencing anomalous warmth over the last century can be explained on the basis of the magnitude of global mean surface temperature (GMST) anomalies. The model performs similarly well when applied to a variety of unforced and forced model simulations and represents a general thermodynamic link between global and local warming on annual timescales. Our model predicts that, independent of the baseline that is chosen, 95% of the globe is expected to experience above‐average annual temperatures at 0.7°C of GMST warming, and 99% at 1.0°C of warming.

Whereas many studies focus on global mean surface temperature

(GMST), it is also critical to understand the spatial patterns of this warming in preparing for a warmer world. Using a simple statistical model, we describe a strong relationship between GMST and the spatial extent of above‐average temperature anomalies each year. This relationship is shown to hold across various model simulations and historical temperature data, capturing both annual and decadal trends. We can use our model alongside projections to gain insight into the extent of warming variations. For instance, 1°C GMST warming from a given baseline climate generally corresponds to a world in which 99% of places across the globe will be anomalously warm every year.

Robinson, 2013). Research using this and other methods documented how the increase of global land area affected by heat extremes has risen with global mean temperatures (Robinson et al., 2021; van der Wiel & Bintanja, 2021). The spatial structure and origins of this relationship, however, have not been explored intensively.

It might seem difficult to relate global and local climates given the heterogeneity of local climate patterns. For example, Arctic amplification is associated with more than twice the warming rate relative to the rest of the globe (Estrada et al., 2021; Rantanen et al., 2022), and land has warmed 50% faster than the ocean surface (Wallace & Joshi, 2018). Given this complexity, it is notable that a first‐order model based on normal statistics is able to predict the spatial extent of extreme temperatures from the global mean (Coumou & Robinson, 2013). One potential explanation could be the presence of an underlying relationship between climate variability and the strength of the local warming trend. Evidence for such a relationship can be seen in recent analysis of instrumental records by Zeng et al. (2021), where temperature trends normalized against internal variability show a higher level of consistency across latitudes compared to the raw trends. In this paper, we adapt the normal statistics of previous work (Coumou & Robinson, 2013; Rahmstorf & Coumou, 2011) and apply it at a global scale to understand the nature of the relationship between global mean temperatures and the spatial extent of above‐average warming.

To quantify the spatial extent of a given climate event, we define a metric that evaluates the proportion of the globe that is warm or cold at any given time. First, annual averages are taken for each temperature record to remove seasonal patterns and allow for easy intercomparison across hemispheres. Then, a local average temperature is defined for every individual gridbox using a common reference period. Each local timeseries is then dichotomized into a binary variable representing whether that area was warm or cold in a given year. If temperature is above its local average, the variable takes a value of 1, representing a warm year, and if it is below local average, the value becomes 0, signifying a cold year. From there, each dichotomized timeseries is averaged to find the percentage of the globe that is warm or cold. The average of dichotomized warming then indicates the fraction of the globe experiencing anomalous warmth in a given year. Although this definition of anomalous warmth is somewhat arbitrary, its simplicity makes it useful for exploring the relationship between magnitude and extent of global warming.

We adopt a simple model relating local and global temperature anomalies,

Tlocal − μlocal = a(Tglobal − μglobal + ϵ), (1)

where T represents local or global temperatures, μ represents local or global means, and ϵ represents local temperature fluctuations that are independent of GMST. The scaling factor a represents the strength of the local warming compared to the global signal. For example, local temperatures in the Arctic, which have warmed at higher rates than GMST since 1960 (Rantanen et al., 2022), would have higher values of a.

We make two assumptions to simplify our analysis. First, we assume that the distribution of temperature anomalies around the global mean follows a normal distribution. Second, we take the ratio of the local warming rate to the standard deviation of ϵ to be homogeneous across the globe, allowing us to prescribe a single parameterization of ϵ globally. These assumptions are supported by analyses presented in Section 4.

The relationship between local and global warming expressed in Equation 1 becomes even simpler after dichotomization, D. Given that D = 0 if T ≤ μ and D = 1 if T > μ, Equation 1 becomes,

where T′ represents temperature anomalies. Note that Equation 2 is independent of a.

The normality assumed for ϵ allows us to write the distribution of local temperature anomalies as,

. Scatterplots of global temperature anomaly versus spatial extent, along with least‐squares fits to Equation 4. The bottom right plot shows the fit for the HadCRUT5 instrumental data set, and the other three plots show fits for forced and unforced runs of three different CMIP6 models.

The distribution is centered about the GMST anomaly and has variance σ2 equal to that of ϵ. Taking advantage of the equality from Equation 2, the probability that D(T′local) = 1 is

derived from the cumulative distribution function of the normal distribution.

We apply Equation 4 to localities across the globe with, as noted, a fixed parameterization of σ. Each annual local record is represented as an independent Bernoulli trial with probability of success, p, given by Equation 4. Summing the results of these Bernoulli trials, we obtain a Binomial distribution with expected value E[B(n, p)] = np, where n is the number of independent Bernoulli trials. The fraction of records experiencing warmer than average temperatures is equal to the Binomial result divided by n, and so has expected value p. This relationship, alongside Equation 4, allows us to draw a direct link between a GMST anomaly and its spatial extent.

Note. Applying our dichotomization procedure to 20 different CMIP6 models and fitting the results with Equation 4, we find a range of values for σ with generally low uncertainty. The average value of σ was 0.480 for unforced models and 0.484 for forced models, with average r2 values of 0.874 and 0.959 respectively. Models are forced according to the historical period of 1850–2014, with anomalies calculated against a 1961–1990 baseline for easy comparison to the HadCRUT5 instrumental data set.

Our simple statistical model predicts a sigmoidal relationship between the global temperature anomaly and the fraction of the globe experiencing anomalous warmth, with one free parameter, σ. In order to explore how well this model captures the behavior of surface temperatures, we apply it to both instrumental and model data sets. First, applying our method to simulated surface temperature data from CMIP6 (Eyring et al., 2016), we

dichotomize the temperature record at each gridsquare and fit the spatial average to Equation 4. For each model, we fit data from a control run without any climate forcing, as well as from a historically forced model run. For comparability with the instrumental data set, we set the temperature baseline for forced model runs as 1961–1990, with the baseline for unforced runs set as the full timeseries average. Examples of the model fits are shown in Figure 1, with the full results in Table 1.

The models show a range of fitted values for σ. Quality of fit across forced simulations is excellent, with r2 values averaging 0.96. Given the lack of global climate forcing influencing the control simulations, as well as the smaller overall range of temperature anomalies, we expect a slightly lower quality of fit for our sigmoidal model in the unforced case. Although this pattern is observed, r2 values for the unforced simulations are still fairly high, averaging 0.87 across the 20 models. Moreover, the values of σ across forced and unforced models have remarkable agreement, with σ averaging

0.484 ± 0.022 in simulations that are forced with greenhouse gases and 0.480 ± 0.024 in control simulations. This regularity implies that the relationship between extent and GMST is similar for internal variations, such as ENSO (El Niño Southern Oscillation) variability, and external forcing from greenhouse gases, consistent with the observed behavior that spatial extent and GMST track each other in both high and low‐frequency variations

(Figure 2).

Next, in order to assess the model's performance on real‐world climate data, we apply the same procedure to the HadCRUT5 instrumental data set (Morice et al., 2021). The data set stretches from 1850 to the modern day, but we exclude the years 1850–1900 from our analysis because of their relative lack of spatial coverage and potential for greater uncertainty, including in interannual trends (Chan et al., 2023).

Fitting the GMST and fractional warm area calculated from the instrumental data to Equation 4, we obtain a value of σ equal to 0.424 ± 0.024 and an r2 value of 0.96. Despite the simplicity of our model, we are able to explain the vast majority of variations in the spatial extent of anomalies as a function of

GMST. Plotting the observed global fraction of anomalous warmth alongside

the predictions, Figure 2 shows that the spatial pattern of climate tracks GMST closely both with respect to underlying trends and higher‐frequency events associated with, for example, El Niño events. As seen in Figure 1, the quality of fit is similar between models and instrumental data. The fitted value of σ for the instrumental data set is slightly lower than the average value of the model simulations but falls within the 95% confidence interval of several of the forced model simulations.

There are several assumptions made in our analysis that warrant further scrutiny involving normality, constant signal‐to‐noise ratios, homogeneous Bernoulli trials, and resolution. There are also larger issues of the degree to which homogeneous thermodynamics can be assumed. We treat each of these issues in turn. First, we have assumed that ϵ is normal across the analysis. An examination of the residuals associated with Equation 1 at the 95% confidence level shows consistency with normality for 89% of local grid squares, such that the rate of rejection of the null is similar to that expected by chance alone.

Second, we have assumed that ϵ in Equation 1 can be described as statistically stationary across the globe. This is equivalent to assuming a constant signal‐to‐noise ratio where the signal is aT′global and the noise is aσ, where σ is the standard deviation of ϵ. It is possible to directly evaluate this assumption of homogeneity by considering a more typical least‐squares problem, whereby Tlocal − μlocal = a(Tglobal − μglobal) + η. Fitting this equation for each HadCRUT5 gridbox and plotting the estimate of a versus the standard deviation of η demonstrates covariance,

Panel (a) shows the observed spatial extent of warming (black) along with the spatial extent predicted by GMST using Equation 4 (brown), compared against the GMST warming relative to a 1961–1990 baseline (red). Data is from the non‐infilled HadCRUT5 data set (Morice et al., 2021). Each local 5° × 5° gridsquare is compared to its 1961–1990 average to determine whether a year is locally warm or cold. The maps in panel (b) show the spatial distribution of warm and cold areas in the years 1900 and 2020.

where the Pearson cross‐correlation is r = 0.62 (Figure 3). If our assumption that ϵ is statistically stationary holds, the slope between values of a and the standard deviation of η should be consistent with the value of σ from our global fit. The empirical slope of this relationship is 0.38, which is 90% the value of σ calculated with Equation 4, indicating consistency given expected levels of regression dilution.

Third, despite overall consistency with normality and an approximately constant signal to noise, there are a distribution of variances (i.e., the data in Figure 3 does not follow an exact line). This implies that, when summing the results of Bernoulli trials to get an estimate on P[D(T′local) = 1] , there will be deviations from a uniform probability of Bernoulli success p. The sum of non‐identical Bernoulli trials is known as the Poisson Binomial Distribution, and, although the classical Central Limit Theorem does not apply to it, a variant known as the Lyapunov Central Limit Theorem does apply for a sufficiently large number of trials (Tang & Tang, 2022). From synthetic experiments, we find that order 10 trials are needed for approximate convergence. Our analysis of the instrumental data averages 2,240 samples per year over the period 1900–2021 from a spatial field estimated to have order 100 degrees of freedom (Jones et al., 1999). As such, we can estimate the result of the sum of n nonidentical Bernoulli trials with a single average probability of success p.

The scale that we represent local conditions at within the observations is the 5° × 5° grid on which the HadCRUT data is provided. It can be asked whether temperature at a finer resolution would behave similarly. This question can partly be addressed by examining model simulations that are conducted at different resolutions. In our analysis across 20 different models, we find no significant covariance between model resolution and the σ value that is inferred for each model (Figure S1 in Supporting Information S1). This lack of sensitivity to finer resolution is expected given that annual surface temperature variability, as noted, contains only order 100 degrees of freedom (Jones et al., 1999). In addition to smaller spatial scales, this method can be applied to shorter temporal scales. On monthly timescales, although the assumption that ϵ is statistically stationary breaks down, it is

Relationship between the sensitivity of local climate to global warming and the variability of local climate. Each blue circle is derived from a timeseries from a different instrumental gridsquare. A linear leastsquares regression is run for each timeseries using the equation

Tlocal − μlocal = a(Tglobal − μglobal) + η. The a parameter gives the sensitivity, and the standard deviation of η is taken as the internal variability. The red dashed line represents the predicted slope from the global fit to Equation 4, with red shading corresponding to the uncertainty (two standard deviations) in σ.

. Projected temperatures from the IPCC's SSP2 (blue) and SSP5 (red) scenarios, along with the thresholds for 95%, 99%, and 99.9% of global area to be warmer than the 1961–1990 climatology. 95% uncertainty intervals are shown in gray. Temperature projections are taken from the CMIP6 models in Table 1, after subsampling to include only those models with transient climate response (TCR) within the “likely” range (Hausfather et al., 2022) and removing the 1961–1990 average. The shaded regions show the upper and lower bounds of the 13 model runs analyzed. Instrumental temperature data from HadCRUT5 are shown in black.

nevertheless possible to relate monthly GMST to fractional anomalous area using our binomial model and the HadCRUT5 data, where we obtain an r2 value of 0.93.

Although the simplicity of the methodology provides an intuitive way to look at spatial climate, there are cases in which this simplicity fails to capture the complex behavior of the climate system. Notably, there has been an observed patch of cooling sea surface temperatures in the North Atlantic, possibly associated with slowing of the Atlantic Meridional Overturning Circulation (Caesar et al., 2018) as well as strengthening ocean heat transport out of the region and cloud feedbacks (Keil et al., 2020). The anti‐correlation between GMST and local sea surface temperatures in this area are unaccounted for in our model. Similarly, another hole in the global warming trend persists in the Southern Ocean, poleward of the Antarctic Circumpolar Current, which has been attributed to upwelling of cold deep water (Armour et al., 2016) and circulation changes induced by glacial meltwater (Rye et al., 2020). Because of this divergence from the global trend, we do not expect internal variability to be correlated with warming sensitivity in either of the warming hole regions.

We can assess the effect these regions have on our model by masking out the regions of delayed warming in the Southern Ocean (below 50°S) and the North Atlantic (between 45–15°W and 45–60°N), and performing the same analysis on the HadCRUT5 instrumental data set. This reduced analysis yields a slightly higher r2 value of 0.97, consistent with temperature variability in these areas diverging from the statistics described by our model. The improvement in skill afforded by omitting the warming holes is most notable in recent decades, as the fraction of the globe above 1961–1990 baseline approaches unity. Equation 4 yields a result that systematically overestimates the proportion of the globe experiencing anomalous warmth over the last 15 years by, on average, 5.8%. When the warming holes are omitted from the calculation of extent, the proportional warm area is within 1.6% of the predictions based on GMST, where GMST is still computed using global coverage (Figure S2 in Supporting Information S1).

The effectiveness of our model at predicting historical variations in the extent of warming suggests it can be used to predict how the spatial extent of warming will change with continued anthropogenic warming. By inverting Equation 4, we find the global temperature anomaly at which a given proportion of global area is warm,

where f is the fraction of global area that is above average temperature. Threshold GMSTs above which various percentages of the globe are expected to be anomalously warm are shown in Figure 4 along with the SSP2 and SSP5 warming scenarios from the IPCC AR 6 (Gulev et al., 2021).

As can be seen in both Figures 4 and 2, the climate system has already reached a point where more than 95% of locations on Earth are warmer than their 1961–1990 average every year. Disregarding the warming holes in the North Atlantic and Southern Ocean, that proportion is likely to exceed 99% in the coming 5 years. Although we utilize a 1961–1990 baseline because it is the default for the HadCRUT5 data set and provides good sensitivity to the entire

References

Armour, K. C., Marshall, J., Scott, J. R., Donohoe, A., & Newsom, E. R. (2016). Southern Ocean warming delayed by circumpolar upwelling and equatorward transport. Nature Geoscience, 9(7), 549–554. https://doi.org/10.1038/ngeo2731

Arnell, N. W., Lowe, J. A., Bernie, D., Nicholls, R. J., Brown, S., Challinor, A. J., & Osborn, T. J. (2019). The global and regional impacts of climate change under representative concentration pathway forcings and shared socioeconomic pathway socioeconomic scenarios. Environmental Research Letters, 14(8), 084046. https://doi.org/10.1088/1748‐9326/ab35a6

Caesar, L., Rahmstorf, S., Robinson, A., Feulner, G., & Saba, V. (2018). Observed fingerprint of a weakening Atlantic Ocean overturning circulation. Nature, 556(7700), 191–196. https://doi.org/10.1038/s41586‐018‐0006‐5

Chan, D., Gebbie, G., & Huybers, P. (2023). Global and regional discrepancies between early‐twentieth‐century coastal air and sea surface temperature detected by a coupled energy‐balance analysis. Journal of Climate, 36(7), 2205–2220. https://doi.org/10.1175/JCLI‐D‐22

Coumou, D., & Robinson, A. (2013). Historic and future increase in the global land area affected by monthly heat extremes. Environmental Research Letters, 8(3), 034018. https://doi.org/10.1088/1748‐9326/8/3/034018

Coumou, D., Robinson, A., & Rahmstorf, S. (2013). Global increase in record‐breaking monthly‐mean temperatures. Climatic Change, 118(3–4), 771–782. https://doi.org/10.1007/s10584‐012‐0668‐1

Diffenbaugh, N. S., Singh, D., Mankin, J. S., Horton, D. E., Swain, D. L., Touma, D., et al. (2017). Quantifying the influence of global warming on unprecedented extreme climate events. Proceedings of the National Academy of Sciences of the United States of America, 114(19), 4881–4886.

https://doi.org/10.1073/pnas.1618082114

Estrada, F., Kim, D., & Perron, P. (2021). Spatial variations in the warming trend and the transition to more severe weather in midlatitudes. Scientific Reports, 11(1), 145. https://doi.org/10.1038/s41598‐020‐80701‐7

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., & Taylor, K. E. (2016). Overview of the coupled model intercomparison Project phase 6 (CMIP6) experimental design and organization. Geoscientific Model Development, 9(5), 1937–1958. https://doi. org/10.5194/gmd‐9‐1937‐2016

Gulev, S. K., Ahn, J., Aðalgeirsdóttir, G., Albani, S., & Allan, R. P. (2021). Changing state of the climate system. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Blair Trewin. https://doi.org/10.1017/9781009157896.004

Hansen, J., Sato, M., & Ruedy, R. (2012). Perception of climate change. Proceedings of the National Academy of Sciences of the United States of America, 109(37), 14726. https://doi.org/10.1073/pnas.1205276109/‐/DCSupplemental

Hausfather, Z., Marvel, K., Schmidt, G. A., Nielsen‐Gammon, J. W., & Zelinka, M. (2022). Climate simulations: Recognize the “hot model” problem setting the agenda in research. Nature, 26.

Jones, P. D., New, M., Parker, D. E., Martin, S., & Rigor, I. G. (1999). Surface air temperature and its changes over the past 150 years. Reviews of Geophysics, 37(2), 173–199. https://doi.org/10.1029/1999rg900002

Keil, P., Mauritsen, T., Jungclaus, J., Hedemann, C., Olonscheck, D., & Ghosh, R. (2020). Multiple drivers of the North Atlantic warming hole. Nature Climate Change, 10(7), 667–671. https://doi.org/10.1038/s41558‐020‐0819‐8

Kim, S. K., Shin, J., An, S. I., Kim, H. J., Im, N., Xie, S. P., et al. (2022). Widespread irreversible changes in surface temperature and precipitation in response to CO2 forcing. Nature Climate Change, 12(9), 834–840. https://doi.org/10.1038/s41558‐022‐01452‐z

Morice, C. P., Kennedy, J. J., Rayner, N. A., Winn, J. P., Hogan, E., Killick, R. E., et al. (2021). An updated assessment of near‐surface temperature change from 1850: The HadCRUT5 data set. Journal of Geophysical Research: Atmospheres, 126(3). https://doi.org/10.1029/ 2019JD032361

Rahmstorf, S., & Coumou, D. (2011). Increase of extreme events in a warming world. Proceedings of the National Academy of Sciences of the United States of America, 108(44), 17905–17909. https://doi.org/10.1073/pnas.1101766108

Rantanen, M., Karpechko, A. Y., Lipponen, A., Nordling, K., Hyvärinen, O., Ruosteenoja, K., et al. (2022). The Arctic has warmed nearly four times faster than the globe since 1979. Communications Earth and Environment, 3(1), 168. https://doi.org/10.1038/s43247‐022‐00498‐3

Robinson, A., Lehmann, J., Barriopedro, D., Rahmstorf, S., & Coumou, D. (2021). Increasing heat and rainfall extremes now far outside the historical climate. npj Climate and Atmospheric Science, 4(1), 45. https://doi.org/10.1038/s41612‐021‐00202‐w

Rye, C. D., Marshall, J., Kelley, M., Russell, G., Nazarenko, L. S., Kostov, Y., et al. (2020). Antarctic glacial melt as a driver of recent southern ocean climate trends. Geophysical Research Letters, 47(11). https://doi.org/10.1029/2019GL086892

Simolo, C., & Corti, S. (2022). Quantifying the role of variability in future intensification of heat extremes. Nature Communications, 13(1), 7930. https://doi.org/10.1038/s41467‐022‐35571‐0

Tang, W., & Tang, F. (2022). The Poisson binomial distribution—Old & new. Statistical Science, 38(1). https://doi.org/10.1214/22‐sts852

Thornton, P. K., Ericksen, P. J., Herrero, M., & Challinor, A. J. (2014). Climate variability and vulnerability to climate change: A review. Global

Change Biology, 20(11), 3313–3328. https://doi.org/10.1111/gcb.12581 van der Wiel, K., & Bintanja, R. (2021). Contribution of climatic changes in mean and variability to monthly temperature and precipitation extremes. Communications Earth and Environment, 2(1), 1. https://doi.org/10.1038/s43247‐020‐00077‐4

Wallace, C. J., & Joshi, M. (2018). Comparison of land‐ocean warming ratios in updated observed records and CMIP5 climate models. Environmental Research Letters, 13(11), 114011. https://doi.org/10.1088/1748‐9326/aae46f

Zeng, X., Eyre, J. E. R., Dixon, R. D., & Arevalo, J. (2021). Quantifying the occurrence of record hot years through normalized warming trends. Geophysical Research Letters, 48(10). https://doi.org/10.1029/2020GL091626