Changes in statistical distributions of sub-daily surface temperatures and wind speed

With the ongoing warming of the globe, it is important to quantify changes in the recent behaviour of extreme events given their impacts on human health, infrastructure and the natural environment. We use the sub-daily, multivariate, station-based HadISD dataset to study the changes in the statistical distributions of temperature, dew point temperature and wind speeds. Firstly, we use zonally averaged quantities to show that the lowest temperatures during both day and night are changing more rapidly than the highest, with the effect more pronounced in the northern high latitudes. Along with increases in the zonally averaged mean temperature, the standard deviation has decreased and the skew increased (increasing positive tail, decreasing negative tail) over the last 45 years, again with a stronger, more robust signal at higher latitudes. Changes in the distribution of dew point temperature are similar to those of temperature. However, changes in the distribution of wind speeds indicate a more rapid change at higher speeds than at lower.

Secondly, to assess in more detail the spatial distribution of changes as well as changes across seasons and hours of the day we study each station individually. For stations which show clear indications of change in the statistical moments, the higher the statistical moment, generally the more spatially heterogenous the patterns of change. The standard deviations of temperatures are increasing in a band stretching from Europe through China but are decreasing across North America and in the high northern latitudes, indicating broadening and narrowing of the distributions, respectively. Large seasonal differences are found in the change of standard deviations of temperatures over North America and eastern China. Temperatures in eastern Asia also have increasing skew in the winter in contrast to the remainder of the year. The dew point temperatures show smaller variation in all of the moments but similar patterns to the temperatures. For wind speeds, apart from the USA, standard deviations are decreasing across the world, indicating a decrease in variability.

Finally, we use quantile regression to show changes in the percentiles of distributions over time. These show an increase in high quantiles of temperature in eastern Europe during the summer and also in northern Europe for low quantiles in the winter, also indicating broadening and narrowing of the distributions, respectively. In North America, the largest changes are at the lower quantiles in northern latitudes for autumn and winter. Quantiles of dew point temperature are changing most in the autumn and winter, especially in the northern parts of Europe.

1 Introduction

The study of changes in the extremes of essential climate variables is vital given their impacts on human health, infrastructure, agriculture and the natural environment. The Intergovernmental Panel on Climate Change (IPCC) released a Special Report on Extremes (SREX; Field et al., 2012) and therein outlined three simple classes of the way changes in extremes occur: shifting the mean, increasing the variability and changing the symmetry. Past studies have indicated that there is uncertainty around how changes in the occurrence or intensity of extreme values arise. Are extremes changing due to changes in the location of the distribution mean with no change in the distribution shape (Griffiths et al., 2005; Simolo et al., 2011), or are changes in the shape of a distribution the primary driver of changes in extremes? If so, is the change in shape only the consequence of a change in the variance, or does it also arise from changes in higher-order moments (Della-Marta et al., 2007; Ballester et al., 2010)?

There are a whole host of ways to study climate extremes and determine how these have changed over the recent past. One common approach is to use climate extreme indices to characterise moderate extremes in time series of station data, for example those developed by the World Meteorological Organization (WMO) Commission for Climatology (CCl), World Climate Research Programme (WCRP), and the Joint Technical Commission for Oceanography and Marine Meteorology (JCOMM) Expert Team in Climate Change Detection and Indices (ETCCDI; Alexander et al., 2006). A number of datasets comprising these indices now exist, allowing for detailed investigations of past changes and comparisons to model and reanalysis fields (e.g. Caesar et al., 2006; Donat et al., 2013a, b). Another route uses extreme value theory to model the tails of the distributions and from these few points characterise the occurrence and intensity of extreme events, including those not yet observed in the modern data record. The “block maximum” approach models a set of e.g. annual maxima with a generalised extreme value distribution (e.g. Christidis et al., 2011). Alternatively, the “peaks over threshold” approach models all peaks over a fixed threshold with the generalised Pareto distribution (e.g. Brown et al., 2008).

A further method characterises the complete distribution using all available data to establish the causes of changes in climate extremes. The advantage of this approach is that it uses all of the available information and can be extended to more or less extremal parts of the distribution. A number of the ETCCDI temperature indices use the 10th and 90th percentile values, which only probe moderate extremes. A recalculation of these indices at 1st and 99th percentile values is possible, but given the diverse sources of data in the HadEX datasets (Donat et al., 2013b), this would be a large undertaking by the international community. Furthermore, as the climate continues to warm, the relative intensity or duration of extreme events against warmer average conditions may be of interest; for example, how warm is the relatively warmest 10 % of days now compared to the same fraction selected between 1961 and 1990?

Hence, a number of recent studies have investigated past changes in the distributions of observed land surface temperatures. Most have used the Global Historical Climate Network Daily (GHCND; Menne et al., 2012) or its gridded version (HadGHCND, temperature only; Caesar et al., 2006) as these provide daily maximum and minimum temperature values. Donat and Alexander (2012) compared two periods to show that the means of both temperature measures have shifted to warmer values. Along with larger changes in the minimum than the maximum temperatures, this has resulted in increases in the skew of the distributions. However, changes in standard deviation were less significant (at 10 %) and more heterogeneous. Over North America, Shen et al. (2011) show that the standard deviation and skewness were decreasing but that the kurtoses of the maximum temperatures were increasing, in contrast to the decrease in the kurtoses of minimum temperatures. A subsequent study investigated seasonal changes in the distribution moments over 1950–2010 (Cavanaugh and Shen, 2014). McKinnon et al. (2016) showed that for most stations between 1980 and 2015, over just the summer months, trends can mostly be explained by a shift in the distribution with no change in shape. Reanalyses (e.g. Huntingford et al., 2013; Gross et al., 2018) and coupled models (e.g. Lewis and King, 2017) have also been used in the study of changes in the shapes of current and future distributions of essential climate variables (ECVs).

In this study, we use sub-daily station data from a single dataset as opposed to the maximum and minimum daily temperatures (either station based or gridded). We note that by using a single dataset, we have not investigated the effect of dataset choice on our findings, something which can have a large effect (Gross et al., 2018). In order to compare with a number of previous studies we perform three different assessments to see how the distributions of temperature, dew point temperature and wind speed have changed in recent years. We outline the dataset and station selection criteria in Sect. 2. The three assessments are outlined in Sects. 3 to 5, with a summary in Sect. 6.

2 Input data

To study the changes in distributions related to climate extremes, we use the sub-daily HadISD dataset (Dunn et al., 2012, 2014, 2016; Dunn, 2019). This is updated annually and now covers the period 1931–2018 inclusive. Improved station selection and merging processes as well as minor changes in the quality control tests are features of version 3.0.0.2018f over version 1.0.x. This version contains 8139 unique station locations, with some of these being composited from individual records in the ISD (Integrated Surface Dataset; Smith et al., 2011). As a result of the long-standing issues surrounding the free sharing of observational climate data (Thorne et al., 2017), the bulk of these stations are in the Northern Hemisphere, concentrated in North America and Europe. In contrast to some of the previous studies outlined in the Introduction (e.g. Donat and Alexander, 2012), we do not create a gridded form of the HadISD for this assessment but retain the individual stations as input to the analyses, similar to e.g. Cavanaugh and Shen (2014) and McKinnon et al. (2016).

HadISD undergoes a homogeneity assessment using the pairwise homogeneity algorithm (PHA; Menne and Williams, 2009) to identify the number and size of inhomogeneities in its four main ECVs, including temperature and wind speed (Dunn et al., 2014). This used monthly averages derived from the sub-daily HadISD observations to compare station–pair differences to identify potential jumps in times eries, a process which is more effective for stations and variables with low variability. Although no adjustments are made, the outputs of the PHA allow stations that have few or small inhomogeneities to be selected, i.e. those with the most homogeneous records (see Dunn et al., 2014). Dunn et al. (2014) found that for temperature the best balance of station network against the exclusion of stations with the worst inhomogeneities occurred between 1 and 2 ^∘C. Therefore, in all of the investigations below, we select those stations that have fewer inhomogeneities, allowing jumps of up to 1 ^∘C (temperature and dew point) or 1 m s⁻¹ (as the distribution of the inhomogeneities in the wind speed data are similar; see Fig. 10 of Dunn et al., 2014). This balances retaining sufficient stations to study the behaviour across the globe with removing those with very large inhomogeneities.

As the selected stations still contain uncorrected inhomogeneities, we note that these will have affected the analysis of the statistical properties of the observations outlined herein. But the distributions of the estimated inhomogeneities as shown in Dunn et al. (2014, 2016) do not have a strongly non-zero mean, and the central portions retained by the approach above have no clear asymmetry. Therefore, as the following analyses use combinations of stations or look for contiguous regions of change that cross national (i.e. observing practice) boundaries, we do not expect large effects arising from the remaining inhomogeneities themselves.

Although HadISD v3.0.x contains data from 1931 onwards, there is a large increase in the number of stations from 1973 onwards (the reason why HadISD v1.0.x started then; see Fig. 2 in Dunn et al., 2016). Therefore, we restrict our analysis to the time period from 1 January 1973 unless otherwise stated. We also remove data from 29 February for ease of analysis.

3 Zonal distribution changes

We firstly look at the changes in the distributions of anomalies in temperatures, dew point temperatures and wind speeds for all stations in a specified latitude band on an annual basis. By combining stations together, the quantity of observations available to define the distribution is increased at the expense of spatial resolution. Donat and Alexander (2012) used the HadGHCND gridded daily temperature dataset (Caesar et al., 2006) to compare two 30-year periods of daily maximum and minimum temperatures (1951–1980 and 1981–2010) globally and for three smaller zonal regions (Northern Hemisphere, the tropics and Southern Hemisphere). Here we develop this method further using the sub-daily data available in HadISD as well as narrower 10^∘ zonal bands.

We select those stations that fall within each 10^∘ band. An approximation for the whole-hour time zone of the station is calculated from the station longitude, enabling the Coordinated Universal Time (UTC) observation time to be converted to local time. We then extract those portions of data that correspond to daytime and night-time, using the time periods 09:00–20:00 and 21:00–08:00, respectively, to ensure that the minimum temperatures are on the whole captured during the night-time portion and the maximum temperatures in the daytime. At high latitudes these definitions of daytime and night-time will vary in accuracy throughout the year but should split the 24 h day into the hours with the most and least insolation. There are few stations at the highest latitudes at which 24 h of day or night occur.

Climatologies for stations for each day of the year for the daytime and night-time data are calculated separately, requiring the equivalent of 15 years of 3-hourly observations over the 1981–2010 period. These climatologies are used to create anomalies from the sub-daily data. The creation of climate anomalies also restricts this analysis to stations that have consistent records within this recent climatology period.

Figure 1 Maps showing the locations of HadISD stations contributing to the zonal analysis along with the 10^∘ latitude bands for (a) temperature, (b) dew points and (c) wind speeds.

To investigate the changes in the distributions of these zonally averaged quantities over time, we split the data up into nine intervals each of 5 years in length (1974–1978, 1979–1983, ..., 2014–2018). By not including 1973, we can assess changes over a 45-year period. Furthermore, we only take those stations that have sufficient observations spread over the entire analysis period. In each interval, we require there to be more observations than the number equivalent to hourly observations for 1∕4 of the 5-year interval. Also, at least 3∕4 of the intervals require sufficient data for the station to be incorporated into the final distribution. This is to ensure that nonuniform distributions of observations within a station record (both within a 5-year interval and across the entire record) do not overly influence the final distributions. The stations which pass these selection criteria are shown in Fig. 1. A more technical summary of these processing and selection criteria is shown in the Supplement.

We determine the values for the first four statistical moments (mean, standard deviation, skew and (excess) kurtosis) of the resulting distributions. These in turn represent the central tendency, the spread, the asymmetry and the “peaked-ness” of the distribution. Distributions with non-zero skew or kurtosis indicate departures from a pure Gaussian shape. We note that the skew and kurtosis measures are not fully orthogonal (McKinnon et al., 2016) but will use these measures in this analysis because of their well-understood nature. Other statistical moments, e.g. linear (or L) moments (Hosking, 1990), have also been used for the analysis of changes in the characteristics of distributions (e.g. Fowler et al., 2005; Simolo et al., 2011). However, as noted above, we use the ordinary statistical moments to enable clearer comparison with previous similar analyses (e.g. Donat and Alexander, 2012).

Table 1 Linear change per decade in fits to parameters in zonal analysis for (a) daytime and (b) night-time temperature observations over 1974–2018. Values in bold show parameters and bands for which the 1σ range of the fitted trend does not include zero.

The linear changes in these ordinary moments, along with the $\mathrm{1}\widehat{\mathit{\sigma }}$ range, have been derived from the median of pairwise slopes method (e.g. Sen, 1968; Lanzante, 1996) and are shown for each latitude band in Tables 1, S1 and S2 for temperature, dew point temperature and wind speed, respectively. We do not expect any changes in the parameters of the distributions to be linear, but it is a useful way to summarise their gross changes over time together with the relative certainty of these changes. Latitude bands for which the range in trend values does not encompass zero are highlighted in bold, as for these we can be more confident that there is a true change in the value over time.

3.1 Temperatures

The results from daytime and night-time temperature observations for the Northern Hemisphere are shown in Fig. 2 and Table 1. We can be more confident of changes that are seen in latitude bands containing a large number of stations than in those which only contain a handful. The Southern Hemisphere is shown in Fig. S1 and Table S1, as it includes considerably fewer stations. It is clearly evident from these figures and tables that there has been a shift in the location of the distributions from the early period to the later period, with the most recent period being on average warmer than the earliest. This behaviour agrees with the results presented in Donat and Alexander (2012), who show increases in the mean of around 1 ^∘C between the early (1951–1980) and late (1981–2010) periods in their study (for all their regions and for both maximum and minimum temperatures). We note that their period of study is longer (60 years) compared to the period presented in this work (45 years). It is also unsurprisingly in agreement with the widespread warming of the Earth's surface (IPCC, 2013).

Figure 2 Dry-bulb temperature distributions in latitudinal bands for daytime (left) and night-time (right) observations for the Northern Hemisphere stations. The y axis has been normalised for ease of comparison. The number of stations contributing to each set of distributions is shown in the top left of each sub-panel, with the latitude range in the top right. The short vertical bars indicate the location of the 5th and 95th percentile values for each distribution.

In both daytime and night-time observations, positive increases in the mean that are clearly different from zero are seen in all latitude bands (Table 1). The strongest increases are seen in the highest northern latitudes, with increases on the whole becoming smaller when heading towards the Equator for both halves of the day. Changes in the Southern Hemisphere show less variation with latitude and are on the whole smaller than in the north, though still showing overall increases. There are few differences in the daytime and night-time increases for each latitude band, and in most cases the ranges overlap.

The northern middle to high latitudes also see a consistent decrease in the standard deviations in both daytime and night-time. However, south of 20^∘ N, the signal is less consistent. In fact, during daytime and night-time most tropical and southern latitude bands have ranges that encompass zero, though most bands suggest decreases. Although Donat and Alexander (2012) also found decreases in standard deviation in the extratropics, between 30^∘ S and 30^∘ N they found an increase when using zonal averages. This decrease suggests that the maximum temperatures are warming less quickly than the minimum temperatures, reducing the spread of the distribution in the Northern Hemisphere.

Consistent changes in the skew and kurtosis for both daytime and night-time or across latitude bands are not as apparent in Table 1. Slight increases in skew (growth in the positive tail and shrinking of the negative tail) are seen for the daytime high latitudes (north and south). In contrast, during night-time, some tropical regions show decreases in skew. The kurtosis decreases in the tropics at both times of day, with a more mixed signal at the high latitudes. Donat and Alexander (2012) also found an increase in the skew in all regions except for T_max in the tropics, which is broadly consistent with the changes shown in Table 1.

3.2 Dew point temperatures

The changes in the dew point temperatures (Table S1 and Figs. S2–S3) show a slightly more complex picture than for the dry-bulb temperatures. In general, the mean again shows an increase in both daytime and night-time values, with larger increases over time in the Northern Hemisphere high latitudes, similar to the dry-bulb temperatures. In the Southern Hemisphere, however, and especially at high latitudes, this increasing signal is less strong with few cases in which the confidence intervals do not span zero. One southern latitude band even shows a clear decrease in mean dew point temperatures over the last 45 years. Also, trends close to the Equator are greater than at latitudes immediately to the north and south of this equatorial belt.

The tendency has been for standard deviations in many northern latitudes to decrease, but there are also a number of increases, especially in the Southern Hemisphere, though these are determined from few stations. For the skew, in the daytime, this increases in the north and decreases in the south, with a similar but not as consistent pattern in the night-time. The changes in kurtosis do not show any clear pattern in comparison to the other moments.

3.3 Wind speeds

We also assess changes in the sub-daily wind speeds (for example, 1 or 2 min average values over the USA as in DeGaetano, 1997, or 10 min average values over the UK), but there are fewer stations where this variable has been sufficiently observed (see Table S2 and Fig. S4). In fact, all latitude bands in the Southern Hemisphere have fewer than 10 stations, and so we do not study this hemisphere. We use the same method as for the temperatures but note that the anomalies will be truncated at the negative end by the zero-bounded nature of this variable.

For the Northern Hemisphere, most latitude bands show reductions in the mean and also the standard deviation over time for both daytime and night-time observations, with an uncertainty range distinct from zero. In the lower mid-latitudes, there are indications for a small increase in standard deviations but at a lower magnitude than other changes. Tropical regions show decreases in the skewness over daytime and night-time, but high latitudes show an increase during the day.

These decreases in the mean wind speed over land have been noted in past studies (termed “stilling”; Roderick et al., 2007; McVicar et al., 2012), and monitoring studies have also shown decreases over time for the high and low wind speed values (e.g. Azorin-Molina et al., 2019). However, these recent regular global assessments also use HadISD, so differences are not expected. A number of possible explanations have been proposed to explain the stilling, which may be reversing in recent years (Azorin-Molina et al., 2019), including amongst others changes in the surface roughness because of land use changes (Vautard et al., 2010), variability of the atmospheric circulation (Azorin-Molina et al., 2016) and instrumental issues with wind sensors (Azorin-Molina et al., 2017, 2018).

3.4 Zonal summary

As can be seen in both Fig. 2 and Table 1 (as well as Fig. S1 in the Supplement) there is a clear increase in the mean temperatures with time across all latitude bands and times of day. However, this is also combined with a decrease in the standard deviation over the Northern Hemisphere. The width of the distribution of temperatures is therefore narrowing. Also, the skew of the distribution is becoming more positive, which indicates that the positive tail is growing and the negative tail is shrinking. Together, these three changes indicate that there has been a greater change in the low tail than in the high tail of the temperature distribution. The changes in the kurtosis are mainly increases (the tropical decreases rely on small samples), indicating a change to more peaked distributions and a reduction in both tails. Although Donat and Alexander (2012) found similar increases in the mean and skew of daily maxima and minima, they showed increases in the standard deviation for tropical regions not found here.

Markers indicating the location of the 5th and 95th percentile values in Fig. 2 track the changing tails of the distributions. Qualitatively for the Northern Hemisphere, the low tails have changed more than the high tails, as the vertical lines cover a greater spread, and more so for the higher latitudes than the tropical regions. These highlight changes in the standard deviation and also the skew on top of changes in the mean. However, there is no striking difference between the daytime and night-time plots for temperature. The number of stations in the Southern Hemisphere is much smaller, but there are no systematic differences between the behaviours of the low and high tails (see the Supplement). A similar pattern is seen for the dew point temperatures, with a greater shift for the lower tail and no striking difference between night-time and daytime values. The wind speeds show little change in the lower tails (Fig. S4) but decrease in the value of the 95th percentile over time. Further investigations into the changes in the 5th and 95th percentile values are presented in Sect. 5.

The intention of investigating the daytime and night-time values separately was to echo analyses conducted in other studies of daily maximum and minimum temperatures. Many of these show that the minimum temperatures are changing faster than the maxima, and so we expected to find a difference between daytime and night-time either for the parameters of the distributions or in the tail markers. However, the clearest difference remains between the upper and lower tails. This analysis shows increases in the mean and skew, as well as decreases in the Northern Hemisphere standard deviations. These changes in the statistical moments of the temperature distributions indicate that the lower tails are shifting more rapidly than the upper tails. Hence, as well as a change in the location of the distribution, the shape has also changed, resulting in changes to the occurrence of extremes. Although fewer stations contribute to the wind speed observations, these also show changes in location and shape, with larger changes in the upper tails caused by decreases in the mean and standard deviations.

Increases in the mean for the temperature and dew point temperatures are expected under a warming world (IPCC, 2013), and so it is the coherent changes in the higher moments which are of greater interest in this analysis. A greater change in the minima than the maxima can drive both the decrease in standard deviation and also the increases in the skew, but interestingly we see no difference between the magnitude of night-time changes compared to daytime ones. Some studies using climate models suggest that changes in extremes in a warming world (from anthropogenic forcings) come from changes in the shape as well as location of the temperature distributions (Clark et al., 2006; Hegerl et al., 2004). More recent work has looked at the effect of moisture on extreme temperatures in the soil (Whan et al., 2015) or the atmosphere (Sherwood and Huber, 2010; Matthews, 2018), where separation into latent and sensible heating when moisture is present can reduce the peak temperatures.

We note that observations in HadISD have a finite numerical resolution arising from the point of recording, through transcription, rounding, conversion and truncation, and in some cases combinations and multiples of these. As in this part of our analysis we blend anomalies calculated using daily climatologies from different stations together, and assuming a random distribution of these effects within a latitude band we do not discuss the impact of all these processes at this point, but we do so in the next section (Sect. 4).

This analysis combines station anomalies together in zonal bands separately for local daytime and night-time. In doing so, it removes any small-scale regional changes in the characteristics of the distributions of the meteorological variables. We also have not split up the analysis into seasons (neither 3-month nor wet–dry), and so large changes occurring in only part of the year will have been diluted in this analysis. In the following section we improve the assessment of regional and seasonal changes but in doing so reduce the number of observations available to characterise the distributions.

4 Station distribution changes

To investigate the geographical patterns of changes in these meteorological parameters in more detail than in Sect. 3, we study each station individually. This reduces the amount of data available to characterise the distribution at the gain of increased spatial information. A similar approach was used by e.g. Cavanaugh and Shen (2014) and McKinnon et al. (2016) in their investigations into the changing shape of temperature distributions on a daily basis from GHCND. However, the sub-daily nature of HadISD means we are able to investigate the distributions at different times throughout the day rather than changes in the maximum and minimum values independently. A large proportion of the stations in HadISD only report every 3 h rather than hourly. Therefore, we use the station longitude to adjust from UTC to local time, collating to the nearest 3-hourly time point (00:00, 03:00, etc.).

For each station a climatology over 1981–2010 is calculated for each day of the year for each 3-hourly time point ($\mathrm{365}\times (\mathrm{24}/\mathrm{3})$), requiring that at least 15 years of data are present. These 3-hourly climatologies are subtracted from the observations to create 3-hourly anomalies. The analysis is performed using observations over the entire year and also for standard 3-month seasons.

The 45 years (1974–2018 inclusive) are split into five 9-year periods, and the anomalies within each period and at each 3-hourly time point (00:00, 03:00, etc.) are combined annually and seasonally. Using fewer periods of longer length than in Sect. 3 allows for a better characterisation of the distributions at the individual station level at the expense of the number of periods. This is in contrast to the previous analysis (Sect. 3) wherein combining stations together increases the number of observations contributing to a distribution. For the 3-month seasons, using periods 9 years in length means there are maximally around 800 observations from which the distribution at a single time point can be characterised if all observations are present (3 months $\times \sim \mathrm{30}$ d ×9 years). The first four moments of the distributions are only calculated if there are more than 300 observations present within a 9-year period (roughly the equivalent of 1 month of observations at a given time point). A more technical summary of these processing and selection criteria is shown in the Supplement.

We again use the median of pairwise slopes method (Sen, 1968; Lanzante, 1996) to characterise the change in the moments over time, but only if there are sufficient data in at least four of the five periods. Again, we do not expect any changes to be linear, but this enables the convenient visual inspection of changes over time. Even if our data completeness criteria result in the inclusion of stations which have a nonuniform distribution of observations across the years and seasons, coherent, robust changes (if they are present) will still stand out across the stations.

4.1 Temperature

In Fig. 3 we show changes in the mean, standard deviation, skew and kurtosis for the temperature data for all seasons from HadISD at 15:00 and at 06:00 local time (LT) in Fig. 4 to correspond approximately with the times of maximum daytime and minimum night-time temperature, respectively. Other times are available in the Supplement. Despite 8139 stations being available in the HadISD v3.0.0.2018f dataset, after selection for homogeneity and length of record, only around 3200 have quantified changes in distributional parameters; the total number of stations is shown in parentheses at the bottom right-hand corner of each panel. Those stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of the fitted trend excludes zero are plotted with a larger symbol than those where this range encompasses zero. As is also clear from these figures, there are many fewer stations in the Southern Hemisphere than in the north. The highest concentration of stations is in Europe and eastern Asia, and this uneven meridional distribution will have affected the zonal average calculations in Sect. 3.

Figure 3 Trend over 1974–2018 in distribution parameters for temperature at 15:00 LT, showing the mean, standard deviation, skew and kurtosis. The stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of the fitted trend excludes zero are plotted with a larger symbol (the number of these is shown in the bottom right-hand corner of each plot). The total number of stations available is shown in parentheses.

Figure 4 Trend over 1974–2018 in distribution parameters for temperature at 06:00 LT, showing the mean (^∘C), standard deviation (^∘C), skew and kurtosis. The stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of the fitted trend excludes zero are plotted with a larger symbol (the number of these is shown in the bottom right-hand corner of each plot). The total number of stations available is shown in parentheses.

For the distributions of daytime maximum temperature (15:00; Fig. 3), it is immediately apparent that there are large and coherent changes in the mean and standard deviation over the period of record of these stations. Almost all changes in the mean are positive, with the strongest signal in Eurasia and lower values in North America. Approximately 95 % (3062 plotted out of a total of 3224) of the stations have trends wherein the uncertainty range does not encompass zero. The majority of changes in standard deviation at high northern latitudes are negative, along with most of central North America. There are increases over Europe and parts of China and southern Russia. Fewer stations (around 45 %, or 1542 out of 3224) have trends wherein the uncertainty range of the trends does not encompass zero.

Changes in skew and kurtosis have a more heterogenous signal, and again around half the stations exhibit changes wherein the uncertainty range does not include zero. However, the changes are not completely random, with consistent regional patterns in some areas, but these are smaller than for the standard deviation. The general pattern is for slight increases in the skew, except in a belt around the tropics and in south-east Europe. The kurtosis shows decreases in north–central Europe and a band extending eastwards from the White Sea across to China and India, with a more mixed signal elsewhere with smaller areas showing coherent changes.

The night-time observations (06:00) show very similar general patterns for changes in the mean and standard deviation, with quasi-global increases in the mean, decreases in the standard deviation at high latitudes, and increases for China, Europe and south-eastern Russia. However, the magnitude of the changes in the mean is lower than for the daytime (15:00) observation, especially over Europe.

There is a stronger change between the skewness for the 06:00 observations and those from 15:00, with larger regions showing a coherent signal. Decreases are visible more widely in Europe and also in India and parts of south-east Asia, and increases are more visible in the high latitudes and North America. Regions with a strong decrease in kurtosis are Midwestern North America, India through China and a band in north-western Russia, but in contrast to 15:00, large parts of Europe show an increase.

4.1.1 Temperature changes across the day

In Fig. 5 we show maps for changes in the mean for all eight of the 3-hourly time stamps. All hours of the day show an increase in the mean, with the strongest signals during the daytime (09:00, 12:00, 15:00). The region with the largest trends is Europe. However, it also has the highest station densities, and so the eye is more drawn to this region than areas which also show large trends but from a sparser station network.

Figure 5 Trend over 1974–2018 in distribution mean for temperature at each 3-hourly interval. The stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of the fitted trend excludes zero are plotted with a larger symbol (the number of these is shown in the bottom right-hand corner of each plot). The total number of stations available is shown in parentheses.

The trends in standard deviation show decreases at all hours in high northern latitudes, North America and southern Asia (Indian subcontinent and surroundings; Fig. 6). Increases in the standard deviation are seen in a band from Europe across to China. The skew from HadISD (Fig. 7) shows increases in most regions except the majority of Europe and parts of Asia. The decreases in skew over Europe are strongest in the mid-morning (09:00) and weakest in the late evening (21:00). Also, changes in standard deviation and skew are the inverse of each other when comparing Figs. 6 and 7.

Figure 6 Trend over 1974–2018 in distribution standard deviation for temperature at each 3-hourly interval. The stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of the fitted trend excludes zero are plotted with a larger symbol (the number of these is shown in the bottom right-hand corner of each plot). The total number of stations available is shown in parentheses.

Figure 7 Trend over 1974–2018 in distribution skewness for temperature at each 3-hourly interval. The stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of the fitted trend excludes zero are plotted with a larger symbol (the number of these is shown in the bottom right-hand corner of each plot). The total number of stations available is shown in parentheses.

4.1.2 Seasonal temperature changes

In Figs. S9–S12 we show the changes for the first four moments from the 15:00 observations across the four seasons of the year (and Figs. S13–S16 for the 06:00 observations), with the standard deviations at 15:00 being reproduced here in Fig. 8. Again, the strongest variations are over Eurasia. Changes in the mean temperature at 15:00 are smallest in boreal winter¹ (DJF) and strongest over a wide area in spring (MAM). However, summer (JJA) temperatures increase strongly with time around the Mediterranean and across continental Europe, and winter (DJF) has strong warming in the western part of Russia and Scandinavia. In North America, the largest trends are seen in SON (especially in higher latitudes), with little variation over the other seasons, especially spring and summer. A similar but less intense pattern is observed at 06:00.

Figure 8 Trend over 1974–2018 in distribution standard deviation for temperature at 15:00 LT across the four seasons. The stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of the fitted trend excludes zero are plotted with a larger symbol (the number of these is shown in the bottom right-hand corner of each plot). The total number of stations available is shown in parentheses.

Although increases in the standard deviation are seen across Asia in springtime and summer at 15:00, a more mixed pattern is present in this region in autumn; decreases are seen in the far east in winter (Fig. 8). At 06:00, the pattern is somewhat similar, though with the Indian sub-continent and surrounding regions showing consistent decreases throughout the year. Over North America, standard deviations increase during spring but mostly decrease during the other seasons, except in summer at 06:00. Further increases in the standard deviation are also found around the Mediterranean and in south-eastern Europe for most of the year, with decreases in the north of the continent (British Isles and Scandinavia). Increases in the standard deviation are generally weaker and decreases more prevalent across all seasons at 06:00 than at 15:00.

Patterns of changes in the skew have some clear features, with decreases in China for all seasons except winter. Europe also shows decreases in the spring but a more mixed pattern for the rest of the year. Increases in the skew at 15:00 are found in northern high latitudes during spring and India for most of the year. Over North America, decreases in skew are observed during spring and autumn in the southern half but in the northern half during summer.

4.1.3 Temperature discussion

The results at 15:00 and 06:00 can be compared to the maximum and minimum temperatures in the study by Donat and Alexander (2012) using HadGHCND between 1951 and 2010. Unsurprisingly, the increase in the mean temperatures is also visible in the HadGHCND study, along with numerous other studies on recent climate change. However, HadGHCND shows stronger increases in the minimum temperatures (analogous to 06:00 here), whereas HadISD shows stronger increases in daytime (15:00) temperatures. The HadGHCND dataset uses inverse distance weighting to convert the daily station observations into a gridded dataset. One aspect of this process is that it smooths the data and can also lead to very few stations contributing to large areas of interpolated values.

Donat and Alexander (2012) find larger regions showing stronger increases in standard deviation around Hudson Bay for both maximum and minimum temperature values than for HadISD. The heterogeneous nature of the changes in the standard deviation, along with the smaller fraction of stations showing clear trends, is in line with the study of Huntingford et al. (2013). For the skew, HadGHCND shows decreasing skew in Greenland, northern Africa, the coasts of eastern Asia and parts of Australia; it shows increases over Europe, in contrast to the results of this study, which in some regions may be the result of low station densities in HadISD (Fig. 1). The strong decreases in kurtosis observed close to the tropics in Sect. 3 are now shown to be concentrated in southern India and south-east Asia and not always linked to an opposite trend in the standard deviation.

Although McKinnon et al. (2016) link the first four statistical moments derived from GHCND between 1980 and 2015 to orthogonal basis functions, these still are closely tied to the original statistical moments. Over Europe, the geographical peaks of the trends in the daytime maximum mean are both in eastern Europe. However, HadISD does not show such a strong feature of increasing standard deviations across Asia, and GHCND does not show the decreasing daytime standard deviations in northern Europe, but those in eastern Asia do agree. The results from HadISD do not show daytime cooling in the eastern half of North America observed in GHCND. However, the general decrease in variance is observed in both. The increase in the central southern region of North America in GHCND is not present here, though the agreement is better by night.

Comparing the seasonal plots with those from Cavanaugh and Shen (2014, 2015) using daily GHCND temperatures over 1950–2010 shows generally good agreement for the changes in the mean. They find the strongest changes over Europe and a clear “warming hole” (Portmann et al., 2009) over the US, especially in JJA and DJF. In HadISD, the strongest changes are over Europe, and the location of the largest values moves between the seasons; but the warming hole is only weakly apparent during JJA over North America. For standard deviation, both Cavanaugh and Shen (2014) and our study find coherent decreases in the Northern Hemisphere and increases in MAM over North America. Over Europe, however, JJA shows the largest differences, with widespread increases in Cavanaugh and Shen (2014, 2015), but changes in both directions are clear in HadISD. McKinnon et al. (2016) also show the decrease in standard deviation over North America, except for a region of increase on the Gulf of Mexico. Also across Eurasia, the patterns are similar in HadISD and McKinnon et al. (2016). The positive trends in skewness found over North America (Cavanaugh and Shen, 2014, 2015) in GHCND are not as prominent in HadISD or in McKinnon et al. (2016) also using GHCND, with contiguous regions of negative trends in all seasons.

Some of these differences may result from inherent differences in the datasets, for example the differing temporal periods used, with HadISD using 1974–2018 compared to GHCND, which uses 1980–2015 and 1950–2010, and HadGHCND using 1951–2010. We note that the number of stations reporting in HadISD prior to 1973 is much lower, restricting the possibility of extending our analyses further into the past. Also, we and McKinnon et al. (2016) do not combine the stations together unlike in Sect. 3, Donat and Alexander (2012), and Cavanaugh and Shen (2014, 2015). The largest difference is that the sub-daily data within HadISD are not the true maximum and minimum values which comprise GHCND and HadGHCND. Although we have used the morning and afternoon values as proxies for the true daily maximum and minimum temperatures, this will not always be the correct match. This does, however, allow for the study of how extremes are changing throughout the day. It is difficult to say which of these analyses is the best at determining the change in the distributions of temperature given the different temporal periods and processing methods, station selections, and level of smoothing. Where a number of studies agree on the sign of the change, results are naturally more reliable than regions where they do not.

Changes in the temperature distributions for stations in China are likely to be the result not only of warming on a global scale, but also from more local effects arising from rapid urbanisation (Sun et al., 2016). Impacts from urbanisation are also likely in other regions undergoing rapid development in recent decades. As only a small fraction of the Chinese station network is unaffected by urbanisation (32 % in 2015; Wang et al., 2015), any changes are more representative of the urban areas (and hence a larger fraction of the population). Our study period coincides with the rapid urbanisation over China, and as noted above we see strong increases in the mean of the temperature distributions at all times of the day. Stronger warming is seen between 06:00 and 12:00 for northern and western parts of China but not for the southern coastal regions. Wang et al. (2017) found that the effect of urbanisation was greater on the minimum (morning) temperatures than the maximum (afternoon) ones, which we find in the northern and western parts but not in the southern coastal ones (Fig. 5). Although China is one example of recent and rapid urbanisation, with resulting effects on observations as well as people's experience of the weather, other parts of the world will also be affected. This may be to a lesser extent because of a station network better suited to climate monitoring or the urbanisation occurring over a longer time period.

A number of different causes for changes in temperatures beyond the anthropogenic warming signal have been suggested. These drivers of change are likely to influence the distributions of observations across hours of the day in different ways. Large-scale circulation changes would change the weather patterns experienced by different regions. Studies across the world have found local and global modes of variability linked to changes in cold and warm extreme temperatures (e.g. Scaife et al., 2008; Kenyon and Hegerl, 2008; Barrucand et al., 2008; Sillmann and Croci-Maspoli, 2009; Zongxing et al., 2012; Ning and Bradley, 2015). Hence, changes in these large-scale circulation patterns and modes will affect the extremes. However, taking the El Niño–Southern Oscillation as an example, there is no consensus on the change under global warming (Stevenson, 2012; Cai et al., 2014; Kim et al., 2014; Vega-Westhoff and Sriver, 2017). That said, as some of the changes in the higher moments shown here have a strong zonal component, especially over Eurasia, this suggests some influence of the poleward expansion of the Hadley cells (Hu and Fu, 2007), which have been linked to anthropogenic forcings (Lu et al., 2007; Lucas et al., 2014; Tao et al., 2016) on changes in the shape of the temperature distributions.

Other drivers of change in observed temperature distributions have been suggested, including enhanced land–atmosphere coupling and soil moisture (e.g. Seneviratne et al., 2006; Jaeger and Seneviratne, 2011; Hirschi et al., 2011) or changes in cloudiness and radiation (e.g. Lenderink et al., 2007). In a modelling study of future summer climate over Europe, Fischer and Schär (2008) find that the seasonal cycle as well as interannual and intraseasonal variability have separate effects on the temperature variability, and each of these are affected to a different level by the other drivers.

4.2 Dew point temperature changes

Comparing the changes in the dew point temperature moments across the day and the seasonal cycle (Figs. S17 to S28) with those from the temperatures shows that on the whole these two variables are quite closely aligned. Despite fewer stations passing the selection criteria for this variable, there are still coherent patterns of change. Over the eight time points during the day, the changes in the mean dew points show less variation than the temperatures, most clearly apparent over Europe in the early afternoon. This suggests that the relative humidity over this region has decreased for this part of the day. A similar but less striking decrease in daytime relative humidity is also seen over North America. This decrease in the relative humidity has been observed regionally and globally in the homogenised HadISDH dataset (Willett et al., 2014, 2019), which is based on HadISD.

The changes in the higher moments of the dew point temperature are also reasonably constant through the day but have a greater spatial heterogeneity and lower magnitude than the temperatures. Notably, there is a clear increase in skew over south-east Asia during the daytime and early evening, which is not seen in the temperature results, indicating a change in shape away from a low-dew-point tail towards a high-dew-point tail.

Seasonally, the greatest change in the dew point mean is in SON and DJF over Europe in both the daytime and night-time, contrasting with the temperatures of MAM and JJA for the day and also DJF at night. Again the standard deviations and kurtoses are very similar for day and night across all seasons between the dew points and temperatures. For the skews, the increase in skew over south-east Asia noted above is most prominent during SON and during the daytime, but it is also present at night and differs from the behaviour of the temperatures for all seasons.

4.3 Wind speed changes

There are fewer stations with sufficient observations to assess changes in wind speed (Figs. S29 to S33) and very few in the Southern Hemisphere. As a result we focus our analysis on North America and Eurasia. The proportion of stations where changes in the distribution could be assessed varies between each time step to a greater extent than for the two temperature variables. Despite this, regions with sufficient station densities show a diurnal cycle in the behaviour of wind speed trends. There are stronger trends for declining wind speeds during the daytime (09:00–18:00) than for the night-time hours. And there are also no large regions that show large increases (at the same magnitude as decreases) in the wind speeds at any time point. However, some individual stations do show some increases, for example in south-eastern Europe. For the other three moments (not shown) there are even smaller differences than for the mean of the annual changes over the day.

We show the seasonal changes at 15:00 for all four moments in Figs. S30 to S33. For the mean there are no large changes across the seasons. The changes in standard deviation are also generally unvarying compared to the temperature measures. There are increases in the USA, with stronger values during spring and autumn. All other regions with sufficient station density have decreases in the standard deviation. Both skew and kurtosis show no strong seasonal variations. Skews are increasing over North America, but there are no strong coherent changes elsewhere. The kurtosis is mostly decreasing in south-eastern Europe and also on the whole across Asia.

4.4 Discussion of station-level distributions

Across all seasons, times of day and locations, there are strong increases in the means of both temperature and dew point temperature with time as expected in a warming world (IPCC, 2013). Although the geographical distribution of stations is not uniform, there appear to be larger increases in Europe, especially during the spring and summer The largest increases in Europe also occur during the daytime rather than night-time, which matches the behaviour seen in Table 1 at 40–60^∘ N (the latitudes corresponding to Europe). The standard deviations show stronger decreases in the northern latitudes but more regions with increases in the middle and tropical latitudes in the Northern Hemisphere. Seasonally, North America shows decreases in the daytime and night-time standard deviation in winter and summer but increases in spring (Figs. 8, S10 and S14), which was also observed in Cavanaugh and Shen (2015). This indicates a change towards less variability in winter and summer but a change towards greater variability in spring for both daytime and night-time temperatures.

We also note that in Sect. 3, the standard deviations of the temperatures and dew point temperatures were decreasing more or less across the whole Northern Hemisphere, whereas in this study there are regions with clear increases. The zonal study combines the stations together, even if individually any change in the statistical moments over time would not show clear changes ($\mathrm{1}\widehat{\mathit{\sigma }}$ range of the trends encompassing zero). In contrast, Figs. 3 and 4 emphasise those stations where the $\mathrm{1}\widehat{\mathit{\sigma }}$ range of fitted trends does not include zero. As noted, this is only around half of the total number of stations, and many of the others show decreasing trends.

Changes in annual skew are different on either side of the North Atlantic, with predominantly increases in North America, decreases in Europe across the day and hence a rather small change in the zonally averaged results (Sect. 3). The seasonal results show increases at higher latitudes (Canada and Scandinavia) and decreases further south (USA and continental Europe) in spring. Other seasons show a more heterogeneous pattern in Europe but increases across central Asia in summer and autumn, similar to the patterns from GHCND in Cavanaugh and Shen (2015). In China, which was not included in their study, all seasons except winter show decreases in the skews. For the kurtoses, Cavanaugh and Shen (2015) found strong decreases in the southern parts of North America in the summer, which are not seen in HadISD. However, the strong decreases over parts of Europe in winter are present in HadISD.

The observations in HadISD have a finite numerical resolution, with more recent observations being more often recorded to 0.1 ^∘C and those from earlier parts of the record to 1 ^∘C. Underlying this is a more complex process of conversion between Fahrenheit and Celsius (or knots and miles per hour to metres per second), possible truncation and rounding (or vice versa), and human-induced errors during transcription, digitisation or recording. In the GHCND, Rhines et al. (2015) found that 65 % of all temperature observations were misaligned as a result of these effects, and so it is likely that a large fraction of the HadISD observations have also been affected.

Rounding has an effect on the mean and standard deviation dependent on the rounding resolution, any offset in the two scales from zero (rounding width and shift parameter), with zero bias for some combinations (Schneeweiß et al., 2010; Sheppard, 1897), and likely therefore on higher moments as well.

These changes are combined with the discrete measurement times, which also are more likely to be hourly in the more recent period, and more stations will be recording only every 3 or 6 h in the earlier period. More frequent recording results in a greater likelihood of values close to the true extremes being recorded. Hence, under a stationary climate, this could result in both an increase in the maximum temperatures and a decrease in the minimum temperatures (for a perfect sinusoidal diurnal cycle). Combined with a warming climate, this could result in an enhancement of the maximum temperatures recorded. There will also have been changes in instruments and recording practices over the period, though the homogeneity assessment should have identified any gross traces these have left in the data record.

It is possible that the effect of these data measurement issues is systematic across the entire global network. We also note that despite the automated quality control and homogeneity assessment of HadISD, there are likely to be data artefacts still present in this global dataset. Therefore, further investigation is needed to determine how these two aspects combine to affect the extremes recorded and their change over time in the sub-daily data. However, in our analysis we have looked for contiguous regions of coherent changes for large numbers of stations spanning geopolitical boundaries. Wholesale changes in national networks and their measurement routines causing systematic biases in the results presented here are likely to stand out as political areas, which are not observed.

5 Quantile regression

Studies of the characteristics of a distribution are useful in describing the overall changes in a meteorological variable. However, they do not make the changes in extremes immediately clear, although they can be calculated through e.g. the extreme value theory. An alternative way to study the changes in the climate that show changes in past extreme events is to investigate the changes in the quantiles of a distribution. This enables the question of how much warmer the warmest 10 % of days are now than they were in the past to be addressed. Quantile regression can be used to calculate trends in specified quantiles of a distribution rather than trends in the mean (Koenker and Bassett, 1978; Koenker and Hallock, 2001; Barbosa et al., 2011; Franzke, 2015). This allows, for example, cases in which there is little or no mean trend but an increased variance over time to be assessed more completely. In the case of climate observations, especially given the impact of extreme events on society and infrastructure, quantile regression allows the changes in these extremes over time to be made apparent.

The observations in HadISD naturally contain both an annual and a diurnal cycle. If run purely on these raw data, then the annual upper quantiles are likely dominated by summer daytime temperatures and the lower by winter night-time temperatures. Although using standard 3-month seasons would enable the more detailed study of winter and summer quantiles, those from the shoulder seasons (spring and autumn) would be dominated by the earlier or later parts of these seasons given the strong annual cycle. For this reason McKinnon et al. (2016) use daily maximum and minimum temperature data from GHCND over July and August only to reduce the effect of the seasonal cycle in their study.

Our approach is to use de-seasonalised data, which removes the seasonal aspect described above and also the issue that quantile regression results can be biased when applied to discrete data (Machado and Silva, 2005; McKinnon et al., 2016). The observations in HadISD are at a finite numerical resolution (e.g. temperatures are reported to the nearest 0.1, 0.5 or 1 ^∘C). We use a daily climatology to create hourly anomalies of the HadISD data, retaining the diurnal cycle. Within each calendar day, the numerical resolution remains unchanged by the subtraction of a daily mean from the hourly observations (within-day differences will still have 0.1, 0.5 or 1 ^∘C resolution). However, between each calendar day, the different climatology value will result in a different offset so that between-day differences would not result in discrete values. Hence, overall, the values passed to the quantile regression algorithm will be less discrete than the raw observations themselves. As discussed in Sects. 3.4 and 4.4, we note that the final resolution of HadISD hides the effects of conversion, rounding and truncation (multiple times in various orders) that could occur between measurements and the current use (Rhines et al., 2015).

Daily means were calculated for 1981–2010, requiring at least eight observations in each 24 h period (equivalent to 3-hourly data). The daily climatology was calculated, requiring at least 15 years of data to obtain an average. Finally, the climatology was smoothed using a five-point binomial filter and subtracted from the sub-daily data for each day to make the climate anomalies. An example of the changes in the percentiles of a time series is shown in Fig. 9.

Figure 9 An example of the quantile regression analysis for a single station, with the trends in the quantiles in red and the ordinary least squares (OLS) fit in yellow. This is for station 154810–99999, Constanţa Mihail Kogălniceanu International Airport, Romania.

Both Europe and North America have high station densities (see Sect. 4) compared with other regions, and so our analysis is focussed on these two regions (see also McKinnon et al., 2016). Furthermore, we show the results using a space-filling method, Voronoi tessellation (derived from Delaunay triangulation; Voronoï, 1908; Delaunay, 1934), which uses the perpendicular bisectors of the lines joining each station to split the land surface up into polygons, with each polygon containing a single station. These have then been coloured to show the results of the quantile regression. Stations for which the slope of the quantile shown is not significantly different from zero (p>0.05 as determined by the quantile regression algorithm) are shown in grey, but we do no further significance testing on these trends. It should be noted that areas of high station density will naturally have smaller polygons, but the human eye is drawn to the larger blocks of colour arising from sparse station networks. However, the advantage of this approach is that the underlying station information is retained as opposed to it being smoothed when applying some form of gridding method.

Figure 10 Trend over 1973–2017 in quantiles of temperature (^∘C per decade) at (a) 0.01, (b) 0.05, (c) 0.10, (d) 0.90, (e) 0.95 and (f) 0.99 over Europe during summer (JJA).

Figure 11 Trend over 1973–2017 in quantiles of temperature (^∘C per decade) at (a) 0.01, (b) 0.05, (c) 0.10, (d) 0.90, (e) 0.95 and (f) 0.99 over Europe during winter (DJF).

We show the changes in the quantiles (0.01, 0.05, 0.1, 0.9, 0.95 and 0.99) over Europe for both JJA and DJF in Figs. 10 and 11, respectively. In some subregions large numbers of stations show similar patterns, suggesting that there is a coherent change in the quantile in that area. However, within these there may be isolated stations and/or cells showing changes in the opposite direction. Although we have selected those stations that do not have many or large breaks in their time series (as determined by PHA; see Sect. 2), and HadISD has undergone a number of quality control checks, there will inevitably be a few data quality issues remaining.

The largest change in the summer (JJA) quantiles is observed in south-eastern Europe and western Asia (Fig. 10), focussed on the area north and west of the Black Sea. The increase observed in these regions becomes progressively less strong towards lower quantiles and also with distance from this region. This indicates that the temperature distribution is broadening over time, with the upper tails changing more rapidly than the lower tails; i.e. the warmest N % of events are becoming warmer as the percentile values change more rapidly than the cool events (where N is 1, 5, 10, etc.). In contrast, the regions close to the North Sea and especially Scandinavia show slight increases for the lower quantiles that become less strong for higher quantiles. For the north and west of the UK in general and Scandinavia, except western Norway, there are even decreases for the highest quantiles.

Conversely, in the winter (DJF), the strongest increases are in the northern regions and the lower quantiles, with Scandinavia, the Baltic states, Benelux, Denmark and northern Germany indicating a narrowing of the temperature distribution over time. In the higher quantiles, there are increases for eastern Europe but not as much as in the summer. Westernmost Europe shows very little change across all quantiles. We show the panels for the shoulder seasons (MAM and SON) in Figs. S34 and S35. There is a general increase for all the percentiles, and the cold quantiles in the north-east of the region shown are warming faster than the rest. But in contrast to Figs. 10 and 11 these increases are much slower and the regions are less coherent. Overall, these results show that warm extremes are becoming more common more quickly during the summer in south-eastern Europe than elsewhere across the continent. Similarly, cool extremes are becoming less common more quickly during winter in northern Europe than elsewhere.

Franzke (2015) performed quantile regression on daily mean, maximum and minimum temperatures between 1950 and 2013. The mean temperatures show strong regions of increase in both the 5th and 95th percentiles in eastern Europe and western Russia, especially north of the Black Sea, with the higher percentile region covering a larger area than for the results presented here. When using summer maxima, Franzke (2015) shows a more uniform increase over all of Europe at the higher percentile, except Scandinavia, which behaves similarly in Fig. 10. The magnitudes of the change in the winter temperature percentiles from HadISD are consistent with the studies of Barbosa et al. (2011) and Franzke (2015). The behaviour over Turkey stands out in Franzke (2015) but not in HadISD, suggesting a difference in the way these data have been processed or effects arising from the difference in temporal coverage or even station selection.

The difference in the behaviour of the lower and upper quantiles indicates a change in the variance of the sub-daily temperatures over Europe, which has already been suggested by the results from Sect. 4. During the summer, the variance in the large region to the north of the Black Sea has increased. In contrast, during winter the variance over Scandinavia and northern Europe has decreased. Both of these changes are also clearly visible in Fig. S10, showing the seasonal change in the standard deviation for each station. Franzke (2015) found a more widespread, uniform decrease in variance in the winter months over Europe but with a stronger decrease over Scandinavia. However, Franzke (2015) found no increase in variance during the summer in eastern Europe.

Over North America (Figs. S36 to S39), the trends in quantiles are on the whole smaller than for Europe for most of the year. In the spring (MAM), there are moderate increases in the south-west, which are stronger for the higher quantiles. For the summer (JJA) the higher quantiles show decreases in the northern and central regions but increases for the lower quantiles in the south-west. In the autumn (SON) and winter (DJF) the strongest changes are in the north of the continent at the lower quantiles, with more moderate increases in the south and the higher quantiles. Hence, the most rapid changes are in the cool extremes during autumn and winter in the north, with less rapid changes in the warm extremes in spring in the south-west. Using GHCND, Rhines et al. (2017) show a similar pattern, with a strong increase in lower percentiles for northern latitudes in winter. In the spring their data show decreases in the lower percentiles in the north and increases in the south, almost the opposite in the autumn, and general warming for the minimum temperatures in the summer. We also find the strongest changes for the lower percentiles during winter, indicating a reduction in the range of temperatures during that season.

We also show changes in the dew point temperatures (Figs. S40 to S43) and wind speeds (Fig. S44) over Europe. The dew point temperatures indicate small increases in the quantiles in most regions, with no strong regional differences for MAM and JJA for any quantile. However, in SON and DJF, there are stronger increases in the lower quantiles for the north of the region (Scandinavia and north-western Russia) but small increases for other regions and quantiles. So, although the summer warm extremes are becoming more common in south-eastern Europe, there is no corresponding change in the dew point temperatures, suggesting a decrease in the relative humidity. However, in the northern parts, the dew point temperatures are changing roughly in step with the dry-bulb temperatures.

Only the annual changes in wind speeds are shown (Fig. S44), with some blank tiles and regions because of the smaller number of stations with sufficient observations. There are no large changes in the lower quantiles but clear and larger decreases in the upper quantiles, which matches the behaviour observed in Sect. 3 and the decrease in variance observed in Sect. 4.3. The larger changes in the upper quantiles also suggest the skews are becoming more negative, with decreases in the number and intensities of extremely high wind speeds.

6 Summary

To investigate the changes in the shape of the distributions of temperatures and wind speeds in the HadISD dataset, we have performed three types of assessment. Firstly, by combining the data from the station-based observations into zonal averages, the resulting distributions are built from large numbers of observations at the expense of some spatial information. Using 5-year bins from 1973 to 2018, we demonstrated the change in the distributions and percentile values. As expected, the mean temperatures increase for both daytime and night-time (IPCC, 2013). However, the standard deviations decrease for the Northern Hemisphere (where there are the most stations). Although changes are seen in both tails of the distribution, they are greater in the lower tail. This behaviour is consistent with both the more rapid change in low temperature extremes (Davy et al., 2017) and the decrease in diurnal temperature range (Alexander et al., 2006; Thorne et al., 2016a, b). The dew point temperatures show similar, if more heterogeneous, behaviour in this analysis. Atmospheric moisture, as well as that in the soil, may be acting as a brake on the change in the maximum temperatures, resulting in the consistent changes in the shape of the distributions observed.

In contrast, the wind speeds showed decreases in the mean but also in the standard deviations, with larger shifts in the upper tail than the lower. The changes in distributions observed on a zonal level suggest this may be driven by changes in the upper tail. A number of causes for this “stilling” have been proposed, from changes in land use and circulation patterns to instrumental degradation.

To improve the spatial discrimination, we assessed changes in the distributions for each station individually. The sub-daily nature of HadISD allows these to be calculated across the hours of the local day as well as annually or seasonally. The temperature means show the expected increase over time, but with a daily cycle showing a more rapid increase in mean temperatures during the early afternoon than in the early hours of the morning. However, the higher-order moments are much less homogeneous but have a smaller variation throughout the day. Trends in the standard deviations are on the whole negative, except for a band from the Mediterranean across to central Asia. Conversely, trends in the skew of the temperatures are mainly positive, except over Europe and across into western Asia, indicating increasing warm extremes (positive tail) and decreasing cool extremes. These patterns are in some regions different from those reported in previous investigations by Donat and Alexander (2012) and McKinnon et al. (2016) obtained using (derivatives of) GHCND (true maximum and minimum values). For the seasonal changes, which can be compared to Cavanaugh and Shen (2014), there are also differences to the patterns obtained from GHCND.

These disparities may result from a variety of differences in the datasets and methods. Firstly, the period used in this study is 1974–2018 compared to 1950–2010 (Donat and Alexander, 2012; Cavanaugh and Shen, 2014) and 1980–2015 (McKinnon et al., 2016). Also, by using all of the sub-daily observations in HadISD, this work incorporates the effect of changes across the entire distribution, even if we are concentrating on either the moments or specific percentile values (e.g. all temperature observations rather than just the maxima and minima). The studies using GHCND use the true maximum and minimum values, which allow for the focussed investigation of the momentary extremes. We also do not perform any spatial smoothing of the underlying stations, which could result in less variable fields with smaller extremes.

A number of different drivers for changes in temperatures have been proposed, including changes in global and large-scale circulation patterns, land–atmosphere coupling, and cloudiness and radiation. And a number of these have also been shown to be changing as a result of anthropogenic climate change. Over Eurasia, the quasi-zonal patterns of change of the higher moments suggest a link to the expanding Hadley cells, but further study is required.

Finally, we perform a quantile regression analysis to determine how the values of the percentiles have changed and therefore how the warmest days and coolest nights feel in comparison to others. This approach is in contrast to the ETCCDI, which use percentiles (e.g. TX90p) that capture the exceedance rate of a fixed percentile value. Here we show that over Europe, the upper percentiles are changing fastest in the large region north-west of the Black Sea in the summer, but the lower percentiles are changing faster in northern Europe and Scandinavia in the winter. Hence, the warmest N % of summer days are becoming warmer more rapidly than the coolest N % in south-eastern Europe, and the reverse is true for winter days and particularly nights in the north. Coherent changes in dew point temperatures are less clear, with increases in the lower percentiles in autumn and winter for the northern parts of Europe. For the wind speeds, as expected from the other analyses, the higher quantiles show a decline over Europe, with stronger decreases in the east.

All the studies presented herein show widespread increases in mean temperature and less rapid rises in the dew points, as expected from numerous other studies on the observed change in global temperatures. However, the changes in the higher moments are of more interest. Although in many cases (seasons and quantities) the patterns of change are very similar to prior studies (e.g. Donat and Alexander, 2012; Cavanaugh and Shen, 2014, 2015; McKinnon et al., 2016), there are still differences in the spatial patterns between all of the studies. These could arise from differences in the input data (daily extrema versus hourly values) as well as the temporal and spatial coverage. Furthermore, we have used a single station-based dataset, and differences between datasets (Gross et al., 2018) or the effect of gridding (Cavanaugh and Shen, 2015) can change the apparent behaviour of higher moments. Nevertheless, in these investigations, we have shown that there are changes in the mean, the standard deviations and the skews of the distributions studied. And, when using a sub-daily dataset of temperatures, changes in the number and intensity of extremes appear to arise from the combination of changes in all three of these moments. Further investigations are, however, needed to determine the underlying drivers for the changes in these distributions and hence the changes in extremes.