A global perspective on past and future change in regional seasonal cycle shape

Ever-worsening climate change increases near-surface air temperatures for almost the entire Earth and threatens living organisms and human society. While annual mean changes are frequently used to quantify past and expected future changes, the increase is rarely uniform throughout the year. In addition, the shape of the annual cycle and its changes can differ considerably between regions around the globe. Therefore, we perform a global analysis resolving the annual cycle and its changes in different regions, focusing on diagnostics that can be evaluated for the variety of existing annual cycle shapes (e.g., single and double waves, different timing of seasons, etc.). Many previous studies relied on parameter-based methods, assuming a sinusoidal shape of the mean annual cycle. We introduce the Functional Data Analysis (FDA) approach, representing the mean annual cycle by a linear combination of Fourier bases. The FDA methodology does not require any prior assumptions about the shape of the temperature seasonal cycle except periodicity and allows to quantitatively assess various aspects of the seasonal cycle shape. The evolution of the mean annual cycle is estimated from daily long-term mean temperature values, which are converted to functional form. We concentrate on diagnostics that evaluate the absolute change in temperature, its seasonal slope, the position of the maximum, and the amplitude of the annual cycle. We analyze two reanalysis datasets (coupled CERA20C and atmospheric ERA5) and a subset of five CMIP6 Earth system models (ESMs). Observed changes in the second half of the 20th century are assessed, and the ability of ESMs to represent them is evaluated. Further, the changes projected for the end of the 21st century under the SSP3-7.0 pathway are analyzed. Among other results, we highlight distinct differences between the two reanalyses, especially over equatorial and polar regions across diagnostics. Our approach also reveals that differences in the historical period between 1951–1980 and 1981–2010 can be negative during (short) parts of the year in many regions. Further, the ESMs future projections show different rates of warming between seasons, resulting in changes in the amplitude. The largest amplitude increase is projected over the Mediterranean region, and the largest decrease over the Arctic Ocean, the latter being due to the considerably stronger warming in the Northern Hemisphere winter. The ESMs also project a delayed maximum near the poles and an earlier maximum in many tropical continental regions. In Europe, the southern and eastern regions experienced a delay of the maximum of up to 10 d, whereas a slightly earlier maximum is found for northern Europe. A similar dipole pattern can be seen between eastern and western regions in North America. Regarding the slope of the annual cycle, higher latitudes detect a higher magnitude of change in the historical period than lower latitudes. The geographical pattern remains the same for future slope changes, with the magnitude twice as high in most regions. The FDA diagnostics introduced here can be tailored for different purposes and applied to different climatic variables, with no need to make any prior assumptions about the annual cycle shape. Potential applications include, e.g., explicitly evaluating the climate model performance or ensemble mean and spread assessment beyond annual or seasonal means.

1 Introduction

Increasing near-surface air temperature is observed and projected for almost the entire Earth (IPCC, 2021), threatening the environment and human society alike. However, this temperature increase is rarely uniform throughout the year, and even in case the annual mean changes only slightly, the annual cycle might change quite dramatically (Marvel et al., 2021; Wang et al., 2021). The changes in seasonal temperature cycle can have potentially large impacts on, e.g., phenological phases of living organisms, agriculture, health, tourism, and other sectors. Widespread expected changes in the annual cycle have even motivated suggestions of new definitions of seasons (Hekmatzadeh et al., 2020; Wang et al., 2021; López-Franca et al., 2022). Moreover, as noted by McKinnon and Huybers (2024), the shape of the temperature annual cycle can be taken as an analogy for temperature changes in general, as it is easily distinguishable from internal variability and can be reliably observed. They emphasize that the seasonality of temperature in the current climate and its changes are strongly related to projected temperature changes, and recent changes in the mean annual cycle can be considered as a proxy for overall future warming. The skill of climate models in depicting correctly the observed shape of the annual cycle and its changes is therefore very informative in terms of confidence in simulated future changes (Lynch et al., 2016).

A large number of previous studies have shown that the observed shape of the temperature annual cycle has already changed in recent decades, including, e.g., a phase shift towards an earlier onset of the seasons over the middle and higher latitudes (evaluated using the sinusoidal approximation of the mean annual cycle shape, the results do not relate to a specific season, Stine et al., 2009), lengthening of summer (Peña-Ortiz et al., 2015; Park et al., 2018) and shortening of all other seasons over Northern Hemisphere midlatitudes (Wang et al., 2021). Wang and Dillon (2014) revealed regionally different changes of annual cycle amplitude over Northern Hemisphere midlatitudes and polar regions, with a prevailing decrease in 1975–2010 in comparison to 1961–1990. In addition to adaptation to recently observed shifts, it is also crucial to investigate the expected future evolution, as the shape of the annual cycle is expected to undergo even more dramatic changes during the upcoming decades. For example, Santer et al. (2018, 2022) found an increase in the temperature amplitude globally and throughout the troposphere in recent observations and future projections and attributed it to anthropogenic forcing. Further, Chen et al. (2019) concluded that the CMIP5 global climate models project increased seasonal amplitudes in low-latitude regions and most global ocean areas. In contrast, the seasonal amplitudes are expected to decrease over the Southern Ocean and high-latitude regions. Ruosteenoja et al. (2020) describe projected lengthening of the summer season in Northern Europe and López de la Franca et al. (2013) show the same for Spain, together with the winter season practically disappearing.

Earth system models (ESMs) are state-of-the-art instruments for assessing possible future climate evolutions and attributing observed and projected climate changes to their potential causes. The multi-model ensemble produced under the Coupled Model Intercomparison Project Phase 6 (CMIP6) initiative, coordinated by the World Climate Research Programme's (WCRP) Working Group on Coupled Modelling (Eyring et al., 2016), represents the newest set of ESM simulations. This ensemble includes simulations of a range of different models under several shared socioeconomic pathways (SSPs, Tebaldi et al., 2021,) enabling the analysis of uncertainties arising from structural model differences (Abramowitz et al., 2019). Indeed, despite indisputable progress in the complexity of the newest generation ESMs, many uncertainties and issues still need to be solved (Shaw and Stevens, 2025; Randall et al., 2019; Bordoni et al., 2025). The issue of the choice of ESMs appropriate for climate change scenarios is a very complex task, and different approaches are still under investigation (e.g., McDonnell et al., 2024; Snyder et al., 2024; Merrifield et al., 2023;  Rahimpour Asenjan et al., 2023).

Table 1 Basic information about the two reanalysis datasets used in the present study.

Table 2 Basic information on the CMIP6 ESMs used in the present study. The values of the equilibrium climate sensitivity are taken from Meehl et al. (2020).

The shape of the annual cycle of air temperature differs significantly among different regions around the globe. Therefore, performing a global analysis requires focusing on quantities that can be evaluated for all these different shapes (e.g., single and double waves, different timing of seasons, etc.). A lot of previous studies relied on Fourier-transform-based methods, assuming a sinusoidal shape of the mean annual cycle and focused on its amplitude and phase (e.g., in Paluš et al., 2005; Stine et al., 2009; Zhao et al., 2021; Marvel et al., 2021; Deng and Fu, 2023; Zhang et al., 2025), which in some cases resulted in omitting certain regions (e.g., Dwyer et al., 2012; Yettella and England, 2018) and a large portion of studies focused on the Northern Hemisphere only. Here, we introduce an innovative approach based on Functional Data Analysis (FDA). The evolution of temperature throughout the year is approximated by daily long-term mean temperature values, which are converted to functional form (see Sect. 3 for details). This approach allows us to assess any existing shape of the temperature seasonal cycle. López-Franca et al. (2022) also employed smoothing of daily temperature values with splines and assessed changes in dates of minimum and maximum of the smoothed annual cycle, as well as the dates of minimum and maximum slope changes. Unlike the methodology presented here, they only concentrated on specific parts of the year and on the midlatitude regions. We previously successfully applied a Functional data analysis approach to investigate the influence of the driving global climate model on nested regional climate simulation within a multi-model ensemble (Holtanová et al., 2019).

2 Data

The present study deals with the annual cycle of daily mean near-surface air temperature. To analyze its recent changes over both land and ocean regions, we use two reanalysis datasets from the European Centre for Medium-Range Weather Forecasts (ECMWF), namely the ERA5 (Hersbach et al., 2020) and CERA20C (Laloyaux et al., 2018). The choice was motivated by the long temporal coverage of these datasets back to the 1950s. Some basic information about these datasets is described in Table 1. One of the main differences between them is that ERA5 is an atmospheric reanalysis; in contrast, CERA20C was created using a coupled modeling system with the representation of not only the atmosphere but also the ocean, land, oceanic waves, and sea ice. The atmospheric modeling system (ECMWF's Integrated Forecast System (IFS) version CY41R2) is the same for both ERA5 and CERA20C (Laloyaux et al., 2018; Hersbach et al., 2020). As the coupling demands large computational costs, CERA20C has a coarser horizontal resolution (Table 1). The CERA20C dataset includes 10 members representing the spread related to the errors in the assimilated observations and the modeling system (Laloyaux et al., 2018). We use the “number0” ensemble member and do not analyze the uncertainty spread here.

Table 3 Overview of time periods investigated in the present study.

Further, we select historical and scenario simulations of five CMIP6 ESMs (Table 2). The model choice is motivated by the different values of the equilibrium climate sensitivity (Meehl et al., 2020) and overall good performance compared to the whole CMIP6 ensemble (Bock et al., 2020). We employ only five models to be able to analyze the individual simulated curves of the mean annual cycle and illustrate the innovative methodology properly. For the scenario period, we analyze outputs for the SSP3-7.0 socio-economic pathway, which represents the medium to high end of the whole range of the SSPs currently considered plausible (Tebaldi et al., 2021).

The analysis focuses on the time periods described in Table 3. The two historical periods are used for the assessment of recent observed changes. The difference between the future and reference periods is referred to as the projected or expected future change. For both the reanalyses and ESMs, the long-term mean values of near-surface air temperature for each day of the year are averaged over the reference regions from Working Group 1 of the IPCC AR6  (Iturbide et al. 2020) directly from the native grids. These daily long-term mean values are then subject to functional data analysis as described in the following section. To enable comparison of our results to global mean, annual mean temperature changes, Table S1 in the Supplement lists these changes for both the historical and future time periods, with additional division into land and ocean areas.

3 Functional data analysis approach

3.1 Construction of the functional data

The modeling of the mean seasonal cycle of temperature uses the techniques of Functional Data Analysis (FDA), a relatively novel statistical approach (Ramsay and Silverman, 2005; Horváth and Kokoszka, 2012; Kokoszka and Reimherr, 2017). Unlike traditional statistics, a single observation of a variable is not a data point but rather a function. This approach is especially suitable for a series of observations with an underlying correlation structure.

Figure 1 (a) Basis functions for the case K=5; (b–d) smoothed temperature with respect to K=5, 15, 55. As the number of coefficients grows more and more variability beyond the mean seasonal cycle is captured by the FDA.

In general, the relation between a covariate x and a response Y can be modeled as a function y=f(x) using the pairs (x_i,Y_i) of the data, $i=\mathrm{1},\mathrm{\dots },n$. In our case, the covariate x is represented by the days of the year (x_i varies from 1 to 365, for leap years, values for 29 February were deleted). The mean seasonal cycle of temperature plays the role of response Y. To account for the periodic nature of the data, the function f(x) is defined as a linear combination of Fourier basis functions:

$$\begin{array}{}\text{(1)}& f\left(x\right)=a_{\mathrm{0}}+{\sum }_{n=\mathrm{1}}^m\left(a_n\text{cos}\text{cos}{\displaystyle \frac{n\mathit{\pi }x}{\mathrm{365}}}+b_n\text{sin}\text{sin}{\displaystyle \frac{n\mathit{\pi }x}{\mathrm{365}}}\right)\end{array}$$

i.e., f(x) depends on $K=\mathrm{2}m+\mathrm{1}$ coefficients $\mathit{\{}{a}_{\mathrm{0}},a_{\mathrm{1}},b_{\mathrm{1}},\mathrm{\dots },a_m,b_m\mathit{\}}$ and basis functions (see Fig. 1a for K=5).

The coefficients $\mathit{\{}{a}_{\mathrm{0}},a_{\mathrm{1}},b_{\mathrm{1}},\mathrm{\dots },a_m,b_m\mathit{\}}$ are chosen to minimize the following functional:

$$\begin{array}{}\text{(2)}& {\sum }_{i=\mathrm{1}}^n[Y_i-f(x_i){]}^{\mathrm{2}}\end{array}$$

Our approach fits a function f with varying degrees of freedom to the data. In contrast to a simple interpolation, this does not necessarily mean that the function just connects all adjacent data points (i.e., in all cases where the degrees of freedom are less than the number of data points, see Fig. 1). In general, the particular values, x_i, of the covariate and the corresponding observed responses, Y_i, are linked by $Y_i=f\left(x_i\right)+\mathit{\epsilon }i$, $i=\mathrm{1},\mathrm{\dots },m$, where ε_i are realizations of the random errors. This corresponds to the situation where the covariate, x_i, is given, and the observed response, Y_i, is the realization of some random variable linked with the value of x_i. The resulting function f balances the size of the errors, ε_i, and the smoothness of the function linking the covariate and the response. The smaller the number of basis functions K is, the less sensitive it is to fluctuations in the data – compare Fig. 1b and Fig. 1d for cases K=5 and K=55. Here we choose K=15. The choice is supported by the fact that for this value, the character of the FDA curve best resembles the 30 d running average. The 30 d average is analogous to the monthly mean, and the length of the month is an intuitive choice in climatology, as generally a lot of climatological analysis is based on monthly mean values. Moreover, even for K=5, the FDA function explains more than 99 % of the variance of the 30 year mean values of temperature, even though the curve does not entirely align with the underlying data (Fig. 1b). On the other hand, for higher K, the curve becomes too fluctuating, resembling high inter-daily variability in the data. Therefore, we consider smoothing based on K=15 appropriate for the current study. However, the results of the analysis are not sensitive to the choice of K (not shown). The FDA-smoothed curves of the mean annual cycle for all the datasets and geographical regions are shown in Figs. S3 and S7 in the Supplement for the historical periods, and in Figs. S4 and S8 in the Supplement for the projections.

Figure 2 The FDA diagnostics interpretation framework. Blue and red lines illustrate an example FDA-smoothed temperature annual cycle in two time periods, except for panel (c), where the lines represent the absolute values of the 1st derivative of the FDA-smoothed curves, i.e., temperature velocity, and the dashed color lines are used to depict negative temperature velocities. (a) The black arrow corresponds to the vertical temperature change between the two periods on a specific day of the year, and the black line represents its values during the whole year. (b) Dashed lines represent the days of temperature maxima, and the black arrow corresponds to the shift of the maximum. (c) The black arrow corresponds to the change of temperature velocity between the two periods on a specific day, and the black line represents its values during the whole year. (d) The blue and red arrows correspond to the amplitudes in each period, and the vertical black arrow illustrates the change in the amplitude between the two periods.

Table 4 Overview of the FDA diagnostics and the corresponding figures showing the results. The “Regional figures” in the Supplement show the FDA results underlying the individual diagnostics.

3.2 FDA diagnostics

Drawing on the FDA representation of the annual cycle in each of the time periods specified in Table 3, we now define diagnostics that evaluate changes in the shape of the annual cycle (Sects. 3.2.1–3.2.4). Figure 2 illustrates the interpretation of the diagnostics on example data. Table 4 provides an overview of the diagnostics and references to the figures showing the results based on them. All the diagnostics are further used to quantify differences between the annual cycle curves in the two historical periods and between the future and reference periods. For the historical periods, we compare the GCMs with ERA5 and CERA20C. For future time periods, we compare individual GCMs with their multi-model mean. We want to emphasize that in the projections, the multimodel mean values are based on the multimodel mean annual cycle, not the multimodel mean of the FDA diagnostics. Therefore, the multimodel mean values of FDA diagnostics can fall outside the range of individual ESMs. Regarding the terminology, we note that we use both terms “seasonal cycle” and “annual cycle” intermittently in the text, with no difference in its meaning.

3.2.1 Annualy integrated temperature change

For each day of the year, we calculate the distances between the smoothed annual cycle curves (see Fig. 2a). We then aggregate these distances in three ways: by calculating the 10th and 90th percentiles and the root mean square of them. The former two diagnostics, hence, represent high and low annual extremes of the temperature changes (allowing both positive and negative values), while the latter diagnostic evaluates the Euclidean distance of the whole annual cycles (positive by definition). Figures S1 and S2 in the Supplement show the occurrence of values below/above 10th/90th percentiles; the red dashed line represents the 10th percentile, and the green dashed line represents the 90th percentile. For all values below/above the 10th/90th percentile threshold, the time periods of the year when these values occur are shown in red/green.

3.2.2 Annual cycle maximum

We define the shift in the annual cycle maximum as the number of days between the maximum of the annual cycle in two periods (black arrow in Fig. 2b). Positive values indicate a delay in the maximum occurrence relative to the reference period, and negative values vice versa. In regions with two (local) maxima in the annual cycle, the “first” and “second” maxima are considered chronologically from 1 January, with no regard to the actual maximum magnitude (the second maximum can potentially have a higher temperature than the first). There are nine regions, where we identify two distinct maxima, see Figs. 5, 6, S3, and S4.

3.2.3 Annual cycle velocity

We calculate 1st derivative of the smoothed curve of the temperature annual cycle. We define temperature velocity as the absolute value of this 1st derivative curve. It gives an indication of the steepness of the annual cycle on individual days (see Fig. 2c). Then we calculate changes in temperature velocities between corresponding days of the year between the two time periods. Positive differences in temperature velocity, hence, indicate days where the annual cycle is getting steeper compared to the reference period, and vice versa for negative values. Note that steepness means faster warming as well as faster cooling because we consider absolute values of the 1st derivative. Similarly to the changes in temperature itself (Sect. 3.2.1), we calculate the 10th and 90th percentiles of the differences and their root mean square.

3.2.4 Annual cycle amplitude

The amplitude of the annual cycle is defined as the difference between the maximum and the minimum value in °C (see Fig. 2d). Here we evaluate the change in the amplitude between two time periods. Consequently, a positive change in the amplitude indicates an increasing temperature range over the year compared to the reference.

Figure 3 (a) Root mean square difference (Euclidean distance, [K]) of the whole FDA-smoothed mean annual cycle curves between the two historical periods 1981–2010 and 1951–1980. (b) 10th percentile and (c) 90th percentile of daily distances [K] between the smoothed annual cycle curves. For each region, the center of the pie plot shows results based on the two reanalysis datasets ERA5 and CERA20C., while the outer part of the pie shows the results for the five CMIP6 ESMs (see Sect. 2 for data description).

4 Results

Here, we discuss changes in the four diagnostics from a high-level perspective; detailed figures for each of the regions can be found in the Supplement (see Table 4 for an overview).

4.1 Annualy integrated changes in the annual cycle

In the often employed annual-mean view, warming is evident almost everywhere on the globe, with land areas and higher latitudes generally warming faster (Gulev et al., 2021). Our approach resolves seasonal differences in the long-term warming signal and reveals that differences in the historical period between 1951–1980 and 1981–2010 can be negative during (short) parts of the year in many regions (Fig. 3b). This is compensated by a strong warming in other parts of the year, which can exceed 2 °C in many Northern Hemisphere land regions (Fig. 3c).

Figure 3a shows the resulting annualy integrated temperature change as the root-mean-square of the daily differences (RMSD, also termed Euclidean distance), which also exceed 1.5 °C in most datasets at northern mid-latitude land regions. In most other parts of the world, except Antarctica, the RMSD remains lower than 1.5 °C for the historical periods. We stress that this diagnostic embraces both negative and positive temperature changes, evaluating the overall change in temperature , unlike simply averaging the changes over the year.

With regard to warming between the historical periods, in general, a stronger signal is seen in the Northern Hemisphere than in the Southern Hemisphere. An exception is the strong warming signal in CERA20C in Antarctica (Fig. 3). Larger disagreement between the reanalyses also occurs over the southern ocean and in some regions near the equator (e.g., SAH and ARP).

Figure 4 Same as Fig. 3, but for the differences between the scenario period 2071–2100 and the reference period 1961-1990. Future model simulations follow the SSP3-7.0 socio-economic pathway.

With regard to the timing, both reanalyses show the highest temperature increase during northern-hemisphere winter, or the changes do not have any distinct maximum/minimum (Fig. S1). Only in New Zealand (NZ), the Southern Ocean (SOO), and Antarctica (WAN and EAN) are the changes larger in the southern-hemisphere winter. In the ESMs, the timing of the largest increase/decrease often does not match the reanalyses  (e.g., over Greenland, the reanalyses show a decrease of temperature in the first three months of the year, whereas the ESMs show the decrease (if any) later in the year, Fig. S1). In the Arctic region (ARO), the ESMs and reanalyses generally agree that the lowest increase in temperature occurs in summer. The timing of the changes projected for the end of the century shows a distinct difference between the regions in the northern middle latitudes and subtropical areas. In the former (e.g., NEU, EEU, NWN, NEN), the highest changes occur during winter, and the latter (e.g. MED, WCA, CNA, ENA) during late summer or autumn (see Fig. S2).

For the warming at the end of the 21st century, the Arctic stands out with temperature increase exceeding 10 °C in all models during the 10 % of strongest warming days (Fig. 4c). Such stronger warming in the polar regions compared to lower latitudes (often referred to as polar amplification) is consistent with theoretical considerations and historical observations (e.g., Stuecker et al., 2018; Previdi et al., 2021). Here, we show that the stronger warming at high latitudes predominantly comes from the upper end of the annual temperature distribution, with the 10th percentile of changes being mostly uniform across latitudes (Figs. 3b and 4b).

Polar amplification has also been reported to be underestimated in CMIP6 models (Casado et al., 2023) and to be weaker in Antarctica than in the Arctic region in both observations and CMIP6 models (e.g., Zhang et al., 2023; Xie et al., 2022). Our results contradict these results to some extent. Mainly, the five ESMs simulate the magnitude of historical warming in the Arctic, higher or comparable to the reanalyses (Fig. 3). Finally, we note that in the historical period, one of the two observation-based reanalyses, CERA20C, shows stronger warming in Antarctica than in the Arctic. This discrepancy might be attributable to high decadal variability in Antarctica (Casado et al., 2023) and large uncertainties of the reanalysis outputs over this remote region with low density of assimilated observations (Laloyaux et al., 2018).

Figure 5 Similar to Fig. 3a, but for (a) shift in the annual cycle maximum and (b) shift of the second maximum in regions with two distinct maxima. Note that the “first” and “second” maxima are considered chronologically from 1 January, with no regard to the actual maximum magnitude (the second maximum can potentially have a higher temperature than the first).

Figure 6 Same as Fig. 5, but for the shifts of the maxima between the scenario and historical periods.

4.2 Shift of the annual cycle maximum

For the end of the 21st century, the five selected ESMs project a delayed maximum near the poles and an earlier maximum in many tropical continental regions (Fig. 6). The shift of the maximum between the two historical periods does not show such distinct pattern (Fig. 5). In Europe, the southern and eastern regions experienced a delay of maximum of up to 10 d, wherease northern Europe a slightly earlier maximum. In North America, a similar dipole pattern of historical changes is seen between eastern and western regions (Fig. 5).

The largest differences between the reanalyses and ESMs are found over southern America and eastern and southern Africa (Fig. 5). Also, in the land regions near the equator, there is a disagreement between the two reanalyses (e.g., WSAF, MDG, ESAF in Africa, and NWS and SAM in South America). This is mainly due to the fact that the annual cycle has no distinct maximum peak and the warm season part of the annual cycle is rather flat; thus, a small temperature change in this season may result in a large shift of the maximum (Fig. S3). In North-Central America (NCA), even though it is farther from the equator, CERA20C gives a large shift of the maximum, but the actual temperature change is small, similar to regions NWS and SAM in South America, which are closer to the equator. In most of the regions further from the equator, the reanalyses agree on the sign of the shift in the maximum. In the oceanic regions near the equator, the reanalyses show a shift to an earlier onset of the maximum. Both between the two historical periods and between the reference and future period, we see a smaller shift of the maximum in Antarctica than in the Arctic (Figs. 5 and 6).

In the regions near the equator, there are two distinct maxima of the annual cycle (Figs. S3 and S4); therefore, we evaluate the shift also for the second maximum (Fig. 5). We stress that the first/second refers to the earlier/later occurrence during the year, not to the magnitude. In some of these regions, the annual cycle has even more “maxima”; it is modulated by at least three peaks (Fig. S3, e.g., CNRM-ESM2 in north-eastern Africa, NEAF). As the amplitude of the annual cycle is generally low in near-equator regions, the whole curves are rather flat, and it is difficult to compare them between the datasets. For example, in the oceanic part of south-eastern Asia (SEA region), the CERA20C reanalysis shows a large shift in the 2nd maximum (Fig. 5). However, in Fig. S3 it is clear that in the first historical period, the annual cycle near the 2nd maximum is very flat, and therefore the large shift rather indicates a clearer emergence of the 2nd maximum. Also, in north-eastern Africa (NEAF), the evaluation of the maximum shift is rather problematic. The maxima in different ESMs and reanalyses are shifted, so it is actually questionable to compare them (Fig. S3). Similarly, in north-west southern America (NWS), the mean annual cycle in the historical periods has, according to the reanalyses, only one distinct maximum (Fig. S3). However, the ESMs show a second maximum. We do not consider it in our analysis, but it is interesting to note that the annual cycle in this region is projected to change in the way that the temperature at this second maximum, not present in reanalyses, becomes higher than the first maximum (Fig. S4). As the (radiation-driven) annual cycle in the near-equator regions is less distinct, other climate system processes, such as the distribution of precipitation, become more important in shaping it. As a result, the shift of the temperature maximum can be indicative of a change in the occurrence of dry and wet seasons at low latitudes. At the same time, we note that in above mentioned regions with rather flat maximum and low amplitude, the ESMs and reanalyses mostly disagree on changes in amplitude, so that overall confidence in the signals is rather low.

Figure 7 Same as Fig. 3, but for temperature velocity, that is the 1st derivative of the smoothed curve of the temperature annual cycle. Note that the color scale for plot in (b) has the same range as plot (c), just reversed, i.e., negative values going from 0 K d⁻¹ (white color) to −0.08 K d⁻¹ (darkest blue). The range is not shown for the sake of better visibility of the plots.

Figure 8 Same as Fig. 7, but for the change in temperature velocity between the scenario and reference periods.

4.3 Annual cycle velocity

Higher latitudes detect a higher magnitude of temperature velocity change than lower latitudes. For future changes, the geographical pattern of projected temperature velocity change remains the same as between the historical periods, with the magnitude of change twice as high in most of the regions (Fig. 8). The velocity change over the oceans is mostly smaller compared to the continents (Figs. 7 and 8). Between the two historical periods, all regions experienced both a decrease and an increase in the slope of the annual cycle, depending on the time of year (see Figs. 7b, c and S5 in the Supplement). Historical changes in temperature velocity are largest in the western-central part of Euroasia (EEU, WSB, and ESB; Fig. 7). Generally, for both historical and projected changes, the regions with larger changes in velocity have a larger range between the 10th and 90th percentiles, which is expected given the definition of the diagnostic.

The temperature velocity changes agree between the reanalyses, except for Antarctica and the RFE (eastern Asia), and CNA (central North America) regions. Further, the ESMs tend to underestimate the reanalysis-based velocity changes in the middle and higher latitudes of the Northern Hemisphere, and largely agree on the smaller changes in the tropics and over the Southern Hemisphere. The temperature velocity changes between the two historical periods are mostly in the interval between −0.1 and +0.1 °C d⁻¹. This change of slope of the annual cycle curve can result in a temperature change of 3 °C per month. It naturally corresponds to changes in the amplitude of the annual cycle and changes in the temperature contrasts between seasons. However, Fig. S5 shows that the changes in the velocity are, in most cases, rather variable during the year, with the sign persisting not for the whole season, but rather for a week up to two months. Still, except for a couple of regions, the two reanalyses have rather similar annual cycle of temperature velocity changes in the historical periods. Unlike the temperature change, the reanalyses agree on the sign and value of annual cycle velocity change over Greenland.

The five ESMs mostly follow the reanalysis-based pattern of change in temperature velocity. If there is any disagreement, the models tend to underestimate the magnitude of changes. This is mainly seen in the Northern Hemisphere's higher latitudes. Generally, over the Northern Hemisphere continents, we mostly see higher fluctuation of velocity changes between negative and positive values in winter than in summer. Regarding the projections, Fig. S6 in the Supplement depicts a distinct annual cycle of velocity changes in the regions where the expected warming is larger in one of the seasons. A nice example is the Arctic, where we can see a decrease in velocity in the spring and the autumn, but near-zero or positive changes in winter and summer (Fig. S6). This stems from projected flattening of the annual cycle and higher warming in winter than in summer. In northern mid-latitude regions, the velocity changes are more variable during winter than in summer (Fig. S6), which is connected to a higher increase of temperature in winter than in summer and shrinking amplitude, as discussed below.

Figure 9 Same as Fig. 3a, but for the change in the amplitude of the annual cycle.

Figure 10 Same as Fig. 4a, but for the change in the amplitude of the annual cycle.

4.4 Annual cycle amplitude

Both reanalyses agree on the prevailing decrease in amplitude, with the largest changes detected in EEU (eastern Europe) and WSB (western Siberia), CNA, NWN (both northern America), the Arctic ocean, and Antarctica (Fig. 9). Unlike the temperature change, both reanalyses show that the amplitude change in the Arctic is larger than in Antarctica. The ESMs, in turn, show diverging changes of amplitude between the historical periods over the globe. Over the oceans, we see small amplitude changes with varying signs in both models and reanalyses, except for the Southern Ocean and the Arctic region, where all the models and reanalyses show an amplitude decrease (ERA5 −0.4 and CERA20C −1.1 °C in SOO, ERA5 −2.9, CERA20C −1.6 °C). Clearly, decreasing annual cycle amplitude arises from a faster increase of temperature in the colder part of the year, in comparison to the warm season (Fig. S7). Such seasonal difference in temperature trends during the 20th century has been reported by Nigam et al. (2017) for the Northern Hemisphere. However, there are several regions where the reanalyses show an increase in amplitude, up to 0.8 °C, e.g., Greenland, the southern part of South America, Madagascar, and interestingly also Siberia (RAR) and some regions in northern America (Fig. 9). The ESMs show decreasing amplitude everywhere. We hypothesize that the discrepancy between ESMs and reanalyses over Greenland could be connected to differences in the evolution of sea ice between simulations. A recent increase in amplitude over Greenland has also been reported by Deng and Fu (2023). Liu et al. (2024) showed an increase in the amplitude of SSTs over most of the ocean basins in recent 40 years. Our study period is longer, and the amplitude increase attributed to the anthropogenic forcings is probably masked to some extent by internal climate variability.

In many regions, the projected future amplitude changes have the sign opposite to the changes between the historical periods (compare Fig. 9 and Fig. 10). The largest increase is projected over the Mediterranean and West-Central Asia regions (due to summer warming being more than winter), and the largest decrease over the Arctic Ocean (due to winter warming being higher in winter than summer). Generally, the amplitude increase is projected for most of the Southern Hemisphere and equatorial areas, whereas most of the middle to higher latitudes of the Northern Hemisphere are projected to experience a decrease in the amplitude. There are a few exceptions: West-Central Europe (WCE), East-Central Asia (ECA), and the western part of the USA (WNA), where we can see an increase in the amplitude of app. 2–3 °C. Over the continents, the increasing amplitude indicates an increase in thermal continentality of climate, with higher contrasts between winter and summer.

To illustrate how the individual FDA diagnostics complement each other, the projected amplitude changes correspond very well to the projected temperature velocity changes; a higher increase in velocity is connected to a higher decrease in amplitude, and the other way around.

5 Discussion

The shape of the mean temperature annual cycle can be considered a very basic feature of climate. Nonetheless, we highlight large observational uncertainty related to its recent changes, i.e., distinct differences between the two reanalyses, especially over equatorial and polar regions. Multiple differences between the reanalyses might be behind the discrepancies. Besides differences in spatial resolution, CERA20C is a coupled reanalysis, whereas ERA5 was produced by the atmospheric model only. Laloyaux et al. (2018) emphasize that the former is expected to be more realistic in terms of ocean heat balance and ocean heat uptake, important for the temporal evolution of near-surface air temperature and its low-frequency variability. Regarding the discrepancies over Antarctica, the assimilated observations are scarce and might be spurious (Laloyaux et al., 2018). Moreover, the number of observations in both CERA20C and ERA5 were increasing during the study period, which might have influenced the results. In case of CERA20C, in which only variables measured over the ocean are assimilated, the data inputs from ships more than doubled, and data from buoys started to be assimilated after 1970 (Laloyaux et al., 2018). For ERA5, the number of assimilated observations increased from 53 000 to 570 000 between 1950 and 1970 (Bell et al., 2021). Furthermore, as pointed out by, e.g., McKinnon et al. (2024), the reanalysis performance is in general questionable over regions that have spurious observations, not only in Antarctica but also over large portions of Africa or Southern America. Furthermore, Yettella and England (2018) emphasized large internal climate variability uncertainty connected to the evolution of annual cycle shape over Northern Hemisphere middle and high latitudes. Brunner et al. (2025) emphasize that the discrepancies between different reanalyses are rather larger in the southern ocean than in other ocean basins, not only for ERA5 and CERA20C, but also for other reanalyses. Moreover, the big difference between CERA20C and ERA5 in southern high latitudes, unlike in northern high latitudes, points to the importance of the coupling between the atmosphere ocean for the southern high latitudes. As demonstrated by, e.g., Kang et al. (2023), there is a strong relationship between tropical and subtropical Pacific and temperature changes in the southern ocean, and the simulation of these features is expected to be different in ERA5 (atmosphere only) than CERA20C (coupled simulation). The coupling does not automatically guarantee a better simulation, naturally.

Even though we analyze only five CMIP6 ESMs, which is admittedly a very small subset of the whole multi-model ensemble, they differ in many aspects, including spatial resolution (Table 2), model family (Merrifield et al., 2023), and climate sensitivity (Table 2, Meehl et al., 2020). It is not thus surprising that they show diverse outcomes. They are not always able to reproduce the reanalysis-based historical changes, and their projections differ in many aspects. The differences in the structure of the models imply different character of internal climate variability, which is certainly behind some of the discrepancies. ESMs with higher climate sensitivity (corresponding to higher globally averaged temperature changes as listed in Table S1) generally project larger annual cycle shape changes (e.g., CanESM5 in the Arctic, Figs. 3 and 4). Even though it has been argued that the higher sensitivities are not plausible (e.g., IPCC, 2021), it is difficult to rule out the hot models, especially in the case of regional impact assessment (Palmer et al., 2023; Swaminathan et al., 2024). This is illustrated in our results by cases where the regional changes in FDA diagnostics do not correspond to the differences in global mean temperature changes. For example, these cases include West Antarctica and regions in south-eastern north America.

Previous studies on changes in the annual cycle mostly concentrated on the amplitude and shift of the maximum. Chen et al. (2019) studied ERA-Interim-based and CMIP5-simulated spatial patterns of seasonal amplitude and phase. They concluded that the seasonal amplitude reduced during the 21st century over high latitudes of both hemispheres because cold-season air temperature increases faster than warm-season air temperature. In contrast, over low latitudes, the expected evolution is exactly the reverse. Further, the maximum of the annual cycle was projected to be delayed by 15–30 d over the high-latitude oceans where the sea ice is expected to shrink significantly (Chen et al., 2019). All these patterns are also obvious in our results, implying consistency between CMIP5 and CMIP6 projections. The gradual decrease of amplitude prevailing over most of the Northern Hemisphere has also been reported in other studies, including Stine and Huybers (2012), Wang and Dillon (2014), Nigam et al. (2017), and Cornes et al. (2017). However, we depict some regions where the reanalyses disagree on the sign of change, and also regions where both ESMs and reanalyses imply an increase in the amplitude. Delayed onset of annual cycle maximum over most of the northern-hemisphere continents was also reported by Deng and Fu (2023).

It is generally expected and observed that the temperature changes over land would be higher than over the ocean (e.g., Sutton et al., 2007, see also Table S1). Our results confirm this expectation in most of cases, especially in middle latitudes, when comparing regions for large ocean basins and surrounding regions (e.g., Figs. 3 and 4), and when comparing the oceanic and land parts of the marginal sea regions (e.g., Mediterranean, Caribbean, south-east Asia). Near the equator, especially in the eastern Pacific Ocean (EPO) region, the position of the annual cycle maximum is expected to change more in the future than over continental regions in South America. This is also a region with relatively high disagreement between individual models. The uncertainty is apparently connected to rather flat maximum and low amplitude of the seasonal cycle, with even a small temperature change in a part of the year potentially causing a large shift of the maximum. As already discussed, this might also be connected to changes in precipitation distribution over the seasons.

In a recent study, Brunner and Voigt (2024) revealed a systematic bias in the definition of percentile-based temperature extremes (Tx90p) when using too long seasonally running windows. One of the pitfalls they revealed was spurious signals of change in Tx90p, as the strength of the bias depends on the shape of the temperature annual cycle (as well as the day-to-day variability). They find two particularly affected regions: a region of increasing bias in oceans north of 45° N, except the very highest latitudes (approximately our NPO, NAO, and MED regions), and a region of decreasing bias in our ARS region (see their Fig. 5a). Connecting to our results, stronger seasonal gradients (corresponding to a higher temperature velocity) favour a stronger bias in Tx90p (Brunner and Voigt, 2024). Indeed, we find a weak (in particular compared to some land regions; see Fig. 8), but clear increase in temperature velocity between 1961–1990 and 2071–2100 in all three regions affected by increasing bias (NPO, NAO, MED), which is particularly pronounced towards and away from the annual maximum (which is located around end of August/beginning of September or day of the year 240, Fig. S4). While the absolute value of the temperature velocity change is considerably higher in other regions, its systematic increase in combination with the low day-to-day temperature variability considerably contributes to the increase in Tx90p bias in these regions.

For the region of decreasing bias in Brunner and Voigt (2024), roughly corresponding to our ARS region, the attribution of the bias change to the temperature velocity is less clear due to a combination of two reasons; first, the decreased bias stems mainly from a very limited number of days surrounding the second annual minimum in the region (July and September; see Fig. 5c in Brunner and Voigt, 2024). For CanESM5, which was used in Brunner and Voigt (2024), we do find a short consistent decrease in the temperature velocity corresponding to those months (Fig. S6). Second, the decrease in bias found in this region is also (at least partly) attributable to an increase in day-to-day variability, which is not evaluated by the FDA diagnostics.

6 Conclusions

This paper presented an innovative method for assessing the shape of the annual cycle. It is applicable for different climatic variables and for various purposes, with no need to make any prior assumptions about the annual cycle shape. The diagnostics we introduced provide important information about different aspects of the seasonal cycle shape and its changes: amplitude, slope, and location of extrema. We analyze annual cycles averaged over 30 year periods, but the method can also serve for analysis of shorter-term variability of seasonal cycle, to study even inter-annual variability of the shape features. Unlike methods based on monthly or seasonal means (e.g., evaluating the amplitude based on monthly values), the FDA diagnostics can capture even slight changes in the shape of the annual cycle, for example, in the timing of the maximum, discussed in Sect. 4.2.

We have illustrated the methodology on the example of the temperature annual cycle and its changes in pre-defined climatological regions. We used it to assess recent and projected changes of the temperature annual cycle in a selection of ESMs and observation-based datasets. Other potential applications include assessing other variables, for example, minimum or maximum air temperature. The changes in the shape of annual cycle of these variables can have implications for the occurrence of extreme cold or heat events. Another possibility is to apply our method for explicitly evaluating the model performance. In that case, differences between models and one or more reference datasets would be investigated. The results can be aggregated, resulting in assessments of the ensemble mean and spread. The diagnostics can be modified to evaluate not only changes between time periods, but also differences between datasets to reveal model biases in the representation of the annual cycle compared to an observational reference. If the evaluation is applied to the outputs for regional climate models or ESMs with finer resolution, the assessment can be done for smaller geographical regions, revealing more details about projected changes and their potential impacts. The definition of FDA diagnostics can thus be tailored for specific interests and applications.