REVISITING THE SPATIOTEMPORAL CHARACTERISTICS OF PAST AND FUTURE GLOBAL WARMING

It is still an open question, which processes lead to the spatiotemporal specifications of observed near-surface temperature changes over recent decades. Here, we contribute to this debate by investigating a large number of theorybased atmospheric fields referring to the radiation and energy budget and to atmospheric dynamics that may serve as predictors for local temperature changes. The predictors are linked to temperature trends from reanalysis and climate model data, using a sophisticated spatial and temporal statistical model. Temperature changes since the mid-20th century exhibit distinct regional and seasonal differences. After 1990, the near-surface warming rate is more enhanced over landmasses rather than oceans and roughly increases with latitude in both hemispheres. While none of the considered predictors solitarily accounts for the spatial heterogeneity of recent temperature trends, their linear combination largely reproduces the observed cooling pattern during the mid-20th century and the enhanced warming pattern after 1990. This excludes high-altitude areas, sea ice margins and upwelling regions where local feedbacks and nonlinear processes prevail. The leading predictors pertain to radiative processes, especially downward longwave radiation, and changes in sensible heat fluxes. In the low latitudes, dynamical processes such as temperature advection and energy flux divergence also play a role. Until the end of the 21st century, the warming rate and its ocean-land contrast steadily increase. The underlying mechanisms are the same as the ones already established in present-day climate, but near-surface temperature follows more straightly the imposed greenhouse gas scenario. Climate models have different skills in reproducing the observed trend pattern but exhibit more or less the same mechanisms of temperature control.


Introduction
Global warming represents a major indicator of man-made climate change and a primary mechanistic explanation for the human interference with climate (IPCC 2007(IPCC , 2013. Therefore, a detailed physical understanding of the causes of observed temperature trends is crucial for creating public confidence in the scientific line of argument from greenhouse gas emissions over radiative forcing to anthropogenic climate change. Yet, the spatiotemporal heterogeneity of temperature changes over recent decades does not suggest a simple link to monotonously and uniformly rising greenhouse gas concentration (Stott et al. 2000;IPCC 2007IPCC , 2013. This raises the question of how increasing greenhouse gas and aerosol concentrations affect near-surface temperature via feedbacks in the radiation budget, cloudiness, surface processes and atmospheric dynamics. In terms of the temporal component, much attention has been paid to the time series of near-surface global-mean temperature. Stott et al. (2000) have been among the first to detect the main drivers of low-frequency variability in the observed global-mean temperature time series by a series of climate model simulations with various combinations of natural and anthropogenic forcings. They concluded that a realistic combination of natural and anthropogenic drivers is necessary to reconstruct the observed time series, with enhanced greenhouse forcing clearly dominating since the 1980s. Such 20 th -century climate hindcasts have also become a standard approach in the latest IPCC reports dedicated to attribute observed climate characteristics, mainly regional-mean temperature time series, to either natural or man-made components of climate variability (IPCC 2007(IPCC , 2013. There is a broad consensus that the early warming period of the 20 th century emanates from a combination of natural and anthropogenic drivers, that the slight cooling during the 1940-1970 period is mainly caused by increasing aerosol emissions from human sources, especially over Northern Hemisphere landmasses, and that the substantial temperature rise since the 1980s is predominated by enhanced greenhouse warming (e.g. PaetH and FeicHter 2006;IPCC 2007;knutSon et al. 2013;IPCC 2013, PaetH 2015, with a presumed hiatus after the turn of the millennium due to internal variations in the coupled ocean-atmosphere system (coHen et al. 2012).
Concerning spatial heterogeneity, the pattern of temperature changes is more difficult to study because some areas of the globe are barely covered by meteorological observations (cf. HarriS et al. 2020). Based on a combination of data from different sources, the 5 th Assessment Report of the IPCC draws a rather comprehensive picture of observed changes in near-surface temperature (IPCC 2013). Supplemented by climate model experiments, it can be concluded that the temperature trends since 1901 feature two systematic spatial structures: an equatorto-pole gradient and an ocean-to-land contrast ( Jain et al. 1999). This is particularly expressed in the Northern Hemisphere. Wallace et al. (1995) have been the first to suggest that the landmasses (oceans) are relatively warm (cold) when the Northern Hemisphere experiences a positive temperature anomaly. They named this pattern cold ocean -warm land (COWL) and explained the contrast by the different heat capacity of land and ocean surfaces, favouring a larger temperature response to advection processes over land. Since then, several authors have shown that, besides near-surface temperature, the COWL pattern leaves a mark in various lower-tropospheric climate features and, hence, represents a leading pattern of climate variability (Wu and StrauS 2004;SHi and BueH 2012). Broccoli et al. (1998) also saw the different thermal inertia of ocean and land surfaces as the major cause for stronger warming trends over the continents. Other explanations pertain to changes in extratropical circulation modes and the related advection of anomalously warm and cold air masses, respectively (Quadrelli and Wallace 2004;rautHe and PaetH 2004;PaetH and Pollinger 2010). More recently, Boer (2011) stated that the argument of different thermal inertia may be too simple to account for the substantial spatial heterogeneity of observed temperature trends. Instead, he referred to the importance of local feedbacks in the context of radiative forcing, surface processes, cloudiness and atmospheric circulation. A crucial mechanism for the maximum warming rates over high-latitude and high-altitude regions is given by the ice-albedo feedback as a result of reduced snow and sea-ice cover and melting glaciers (IPCC 2007(IPCC , 2013. Although Wallace et al. (1995) originally interpreted the COWL pattern as being an expression of internal climate dynamics, the spatial pattern of observed and simulated warming trends is often referred to as a fingerprint of anthropogenic climate change and used for detection and attribution studies (e.g. Broccoli et al. 1998;Jain et al. 1999;PaetH and HenSe 2001;van oldenBorgH et al. 2013;PaetH et al. 2017). In addition, the differential warming of the Earth is supposed to af-fect climate phenomena that arise from spatial temperature gradients such as jet streams, atmospheric waves and monsoons (Sutton et al. 2007;Jain et al. 1999, PaetH andPollinger 2019).
While the pattern of past and future temperature changes plays a crucial role for ongoing climatological research, but also for regional climate impacts and necessary adaptations, its spatiotemporal heterogeneity still lacks a comprehensive explanation (Sutton et al. 2007). In the current study, we pursue three major objectives: (1) we aim at decoding the trend patterns of near-surface temperature from observational, reanalysis and climate model data in the light of mechanisms involved in the process chains between radiation budget, surface energy fluxes, atmospheric dynamics, land surface and cloud feedbacks. These mechanisms are represented by meteorological fields that are selected according to thermodynamic theory (see section 2) and assessed by means of a statistical model approach that allows for identifying the relative importance of these mechanisms in different regions. (2) We compare the prevailing mechanisms in the past and future in order to gain insight into the temporal signature of the processes that are related to the warming of our planet. (3) We evaluate climate models in terms of their ability to reproduce the observed temperature trend patterns and the underlying processes as derived from reanalyses.
The research approach is guided by two scientific hypotheses: (1) the spatiotemporal heterogeneity of observed and simulated temperature changes cannot be satisfactorily accounted for by a single process but by a linear combination of a limited number of well-founded processes. (2) The underlying processes of global warming differ between past and future periods.
In the next section, we delineate the theoretical background of temperature changes, the resulting statistical methodology and the analysed model and observational data sets. Section 3 is dedicated to the results. The findings of the present study are discussed in section 4 against the background of previously published work. The main conclusions are drawn in section 5.

Considered data sets
The main part of this study refers to the NCEP/ NCAR reanalysis (hereafter NCEP) that is globally available in 2.5° resolution and extends from 1948 until 2019 (kalnay et al. 1995). The analysis is based on the output variables 2 m temperature, 10 m zonal and meridional wind velocity, surface sensible and latent heat flux, surface downward short-wave and long-wave radiation, surface upward short-wave and long-wave radiation (Tab. 1). Reanalyses provide complete and multivariate meteorological data sets with a high level of inherent physical consistency, nudged to available observational data. Yet, we are aware that reanalyses are not observations in a stricter sense, but the required meteorological fields are not globally available from station data or gridded data sets. For this reason, the temperature trend patterns from the NCEP reanalysis are compared with the Climatic Tab. 1: Considered data sets and atmospheric fields (T = temperature at 2 m height, u = zonal wind at 10 m height, v = meridional wind at 10 m height, for additional variable symbols see Tab. 3).  (HarriS et al. 2020). The CRU data set represents a gridded temperature data set in 0.5° resolution that is based on a large number of meteorological station observations (Tab. 1). Despite the lower horizontal resolution compared with the most recent ERA5 reanalysis, the NCEP reanalysis has been chosen because it provides continuous data since 1948, capturing the major cooling and warming phases of the post-war period. Trend periods in the early 20 th century would also be interesting to study in this context, e.g. from 20 th -century reanalyses. However, temperature data during the first half of the last century are characterized by substantially lower quality (Brönnimann 2009) and, hence, may imply spurious effects in the temperature trends. The second part of our study includes the same meteorological fields from coupled climate models (Tab. 1). As climate change detection in a multi-model ensemble framework is not a primary goal of this analysis, we rely on a small number of selected climate model simulations that are presumably characterized by a different ability to reproduce observed temperature trend patterns. This selection goes back to the work from ring et al. (2019) who have ranked climate models from Coupled Model Intercomparison Project (CMIP) version 3 and 5 with respect to exactly this feature.

Data set
The main focus is on the German MPI-ESM-LR Earth system model from the Max-Planck Institute for Meteorology (giorgetta et al. 2013). This global model has a T63 spectral resolution (~2°) that has been interpolated to the 2.5° grid of the NCEP reanalysis. The simulation covers the period 1850-2100, using observed natural (Earth orbit, solar variability, natural and volcanic aerosols) and anthropogenic (well-mixed greenhouse gases, sulphate aerosols, land For model inter-comparison, two additional coupled climate models are considered, one from the CMIP3 framework (IPCC 2007) and one from CMIP5 (IPCC 2013). According to ring et al. (2019), they are assumed to have a lower skill to reproduce the observed temperature trend patterns. In contrast to MPI-ESM-LR, the CMIP5 model provides SO 4 aerosols on a constant yearly cycle, but no other aerosols. The CMIP3 model uses the 20C3m scenario until 2000 and the SRES A1b scenario from there on. Land use changes and volcanic forcings are not included and solar forcing exclusively for 20C3m. Only SO 4 aerosols are prescribed. They serve as a negative example and remain incognito. The interesting question is whether the investigated climate models also differ in terms of the mechanisms leading to temperature changes or just have different sensitivities to the imposed forcings.

Theoretical background and predictor selection
A key element of our approach is the selection of adequate meteorological fields that serve as predictor variables for the spatial and temporal characteristics of near-surface temperature changes, representing the predictand in the statistical model. The predictor selection must satisfy two constraints: first, the link between predictors and temperature changes must be physically based in accordance with thermodynamic theory and, second, the number of predictors must be all-embracing with respect to the required information, but at the same time small in order to avoid overfitting.
The total time derivative of temperature T emanates from the first law of thermodynamics, excluding chemical reactions (Wallace and HoBBS 2006). When applying this to near-surface temperature T 0 in our case the temperature at 2 m above ground, vertical wind velocity vanishes. Thus, for a local change in time, e.g. in a model grid box, temperature variations are caused by horizontal advection, local changes in diabatic heating and the divergence of energy fluxes in the diabatic term with horizontal wind vector w H and heat capacity cp. The diabatic term refers to warming and cooling effects arising from radiative processes, sensible heat fluxes and phase transitions of water with the mass-specific energy fluxes for radiation q, sensible heat h and latent heat l (with unit J). The negative signs occur because energy fluxes are typically defined from the surface into the atmosphere, with a positive sign marking a loss of energy at the surface.
Concerning the third term on the right hand side, it must be distinguished between energy fluxes that simply pass through the considered atmospheric layer without heat release, and energy fluxes that transfer a part of their energy to the intrinsic energy of the atmospheric layer and, hence, lead to temperature rise -or decrease when the incoming heat flux is smaller than the outgoing (Hantel and HaimBerger 2016). In the prognostic equation for temperature this is taken into account by considering the divergence of the energy fluxes. Note that the divergence of latent heat fluxes l within an atmospheric layer directly affects the atmospheric water vapour content but not temperature, as long as no phase transitions of water occur. As condensation and evaporation barely prevail at 2 m height, we do not expect a strong impact of latent heat flux divergence on near-surface temperature.
In high-resolution climate models with time steps of 1-2 minutes, local time changes in the second term are very small and usually not taken into account (cf. Bott 2012). In the present study, we pursue a climatological approach on the basis of multi-decadal temperature variations and, hence, include long-term changes of surface energy fluxes H and L into the atmosphere and surface net radiation Q (now with standard output unit Z). In order to improve our process understanding, Q is further differentiated into surface global radiation I, surface albedo A, downward long-wave radiation at surface G (a.k.a. greenhouse effect) and outgoing surface long-wave radiation E.
In total, we use 11 potential predictor fields to explain the spatiotemporal characteristics of observed and simulated temperature changes: horizontal advection, 3 terms of energy flux divergence, sensible and latent heat fluxes, and 5 variables referring to the surface radiation balance. The predictor fields themselves arise from a large number of processes in the climate system and, hence, stand for the complexity of mechanisms leading to warming and cooling tendencies at the surface: the first four predictors involve changes in atmospheric circulation. Sensible and latent heat fluxes are affected by changes in surface net radiation and skin temperature, vegetation cover, soil hydrology and atmospheric turbulence. The four components of surface net radiation reflect a wide variety of processes involved in climate variability and change: global radiation includes signals of solar variability, volcanic eruptions and variations in cloudiness and aerosol absorption. Surface albedo relates to land cover changes and to the spatial ex-tent of sea ice, snow and glaciers with feedbacks from temperature. Downward long-wave radiation is influenced by greenhouse gas concentrations, atmospheric water vapour content and cloudiness, while outgoing long-wave radiation is a direct nonlinear function of near-surface temperature. For the temporal statistical model, the time series of atmospheric greenhouse gas concentrations, expressed as CO 2 -equivalents, is used additionally as a global-mean predictor, based on observations until 2005 and the RCP8.5 emission scenario until 2100 (van vuuren et al. 2011).

Statistical model
While temperature variations at the scale of temporal discretization in climate models and reanalyses (i.e. minutes) emanate from the equation described in subsection 2.1, we are interested in the statistical relationships between the selected theory-based predictors and spatial and temporal specifications of temperature changes at the climatological time scale (i.e. years and decades). The statistical transfer functions between the predictor fields and temperature changes are determined by means of a stepwise linear multiple regression model with a bootstrapping approach for cross validation (WilkS 2006). The statistical model provides information on the explained variance and statistical significance of each predictor and of the entire model, the relative importance (ranking) of the predictors, and the number of robust predictors as approved by cross validation. Especially this latter information is a strong argument for such a classical statistical model rather than using machine-learning approaches which represent efficient optimization algorithms but do not provide detailed information on the role of each individual predictor (cf. Pollinger et al. 2017). A detailed description of the statistical model is given by PaetH (2011). The assumed linearity is certainly a major constraint that needs to be scrutinized in the light of the gained results. We are aware that correlation and regression analysis is not dedicated to detect physical mechanisms nor causality between the considered variables. While the physical processes can be assumed to be reasonably represented in the analyzed reanalysis and climate model data, our main goal is to generate process understanding at the climatological time scale by explaining the space and time structure of global warming in the light of meteorological fields that, from the per-spective of theory, are directly related to near-surface warming. In this context, it is not obstructive that temperature is not only affected by the predictor variables but also feeds back on most of them.
The statistical model is applied to spatial patterns (spatial approach) and to time series at individual grid boxes (temporal approach). Both approaches use various time windows and seasons in the past and future in order to assess the many-sided facets of observed and simulated temperature changes. The spatial approach is based on a large sample of 10,244 grid boxes. As the number of predictors is much smaller (= 11), 1,000 grid boxes can be retained for cross validation. The temporal approach uses time windows of 30 years or longer that are suitable to distinguish between noise and potential signals of natural or anthropogenic climate drivers (Santer et al. 2011;Paxian et al. 2013). Here, we reserve only 5 years for cross validation. The bootstrap selection of retained grid boxes and years, respectively, is randomized and repeated 50 times. It has been found that this is by far sufficient to obtain a stable solution. The results illustrated in section 3 represent either averages or sums over the 50 iterations of the statistical model.
When spatial means and correlations are computed the grid boxes are weighted by their spatial extent using a trigonometric function. All predictors are standardized to equilibrate differences in spatial and temporal variability. For the temporal approach of the statistical model, all time series are de-trended with the aim of avoiding spurious correlation due to long-term trends. The assessment of statistical significance is based on standard t-test (correlation coefficients and trends) and F-test (multiple regression model) statistics (WilkS 2006). For the spatial approach, the degrees of freedom are reduced due to the substantial spatial autocorrelation of temperature and, to a lower extent, of the various predictor variables. The reduced degrees of freedom Φ r are computed by with Φ = n -2 = 10,242 arising from the original grid box resolution and q representing the global-mean autocorrelation of de-trended temperature time series between neighboring grid boxes (WilkS 2006). The largest autocorrelation is found for annual-mean temperature ( q ), leading to Φ r = 600. For a conservative test design, this applied to all considered variables and seasons.

Spatiotemporal characteristics of global warming from NCEP
The spatial pattern of temperature trends (specified as linear temperature change per decade) during the NCEP periods 1948-2019, 1948-1977 and 1990-2019 are shown in Fig. 1. Concerning the long period (top panel), increasing temperatures clearly dominate over the globe. The strongest warming has occurred over the polar regions with more than 1°C per decade in the Barents Sea. There are also several spots with minor cooling tendency, especially over the mid-latitude oceans and over some land areas such as the Sahara, the Andes and parts of China. The trend patterns in boreal spring and autumn are very close to the one for annual-means, whereas Antarctica has experienced a clear cooling in austral summer since 1948 of up to 0.4°C per decade (not shown). The temperature trend patterns from the CRU data set are quite similar but less heterogeneous in space (not shown). The spatial correlation coefficient over the landmasses covered by NCEP and CRU amounts to 0.6 which is quite high given the large sample of more than 2,700 grid boxes.
The mid-20 th century (middle panel) was marked by a hemispherically asymmetric pattern with cooling tendency over most parts of the Northern Hemisphere mid-latitudes, especially over North America, Europe and Asia. At the same time, most areas of the Southern Hemisphere underwent a noticeable warming trend. Over and around Antarctica, the trend pattern is very noisy with positive and negative trends in close vicinity. During the last 30 years (bottom panel), a massive warming prevails, revealing two basic structures: an equator-to-pole gradient with maximum warming over the polar regions, and an ocean-to-land contrast (cf. Jain et al. 1999). Nonetheless, some regions still experienced decreasing temperatures. This mainly holds for parts of the South American highlands, some desert regions and mid-latitude ocean basins. Temperatures from the CRU data set confirm the trend patterns from NCEP, but the cooling before 1977 is less pronounced in CRU (not shown).
In order to assess systematic differences in temperature trends between land and ocean areas, Tab. 2 lists the spatial-mean trends from NCEP for different trend periods and seasons. Over the whole period from 1948 to 2019, only positive trends have occurred. Most trends are statistically significant and larger over the continents rather than oceans. Yet, the difference is small. Between 1948 and 1977, negative trends prevail but they are close to zero and not at all significant. The cooling has been a clearly continental phenomenon. Since 1990, our planet has warmed substantially and much more over the landmasses than over the oceans, supporting the postulation of the COWL pattern. The highest warming rate is found over the land areas in boreal autumn with 0.44°C per decade. Anyway, not all temperature trends since 1990 are statistically significant.

Predictor fields from NCEP
The long-term climatological mean of all predictor fields is depicted in Fig. 2, along with the climatological pattern of near-surface temperature as the predictand in our statistical model approach. Most annual-mean patterns are well-known and easily interpretable: the temperature pattern (T2M) mirrors the equator-to-pole gradient of surface net radiation and topography over the landmasses. The strongest sensible heat fluxes (H) into the atmosphere occur over dry land areas, negative fluxes from the atmosphere to the ground prevail in polar regions. Oceanic regions, the inner tropics and parts of the mid-latitudes barely exhibit any sensible heat fluxes. Instead, available radiative energy is almost entirely transformed to latent heat fluxes (L) because ample sources of water vapor are given from open water surfaces or dense vegetation cover. Solar irradiation (I) peaks over the low latitudes with a relative minimum in the vicinity of the intertropical convergence zone (ITCZ). The minimum occurs over the poleward mid-latitudes in both hemispheres with a slight increase towards the poles themselves due to the polar day phenomenon. Albedo (A) is largest over the regions covered by snow and ice, i.e. polar caps and high-mountain areas, and over deserts. Terrestrial emission (E) is a direct function of surface temperature and, hence, largest in the low latitudes and close to zero over Antarctica. A very similar pattern is given for the greenhouse effect (G), i.e. the downward long-wave radiation that depends on atmospheric greenhouse gas concentrations, water vapor content and cloudiness. The surface net radiation (Q) results from the interaction of the components of the surface radiation balance and is positive (negative) over the low (high) latitudes. The remaining four predictor fields are tied to near-surface atmospheric circulation. As zonal and meridional wind velocity are equally included in the calculation of the divergence terms of sensible heat flux (-div(H flux)), latent heat flux (-div(L flux)) and radiation flux (-div(Q flux)), these patterns look quite similar, but exhibit different amplitudes. The negative sign implies that a positive value in the maps indicates a net warming effect in the given region. Heat sources emerge in the ITCZ region, along the northern and southern equatorial current and over warm extratropical ocean currents, particularly the Gulf Stream and Kuroshio. The divergence of sensible heat fluxes is minor compared with the other energy fluxes, except for the polar regions. The last predictor field pertains to temperature advection (-adv(T)). Warm air advection prevails in the southern midlatitudes, especially on the western sides of the continents, and in monsoon regions. Cold air advection is a typical characteristic over cold ocean currents, oceanic upwelling regions and the polar caps.

Period
The changes per decade over the 1990-2019 period are displayed in Fig. 3 for the eleven predictor fields and 2 m temperature, using the same arrangement as in Fig. 2. The temperature trend pattern has already been dealt with in the previous subsection. Sensible heat fluxes have decreased since 1990 in most regions, while latent heat fluxes have increased over many oceanic regions, yet the trend pattern is quite noisy. The trend patterns of solar irradiation and greenhouse effect are negatively correlated, probably due to enhanced cloudiness that reduces short-wave irradiance but favors downward longwave radiation at the surface. Albedo has become weaker along sea-ice margins and in some highmountain regions. Outgoing long-wave radiation directly reflects the trend pattern of near-surface temperature. The trend patterns of surface net radiation and of the four predictor fields that are related to atmospheric circulation are spatially quite incoherent. There is a tendency towards stronger warm air advection in higher latitudes and more heating from convergent energy fluxes in the tropical oceans.

Spatial statistical model applied to NCEP
As a first glimpse on the spatial link between annual-mean temperature trend patterns and predictor fields, Tab. 3 shows the spatial correlation coefficients with respect to all three trend periods. Note that the sample comprises 10,224 grid boxes and, hence, the critical value for statistical significance is very small, even when accounting for the spatial autocorrelation in the temperature field (cf. subsection 2.3). The closest link is found between the trend patterns of temperature and surface long-wave emission, with correlation coefficients of about 0.9. However, this statistical relationship is trivial because emission is a direct nonlinear function of temperature (E = σ · T 4 , with Stefan-Boltzmann constant σ). For this reason, upward long-wave radiation is omitted from the list of potential predictors in the multiple regression model, noting that there are also some more indirect feedbacks between temperature and the other predictors. Quite high positive (negative) correlations exist with changes in downward long-wave radiation (albedo). This is in perfect agreement with theoretical considerations. Variations in sensible heat fluxes, surface net radiation and divergence of Q and H flux vector account for up to 10 % of the total variability of temperature trend patterns. The remaining four predictor fields only play a minor role. This picture is identical for all seasons (not shown). In summary, none of the selected predictors alone provides a satisfactory explanation for the spatial heterogeneity of temperature trends from NCEP (nor CRU).
In the next step, we deploy a spatial multiple regression model that is based on the ten remaining predictor fields and the temperature trend pattern as predictand. The main results are summarized in Tab. 4 for different trend periods and seasons. The total explained variance of the statistical model amounts to 80-95 % of the spatial variability of the cooling and warming patterns since 1948, at an error level of α = 1%. The explained variance is systematically higher in boreal winter rather than summer. The step-wise multiple regression model is combined with a bootstrapping approach for cross validation, dedicated to identify the ranking of predictors and the number of robust predictors. For this purpose, the root mean square error (RMSE) between the predictand and the predicted values from the model is computed for the retained independent data (here 1,000 out of 10,224 grid boxes, cf. subsection 2.3). The maximum number of robust predictors is given by the RMSE minimum in the spectrum over all predictors. Tab. 4 reveals that all ten predictor fields are identified as robust predictors. On the one hand, this is surprising because some of the predictor fields have similar spatial structures (cf. Fig. 3) and, hence, are subject to multi-collinearity. On the other hand, the total number of predictors is very small compared with the sample size and heterogeneity of the patterns. Therefore, every predictor contributes valuable information, at least for a sub-region of the

1990-2019
Sensible predictand. In most analyzed periods and seasons, downward long-wave radiation represents the leading predictor, sometimes it is albedo or surface net radiation. In summary, the total explained variance of the multiple regression model is substantially higher than the individual contributions of every single predictor (cf. Tab. 3). This implies that either more than one mechanism accounts for near-surface temperature changes or that different regions are affected by different mechanisms. This question will be picked up in the next subsection. Fig. 4 demonstrates that the 1948-1977 cooling pattern and the 1990-2019 warming pattern can indeed be reproduced by the spatial statistical model with reasonable accuracy: the original temperature trend patterns from NCEP (Fig. 1) and the ones predicted by the model (top panels in Fig. 4) are largely coherent. The model slightly overestimated the Northern Hemisphere cooling in the mid-20 th century and the recent warming over most landmasses (middle panels). For the overall ranking of the predictor fields, a counting statistic is used over all five seasons and 50 iterations of the model (bottom panels). It is obvious that the greenhouse effect prevails as leading predictor in both periods. During the cooling phase, the second rank is occupied by either sensible heat fluxes or solar radiation, while the latter dominates during the recent warming phase. The subsequent ranks are assigned to the components of the surface radiation balance. Since 1990, temperature advection has played a more important role than during the decades before. Changes in latent heat flux usually occupy the last rank, as to be expected (cf. subsection 2.2).

Temporal statistical model applied to NCEP
With the temporal statistical model, we aim at tailoring the model more to the local characteristics of temperature changes. This model is based on detrended time series (not patterns) of the predictand and the ten predictors at every grid box plus the global-mean time series of increasing CO 2 -equivalent greenhouse gas concentrations. The results for annual-means are mapped in Fig. 5 for the mid-20 th century cooling phase and the recent warming period. First of all, it must be stated that both time windows exhibit very similar results: the temporal statistical approach accounts for almost 100 % of total temperature variability in most grid boxes (top panels), implying that the predictor selection is adequate and the assumed linearity does not represent a major constraint. Interesting insight is given by the regions where the explained variance partly goes down to 60 % minimum: in some high-mountain areas, along sea-ice borders and in some oceanic sectors, especially in upwelling regions, the model is less successful because additional mechanisms or nonlinearities play a larger role in the temporal evolution of nearsurface temperature.
The number of robust predictors varies between 4 and 11 (middle panels). Note that the 11 predictors are numbered from 2 to 12 because the long-term mean as the simplest prediction is assigned number 1. On average, the number of robust predictors is higher in the low latitudes compared with extratropical landmasses and polar regions. As the explained variance is still high in the latter areas, it means that the identified predictors make large contributions to total temperature variability. Downward long-wave radiation is the main driver of cooling and warming in most regions across the planet (bottom panels). However, in some oceanic sectors sensible heat fluxes are more relevant. Dynamical predictors, i.e. divergent energy fluxes, play an important role in the low latitudes, particularly in monsoonal regions. The time series of increasing greenhouse gases prevails in a very limited number of grid boxes in the tropics and subtropics, but is more present since 1990 than in the mid-20 th century. Different seasons barely differ in terms of the findings illustrated in Fig. 5 (not shown). Fig. 6 displays the ranking of the predictors over 50 model iterations and their individual contributions to total explained variance of the statistical model. It has been found that the predictor ranking differs more from region to region than between time periods or seasons. Therefore, Fig. 6 distinguishes between regional averages over the Northern Hemisphere polar region (60°-90°N), the mid-latitudes (30°-60°N) and the low latitudes (0°-30°N). North of 60°N, the greenhouse effect clearly stands out as the leading predictor, followed by changes in sensible heat fluxes and the remaining components of the surface radiation balance. The dynamical predictors occupy the last ranks. Most of the explained variance accounted for by the statistical model is assigned to the leading predictor (no. 2, almost 80 %), while the climatological mean (no. 1) has no predictive skill at all. In the Northern Hemisphere tropics and subtropics, changes in surface radiation still play a dominating role but dynamical changes are much more prevailing. The time series of increasing greenhouse-gas concentrations (GHG in Fig. 6)  contributions to total explained variance are more equally distributed among the predictors. The midlatitudes represent a transition zone with respect to the predictor ranking. We also repeated the statistical model approach without de-trending the time series: in this case, the global-mean time series of enhanced greenhouse-gas concentrations is always one of the leading four predictors (not shown).

Temporal statistical model applied to climate models
In the last step of this study, the temporal statistical model is applied to output from coupled climate models, covering the 1948-2100 period. The main focus is on the MPI-ESM model because it is state-ofthe-art in terms of the included processes and forc- ing components in the past and future (cf. subsection 2.1). Thus, it is assumed to reproduce the patterns and mechanisms of temperature changes as derived from the NCEP reanalysis. In order to assess differences between climate models, we also investigate two other climate models that are characterized by lower complexity and have been found to have less skill with respect to the observed temperature trend patterns (cf. ring et al. 2019).
For the MPI-ESM model, the annual-mean decadal temperature trends over different periods in the past and future are listed in Tab. 5. The climate model produces slightly higher warming rates during the 1948-2019 and 1990-2019 periods compared with NCEP. The mid-20 th century cooling is not simulated, but the warming is much smaller than during any other time window and statistically not significant. Until the end of the 21 st century, the temperature rise accelerates continuously, reaching up to 0.71°C per decade over land areas. The ocean-to-land contrast in the global warming pattern also increases until 2100. Fig. 7 shows the global-mean time series of nearsurface temperature from MPI-ESM and the trend patterns over two time spans in the past and two in the future. Compared with the NCEP reanalysis (cf. Fig. 2), it is obvious that the cooling phase during the 1948-1977 period looks quite different: there are several areas with decreasing temperature but the hemispheric asymmetry is not reproduced by the climate model. In contrast, the warming patterns since 1990 look very much the same in NCEP and MPI-ESM, including the equator-to-pole gradient and the oceanto-land contrast. In the middle and late 21 st century, this basic structure is more and more apparent and near-surface heating becomes a circumglobal phenomenon. The comparative climate models with lower complexity reproduce only one aspect of the observed temperature changes, i.e. the maximum heating of the Arctic Ocean (not shown). Other main structures such as the ocean-to-land contrast are not simulated, and the climate model from the CMIP3 family still projects major cooling areas over Northern Hemisphere landmasses during the 2071-2100 period.
The results of the temporal statistical model that have been derived from the NCEP reanalysis (cf. Fig.  5), are also found for MPI-ESM and for the other two climate models (not shown). The only difference is that the global-mean time series of enhanced atmospheric greenhouse conditions prevails more and more as leading predictor towards the end of the 21 st century, especially in the low latitudes. This is also illustrated by the predictor matrices for MPI-ESM ( Fig. 8): during the last 30 years of our century, the mechanisms of surface warming reveal the same basic structure as found in NCEP during the 1948-2019 period (cf. Fig. 8), with the components of surface radiation balance being the most relevant drivers and dynamical processes being more important in low rather than high latitudes. However, future temperature time series follow more straightly the monotonous forcing by increasing atmospheric greenhousegas concentrations. The other two climate models exhibit basically the same predictor ranking as in NCEP and MPI-ESM (not shown).

Discussion of results
During the 1948-1990 period, NCEP exhibits the well-known temporal features of global-mean near-surface temperature changes: a slight cooling phase until the 1970s and a considerable temperature increase since the 1980s (cf. IPCC 2007(cf. IPCC , 2013. Both phenomena have been interpreted by many authors in the light of natural climate drivers and enhanced anthropogenic greenhouse gas and aerosol concentrations (e.g. Stott et al. 2000;PaetH 2015). When comparing the NCEP temperature time series averaged over the landmasses with the one from the station-based CRU data set, it becomes evident that NCEP and CRU are in excellent agreement with each other in terms of the long-term trend and even interannual variability, but NCEP overestimates the cooling during the mid-20 th century due to a warm bias in the 1950s (cf. SimmonS et al. 2004).
Over the landmasses, the trend patterns from NCEP and CRU are quite coherent (cf. SimmonS et al. 2004). In addition, NCEP mainly reproduces the trend patterns illustrated by IPCC (2013, Fig. 2.22) Tab. 5: Annual-mean temperature trends from the MPI-ESM climate model in different periods averaged over the globe, the land area and the ocean area (unit is °C/decade). Bold numbers indicate linear trends significant at the 5% level.

Period
Globe with mid-20 th century cooling over large parts of the Northern Hemisphere and a warming tendency prevailing since 1990. The recent warming rates are characterized by an equator-to-pole gradient and an ocean-to-land-contrast as reported by several previous studies (e.g. Broccoli et al. 1998;Jain et al. 1999;Boer 2011;SHi and BueH 2012). NCEP also reveals some regional and seasonal details of observed temperature changes, e.g. the cooling of Antarctica in austral summer that is probably related to enhanced cyclone activity in the Weddell Sea, a stronger southern annular mode and ozone dynamics (PaetH and Ocean (cf. SimmonS et al. 2004;IPCC 2013;JoneS et al. 2013). In contrast, NCEP appears to produce spurious cooling trends in the Andes. The leading predictor in the cooling and warming period is downward long-wave radiation at the Earth's surface. Although this predictor is also named the atmospheric greenhouse effect, it is much more complex than the time series of monotonously increasing greenhouse-gas concentrations. It may be affected by enhanced radiative forcing, but also includes feedbacks via surface energy fluxes, atmospheric water vapour content and cloudiness as well as signals from natural forcings (solar and volcanic activity) and internal variability (Wallace and HoBBS 2006;Hantel and HaimBerger 2016). Thus, a possible mechanism is that increasing greenhouse gas concentrations induce stronger sensible and latent heat fluxes into the atmosphere, higher water vapour content as the most efficient natural greenhouse gas (cf. FloHn et al. 1992), deeper convection and enhanced cloudiness at greenhouse-sensitive atmospheric levels. This mechanism is obviously modulated by the land-sea contrast, topography, land surface conditions, climate seasonality and latituderelated conditions such as solar incidence and tropopause height. In addition, warming in high latitudes and high altitudes is strongly affected by reduced surface albedo, indicating an ice-albedo feedback in regions with snow and ice cover (cf. Jain et al 1999). In the low latitudes, especially in monsoon regions and along the ITCZ, changes in wind divergence and related energy fluxes play a noticeable role. In summary, our analysis gives support to the argumentation by Boer (2011), suggesting that complex timedependent local feedbacks are more relevant to the patterns of observed temperature changes than the static contrast of thermal inertia between ocean and land surfaces (Broccoli et al. 1998).
Advection processes occupy a medium rank among all considered predictors. Changes in the advection term have been found to be quite noisy at the grid box scale. Therefore, the large-scale effect of advection may be better represented by index time series of circulation phenomena such as the North Atlantic Oscillation or the northern and southern annular modes that govern these processes at the climatological time scale (rautHe and PaetH 2004; PaetH and Pollinger 2010; IPCC 2013).
A linear combination of up to 11 predictors achieves around 90 % of total variance in the spatial modelling approach and almost 100 % in the temporal approach, even when using de-trended time series. The fact that the temporal approach outperforms the spatial one may be interpreted as an indicator for the spatial heterogeneity of the relevant mechanisms. However, it may also emanate from the very different ratios of number of predictors to sample size (10 to 10,224 for the spatial model versus 11 to 30 in the temporal approach). In some high-mountain areas, along sea-ice margins and in oceanic upwelling regions, the residual part of temperature variability not accounted for by our statistical model amounts to 10-40 %. In these regions, additional mechanisms, subgrid-scale and/or nonlinear processes play a more important role (cf. Boer 2011). While the trend patterns of near-surface temperature clearly differ from period to period and season to season, the underlying mechanisms, i.e. the predictor ranking, are virtually the same.
The amplitude and spatial structure of the observed mid-20 th century cooling phase is barely reproduced by the MPI-ESM Earth system model although it is state-of-the-art in terms of the complexity of the simulated climate system and the imposed natural and anthropogenic forcings. It has been argued that some of the 20 th -century climate variations are partly related to internal variability and, hence, cannot be replicated by climate model simulations with randomized initial conditions (korHonen et al. 2010). In contrast, the warming pattern since 1990 is very coherent between climate model, reanalysis and CRU station data, revealing the well-known equatorto-pole and ocean-to-land gradients (cf. Jain et al. 1999) and even including the slight cooling tendency in the North Atlantic (driJtHout et al. 2012). van oldenBorgH et al. (2013) have reported on the ability of CMIP5 models to simulate the observed patterns of recent warming when all natural and anthropogenic climate drivers are correctly implemented. JoneS et al. (2013) have defined this a major improvement of CMIP5 compared with CMIP3 models. However, they also stated that CMIP5 models do not adequately reproduce the observed cooling patches in the tropical Pacific and Southern Ocean. To a certain extent, this is also found here when comparing MPI-ESM and NCEP. It is likely due to missing processes in the climate model such as changes in seasalt aerosols (korHonen et al. 2010) and melting ice shelfs (BintanJa et al. 2013).
Until the end of the 21 st century, the warming rates steadily increase, culminating in 0.71°C per decade over the landmasses. knutSon et al. (2013) have demonstrated that the CMIP5 multi-model ensemble mean exhibits weaker warming rates than observed. For MPI-ESM the opposite is true, implying that this model has an above-average climate sensitivity in the CMIP5 framework. The ocean-to-land contrast in terms of the spatial-mean temperature rise per decade enhances as well from 0.07°C in the mid-20 th century to 0.27°C after 2071. Sutton et al. (2007) have drawn the same conclusion from CMIP3 models. The NCEP reanalysis foreshadows this tendency when comparing the warming rates over the 1948-2019 period. The intensification of the equatorto-pole gradient and ocean-to-land contrast in the global warming pattern will affect crucial phenomena in the climate system that arise from horizontal temperature gradients, i.e. monsoons, land-sea breeze systems, jet streams, atmospheric waves and weather variability (e.g. Jain et al. 1999;Sutton et al. 2007; PaetH and Pollinger 2019).
The statistical model has identified virtually the same underlying mechanisms in present-day and future warming patterns. In addition, the ranking of predictors is more or less identical in climate models that do not reproduce the observed warming pattern since 1990. Thus, this model failure is not caused by different thermodynamic principles and involved feedbacks relevant to temperature simulation, but by different forcings and/or different climate sensitivity to the imposed forcings.

Conclusions
Returning to the first scientific hypothesis formulated in section 1, it can be concluded that none of the selected predictors alone can account for the spatiotemporal heterogeneity of observed and simulated temperature changes. At every grid box, temperature variations arise from a combination of several mechanisms, mainly related to the surface radiation balance. In addition, these mechanisms vary from region to region, particularly between low and high latitudes. The limited number of predictor fields and the assumed linearity of the statistical transfer functions do not represent major constraints to the explanation of the heterogeneous warming pattern.
Climate models are indeed able to reproduce the observed warming pattern with high accuracy, but only if they are state-of-the-art with respect to the complexity of the simulated climate system and the imposed natural and anthropogenic forcings. Yet, the second hypothesis must be rejected: the underlying processes of global warming barely differ between past and future periods. The same predictor fields are responsible for the mid-20 th century cooling phase, for the recently observed warming pattern and for future surface heating under the assumption of a business-as-usual emission scenario. Strictly speaking, this is not surprising because the selected predictors are more or less complete regarding the theory of temperature changes (see subsection 2.2). However, it is an unexpected result that the predictor ranking does not noticeably change between different time periods, except for an increasing impact of monotonously enhanced greenhouse gas conditions over the 21 st century. In fact, downward long-wave radiation is the main driver for the 1948-1977 cooling and, at the same time, for the subsequent warming era between 1990 and 2100. Thus, the processes associated with the atmospheric greenhouse effect, i.e. changes in evapotranspiration, water vapour content, cloudiness and greenhouse gas concentrations, are more responsible for near-surface temperature variations in most regions of the globe than changes in surface albedo, solar radiation, energy flux divergence or temperature advection.