Characterization of seasons over the extratropics based on the annual daily mean temperature cycle

A proposal to characterize seasons based on the annual cycle of daily mean temperature over extratropical regions is presented. Four metrics are computed, based on the dates along the year when the maximum and minimum temperature values and the maximum warming and cooling (when the rate of temperatura change is maximum or minimum) are obtained. These four metrics are computed for 1980–2014 period for ERA‐Interim, JRA‐55 and NCEP2 reanalysis products. Results indicate that these four metrics are able to represent, with geographical coherence over any extratropical region, the start and end of subperiods along the year and that can be seen and used as a consistent procedure to define seasons. The two more extreme metrics (dates for the maximum and minimum temperature) are closely influenced by the variability of the net radiation and near‐surface circulation. Regions with different climatic regimes, defined in terms of land–atmosphere coupling conditions, have different contributing factors influencing the interannual variability of the metrics. In humid regions the net surface radiation is critical for the timing of the maximum temperature date, while the variability of the atmospheric circulation is relevant for the minimum temperature metric. In transitional regions, land–atmosphere interaction plays a key role in the timing of the annual cycle, both in summer and winter. Variations in net radiation are relevant for dry regions. Therefore, a robust methodology to establish seasons based on the annual temperature cycle is presented. We also inspect how the metrics timing responds in terms of the interannual variations of different relevant atmospheric variables.

In humid regions the net surface radiation is critical for the timing of the maximum temperature date, while the variability of the atmospheric circulation is relevant for the minimum temperature metric. In transitional regions, landatmosphere interaction plays a key role in the timing of the annual cycle, both in summer and winter. Variations in net radiation are relevant for dry regions.
Therefore, a robust methodology to establish seasons based on the annual temperature cycle is presented. We also inspect how the metrics timing responds in terms of the interannual variations of different relevant atmospheric variables.

K E Y W O R D S
climate processes, extratropics, flexible threshold definition, reanalyses, seasons, temperature annual cycle

| INTRODUCTION
Temperature annual cycle is one of clearest variability patterns of land-atmosphere system that controls phenological and climatic processes over most of the world areas, in particular over extratropical regions. From a phenological perspective, temperature is considered as one of the main factors that determine vegetal growing season (Jeong et al., 2011;Piao et al., 2015;Tan et al., 2015). This fact determines the dynamics of natural ecosystems and their services (Fisher et al., 2009;Asseng et al., 2011;Brouček et al., 2011;Iizumi and Ramankutty, 2015). From a climatic perspective, temperature variability encompasses several peak frequencies, from the diurnal to the Milankovitch cycles (Milankovitch, 1941). The annual temperature cycle is dominant in extratropical regions, influencing other meteorological variables like pressure, wind, that is, circulation patterns (Peixoto and Oort, 1992;Wallace and Osborn, 2002), and therefore affecting continental climates both at regional and local scales. There are two main methods to define seasonality based on temperature: by means of the sinusoidal properties of the annual cycle or through thresholds. The first one is able to capture phase changes and annual cycle amplitude at local scale, allowing for regions comparison (Eliseev and Mokhov, 2003;McKinnon et al., 2013). It can show cold or hot periods, but it is unable to distinguish the start and end of seasons, and so their relevance to understand subseasonal temperature changes and their impacts on climatic and biological systems. Characterizing the seasons separately would allow analysis of the climatic processes involved in their occurrence and understanding the year-to-year variations related to the advance or delay of the different seasons and their impact on climatic and biological systems. On the other side, threshold-based methods allow for the identification of yearly subperiods, and so they are commonly proposed to define seasons both climatically (Jaagus and Ahas, 2000;Dong et al., 2010;L opez-Franca et al., 2013) or phenologically (Zhu et al., 2012;Wang et al., 2014).
Besides describing different features of the temperature annual cycle (as it can be the case of heating or cooling subperiods), it is of high importance to know which processes and mechanisms control those seasonal characteristics. Incoming solar radiation at surface is clearly the main temperature modulator along the year. The amount of solar energy arriving at surface varies during the annual cycle and at each latitude, leading to the astronomical seasons. This cycle is more pronounced over the extratropical regions (Ahrens and Henson, 2016).
Extreme maximum/minimum insolation days (solstices) do not coincide with days of maximum/minimum temperature, being delayed around 30 days on average (Trenberth, 1983;Douglass et al., 2004). The magnitude of this delay is locally related to surface characteristics, microclimates, the season of the year, degree of continentality, or orographic features (Seco et al., 1993). Coastal regions exhibit larger delays than continental ones (Prescott and Collins, 1951;El-Hussainy and Essa, 1997;Donohoe et al., 2020).
On regional scales, surface properties and atmospheric circulation seem to be among the main factors that modulate annual temperature cycle (Stine et al., 2009;Ahrens and Henson, 2016;Nigam et al., 2017). Net surface energy is closely related to sensible (SHF) and latent (LHF) heat fluxes. SHF, proportional to wind and vertical temperature gradient near the surface, is typically positive during daytime. Meanwhile, LHF is related to evaporation or condensation of surface water so the heat transfer is then indirect, through subsequent phase changes of the water in the atmosphere (Peixoto and Oort, 1992). Both fluxes are related through Bowen ratio (BW = SHF/LHF), that is usually employed for energy budgets for ecosystems (Bonan, 2016;Morwal et al., 2017;Ping et al., 2018).
Moreover, both energy surface fluxes are related to atmospheric conditions and soil moisture (Bonan, 2016;May, 2017) for intra-annual and interannual scales (Cho et al., 2012;Tang et al., 2014;Schwingshackl et al., 2017;Ping et al., 2018). The role of soil moisture is especially relevant over land-atmosphere coupled regions, where transition between dry and humid climates takes place (Seneviratne et al., 2006;. Over those regions, negative soil moisture anomalies are related to temperature increases through a negative evapotranspiration anomaly, as SHF to LHF fraction is increased. On the other side, positive soil moisture anomalies led to a surface cooling due to LHF increases. Noncoupled land-atmosphere regions are those with high soil water supply. There, net radiation (NR) is the limiting factor of LHF, and so atmospheric conditions impact on soil moisture (Schwingshackl et al., 2017). It must be noticed that evapotranspiration, and consequently LHF, is the only variable that is both present on the water and energy budget equations (Seneviratne et al., 2010).
Main global circulation modes of variability modulate the annual temperature cycle on global and regional interannual scales (Stine et al., 2009). Studies over Northern Hemisphere (NH) extratropical regions point to the influence of the Arctic Annular Oscillation (AO), the North Atlantic Oscillation (NAO), and the Pacific-North American teleconnection pattern (PNA) over the temperature field (Leathers and Palecki, 1992;Hurrell, 1995;Luterbacher et al., 2004;Stine et al., 2009;Kryzhov and Gorelits, 2015). Over the Southern Hemisphere (SH), studies mainly point to the effect of the Southern Annular Mode (SAM) over temperature (Thompson and Wallace, 2000;Thompson and Solomon, 2002;Gupta and England, 2006;Menéndez and Carril, 2010;Fierro and Leslie, 2014).
The complexity of the involved processes on the characterization of continental temperature seasonality, together with the uneven geographical availability of quality observational data over long periods, especially those related to surface energy budgets, make difficult the analysis based on observations (Nigam et al., 2017). Thus, databases coming from numerical models can be quite helpful, not only to understand global processes, due to their physical consistency, but also to be able to analyse regions with scarce data. Reanalyses are numerical systems based on assimilation of meteorological data and forecast modelling over regular grids, developed with several parameterizations and spatial and time resolutions. Ensembles of climate models or reanalyses seem to behave better than any individual model/reanalysis as a tool to represent near-surface air temperature and land-atmosphere interaction (Dirmeyer et al., 2006;Zaninelli et al., 2015).

| OBJECTIVE
The main purpose is to present a novel and comprehensive global methodology to quantify seasonal features of extratropical continental zones (over both hemispheres from the Tropics to 60 ) based on the mathematical analysis of the daily temperature annual cycle. Four seasonal metrics are defined, two based on the dates of occurrence of both extreme maximum and minimum mean daily temperature along the year, and two related to when temperature changes are larger, that is, the maximum warming and cooling dates. The dates related to both extremes of temperature, that are likely to be related to main climatic processes will be analysed, with the aim to add understanding to the obtained results from the annual temperature cycle phases (as amplitude of the cycle is not taken into account here). Therefore, seasonal averages of the more relevant surface temperature climate modulators (net surface radiation, sensible and latent heat fluxes, and 10 m wind, to take into account of atmospheric circulation) will be studied with that purpose. The following section describes the data and methods applied to define and analyse the proposed seasonal metrics. Then, results section starts with the main results related to the four proposed metrics for the whole extratropical areas of the globe. Then, mechanisms related to the two extreme temperature dates are inspected. Finally, discussion and conclusions are presented.

| Data
Daily mean air temperature (T2M hereafter) is used to define the four seasonal metrics of the annual cycle of temperature. Net radiation at surface, sensible and latent heat fluxes, zonal and meridional wind components at 10 m and total precipitation were the variables chosen to analyse mechanisms involved in the occurrence of these metrics. All variables were taken from three global reanalysis: ERA-Interim (Dee et al., 2011), JRA-55 (Kobayashi et al., 2015) and NCEP2 (Kanamitsu et al., 2002). The time period is the common one to all of them .

| Definition of metrics
A cubic spline function of grade D = 3 at each grid cell (Maindonald and Braun, 2010) has been used to smooth the natural variability of T2M, with the aim to capture, at the same time, as much as possible of the structure of the annual cycle (Dong et al., 2010;Buitenwerf et al., 2015;Donohoe et al., 2020). The function smooth. spline() from stats R package (R Core Team, 2016) was employed for this task. Each year was split into four splines with equidistant knots. Figure 1 shows an idealized description of the proposed metrics: the day of the year of maximum (DXT) and minimum (DNT) T2M, and the day of maximum warming (DXWARM) and cooling (DXCOOL). A similar smoothing procedure was proposed in Donohoe et al. (2020).
These last two metrics were computed as those days on which the instantaneous rate of change (first derivative respect to time) of T2M at time t from the spline function is maximized and minimized, taking into account leap and nonleap years:

| Statistical analysis
Each metric (DXT, DNT, DXWARM, DXCOOL) was calculated at each pixel from the daily T2M values, for each year and for each of the reanalyses covering the period 1981-2013. Daily data are used, as they is more accurate than the approach of using monthly data seen in several previous works (Stine et al., 2009;Dwyer et al., 2012;Stine and Huybers, 2012). First and last years were removed from computations to avoid unrealistic or spurious behaviour of the spline function at their extreme control points. The mean value of each metric was calculated by averaging the annual values on the original grid of each reanalysis. Subsequently, the multireanalysis ensemble mean (ENS) was computed, interpolating to a common 0.5 grid. Previous studies indicate that a multimodel ensemble mean instead of individual models usually shows a better performance for any analysed variable (Hagedorn et al., 2005;Carril et al., 2012;Flato et al., 2013), properly representing the feedback between land and atmosphere (Dirmeyer et al., 2006). As the annual cycle of temperature is strongly and primarily related with the radiative forcing (Prescott and Collins, 1951;Stine et al., 2009;Donohoe and Battisti, 2013;Ahrens and Henson, 2016), metrics will be expressed as days of delay with respect to the start of the corresponding astronomical season at each hemisphere (Song et al., 2009). Thus, DXT is related to June 21 (December 21) in the Northern (Southern) Hemisphere while DXWARM is related to March 20 (September 22). The opposite is in the case of DNT and DXCOOL, respectively. This idea has been proposed by the authors (L opez-Franca et al., 2016; Sanchez et al., 2018), and also recently seen in Donohoe et al. (2020). Moreover, because the reanalyses are affected by temporal and spatial inhomogeneities depending on the assimilated observations and also on errors from different sources, the dispersion of the metrics among the three reanalysis was computed in order to analyse the uncertainties between them as max(ENS)min(ENS) for each seasonal metric. These results are presented as Supporting Information.

| Processes related to DXT and DNT metrics
The analysis of physical mechanisms is focused on DXT and DNT, the two extreme metrics related to temperature. It seems reasonable to expect that a close relation should exist with the main climatic processes that explain energy budgets and near-surface atmospheric circulation. Thus, BW (Bowen ratio = SHF/LHF), surface net solar radiation (NR) and zonal (u10) and meridional (v10) wind components and wind intensity (W 10= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u10 2 +v10 2 p ) at 10 m are employed to perform this analysis. They are related to several atmospheric climatic processes, such as cloudiness and precipitation, land-atmosphere coupling, and heat transport. These magnitudes are thought as a first-step procedure to relate these climatic mechanisms with the obtained metrics, and so we are aware that some limitations and overlapping could exist when using them. This is the case, for example, of the 10-m wind to properly describe the impact of atmospheric circulation on the variability of the annual temperature cycle, but we consider that the approach can be good enough for this first approximation. Partial correlations coefficients (ρ), using nonparametric Spearman's rank test R package (Kim, 2015) were computed, with the aim to quantify the contribution of each variable (NR, BW and W10) to each metric. It is important to point out that correlations do not mean necessarily causality, but it can be considered as indirect estimations of some relation between both magnitudes, if physical coherence exists (Seneviratne et al., 2010), taking into account that the obtained correlation signs can be relevant (Karl et al., 1993). This proposed statistical procedure has been widely used (Karl et al., 1993;Seneviratne et al., 2010;Peng et al., 2013;Du et al., 2017).
For each year, NR, BW and W10 were seasonally averaged following the classical meteorological season definition. Thus, meteorological summer comprises June-July-August (JJA) over NH and December-January-February (DJF) over SH, while winter encompasses DJF over NH and JJA over SH. Variables averaged in meteorological summer were correlated with DXT and winter with DNT. Significant correlations were considered for a 90% confidence interval.
As a measure of the influence of land-atmosphere coupling, the simple linear correlation (r) between T2M and LHF was computed (Seneviratne et al., 2006;Menéndez et al., 2016). Note that LHF and evapotranspiration ET differ from each other by an approximately constant factor, ET = LHF/λ, being λ the latent heat of vaporization (λ = 2.45 × 10 6 JÁkg −1 ) (Peixoto and Oort, 1992). Regions where evapotranspiration is limited by soil moisture are associated to r(T2M, LHF) < 0, meanwhile regions where evapotranspiration is limited by radiation show r(T2M, LHF) > 0 (Seneviratne et al., 2010). This description of land-atmosphere coupling allows the classification of climate regimes related to the limiting factors of evapotranspiration. Each grid could be classified as dry, transitional or humid, based on the theoretical framework of Seneviratne et al. (2010), and the climatic regimes from Arora (2002), which are based on the radiative aridity index (D) defined in Budyko and Miller (1974). This index is defined as D = R/(P × λ), where R is the net annual averaged radiation, P is the annual averaged precipitation and λ the vaporization latent heat. Thus, humid regions are those where there is no coupling between atmosphere and land (r(T2M, LHF) > 0); meanwhile, dry areas are those where LHF is limited by soil moisture (r(T2M, LHF) < 0) and also D ≥ 5. (This last threshold is established by Arora (2002) to define dry regions.) Transitional regions are those where landatmosphere coupling exists, but D < 5.
Moreover, in order to further analyse the effects of interannual wind variability over seasonal metrics, the interannual linear correlations between zonal (u10) and meridional (v10) wind components during summer (NH: JJA, SH: DJF) and winter with DXT and DNT, respectively, are also computed.
Both interannual partial and linear correlations were computed at each grid cell and for each reanalysis, and then averaged to obtain the ensemble value. The uncertainty of the ensemble in the interannual partial correlations was addressed as the significance agreement among reanalyses. Thus, areas where at least two out of the three reanalyses present significant correlation (α = .1) were considered as statistically significant.  Figure 2a shows the delay (in days) with respect to summer hemisphere astronomical solstice (DXT), as defined on the methodology section. In relation to an overall delay of about 30 days, some specific features are observed. Coastal regions show larger delays values than continental parts, that is, coastal areas take longer to arrive to maximum T2M that the inland parts of the continents. This delay is larger over the NH (coastal areas above 45 days, continental ones around 30 days) than over the SH (coastal areas with 30 days, continental ones around 20 days).
Day of minimum temperature (DNT, Figure 2b) looks quite homogeneous over the whole domain, with delays with respect to the winter hemisphere astronomical solstice of around 20-30 days. In that sense, it is more uniform and closer to those values than the DXT case. Nevertheless, some clear similarities also appear with DXT, like the larger delays for coastal areas when compared with continental counterparts. This difference between coastal and continental zones is larger for higher latitudes (>40 ), reaching up to 15-20 days.
The delay in the day of maximum warming (DXWARM, Figure 3a) is given with respect to the astronomical equinox day of the spring hemisphere, in consistency with the other metrics. Although large areas of delays around 30 days are obtained, the patterns and values are less homogeneous than for the first two metrics (DXT and DNT). There is a strong latitudinal pattern in both hemispheres, with smaller delays over areas close to the Tropics, and increasing as latitudes are higher, with values up to 50 days over latitudes above 50 N. This is the case also over the Mediterranean basin.
Finally, the delay of the maximum cooling day (DXCOOL, Figure 3b), referred to the day of the autumn hemisphere astronomical equinox, shows a roughly reversed latitudinal pattern compared with DXWARM ( Figure 3a). Largest delays are obtained near the subtropical areas (≈50 days), that is, the day of maximum cooling generally occurs first in the higher latitudes and later in the subtropical regions. Delays around 30 days are seen less frequently than for DXWARM, but on that reversed latitudinal behaviour, also tend to occur over mid-latitudes, and in particular, over continental regions of Eurasia or South America, with smaller delays (around 25 days) than their surrounding areas (with 30-40 days).

| Processes related to DXT and DNT metrics
Once metrics results are analysed, and their capability to describe in a realistic and consistent way the main dates of the annual cycle that can define subperiods (seasons) along the year, the climatic variables chosen (NR, BW F I G U R E 2 Extratropical spatial patterns of metrics proposed in Figure 1, averaged for the 1981-2013 period and for the ensemble of three reanalyses (ENS). Values represent the difference in days between the metric date occurrence and the corresponding astronomical season onset at each hemisphere: (a) day of maximum T2M (DXT) with respect to summer onset, (b) day of minimum T2M with respect to winter onset [Colour figure can be viewed at wileyonlinelibrary.com] and W10) that can help to explain such results are now analysed. The focus is set on the days of the maximum and minimum T2M (DXT and DNT), which probably show a strong relationship with atmospheric circulation and with the main climatic variables that determine the surface energy balance. Partial linear correlations, as described in the methodology section are now displayed, computed as ρ(A, B), here A is the chosen variable removing the effect of the other variables (NR, BW, W10) and B is each of the metrics (DXT, DNT). These interannual correlations were computed at each grid cell for each reanalysis, and then averaged to obtain the ensemble value. For these statistics, the absolute DXT and DNT values are used, not the delays with respect to astronomical dates. Figure 4 shows the interannual partial correlations between DXT and the variables selected as representative of processes during the corresponding meteorological summer (JJA over NH, DJF over SH) for the period 1981-2013. Note that areas where the partial correlations were considered as statistically significant are outlined. Larger partial correlations are obtained for NR (Figure 4a), compared with the other two magnitudes (BW and W10). In areas where these partial correlations are negative, increases of NR are related to decreases of DXT, that is, larger net surface radiation correspond to maximum annual temperature earlier days. This effect is usually seen in humid areas, where no land-atmosphere coupling exists (Figure 4d, r (T2M, LHF) > 0) and so net radiation controls surface air temperature. This is especially for latitudes over 40 N where significant partial correlation areas are found in North America and Eurasia. Positive correlations between DXT and NR (later/advanced occurrence of DXT related to increased/decreased NR) generally coincide with dry or transitional regions, where summer land-atmosphere coupling exists (Figure 4d, r(T2M, LHF) < 0). This is the case of the southern part of North America, inland part of F I G U R E 3 As Figure 2, but for (a) day of the maximum warming of T2M (DXWARM) relative to spring onset and (b) day of the maximum cooling (DXCOOL) relative to autumn onset [Colour figure can be viewed at wileyonlinelibrary.com] Eurasia, Sahara desert and western Australia that show areas with significant partial correlations, but also, Arabian Peninsula, southern South America and western South Africa. Regarding the relationship between DXT and BW (Figure 4b), significant partial correlations are mainly negative, especially in the transition regimes ( Figure 4d). This suggests that during meteorological summer, when BW is increased (related to large SHF and relatively dry soils), DXT value is decreased, that is, the day where maximum T2M occurs earlier in the year. Figure 4c shows the interannual correlation between DXT and wind intensity at 10 m (regardless the wind direction). The areas with positive (negative) correlations are characterized by delays (advances) in the day of occurrence of DXT associated with more windy summers. However, there are few regions where these relationships are significant (e.g., some areas of northern North America and Eurasia).

| DXT-related processes
To establish more clearly the impact of interannual wind variability, the interannual linear correlations of its zonal (u10) and meridional (v10) components with DXT are shown in Figure 5. In this case, the sign of the wind components is considered, that is, u10 > 0 means westerly wind and v10 > 0 means southerly wind in both hemispheres. Variations in the zonal component are associated with oceanic (in areas near the windward coasts) or continental (in the interior of the large continents of the NH) influence. Variations in the meridional component are linked to seasons characterized by cold (stronger northerly winds in NH or southerly winds in SH) or warm (stronger southerly winds in NH or northerly winds in SH) flows. The timing of the maximum summer temperature in some regions is related to the interannual variations of the zonal wind and, in others, of the meridional wind. For example, negative significant correlation between u10 and DXT ( Figure 5a) over central Canada indicates that an intensification (weakening) of the westerlies is associated with an advance (delay) of DXT. In contrast, in the central United States the relationship is reversed: positive significant correlations (intensification of westerlies are related to DXT delay). The influence of variations in meridional flow is evident (Figure 5b), for example, in southeastern South America, where significant r(v10, DXT) > 0 suggests that in those summers with stronger warm northerlies DXT is advanced. Summer timing in western Russia is also sensitive to the variability of the meridional flow: the significant negative correlation in that area indicates that in summers with stronger southerly wind DXT is advanced while the opposite occurs if the cold northerly flow increases.  Figure 6 shows interannual partial correlations between the DNT metric and the variables NR, BW and W10 during the meteorological winter (DJF over NH, JJA over SH). As in the previous figure, note that areas where the partial correlations were considered as statistically significant are outlined. When NR is analysed (Figure 6a) significant positive partial correlations are observed over central and northeastern of North America and over large parts of Eurasia. These regions tend to coincide with humid regions where land and atmosphere are not coupled in winter (r(T2M, LHF) > 0, Figure 6d). Meanwhile, areas over southwestern North America, northeast of Asia, south of Africa and Australia show significant negative correlations between DNT and NR, that is, increases in NR are related to smaller DNT values (earlier days). These zones coincide with transitional and dry regions, where land-atmosphere coupling exists (r(T2M, LHF) < 0, Figure 6d). Regarding BW (Figure 6b) there is no clear spatial pattern, with both positive and negative values. However, there are some areas with more homogeneous negative correlations (e.g., southeastern United States and southeastern China) suggesting some modulation in DNT timing due to the partitioning between SHF and LHF. The partial correlation between DNT and W10 (Figure 6c) shows some well-defined regional patterns with significant r(W10, DNT) > 0 over western Europe, northern Eurasia, western North America and the Patagonian Andes. In these regions, more (less) windy winters correspond to delays (advances) in the coldest day date. Since these regions are exposed to the flow from the oceans, this result may suggest a possible ocean-continent coupling mechanism that contributes to modulate the timing of winter temperature.

| DNT-related processes
The interannual linear correlations between DNT and zonal (u10) and meridional (v10)  westerlies delay/advance DNT), and negative ones for example, in Alaska and eastern North America. Whereas the meridional component variations show significant correlations with DNT (Figure 7b), especially in the NH. In general, r(W10, DNT) > 0 over large parts of North America, Eurasia and around the Mediterranean Sea, indicating that southerly (northerly) winds delay (advance) the occurrence of the minimum temperature.

| DISCUSSION
This work has analysed the main characteristics of the annual cycle timing as obtained from T2M proposed metrics, computed from the days of maximum and minimum temperature, and those related to the dates of maximum and minimum rhythm change (fastest warm-up and cooldown): DXT, DXN, DXWARM, DXCOOL. Regions of analysis span from the Tropics to 60 latitude, on both hemispheres. This proposal was applied to three independent reanalysis (ERA-Interim, JRA-55 and NCEP2), for 1980-2014 period. Three processes have been studied with the aim to understand the interannual variations of such metrics.
Extreme metrics (DXT and DNT) showed similar patterns with around 20 to 30 days on average lagged the respective astronomical solstice. This overall result, although somewhat widely known, has been just briefly documented on scientific literature (Trenberth, 1983;Douglass et al., 2004), and not based on the more used reanalysis databases. Delays of around 30 days between solar incoming radiation and temperature extremes over extratropical zones are found. Therefore, this can be considered a remarkable result from the proposed methodology. Smaller delays (day of extreme temperature closer to astronomical value) over inland is obtained over wide regions, when compared with coastal values, which is also in agreement with previous studies (Prescott and Collins, 1951;El-Hussainy and Essa, 1997;Donohoe et al., 2020). This spatial contrast is clearer for the NH than for the SH, which could be related to the oceanic influence over the metrics (related to the water to land heat capacities; Ahrens and Henson, 2016;Donohoe and Battisti, 2013;Trenberth, 1983). Thus, continental regions, where heat capacity is lower than in the ocean, energy income/outcome leads to earlier day of maximum/ minimum temperature (DXT/DNT). Over the SH, the oceanic influence smooths the temperatures, provoking a more homogeneous spatial pattern for the day of temperature extremes, and so the computed metrics, reducing the inland/coast contrast. Meanwhile, DXWARM and DXCOOL metrics (Figure 3) show reversed meridional gradients between each other, from the Tropics to mid-latitudes. These results are related to the meridional contrast of insolation at the spring and autumn equinoxes. As a measure of the uncertainty of the reanalysis ensemble, Figure S1, Supporting Information shows the intermodel range computed as max(ENS)-min(ENS) for each seasonal metric. In general, the spread between reanalyses remains reasonably narrow over most of the continental areas, especially for DXT. Eighty percent of the pixels have an inter-reanalysis range below 4.3 days (DXT), 9.8 days (DNT), 8.3 days (DXWARM) and 13.3 days (DXCOOL). The largest uncertainties are usually found in areas with complex orography and in high latitudes.
It is reasonable to expect that the variables studied, relevant in the land-atmosphere interface, related to energy budget (NR, BW) and atmospheric circulation (nearsurface wind) may contribute to understanding the timing of the metrics. However, at the same time, some shortcomings are inherent to this analysis. First, the limitations of the methodology, for example, when it does not take into account the intra-seasonal variability of the chosen variables, which could be important (Cornes et al., 2017). Second, the limitations of the reanalysis products (the global models used, the assimilation methodologies with the available observations). Third, the intrinsic internal variability among the complex feedbacks between land, atmospheric boundary layer and large-scale circulation. Yettella and England (2018) indicates that internal variability makes complex the capability to detect changes on the seasonal cycle over many regions. These difficulties are manifested through low or not significant correlations between seasonal metrics and climate variables. However, we consider that our results are generally consistent and reflect different processes that affect the timing of the annual temperature cycle. These processes (related to the variability of net surface radiation, land-temperature coupling or atmospheric flux) vary according to climate regime and season. A more detailed analysis focused on climate regimes is now presented.

| Humid regions
In humid regions, where there is no coupling between atmosphere and land, the summer partial correlations between DXT and NR have greater spatial homogeneity and significance than the correlations of DXT with either BW or W10 (Table 1). This suggests that variations in net surface radiation are the main driver of variations in DXT timing. When looking at DNT, the relationship with W10 shows the highest percentage of cells with statistically significant partial correlations respect to RN and BW. This means that the interannual variability of circulation play a major role. The significant correlation areas between NR and DXT ( Figure 4a) are mostly negative over humid regions, that is, increased NR is related to earlier day of maximum temperature (DXT). As no landatmosphere coupling exists, air temperature is just controlled by the atmosphere (net radiation). For DNT, most of significant correlations with W10 ( Figure 6c) are positive. Variations in wind intensity and direction are related to the strength of the oceanic influence on the continents and to variations in cold or warm flows, mainly affecting the timing of the minimum temperature in winter. The pattern of opposite correlation sign obtained over north of Eurasia and north of Africa, respectively, is reminiscent of the NAO mode of variability. This fact suggests the hypothesis (to be explored in future research) that the NAO could modulate the timing of the annual temperature cycle over the continents. The hypothesis is consistent with previous studies on the impact of the modes of variability on NH winter climate (Thompson and Wallace, 1998;Hurrell et al., 2001). The positive phase of the NAO is related to anomalously intense westerly winds, transporting heat and moisture from the ocean to northern Europe. This results in increased temperature and precipitation (Hurrell, 1995;Rodwell et al., 1999), which could lead to later DNT values. In addition, this mode of variability is also related to the decrease in temperature over North Africa, which could be consistent with earlier DNT values. Similarly, a possible influence of the SAM mode of variability (Thompson and Wallace, 2000;Thompson and Solomon, 2002) on the timing of temperature for SH could be hypothesized. The high correlations between DNT and wind observed over southern South America, New Zealand and Tasmania could be associated with the SAM driving westerly winds. In particular, the SAM positive phase is associated with warm temperature anomalies over Patagonia (Garreaud et al., 2009), and so a later day of minimum temperature occurrence (DNT).

| Transitional regions
DNT timing is mainly related to NR, while for DXT the variability of W10 is the most influential contributor (Table 1). At the same time, the land-atmosphere interaction (represented by the BW magnitude) is also relevant for DXT and DNT values. These results would point to a relevant role of the partitioning of sensible and latent heat fluxes for the timing of the annual temperature cycle in transition regions. These regions are characterized by variability in soil moisture conditions, alternating dry and humid years. Thus, for example, a humid summer year leads to larger fraction of net radiation used to evaporate soil moisture compared to a dry summer, causing a delay on DXT date (consistent with the negative correlations between BW and DXT in Figure 4b) due to (a) a decrease on the percentage of available energy to increase air temperature (lower sensible heating) and (b) the evaporative cooling process. These results are coherent with those of Menéndez et al. (2019Menéndez et al. ( , 2016, Seneviratne et al. (2010Seneviratne et al. ( , 2006 and Spennemann et al. (2018), in which the importance of soil moisture-atmosphere coupling processes for near-surface air temperature during the summer season over transition regions is found.

| Dry regions
DXT and DNT are related mainly to net radiation, in terms of the largest areas with significant partial correlation between NR and any of them (Table 1), indicating that surface radiation variability is relevant to these seasonal metrics. Over these dry regions, DXT is positively correlated to NR, suggesting that summers with positive net radiation anomalies are related to delays on DXT days. Under dry soil conditions, evapotranspiration and temperature variability is affected by precipitation variability . Net radiation is positively correlated with precipitation over these dry areas ( Figure S2). Therefore, a possible chain of interactions is as follows: more NR, more precipitation, more soil moisture, more evapotranspiration and less sensible heating, less temperature, which in turn induces a later DXT. DNT timing in dry regions is influenced by interannual T A B L E 1 Percentage of cells with statistically significant partial correlations (α = .1) between metrics (DXT, DNT, DXWARM, DXCOOL) and selected variables (NR, BW, W10), for at least two of the three reanalysis (ERA-Interim, JRA-55 and NCEP2) variations of NR, BW and W10, with similar relative importance for the three variables.

| CONCLUSIONS
A proposal to characterize the annual daily temperature cycle over extratropical land regions based on the computation of four critical dates (metrics) from the annual cycle characteristics has been presented. The results indicate that it can be considered a valuable procedure to quantify the subperiods (seasons) of the annual cycle and to describe their interannual variability compared to different variables representative of physical processes. The four dates show values with a coherent geographical structure, related to the astronomical dates references. The main mechanisms that could explain the interannual variations of these metrics are those affecting surface energy fluxes, including land-atmosphere coupling processes and circulation variability. Near-surface atmospheric circulation (which introduces oceanic influence, as well as cold or warm flows) and net surface radiation are clear influences on the DXT and DNT metrics. It is useful to consider the different climatic regimes related to land-atmosphere coupling. Thus, in humid regions, where there is no land-atmosphere coupling, summer radiation variability is essential for the timing of DXT, while winter circulation variability is most important for DNT. In transition regions, where landatmosphere coupling plays an important role in the energy balance, the partitioning of energy between latent and sensible heat fluxes is a dominant factor for the timing of metrics, although in summer circulation variability is also relevant, causing DXT dates to be advanced or delayed. Lastly, the annual temperature cycle in dry regions is mainly determined by radiative processes. The fact that reanalyses provide comprehensive and physically consistent information on atmospheric conditions at regular intervals and over periods of several decades made possible the present study on the timing of the annual temperature cycle and possible coadjuvant factors, on global and interannual scales. However, while reanalyses may have shortcomings, the major drawback is that their errors and uncertainties are not quantified and are often insufficiently understood (Parker, 2016). Our work, by comparing reanalysis datasets with each other, quantifying the spread between the metrics calculated with the different reanalyses, and analysing the significance of the results from the ensemble, is a contribution towards a better description of this uncertainty. Our approach can be used with other gridded temperature databases from both observational products and climate models. It can also be adapted to analyse the annual cycle of other climatic variables and vegetation indices for the study of phenological seasons. The proposed method to compute seasons presented here is somewhat complementary to a proposal already presented by the authors, based on thresholds associated to 25th and 75th percentiles of maximum and minimum daily temperature (L opez-Franca et al., 2013), with the aim to inspect different approaches to characterize the annual cycle of temperature.