Detection of non‐climatic biases in land surface temperature records by comparing climatic data and their model simulations

Abstract

The 0.6 °C warming observed in global temperature datasets from 1940 to 1960 to 2000–2020 can be partially due to urban heat island (UHI) and other non-climatic biases in the underlying data, although several previous studies have argued to the contrary. Here we identify land regions where such biases could be present by locally evaluating their diurnal temperature range (DTR = T_Max − T_Min trends between the decades 1945–1954 and 2005–2014 and between the decades 1951–1960 and 1991–2000 versus their synthetic hindcasts produced by the CMIP5 models. Vast regions of Asia (in particular Russia and China) and North America, a significant part of Europe, part of Oceania, and relatively small parts of South America (in particular Colombia and Venezuela) and Africa show DTR reductions up to 0.5–1.5 °C larger than the hindcasted ones, mostly where fast urbanization has occurred, such as in central-east China. Besides, it is found: (1) from May to October, T_Max globally warmed 40% less than the hindcast; (2) in Greenland, which appears nearly free of any non-climatic contamination, T_Mean warmed about 50% less than the hindcast; (3) the world macro-regions with, on average, the lowest DTR reductions and with low urbanization (60S-30N:120 W–90 E and 60 S–10 N:90 E–180 E: Central and South America, Africa, and Oceania) warmed about 20–30% less than the models’ hindcast. Yet, the world macro-region with, on average, the largest DTR reductions and with high urbanization (30 N–80 N:180 W–180 E: most of North America, Europe, and Central Asia) warmed just a little bit more (5%) than the hindcast, which indicates that the models well agree only with potentially problematic temperature records. Indeed, also tree-based proxy temperature reconstructions covering the 30°N–70°N land area produce significantly less warming than the correspondent instrumentally-based temperature record since 1980. Finally, we compare land and sea surface temperature data versus their CMIP5 simulations and find that 25–45% of the 1 °C land warming from 1940–1960 to 2000–2020 could be due to non-climatic biases. By merging the sea surface temperature record (assumed to be correct) and an adjusted land temperature record based on the model prediction, the global warming during the same period is found to be 15–25% lower than reported. The corrected warming is compatible with that shown by the satellite UAH MSU v6.0 low troposphere global temperature record since 1979. Implications for climate model evaluation and future global warming estimates are briefly addressed.

1 Introduction

There has been considerable debate in the literature over the extent (if any) to which non-climatic biases – in particular, those due to urbanization elements such as the urban heat, aerosol emissions, and other factors—have contaminated regional, hemispheric and global temperature records, and artificially raised, or in any way altered global warming trends since the 1940s (e.g.: Freitas et al. 2013; Hansen et al. 2010; McKitrick and Michaels 2007; Menne et al. 2018; Pielke et al. 2007a, b; Scafetta and Ouyang 2019; Soon et al. 2015, 2018; and many others).

The issue is of great concern because from 1950 to 2020 the world population increased from 2.5 billion to 7.5 billion (Population Division of the Department of Economic and Social Affairs of the United Nations: https://population.un.org/wup/). Urbanization did not develop uniformly. In fact, Fig. 1a shows that the world population density is not uniformly distributed on the globe, and Fig. 1b shows how the world city population has increased reporting data for the year 1950, 1990, 2015 and the projected values for 2030 (https://population.un.org/wup/). For example, from 1950 to 2015 in China the urban population increased 12 times, from 65 million to 775 million whereas in Europe it just doubled by increasing from 284 million (1950) to 547 million (2015).

Fig. 1

a World population density in 2015. (GPWv4 Population Density, v4.11_2015; https://sedac.ciesin.columbia.edu/data/collection/gpw-v4). b World city population in 1950, 1990, 2015 and the one projected in 2030 (https://population.un.org/wup/)

Urban areas are usually warmer than the surrounding rural ones (Mitchell 1953, 1961; Landsberg 1981). This is commonly known as the urban heat island (UHI) effect (Oke 1987; Stull 1988; Kershaw 2017). Estimating its relevance is still very challenging.

For example, using a few historical stations located inside and outside London, some authors found that the city’s UHI could be around 1.1 °C (Jones et al. 2009); Wolf and Mc Gregor (2013) found a UHI of 2.8 °C in the summer temperature from 1990 to 2006; Wilby et al. (2011) found that UHI intensity depends on weather conditions. The above values vary from each other because different periods were used, and few and different historical meteorological stations located in places qualitatively labeled as “urban” or “rural” were compared. More extended and detailed studies of London’s thermal distribution using satellite temperature images and hourly air temperature data from 77 stations optimally displaced in the city area found that, in specific hot weather and low wind conditions, the air temperature of the city’s center can also become 6–9 °C warmer than the surrounding rural area (Kershaw 2017; Kolokotroni et al. 2006; Kolokotroni and Giridharan 2008; Watkins et al. 2002a, b; Wolf and McGregor 2013; Holderness et al. 2013).

UHI generally occurs because of the direct release of urban anthropogenic heat and of the cities’ dense concentrations of materials like asphalt, concrete, and buildings (Rizwan et al. 2008). Urban materials absorb more heat during the daytime and release it more slowly at night than the soil and vegetation that characterize the rural areas. Besides, wind ventilation is lower in urban canyons, which contributes to keeping the cities warm (Dimoudi and Nikolopoulou 2003).

UHIs create a genuine localized climatic effect. Thus, the individual meteorological records are not technically biased; as a matter of fact, they are recording genuine localized warming. However, although in stationary conditions UHI should not change climatic averages and trends, a warming trend occurs in meteorological records when cooler rural landscapes are progressively transformed and included into warmer urban areas (e.g.: Saaroni et al. 2000; Kim and Baik 2004; Gaffin et al. 2008; Emmanuel and Krüger 2012).

Localized warming trends are of concern when these records are used for calculating regional, hemispheric, or global temperature climatic trends. Moreover, urban areas only comprise a few percent of the land area, but urban stations comprise a much larger fraction of the used ones, and this fraction is typically greater for the longest station records. Therefore, urbanization bias is essentially an over-representation sampling problem.

UHI effects have been often studied by attempting to distinguish urban from rural land regions from satellite (Wickham et al. 2013), reanalysis (Compo et al. 2013), and homogeneity algorithms (Menne et al. 2018). More specifically, some studies attempted to determine any UHI bias by comparing urban temperature records against those of the surrounding rural areas (e.g.: IPCC AR5 2013; Hansen et al. 2010; Hausfather et al. 2013; Parker 2006). Stewart and Oke (2012) discussed in detail the challenges in estimating UHI biases by comparing urban and rural stations, which are often all partially affected by urbanization. Another approach tried to separate urban and rural regions using nightlight brightness metadata (Peterson 2003). The main conclusion was that the climatic series obtained after homogenization could be nearly free of UHI bias (e.g., Hausfather et al. 2013; Venema et al. 2012; Mestre et al. 2013; Freitas et al. 2013). Yet, these mathematical approaches are specifically designed for reducing different kinds of non-climatic biases such as step-change biases linked to station moves, changes in instrumentation, etc. Trend biases, such as those expected from urbanization development, are much more challenging to correct.

Menne et al. (2018) acknowledged the problem and adopted a specific algorithm claimed to remove UHI by a series of steps that was supposed to work at least for those sites that were considered to be more affected, e.g. Reno, Phoenix in the USA, and Shenzhen in China.

Yet, all the above methodologies have serious limits and the issue of whether the processed climatic records are still contaminated by UHI or other non-climatic biases remains. In fact, when entire extended regions are homogenized, UHIs can influence the same rural regions used to detect and quantify the bias of the nearby urban centers (e.g.: de Gaetano 2006; Pielke et al. 2007a, b; Soon et al. 2015, 2018).

Essentially, the current temperature homogenization algorithms (e.g.: Menne and Williams 2009; Mestre et al. 2013) effectively use the trends of neighboring stations to determine the climatic trends whenever an apparent non-climatic bias is identified (e.g., from a station move). If the station records for several of these neighbors are affected by urbanization bias, then some of that urbanization bias will be aliased onto the target station’s record, even if the target station is still rural and had originally been unaffected by urbanization bias. For example, such a homogenization process tends to blend the trends of all stations towards the average of the neighbors (D’Aleo 2016; Soon 2018). In general, if the neighboring stations have each been affected by urbanization bias (or other non-climatic biases) to different extents, then the most urbanized stations will have their urbanization bias partially reduced. But, simultaneously, the most rural stations will have warming trends artificially added: that is, the homogenization process tends towards “urban blending”.

To make matters worse, these blended biases can no longer be identified by the standard approaches of comparing rural and urban stations, since the homogenized station records of the rural stations now contain this blended urban bias. Soon et al. (2015, 2018) argued that this has contributed to many of the studies using homogenized temperature datasets failing to identify evidence of major non-climatic biases. The consequence is that land climatic records should show warming trends that can only partially be associated with real climatic changes.

Moreover, any criteria attempting to discriminate some regions as “urban” and others as “rural”, where only the former ones could be affected by UHI, could be unsatisfactory. Indeed, the heat produced by the urban areas can also be transported toward the surrounding areas by wind. For example, in London UHI varies from 9 to 1 °C as wind speed increases from 1 to 7 m/s (Kolokotroni et al. 2006; Kolokotroni and Giridharan 2008), which indicates that wind extract and transport a large amount of the urban heat toward the nearby surrounding rural areas. Therefore, also the so-called “rural” areas can be directly contaminated by the surrounding UHIs. This has likely happened in extended and densely populated regions such as those found in central-east China (cf. Scafetta and Ouyang 2019). Besides, as explained above, homogenization approaches interpolate the temperature data among far stations and contaminate large regions by artificially enlarging the UHIs’ area of influence even for hundreds of kilometers (cf.: de Gaetano 2006; Pielke et al. 2007a, b; Soon et al. 2015, 2018).

Solving the urbanization problem is challenging also because the definitions of “rural” and “urban” often differ substantially between studies, and a “rural” station in one study might be classified as “urban” according to another study (and vice versa) (Stewart 2011; Stewart and Oke 2012). Moreover, UHI can develop in small urban settlements and villages as well (e.g.: Hinkel et al. 2003; Szegedi et al. 2013; Dienst et al. 2018). Thus, there is no guarantee that UHI or other local biases can be efficiently removed from the data using some mathematical algorithm. Yet, if the net magnitude of the biases in the adopted climatic records, which are already processed by homogenization, is greater than currently assumed, this could have important implications for climate change estimation, attributions, and climate model validation.

Uncertainty and underestimation of urbanization biases can artificially increase the warming trends in several Earth’s regions (e.g.: Balling Jr. and Idso 1989; Kato 1996; Wang and Ge 2012). McKitrick and Michaels (2007) claimed that UHI effects could reduce the estimated 1980–2002 global average temperature trend over land by about half. Soon et al. (2018) reviewed the literature stressing the current large uncertainty regarding the surface air temperature warming trends over China during the last century due to non-climatic urbanization biases (cf.: Wang et al. 2001, 2004; Tang and Ren 2005; Tang et al. 2010; Ren et al. 2012, 2017; Cao et al. 2013; Ding et al. 2014; Wang et al. 2014; Soon et al. 2015; Li et al. 2017). Temperature reconstructions from tree-rings show significantly lower warming trends than the corresponding land instrumental temperature averages since the 1970s (Esper et al. 2012, 2018).

There are also certain specific cases – chiefly in arid regions – where the process of urbanization could lead to a net cooling, e.g., through increased shading caused by buildings and the presence of more water: oases are obvious examples (cf.: Stohlgren et al. 1998; Rasul et al. 2015). Non-climatic biases are also controlled by variations in the capacity of urban and rural areas to evaporate water (Li et al. 2019), by local microclimates (Fall et al. 2011) and by changes in instrumentation (Hubbard and Lin 2006) and in the observation time (Karl et al. 1986). Finally, non-climatic biases in temperature records could be induced by alternative factors including land-use change, aerosols, and others (cf.: Robinson et al. 1995; Gallo et al. 1996; Dai et al. 1999; Braganza et al. 2004; Makowski et al. 2008; Lim et al. 2005).

The purpose of the current work is to propose a methodology aimed to identify the regions where possible urbanization or other unresolved local non-climatic biases could be still present in the available climatic temperature records used to evaluate global warming trends.

We propose a workaround methodology by extending the analysis proposed in Scafetta and Ouyang (2019) for China. The idea is that non-climatic biases could be localized using some property of the atmospheric boundary layer physics. This physics predicts that UHI effects are not constant in time but vary according to weather conditions, season, and, above all, the time of day.

The strategy is to compare various temperature estimates (T_Mean, T_Max, T_Min) and their climate model simulations, and to process the found divergences at global, macro, and local scales to empirically quantify how much the land and the global surface temperature warming could be fictitious because of non-climatic biases affecting the records.

2 Physical background

Although temperature records could be affected by several non-climatic biases, in this section we focus mainly on possible urban contaminations such as those due to UHI and aerosol emission.

UHI is generally most pronounced at nighttime than during daytime; that is, urbanization makes the minimum temperature (T_Min) to warm more than the maximum one (T_Max) (cf.: Oke 1987; Stull 1988; Kershaw 2017). The difference between T_Max and T_Min is referred to as the diurnal temperature range (DTR). A larger warming trend at nighttime than during daytime leads to a DTR negative trend in time.

There is a long history of studies on DTR changes. One of the earliest papers looking at this was Karl et al. (1993) and others followed (Easterling et al. 1997; Vose et al. 2004; Thorne et al. 2016a, b; and others). These papers also discussed some ideas as to why T_Min could warm more than T_Max in specific situations. Figure 2 schematically summarizes the physics involved in the process.

Fig. 2

(Adapted from Kershaw 2017; from Scafetta and Ouyang 2019)

Boundary layer structure over a city and its surrounding area. Top, daytime boundary layers; middle, nocturnal boundary layers; bottom, microscale structures showing how the wind extracts heat from an urban area and transports it away from it toward the rural areas. Planetary boundary layer (PBL); Urban boundary layer (UBL); Urban canopy layer (UCL).

The atmospheric boundary layer physics predicts that during daytime the surface, including the urban areas, gets warmed by sunlight, which activates a convective instability resulting in a mostly vertical flux of warm air parcels. In other words, during the daytime, the atmospheric capping inversion layer rises and the warm air produced by the city more easily diffuses vertically than horizontally, which reduced the possibility that the surrounding rural areas could get warmed by it. Essentially, vertical air-convection processes favor heat dissipation.

On the contrary, during the nighttime, the UHI area of influence spreads around becoming larger because the surface cools fast forming a low capping inversion layer that reduces vertical air convection flux. Consequently, nocturnal surface winds can more easily diffuse a city’s warm air into the surrounding region and warm it. Rural and suburban areas are, therefore, more exposed to the heat produced by their surrounding UHIs during nights and, consequently, cool less than what would be expected. Thus, urbanization can make warming trends of station placed nearby urban centers anomalously higher during the night than during the day. Also, McNider et al. (2012) showed that a long-term increase in minimum temperatures could derive from a redistribution of heat by changes in wind and air turbulence and not just by an accumulation of heat in the boundary layer. Consequently, these authors concluded that minimum temperatures could be a misleading global warming metric, which implies that also mean temperatures could be misleading.

DTR biases could emerge also for other reasons such as the building-up of sulfate or other aerosols from fossil fuel combustion within the cities in the lower troposphere, which should reduce daytime warming since the early 1950s. For example, Thorne et al. (2016b) showed time series (for global land, North America, Europe, and Australia) with significant year-to-year variability: from 1900 to 1950 DTR values were relatively constant, then they decrease from 1950 to the 1980s, and then remained at their newer level since. The proposed reason for the drop was that power station emissions of sulfate aerosols were reduced leading to a cleaner lower troposphere. However, Zdunkowski et al. (1976) also claimed that aerosols should cool both T_Max and T_Min. In fact, aerosols reduce the amount of solar radiation reaching the surface and this makes a city cooler during the daytime. In such a situation, a city should be cooler also during nighttime as its warmth is mostly stored during the sunny hours. UHI intensity is also controlled by variations in the capacity of urban and rural areas to evaporate water (Li et al. 2019) and by local microclimates (Fall et al. 2011). DTR could be also influenced by changes in instrumentation too (Hubbard and Lin 2006) and in the observation time (Karl et al. 1986), but these biases are supposed to be properly handled by the homogenization algorithms.

In any case, Ren and Zhou (2014) and Jiang et al. (2020) have confirmed that there is a significant asymmetric warming rate between T_min and T_max around most of the meteorological stations in China that experienced rapid urbanization, and a similar situation could have occurred in similar conditions everywhere in the world.

3 Data

We use the following data:

• 1. The Climate Research Unit (CRU) (University of East Anglia) T_Max and T_Min global land surface temperature datasets (CRU-TS4.04) (Harris et al. 2014, 2020; Jones et al. 2012) (https://crudata.uea.ac.uk/cru/data/hrg). These records cover the period from 1901 to 1919 with a spatial resolution of up to 0.5° × 0.5°. The same records are also available at Climate Change Atlas at KNMI Climate Explorer (https://climexp.knmi.nl/plot_atlas_form.py). The regions with missing or too few temperature data are about 5% of the world’s land surface. These places are left in white on the maps and include most of India and Pakistan, most of the Amazonian region in Brazil, and extended parts of Africa. In these regions the T_Min and T_Max records are missing or, as it happens for example for India, they are present but are characterized by identical temperature anomalies between 1945 and 1954 and 2005–2014 as if they were derived from mean temperature values so that DTR values are zero. We also use the HadCRUT4.6 record for the land and ocean surface temperature and for the global surface temperature (Morice et al. 2012), and the UAH MSU v6.0 low troposphere satellite global temperature record (Spencer et al. 2017).

• 2. The Coupled Model Intercomparison Project Phase 5 (CMIP5: https://cmip.llnl.gov/cmip5/) (full set) GCM mean simulations for the T_max and T_min temperature at the surface available at Climate Change Atlas at KNMI Climate Explorer (https://climexp.knmi.nl/plot_atlas_form.py) with a spatial resolution of 2.5° × 2.5°. The records were interpolated to fit a 0.5° × 0.5° grid. The chosen computer simulations adopt the full known set of historical radiative forcings from 1860 to 2005. From 2006 to 2020 the simulations available for alternative RCP (2.6, 4.5, 6.0, 8.5) are statistically equivalent (cf.: Bindoff et al. 2013).

Herein, the ensemble mean rather than the envelope of the full range of models or individual model runs is used because the latter appear too random at the required decadal scale. This random “internal variability” tends to cancel out in the ensemble mean (cf.: Scafetta 2013; Connolly et al. 2019). The CMIP5 ensemble mean could be considered as an ideal and optimized syntetic climate record that combines the performance of all models. Its statistical error is supposed to be very small, and it is herein ignored, because it would be given by the entire dispersion among the models’ single runs divided by the squared root of the number of all used simulations, which are more than 100.

4 Method

The proposed analysis is based on the study of the negative change of the DTR (− DTR = T_Min −  T_Max) metric over a 60-year period. This measure positively correlates with the warming bias that UHI would usually induce and, therefore, it could be more intuitive for representing UHI warming biases.

Our methodological approach aims to determine whether some urbanization or other non-climatic biases are still present in the climatic surface land records obtained after homogenization of the individual station records. The assumption to be tested is whether the climatic processed global surface land temperature databases are free or at least little affected by local thermal contaminations such as UHI. To check the situation we determine the warming divergence observed between nocturnal and diurnal hours between the decades 1945–1954 and 2005–2014 in each region where the data are available by analyzing and comparing the available T_Min and T_Max local climatic records throughout the world.

The physical complexity of the process makes very challenging any interpretation of the data. We propose to compare the estimated DTR changes against their predictions made by the GCM ensemble mean simulations (Bindoff et al. 2013). The reason why such a comparison could be relevant is that the DTR changes are known to have multiple possible causes (e.g.: cloud cover, wind, urban heat, land-use change, aerosols, water vapor, and greenhouse gases) and different regions can be affected by different factors (cf.: Robinson et al. 1995; Gallo et al. 1996; Dai et al. 1999; Braganza et al. 2004; Makowski et al. 2008; Lim et al. 2005). These climate models are supposed, at least theoretically, to properly simulate the global and local climatic variations induced by cloud cover, land-use change, aerosols, water vapor, greenhouse gases, etc. A good correlation has also been found between the DTR reduction from 1940 to 2014 and the precipitation series for the northern hemisphere (Sun et al. 2019), but once again the models are supposed, at least in principle, to reproduce precipitation trends on such a long time scale. Thus, the climate models should, at least teoretically, predict local and global changes in DTR due to real climatic changes, at least on sufficiently long periods such as 60 years. However, the climate models do not parameterize the cities so that significant differences between the observed DTR climatic changes and those predicted by the models could highlight possible urbanization biases, particularly in highly densely populated regions. In any case, other non-climatic biases could be involved as well, so that our results could have a wider interpretation as well.

For each of the four analyzed records and every grid cell, the T_Min and T_Max mean values for the decades 1945–1954 and 2005–2014 are evaluated together with the respective temperature variation ΔT_Max and ΔT_Min between the same decades: this gives the warming observed in T_Max and T_Min in each place from the decadal averages centered on 1950 and 2010.

The above 60-year period was chosen for the following reasons: (1) it covers the post World War II period when the world has experienced a large population and urbanization increase, and anthropogenic global warming; (2) it avoids possible statistical biases due to the likely presence of a quasi 60-year natural climatic cycle that is detected across the world in several climatic records (e.g.: Gervais 2016; Scafetta 2013, 2014; Scafetta et al. 2020; Wyatt and Curry 2014; and many others); (3) it avoids other possible statistical biases due to the spatial inhomogeneity of the climatic effects of the strong El-Nino warming event that occurred between 2015 and 2016 (e.g. Scafetta et al. 2017a, b).

The proposed DTR analysis uses monthly averages since January 1945 when the CRU-TS4 record covers most of the land surface. In any case, only the 0.5° × 0.5° cells that contain at least 60 months of data for each decade (e.g. 1945–1954 and 2005–2014) are taken into account to avoid biases due to data deficiency.

The T_Min to T_Max warming divergence from 1945 to 1954 to 2005–2014 is evaluated for each grid cell using the equations:

$${\text{Div}}_{{{\text{data}}}} = (\Delta {\text{T}}_{{{\text{Min}}}} {\text{ }} - {\text{ }}\Delta {\text{T}}_{{{\text{Max}}}} {\text{)}_{\text{data}}}$$

(1)

$${\text{Div}}_{{{\text{model}}}} {\text{ = (}}\Delta {\text{T}}_{{{\text{Min}}}} {\text{ }} - {\text{ }}\Delta {\text{T}}_{{{\text{Max}}}} {\text{)}_{\text{model}}}$$

(2)

Finally, the difference between Eqs. 1 and 2 is calculated for each grid cell. This gives the desired local temperature bias estimate:

$${\text{Bias = Div}}_{{{\text{data}}}} {\text{ }} - {\text{ Div}}_{{{\text{model}}}}$$

(3)

This operation is repeated for the CRU-TS4 data and the CMIP5 GCMs ensemble mean simulations, respectively. The analysis is then repeated also for the period between the decades 1951–1960 and 1991–2000, and it is repeated for both data and the CMIP5 ensemble means.

Finally, additional analyses are proposed by comparing mean instrumental temperature records for land and sea surface, and temperature reconstructions from tree-rings to attempt to evaluate the bias induced in by non-climatic mechanisms.

5 DTR trend analysis

Figure 3a, b shows the T_Min and T_Max global surface records for both the CRU-TS4 data and the CMIP5 GCM ensemble mean simulations overland alone. The records are anomalies relative to the decade 1945–1954. The bottom panels (Fig. 3c, d) show the T_Min minus T_Max, namely the negative DTR records, in both cases.

Fig. 3

T_Min and T_Max global surface records for both the CRU-TS4 data (a) and the CMIP5 GCM ensemble mean simulations overland alone (b) and (bottom) their differences (c and d)

From the decades 1945–1954 to 2005–2014 in the CRU-TS4 data T_Min warmed by about 1.0 °C while T_Max warmed by about 0.75 °C. On the contrary, in the computer simulation T_Min warmed by about 1.05 °C while T_Max warmed by about 0.95 °C. In both cases, T_Min warmed more than T_Max. However, regarding the temperature data, from the decade 1945–1954 to 2005–2014 T_Min warmed by 0.25 °C more than T_Max, while for the computer simulations T_Min warmed by 0.10 °C more than T_Max. The data also show that most of such 0.25 °C warming occurred between 1950 and 1990, while during the same period the computer simulations show that T_Min warmed just about 0.05 °C more than T_Max. This anomaly could have been caused because from 1950 to 1990 a significant industrialization and urbanization occurred in vast regions of the Earth such as in Europe, in the USA, and other places. From 1900 to 1950 T_Min and T_Max trends are nearly equivalent.

Figure 4 repeats the same analysis shown in Fig. 3 but now the Earth’s surface is divided into three regions: a northern region from latitude 25° N to 90° N; an intertropical region with latitude between 25° S to 25° N; and a southern region with latitude between 90°S and 25°S. The results are also qualitatively similar. In all regions and from 1950 to 1990, a strong divergence of T_Min relative to T_Max is observed: Figure 4a, b, and it is larger for the northern and southern regions (about 0.35 °C) than for the intertropical region (about 1.5 °C). However, after 2000 the southern region shows nearly compatible T_Min and T_Max values. On the contrary, the CMIP5 GCM ensemble mean simulations show T_Min values slightly larger than the T_Max ones: Figure 4c, d. From 1950 to 1990 the northern and intertropical regions are similar (the divergence is about 0.1 °C) while for the southern region the divergence is about 0.05 °C.

Fig. 4

As in Fig. 3 but the Earth’s surface has been divided into three regions: (left) latitude 25° N to 90° N; (center) latitude between 25° S to 25° N; (right) latitude between 90°S and 25°S. a Temperature data; b DTR of the data; c Synthetic data; d DTR of the synthetic data

Table 1 collects the statistics derived from the records depicted in Figs. 3 and 4 for the decades 1945–1954 to 2005–2014.

Table 1 Temperature variation ΔT_Min and ΔT_Max, and their difference between the decades 1945–1954 and 2005–2014 for the CRU-TS4 data and the CMIP5 GCM ensemble mean simulations. See Figs. 3 and 4

Figure 5a, b depict the world distribution of the divergence between T_Min and T_Max between 1945 and 1954 and 2005–2014, and between 1951 and 1960 and 1991–2000 (using Eqs. 1 and 2) using the CRU-TS4 data. Figure 6a, b do the same using the CMIP5 GCM ensemble mean simulation. The two figures show quite different patterns.

Fig. 5

World distribution of the divergence between T_Min and T_Max between a 1945–1954 and 2005–2014 and b 1951–1960 and 1991–2000 using Eqs.1 and 2

Fig. 6

Same as Fig. 5 using the CMIP5 GCM ensemble mean simulation a 1945–1954 and 2005–2014 and b 1951–1960 and 1991–2000 using Eqs. 1 and 2

Both figures show that T_Min has usually warmed more than T_Max in most regions. However, Fig. 5 highlights that the temperature records are characterized by a strong regional variability spanning between − 2 and 2 °C. On the contrary, Fig. 6 shows that the CMIP5 GCM ensemble mean temperature simulations produce, as expected, nearly spatially homogeneous results. These are characterized by a much less regional variability spanning between − 0.3 and 0.5 °C.

Figure 5 shows that vast regions of Asia (in particular Russia and China) and North America, a significant part of Europe and Australia, and part of South America (in particular Colombia and Venezuela) are characterized by a T_Min-T_Max warming divergence up to 1–2 °C (red-dark red color). Africa too shows a significant variability but overall there is a balance between positive (red) and negative (blue) values. Across the world, the only negative values (blue) are found in Algeria, between Niger and Chad, north-east Ethiopia, west Yemen, and some regions of Kazakhstan and South Australia. These places are mostly arid and in such environments urbanization may have an inverted effect by generating cool islands (e.g. oases) instead of warm ones (cf. Rasul et al. 2015). Central Bolivia too is blue, but there the inverted effect could have been induced by the intense deforestation of the place which increased the local albedo. This extended land transformation is evident from Google Earth satellite photographs (cf. Killeen et al. 2008; Müller et al. 2012).

Over the ocean, Fig. 5a shows a white color because the CRU-TS4 T_Min and T_Max data cover only the land. The same happens for those land regions with missing data. The data coverage improves in Fig. 5b. Investigating this issue goes beyond the purpose of the present work.

Figure 5 shows clear patterns whose spatial resolution may also depend on the land distribution of the weather stations used to build the CRU-TS4 database, which does not cover uniformly the entire world surface and varied in time (Harris et al. 2020). However, this sampling problem is herein ignored because the issue here is to study the temperature coverage as contained in the CRU-TS4 database, not the performance of single meteorological stations.

Figure 6 shows that the climate models show that T_Min and T_Max warmed in the same way (nearly white color) over the oceans. However, from 1945 to 1954 and 2005–2014, in the polar regions, T_Min warmed slightly more than T_Max up to 0.5 °C (from coral to light violet color) as in most land regions. For most of Europe, Turkey, north Arabia, south of United States of America, Mexico, Brazil, south Argentina, south Chile and South Africa, the CMIP5 model ensemble mean predicts that T_Max had to warm just a little bit more (up to 0.2 °C) than T_Min (yellow color).

Using single model runs, a larger regional DTR variability during the same periods can be found. However, this variability does not correlate with that observed in the data. In fact, it appears quite random, as indirectly demonstrated by the very low regional DTR variability found in the mean ensemble simulations shown in Fig. 6, which indicates that the various models show very different local patterns that contradict each other while they are supposed to reproduce climatic patterns on a 60-year period also locally.

Figure 7 shows the difference between the results depicted in Figs. 5 and 6 using Eq. 3. The figure highlights the spatial distribution of the discrepancy between the results obtained using the T_Min and T_Max records and the synthetic ones obtained with the CMIP5 GCMs. The results are qualitatively similar to those commented above in Fig. 5.

Fig. 7

World distribution of the discrepancy between the CRU-TS4 temperature data (Fig. 5) and the CMIP5 ensemble mean simulation (Fig. 6) results using Eq. 3. a 1945–1954 and 2005–2014 and b 1951–1960 and 1991–2000 using Eq. 3

Figure 7 reveals that extended land regions at 30° N–50° N and 65° N–75° N (mostly Siberia) are characterized by a strong positive bias that is on average up to 0.5 °C. In general, the 20°-80° N region, which is also the most populated one, has a strong positive bias (between 0.0 and 0.5 °C). The 20° S–20° N region has a moderate positive bias (between 0.0 and 0.2 °C). Finally, the 20° S–60° S region is more balanced (between − 0.3 and 0.3 °C). Data from Antarctica are missing.

Figure 8 presents a seasonal analysis of the records. It shows T_Min and T_Max global surface records for both the CRU-TS4 data (left) and the CMIP5 GCM ensemble mean simulations overland alone (right) divided into summer records (from May to October) and winter records (from November to April). The first raw refers to the entire world, whereas the second, third, and fourth raw refers to the latitude bands shown in Fig. 1a: 25° N–90° N, 25° S–25° N, and 90° S–25° S.

Fig. 8

T_Min and T_Max global surface records for both the CRU-TS4 data (left) and the CMIP5 GCM ensemble mean simulations overland alone (right) divided into summer records (from May to October) and winter records (from November to April). First raw, world data; second, third, and fourth raw, 25° N–90° N, 25° S–25° N and 90° S–25° S latitude regions

Table 2 reports the temperature variation ΔT_Min and ΔT_Max and their difference between the decades 2005–2014 and 1945–1954 of the records depicted in Fig. 8.

Table 2 Temperature variation ΔT_Min and ΔT_Max and their difference between the decades 2005–2014 and 1945–1954 of the records depicted in Fig. 8

The result confirms again that there is a significant divergence between data and simulations which is stronger in the northern regions (25° N–90° N) and it is more stressed during the cold months from November to April. In general T_Max during the warm months from May to October not only warmed significantly less than T_Min but also warmed significantly less than the model predictions. For example, for the northern region 25° N–90° N, T_Max warmed 0.34 ± 0.06 °C less than T_Min, while the models’ prediction was 0.09 ± 0.08 °C. Moreover, during the warm months, T_Max warmed about 0.41 ± 0.07 °C less than the CMIP5 model prediction. Similar results were also found for China (Scafetta and Ouyang 2019).

These seasonal results are consistent with a UHI contamination since such biases would be stressed during the cold months more than during the warm ones because, according to the atmospheric boundary layer physics, during the warm season atmospheric convective instability and vertical flux are larger than during the cold months.

6 Visual evidence for UHI biases in macro‐regions

Figure 9 enlarges four macro-regions depicted in Fig. 5a. The depicted areas include America, central and north Africa, Europe, and south–east Asia, and Australia. The maps also report the position of the cities with more than 50,000 people (cyan dots) and, in the case of the bottom-right panel regarding Australia, we added the cities with more than 10,000 people (black dots) (Free World City Database, worldcitiespop.txt 2012). The figure confirms that the most populated regions, which are characterized by a high concentration of urban centers, are often also the ones with the largest DTR negative trend.

Fig. 9

Zooms of Fig. 5a. The cyan dots are cities with more than 50,000 people, the black dots in the bottom-right panel are cities with more than 10,000 people

Figure 9a shows reddish color areas in correspondence with the highly populated central-east and west USA areas, south Alaska around its largest city Anchorage, north Columbia, and northern Venezuela.

Figure 9b shows a reddish area in the highly populated region from Italy up to the United Kingdom and of south-east Europe from Greece to Russia compared to the less populated areas of Spain, Ireland, and north-east Europe which are mostly yellow. Africa and West Asia show reddish areas from western Turkey to Israel and Syria, in the eastern regions of the Arabia peninsula where, on the coast, large cities are present such as Manama, Dammam, Dubai. The same is noted in the vast highly populated region of Nigeria, and of Egypt and North Sudan along with the populated areas where the Nile river runs (cf. Fig. 1). Also, the blue desert regions in Algeria, east Ethiopia, and south-west Arabia seem influenced by the presence of some cities.

Figure 9c highlights the reddish color around the highly populated area of south-east Brazil, northern Argentina, and Chile. The bluish area of Bolivia appears to be due to the deforestation of the area.

Figure 9d shows a reddish color in the highly populated area of Vietnam, Cambodia, the Philippines, and Indonesia. New Zealand is violet while Australia, which is poorly populated and vastly desertic, still presents some reddish areas in the north–east coastal regions (Queensland) while the greenest area of the south–east (Victoria and Canberra) are less affected. A reddish color is found around some central Australia minor cities such as Alice Spring, which is located in a darker area relative to the surrounding desertic one.

In general, the salmon–yellow–green areas of Fig. 9 do not present large urban agglomerates.

Similar visual evidence of a diffused UHI effects captured by the DTR local trends and distribution regarding the most populated regions of China is reported and studied in detail in Scafetta and Ouyang (2019). In any case, again, the reddish areas are usually the most populated ones.

7 Comparisons between T _Mean trends and their model simulations in macro-regions.

If some regions of the Earth are particularly contaminated by UHI, it would affect both T_Max and T_Min trends and, therefore, their mean temperature records should present a warming trend larger than what happened. However, the models are optimized to reconstruct the warming shown in the actual global mean surface temperature record. Thus, if there is a warming bias in the temperature record, the CMIP5 ensemble mean temperature simulations should overestimate the warming where such biases are absent or minimally present, and underestimate it where the UHI biases are expected to be larger. To check whether such an expectation is correct, we run a few tests by comparing mean temperature records and their CMIP5 simulations in specific macro-regions that can be discriminated by urbanization levels.

7.1  Greenland

Figures 5 and 7 suggest that Greenland should be one of the few macro-regions of the world that is mostly free of any non-climatic bias. In fact, in Fig. 7 the area color varies from lite yellow (– 0.2 °C) to lite coral (0.2 °C). Indeed, Greenland has a surface of about 2,200,000 km² (approximately as West Europe) and a total population of just a small European city, that is only 56,000 people (by 2018) that did not grow significantly since 1980 (source: https://en.wikipedia.org/wiki/Greenland). Moreover, only one town has a population larger than 10,000 people (Nuuk, 18,000) and the other four largest towns have a population between 3000 and 6000. Only Southern Greenland, where its largest towns are located, appears modestly affected by a non-climatic bias (slight coral color in the figures). Indeed, UHI effects have been observed also for small arctic towns of comparable size such as Barrow (AK, USA, actual population of about 4500) (Hinkel et al. 2003). Thus, non-climatic anthropogenic biases (such as UHI, sulfate aerosol, and others) should be nearly negligible in Greenland.

The climate record of Greenland was carefully reconstructed since 1840 (Box et al. 2009). This record shows that the Greenland mean temperature in the 1940s was just moderately colder than the post-2000 warming period. The same is true for the entire polar arctic region from 70°–90° N.

Figure 10a compares the CRU-TS4 mean temperature record over Greenland and its corresponding GMIP5 GCM ensemble mean reconstruction. Both records are available at Climate Change Atlas at KNMI Climate Explorer (https://climexp.knmi.nl/plot_atlas_form.py). It is observed that from 1930 to 1950 to 2005–2017, Greenland warmed by about 0.7 ± 0.2 °C. On the contrary, the CMIP5 ensemble mean temperature record predicts a warming of about 1.5 ± 0.1 °C between the same two periods.

Fig. 10

a CRU-TS4 observed and CMIP5 GCM simulated ensemble mean temperature anomaly over Greenland relative to the 1930–1950 mean value. b Locations of towns in Greenland with population > 10,000 (red), > 3000 (orange), > 1000 (green), and > 300 (blue) (https://en.wikipedia.org/wiki/List_of_cities_and_towns_in_Greenland)

Thus, the models (which are calibrated to reproduce, on average, the observed global warming trend since 1900) hindcast almost twice the warming than what has been recorded in Greenland since 1930–1950.

7.2 Three continental macro‐regions

The color patterns of Figs. 5 and 7 suggest that most of the world can be approximately divided into three macro-regions with different mean DTR biases. We analyze the following ones.

• 60 S–30 N:120 W–90 E: it comprehends Central and South America and Africa and appears to be on average less affected by DTR biases and is also modestly populated.

• 60 S–10 N:90 E–180 E: it comprehends South-East Asia and Oceania and appears to be more affected than the previous one, but less affected than the following one.

• 30 N–80 N:180 W–180 E: it comprehends most of North America, Europe, and Central and North Asia. This macro-region is the most urbanized and populated one and also appears to be the most DTR biased one.

For each of these macro-regions, Fig. 11 compares the CRU-TS4 mean records and the CMIP5 mean surface temperature hindcasts on land alone.

Fig. 11

a–c Comparison between the CRU-TS4 records and the CMIP5 mean prediction in the macro-regions 60 S–30 N:120 W–90 E, 60 S–10 N:90 E–180 E, and 30 N–80 N:180 W–180 E. The temperature anomalies are relative to the 1900–1950 period. d–f Residual functions between predictions and data are fit with parabolic curves

The upper panels a–c show the actual mean temperature records, while the bottom panels d– f show the discrepancy between the synthetic records and the reported temperature records. These residual functions are fit with parabolic curves to better highlight systematic biases.

This comparison confirms that from 1900 to 2020 in the first two and less DTR-affected macro-regions (60 S–30 N:120 W–90 E and 60 S–10 N:90 E–180 E), the computer simulations significantly overestimate the observed warming by about 0.3 °C and 0.2 °C, respectively, or about 20–30% of the observed warming of about 1 °C. On the contrary, for the third and most DTR-affected macro-region (30 N–80 N:180 W–180 E) the models underestimate the reported warming by about 0.08 °C, or about 5% of the observed warming of about 1.6 °C.

7.3 Instrumental and tree‐based land temperature reconstructions

A very large number of temperature proxy reconstructions have been proposed in the literature. Regarding the proxy temperature records for the late 20th century a so-called “divergence problem” has emerged (Taubes 1995; Jacoby and D’Arrigo 1995; Briffa et al. 1998; Esper et al. 2012, 2018). This is the observed disagreement between the reported land instrumental temperatures warming of the last decades and the one reconstructed from the latewood densities and tree-rings widths: the latter show much less warming.

For example, Esper et al. (2018) proposed six tree-ring temperature reconstructions for the 30° N–70° N land area which approximately agrees with the third macro-region analyzed above that is suspected to be seriously affected by non-climatic warming biases. Four of the proposed tree-based reconstructions last until 2002 and were compared against the summer (JJA) instrumental temperatures averaged over 30° N–70° N land area since 1880. The summer temperature data were used because the trees are supposed to grow mostly during such a season.

Figure 12 shows the average of Esper’s four proxy temperature records against the corresponding JJA instrumental temperature record over the 30° N–70° N land areas (relative to the 1930–1960 mean period). The two records agree sufficiently well from 1880 to the 1970s, but significantly diverge afterward. Relative to the 1930–1960 period, the instrumental records show that the 1990–2002 period warmed by 0.4 ± 0.1 °C. However, during the same period, the tree-based temperature reconstruction warms by just 0.05 ± 0.05 °C.

Fig. 12

Tree-based mean temperature reconstructions (red) against the JJA instrumental temperatures averaged over 30–70 °N land areas (blue). The curves are scaled over the 1930–1960 period (Esper et al. 2018)

The explanation for the observed strong divergence (0.35 °C from 1980 to 2002, on a recorded warming of 0.4 °C) is still unclear. However, a possibility is that the instrumental temperature records are affected by non-climatic warming biases, such as those due to urbanization, which cannot influence trees that grow in uninhabited large forests which are unambiguously genuine rural regions.

In general, to properly reproduce the temperature trends in Greenland, as well as in other macro-regions or even globally, the models might need to simulate correctly the climatic oscillations including the North Atlantic Oscillation and the amount of sea-ice around Greenland, which they appear to fail for several reasons (cf.: Scafetta 2013; Connolly et al. 2017; Scafetta et al. 2020).

8 A tentative estimate of the non-climatic land and global warming bias

The above results suggest that the land climatic record is affected by significant non-climatic warming biases. In this section, we propose a methodology to estimate the global non-climatic warming bias affecting both the land and the global surface temperature records. It uses a comparison between ocean and land global mean temperature data against their CMIP5 predictions: see Fig. 13. In this test, we use the HadCRUT4.6 land and ocean surface temperature database.

Fig. 13

Comparison between Land and SST records using the HadCRUT4 records and the CMIP5 ensemble mean values. The temperature anomalies are relative to the 1940–1960 period

The rationale is that sea surface temperatures (SST) are expected to be free of urbanization and other non-climatic biases that may be affecting the land climatic records. These records may be affected by other forms of biases, particularly for the early 20th century (e.g.: Kennedy et al. 2011; Kennedy 2014; D’Aleo 2016). Davis et al. (2019) even found that the implied 1950–1975 global SST trends were of opposite sign depending on which type of SST data they used. However, in the following, we assume that the global SST record is sufficiently accurate at least since 1940–1960 and use it as a reference.

The proposed methodology takes into account that SST is expected to warm less than the land because of its larger heat capacity: as the global land temperature and SST records depicted in Fig. 13a show. The same qualitative behavior is reproduced by the CMIP5 GCM ensemble mean predictions because these models take into account the different thermodynamical properties of the two regions (Fig. 13b).

However, when the land and SST records are compared directly against their CMIP5 GCM ensemble simulations, the simulation almost exactly reproduces the land record (Fig. 13c), but significantly overestimates the SST warming (Fig. 13d).

Using the 1940–1960 period as the anomaly reference, in the 2000–2020 period the land temperature was 0.97 ± 0.05 °C against the modeled 1.03 ± 0.04 °C. However, in 2000–2020, the SST temperature was 0.41 ± 0.03 °C against the CMIP5 ensemble mean simulation show a warming of 0.69 ± 0.03 °C. Thus, there is a significant observation versus simulation divergence regarding the ocean region, which covers 70% of the entire world surface.

By assuming the SST warming trend accurate since 1940–1960 the above result once again confirms that the land temperature record significantly exaggerates the warming. In fact, the models likely overestimate the SST warming trend simply because they were calibrated to reproduce the global mean temperature warming of the 20th century, as also discussed and found in Sect. 7.

By assuming the SST warming trend reliable, the above results can be used to calibrate the land temperature record using adjusted model estimates. This operation assumes that (1) the synthetic records could correctly evaluate the different climatic warming expected by the land versus the sea surface and (2) they can be adjusted by modifying their sensitivity to radiative forcings within a large but certain range. The first assumption could be justified by the claim that these models simulate, at least theoretically, all known physical aspects of the phenomenon; the second assumption is supported by the finding that the equilibrium climate sensitivity (ECS) to radiative forcing is still extremely uncertain because ranging from about 1 to 6 °C for a CO₂ doubling (Knutti et al. 2017). Since the mean CMIP5 ECS is about 3 °C (IPCC 2013), the CMIP5 ensemble mean simulation could theoretically be scaled down by even a factor of 3 and still be consistent with the ECS values determined in the scientific literature. The proposed adjustment can be made, in the first approximation, by applying appropriate scaling factors to the temperature reconstructions as explained below.

First, we scale the CMIP5 synthetic ocean record by a factor equal to 0.41/0.69 to make its warming from the twenty-year period 1940–1960 to 2000–2020 compatible with the SST record (which has been assumed to be accurate). The same scaling must then be applied to the synthetic land record because the CMIP5 models process the sea surface temperature together with the land one. With this scaling, in 2000–2020 the modeled land mean warming decreases from 1.03 ± 0.04 °C to 1.03*0.41/0.69 = 0.61 ± 0.06 °C.

At this point, a second scaling equal to 0.61/0.97 = 0.63 is applied to the instrumental land record to make its 2000–2020 warming relative to the 1940–1960 period compatible with that hindcast by the scaled model expectation. This second scaling reduces the modeled land warming from 1940 to 1960 to 2000–2020 from 0.97 ± 0.05 °C to 0.61 ± 0.07 °C.

Thus, under the above assumptions, the 0.97 ± 0.05 °C instrumental land warming observed from 1940 to 1960 to 2000–2020 is likely made of a 0.61 ± 0.07 °C possibly autentic climatic warming plus a 0.36 ± 0.04 °C non-climatic warming bias. This means that 25–45%, that is about a third of the recorded global surface land warming from 1940 to 1960 to 2000–2020 could have been due to urbanization and other unidentified non-climatic biases.

The above results can be used to evaluate a corrected global surface temperature record using a combination of the original SST record and the land record once corrected of the estimated non-climatic bias using the above-calculated scaling. These two records are combined by weighting the land one by 30% and the SST one by 70%, that is: T_corr = 0.3 T_land,corr + 0.7 SST. In fact, using linear regression, we determined that with such weights the original land and SST records reconstruct the HadCRUT mean surface temperature record.

Figure 14a shows the original HadCRUT4.6 record (black), the proposed corrected one (red) versus with the CMIP5 mean ensemble temperature hindcast (yellow) together with 106 independent model runs. Figure 14b shows the same versus the UAH MSU v6.0 low troposphere satellite global temperature record (Spencer et al. 2017) using the same 1979–1985 baseline. The data are anomalies relative to the 1940–1960 period. Relative to the 1940–1960 period, the 2000–2020 period was 0.59 ± 0.03 °C using the original HadCRUT record, 0.44 ± 0.03 °C using the UAH MSU record, 0.48 ± 0.03 °C using the proposed corrected global surface temperature record, and 0.78 ± 0.03 °C using the CMIP5 GCMs ensemble mean record.

Fig. 14

The original HadCRUT4.6 global surface temperature record (black), the proposed one corrected of the land non-climatic warming biases (red). Together with: a the CMIP5 mean ensemble temperature prediction (yellow) with 106 single different model simulations (green); b the UAH MSU v6.0 low troposphere global temperature record (blue) using the same 1979–1985 baseline. The data are anomalies relative to the 1940–1960 period

Thus, according to the above assumptions, the non-climatic biases may have contributed between 15 and 25%, that is about a fifth of the global warning from 1940 to 1960 to 2000–2020 recorded in the official global mean surface temperature record, which is not a negligible error. It is also relevant that the corrected global temperature record is statistically compatible with the satellite UAH MSU v6.0 low troposphere global temperature record.

Moreover, from 1940 to 1960 to 2000–2020, the CMIP5 GCMs may have on average overestimated the climatic warming by about 40%. Even assuming the 106 CMIP5 single model runs, Fig. 14 shows that they appear to produce almost always warmer trends than the proposed corrected global surface record.

9 Discussion and conclusion

The warming observed globally from the twenty-year period 1940–1960 to 2000–2020 (about 0.6 °C) may be partially due to UHI and other non-climatic mechanisms even though several previous studies have argued to the contrary. Herein we have addressed this issue by comparing land and SST climatic local and global records and their CMIP5 model predictions. We have looked at warming trends over periods of about 60 years. This time scale coincides with one of the largest well-known climatic oscillations observed in climatic records for centuries and was chosen to overcome possible statistical and trend artifacts that could affect the analysis by not considering it (Scafetta 2014; Scafetta et al. 2020; Wyatt and Curry 2014).

The DTR trends were studied because UHI area enlarges during nighttime and nocturnal near-surface winds diffuse more urban heat upon suburban and rural areas, warming them. Thus, in presence of urbanization development, DTR is expected to decrease.

In general, land use and management changes such as irrigation, deforestation, overgrazing of grasslands, etc., like urbanization, also produce non-global warming trends (e.g.: Liu et al. 2019; Pielke et al. 2016). In this regard, Zipper et al. (2019) showed that land surface effects can be regionally teleconnected to locations without a landscape change. In specific environments, e.g. in some arid regions, urbanization, which develops in the proximity of oases, could induce an opposite effect because their urban areas should form cold islands. The same could happen when a region is transformed by intense deforestation, as happened in Bolivia, which causes a local albedo change.

Also anthropogenic aerosols (such as sulfate emitted mostly by the combustion of fossil fuel) cools daytime temperature (Cai et al. 2016). However, according to Zdunkowski et al. (1976), aerosols should cool both T_Max and T_Min. Yet, what was found here is that often T_Max warmed less than what the models hindcast while T_Min warmed more than the synthetic records. Moreover, GCMs include aerosol, land-use changes, and other forcings. Thus, if these alternative factors were responsible for the observed DTR changes on global and local scales, the models would be supposed to reproduce the observations, which, according to the above figures, they cannot do. Thus, although a localized aerosol bias may have contributed to the observed patterns, it cannot explain them fully, as already noted for China (Scafetta and Ouyang, 2019; Jiang et al., 2020).

Figures 3 and 4 showed that DTR decreased from 1950 to 2000 (when the urbanization greatly increased in the world significantly) more than what the CMIP5 GCMs hindcast. Yet, in the inter-tropical regions, this bias is less relevant. This result reveals possible UHI contaminations as well as aerosol effects, land use, and/or other non-climatic biases.

Figures 4, 5, 6, 7, 8 and 9 confirm that in vast regions of the world, from the decade 1945–1954 to 2005–2014, DTR values decreased more than the model predictions. Although part of the observed data-model discrepancy could be due to an inability of the GCMs to properly reconstruct climatic patterns in specific environments, UHI, and aerosol effects could still explain a significant portion of the results because many urbanized areas coincide with areas characterized by significant DTR decrease not modeled by the CMIP5 models. The same result was found for China (Scafetta and Ouyang 2019; Jiang et al. 2020). In general, the complexity of these phenomena makes it difficult to correctly interpret what happens in each region. In any case, the results demonstrate that either the data are still significantly affected by local non-climatic biases such as urbanization, or that the CMIP5 models poorly interpret the data patterns at both local and global scales, and must be greatly improved, or both.

It may be possible to wonder whether such a result emerged just because the data were compared directly against the surface temperature ensemble mean model. Yet, the models should be expected to reproduce the DTR trends on 60-year periods at both the local and global time scales. Thus, if the data were free from local non-climatic biases, a failure to reproduce the observed patterns in a consistent way among the various models, would still argue against the latter. In any case, the results suggest that the climatic global land surface temperature is likely contaminated by regional non-climatic biases that the models are not able, or not supposed to simulate for several reasons.

In any case, many regions showing large DTR biases (reddish color in Figs. 5 and 7) are also densely populated and have experienced significant urbanization (cf.: Fig. 1 and 9). These regions include east China, south-east Asia, east and western USA, north Colombia and Venezuela. Even localized regions such as the area around Rio de Janeiro and San Paolo in south–east Brazil are affected by UHI. The urbanization of these areas has on average increased 10 times from 1950 to 2015: see Fig. 1b (https://population.un.org/wup/). A moderate warming bias is found also in Europe (United Kingdom, France, Italy, Greece, and Slavic countries up to Ukraine and south Russia), and central and north-east Australia. North Siberia (60°–80° N, 75°–120° E) is also characterized by a strong hot spot but this region is sparsely inhabited. This bias may be due to a lack of weather stations in the region and to the development of a very localized UHI effect (cf. Konstantinov et al. 2015), or alternatively it may be due to the melting of an increasing mass of permafrost, which could cool daytime temperatures by melting and absorbing latent heat from the environment and warm nighttime temperatures by freezing again and releasing latent heat (cf. Biskaborn et al. 2019), which could have an effect on local instrumental temperature records similar to increasing UHI.

The seasonal analysis shown in Fig. 8 reinforces the interpretation that the land temperature climatic record contains a significant non-climatic UHI contamination component. In fact, the warming observed during the summer months in the highly populated north hemisphere is significantly lower than the model hindcasts. This result may not be fully explained by the addition of atmospheric aerosols such as sulfate since air pollution is lower during the warm season (Lin et al. 2008). On the contrary, it may be explained by the atmospheric boundary layer physics nearby UHIs and by the seasonal changes in wind strength, which in summer are usually stronger.

To further investigate the issue, we analyzed the divergence between global and local mean temperature records relative to the model simulations. In fact, the models are calibrated to reproduce the global warming observed during the 20th century. Thus, if some regions present a non-climatic warming bias, the models are on average expected to underestimate the warming where such a bias is more prominent and overestimate the warming where the non-climatic warming biases are less relevant.

The above expectation was confirmed in several ways. We found that the CMIP5 ensemble model simulation significantly overestimates the warming observed in Greenland and in other macro-regions expected to be little affected by urbanization while underestimating the warming observed in more urbanized macro-region such as in the 30 N–80 N:180 W–180 E macro-region, which includes the most urbanized areas of North America, Europe, and Central and North Asia. That this macro-region is problematic is also indirectly confirmed by the late 20th century tree-based proxy temperature reconstructions covering the 30° N–70° N land area, which produce significantly less warming than what was shown by the instrumental temperature record covering the same region at least since 1980.

Finally, we compared global land and sea surface temperature records and found that, while the models slightly underestimate the land warming observed from 1940 to 1960 to 2000–2020, they significantly overestimate the SST warming during the same period. In the light of the above findings, and under the assumption that the SST warming since 1940–1960 is accurate, the models can be scaled on the SST record and used to estimate an expected land warming. Corrected in such a way, we determined that 25–45% of the recorded 0.97 ± 0.05 °C land warming from 1940 to 1960 to 2000–2020 is likely due to urbanization and other unidentified non-climatic factors biasing the available climatic records (cf. McKitrick and Michaels 2007).

Finally, we produced a corrected global surface temperature record using a combination of the original SST record and the corrected land temperature record. We found that about 20% of the reported global surface warming from 1940 to 1960 to 2000–2020 could be due to non-climatic biases such as urbanization. Our result appears also experimentally confirmed by the fact that the warming observed since 1979 in our corrected global surface temperature record is statistically compatible with that shown by the satellite UAH MSU v6.0 low troposphere global temperature record, which should not be significantly affected by UHI and other surface non-climatic biases.

Thus, the present findings stress the importance of better investigating possible urbanization and other non-climatic biases throughout the world that could have exaggerated the warming so that the climatic records could be properly corrected. Alternatively, the UAH MSU temperature record could be used to better characterize global warming and climate changes since 1979.

The result is important also to evaluate and improve the climate models. In fact, Fig. 14 suggests that, on average, from the twenty-year period 1940–1960 to 2000–2020 the CMIP5 GCMs have overestimated the anthropogenic warming trend by on average about 40%.

The result has obvious consequences also for the models’ warming expectations for the 21st century because, to make these models consistent with our proposed adjusted temperature record, their projecterd warming should be reduced by about 40% for all emission scenarios.

The latter finding also supports the conclusion of several studies suggesting that these models significantly overestimated the ECS by about a factor of two (cf.: Scafetta 2013; Gervais 2016; Lewis and Curry 2018; Bates 2016; Lindzen and Choi 2011; and others). The models predict ECS values spaning within the range from 2.1 °C to 4.7 °C, which is usually larger than what is determined by several recent empirical studies (Knutti et al. 2017). Moreover, the situation may be worsening because the recent CMIP6 models even predict on average larger ECS values (spanning 1.8–5.6 °C) than the CMIP5 models (Zelinka et al. 2020). Finally, Scafetta (2013) argued that ECS should be low because the temperature reconstructions of the last several millennia also show a large quasi-millennial cycle responsible, for example, of very warm Roman and Medieval climatic anomalies (cf.: Alley 2004; Ljungqvist 2010; Esper et al. 2012; Kutschera et al. 2017; Hao et al. 2020; Margaritelli et al. 2020). This large millennial oscillation is expected to peak in the second half of the 21st century (Scafetta 2013) and appears to be induced by some kind of solar-astronomical forcings (Scafetta 2020; Scafetta et al. 2020). A large natural millennial oscillation is not reproduced by the models, but it should have contributed to the warming observed in the 20th century despite the CMIP5 models attribute the post industrialization global warming intirely to anthropogenic forcings (cf.: IPCC 2013; Scafetta 2019).

Appendix

Colour-blind readers could have some difficulty in correctly interpret Figs. 5, 6, 7 and 9 because the multicolor images also include red and green colors. Using specific software, I simulate how these figures would appear to various types of color-blindeness and it seems to me that the different marked regions could be still recognizable although they would appear in different colors. However, specific software tools are freely available to transform image colors into those discriminable for various people with color vision deficiency: e.g. Visolve, http://www.ryobi-sol.co.jp/visolve/en/.