Quantifying changes in spatial patterns of surface air temperature dynamics over several decades

We study daily surface air temperature (SAT) reanalysis in a grid over the Earth's surface to identify and quantify changes in SAT dynamics during the period 1979–2016. By analysing the Hilbert amplitude and frequency we identify the regions where relative variations are most pronounced (larger than ±50 % for the amplitude and ±100 % for the frequency). Amplitude variations are interpreted as due to changes in precipitation or ice melting, while frequency variations are interpreted as due to a northward shift of the inter-tropical convergence zone (ITCZ) and to a widening of the rainfall band in the western Pacific Ocean. The ITCZ is the ascending branch of the Hadley cell, and thus by affecting the tropical atmospheric circulation, ITCZ migration has far-reaching climatic consequences. As the methodology proposed here can be applied to many other geophysical time series, our work will stimulate new research that will advance the understanding of climate change impacts.

1 Introduction

The unprecedented intensification of weather extremes is motivating research aimed at understanding long-term climatic variations (Barreiro et al., 2008; Cai et al., 2014; Coumou and Rahmstorf, 2012; England et al., 2014; Turco et al., 2015) that can have profound socio-economic impacts (Ghil et al., 2011) and trigger complex ecological adaptation mechanisms (Beaumont et al., 2011; Bordeu et al., 2016; Gottfried et al., 2012; Lejeune et al., 2002).

Figure 1 Relative change in the time-averaged Hilbert amplitude. (a) Amplitude averaged over the first 10 years (January 1979 to December 1988). (b) Amplitude averaged over the last 10 years (July 2007 to June 2016). (c) Relative change in the Hilbert amplitude, $\mathrm{\Delta }a/<a>$. (d) Relative change in amplitude of the seasonal cycle computed from the amplitude of the climatology, Δa^(clim)∕a^(clim). A good qualitative agreement is seen in the spatial structures in (c) and (d). Importantly, the structures uncovered by the Hilbert amplitude are well defined in comparison with those uncovered by the analysis of the climatology amplitude, which look noisier.

Quantifying variations in surface air temperature (SAT) dynamics over several decades is a challenging problem because of non-stationarity and the presence of trends, measurement noise, multiple timescales, memory, and correlations in the data (Franzke, 2012; Massah and Kantz, 2016); in addition, reanalysis data can be unreliable (due to the lack of observational constraints in many geographical regions), and reanalysis time series are insufficiently long (as reanalysis starts at the beginning of the satellite era). These challenges have motivated the use, for climate data analysis, of data-driven approaches that have been commonly used for investigating observed complex signals in other fields of science (e.g. neurological, physiological, financial, etc.). Univariate analysis tools that have been used to analyse SAT time series include detrended fluctuation analysis, fractional analysis, and wavelet analysis. Bivariate analysis and the complex network approach has also allowed us to uncover inter-relations between SAT anomalies in different regions (Barreiro et al., 2011; Donges et al., 2009; Tsonis and Swanson, 2008). In this approach the seasonal cycle is removed to eliminate the influence of solar forcing, and the links represent correlations (linear or non-linear) or statistical similarities between SAT dynamics in different areas (Tirabassi and Masoller, 2016). On the other hand, changes in the SAT seasonal cycle have also been investigated, and a trend toward reduced cycle amplitude has been detected in many regions (Chambers et al., 2013; Duan et al., 2017; Dwyer et al., 2012; Qian et al., 2011; Stine and Huybers, 2012; Stine et al., 2009; Wang and Dillon, 2014). However, changes in SAT dynamics over several decades (such as those observed in Fig. 1) have not yet been investigated at a global scale. In order to fill this gap, we use Hilbert analysis (described in the Supplement) to investigate SAT time series with daily resolution (reanalysis covering the Earth's surface in the period 1979–2016). Our goal is to detect the most sensitive regions (“hotspots”) where variations in SAT dynamics over the last decades are more pronounced.

The Hilbert transform (HT) provides, for a real oscillatory time series, x(t) with t ∈ [1, T], an instantaneous amplitude, a(t), and an instantaneous frequency, ω(t), for each data point of the time series and thus allows us to characterise how the amplitude and the frequency of a signal vary in time. If a signal does not have a sufficiently narrow frequency band, a(t) and ω(t) will not have a clear physical meaning (Pikovsky et al., 2001). The usual solution is based on band-pass filtering to isolate a narrow frequency band; however, HT directly applied to the signal can still yield useful information. An alternative solution is based on the Hilbert–Huang transform (Huang et al., 1998) that combines Hilbert analysis with empirical mode decomposition that decomposes an arbitrary real time series into components, each having the physical meaning of a rotation in the complex plane.

Because many natural geophysical time series have a seasonal periodicity, this has motivated the use of Hilbert analysis to characterise the time-varying oscillation amplitude and to investigate phase shifts and phase–amplitude couplings. Applications in various geophysical fields are discussed in Huang and Wu (2008). As more recent examples, Massei and Fournier (2012) used Hilbert analysis to characterise the daily variability of the Seine river flow from 1950 to 2008, uncovering linkages between river flow variability and global climate oscillations (the North Atlantic Oscillation and the Madden–Julian Oscillation). Sun (2015) used Hilbert analysis to compute the daily phase shift between temperature signals recorded at the ground surface and at a depth of 5 m in two meteorology stations in Taiwan from 1952 to 2008. Significant reductions in the phase shift from the 1980s to 1990s were found, which was interpreted to be related to the warming of the Pacific Decadal Oscillation. Reddy and Adarsh (2016) applied Hilbert analysis to rainfall time series in India and found that the multi-scale components of rainfall series have a similar periodic structure as global climate oscillations (the Quasi-biennial Oscillation, El Niño Southern Oscillation, etc.).

We have recently applied the Hilbert transform to unfiltered daily SAT reanalysis (Zappalà et al., 2016). We have shown that the maps of time-averaged Hilbert frequency, $<\mathit{\omega }>$, and standard deviation, σ_ω, revealed well-defined large-scale structures which were consistent with known dynamical processes.

Here we use a(t) and ω(t) to quantify SAT variations. Our hypothesis is that changes in a(t) and ω(t) can yield information about variations in SAT dynamics. Specifically, we are interested in addressing the following questions: which properties of a(t) and ω(t) display relevant variations? Where are the regions in which these variations are more pronounced? Which processes can be responsible for these variations? Can these variations be used as a quantitative measure of regional climate change?

2 Data

In the main text we present results from an ERA-Interim daily SAT reanalysis (Dee et al., 2011) that covers the period from January 1979 to June 2016 with a spatial resolution of 2.5^∘, both in latitude and in longitude. Thus, there are N = 73 × 144 = 10 512 geographical sites and in each site the SAT time series has T=13696 days. In the Supplement we compare ERA-Interim with NCEP-DOE Reanalysis 2, which is an improved version of the NCEP Reanalysis I model (Kistler et al., 2001). It covers a longer time interval and has 94 × 192 = 18 048 geographical sites. In order to perform a precise comparison between the results of the two datasets, in the NCEP-DOE Reanalysis 2 we consider the same time interval as the ERA-Interim dataset.

3 Methods

3.1 Hilbert analysis

To apply the Hilbert transform (described in the Supplement) we first pre-process each raw SAT time series, r_j(t) (where j ∈ [1, N] represents the geographical site and t ∈ [1, T] represents the day): we eliminate the linear trend and normalise to zero mean and unit variance, obtaining x_j(t). The Hilbert transform is then applied to x_j(t), obtaining y_j(t) = HT[x_j(t)]. From x_j(t) and y_j(t), the amplitude a_j(t) and the phase φ_j(t) were calculated as a_j(t) = $\sqrt{\left[x_j\right(t){]}^{\mathrm{2}}+[y_j\left(t\right){]}^{\mathrm{2}}}$ and φ_j(t) = arctan[y_j(t)∕x_j(t)]. Taking into account the signs of x_j(t) and y_j(t), the phase is constrained to the interval [−π, π). Whenever an extreme value of the interval is reached, the phase jumps to the other end of the interval. By eliminating these sudden jumps (using a standard library function that appropriately adds ±2π) we “unwrapped” the phase obtaining a continuous variation in time from which the frequency time series, ω_j(t), was obtained by calculating the derivative. Since the Hilbert algorithm (Bilato et al., 2014) gives deviations from the true values of the amplitude, phase, and frequency near the extremes, in each time series (a_j(t), φ_j(t), and ω_j(t)) we disregarded the initial and final 5% (see the Supplement). This way, we have time series of length T = 12 328.

3.2 Measures used to quantify variations in SAT dynamics

Variations in the Hilbert amplitude were quantified by the relative change, $\mathrm{\Delta }a/<a>$ = ($<a{>}_{\mathrm{l}}$ − $<a{>}_{\mathrm{f}})/<a>$, where $<a{>}_{\mathrm{f}}$ is the average value of the amplitude during the first 10 years of the time series (January 1979 to December 1988), and $<a{>}_{\mathrm{l}}$ during the last 10 years (July 2007 to June 2016). Analogously, we calculated the relative change in amplitude variance, $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}$, of average frequency, $\mathrm{\Delta }\mathit{\omega }/<\mathit{\omega }>$, and of frequency variance, $\mathrm{\Delta }{\mathit{\sigma }}_{\mathit{\omega }}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathit{\omega }}^{\mathrm{2}}$. In the Supplement we analyse how the spatial structures uncovered depend on the time intervals used to calculate the relative variations: we compare with relative variations during the first and final 5 years of the reanalysis and also during the first half-period and the second half-period of the reanalysis. While the values of the relative variations vary with the time interval considered, the spatial maps are remarkably robust as the same structures are found with the three time intervals considered.

A similar analysis was performed to detect changes directly from the raw SAT time series, r_j(t), by computing the amplitude of the climatology (or seasonal cycle), c_j(t), and the variance of anomaly time series, z_j(t).

Specifically, the amplitude of the climatology was calculated as $a_j^{\left(\mathrm{clim}\right)}\left(I\right)$ = max[$c_j^I\left(t\right)$] − min[$c_j^I\left(t\right)$], where $c_j^I\left(t\right)$ is the climatology series calculated only in the time interval I. We remark that the climatology amplitude $a_j^{\left(\mathrm{clim}\right)}\left(I\right)$ is a scalar number that depends on the choice of the time interval I. We calculated the climatology amplitude in the first and last decade, as well as in the whole series. As before, we used these values to calculate the relative change Δa^(clim)∕a^(clim). Also, the variance of the anomaly time series z_j(t) was calculated and then used to find the relative change, $\mathrm{\Delta }{\mathit{\sigma }}_z^{\mathrm{2}}/{\mathit{\sigma }}_z^{\mathrm{2}}$.

With the goal of relating changes in Hilbert frequency with changes in the statistical properties of SAT time series, an analysis of the number of zero crossings was performed: for each x_j(t) we counted the number of crossings through the mean value, x = 0. As with other quantities, we then calculated the relative change.

Figure 2 Surface air temperature in two regions where a clear change in the oscillation amplitude in the last 10 years is observed with respect to the first 10 years. (a) Site of coordinates (7.5^∘ S, 307.5^∘ E) marked with a triangle in Fig. 1c. (b) Site of coordinates (75^∘ N, 40^∘ E) marked with a circle in Fig. 1c.

3.3 Significance analysis

A statistical significance analysis was performed by surrogating Hilbert series. For each amplitude time series (i.e. in each grid point) 100 shuffle surrogates were generated and for each surrogate the relative change, $\mathrm{\Delta }{a}^{\mathrm{s}}/<a^{\mathrm{s}}>$, was calculated. Then, the average over the 100 surrogates, $<\mathrm{\Delta }{a}^{\mathrm{s}}/<a^{\mathrm{s}}{\mathit{>>}}_{\mathrm{s}}$, and its standard deviation, σ_s, were used to define the significance threshold: the relative change computed from the original data was considered significant if it was higher than $<\mathrm{\Delta }{a}^{\mathrm{s}}/<a^{\mathrm{s}}{\mathit{>>}}_{\mathrm{s}}$ + 2σ_s or lower than $<\mathrm{\Delta }{a}^{\mathrm{s}}/<a^{\mathrm{s}}{\mathit{>>}}_{\mathrm{s}}$ − 2σ_s. In the colour maps, regions where variations are not significant are displayed in white. The same test was applied to frequency variations and the other quantities, except for the climatology for which a surrogate test is not applicable. In the Supplement various thresholds are considered and, in addition, a non-parametric significance test is used. Here we present only the maps obtained with threshold ±2σ because it is a compromise between uncovering the spatial regions where SAT changes are pronounced and disregarding the areas where the variations are small.

4 Results

We analyse the maps of $<\mathit{\omega }>$, $<a>$, ${\mathit{\sigma }}_{\mathit{\omega }}^{\mathrm{2}}$, and ${\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}$ in the first 10 years and in the last 10 years of the period covered by the reanalysis, as well as the relative change between the two decades.

4.1 Analysis of amplitude variations

Figure 1a and b display $<a>$ in the first and in the last 10 years, respectively, and Fig. 1c displays the relative difference (see Sect. 3 for details). In Fig. 1c we see an area of large increase (more than 50 %) in average amplitude located in South America (red spot marked by a triangle) and an area of large decrease (again, more than 50 %) located in the Arctic (blue spot marked by a circle). The raw SAT time series in these regions are displayed in Fig. 2.

In both time series we clearly observe a change in the amplitude of the oscillations in the last 10 years with respect to the first 10 years, having a visual confirmation of the changes detected by the Hilbert amplitude. The red spot in Amazonia, whose SAT series shown in Fig. 2a has an increasing amplitude, can be interpreted in terms of changes in precipitation. In particular, the increase in the Hilbert amplitude is linked to the decrease in precipitation and to the lengthening of the dry season (as reported in Fu et al., 2013; Gu et al., 2016; Liebmann et al., 2004). This is due to the fact that when precipitation decreases, the fraction of solar radiation that is not used for evaporation is used to heat the ground, which in turn heats the surface air. This leads to higher extreme temperatures during the dry seasons, as can be observed in Fig. 2a. Regarding the blue spot in the Arctic region where the SAT series shown in Fig. 2b has a decreasing amplitude, it can be interpreted as due to the melting of sea ice. In fact, when ice is present at the surface of the sea, it acts as an insulator preventing heat exchange between sea and air. This causes a large amplitude in the SAT cycle. On the other hand, if the ice melts, the air–sea heat exchange reduces the amplitude of the cycle. In particular, during winter the air temperature is mitigated by the sea and tends to have more moderated values. It is important to take into account that this blue spot is in a region for which the observational constraints from satellites on the reanalysis are scarce, which decreases the quality of the reanalysis in the region. Therefore, in order to check whether the detected changes are robust, we performed the same analysis using the NCEP-DOE reanalysis dataset. The results, presented in the Supplement, confirm the presence of the blue spot in the Arctic.

Next, we compare the changes detected by the Hilbert amplitude with those computed directly from SAT (by decomposing the SAT time series into climatology and anomaly, as explained in Sect. 3). Since the climatology term retains the seasonal variation, we expect its amplitude change to give similar indications as the Hilbert amplitude change. On the other hand, the anomaly term contains all the rapid variability, so we expect its variance to give similar results as the variance of the Hilbert amplitude.

Figure 3 Relative change in amplitude fluctuations computed from the variance of the (a) Hilbert amplitude, $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}$; (b) the anomaly time series, $\mathrm{\Delta }{\mathit{\sigma }}_z^{\mathrm{2}}/{\mathit{\sigma }}_z^{\mathrm{2}}$.

Figure 4 Relative change in the time-averaged Hilbert frequency (in units of oscillations per year). (a) Average in the first 10 years (1979–1988). (b) Average in the last 10 years (2007–2016). (c) Relative change in Hilbert frequency, $\mathrm{\Delta }\mathit{\omega }/<\mathit{\omega }>$. (d) Relative change in the number of zero crossings of the normalised SAT time series. In (a) and (b) the colour scale is adjusted to represent in white the regions where the average frequency is one oscillation per year. In (c) and (d), a good qualitative agreement of spatial structures is seen; however, we note that the Hilbert frequency detects stronger variations than those measured by the number of zero crossings.

Figure 1c and d, which respectively display the relative change in Hilbert amplitude and in climatology amplitude, and Fig. 3a and b, which respectively display the relative change in Hilbert amplitude variance and in anomaly variance, confirm these expectations.

The good qualitative agreement seen in the spatial structures in these maps confirms that Hilbert analysis directly applied to unfiltered SAT indeed gives a physically meaningful instantaneous amplitude, with average and variance values that are consistent with those computed from SAT.

In Fig. 3a and b, however, there is a difference in the eastern Pacific Ocean in the area marked with a circle. In particular, in Fig. 3b there is an area with large decrease in variance (dark blue, around −100 %), while in Fig. 3a the decrease is less pronounced (light blue, around −65 %) and extended over a smaller area. In addition, in Fig. 3a there is a reddish orange area that indicates a moderate increase in variance (around 45 %), while in Fig. 3b such an area is absent. The reasons underlying these differences will be discussed later.

Figure 5 Normalised SAT time series and number of zero crossings in the regions indicated with a circle in Fig. 4c. In the red region (2.5^∘ N, 245^∘ E) (a), the number of zero crossings increases in the last 10 years with respect to the first 10 years (289 and 202, respectively), while in the blue region (7.5^∘ S, 250^∘ E) (b), it decreases (128 and 258 in the last and first 10 years).

4.2 Analysis of frequency variations

Figure 4a displays the average frequency $<\mathit{\omega }>$ in the first 10 years, Fig. 4b in the last 10 years, and Fig. 4c displays the relative change, $\mathrm{\Delta }\mathit{\omega }/<\mathit{\omega }>$. In Fig. 4c we note that in the eastern Pacific Ocean there are two small areas, enclosed by the circle, of intense increase (red) and decrease (blue) in frequency. They both represent frequency changes whose absolute values are larger than 100 % and correspond to the same region where differences were detected in Fig. 3.

These two areas of opposite signs suggest that, between the initial and the final decade, there is a shift of the inter-tropical convergence zone (ITCZ) toward the north. The ITCZ involves strong convective activity, which causes rapid fluctuations of SAT, thus giving high values of instantaneous frequency, as shown in Fig. 4a and b. Therefore, in the relative change in frequency, in regions corresponding to the initial position of the ITCZ we see a decrease, while in regions corresponding to the present position of the ITCZ we see an increase. For the same reason, the two red areas in the western Pacific Ocean (indicated by two squares) suggest an expansion of the tropical convective regions. This interpretation is in agreement with previous works that have related a northward shift of the ITCZ to an inter-hemispheric temperature gradient, such as the one experienced during the last decades (Frierson and Hwang, 2012; Kang et al., 2009; Schneider et al., 2014; Talento and Barreiro, 2016; Yoshimori and Broccoli, 2008). Regarding the red areas in the north Atlantic, the north Pacific, and the south Pacific, they are consistent with an increase in the occurrence of fronts, which cause large daily fluctuations of temperature and thus an increase in the Hilbert frequency.

To gain insight into the physical meaning of the changes that are captured by the Hilbert frequency, we use an alternative approach to estimate frequency variations: we define “events” as the zero crossings of SAT time series (Pikovsky et al., 2001) (detrended and normalised to zero mean as described in Sect. 3). Then, we count the number of events in the first 10 years and in the last 10 years and calculate the relative variation.

Figure 4d displays the map of relative change in zero crossings. We see that there is a qualitative good agreement with the spatial structures seen in Fig. 4c, thus providing a physical interpretation for the observed variation in the Hilbert frequency: the areas where the frequency increases (decreases) correspond to areas where the number of zero crossings increases (decreases). We note that the relative variations in the Hilbert frequency are more pronounced than those in the number of crossings, and this specifically holds in the regions where frequency variations are interpreted in terms of ITCZ migration.

Figure 5a and b display SAT time series in the dipole region indicated with the circle in Fig. 4c and also indicate (in red) the zero crossings. We can understand the difference that was detected in this region between the variance of the Hilbert amplitude (Fig. 3a) and the variance of the anomaly (Fig. 3b). This difference is explained in the following terms: in the first decade the seasonal cycle is more irregular than in the last decade, probably a consequence of an El Niño event in 1982–1983. The anomaly series contains these slow fluctuations as well as the rapid ones, and thus its variance is affected by both effects. In contrast, the Hilbert amplitude is less affected by the slow fluctuations as its variance captures mainly the rapid fluctuations of SAT.

To demonstrate the robustness of our findings, in the Supplement we compare the results obtained from ERA-Interim with those obtained from another reanalysis dataset, NCEP-DOE. We find a good qualitative agreement in the spatial structures in the maps of $<\mathit{\omega }>$, $<a>$, ${\mathit{\sigma }}_{\mathit{\omega }}^{\mathrm{2}}$, and ${\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}$, but we also discuss some relevant differences. In addition, to further understand the relationship between the statistical properties of SAT and those of the Hilbert amplitude and frequency, in the Supplement we apply Hilbert analysis to synthetic data generated by an autoregressive AR(1) process. We chose an AR(1) process because it is commonly used in the literature to model climate data. We find that, when increasing the noise intensity in the synthetic series, the Hilbert amplitude decreases while the frequency increases, which shows that this trend is also observed in real SAT time series.

5 Conclusions

We have used Hilbert analysis to quantify the changes in SAT dynamics, on a global scale, that have occurred over the last 3 decades. From the SAT time series with daily resolution we derived the amplitude and the frequency time series and then calculated the relative change (between the first and the last decade) in the average and variance of these series. Large variations in the Hilbert amplitude (more than 50 %) in the Arctic and in Amazonia were respectively interpreted as due to ice melting and precipitation decrease. The analysis of the Hilbert frequency also uncovered areas of large changes. In particular, two areas of opposite changes in the eastern Pacific Ocean and two areas of increase in the western Pacific Ocean suggest a shift towards the north and a widening of the ITCZ. While there is evidence that the ITCZ has moved north–south in the past, to the best of our knowledge our work is the first to confirm this migration in the last decades. Our findings have important implications because, as the ITCZ is the ascending branch of the Hadley cell, its migration affects both the Earth's radiative balance and the release of latent heat that drives the tropical atmospheric circulation. Taken together, these effects have not only local but also far-reaching climatic consequences. Additional analysis provided in the Supplement confirms the robustness of these observations.

As the methodology used here can be applied to many other climatological time series that exhibit well-defined oscillatory behaviour, we believe that our work will stimulate new research to identify and quantify the impacts of climate change directly from observed data.