Anthropogenic influence in observed regional warming trends and the implied social time of emergence

The attribution of climate change allows for the evaluation of the contribution of human drivers to observed warming. At the global and hemispheric scales, many physical and observation-based methods have shown a dominant anthropogenic signal, in contrast, regional attribution of climate change relies on physically based numerical climate models. Here we show, using state-of-the-art statistical tests, the existence of a common nonlinear trend in observed regional air surface temperatures largely imparted by anthropogenic forcing. All regions, continents and countries considered have experienced warming during the past century due to increasing anthropogenic radiative forcing. The results show that we now experience mean temperatures that would have been considered extreme values during the mid-20th century. The adaptation window has been getting shorter and is projected to markedly decrease in the next few decades. Our findings provide independent empirical evidence about the anthropogenic influence on the observed warming trend in different regions of the world. Most regions are experiencing mean temperatures which would have been considered extreme in 1880, while what is currently considered extreme may become normal within 10–20 years, according to statistical analyses of regional observational datasets.

T he attribution of the observed warming to natural and anthropogenic drivers is fundamental for the study of climate change [1][2][3][4][5] . Estimates of the contribution of natural and anthropogenic forcing to the observed changes in regional temperatures provide relevant information for policymakers to propose portfolios of mitigation and adaptation strategies 6,7 . These can also create more accurate perceptions of the risks of climate change and support pertinent climate policies. The Paris Agreement illustrates the international consensus about the anthropogenic contribution to warming and the need to curb greenhouse gas emissions 6 . However, in some political circles, (e.g., in the US), there seems to be an increasing skepticism that could hamper the success of international agreements 8 . This is important since even short delays in the participation of key countries in the Paris Agreement could render the climate targets unattainable 9 . Such a situation would require the adoption of prompter and stronger adaptation and risk reduction strategies. The current global health crisis 10,11 and the expected economic downturn could significantly reduce international efforts to mitigate climate change.
A clear human fingerprint on global and hemispheric warming has been corroborated by a variety of physical and statistical methods, sometimes based on conflicting assumptions and approaches [12][13][14][15][16][17] . In contrast, regional attribution studies of longterm trends have been largely based on climate models 1,18,19 , which depend on the models' performance to reproduce the observed climate and, despite their improvements, are less reliable at this scale 5,20 . Furthermore, model performance evaluation, and thus attribution results, may imply circular reasoning 21 .
Observation-based attribution studies do not depend on climate models' performance, tuning or calibration and can provide independent evidence (see "Methods").
The attribution of regional changes in climate can influence the perception and participation in risk management efforts due to a sense of information sufficiency and responsibility [22][23][24] . The time of emergence (ToE) 25 is a property of the climate system and provides an estimate of the date when the projected anthropogenic climate signal would exceed some measure of long-term natural variability. The ToE has also been used to estimate risk measures in ecological systems 26 as well as the benefits of mitigation policy for reducing climate risk 27 . It can further be extended to represent a risk in the socio-environmental realm as a property of the interaction of the climate and human systems. For this, a more relevant metric for decision-makers consists in determining when the climate signal becomes larger than the variability societies have experienced in the recent past and what they are presently prepared for. Current thermal comfort standards, infrastructure, services, productive, and recreational activities, energy consumption, disease distribution, among many others are more a function of observed rather than of long-term climate natural variability [28][29][30][31] . We define this metric as the social time of emergence (SToE) and it consists of the signal to noise ratio of the climate response to changes in radiative forcing relative to a measure of the observed natural variability to which a society is able to adapt ("Methods"). A related metric is defined as the Time to Adapt (TtA) which expresses the expected number of years to reach a novel climate.
Based on recent statistical methods, we analyze annual mean surface air temperature from a variety of sources to evaluate the anthropogenic contribution to the observed warming trends. The analysis is performed for three different spatial scales: latitude belts, continent, and country levels. A set of 13 countries was selected based on a combination of spatial extension, their aggregate emissions and fossil fuels reserves, and relevance for international climate policy negotiations. Detection of climate change based on analyzing observed trends is not addressed here, instead, we focus on inferring the underlying warming trends produced by global aggregates of external forcing 12,32,33 . Natural variability and external forcing factors with marked influence on regional scales modulate the underlying warming trends imparted by the globally aggregated radiative forcing 34,35 and can make observed and inferred warming trends differ in magnitudes and features such as breaks and their associated break dates 32,33 . We apply statistical methods to test the existence of common trends in nonstationary variables [36][37][38] . The tests allow us to evaluate the influence of anthropogenic forcing on the observed warming and whether the trend imparted by anthropogenic greenhouse gas forcing is present, accounting for factors such as the regional effects of forcings and natural variability that may have altered some relevant features (e.g., breaks in trend) of temperature series 12 .

Results
Common trends, the observed contributions of external forcing, and the evaluation of SToE and TtA. Cotrending tests [36][37][38] were applied to sets of variables that include the radiative forcing of well-mixed greenhouse gases (WMGHG), the sum of all natural and anthropogenic radiative forcing (TRF), and the regional temperatures of interest (see "Methods"). This allows us to investigate the attribution of the observed warming in regional annual temperatures to anthropogenic influence. The selection of variables allows us to test for a common trend in WMGHG, TRF, and regional temperature. Finding a common trend implies that WMGHG, which has an anthropogenic origin, imparts the secular movement in TRF and regional temperatures, allowing other factors such as low-frequency natural variability and the regional effects of external forcing to have an important modulating influence.
For every region, continent and country analyzed there is strong evidence of a common trend between forcing variables and annual mean temperatures due to anthropogenic forcing (Supplementary Tables S1-S18). Across tests and datasets, the existence of cotrending is seldomly rejected and in all cases at least two out of the three tests support the conclusion of a common trend (Supplementary Table S19). These results strongly suggest a warming trend attributed to anthropogenic forcing at all spatial scales: latitude belts, continent, and country levels. The warming has been widespread and followed a common secular movement, modulated by the effects of regional climate variability, and the regional effects of forcings. These conclusions are robust to the use of different radiative forcing, temperature datasets, and cotrending tests (Supplementary Tables S1-S19).
The existence of a common long-term trend allows the use of ordinary least squares regression to estimate the transient climate response (TCR) for each temperature series, which measures the response of regional temperatures to changes in external radiative forcing (Supplementary Tables S20-S22). This is in accordance with the established linear relationship between radiative forcing and perturbative temperature changes, including regional temperatures with the exception of areas in which strong regional effects of external forcing occur 1,39,40 . At the most aggregated level, only the TCR of the extratropical region in the northern hemisphere (24°N -90°N) is statistically larger (±2SE) than those of the southern hemisphere and the tropics. The results for finer latitude belts reveal that, in general, only the values for the high latitudes in the northern hemisphere are statistically larger than those for other latitude belts ( Supplementary Fig. S1). This is consistent with the results from climate models, in which the high latitudes warm faster than the rest of the world due to the arctic amplification phenomenon 41,42 . At the continental and country levels, the TCR tends to be larger at higher latitudes ( Supplementary Figs. S2 and S3). With the BEST dataset, Asia and North America have the largest TCR values (0.79 and 0.76°C per W m −2 , respectively) followed by Europe, Africa, and South America (0.70, 0.63, 0.62°C per W m −2 , respectively). The lowest value is for Oceania (0.52°C per W m −2 ), probably because of the ocean influence. These estimates are not statistically different from each other except for Asia and Oceania: for Asia, the TCR value is larger than those of Africa, South America, and Oceania; for Oceania, it is smaller than those of North America and Asia. Canada and Russia show the largest TCR of all countries, regardless of the temperature dataset used.
Estimates of the contributions of total and specific radiative forcing components to the secular trend at the country level can help decision-makers and the society to have a clearer picture of the magnitude of climate change and its causes. Using the estimated values of TCR, these contributions to the warming trends can be approximated (Fig. 1). To estimate the increase in temperature due to changes in both natural and anthropogenic radiative forcing, the TCR values are multiplied by TRF, while the TCR values are multiplied by WMGHG and solar forcing (SOLAR) to approximate the changes in temperatures due to these specific components of anthropogenic and natural forcings, respectively. Note that the estimates based on WMGHG and SOLAR do not represent the full response to changes in anthropogenic and natural forcings. The increase in anthropogenic radiative forcing due to WMGHG with respect to 1880 has produced significant warming in all selected countries, while non-WMGHG forcing had a modulating effect. By the year 2011, TRF contributed to an average increase of about 1°C, with warming near 1.2-1.3°C in France, Germany, and China and up to 1.6°C and 1.7°C in Canada and Russia ( Fig. 1a and Supplementary Fig.  S4a). Without the modulating effect of other forcing factors ( Fig. 1b and Supplementary Fig. S4b), WMGHG would have led to a warming of at least 1.5°C in most countries and more than 2.8°C in Canada and Russia. Solar forcing alone would have produced slight warming of about 0.1-0.2°C ( Fig. 1c and Supplementary Fig. S4c). However, these warming estimates do not provide a clear perspective about how extreme these changes are. The SToE and TtA metrics offer a simple and easy way to quantify the relative magnitude of the changes and how fast they are occurring.
The SToE was calculated for four reference years: 1880 and 1960, using the historical TRF 43 ; 2010 and 2020 using the projections of TRF from the Representative Concentration Pathways (RCP) RCP8.5, RCP4.5, and RCP2.6 scenarios 44 . For each case, the warming signal was obtained as ΔT ref ;i;t ¼ TCR i Ã ðTRF t À TRF ref Þ, where i and t denote region and time, respectively, and ref is the reference year ("Methods"). The estimates of the temperature variability experienced by societies (σ i ) are the standard deviation of the residuals when regressing temperatures on TRF. The social time of emergence is defined as is the indicator function and year is the calendar date ("Methods"). Note that for trend stationary variables such as temperatures 12,45 , a warming signal (ΔT) equal to Bσ i implies that the whole distribution of temperatures has shifted toward higher temperature values, not only the mean value. For example, if changes in temperatures are normally distributed and B = 2, the reference climate as experienced by the affected population at date SToE2 would have a very small overlap. The new mean temperature societies would experience at date SToE2 corresponds to the 97.5th percentile of the previous temperature distribution. For a SToE based on three standard deviations (SToE3), there would be no discernable overlap between the reference climate and the SToE.
The results show that most groups (whether by countries, continents, and latitude belts) are already experiencing annual mean temperatures that would have been considered extreme values with respect to the temperature distributions at the end of the 19th century or even during the second part of the 20th Fig. 1 Contributions to the warming trend from total and specific radiative forcing factors in the year 2011 with respect to 1880. Panel a shows the response of regional temperature trends to changes in TRF. The response of regional temperature trends to changes in WMGHG is shown in panel b. Panel c shows the response of temperature trends to changes in SOLAR. Estimates of radiative forcing for 2011 are from Hansen's dataset (see "Methods"). Units are°C.
century. For the reference year 1880, the SToE2 dates are similar for all groups (Tables 1-3 and Supplementary Tables S23-S25), with estimated values in the late 1970s. When considering 1960 as the reference year, the SToE2 values increase slightly to the late 1980s. Such small variations when changing the reference year indicate that the climate is warming at an accelerated pace. The general spatial patterns of SToE2 values are similar to those of ToE reported in the literature 25,46,47 . The warming signal emerges from natural variability sooner in the tropics than in high latitudes where SToE occurs decades later. We also find that midlatitude regions in the southern hemisphere may also experience novel climate conditions sooner than other latitudes. These regions are commonly identified as particularly vulnerable and with low financial and adaptation capacities 46 . The uncertainty in the estimates of StoE and TtA was evaluated using a bootstrap procedure (see "Methods"; Supplementary Tables S33-S35). The confidence intervals show the uncertainty in SToE to be highly skewed towards higher values and, in agreement with previous studies based on climate models' simulations, the uncertainty is larger for higher latitudes than for near-tropical regions 48 .
The TtA metric ( measures the transition time to reach a novel climate, i.e., the time available to adapt to a new climate. Taking 1880 as the reference year, TtA 2,1880 is close to 100 years for most latitude belts and continents (Tables 1 and 2 and Supplementary Tables S23-S25), and for some countries it can be as large as 130 years (Table 3 and  Supplementary Table S24). The TtA value decreased rapidly during the 20th century. By the mid-20th century (circa 1960), for most latitude belts, transitioning to a new climate (TtA 2,1960 ) took only about 2-3 decades and at the continental level about 3 to 4 decades, with few exceptions such as Europe and North America, where this threshold is reached a decade later. The SToE2 dates and TtA are slightly more heterogeneous at the country level but follow a similar pattern (Table 3 and Supplementary Table S25). The time to adapt was more than halved in the mid-20th century, with TtA 2,1960 values ranging from 30 to 60 years; e.g., the TtA 2,1960 value is 57 years for the UK and Germany, while it is as low as 28 years for China and Brazil.
To illustrate the pace of change in climate and the reduction in the time societies may have to adapt given the change in climate, three RCP scenarios were used to project TtA and StoE using 2020 as the reference year (Tables 1-3 and Supplementary Tables  S26-S32). These were selected to represent different levels of international mitigation policy: RCP8.5 represents a high warming scenario with no mitigation actions, while RCP4.5 and RCP2.6 can be viewed as approximations to the strict compliance of the Nationally Determined Contributions (NDC) commitments and the goals of the Paris Agreement 49 , respectively. Under the RCP8.5, the time available to adapt during this century becomes markedly shorter. The TtA values with respect to 2020 are about 12 years for Asia, Africa, and South America and 20 years for Europe and North America (Table 2 and Supplementary  Tables S28-S29). At the country level, the TtA 2,2020 values are around 20 years, with countries such as China and Brazil having values as low as 13 years (Table 3). Hence, in general, it would take about 10-20 years to experience mean temperatures that would be considered extreme by today's standards. It is also notable that the confidence intervals for TtA and SToE2 become increasingly narrower for higher warming trajectories (Supplementary Tables S33-S35). This illustrates the magnitude of the adaptation challenges that societies would experience under high warming and no climate policy scenario. Moreover, the TtA values get smaller through the 20th century as the reference period is updated, making it unlikely that adaptation processes could keep up with the pace of change in climate conditions. Achieving strict compliance of the NDC would provide about 10-15 years more for most countries to adapt, although, for countries like Germany and the UK, the TtA 2,2020 values could be about 25 years (Table 3). In contrast, under an emission scenario consistent with the goal of the Paris Agreement, the pressure to adapt would be much lower as the future and current climate would remain similar to that of 2020 during this century and the SToE2 2020 would take values beyond the present century (Table 3  and Supplementary Tables S31-S32).

Conclusions
We provided an observation-based regional attribution study that shows the existence of a dominant warming trend of anthropogenic origin in annual mean temperature series at different spatial scales, including continent and country level.   latitudes of the northern hemisphere. The results are consistent with estimates of regional warming produced by climate models and with attribution studies based on models' simulations 1,42,50 , providing an independent confirmation. The concept of time of emergence was modified to make it more informative for decision-makers. It is based on the natural temperature variability over the observed period, which is better able to reflect what societies are prepared for, than long-term natural variability used in the standard construction of the ToE metric. Our study uses the attribution results to compute all the necessary ingredients to obtain SToE and TtA, which are interrelated metrics that express: (1) the date when the climate signal would exceed a threshold of B times the standard deviation of experienced natural temperature variability; and (2) the number of years left to reach this new climate. The results show that the time to adapt to a novel climate has been rapidly decreasing since the late 19th century. While TtA had values near 100 years at the beginning of the 20th century, by the mid-20th century it decreased to about 3 or 4 decades. Then, in only 3 or 4 decades most parts of the world have experienced mean temperature values that would have been considered extreme realizations from the right tail of temperature distributions with respect to the 1960s. Projections for this century indicate that taking 2020 as the reference year, most parts of the world would experience in about 10-20 years a novel climate that would now be considered extreme if warming is not controlled by relevant climate policies. Moreover, TtA is projected to rapidly decrease during this century. These results have important implications about the time available to adapt and whether successful adaptation is feasible. The estimates based on the RCP4.5 scenario suggest that an intermediate international mitigation effort could help by providing about 10-15 additional years for adaptation. Under an emissions scenario consistent with the Paris Agreement, most countries, continents, and regions would not experience temperature conditions much different than current ones.

Methods
Statistical methods for detection and attribution of climate change. The use of time-series methods in detection and attribution studies has a long history in climate change literature 12,33,45,[51][52][53] . While these methods offer the advantages of not depending on the accuracy and performance of complex climate models, attribution is not possible without invoking a physical model that may be implicit in the statistical framework used and without assessing its physical consistency 1 . As described below, our attribution analysis is grounded on a zero-dimensional energy balance model and on the well-established result of stationary spatial patterns found in different generations of complex global atmosphere-ocean general circulation models [54][55][56] . Another challenge of attribution studies based on statistical techniques has been the correct identification of the time-series properties of climate and forcing variables. One of the main problems for obtaining empirical evidence about the influence of anthropogenic forcing using observed temperature records was the lack of adequate statistical tests to relate trending variables. While this was apparently solved with the introduction of cointegration techniques 51,52 , they impose important restrictions about the type of data-generating processes for temperature and forcing time series 57 . The analysis of the time-series properties of these variables using modern econometric tests indicates that cointegration techniques may find spurious relationships which generated a debate 12,58-60 . Here, we use an encompassing strategy based on a selection of cotrending tests that (1) makes results robust to the possible misidentification of the integration order and; (2) implies the existence of common breaks between the group of time series when a cotrending relationship is found. This is of importance because it provides strong evidence of causality in the relationships as recently argued 61 .
Zero-dimensional energy balance model and analyzing the influence of aggregate anthropogenic forcing at regional scales. The cotrending tests used here to analyze the influence of anthropogenic forcing are interpreted in terms of zero-dimensional energy balance model 12,62,63 and the assumption of stationarity in the spatial patterns of change frequently used in the literature [54][55][56] .
The foundations of time-series models to analyze the relationships between radiative forcing and temperature series have been established in previous publications 62,63 . The time-series models used in our study can be represented by where T t is global temperature, F t is a measure of the change in radiative forcing, α and γ are the intercept and slope parameters, respectively, and ε t is a stochastic noise process that represents high to low-frequency natural variability 12 . The structural model supporting this statistical model can be described by a simple twocompartment climate model 12,62,64,65 . It consists of an upper compartment (U) which represents mainly the atmosphere and the upper ocean, and a lower compartment (L) associated with the deep ocean. These components are thermally coupled and can be described by the following equations 62 : where C U and C L represent the heat capacity of the upper and lower compartments, respectively, and ΔT U and ΔT L are the changes in temperature in the respective compartments. F is the external forcing, and λ and β are the climate response and heat exchange coefficients. The heat capacity is much larger in the lower than in the upper compartment. The time constant of the response of the upper compartment is~4-9 years, while for the lower compartment this value ranges from 400 to 580 years 62,64 . The analysis in this paper relates to the transient climate response (TCR) which characterizes how the upper compartment responds to sustained increases in the external radiative forcing. The TCR is defined by S tr ¼ κ þ λ ð Þ À1 were κ is the heat uptake coefficient of the climate system. This coefficient relates the timedependent changes in surface temperature and external forcing via ΔT(t) = S tr F(t) and it is represented as γ in Eq. (1) 12 . The response of surface temperatures to external radiative forcing is dominated in the observed period by the short time constant of the upper compartment and the TCR. This provides a physical explanation of why temperatures and external forcing share a common nonlinear trend and common features such as co-breaks. By construction and following the physical model described above, the selection of variables used for two of the three cotrending tests 36,38 (TRF, WMGHG, and temperature) implies evaluating causality and the importance of anthropogenic greenhouse gas forcing. The simplest radiative forcing series included in the cotrending tests is WMGHG and if its nonlinear trend is common to TRF and temperature series, this indicates that it is the dominant driver of the secular movement in TRF and temperatures. Moreover, co-breaking is a necessary condition for cotrending in one of the tests used and the existence of common breaks provides strong evidence for causality 61 . The third cotrending test used was developed specifically to test for a linear stable relationship between temperature and TRF and complements the other cotrending results 37 .
The analysis of simulations from different generations of global coupled ocean-atmosphere general circulation models have shown strong support for stationary spatial patterns of change in variables such as surface temperature 54 . This result has been used for a variety of purposes such as producing scaling patterns or climate scenarios 56,66 , time of emergence 25 , and regional TCR estimates based on aggregate quantities 55 . The regional estimates of TCR can be estimated regressing regional temperatures on external forcing series: where β = φγ is the regional TCR, γ is the global TCR from Eq. (1), c is the intercept, F t is a measure of TRF, and T t,r is surface temperature in region r. It is important to note that if regional forcing (such as aerosols) significantly distorts the underlying global forcing trend, instead of only modulating it, there would be no cotrending between WMGHG, TRF, and regional temperature. Our results show that at the spatial scales and regions chosen for the analysis, regional forcing factors have not altered the secular movement imparted by global radiative forcing enough to preclude finding a common nonlinear trend between WMGHG and/or TRF and regional temperatures. Finding a cotrending relationship between TRF and regional temperatures helps to identify physically meaningful TCR, SToE, and TtA estimates. It also provides a rationale consistent with climate physics to use TRF to decompose regional temperature series into their trend and noise components, the latter including both natural oscillations and the modulation effects produced by regional external forcing.
Carrion-i-Silvestre and Kim cotrending and cointegration tests. Carrion-i-Silvestre and Kim 37 proposed joint tests of cotrending and cointegration, which were developed specifically for a linear stable relationship between temperature and TRF. These tests are quasi-likelihood ratio tests in a bivariate system of temperature and TRF. Under the null hypothesis, TRF and temperature are linearly related; under the alternative, they are not. These tests aim at offering valid inferences regardless of the true type of the trends in the data. TRF is allowed, but not assumed, to have a stochastic trend in addition to an obvious trend due to the well-mixed greenhouse gas forcing. Temperature is also allowed, though not assumed, to have any of a stochastic trend, a broken linear trend, a trend due to the well-mixed-greenhouse-gases forcing. We use theŜ 2 test, which tests for an exact linear relationship. It is constructed as follows. Suppose that y, x, and m are observed vectors of temperature, total radiative forcing, and well-mixed greenhouse-gases series: y = (y 1, …, y T )′, x = (x 1 , …, x T )′, and m = (m 1, …, m T )′. In addition, define regressors such as Þ¼ ð0; Á Á Á ; 0; 1; 2; Á Á Á ; T À ½TπÞ 0 with π 2 ð0; 1Þ. TheŜ 2 statistic combineŝ S þ 2 andŜ À 2 viaŜ 2 ¼Ŝ þ 2 Á 1ðα s ¼ 1Þ þŜ À 2 Á 1 À 1α s ¼ 1 ð Þ ð Þ , where 1(·) is an indicator function;α s ¼α, ifα À 1 j j> 1:96T À1=2 , and 1, otherwise; andα is a consistent estimator for the sum of autoregressive coefficients ofx 0 ¼ x-D x ðD 0 x D x Þ À1 D x x. TheŜ þ 2 statistic is computed as follows. Letβ andδ y be the ordinary least squares estimate from a regression of y on x and D x and letδ x be the ordinary least squares estimate from a regression of x on D x . Obtainû ¼ ½y À xβ À D xδy ;û 0 being a lower triangular matrix of ones. Then, estimate the long-run variance ofû t , the transpose of the tth row ofû byΣ ¼ T À1 Pû tþj , andΩ ¼Σ þΓ þΓ 0 , where κ(·) is a kernel function and b T is a bandwidth. PartitionΩ ¼ω yyωxŷ ω yxωxx ,Σ ¼σ yyσxŷ σ yxσxx , andΓ ¼γ yyγxŷ γ yxγxx .
Bierens nonparametric nonlinear cotrending test. The advantage of the test proposed by Bierens 36 is that the nonlinear trend does not have to be parameterized, it can be modeled as z t = g(t) + μ t where g t ð Þ ¼ β 0 þ β 1 t þ f t ð Þ; z t is a k-variate time series, μ t is a k-variate zero-mean stationary process and f(t) is a deterministic k-variate general nonlinear trend function that allow for structural changes. Nonlinear cotrending occurs when there exists a non-zero vector θ such that θ′f(t) = 0 where a'b is the inner product of two vectors a and b. Hence, the null hypothesis of this test is that the multivariate time series z t is nonlinear cotrended, and one or more linear combinations of the time series are stationary around a constant or a linear trend. This test is a cointegration test in the case when it is applied to series that contain unit roots. The nonparametric test for nonlinear cotrending is based on the generalized eigenvalues of the matrices M 1 and M 2 defined by ½ T À1 Þ, withβ 0 andβ 1 being the estimates of the vectors of intercepts and slope parameters in a regression of z t on a constant and a time trend. M 2 is defined as M 2 ¼ T À1 P T t¼m ½m À1 P mÀ1 j¼0 ðz tÀj Àβ 0 Àβ 1 ðt À jÞÞ½m À1 P mÀ1 j¼0 ðz t À j Àβ 0 Àβ 1 ðt À jÞÞ 0 , where m = T a with T the number of observations and α = 0.5. Solving M 1 À λM 2 j j¼0 and denoting the solution of the largest rth eigenvalue byλ r , the test statistic is T 1Àαλ r . The null hypothesis is that there are r cotrending vectors against the alternative of r − 1 cotrending vectors. The test has a non-standard distribution and the critical values are tabulated in Bierens 36 . The existence of r cotrending vectors in r + 1 series indicates the presence of r linear combinations of the series that are stationary around a linear trend and that these series share a single common nonlinear deterministic trend. Such a result indicates a strong secular co-movement in the r + 1 series.
Guo and Shintani consistent cotrending rank selection when both stochastic and nonlinear deterministic trends are present. Guo and Shintani 38,67 propose a model-free consistent cotrending rank selection procedure when both stochastic and nonlinear deterministic trends are present in a multivariate system. Their procedure defines two cotrending ranks: r 1 , which is the total number of linearly independent vectors in an m-variables system that can eliminate both stochastic and deterministic trends at the same time; and r 2 , which is the total number of linearly independent vectors that can eliminate the deterministic trend, regardless of whether they eliminate the stochastic trend at the same time. r 2 is called the weak cotrending rank and m − r 2 is the total number of common deterministic trends. These cotrending ranks can be estimated in paired or joint ways, but for large enough samples (T = 100) the two methods yield similar results. The proposed cotrending rank selection procedure allows for a wide variety of nonlinear trends, including breaks in the trend function and smooth transitions. In the case of breaks in the trend function, no common deterministic trends will be reported by this procedure if either: (1) there is no co-breaking (i.e., the breaks in the trend functions of the variables in the system are not common) or; (2) not all of the piece-wise trend slope coefficients are proportional between the trend functions. For cotrending to hold, co-breaking is a necessary, but not sufficient, condition. Here we report only the weak cotrending rank r 2 as the literature has shown temperature and radiative forcing series to be better represented as trend stationary processes with breaks in their trend function 3,39,43 . Let Z t denote an m-variate process for t = 1, …, T. Let λ 1 ≥λ 2 ≥ Á Á Á ≥λ m denote the eigenvalues of S À1 11 S 00 where S 11 ¼ T À1 P T t¼1 z t z 0 t and S 00 ¼ T À1 P T t¼2 ðz t À z tÀ1 Þðz t À z tÀ1 Þ 0 . The paired procedure selects r 1 and r 2 independently by minimizing each of VN 1 r 1 ð Þ ¼ À P r 1 i¼1λi þ f r 1 ð ÞC T =T and VN 2 r 2 ð Þ ¼ À P r 2 i¼1λi þ f r 2 ð ÞC 0 T =T 2 . The joint procedure selects both r 1 and r 2 by minimizing VN r 1 ; r 2 ð Þ¼ÀT a P r 1 i¼1λi À P r 2 i¼r 1 þ1λi þ f r 1 ð ÞC T =T þ f r 2 ð ÞC 0 T =T 2 with 0 < a < 1. C T = In(T), 2 In (In(T)), and 2 for the Bayesian information criterion (BIC), Hannan-Quinn criterion (HQ), and Akaike information criterion (AIC), respectively. f(s) is an increasing function in s. The BIC was used for the results presented here.
Regression-based method to estimate the TCR and natural variability around the warming trend. The results from cotrending tests and of the procedure to estimate the cotrending rank indicate that temperature series from all spatial scales considered share a single common trend which is imparted by TRF and dominated by WMGHG. Deviations from such a trend are consistent with a second-order stationary process and thus natural variability is deemed to remain stable during the sample period. Both the estimation of TCR and the natural variability around the trend use the results of the cotrending tests. First, we use ordinary least square to estimate the regression T t = α + βTRF t + ε t , in which β represents the TCR. The residuals ε t represent the variations around the warming trend produced by natural variability and other factors such as regional/local forcing. We use ε t to produce estimates of the natural variability around the warming trend (σ) required to calculate SToE and TtA. Given the time-series properties of ε t , σ is assumed to be invariant during the observational period and thus we use the full sample to estimate the value of σ to provide a measure of the variations in temperature for which societies are assumed to be prepared for and able to cope with.
ToE and the procedure to compute SToE and TtA. ToE has been commonly used in the climate change literature to estimate when the warming signal will become larger than the background climate variability 25 . ToE helps to put into context how extreme warming projections are with respect to climate variability and thus seek to identify regions where risk could be higher and faced sooner. Examples of the application of this concept are found in ecology 26 , climate policy, and reducing the risk of climate extremes 27 . Its application covers diverse spatial and temporal scales as well as climate variables 25,46,47,68 . These studies typically involve using simulations from complex physical models to approximate the climate signal and a measure of long-term climate variability.
The social time of emergence SToE is evaluated using the estimates of TCR and σ via the following steps: (1) obtain estimates of the warming trend for a given reference period; (2)  Bootstrap procedure for constructing SToE and TtA confidence intervals. The bootstrap procedure we employ consists of two steps: (i) we estimate the regression, T t ¼ α þ βF t þ ϵ t to obtainε t ¼ T t Àα ÀβF t . Then, we fit an auto-regression of order 2 onε t , which isε t ¼ c 0 þ c 1εtÀ1 þ c 2εtÀ2 þ η t . Since the ordinary least squares estimator used for an auto-regression is subject to bias, we correct it by a bootstrap procedure. That is, we subtract from the ordinary least squares estimate of (c 0 , c 1 , c 2 ) the bias computed as the difference between the average of bootstrap estimates and the original estimate. If the bias adjusted estimates imply a nonstationary autoregressive process, then we reduce the magnitude of bias adjustment until the implied autoregressive process becomes stationary; (ii) Using the biascorrected estimate of (c 0 , c 1 , c 2 ), we re-construct bootstrap samples ofε t and combine them withα þβF t to generate bootstrap samples of T t . For each bootstrap sample of T t , we compute and store SToE and TtA measures. The confidence intervals reported in Supplementary Tables S33-S35 are percentiles of the bootstrap distributions of SToE and TtA.
Data description and sources. The annual mean near-surface temperature anomalies for zonal means (latitude belts) were obtained from the GISTEMP dataset 69  The continent-level annual mean near-surface temperature anomalies were obtained from the Berkeley Earth Surface Temperatures (BEST; http:// berkeleyearth.org/) and; the NOAA GlobalTemp dataset (https://www.ncdc.noaa. gov/cag/global/time-series). The country-level annual mean surface temperature anomalies were obtained from the Berkeley Earth Surface Temperatures (BEST; http://berkeleyearth.org/) and the CRU TS v. 4.03 72 and the derived CY v. 4.03 dataset (https://crudata.uea.ac.uk/cru/data/hrg/). We use two radiative forcing datasets from the Goddard Institute for Space Studies of NASA (https://data.giss.nasa.gov/modelforce/). The main one used is that of Hansen 43 and for sensitivity analyses we also consider the Miller dataset 73 . The aggregated radiative forcing variables used are the well-mixed greenhouse gases (WMGHGs; carbon dioxide (CO 2 ), methane (NH 4 ), nitrous oxide (N 2 O), and CFCs), and the total radiative forcing (TRF), which includes WMGHG plus ozone (O 3 ), stratospheric water vapor, solar irradiance, land-use change, snow albedo, black carbon, reflective tropospheric aerosols, and the indirect effect of aerosols.
We use simulations from climate models included in the Atlas subset of the CMIP5 ensemble 50 to compare our regional TCR estimates with those that can be inferred from current climate projections. The Atlas subset includes a single realization per model in the CMIP5 ensemble. We use the historical simulations and those produced under the RCP4.5 scenario. The complete list of models is given in Table 1 of the Annex of the Working Group I AR5 report 50 Fig. 1 are reported) and obtained the TRF value for 2011 from the RCP emission scenario database (http://www.pik-potsdam.de/~mmalte/rcps/index.htm) 44 . TCR is calculated as TCR ¼ Δ . Supplementary Tables S20 to S22 provide the TCR estimates (ensemble mean and 95% confidence intervals) for each region based on the IPCC's Atlas subset. The results show that the estimates of TCR based on the Atlas subset of the CMIP5 and those presented here are similar for all regions. Moreover, in all cases, our estimates fall within the confidence intervals of the models' projections.

Data availability
No datasets were generated during this study. All data used in this study are publicly available from the sources cited in "Methods" and at https://doi.org/10.6084/m9. figshare.13380779.v1.