Abstract

The first four statistical moments and their trends are calculated for the average daily surface air temperature (SAT) from 1950 to 2010 using the Global Historical Climatology Network–Daily station data for each season relative to the 1961–90 climatology over the Northern Hemisphere. Temporal variation of daily SAT probability distributions are represented as generalized linear regression coefficients on the mean, standard deviation, skewness, and kurtosis calculated for each 10-yr moving time window from 1950–59 to 2001–10. The climatology and trends of these statistical moments suggest that daily SAT probability distributions are non-Gaussian and are changing in time. The climatology of the first four statistical moments has distinct spatial patterns with large coherent structure for mean and standard deviation and relatively smaller and more regionalized patterns for skewness and kurtosis. The linear temporal trends from 1950 to 2010 of the first four moments also have coherent spatial patterns. The linear temporal trends in the characterizing statistical moments are statistically significant at most locations and have differing spatial patterns for different moments. The regionalized variations specific to higher moments may be related to the climate dynamics that contribute to extremes. The nonzero skewness and kurtosis makes this detailed documentation on the higher statistical moments useful for quantifying climate changes and assessing climate model uncertainties.

1. Introduction

Weather-related damage often occurs on the time scales of days or shorter. As such, understanding and quantifying climate risk requires consideration of daily variability as well as variability on longer time scales. For temperature, damages at (sub)daily time scales often concern human health (e.g., heat stress; Sherwood and Huber 2010) and agriculture (e.g., frost days; Meehl et al. 2004). Historically, most climate variability studies focus on either time scales of months or longer (e.g., Parker et al. 1994; Michaels et al. 1998; Solomon et al. 2007; Hansen et al. 2012), or spatial averages over regions, such as individual countries or over an entire hemisphere (e.g., Karl et al. 1995; Gong and Ho 2004; Shen et al. 2011). While these analyses deduce useful results that can aid in the understanding of global climate changes, finer space–time scale details for not only the mean but also higher statistical moments relevant to high-frequency geophysical dynamics and regional climate risk still require careful documentation as climate data and models improve. The purpose of this paper is to document the higher statistical moments, their seasonal climatology, their regional structure of temporal trends, and their associated non-Gaussian probability density functions (PDFs).

It is known that statistical moments can determine a PDF via a Taylor series expansion and Fourier transform (Wackerly et al. 2008). Perron and Sura (2013) used the third and fourth statistical moments to capture the non-Gaussian nature of daily atmospheric variables. They compiled a climatology of non-Gaussian atmospheric statistics from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) Reanalysis 1 dataset, which includes 925-hPa temperature. Petoukhov et al. (2008) examined the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) dataset, focusing on skewness and mixed third-order moments for a variety of atmospheric variables, also including temperature. Shen et al. (2011) examined changes in the higher-order moments of surface air temperature (SAT) of the contiguous United States from 1901 to 2000 and found decreasing second- and third-order moments and increasing fourth-order moments. Donat and Alexander (2012) calculated the differences between the moment statistics (up to skewness) of two 30-yr climatologies, 1951–80 and 1981–2010, for the daily maximum and minimum temperatures in the Hadley Centre Global Historical Climatology Network–Daily (HadGHCND), a gridded dataset interpolated from the station data used here.

Quite often, approximate Gaussianity is assumed for simplicity in studies of geophysical systems (e.g., in independent EOF analysis for physical interpretations and the independent data and model errors in Kalman filter data assimilation). However, non-Gaussian variability in observations and its prediction from ensembles of nonlinear dynamical systems perturbed with Gaussian white noise (Palmer 1993, and references therein) has been known for quite some time. Principally, non-Gaussian ensemble statistics necessarily arise from nonlinear deterministic dynamics or a non-Gaussian noise term when considering a stochastic forcing. Recent theory in empirical stochastic–dynamic modeling has adapted to capture non-Gaussianity through the form of multiplicative noise (Berner 2005; Sura et al. 2005) and has been shown to spatially approximate statistics from idealized atmospheric models (Sardeshmukh and Sura 2009). Nonzero skewness measures the asymmetry of a PDF and is sufficient to demonstrate non-Gaussianity, and the non-Gaussian structure of climate data has been documented for several climate datasets as aforementioned.

Two major pitfalls are detrimental to the study of high-frequency probabilistic climate variability. First, temporal averages result in Gaussian anomaly distributions asymptotically due to the central limit theorem. For many variables and locations, monthly or longer time averages are approximately Gaussian (Stephenson et al. 2004). Second, spatial correlation scales vary globally, necessitating that close attention be paid to the aggregation of data over homogeneous regions. Averaging or aggregating data over large scales can sully the underlying statistics.

The current paper studies not only the first two statistical moments of SAT variability, but also the higher-order moments, whose changes may imply deviations in the structure of the governing dynamical system. These changes of higher statistical moments, which impart PDF changes, can help diagnose dynamic variability and quantify variations of many applications-oriented indices including frost days, heating and cooling degree days, corn-heat-units, and heat stress index. Our work provides a comprehensive documentation of non-Gaussian structure for Northern Hemisphere SAT using the ground station observations of the Global Historical Climatology Network–Daily (GHCND) from 1 January 1950 to 31 December 2010. Information on non-Gaussianity can help explain why certain EOF modes are not independent and/or have little physical meaning (Monahan et al. 2009). The climatology may practically be used as benchmark distribution characteristics for validating model output mean, variance, and extremes. They can also be used for a method of moments estimation for chosen distributions and to derive inferences about the probabilities of specific climatic extremes (e.g., Rahmstorf and Coumou 2011). While this dataset suffers from incomplete spatial and temporal coverage, no a priori spatial or temporal averaging has been performed, which reduces the chance of systematic statistical bias and increases the chance of the results’ fidelity to truth. No attempt at assessing field significance about the spatial average of a large region is made since any demarcation would be inherently arbitrary without a more complete assessment of climatic regions based on daily PDFs or other high-frequency regime measures.

This paper is organized as follows. Section 2 discusses the observational dataset and its preprocessing steps, followed by a thorough explanation of further processing methodology. Section 3 describes data quality measures, statistical characterizations, regression methods, and PDF approximations. Section 4 includes topical regional analysis of results derived directly from GHCND data for both the climatology and trends. A discussion regarding the spatial patterns of climate change is supplemented with two examples to facilitate understanding and to demonstrate the use of these metrics for regional risk analysis. Section 5 concludes with a summary of findings and suggestions for future work.

2. GHCND data and preprocessing

GHCND contains in situ observational data from over 80 000 stations in 180 countries and territories worldwide (Menne et al. 2012). Recorded variables commonly include daily maximum Tmax and minimum temperature Tmin, total daily precipitation, snowfall, and snow depth. In this study, Tmax and Tmin records from approximately 4000 stations distributed across the Northern Hemisphere are analyzed. Each of these station records must pass both the GHCND consistency and quality assurance (QA) measures, outlined below, as well as record completeness measures from 1950 to 2010, outlined in section 3.

The GHCND source data have been compiled with the goal of maximizing spatial coverage. Records from numerous agencies have been collected and quality controlled, but not homogenized. Before new data are added to the dataset, they are first cross-checked with included data to confirm record uniqueness, spatial consistency, and temporal consistency, as outlined in Menne et al. (2012), which reduces the chances of redundancy. Each record is then subject to 19 QA measures and flagged with an appropriate tag so that data can be utilized or discarded as necessary by a user (Durre et al. 2010). Toward the goal of capturing trends in the moments of daily anomalies aggregated over many seasons, each flagged temperature datum is viewed as questionable and removed from the analysis. It is assumed that each seasonal daily distribution is not undersampled over a 10-yr period, obviating the need to retain questionable data for the sake of robust statistical moment estimation.

The daily mean SAT anomaly, denoted by Tavg, for this study is defined as the average of Tmax and Tmin. First, a running average of five days is performed on Tmax and Tmin and seasonal cycles are computed from the 1961–90 smoothed values, producing the 30-yr daily Tmax and Tmin climatologies, which are then removed from the GHCND Tmax and Tmin data to produce Tmax and Tmin anomalies. Only stations that include at least 20 years of data (7300 days) over the 1961–90 period are included in this study. This anomaly calculation procedure follows Alexander et al. (2006) and Caesar et al. (2006) to facilitate intercomparison with other works. The station density per 3.75° longitude by 2.5° latitude grid cell is shown in Fig. 1 for stations that are at least 75% complete (i.e., 15 yr of data in each 20-yr period).

Fig. 1.

GHCND station density (number of stations per grid cell) over the regression period, separated into 20-yr intervals. Stations accounted for are those that are at least 75% complete over the 20-yr interval.

Fig. 1.

GHCND station density (number of stations per grid cell) over the regression period, separated into 20-yr intervals. Stations accounted for are those that are at least 75% complete over the 20-yr interval.

3. Methodology

a. Record completeness measures

Before each station’s Tavg anomaly record is analyzed, it is subject to record completeness checks to ensure robust statistical moment estimation for each 10-yr period as well as for completeness over the 1950–2010 regression period. As a preliminary measure, each station record is required to span the entire regression period; in general, for a given station, data quality and completeness (inferred by number of missing days) improves nearer to the present, so each record is required to include data for 1 January 1950 or before, but is not explicitly required to include data for 31 December 2010 or after, since many stations lack the most recent values (likely due to unreporting) but are otherwise of high quality for the most recent decade. To ensure an adequate sampling of the regression period, each station record that is less than 75% complete is removed.

It is assumed here that SAT anomalies are approximately interseasonally piecewise stationary, but maintain interannual nonstationarity (following Shen et al. 2011). The data from each station Tavg anomaly record are thus divided into four seasonal records, for December–February (DJF), March–May (MAM), June–August (JJA), and September–November (SON). Each seasonal record is then subject to a running 10-yr data-aggregating moving time window so that the first bin contains the 1950–59 days, the second bin the 1951–60 days, and so forth. Thus, we have used 52 bins, beginning with 1950–59 and ending with 2001–10. A 10-yr bin is selected to smooth the fluctuations of yearly anomaly statistics arising from the El Niño–Southern Oscillation (ENSO; see Stefanova et al. 2013), to yield an adequate sample size for robust estimation of summary statistics, to enable a robust trend detection over the relatively short regression period, and to help satisfy the statistical assumptions associated with the regression method.

The first four statistical moments for each 10-yr data bin that contains at least a 75% complete record, for each season and each station, are calculated. Out of these seasonal records, only records that contain summary statistics for 75% of the regression period (i.e., at least 39 out of the 52 binned 10-yr periods) are included in the trend analysis. This regression period completeness measure is in addition to the completeness measure applied to the daily data. The 75% completeness measures were selected as an intermediate level of strictness and for simplicity. The results are insensitive to this choice with values from 50% to 90% tested.

b. Statistical moments

Data in the binned 10-yr periods are used to calculate the first four moments, as follows. The expectation of each 10-yr seasonal distribution is characterized by the biased arithmetic sample mean

 
formula

for the jth station for time index t with nonmissing Xj,t contained within the 10-yr distribution. The biased estimator, a simple average, is used since the sample size τ is over 600, large enough so that the bias is effectively reduced. In comparison, the sample size for monthly or yearly data may be small and the unbiased estimator may be more appropriate.

To assess changes in the second moment, the sample standard deviation

 
formula

for the jth station for time index t contained within the 10-yr distribution is used. The standard deviation is chosen so that the units of the first two moments are comparable.

The third and fourth moments are characterized by the sample skewness and sample (excess) kurtosis, respectively:

 
formula
 
formula

For a Gaussian distribution, both the skewness and excess kurtosis are zero. It is customary to leave out “excess” in discussions of kurtosis (as will be done in this paper).

Skewness and kurtosis are considered descriptors of extreme probabilistic behavior in that they help define the tail behavior of a PDF. Skewness is a measure of the degree of asymmetry in a distribution. Positive (negative) skewness indicates the presence of a long right (left) tail; for a unimodal distribution this manifests itself as a more positive (negative) mean than mode, and generally indicates that equally probable high-valued (low-valued) events are more extreme than low-valued (high-valued) events. Conversely, this implies that equally extreme high-valued (low-valued) events are more (less) probable than their low-valued (high-valued) counterparts.

Kurtosis is a measure of the peakedness of a distribution relative to a Gaussian distribution. Simultaneously, kurtosis is also a measure of the tail thickness relative to a Gaussian distribution, since the area under a PDF must integrate to unity. Geometrically, distributions with positive kurtosis (and zero skewness in this case) will show greater density around the mean and tails than a Gaussian and are commonly considered “fat”-tailed, indicating higher probabilities for very large anomalies.

c. Generalized linear regression model

For dynamic systems, nonzero autocorrelation structure leads to a reduced effective number of observations. In this case, midlatitude atmospheric degrees of freedom can be considered largely controlled by the propagation of atmospheric Rossby waves (Perron and Sura 2013). Oceanic variables and climate oscillations are more complex with decorrelation times on the order of years to decades. The methodology described above introduces another form of autocorrelation in that each binned moment observation is a function of each seasonal/daily observation of the previous 10 years. For the mean, this autocorrelation is analogous to that induced by a low-pass running mean filter. Within the context of the simple linear regression model, this leads to overconfident estimates of the slope β (Wilks 2011). Common methods to estimate confidence intervals in autocorrelated time series trend analysis involve the estimation of a decorrelation time and subsequent adjustment of confidence intervals, bootstrap sampling methods, nonparametrics, and/or assumptions of a particular autocorrelation structure such as an autoregressive (AR) model. In this analysis, moments are assumed to have an AR(1) structure of the form

 
formula
 
formula
 
formula

for the seasonal moment M at station j for year t, with y-intercept α, slope β, error term u, AR(1) coefficient ρ, and random normal noise ε of zero mean and σ standard deviation. Further details are discussed in Box et al. (1994).

4. Results

a. Seasonal climatology of non-Gaussian Tavg statistics

The left panels of Figs. 25 show each station’s 1961–90 moment climatology for MAM, JJA, SON, and DJF, respectively. The mean is presented as the seasonal mean over the climatology period, and the higher moments are calculated from the distribution of daily anomalies. Each moment subplot from each season can be analyzed individually to assess the effect of regional geography on that particular distribution characteristic. Seasonal dependence of regional distribution characteristics can be assessed by comparing the appropriate subplots from each season, and the four-moment subplots from each individual season can be taken in conjunction to draw conclusions about the behavior of Tavg variability relative to regional geography. Grid cells corresponding to 3.75° longitude × 2.5° latitude bins on skewness and kurtosis plots are hatched if at least 50% of the stations contained within that particular grid cell are significantly different from 0 at a 5% significance level, based on the standard errors (SE) of skewness and kurtosis (Brooks and Carruthers 1953):

 
formula
 
formula

where Neff is the effective number of sampled days. Perron and Sura (2013) show that these standard errors are approximately valid for the weakly non-Gaussian distributions observed in the atmosphere. We assume here that Neff = N/7 based on the reasoning put forth by Perron and Sura (2013), namely that midlatitude Rossby waves traverse the globe in approximately one week. For all stations, we take N equal to 20 years’ worth of seasonal data, the minimum to be considered in the climatology.

Fig. 2.

(a),(c),(e),(g) Climatological moment statistics and (b),(d),(f),(h) trends over the 1950–2010 regression period for MAM. The mean in (a) is calculated from the seasonal climatology, whereby the variability statistics in (c), (e), and (f) are calculated from the anomalies. Individual subplots for each season are plotted on the same color scale to facilitate seasonal intercomparison. Skewness and kurtosis grid cells are hatched if at least 50% of the contained stations are significantly different from 0 at a 5% significance level using standard error formulas [Eqs. (6a) and (6b)]. Trend grid cells are hatched if at least 50% of the contained stations show a significant trend at 5% one-tailed significance level, taken here to be “robust.” Contours in (a) and (c) show climatological 250-hPa geopotential height with 200-m spacing; contours in (e) show climatological sea level pressure with 5-hPa spacing, both from the NCEP–NCAR Reanalysis 1 (Kalnay et al. 1996). Extremes and outliers are truncated slightly to enhance interpretability.

Fig. 2.

(a),(c),(e),(g) Climatological moment statistics and (b),(d),(f),(h) trends over the 1950–2010 regression period for MAM. The mean in (a) is calculated from the seasonal climatology, whereby the variability statistics in (c), (e), and (f) are calculated from the anomalies. Individual subplots for each season are plotted on the same color scale to facilitate seasonal intercomparison. Skewness and kurtosis grid cells are hatched if at least 50% of the contained stations are significantly different from 0 at a 5% significance level using standard error formulas [Eqs. (6a) and (6b)]. Trend grid cells are hatched if at least 50% of the contained stations show a significant trend at 5% one-tailed significance level, taken here to be “robust.” Contours in (a) and (c) show climatological 250-hPa geopotential height with 200-m spacing; contours in (e) show climatological sea level pressure with 5-hPa spacing, both from the NCEP–NCAR Reanalysis 1 (Kalnay et al. 1996). Extremes and outliers are truncated slightly to enhance interpretability.

Fig. 3.

As in Fig. 2, but for JJA.

Fig. 3.

As in Fig. 2, but for JJA.

Fig. 4.

As in Fig. 2, but for SON.

Fig. 4.

As in Fig. 2, but for SON.

Fig. 5.

As in Fig. 2, but for DJF.

Fig. 5.

As in Fig. 2, but for DJF.

The most profound (and possibly most subtle) observation is that for each moment, which characterizes an aspect of probability distributions in each individual season, the spatial structure is smooth at the synoptic scale, but the emergent patterns over the Northern Hemisphere in each of the moments are different. The implications of this are twofold. First, for any particular region the contained stations’ anomaly PDFs are the same or similar regardless of local elevation or other geography. This is another justification of using anomalies to assess climate changes. Second, each climatological statistical moment has its own coherent spatial pattern and the direction and magnitude of the trend of each moment also has a coherent spatial pattern. These coherent spatial patterns provide the potential for skill in both spatial and temporal interpolation (and perhaps temporal extrapolation). Moment interpolation, as opposed to time series interpolation, has the benefit of preserving probability distribution characteristics and hence extremes. Comparing the patterns in moments contained within each season, the negative linear relationship between mean and standard deviation noted earlier by Robeson (2002) is evident (see Fig. 6). Cloud cover, cloud type, humidity, and wind speed are key controllers of SAT, and their variability in time is impacted by the variability of synoptic conditions. Robeson (2002) suggests generally that regions with higher average temperatures are less likely to experience highly variable synoptic-scale conditions, such as the frequent passage of cold fronts, which lead to increased standard deviations in SAT. While a parabolic relationship between skewness and kurtosis is expected for a diffusive system, a quite precise inequality for a random variable governed by a Fokker–Planck equation forced by multiplicative noise is kurt ≥ (3/2)skew2 + constant (Sura and Sardeshmukh 2008), which is also valid for our data as shown in Fig. 7. This suggests that Tavg might also be modeled by such a stochastic system. The offset by a small negative constant is discussed in detail in Sura and Perron (2010). The dynamical generator of these higher-order moments in the climate system is thought to be local adiabatic turbulence (Sardeshmukh and Sura 2009). These relationships do not hold over coastal regions where near-shore stations have a strong marine influence.

Fig. 6.

The 1961–90 climatological mean and standard deviation (std) for DJF for GHCND stations included in the climatology. Other seasons are qualitatively similar. An empirical linear fit is plotted in red.

Fig. 6.

The 1961–90 climatological mean and standard deviation (std) for DJF for GHCND stations included in the climatology. Other seasons are qualitatively similar. An empirical linear fit is plotted in red.

Fig. 7.

The 1961–90 climatological skewness and kurtosis for DJF for GHCND stations included in the climatology. Other seasons are qualitatively similar, except with varying levels of vertical shift. The kurtosis–skewness inequality is plotted in blue. An empirical parabolic fit is plotted in red.

Fig. 7.

The 1961–90 climatological skewness and kurtosis for DJF for GHCND stations included in the climatology. Other seasons are qualitatively similar, except with varying levels of vertical shift. The kurtosis–skewness inequality is plotted in blue. An empirical parabolic fit is plotted in red.

Unsurprisingly, climatological means of the raw field (Figs. 2a, 3a, 4a, and 5a), as opposed to the anomaly, exhibit a strong meridional gradient where higher (lower) latitudes are colder (warmer) which can be explained simply by top-of-atmosphere (TOA) insolation; deviations from the TOA insolation pattern are dynamically induced, and mirror to a large extent the climatological jet stream path. Winter months produce colder means nearer to the pole, whereas the tropics remain relatively constant year-round. Regions with a more maritime influence (e.g., western Europe) whose temperature fluctuations are strongly damped exhibit much lower standard deviations relative to other areas within the same zone. MAM, SON, and DJF (Figs. 2c, 4c, and 5c) show a similar standard deviation pattern with slightly elevated values above 45°N during DJF. Seasonal cycles of standard deviations are evident over the entire NH with maxima in DJF (Fig. 5c) and minima in JJA months (Fig. 3c), which is consistent with our understanding of the mean and standard deviation relationship.

Although the large spatial patterns of the mean and standard deviation in Figs. 2a,c, 3a,c, 4a,c, and 5a,c are not surprising, the smaller spatial scales of skewness and kurtosis seem striking (Figs. 2f,h, 3f,h, 4f,h, and 5f,h). Skewness is phase shifted from the locations of semipermanent low-level high and low pressure systems, as originally proposed by Loikith and Broccoli (2012), suggesting that large-scale circulation patterns may play a role in modulating the likelihood or magnitude of asymmetrical temperature extremes. For example, during DJF (Fig. 5e), the areas east (west) of the Canadian and Siberian highs show positive (negative) skewness, and areas east (west) of the Aleutian and Icelandic lows show negative (positive) skewness. These results are largely consistent with those of Perron and Sura (2013). In other areas and at different times, more local effects appear to be at play. These effects include modulation of extreme distribution characteristics by large bodies of water and, perhaps in some places, atmospheric conditions at elevation. Coastal areas, particularly near the west and north sides of continents, exhibit positive skewness, which becomes both more prominent and extends farther poleward during the warmer months. Skewness is most prominently negative along the Great Continental Divide in the Americas during cold seasons (Fig. 5e). Similarly, kurtosis values are most prominently positive in these areas at the same times. It is unknown exactly how the synoptic-scale wavelike patterns of skewness are dynamically linked to nonlinear atmospheric dynamics [see, e.g., chapters 6 and 7 in Holton (2004)]. Northern Hemisphere kurtosis reaches a maximum in MAM and SON months (Figs. 2g and 4g) and is predominantly positive below 30°N year-round. Soil moisture levels regulate surface air temperature extremes (Durre et al. 2010, and references therein), and it qualitatively appears that seasonal precipitation [see the Global Precipitation Climatology Project (GPCP) monthly climatology: http://jisao.washington.edu/data/gpcp/#climatology] plays an important role in the extremes and contribute to skewness and kurtosis values (particularly in the American northwest and Southeast Asia), but the exact relationship is unclear and left to future research.

b. The 1950–2010 temporal trends of statistical moments

The right panels of Figs. 25 show trends of the first four moments for each season in a 61-yr period from 1950 to 2010. Each station’s trend is plotted individually. Grid cells corresponding to 3.75° longitude × 2.5° latitude bins are hatched if at least 50% of the stations contained within that particular grid cell are significant at a 5% significance level. At this p value, hatched regions will be termed “robust.”

Positive trends in the mean for the 61-yr regression period are widespread over the Northern Hemisphere, with many regions exhibiting coherent robust trends, particularly over Eurasia during MAM (Fig. 2b). One notable exception is the so called “warming hole” located over the central southwest United States (Portmann et al. 2009), which shows robust weak cooling and is observable predominantly during JJA and DJF (Figs. 3b and 5b). Warming trends are strongest over Eurasia during MAM and DJF at about 0.4°C decade−1 (Figs. 2b and 5b). During MAM (Fig. 2b), there is also a weak cooling trend over northern Mexico.

Trends in standard deviation are seasonally variable, although broadly most of the robust and regionally coherent trends tend to be negative over the Northern Hemisphere. During MAM (Fig. 2d), the dominant robust features include positive trends over the Iberian Peninsula and negative trends over the southeastern United States, the Pacific Northwest region of the United States and Canada, Europe, western Russia, Central Asia, and noncoastal Siberia. During JJA (Fig. 3d), robust trends include positive trends over western Europe and negative trends over the United States (except in the northwest), western Canada, Alaska, and inland Siberia. During SON (Fig. 4d), there are no coherent, robust, regional, positive trends in standard deviation, but negative robust trends are observed over the contiguous United States (except the Rocky Mountains region), Canada, Alaska, Great Britain, eastern Europe, and Central Asia. During DJF (Fig. 5d), robust positive trends are over the Ural Mountains and robust negative trends are over most of the contiguous United States (except the Midwest region), Canada, Europe, Siberia, and Central Asia.

During MAM (Fig. 2f), positive robust trends of skewness are observed over the Ural Mountains and Central Asia, with negative robust trends over the Iberian Peninsula, eastern Europe, much of Siberia, and Mexico. During JJA (Fig. 3f), robust positive trends exist over parts of eastern Europe, with widespread robust negative trends over Alaska, central Canada, the Pacific Northwest, Appalachia, western Europe, the Ural Mountains, Central Asia, and most of Siberia. During SON (Fig. 4f), robust trends include localized positive trends over much of Russia and negative trends over the west coast of the United States and Canada and in eastern Siberia. During DJF (Fig. 5f), positive robust trends are observed over western Europe, eastern Siberia, and the Ural Mountains, with negative robust trends over Central Asia, central Siberia, northern Canada, Alaska, California, and Mexico. It remains to be investigated whether these trend patterns can be linked to changing climate dynamics.

Robust trends for kurtosis over each season also tend to be predominantly negative. During MAM (Fig. 2h), negative robust trends are observed over the western United States and Canada, Alaska, the Iberian Peninsula, eastern Europe, Russia and Siberia, and Central Asia. During JJA (Fig. 3h) negative trend regions include central Canada, the United States (except for the West Coast), western Europe, the Balkans, Central Asia, the Ural Mountains, and eastern Siberia. During SON (Fig. 4h), negative trends are observed over the United States (except the Midwest), eastern Europe, Russia and Siberia, and Central Asia. There are no robust positive trend regions over MAM, JJA, or SON. During DJF (Fig. 5h), robust positive trend regions are observed over central Canada and Siberia, with negative trends over the United States (except the Midwest), mainland Europe, the Ural Mountains, and western Russia.

c. Discussion of the results

As shown previously in section 4a, the climatological moments, and thus the PDFs that they describe, are spatially coherent and smoothly transition at synoptic scales (left-hand side of Figs. 25). This suggests that after removing a seasonal cycle, the underlying anomaly PDFs of individual GHCND stations over large areas are remarkably similar. The spatial patterns of each moment for each season, however, are different. There is a markedly negative linear relationship between the mean and standard deviation and a quadratic (parabolic) relationship between skewness and kurtosis, as shown in Figs. 6 and 7. This observed skewness–kurtosis relationship is predicted through moment closure in stochastic–dynamic theory (Sura 2011) and suggests that surface temperature may be appropriately modeled with a linear stochastic–dynamic system with correlated additive and multiplicative noise.

Meridionally, the climatological mean and standard deviation exhibit a strong gradient, which has been widely discussed in existing literature. Contour lines of 250-hPa geopotential height Z250 are shown in Figs. 2a,c, 3a,c, 4a,c, and 5a,c and appear to be conspicuously aligned with zonal waves of mean and standard deviation. These zonal waves are particularly prominent during DJF (Figs. 5a,c) and are coincident with the locations of the greatest zonal circulation gradients. Northern Hemisphere extratropical circulation gradients are smaller and exhibit weaker wave formation during JJA (Figs. 3a,c), resulting in reduced zonal fluctuations in mean and standard deviation during these months.

Recent studies have investigated the linkage between 500-hPa geopotential height Z500 and surface temperature extremes and suggest that persistent Z500 anomalies project strongly (Meehl and Tebaldi 2004) and, particularly in DJF, linearly (Loikith and Broccoli 2012) on SAT extremes. Loikith and Broccoli (2012) also suggest that regional variations in skewness may point to regions whose temperature extremes are dominated by local or regional circulation dynamics. During DJF, when semipermanent high and low pressure systems are strongest, geographical skew patterns are coherent on synoptic scales, suggesting that DJF extremes are dominated by large-scale dynamic regimes, consistent with Loikith and Broccoli (2012). During JJA, geographical skew patterns are much more regional and resemble climatological relative humidity (Dai 2006), suggesting that JJA temperature extremes are dominated by local conditions such as soil moisture content (e.g., Hirschi et al. 2011), also consistent with Loikith and Broccoli (2012). During transition months, it is clear that there is a mixture of both synoptic and regional impacts on Tavg. Patterns in kurtosis may result by the linear mechanism suggested by Loikith and Broccoli (2012) from distributions of Z500 whose DJF kurtosis climatologies (Perron and Sura 2013) resemble the kurtosis climatologies here or may also result from local nonlinear interactions, exemplified by the strong positive kurtosis region in the northwestern United States during DJF.

We additionally investigate the trends of these moments over a 1950–2010 regression period against an AR(1) model. Trends in the mean show familiar global warming trend patterns for each season that are compatible with estimates of recent climate change (Solomon et al. 2007). Warming trends during MAM are the most suggestive, particular over Eurasia. Robust trends in the standard deviation are predominantly negative globally, which is also consistent with recent research (Karl et al. 1995; Michaels et al. 1998; Shen et al. 2011). Although standard deviation trends during MAM and JJA are more positive overall relative to the other seasons and show spatially coherent positive trends over large areas, the positive standard deviation trends tend not to be robust at grid scale whereas the negative trends are. Trends in skewness and kurtosis show robustness at grid scale and exhibit spatial coherence as well. Over the United States, trends in skewness are predominately positive whereas kurtosis trends are predominantly negative, which is in agreement with the literature (Shen et al. 2011). Here we include regional details, which make the results more useful as benchmark data for climate change assessment via other means, such as remote sensing and mathematical modeling. There is no other comparable literature for other regions.

Karl and Katz (2012) briefly discussed how sampling Gaussian distributed variability over nonstationary time periods when means and standard deviations have changed can in fact lead to nonzero, biased estimates of skewness. If such nonstationary time series have nonlinear trends in the first two moments, it could also lead to biases in the estimates of higher-moment trends as well. We investigated this possibility briefly through Monte Carlo simulation (not shown). Biases in skewness and kurtosis estimates arising from linear trends in the mean and standard deviation of the order of magnitude as those observed in the atmosphere are at least an order of magnitude smaller than observed Tavg skewness and kurtosis. Additionally, nonlinear trends in the mean and standard deviation lead to trends in skewness and kurtosis that are also at least an order of magnitude smaller than those observed for all reasonable mean and standard deviation temporal trends.

No attempt to quantitatively assess trend significance at larger than grid scale is made here, since it is unclear which stations can be appropriately considered homogeneous representations of a climate (trend) region. For example, a previous study by Shen et al. (2011) estimated higher-moment trends over the contiguous United States; however, it is apparent from this station-scale analysis that both the climatological PDFs and moment trends for stations within the United States are not a homogeneous sample. It is also clear, though, that taken individually there is strong regionalization for trends in each moment (whose regions may or may not coincide with the trend regions for other moments), suggesting that the underlying dynamic mechanisms responsible for these PDF characteristics also have changes on regional scales, globally. One exceptional note is that, in general, during winter months, skewness trends tend to oppose their climatological values, which is perhaps an indication of weakening or shifting climatological high and low pressure systems that manifest themselves in skewness statistics.

Information regarding entire distributions, both the bulk of observations and the extremes, is critical for understanding and adapting to regional climate change, and high-frequency variability (e.g., daily statistics) driven by low-frequency climate modes (e.g., ENSO phase) is increasingly in the public interest. Recent climate extremes, such as the Russian heat wave and European blocking event of 2010, has piqued interest in regional climate extremes and global teleconnections (e.g., Trenberth and Fasullo 2012; Rahmstorf and Coumou 2011), as well as the role of natural variability in the context of high-frequency extreme events (e.g., Dole et al. 2011; Matsueda 2011). In Fig. 8, we choose weather stations located near Moscow, Russia, during JJA to illustrate how changes in moment statistics are manifested through shifts in entire PDFs. Figure 8 shows that the PDFs of Tavg from each station located within the grid cell for the climatological period are nearly identical, as suggested by the smoothness of climatological moments in space. During JJA, the climatological anomaly PDFs are nearly Gaussian, as indicated by skewness and kurtosis near zero (Figs. 3e,g). Over the 1950–2010 regression period, this region shows positive trends in mean and skewness (Figs. 3b,f), which are indicated on Fig. 8 as a positive shift of the bulk of the distribution as well as a lengthening and thickening of the positive anomaly tail, which is particularly noticeable at the extremes. The cumulative probability density of large positive anomalies over JJA has noticeably increased each decade for the last four decades, indicating that the probability of large positive anomalies has increased, particularly during the 2001–10 period.

Fig. 8.

The 1961–90 climatological PDFs for each station contained in the grid cell covering Moscow (gray) overlaid with decadal PDFs determined by aggregating data from each station in the contained grid cell for each interval (colored). Decades are marked as the last year in the decadal period (e.g., 2010 represents the 2001–10 interval).

Fig. 8.

The 1961–90 climatological PDFs for each station contained in the grid cell covering Moscow (gray) overlaid with decadal PDFs determined by aggregating data from each station in the contained grid cell for each interval (colored). Decades are marked as the last year in the decadal period (e.g., 2010 represents the 2001–10 interval).

Similarly, a recent reduction in water levels in the Colorado River basin has prompted concern regarding recorded and projected climate change over the western United States. Studies have shown that while “business as usual” climate projections show slightly decreased precipitation levels (Christensen et al. 2004), temperature changes and the early onset of spring conditions, which result in an increased rain–snow ratio and contribute to early season snowmelt, have the largest impact on the basin (Cayan et al. 2001; Stewart et al. 2004). Examining the climatological PDFs of stations in the region (Fig. 9) again reveals similar PDFs for each of the contained stations. The Colorado stations have a slightly elongated negative anomaly tail, which is manifested in the climatological statistics as negative skewness over the region (Fig. 3e). Over the 1950–2010 regression period, the Colorado region exhibits a nonrobust warming trend, with a robust negative trend in kurtosis. Neither trend is particularly distinguishable in decadal PDFs (Fig. 9) as shifts in kurtosis are subtle and difficult to discern; however, the most recent decadal PDF does display a distinct linear, positive shift, and PDFs for the last two decades have become less peaked, which is consistent with a negative shift in kurtosis.

Fig. 9.

As in Fig. 8, but for the region over central Colorado.

Fig. 9.

As in Fig. 8, but for the region over central Colorado.

5. Conclusions

We present non-Gaussian 1961–90 climatology statistics for GHCND station Tavg data in the Northern Hemisphere for each season, as well as trends in these statistics over a 1950–2010 regression period. The relevant climatology results for gridded data mentioned in the introduction resemble those derived from station data in this paper. Non-Gaussian probability functions of SAT reflect strongly on the probabilities of both high and low temperature extremes (Ruff and Neelin 2012), and proper assessment of climate distribution tail behavior is critical for proper model and reanalysis verification and climate risk assessment under climate change. Climatological moments show coherence on regional–synoptic scales and are seasonally dependent, suggesting that anomaly PDFs are quite similar at regional scales but exhibit a strong seasonal dependence. The mean and standard deviation show synoptic coherence and a strong meridional gradient and are modulated zonally by circulation patterns, while skewness and kurtosis fluctuate on smaller scales, implicating regional dynamics in the modulation of temperature extremes, particularly for summer months. Trends in moments show large regionalization with robust changes at grid scale for all moments and all seasons. For mean and standard deviation, trends are coherent on synoptic scales and show robust warming in most areas during each season and reduced variability seen in most areas during most seasons. Trends in skewness and kurtosis vary on smaller regional scales and suggest changes to the regional dynamics that modulate higher-order statistics at those scales. We further present as case studies the climatological and decadal PDFs for the regions around Moscow and Colorado to illustrate the homogeneity of regional station PDFs and how shifts in moments are represented in entire PDFs over the regression period.

Further work on the mechanisms of these shifts is needed to gain a physical understanding of the (changing) processes that manifest themselves in the (changing) statistical moments. If explicit dynamical links can be made to climatological moments, either through stochastic dynamics (e.g., multiplicative noise discussed previously) or through numerical models (e.g., the specific contribution of extreme weather events to variability statistics), trends in these moments might provide further insight regarding local, regional, and global climate change. Additionally, a greater understanding of the spatial correlation structure of station time series is needed so that biases in the moments arising from averaging can be anticipated and potentially corrected or avoided. Spatial averaging and/or interpolation results in reduced standard deviations, as well as skewness and kurtosis values closer to zero, as is expected by the central limit theorem. Finally, further work on assigning regions of climate change with respect to daily variability (e.g., through cluster analysis; Loikith et al. 2013) is needed so that spatial correlation in moment trends can be appropriately accounted for to provide regional estimates of climate change statistical significance.

Our moment method can be used to assess climate changes, diagnose climate dynamics, and validate models in terms of probability distributions. Toward the specific ends of examining extremes in observational data, Donat et al. (2013) have compiled the GHCNDEX, a global database of extremes indices compiled from GHCND. As for daily-observed data, we are also aware of the existence of the gridded version of GHCND data, HadGHCND. A detailed assessment of the differences between the station GHCND and gridded GHCND will be made in another study, which will particularly facilitate climate model validations when regarding climate model output at a grid point as a gridbox average.

Acknowledgments

This study was supported in part by the National Science Foundation (Award AGS-1015926). Shen was also supported by the National Oceanic and Atmospheric Administration (Award EL133E09SE4048), National Science Foundation (Award AGS-1015957), and U.S. Department of Energy (Award DE-SC002763). We greatly appreciate the efforts of two anonymous reviewers whose comments greatly improved the clarity of this paper. Cavanaugh would also like to acknowledge funding received from the Scripps Institution of Oceanography and thank Alexander Gershunov for his thoughtful contributions to the paper.

REFERENCES

REFERENCES
Alexander
,
L. V.
, and Coauthors
,
2006
:
Global observed changes in daily climate extremes of temperature and precipitation
.
J. Geophys. Res.
,
111
, D05109, doi:.
Berner
,
J.
,
2005
:
Linking nonlinearity and non-Gaussianity of planetary wave behavior by the Fokker–Planck equation
.
J. Atmos. Sci.
,
62
,
2098
2117
, doi:.
Box
,
G. E. P.
,
G. M.
Jenkins
, and
G. C.
Reinsel
,
1994
:
Time Series Analysis: Forecasting and Control
. 3rd ed. Prentice Hall, 598 pp.
Brooks
,
C. E. P.
, and
N.
Carruthers
,
1953
:
Handbook of Statistical Methods in Meteorology
. Her Majesty’s Stationery Office, 412 pp.
Caesar
,
J.
,
L.
Alexander
, and
R.
Vose
,
2006
:
Large-scale changes in observed daily maximum and minimum temperatures: Creation and analysis of a new gridded data set
.
J. Geophys. Res.
,
111
,
D05101
, doi:.
Cayan
,
D. R.
,
M. D.
Dettinger
,
S. A.
Kammerdiener
,,
J. M.
Caprio
, and
D. H.
Peterson
,
2001
:
Changes in the onset of spring in the western United States
.
Bull. Amer. Meteor. Soc.
,
82
,
399
415
, doi:.
Christensen
,
N. S.
,
A. W.
Wood
,
N.
Voisin
,
D. P.
Lettenmaier
, and
R. N.
Palmer
,
2004
:
The effects of climate change on the hydrology and water resources of the Colorado River basin
.
Climatic Change
,
62
,
337
363
, doi:.
Dai
,
A.
,
2006
:
Recent climatology, variability, and trends in global surface humidity
.
J. Climate
,
19
,
3589
3606
, doi:.
Dole
,
R.
, and Coauthors
,
2011
:
Was there a basis for anticipating the 2010 Russian heat wave?
Geophys. Res. Lett.
,
38
,
L06702
, doi:.
Donat
,
M. G.
, and
L. V.
Alexander
,
2012
:
The shifting probability distribution of global daytime and night-time temperatures
.
Geophys. Res. Lett.
,
39
,
L14707
, doi:.
Donat
,
M. G.
,
L. V.
Alexander
,
H.
Yang
,
I.
Durre
,
R.
Vose
, and
J.
Caesar
,
2013
:
Global land-based datasets for monitoring climatic extremes
.
Bull. Amer. Meteor. Soc.
,
94
,
997
1006
, doi:.
Durre
,
I.
,
M. J.
Menne
,
B. E.
Gleason
,
T. G.
Houston
, and
R. S.
Vose
,
2010
:
Comprehensive automated quality assurance of daily surface observations
.
J. Appl. Meteor. Climatol.
,
49
,
1615
1633
, doi:.
Gong
,
D.-Y.
, and
C.-H.
Ho
,
2004
:
Intra-seasonal variability of wintertime temperature over East Asia
.
Int. J. Climatol.
,
24
,
131
144
, doi:.
Hansen
,
J.
,
M.
Sato
, and
R.
Ruedy
,
2012
:
Perception of climate change
.
Proc. Natl. Acad. Sci. USA
,
109
,
E2415
E2423
, doi:.
Hirschi
,
M.
, and Coauthors
,
2011
:
Observational evidence for soil-moisture impact on hot extremes in southeastern Europe
.
Nat. Geosci.
,
4
,
17
21
, doi:.
Holton
,
J. R.
,
2004
:
An Introduction to Dynamic Meteorology
. 4th ed. Academic Press, 535 pp.
Kalnay
,
E.
, and Coauthors
,
1996
:
The NCEP/NCAR 40-Year Reanalysis Project
.
Bull. Amer. Meteor. Soc.
,
77
,
437
471
, doi:.
Karl
,
T. R.
, and
R. W.
Katz
,
2012
:
A new face for climate dice
.
Proc. Natl. Acad. Sci. USA
,
109
,
14 720
14 721
, doi:.
Karl
,
T. R.
,
R. W.
Knight
, and
N.
Plummer
,
1995
:
Trends in high-frequency climate variability in the twentieth century
.
Nature
,
377
,
217
220
, doi:.
Loikith
,
P. C.
, and
A. J.
Broccoli
,
2012
:
Characteristics of observed atmospheric circulation patterns associated with temperature extremes over North America
.
J. Climate
,
25
,
7266
7281
, doi:.
Loikith
,
P. C.
,
B. R.
Lintner
,
J.
Kim
,
H.
Lee
,
J. D.
Neelin
, and
D. E.
Waliser
,
2013
:
Classifying reanalysis surface temperature probability density functions (PDFs) over North America with cluster analysis
.
Geophys. Res. Lett.
,
40
,
3710
3714
, doi:.
Matsueda
,
M.
,
2011
:
Predictability of Euro-Russian blocking in summer of 2010
.
Geophys. Res. Lett.
,
38
,
L06801
, doi:.
Meehl
,
G. A.
, and
C.
Tebaldi
,
2004
:
More intense, more frequent, and longer lasting heat waves in the 21st century
.
Science
,
305
,
994
997
, doi:.
Meehl
,
G. A.
,
C.
Tebaldi
, and
D.
Nychka
,
2004
:
Changes in frost days in simulations of twentyfirst century climate
.
Climate Dyn.
,
23
,
495
511
, doi:.
Menne
,
M. J.
,
I.
Durre
,
R. S.
Vose
,
B. E.
Gleason
, and
T. G.
Houston
,
2012
:
An overview of the Global Historical Climatology Network-Daily database
.
J. Atmos. Oceanic Technol.
,
29
,
897
910
, doi:.
Michaels
,
P. J.
,
R. C.
Balling
Jr.
,
R. S.
Vose
, and
P. C.
Knappenberger
,
1998
:
Analysis of trends in the variability of daily and monthly historical temperature measurements
.
Climate Res.
,
10
,
27
33
, doi:.
Monahan
,
A. H.
,
J. C.
Fyfe
,
M. H. P.
Ambaum
,
D. B.
Stephenson
, and
G. R.
North
,
2009
:
Empirical orthogonal functions: The medium is the message
.
J. Climate
,
22
,
6501
6514
, doi:.
Palmer
,
T. N.
,
1993
:
Extended-range atmospheric prediction and the Lorenz model
.
Bull. Amer. Meteor. Soc.
,
74
,
49
65
, doi:.
Parker
,
D. E.
,
P. D.
Jones
,
C. K.
Folland
, and
A.
Bevan
,
1994
:
Interdecadal changes of surface temperature since the late nineteenth century
.
J. Geophys. Res.
,
99
,
14 373
14 399
, doi:.
Perron
,
M.
, and
P.
Sura
,
2013
:
Climatology of non-Gaussian atmospheric statistics
.
J. Climate
,
26
,
1063
1083
, doi:.
Petoukhov
,
V.
,
A. V.
Eliseev
,
R.
Klein
, and
H.
Oesterle
,
2008
:
On statistics of the free-troposphere synoptic component: An evaluation of skewnesses and mixed third-order moments contribution to the synoptic-scale dynamics and fluxes of heat and humidity
.
Tellus
,
60A
,
11
31
.
Portmann
,
R. W.
,
S.
Solomon
, and
G. C.
Hegerl
,
2009
:
Spatial and seasonal patterns in climate change, temperatures, and precipitation across the United States
.
Proc. Natl. Acad. Sci. USA
,
106
,
7324
7329
, doi:.
Rahmstorf
,
S.
, and
D.
Coumou
,
2011
:
Increase of extreme events in a warming world
.
Proc. Natl. Acad. Sci. USA
,
108
,
17 905
17 909
, doi:.
Robeson
,
S. M.
,
2002
:
Relationships between mean and standard deviation of air temperature: Implications for global warming
.
Climate Res.
,
22
,
205
213
, doi:.
Ruff
,
T. W.
, and
J. D.
Neelin
,
2012
:
Long tails in regional surface temperature probability distributions with implications for extremes under global warming
.
Geophys. Res. Lett.
,
39
, L04704, doi:.
Sardeshmukh
,
P. D.
, and
P.
Sura
,
2009
:
Reconciling non-Gaussian climate statistics with linear dynamics
.
J. Climate
,
22
,
1193
1207
, doi:.
Shen
,
S. S. P.
,
A. B.
Gurung
,
H.-S.
Oh
,
T.
Shu
, and
D. R.
Easterling
,
2011
:
The twentieth century contiguous US temperature changes indicated by daily data and higher statistical moments
.
Climatic Change
,
109
,
287
317
, doi:.
Sherwood
,
S. C.
, and
M.
Huber
,
2010
:
An adaptability limit to climate change due to heat stress
.
Proc. Natl. Acad. Sci. USA
,
107
,
9552
9555
, doi:.
Solomon
,
S. D.
,
M.
Qin
,
M.
Manning
,
Z.
Chen
,
M.
Marquis
,
K. B.
Averyt
,
M.
Tignor
, and
H. L.
Miller
, Eds.,
2007
: Climate Change 2007: The Physical Science Basis. Cambridge University Press, 996 pp.
Stefanova
,
L.
,
P.
Sura
, and
M.
Griffin
,
2013
:
Quantifying the non-Gaussianity of wintertime daily maximum and minimum temperatures in the Southeast
.
J. Climate
,
26
,
838
850
, doi:.
Stephenson
,
D. B.
,
A.
Hannachi
, and
A.
O’Neill
,
2004
:
On the existence of multiple climate regimes
.
Quart. J. Roy. Meteor. Soc.
,
130
,
583
605
, doi:.
Stewart
,
I. T.
,
D. R.
Cayan
, and
M. D.
Dettinger
,
2004
:
Changes in snowmelt runoff timing in western North America under a ‘business as usual’ climate change scenario
.
Climatic Change
,
62
,
217
232
, doi:.
Sura
,
P.
,
2011
:
A general perspective of extreme events in weather and climate
.
Atmos. Res.
,
101
,
1
21
, doi:.
Sura
,
P.
, and
P. D.
Sardeshmukh
,
2008
:
A global view of non-Gaussian SST variability
.
J. Phys. Oceanogr.
,
38
,
639
647
, doi:.
Sura
,
P.
, and
M.
Perron
,
2010
:
Extreme events and the general circulation: Observations and stochastic model dynamics
.
J. Atmos. Sci.
,
67
,
2785
2804
, doi:.
Sura
,
P.
,
M.
Newman
,
C.
Penland
, and
P.
Sardeshmukh
,
2005
:
Multiplicative noise and non-Gaussianity: A paradigm for atmospheric regimes?
J. Atmos. Sci.
,
62
,
1391
1409
, doi:.
Trenberth
,
K. E.
, and
J. T.
Fasullo
,
2012
:
Climate extremes and climate change: The Russian heat wave and other climate extremes of 2010
.
J. Geophys. Res.
,
117
, D17103, doi:.
Wackerly
,
D.
,
W.
Mendenhall
, and
R. L.
Scheaffer
,
2008
:
Mathematical Statistics with Applications
. 7th ed. Cengage Learning, 914 pp.
Wilks
,
D. S.
,
2011
:
Statistical Methods in the Atmospheric Sciences
. 2nd ed. Academic Press, 649 pp.