Increases in atmospheric greenhouse gases will not only raise Earth’s temperature but may also change its variability and seasonal cycle. Here CMIP5 model data are analyzed to quantify these changes in surface air temperature (Tas) and investigate the underlying processes. The models capture well the mean Tas seasonal cycle and variability and their changes in reanalysis, which shows decreasing Tas seasonal amplitudes and variability over the Arctic and Southern Ocean from 1979 to 2017. Daily Tas variability and seasonal amplitude are projected to decrease in the twenty-first century at high latitudes (except for boreal summer when Tas variability increases) but increase at low latitudes. The day of the maximum or minimum Tas shows large delays over high-latitude oceans, while it changes little at low latitudes. These Tas changes at high latitudes are linked to the polar amplification of warming and sea ice loss, which cause larger warming in winter than summer due to extra heating from the ocean during the cold season. Reduced sea ice cover also decreases its ability to cause Tas variations, contributing to the decreased Tas variability at high latitudes. Over low–midlatitude oceans, larger increases in surface evaporation in winter than summer (due to strong winter winds, strengthened winter winds in the Southern Hemisphere, and increased winter surface humidity gradients over the Northern Hemisphere low latitudes), coupled with strong ocean mixing in winter, lead to smaller surface warming in winter than summer and thus increased seasonal amplitudes there. These changes result in narrower (wider) Tas distributions over the high (low) latitudes, which may have important implications for other related fields.
Increases in atmospheric CO2 and other greenhouse gases (GHGs) will not only cause global warming, but may also change the variability in surface temperature and other fields, which have important implications for human society, especially when considering extremes (Katz and Brown 1992; Schär et al. 2004; Fischer and Schär 2009; Hansen et al. 2012; Huntingford et al. 2013). Changes in the probability density functions (PDFs) of surface air temperature (Tas) and other variables have been widely used to quantify variability changes (e.g., Alexander et al. 2006; Donat and Alexander 2012; Fischer et al. 2013; Raghavendra et al. 2019). A shift of the mean of a PDF would lead to a large percentage change in the extremes on the two tails of the PDF, but a change in the shape of the PDF (Hansen et al. 2012), such as a flattening in the PDFs of soil moisture content and runoff under GHG-induced warming (Zhao and Dai 2015), could also induce large changes to the extremes. Changes in the shape of a PDF would imply a change in the variability of the underlying variable. Thus, studying changes in the variability of Tas has major implications for the occurrence of extreme temperature events and other closely related fields, such as atmospheric water vapor (Dai 2016).
Some observational studies have revealed changes in Tas variability over different time periods. For example, Huntingford et al. (2013) found substantial changes in Tas variability over the globe, with increased variability at midlatitudes in both hemispheres from the period of 1963–80 to 1981–96. Screen (2014) revealed a significant decrease in Tas variability during the cold season over the northern mid–high latitudes during 2004–13. Many studies have examined model-projected future changes in Tas variability as the climate warms. For example, Schneider et al. (2015) found a reduction in Tas variability over Northern Hemisphere midlatitudes in climate model projections with increased GHGs, especially in winter, and they attributed it to the reduced meridional thermal advection associated with enhanced Arctic warming. Brown et al. (2017) found reduced variability of the global-mean Tas in the GFDL CM3 2 × CO2 runs. Other studies (Räisänen 2002; Boer 2009; Holmes et al. 2016) showed that Tas interannual variability would decrease over the extratropics (except for northern midlatitudes in summer) but increase in the tropics and subtropics under future warming scenarios. Ylhäisi and Räisänen (2014) also found that future variability of daily Tas would decrease over northern mid–high-latitude land in CMIP3 model simulations. In Europe, daily Tas variability is projected to decrease in winter (de Vries et al. 2012) but increase in summer (Fischer and Schär 2009) during the twenty-first century with increasing GHGs.
The underlying mechanisms for the Tas variability changes have also been examined. Gregory and Mitchell (1995) found that in Europe the winter Tas variability reduction is related to a reduced land–sea temperature contrast while the summer variability increase is caused by reduced evaporative cooling because of the decreased soil moisture in their 2 × CO2 climate. Tas variability decreases in mid–high latitudes (except for northern midlatitudes in summer) are linked to sea ice loss and the associated amplified Arctic warming (Stouffer and Wetherald 2007; Screen 2014; Frankcombe et al. 2018), in conjunction with a decreased latitudinal temperature gradient (Rind et al. 1989; Schneider et al. 2015) and a weaker albedo feedback (Brown et al. 2017). The decreased subseasonal Tas variability in cold-season northern mid–high latitudes is related to the Arctic amplification (AA)—the enhanced warming of the Arctic region compared with lower latitudes (Screen 2014; Screen et al. 2015). AA is seen mainly in the cold season but absent in summer in ERA-Interim and phase 5 of the Coupled Model Intercomparison Project (CMIP5) model simulations (Screen and Simmonds 2010; Dai et al. 2019). AA reduces meridional temperature gradients in northern mid–high latitudes and thus may weaken thermal advection and contribute to Tas variability reduction at northern midlatitudes (Schneider et al. 2015). It may also weaken westerly winds and cause other atmospheric circulation changes (Barnes and Screen 2015; Francis and Vavrus 2015; Luo et al. 2018), which may affect midlatitude weather, including temperature variability and extremes (de Vries et al. 2012; Screen 2014; Holmes et al. 2016; Yao et al. 2017). Qian et al. (2011b) analyzed the spatial patterns of the changes in subseasonal temperature variability over China and suggested that they may be associated with changes in East Asian winter and summer monsoons. In addition, temperature variability changes in summer were discussed in several studies and found to be related to the radiative balance and land surface processes (Gregory and Mitchell 1995; Lenderink et al. 2007; Vidale et al. 2007; Fischer and Schär 2009; Fischer et al. 2012). The increased tropical interannual variability may be related to changes in ENSO and the South Asian monsoon (Meehl et al. 1994).
Seasonal variations often dominate the total variance in daily data for many climate variables outside the tropics (Thomson 1995; Stine et al. 2009; Qian et al. 2011a); thus, studying the changes in the seasonal phase and amplitude is important for examining Tas variability change and for detecting human-caused climate changes (Santer et al. 1991; Thomson 2009; Qian and Zhang 2015; Marvel et al. 2017; Santer et al. 2018). Previous studies have shown substantial changes in the seasonal cycle of Tas over many regions. For example, a trend toward decreased amplitude and advanced phase of the seasonal cycle in Northern Hemisphere mean Tas during the twentieth century was found by Mann and Park (1996) and Wallace and Osborn (2002). Using monthly data only, Stine et al. (2009) found that the phase of the Tas seasonal cycle over extratropical land shifted toward an earlier time by 1.7 days and the amplitude also decreased by 3% from 1954 to 2007. Cassou and Cattiaux (2016) also showed an earlier onset of the summer date over Europe by about 10 days from the 1960s to the 2000s. Besides the changes in Tas, Santer et al. (2018) also analyzed tropospheric air temperature and found that its seasonal amplitude increased at midlatitudes in both hemispheres, but decreased at high latitudes in the Southern Hemisphere, and changed little in the tropics from 1979 to 2016.
Besides the observational studies, model-projected changes in Tas seasonal cycles have also been analyzed. For example, Dwyer et al. (2012) analyzed the future changes in Tas seasonal cycles under increasing GHGs using monthly Tas data. They found a phase delay and an amplitude reduction at high latitudes, and a small phase delay but with increased amplitude over the tropics and midlatitudes during the twenty-first century. They attributed the high-latitude changes to increased surface effective heat capacity associated with sea ice loss at high latitudes and linked the low-latitude changes to changes in surface heat fluxes, but they did not reveal how the surface flux changes may impact the Tas seasonal cycle at low latitudes. Because it is often the change in surface energy fluxes that determine future Tas while the change in the heat capacity of the surface layer is small for most places (e.g., over low–midlatitude oceans), in this study we focus on the surface flux changes rather than the heat capacity changes for explaining Tas changes.
Other studies also analyzed model-projected seasonal cycle changes. For example, McKinnon et al. (2013) used a Lagrangian trajectory model and a simple energy balance model to explain some features of the Tas seasonal cycle using the coupling of land and ocean through atmospheric circulation. Yettella and England (2018) found that large internal variability could delay the detectability of the GHG-forced changes in Tas seasonal cycle over Northern Hemisphere land. Besides Tas, Donohoe and Battisti (2013) also examined the change in seasonal atmospheric heating due to CO2 doubling and found that the seasonal amplitude of air temperature would decrease at the surface in polar regions with sea ice loss but increase in the upper troposphere with increased atmospheric shortwave absorption. Delayed tropical precipitation was also found in a future warmer climate by Biasutti and Sobel (2009) in CMIP3 models that was attributed to a delayed sea surface temperature (SST) response, and by Song et al. (2018) in CMIP5 models that was attributed to increased atmospheric inertia associated with more water vapor. A seasonality change in northeast Pacific atmospheric rivers is implied by a shift of about one month earlier in their peak frequency from the period 1970–99 to 2070–99 under the RCP8.5 scenario (Warner and Mass 2017). A sharpened seasonal cycle with an enhanced amplitude of the monthly mean precipitation over California at the end of the twenty-first century has been revealed by Swain et al. (2018), which may result from the seasonally dependent changes in land–sea moisture contrast (Dong et al. 2019).
Some of these studies found that the seasonal cycle changes are related to changes in Northern Hemisphere atmospheric circulations like the northern annular mode (NAM) and the Pacific–North American (PNA) pattern (Stine et al. 2009; Ault et al. 2011; Stine and Huybers 2012). Other mechanisms are related to changes in surface radiation (Fischer and Schär 2009; Qian et al. 2011a), sea ice and snow cover (Bye et al. 2013; Donohoe and Battisti 2013; Cassou and Cattiaux 2016), land–ocean coupling (McKinnon et al. 2013), and synoptic variability over different regions (Cornes et al. 2017). A human influence is detected in the weakening in Tas seasonal cycle during 1950–2005 over the Northern Hemisphere mid–high-latitude land (Qian and Zhang 2015).
These previous studies often used monthly data to study the changes in Tas variability (e.g., Holmes et al. 2016) and its seasonal cycle (e.g., Stine et al. 2009; Stine and Huybers 2012; Dwyer et al. 2012). The changes in Tas variability and seasonal cycle using daily data remain understudied in comparison with the large number of studies on global warming (i.e., increases) of surface temperature. This is especially true given the potential implications and links of the Tas variability change to other variables. For example, the Tas variability change can directly affect extremes in temperature, water vapor, and precipitation, and may be linked the flattening of the PDFs in soil moisture and runoff (Zhao and Dai 2015). Furthermore, the exact physical mechanisms underlying the Tas variability and seasonal cycle changes are not well understood. For example, while many studies (e.g., Stouffer and Wetherald 2007; Dwyer et al. 2012; Screen 2014; Frankcombe et al. 2018) linked the high-latitude changes to sea ice loss, it is still unclear how sea ice loss would weaken the high-latitude Tas variability and seasonal cycle. While the increased seasonal amplitude at low latitudes is linked to changes in surface fluxes (Dwyer et al. 2012), how the surface flux changes would increase the seasonal amplitude requires further analyses. The previous studies also did not address whether the projected changes over the globe are already evident in historical changes during the recent decades, although some studies (e.g., Qian and Zhang 2015) found a detectable human influence in the weakening of the Tas seasonality over Northern Hemisphere mid–high-latitude land while others (e.g., Yettella and England 2018) suggested that internal variability will delay the detectability of the Tas seasonality change over many Northern Hemisphere regions into the mid-twenty-first century. Furthermore, while the projected change patterns for the Tas variability and seasonal cycle are broadly consistent among the various cited studies, they are often attributed to different processes by different authors and a consistent explanation of the global changes in Tas variability and seasonal cycle is still lacking. Here we attempt to address many of these issues.
Our study provides more detailed (e.g., related to Arctic energy fluxes) and new (e.g., related to evaporative cooling over low-latitude oceans) physical explanations of the changes in Tas seasonal cycle and variability. It differs from the previous studies in many important ways. For example, we averaged daily Tas data to define a composite seasonal cycle instead of Fourier fitting of monthly data in Stine et al. (2009), Stine and Huybers (2012), and Dwyer et al. (2012). We also used a novel empirical orthogonal function (EOF) analysis of the combined monthly Tas data series to graphically show that the seasonal amplitude changes as a result of different warming rates among the cold and warm months. Another new aspect is that we compared the changes in the seasonal amplitude and the standard deviation (SD) of daily Tas anomalies. We also examined the histogram changes for daily Tas anomalies, which are not done in most of the cited studies.
We first evaluate CMIP5 models’ ability in simulating the observed spatial patterns of daily Tas variability and seasonal amplitude and phase and quantify their changes during 1979–2017 using ERA-Interim reanalysis data and CMIP5 historical simulations. We then analyze the CMIP5 model-projected multimodel mean changes in Tas variability and seasonal cycle over the globe under the RCP8.5 high emissions scenario during the twenty-first century, and further investigate the underlying mechanisms. The focus here is on the externally forced change, not variations due to internal variability. The data and methods used in this study are described in section 2. Model evaluations and historical changes are discussed in section 3. Section 4 describes the model-projected future changes in Tas variability and seasonal cycle, and the underlying mechanisms are discussed in section 5. Conclusions and some discussions are presented in section 6.
2. Data and methods
We used daily 2-m air temperature (T2m, used as Tas) data from ERA-Interim reanalysis on a 2.5° grid from 1979 to 2017 (Dee et al. 2011; ECMWF 2011). Unlike most other reanalyses, the ERA-Interim utilized surface observations of air temperatures and humidity; as a result, its Tas data are close to observations (Simmons et al. 2010). As observational global datasets of daily Tas are not readily available, we used the ERA-Interim T2m daily data for evaluating the CMIP5 models and for examining recent changes in Tas variability and seasonal cycle. We also used daily Tas data from historical and future simulations by 25 fully coupled climate models (Table 1) included in CMIP5 (Taylor et al. 2012). The period from 1979 to 1998 was used for evaluating the models’ performance in simulating the current variability and seasonal cycle of Tas, and the differences between the periods of 1999–2017 and 1979–98 were used to quantify recent changes. To investigate the future projected changes in the Tas variability and seasonal cycle, three 30-yr periods 1950–79, 2010–39, and 2070–99 were compared with the model data from the twentieth-century historical simulations for 1950–79 and twenty-first-century simulations under the RCP8.5 high emissions scenario for the two future periods. We used only one realization from each model for the daily data analysis.
For each period (e.g., 1979–98 or 1950–79), the mean seasonal cycle was defined using the time period averaged daily Tas for each calendar day, and the daily Tas anomalies were obtained by subtracting its own mean seasonal cycle (e.g., for 2070–99) from the original daily Tas. For each model, we calculated the mean seasonal cycle and obtained the daily Tas anomalies at each grid point on the original model grid. The PDF of the Tas anomalies, which is approximately Gaussian, and its SD for each season and annual mean were estimated to quantify the Tas variability (excluding the mean seasonal cycle) for each model. Thus, the SD includes variations and changes from both internal variability (IV) and external forcing, similar to that estimated using ERA-Interim. We then interpolated the variability (viz., SD) and the mean seasonal cycle on the original model grid onto the T42 grid (of 2.8125° longitude × ~2.8125° latitude) that is used by the BCC_CSM1.1 and four other models (Table 1) using bilinear interpolation. The multimodel ensemble mean results were obtained from the values on the T42 grid. For each model, we used the seasonal amplitude (defined as the maximum minus minimum Tas in period-averaged daily Tas) and the day of the maximum and minimum Tas to quantify the mean seasonal cycle for each time period. Note that the averaging over a 20- or 30-yr period should remove most of the synoptic noise and produce a stable composite seasonal cycle, which is confirmed in our examination of the Tas seasonal cycles at the select sites indicated in Fig. 1. Some smaller short-term variations still remain in the composite seasonal cycle, but they only have a minor effect on our results. Some of the models had a 360-day calendar and we simply scaled it by multiplying its time with 365/360 to the 365-day calendar, with 29 February being ignored in all the calculations (most models do not have this day). The phase angle θ (=calendar day/365 × 360°) of the maximum or minimum Tas of the seasonal cycle was first projected onto the x = cos(θ) and y = sin(θ) components, which were then averaged over the models, and the averaged x and y components were then converted back to a phase angle that was used as the multimodel mean phase.
We selected 12 locations (Fig. 1) over the world to further show the 30-yr-mean seasonal cycle and the PDF of daily Tas anomalies for the three 30-yr periods. For each location, data from the nearest grid point of each model were used for calculating the ensemble mean to avoid interpolation errors. The range and bin size of the Tas used for its PDF differ among the locations, as daily Tas has a smaller range in the tropics than the extratropics. However, we used a fixed number (52) of bins in estimating the histograms of the PDFs.
We also used CMIP5 monthly data for Tas, sea ice concentration (SIC), surface upwelling longwave radiation, surface upward latent and sensible heat fluxes, 10-m-height wind speed, 2-m-height specific humidity, surface temperature and pressure to analyze possible mechanisms for the changes in the seasonal cycle. An EOF analysis of the monthly anomaly time series (with all months concatenated) was performed to reveal the different change rates among the months that contribute to the seasonality changes. For the monthly data analysis, we averaged the data over several ensemble runs for some of the models before averaging over all the models (after remapping to a 2.5° grid) to produce the multimodel ensemble mean, on which the EOF analysis was done, with the total number of model runs indicated in related figure captions.
3. Recent climatology and change of Tas variability and seasonal cycle
Before analyzing the projected future changes, it is necessary to evaluate the models’ ability in simulating the distributions and changes during the historical period as seen in ERA-Interim. Two periods were used in this analysis here: 1979–98 and 1999–2017. We emphasize that ERA-Interim contains changes resulting from both internal variability and external forcing, while the CMIP5 ensemble mean contains mostly forced response. Thus, they are not entirely comparable. Nevertheless, such a comparison can still help us identify whether a regional change pattern in ERA-Interim is similar to that in the CMIP5 ensemble mean and thus is likely caused by historical external forcing, or if it differs from the model ensemble mean and thus is likely caused by internal variability.
The amplitude of the mean seasonal cycle from both ERA-Interim (Fig. 2a) and CMIP5 models (Fig. 2b) is larger over land than over ocean and increases from the tropics to polar regions, as the seasonal amplitude of solar radiation increases with latitudes. The seasonal amplitude increases from west to east over the midlatitude continents as the oceanic influence weakens downwind, but it decreases from west to east over midlatitude ocean basins because the continental influence weakens as the westerlies leave the coastal lines. The CMIP5 models capture well the seasonal amplitude distribution seen in ERA-Interim, with maximum amplitudes over northeastern Asia and northern North America and minimum amplitudes over the western and central tropical Pacific (Figs. 2a,b).
The Tas variability excluding the mean seasonal cycle during 1979–98, expressed as the SD of daily Tas anomalies (Figs. 2c,d), is also larger over land than over ocean and larger in the high latitudes than in the low latitudes. Maximum variability with SD > 6°C is seen in north-central Eurasia and northwestern North America, while the Tas variability is lowest over equatorial oceans with SD < 0.8°C (Fig. 2c). The Tas over ice-free ocean is tightly coupled with SSTs, which are less variable due to ocean’s large thermal inertia. This helps explain the small Tas variability over ocean in comparison with land. The multimodel-averaged SD (Fig. 2d) shows very similar spatial variations as seen in ERA-Interim (Fig. 2c), with some small differences in magnitude. For example, the SD of the Tas anomalies is slightly lower over tropical oceans in the models than in ERA-Interim, while it is the opposite in the Arctic region. Overall, the CMIP5 multimodel ensemble mean captures well the Tas variability patterns revealed by ERA-Interim.
Figures 3a and 3b show the percentage change of the seasonal Tas amplitude from the period of 1979–98 to 1999–2017 in ERA-Interim and CMIP5 multimodel ensemble mean. Note that the amplitudes in Figs. 3a and 3b (and the SD in Figs. 3c and 3d) both include internal variability and forced changes, but the internal variability (IV)-induced changes in the amplitude and SD are uncorrelated among the models and thus should be canceled out during the averaging over the models, leaving mainly the forced changes in Figs. 3b and 3d. The ensemble averaging over the models should also lead to smoother patterns. Nevertheless, about 64% of the grids over the globe have the same sign of change between the reanalysis and the model ensemble mean. Thus, the overall change patterns are broadly similar between ERA-Interim and the models, with decreases at high latitudes and increases over many oceans and low-latitude land areas (Figs. 3a,b). This implies that these changes are likely externally forced because the chance for the IV in ERA-Interim to generate amplitude changes similar to those of the CMIP5 models is very low.
The percentage changes in Tas SD from the period 1979–98 to 1999–2017 from ERA-Interim and CMIP5 multimodels are shown in Figs. 3c and 3d. Again, the SD change for ERA-Interim include IV-induced decadal variations, while the model-averaged change represents mainly forced response. Nevertheless, about 57% of the global grids have the same sign of change. Both ERA-Interim and CMIP5 models show decreased Tas variability at high latitudes and increased variability over many low-latitude areas (Figs. 3c,d). In particular, a large reduction in Tas variability over the Barents–Kara Sea (BKS) is seen in both ERA-Interim and CMIP5 models (Figs. 3c,d), which is likely linked to the decreasing sea ice cover over the region (Dai et al. 2019).
The models show large increases in Tas variability (mainly in response to historical external forcing) over the central and eastern tropical Pacific (Fig. 3d). In contrast, ERA-Interim shows decreasing Tas variability over the eastern equatorial Pacific and tropical Atlantic (Fig. 3c). This suggests that internal variability, such as that associated with the recent phase change in the interdecadal Pacific oscillation (IPO; Dai et al. 2015), still contributes significantly to the Tas variability change in low latitudes in individual realizations such as ERA-Interim, which also shows noisier change patterns than the multimodel ensemble mean. Specifically, the phase transition of IPO around 1998/99 led to the recent cooling and reduced interannual variability (Hu et al. 2013) in the eastern tropical Pacific Ocean (models are unable to capture this unforced change), which overwhelms the GHG-induced global warming and the effect of the Atlantic multidecadal oscillation (AMO) (Dong and Zhou 2014). The decreasing trend in Tas variability over the eastern tropical Pacific shown in Fig. 3c closely matches the decreased SST variance during 2000–11 shown by Hu et al. (2013). Hu et al. (2013) proposed that stronger surface trade winds and a steeper thermocline slope hamper the eastward migration of the warm water along the equatorial Pacific during 2000–11 and therefore reduce the temperature (and other) variability in the eastern tropical Pacific. Thus, the recent changes in Tas variability from the period 1970–98 to 1999–2017 show robust decreases at high latitudes (especially over the BKS with large sea ice loss) mainly in response to external forcing, while the changes at low latitudes are noisy and are still heavily influenced by internal variability.
We define the date of the maximum and minimum Tas in the 20-yr-mean seasonal cycle to represent the seasonal phase during 1979–98. Substantial contrasts between land and ocean (except the Arctic Ocean) are observed with later dates of maximum (minimum) Tas over the oceans in the Northern (Southern) Hemisphere (Fig. 4a). The daily Tas peaks around the 180–200th day of the year over the Arctic Ocean, around the 190–210th day of the year over most Northern Hemisphere land except South Asia and tropical Africa where it peaks around the 170–180th day of the year. The Tas peak is around the 220–260th day of the year in the northern North Atlantic and northern Pacific, and around the 260–280th day of the year in their tropical regions (Fig. 4a). The daily minimum Tas occurs around the 30–60th days of the year over the Arctic Ocean, around the 350–365th or 1–20th day of the year over most Northern Hemisphere land, and around the 50–90th day of the year over the North Atlantic and North Pacific (Fig. 4c). The daily minimum (maximum) Tas over the Southern Hemisphere (except for Antarctica) occurs around the time of the maximum (minimum) over the Northern Hemisphere (Figs. 4a,c). CMIP5 models capture the general patterns of the maximum (minimum) date distributions but with some regional biases (e.g., over the western South Pacific, Figs. 4b and 4d). The changes in the time of the maximum (minimum) Tas from ERA-Interim (not shown) are noisy without coherent regional patterns, and the models also do not show robust change patterns. This differs from Stine et al. (2009), who found that the seasonal cycle in extratropical land Tas shifted earlier by about 1.7 days from 1954 to 2007 based on their Fourier fitting of monthly data. However, we did find some seasonal phase shifts (toward earlier spring) by 0–5 days during 1979–2017 over most extratropical land areas (except for Canada) when a Fourier transform was applied to the ERA-Interim daily Tas data (not shown).
In summary, the CMIP5 models can capture the spatial patterns of daily Tas variability and seasonal cycle, but not all the changes from the period 1979–98 to 1999–2017. The only robust change pattern seen in both ERA-Interim and CMIP5 models appears to be the decreasing Tas variability and seasonal amplitude over the Arctic and Southern Ocean, which is likely related to decreasing sea ice there (Screen 2014; Frankcombe et al. 2018). Over most other regions, changes induced by internal variability appear to dominate over the forced changes during the recent period in both ERA-Interim and the CMIP5 models, as reflected by the low consistency even for the sign of the change among the CMIP5 models over these regions.
Previous studies have also found that climate models can capture the observed mean seasonal cycle but failed to reproduce many of the recent changes in seasonal phase and amplitude (Mann and Park 1996; Wallace and Osborn 2002; Stine et al. 2009; Ylhäisi and Räisänen 2014). For example, Stine et al. (2009) found that few of the CMIP3 models could reproduce the observed amplitude decrease and none could reproduce the shift toward earlier seasons over extratropical land during 1954–2007. This failure to reproduce the observed changes may be due to the fact that these changes are caused primarily by realization-dependent internal variability (such as the northern annular mode and Pacific–North American pattern; Stine and Huybers 2012) and thus they are not comparable with the randomly initialized model simulations. Comparisons of the recent changes need to account for the realization-dependent internal variations that can dominate over the externally forced changes up to present over many regions (Dai and Bloecker 2019). Since these unforced internal variations are expected to differ among different realizations (such as the observations and individual model runs), substantial differences in the observed and model-simulated historical changes in any given field should be expected over regions where the forced change is still not much larger than internal variability (Dai and Bloecker 2019). Many studies have made misleading conclusions regarding climate models’ ability in simulating historical temperature, precipitation and other changes because the authors incorrectly expected the models to reproduce the observed changes without realizing that these changes may be dominated by internally generated decadal–multidecadal variations that are not reproducible by the CMIP3 or CMIP5 models. These model simulations were designed to simulate the externally forced changes, not to reproduce the observed internal variations, as they were initialized with random initial conditions for 1850.
4. Model-projected changes in Tas daily variability and seasonal cycle
The spatial patterns of the amplitude of the seasonal cycle for these periods (not shown) are similar to Fig. 2b. Here we only show its percentage changes. Increasing (decreasing) amplitudes over low-latitude land and low–midlatitude oceans (northern midlatitude land and the Arctic and Southern Ocean) are projected in the future periods (Figs. 5a,c). These amplitude changes are robust among the models with more than 80% of them agree on the sign of the change, especially by 2070–99 (Fig. 5c). The spatial patterns of the SD in Tas anomalies for these periods (not shown) are similar to Fig. 2d. From the period 1950–79 to 2010–39, the Tas variability increases by 5%–20% over low-latitude land and most oceans, except for the northern high latitudes and Southern Ocean where it decreases by 5%–20% (Fig. 5b). By 2070–99, these change patterns become more robust than during 2010–39, with 10%–30% increases (decreases) in the respective regions (Fig. 5d). The projected decreasing variability at high latitudes and increasing variability at low latitudes are consistent among the models as indicated by the stippling in Figs. 5b and 5d.
The spatial patterns of the dates of the maximum (minimum) of the seasonal cycle for these periods (not shown) are similar to Figs. 4b and 4d. For the future periods, the dates for both the maximum and minimum Tas are delayed by 15–30 days over the Arctic Ocean and the Southern Ocean (for date of minimum only) where sea ice loss is significant, while the date changes are much smaller with noisy patterns over other regions (Fig. 6).
Figure 7 shows the 30-yr-mean seasonal cycle and the PDFs of the Tas anomalies for 12 selected locations with relatively large changes (Fig. 1). Figure 7 confirms that the changes shown in Figs. 5a and 5c, with reduced seasonal amplitudes in the late twenty-first century at the high-latitude locations (i.e., A1, A2, A3, A12) and increased amplitudes at oceanic locations (i.e., A6, A7, A9) and low-latitude land locations (i.e., A10, A11). The seasonal cycle over midlatitude land locations (i.e., A4, A5) changes little. Figure 7 also shows that the mean seasonal cycle may differ substantially from a sine curve, making a Fourier fit imperfect.
The PDFs of the Tas anomalies show Gaussian-like distributions (Fig. 7), as noticed previously (e.g., Schneider et al. 2015). This indicates that the SD changes can well represent the variability changes in daily Tas anomalies under global warming. The PDFs at some selected high-latitude locations (i.e., A1, A2, A3, A4, A12) show higher peaks with narrower distributions by the end of the twenty-first century than the 1950–79 period (Fig. 7), which indicates reduced variability. However, the projected PDFs become broader with lower peaks for some low-latitude locations (i.e., A8, A9, A10, A11) (Fig. 7), which suggests increased variability with warmer temperatures in the future (the increase in the mean Tas would make the low Tas extremes less frequent). Other locations (i.e., A5, A6, A7) show small changes. There is a transition zone between the high and low latitudes where the seasonal or PDF changes are small.
We also calculated the SD of the daily Tas anomalies (i.e., without the mean seasonal cycle) in boreal winter [December–February (DJF)] and summer [June–August (JJA)] separately (Fig. 8). The Tas anomaly variability is larger in the winter than in the summer over the northern extratropics and the Antarctic region (Figs. 8a,b), which is consistent with the interannual variability related to the large latitudinal temperature gradient in winter (Hansen et al. 2012). The percentage changes in DJF (Figs. 8c,e) show results similar to the annual case (Figs. 5b,d). But in JJA, the variability shows small increases in northern mid–high latitudes (except for the polar region) from the period 1950–79 to 2010–39 and the increases become larger from the period 1950–79 to 2070–99 (Figs. 8d,f). The patterns of the Tas variability and projected percentage changes for March–May and September–November (not shown) are similar to the annual case (Figs. 5b,d) with some differences in magnitude. The seasonal PDFs of the Tas anomalies (DJF and JJA in Fig. 9; MAM and SON not shown) also show differences mainly between JJA and the other seasons for locations in northern mid–high latitudes (i.e., A1, A2, A3, A5), with projected lower (higher) peaks and broader (narrower) distributions in JJA (other seasons). The PDFs for other locations generally show consistent distribution changes across the seasons but with some difference in magnitude.
To further investigate the changes in Tas variability, we calculated the mean Tas anomalies for the 5% coldest and 5% warmest days during 1950–79 and 2070–99 for each season. Figure 10a shows that the coldest days in DJF would become warmer at high latitudes around 60°N and 60°S but become colder at low latitudes. For the warmest DJF days (Fig. 10c), roughly the opposite occurs. Thus, these changes would lead to decreased variability at high latitudes, but increased variability at low latitudes in DJF, consistent with the SD change (Fig. 8e). In JJA (Figs. 10b,d), the coldest days become colder while the warmest days become warmer over land and most low-latitude oceans. This would lead to increased variability over these areas as seen in Fig. 8f, while it is the opposite around 60°S, leading to decreased variability there (Fig. 8f). The Tas changes for the coldest and warmest days in MAM and SON (not shown) are similar to DJF. Note that the increased variability in JJA is masked by the decreased variability in other three seasons, leading to decreased variability for the annual mean over northern mid–high latitudes (Fig. 5d). Therefore, there exist asymmetric warming rates for the coldest and warmest days and they can help explain the variability changes in daily Tas anomalies.
5. Underlying mechanisms for the seasonality changes
Observational studies found an overall larger warming in winter than in summer during the twentieth century (Wallace et al. 1995; Balling et al. 1998; Trenberth et al. 2007; Hartmann et al. 2013; Nigam et al. 2017), which is also the case for model projections (Collins et al. 2013), although the recent decades since the 1990s show cooling trends over many Northern Hemisphere land areas in winter (Cohen et al. 2012). This motivated us to focus on varying warming rates among different months over different regions. To do so, we performed an EOF analysis of the monthly Tas anomaly data (with all months concatenated) from the multimodel ensemble mean, as well as other relevant fields, such as polar SIC. The first leading EOF mode (EOF1) of Tas over the northern high latitudes (50°–90°N) is characterized by large values from the eastern Siberian Sea to Beaufort Sea and the Hudson Bay, and over the BKS (Fig. 11c). The corresponding time series show upward trends for Tas (Fig. 11a) and downward trends for SIC (Fig. 11b). Furthermore, they show varying change rates among the months, with each colored line in Figs. 11a and 11b (and Figs. 11e and 11f) connecting the data points corresponding to a specific month. For example, Fig. 11a shows that the warming at the northern high latitudes is fastest in December but slowest in June, with all other months falling between them. Similarly, Fig. 11e shows that the warming rate at the southern high latitudes (50°–90°S) is largest in July but smallest in December. Thus, the warming is larger during the cold season than the warm season at high latitudes in both hemispheres, which would dampen the seasonal amplitude at mid–high latitudes as shown in Figs. 5a and 5c.
The warming patterns revealed by EOF1 (Figs. 11c,g) are linked to polar sea ice loss (Figs. 11d,h), with the largest warming over the areas with the largest sea ice loss in both hemispheres. The pattern correlation coefficients between the EOF1s of Tas and SIC are 0.70 and 0.42 for the Northern Hemisphere and Southern Hemisphere, respectively. Previous studies (e.g., Screen and Simmonds 2010; Dai et al. 2019) have shown that sea ice loss during the cold season exposes the polar warm waters to the frigid air and thus provides strong heating of the lower troposphere with large upward longwave radiation and sensible and latent heat fluxes into the air. This extra heating greatly enhances near-surface and lower-tropospheric warming during the cold season at high latitudes. Note that while the ocean is losing energy to heat the air (mainly using the energy stored during the warm season), the ocean surface temperature (Ts) is still warmer than the frigid air above it during the cold season (Dai et al. 2019). Eventually, the seasonal cooling would lead to formation of seasonal ice over most polar regions, which would stop this oceanic heating of the atmosphere and result in small Arctic amplification for February–March when sea ice cover is at maximum (Dai et al. 2019). Thus, the ocean surface cooling (due to energy loss to the air) would not cause a reduction in Tas. This is different from many other locations, where the Ts − Tas gradient is small and any changes in Ts are often reflected quickly in Tas. Our analysis of the surface fluxes (Fig. 12) confirms this mechanism for both the northern and southern high latitudes, that is, sea ice loss exposes the relatively warm water to the frigid polar air during the cold season, causing large releases of latent and sensible heat fluxes and upwelling longwave radiation into the air, which enhance the surface and lower-tropospheric warming during the cold season. The enhanced warming patterns shown in Figs. 11c and 11g come mainly from the cold season, because during the warm season, most of the extra heating (from increased absorption of solar radiation) is stored in the polar ocean mixed layer, which are a heat sink to the polar atmosphere during the warm season, and thus does not lead to enhanced surface warming during the summer (Dai et al. 2019).
Figure 11b shows that sea ice loss is largest during October but smallest during March during 1950–2100 in the Arctic, leading to an increased seasonal amplitude for SIC; while Antarctic sea ice loss is largest in July and smallest in February (Fig. 11f), leading to a reduced seasonal amplitude. Thus, the large polar amplification of surface warming during the cold season caused primarily by sea ice loss (Screen and Simmonds 2010; Kirtman et al. 2013; Bintanja and Van der Linden 2013; Dai et al. 2019) is the main reason for the reduced seasonal amplitude of Tas at high latitudes under GHG-induced warming in the twenty-first century. The heat capacity argument of Dwyer et al. (2012) cannot explain this seasonal change in the polar warming rates and thus cannot explain the seasonal amplitude change there. This is because it is the season-dependent change in the surface energy fluxes, not the change in the heat capacity of the surface layer, that affects the Tas change in the Arctic (Dai et al. 2019). In fact, the heat capacity of the surface ice layer during the cold season should be larger than that of the surface water layer during the warm season, and this should slow down the surface warming rate in the cold season if it has any effect at all.
The seasonal cycle and its changes at low–midlatitudes are much weaker than at high latitudes. However, the percentage changes are comparable to the high latitudes and show amplitude increases over the low–midlatitude oceans (Figs. 5a,c). The low–midlatitudes have a faster warming rate in the warm season than in the cold season for both the hemispheres (not shown), resulting in the amplitude increase. These changes at low–midlatitudes may be linked to changes in surface heat fluxes, as suggested previously (Dwyer et al. 2012). Here we focus on the surface latent heat flux (LH) because the surface net radiation changes little and the sensible heat flux changes show weak spatial correlation with the seasonal amplitude changes (not shown).
Figure 13 shows that LH increases faster during the cold season than the warm season over low–midlatitude oceans and some land areas (Fig. 13). Thus, more of the extra surface heating from increased GHGs would be used for evaporation rather than for raising surface temperatures during the cold season than the warm season, leading to a slower surface warming rate during the cold season and thus contributing to the seasonal amplitude increases at low–midlatitudes (Figs. 5a,c). However, the increased seasonal amplitude over parts of South America and Africa cannot be explained by this process. Another factor is the stronger ocean mixing during the cold season (due to stronger winds) than the warm season, which would slow down ocean surface warming more during the cold season than the warm season. Note that the LH may increase atmospheric water vapor content initially, but atmospheric mixing, advection, and condensation would quickly remove any local water vapor anomalies. Thus, the increased LH may not necessarily lead to extra heating of the surface through downward longwave radiation, while its evaporative cooling effect can directly slow down the surface warming rate. The evaporative cooling and ocean mixing also are the main reasons for the smaller warming over ocean than over land under increasing GHGs (Collins et al. 2013; Dai 2016).
Based on the aerodynamic formula for LH widely used in models: LH = LρCLU(qs − qa), where U is the 10 m wind speed and W ≡ (qs − qa) is the specific humidity difference between the surface and 2-m height, we found that the LH change (dLH) is proportional to U0 × dW + W0 × dU = term1 + term2 = sum, where U0 and W0 are the current climate values and dU and dW are their future change. To diagnose the LH change, we simply used the monthly data from the multimodel ensemble mean to calculate these terms. Figure 14 shows the 1950–79 30-yr-mean seasonal cycle of the U0 and W0 averaged over the low–midlatitude oceans for the Northern and Southern Hemispheres separately and their changes from the period 1950–79 to 2070–99 under the RCP8.5 scenario, and the seasonal LH changes. Note the formula we used for calculating qs (in g g−1) is qs = 0.622es/(ps − 0.378es), es = 6.112 exp[17.67(Ts − 273.16)/(Ts − 29.65)], where es (in hPa) is saturation vapor pressure based on Bolton (1980), ps (in hPa) is surface pressure and Ts (in K) is surface temperature. For the Northern Hemisphere, U0 is stronger in the cold season than in the warm season (Fig. 14a) and the positive dW is also larger in the cold season than in the warm season (Fig. 14c) with a slightly earlier phase compared with U0; thus term1 is positive with larger values in the cold season than in the warm season (Fig. 14e), which favors a faster LH increase during the cold season than the warm season. However, the positive W0 (Fig. 14a) and negative dU (Fig. 14c) are generally out of phase, leading to negative term2 (Fig. 14e), which slows down the LH increase mainly during the cold season. Since term1 is larger than term2, their sum results in stronger LH increase during the cold season than the warm season in the Northern Hemisphere (Fig. 14e). For the Southern Hemisphere, the positive U0 (Fig. 14b) and positive dW (Fig. 14d) result in a positive term1 (Fig. 14f), which leads to stronger LH increases around May than around October. Like the Northern Hemisphere, the positive W0 (Fig. 14b) and negative dU (Fig. 14d) lead to negative term2 (except for August, Fig. 14f) and damps the LH increase mainly during the warm season. The term1 and term2 together contribute to faster LH increases during the cold season than the warm season in the Southern Hemisphere (Fig. 14f). Furthermore, the similarity of the seasonal cycle between the sum of the decomposed two terms (Figs. 14e,f) and the LH flux directly from the models (Figs. 14g,h) suggests that our decomposition method works well.
Figure 15 further shows the zonal-mean results with more details of meridional distributions for the LH decomposition. It shows that the seasonal cycle of term1 (Fig. 15e) is influenced by both U0 (Fig. 15a) and dW (Fig. 15d) in the Northern Hemisphere but is dominated by U0 in the Southern Hemisphere, while the seasonal cycle of term2 (Fig. 15f) is mainly determined by dU (Fig. 15b). Figure 15e shows positive term1 (which leads to LH increases) at all latitudes with larger values in the cold season than in the warm season, while term2 (Fig. 15f) shows both negative and positive values. The sum of term1 and term2 (Fig. 15g) leads to larger LH increases during the cold season than the warm season, which is similar to the seasonal cycle change revealed by the latent heat flux data from the models (Fig. 15h).
The above simple analyses help understand the seasonality of the LH change as revealed by the EOF results shown Fig. 13, which in turn helps explain the faster warming rate in summer than winter over the mid–low-latitude oceans. Sobel and Camargo (2011) suggested that the strengthening (weakening) surface winds in the winter (summer) subtropics contribute to the seasonal amplitude increases of surface latent heat fluxes. Our results show that the strong mean surface winds in the cold season, the strengthening of the surface winds in the cold season in the Southern Hemisphere, and the large increase of surface humidity gradient in the cold season over the Northern Hemisphere low latitudes all contribute to the larger LH increases during the cold season than the warm season, which then damps the warming rate more in winter than in summer.
6. Conclusions and discussion
We first evaluated CMIP5 models’ ability in simulating the current spatial patterns of daily surface air temperature (Tas) variability and seasonal cycle and their recent changes during 1979–2017, as represented by ERA-Interim. We then analyzed CMIP5 multimodel-simulated changes in the Tas variability and seasonal cycle during the twenty-first century under the RCP8.5 high emissions scenario and investigated the underlying causes for these changes. Results show that the models can well capture the spatial patterns of Tas variability and seasonal cycle. ERA-Interim data show decreasing seasonal amplitude from the period 1979–98 to 1999–2017 over the Arctic and Southern Ocean, especially the areas with large sea ice loss, which also reduces daily Tas variability over these areas. Most CMIP5 models reproduce these changes, suggesting that they are likely to be a response to recent GHG and other external forcing. On the other hand, changes in Tas variability and seasonal cycle over the low–midlatitudes are noisy in ERA-Interim and less consistent among the CMIP5 models, which suggests that random internal variations still dominate over the forced change over the low–midlatitudes.
The Tas variability and seasonality changes are examined by comparing the 2010–39 and 2070–99 periods under the RCP8.5 scenario to the 1950–79 period in the historical simulation using daily Tas data from 25 CMIP5 models. Results show increasing seasonal amplitudes over low-latitude land and most oceans (except for the Arctic and Southern Ocean) and decreasing amplitudes over northern high latitudes and the Southern Ocean. For daily Tas without the mean seasonal cycle, its variance would increase at low latitudes but decrease at high latitudes. The variance decrease at northern high latitudes occurs mostly in the cold season, as the Tas variance increases slightly there in boreal summer. These changes in Tas seasonal amplitude and variability are robust and consistent among the CMIP5 models, with more than 80% of the models agreeing on the sign of change over most areas except for the transitional zone. The shifts in the date of maximum or minimum of the seasonal cycle suggest large delays by 15–30 days over the high-latitude ocean areas with significant sea ice loss, especially over the Arctic Ocean. The seasonal phase changes in the subtropics and tropics are small and inconsistent among the models.
These changes in Tas seasonal cycle and variability are confirmed by the composite curves of the seasonal cycle and PDFs of daily Tas anomalies for selected locations. The Tas variability change is also reflected in the asymmetric warming rates for the coldest (bottom 5 percentiles) and warmest (top 5 percentiles) days, which are located on the two tails of the PDF of daily Tas anomalies. For example, during DJF the coldest days warm faster than the warmest days at high latitudes, leading to decreased variability there, whereas it is the opposite at low latitudes, leading to increased variability. Thus, the PDFs of the Tas anomalies tend to narrow at high latitudes but broaden at low latitudes.
The seasonal amplitude changes are also revealed by an EOF analysis of monthly CMIP5 multimodel ensemble mean Tas. Winter months warm faster than summer months at high latitudes in both hemispheres, resulting in reduced seasonal amplitude there; while it is the opposite at low latitudes. At high latitudes, the enhanced warming patterns during the cold season are linked to sea ice loss, which causes large polar amplification of surface warming during the cold season and contribute to the reduced seasonal amplitude. The polar amplification results from the large release of longwave radiation, sensible and latent heat fluxes from the newly exposed ocean waters during the cold season (Dai et al. 2019). Over low–midlatitude oceans, surface latent heat (LH) flux increases faster during the cold season than the warm season. Coupled with stronger ocean mixing during the cold season, this leads to a slower surface warming rate during the cold season than the warm season. The larger increases in winter LH fluxes result from strong winter surface winds in both hemispheres, strengthened winter surface winds in the Southern Hemisphere and large increases in winter near-surface gradients of specific humidity over the Northern Hemisphere low latitudes.
The projected changes in daily Tas variability revealed in this study are consistent with previous results obtained using daily (Ylhäisi and Räisänen 2014; Screen 2014), monthly (Holmes et al. 2016) and annual-mean Tas data (Olonscheck and Notz 2017; Frankcombe et al. 2018). The results suggest that the Tas variability would decrease over mid–high latitudes and increase in the tropics in a warmer climate. The projected seasonal phase delay and amplitude decrease at high latitudes and amplitude increase at low latitudes are also consistent with the results based on CMIP3 model simulations (Dwyer et al. 2012). In particular, the decreasing Tas variability and decreasing seasonal amplitude over the Arctic regions and around Antarctica are very robust responses to GHG-induced warming resulting mainly from sea ice loss, mainly through the heat release by the newly opened waters, rather than changes in surface heat capacity as suggested by Dwyer et al. (2012). They are seen in most CMIP5 models and are already evident during 1979–2017 as seen in ERA-Interim and CMIP5 historical simulations. Decreasing sea ice cover not only enhances the GHG-induced warming during the cold season (Dai et al. 2019) and thus reduces the Tas seasonal amplitude, but also decreases the variability of sea ice and thus the Tas variability associated with it. However, how the reduced sea ice exactly causes the Tas variability to decrease needs further investigation. The decreased Tas variability at northern mid–high latitudes has also been linked to reduced meridional temperature gradient (Schneider et al. 2015) and a weaker albedo feedback (Brown et al. 2017). In this study, we also did not explore the causes for the Tas variability increase at low latitudes (Figs. 5d, 8e,f), as we only showed the different warming rates for the coldest and warmest days at low latitudes (Fig. 10) that confirm the Tas variability increase. What causes these different warming rates and thus the Tas PDF to widen at low latitudes requires further investigation.
These Tas variability and seasonality changes likely influence other climate fields. For example, CMIP5 models project an amplification and a phase delay of the seasonal cycle for tropical sea surface temperatures and precipitation under RCP8.5 scenario in the twenty-first century (Dwyer et al. 2014). A seasonal delay in Sahel rainfall is another regional manifestation of the projected change in the seasonal cycle (Biasutti and Sobel 2009), while the seasonal delay of tropical rainfall may be related to interhemispheric energy contrasts (Song et al. 2018). Since air temperature is closely coupled with many other climate fields, such as water vapor content and precipitation, the Tas variability and seasonality changes, such as the widening (narrowing) of its PDFs at low (high) latitudes revealed here, may have important implications for variability and seasonality changes in other climate fields, such as the flattening of the PDF for soil moisture and runoff (Zhao and Dai 2015). This issue requires further investigation.
This work is partly supported by the National Natural Science Foundation of China (41621005, 41475092). Chen also acknowledges the support from the China Scholarship Council to sponsor her visit in Albany. Dai acknowledges the funding support from the U.S. National Science Foundation (Grant AGS-1353740 and OISE-1743738), the U.S. Department of Energy's Office of Science (Award DE-SC0012602), and the U.S. National Oceanic and Atmospheric Administration (Award NA15OAR4310086 and NA18OAR4310425).