Observed Changes in One-in-20 Year Extremes of Canadian Surface Air Temperatures

Weather and climate extremes are often associated with substantial adverse impacts on society and the environment. Assessment of changes in extremes is of great and broad interest. This study first homogenizes daily minimum and maximum surface air temperatures recorded at 146 stations in Canada. In order to assess changes in one-in-20 year extremes (i.e., extremes with a 20-year return period) in temperature, annual maxima and minima of both daily minimum temperatures and daily maximum temperatures are derived from the homogenized daily temperature series and analyzed with a recently developed extreme value analysis approach based on a tree of generalized extreme value distributions (including stationary and non-stationary cases). The procedure is applied to estimate the changes over the period 1911 to 2010 at 115 stations, located mainly in southern Canada, and over the period 1961 to 2010 at 146 stations across Canada (including 37 stations in the North). The results show that warming is strongest for extreme low temperature and weakest for extreme high temperature and is much stronger in the Canadian Arctic than in southern Canada. Warming is stronger in winter than in summer and stronger during nighttime than daytime of the same season.


Introduction
Weather and climate extremes can have far-reaching consequences for our society, economy, and environment. Assessment of changes in climate extremes has become an important topic in current climate research. Thus, there have been a large number of studies on extremes, as summarized and reported by the Intergovernmental Panel on Climate Change (IPCC, 2012).
There are two fundamental approaches to the analysis of extremes. One is based on extreme indices, such as analyzing time series of the annual count of days above or below a percentile-based threshold (e.g., Alexander et al., 2006;Vincent & Mekis, 2006;Zhang et al., 2011), which represent moderately extreme events with re-occurrence times (i.e., return periods) of one year or shorter. Another is based on fitting an asymptotic extreme value distribution, such as the generalized extreme value (GEV) distribution, which provides insight into the behaviour of extremes with return periods of several years to decades that are important for engineering design and risk management (e.g., Kharin, Zwiers, Zhang, & Hegerl, 2007;Wang, Trewin, Feng, & Jones, 2013;Zwiers, Zhang, & Feng, 2011).
There have been a few studies on extreme indices of temperature in Canada (e.g., Bonsal, Zhang, Vincent, & Hogg, 2001;Vincent & Mekis, 2006), generally finding that there are fewer cold nights, cold days, frost days, and more warm nights, warm days, and summer days across Canada (Vincent & Mekis, 2006). An increase of 1.5°C in the national average annual mean temperature over the 1950-2010 period has been reported (Vincent et al., 2012).
However, there has not been a study of Canada-wide temperature extremes with return periods of several years to decades nor any assessment of changes in temperature extremes as inferred from extreme value distribution analysis. This study is the first attempt in this regard. Annual extremes of daily minimum (T min ) and daily maximum (T max ) temperatures recorded at 146 stations in Canada are analyzed using a GEV tree approach to infer changes in one-in-20 year extremes (i.e., 20-year return values), as detailed in Section 2. As in Wang et al. (2013), here we choose to focus on one-in-20 year extremes. This choice is a compromise between the rarity of the event of interest and the uncertainty in estimating the return values. Here, the objective is to describe systematic changes that have occurred in extremes over time without determining the underlying physical causes of those changes. Thus, changes could be a result of both external forcing and persistent long-term natural variations.

Data and methodology
The stations analyzed include 115 stations that have daily temperature observations for at least 90 years during the 1910-2011 period, as well as 31 northern stations that have at least 44 years of data during the 1951-2011 period. The 31 stations with shorter records are included to represent northern Canada (north of 56°N) because only six of the 115 long-term stations are located in northern Canada and none of these six stations is located north of the Arctic circle or in Nunavut.
We first homogenize the T min and T max daily data time series for changepoints (points in time when a sudden change occurs) identified in Vincent et al. (2012) by applying the RHtestsV3 software  with reference series to the corresponding monthly mean data series and using available metadata. The quantile-matching (QM) adjustment method (Wang, Chen, Wu, Feng, & Pu, 2010) is used to homogenize the daily data series. The QM method aims to diminish the difference between the empirical distributions of the daily temperatures before and after a changepoint. Here, the QM adjustments for daily data are estimated using the same reference series that was selected and used in Vincent et al. (2012). The reference series is the best correlated neighbour series that has a long period of homogeneous data surrounding the changepoint to be adjusted. As in Vincent et al. (2012), up to 10 years of data before or after a changepoint are used to estimate the QM adjustments with the cumulative distribution being evaluated for 12 quantile categories whenever there is enough data to do so . Figures 1a and 1b show the QM adjustments that were used to homogenize the Mont Joli daily temperature series. These were estimated using the Causapscal daily temperature series for the 1933-52 period as a reference. The Mont Joli temperature series has a shift at the end of 1942 because of the joining of records from two sites. Before 1943, the observations were taken at Pointe-au-Père, which is located near the St. Lawrence river. From 1943 onwards, the observations were taken at Mont Joli, an inland site near the airport. For both T max and T min , as shown in Figs 1a and 1b, the QM adjustments show a clear seasonal cycle, with negative adjustments for below-average temperatures (lower quantiles) and positive adjustments for above-average temperatures (upper quantiles). These adjustments have noticeable effects on the estimate of trends for extreme high and extreme low temperatures, as shown in Figs 1c and 1d. Note that the inland site could be much warmer than the site near the river in summer (especially in summer daytime), and the site near the river is slightly warmer in winter nighttime. Thus, the site change has a much greater effect on extreme high temperatures than on extreme low temperatures (Figs 1c and 1d).
For each station, we derive the annual minimum and maximum of T min , denoted as T Nn and T Nx , respectively, for each year that has at least 310 valid daily T min (i.e., at least 85% of daily observations in the year are not missing). Similarly, we derive the annual minimum and maximum of T max , denoted as T Xn and T Xx , respectively. A few stations that satisfy the above requirements but have too many missing values in the same season for several years or have missing values for a sub-period of a few years are also excluded.
Return values can be estimated by fitting a GEV distribution to annual maxima or by fitting a generalized Pareto distribution (GPD) to excesses of peaks-over-threshold (Coles, 2001). In this study, we take the GEV approach because it has been shown to be comparable to the GPD approach in terms of stability of the estimates when the sample size is more than twice the return period of the extreme to be estimated (Wang et al., 2013), which is the case here. We have 90-102 years of data in the 1910-2011 period (or 44-61 years in the 1951-2011 period) for estimating 20-year return values (i.e., one-in-20 year extremes; the expression "onein-20 year" is often omitted hereafter). We apply a GEV tree analysis to data for the period from 1910 or earlier to 2011 at the 115 long-term stations and also to data for the period from 1951 to 2011 at all 146 stations.
Here, the GEV tree analysis procedure is similar to that proposed in Wang et al. (2013), except that the shape parameter is kept constant for all stations. This is because we find that, for Canadian temperature extremes, the shape parameter is often constant and the occasional inclusion of any trend in the shape parameter could result in a sub-optimal fit because Observed Changes in One-in-20 Year Extreme Temperatures in Canada / 223 return value estimates are extremely sensitive to shape estimation error. Thus, the GEV tree has a stationary GEV model at its root and different non-stationary GEV models at different upper levels. The model with a non-linear trend in both the location and scale parameters is at the top of the tree. The most suitable model for each time series of annual extremes is chosen from the GEV tree using likelihood ratio tests. When the most suitable model is non-stationary, it indicates that the time series in question has experienced changes that are significant (increase or decrease) at the 5% level or higher; when it is stationary, it indicates that the time series in question has not experienced a significant change. Details of the GEV tree analysis approach are described in Wang et al. (2013).
The non-stationary models (i.e., models with one or more parameters that vary over time) can have linear or non-linear temporal trends. Here, the non-linear trend is of the form a + bt + ct k , where k [ {1.1, 1.2, …, 2.0} (this range is determined by our fitting experiments with the data). We set k = 1.5 at the beginning of the GEV tree analysis, and we refine the k of T Xx (annual maxima of T max ) and T Nn (annual minima of T min ) before and after adjustments (red and black curves), and their 20-year return values before and after adjustments (red and blue lines) as estimated using the GEV tree approach for the Mont Joli station. The grey shading (the interval between the red dashed lines) represents the 95% confidence interval of the blue (solid red) lines.
value (namely, choose the k value from {1.1, 1.2, …, 2.0} that provides the best fit) when the most suitable model involves a non-linear trend (i.e., when the fit with a polynomial trend of order k = 1.5 is significantly better than the fit with a linear trend). For each parameter, we first test if it is necessary to include a linear trend before testing if it is better to include a non-linear trend with k = 1.5. A trend in the location parameter is introduced when it is significant at the 5% level, but a trend in the log of the scale parameter is introduced only if it is significant at the 0.1% level.
In summary, the model selection procedure allows a trend that is more likely to be included in the location than in the scale parameter, while the shape is kept constant. This has been shown to provide much better fits than allowing a trend to be equally likely in any of the three parameters (Wang et al., 2013). The resulting k values for T Nn and T Xx are reported in supplementary Tables S1 and S2 (see online supplementary material at http://dx.doi.org/10.1080/07055900.2013.818526). For T Nn , T Nx , T Xn , and T Xx , a polynomial trend is identified at 5, 17, 5, and 13 out of the 115 stations over the century-long period, respectively (see supplementary Table S1 for T Nn and T Xx ), or at 11, 7, 8, and 6 out of the 146 stations over the 1951-2011 period (see supplementary Table S2 for T Nn and T Xx ). The resulting most suitable model is also listed in supplementary Tables S1 and S2. In general, changes in the location parameter dominate the change in Canadian temperature extremes. For T Nn and T Xx at all 146 stations, changes are identified only in the location parameter, with no significant trends in the scale and/or shape parameters (supplementary Tables S1 and S2).
Finally, we use the fitted most suitable GEV model to estimate time series of 20-year return values of the respective extreme variable, which are shown as the red curves for a few selected stations in Fig. 2. Also, we use the asymptotic variance-covariance of the parameter estimates (Coles, 2001) to estimate the 95% confidence interval of the one-in-20 year extreme value (Fig. 2, grey shading).
A 20-yr return value (RV20y) is estimated for each year in the period of the data that were used in the GEV model fitting. For the 115 stations, the RV20y time series starts in 1910 or earlier and ends in 2011, and these are used to estimate changes in RV20y over the period from 1911 to 2010. In order to estimate changes over the period from 1961 to 2010, the RV20y time series are also derived for all 146 stations for the period from 1951 to 2011, using the most suitable model chosen for the data series over the period 1951 to 2011 (supplementary Table S2). More specifically, the RV20y for the climate condition of year 2010 (i.e., the time-dependent GEV parameters take their values at the time that corresponds to year 2010, shown in Fig. 3 for T Nn and T Xx ) is estimated, as is the RV20y for the climate condition of year 1910. The difference between these two RV20y values (RV20y 2010 minus RV20y 1910) represents the change in RV20y over the period 1911 to 2010, which is shown in Fig. 4 (also reported in supplementary Table S1 for T Nn and T Xx ). Similarly, changes in RV20y over the period 1961 to 2011 are also estimated and shown in Fig. 5 (also reported in supplementary Table S2 for T Nn and T Xx ).
As shown in Fig. 2 (red curves), the 20-year return value changes over time when the most suitable GEV model is non-stationary (i.e., when there is a temporal trend in the location and/or scale parameters). For each time series of 20-year return values, we first calculate the 10 yr means for each of the decades 1911-20, 1921-30, …, and 2001-10, and then derive the differences between each of these 10 yr means and the 10 yr mean for the reference decade 1911-20 or 1951-60. Such differences, which represent changes in 10 yr mean RV20y since the reference decade, are then arithmetically averaged over all the stations analyzed or over the southern and northern stations, separately, as shown in Fig. 6. An estimate of the rate of regional mean change in the 10 yr means is estimated by a linear fit to the mean changes in each of the four extreme variables and also reported in the legend of Fig. 6. The linear fit seems reasonable for most cases, except for the extreme T Nn and T Xn , and the T mean for the North, which appear to have a polynomial trend with an accelerating rate of warming over the last few decades (Figs 6c to 6e; changes in the earlier decades lie below the linear trend line, and changes in the last decade lie above the trend line).

Observed changes
Over the past century, warming has been much stronger and more extensive for the extreme low temperature than for the extreme high temperature, as shown in Fig. 4. Warming is identified at 88 (76%) of the 115 long-term stations for the extreme low temperature (Fig. 4a), and at 38 (33%) of the 115 stations for the extreme high temperature (Fig. 4b).
Warming is also identified at 58 stations (50%) for the extreme low daily maximum temperature (Fig. 4c) and at 73 stations (63%) for the extreme high daily minimum temperature (Fig. 4d). The regional mean warming rates (arithmetically averaged over the 115 stations) are estimated to be 3.5°C, 1.9°C, 1.2°C, and 0.5°C per century for the extreme low daily minimum (T Nn ), extreme low daily maximum (T Xn ), extreme high daily minimum (T Nx ), and extreme high daily maximum (T Xx ), respectively (Fig. 6a). The rates are slightly lower when averaged over the 109 southern stations (Fig. 6b). Because there are few stations with century-long data records in the Canadian Arctic, the above regional mean warming rates mainly represent the warming in southern Canada.
Warming is found to be much stronger and more extensive in the Canadian Arctic over the 1961-2010 period than in southern Canada. As shown in Fig. 5a, warming is identified at 71 (48%) of the 146 stations and at 31 (83%) of the 37 northern stations, for the extreme low temperature (extreme T Nn ). The results also show more significant trends in western Canada than in eastern Canada (Figs 5a and 5c). For the extreme high temperature (extreme T Xx ), warming is seen mainly in the Canadian Arctic and in the eastern maritime area, with 14 (37%) of the 37 northern stations, or 28 (19%) of Observed Changes in One-in-20 Year Extreme Temperatures in Canada / 225 Fig. 2 Time series of T Xx , T Nx , T Xn , and T Nn (black curves) and their 20-year return values (red curves) estimated using the GEV tree approach for the indicated stations. The horizontal axis is year. The grey shading represents the 95% confidence interval of the red curve.
Warming of the extreme high nighttime temperature is also more extensive than for the extreme high daytime temperature (Figs 5b and 5d).
The trend significance test is performed separately for each time series at each station. In this case, one needs to evaluate the joint statistical significance of the multiple tests (i.e., field significance). This is because, statistically, 5% of the stations Fig. 3 The 20-year return value (RV20yr ;°C) of annual minima T Nn and annual maxima T Xx , respectively, estimated for the 2010 climate using the GEV tree approach (see Section 2).
Observed Changes in One-in-20 Year Extreme Temperatures in Canada / 227 analyzed can be identified by chance to have a significant trend (or no trend) when these stations are independent of each other, and this percentage could be a little higher when these stations are not independent of each other. In this study, field significance is assessed by performing Walker's test because Wilks (2006) has shown that this test is more powerful than the traditional method of Livezey and Chen (1983) and relatively insensitive to non-independence of the local test results. The Walker's test results suggest that the changes in each of the four extreme variables are field significant at the 5% level, in the field of 115 or 146 stations. Nevertheless, this does not eliminate the small possibility that two neighbouring stations are identified by chance to have different trends. As shown in Fig. 6, warming is strongest for extreme low temperature and weakest for extreme high temperature; it is generally much stronger for winter extremes (i.e., T Nn and T Xn ) than summer extremes (i.e., T Nx and T Xx ). For the onein-20 year extremes of T Nn , T Xn , T Nx , and T Xx , the regional mean rate of warming over the 1951-2010 period is estimated to be 3.0°C, 2.1°C, 1.4°C, and 0.5°C per century, respectively, for southern Canada (Fig. 6c); it is as high as 6.8°C, 6.2°C, 1.7°C, and 1.7°C per century, respectively, for northern Canada (Fig. 6d), and it is 3.9°C, 3.1°C, 1.5°C, and 0.8°C per century, respectively, for the national average (Fig. 6e). For the extreme T Nn over southern Canada, the regional mean warming rate is estimated to be 3.4°C and 3.0°C per century over the last century and over the last six decades, respectively (Figs 6b and 6c), indicating a slight slowdown in warming of extreme low minima in the last half of the twentieth century. However, these are linear estimates; the polynomial fit to the mean changes in the extreme T Nn over the last six decades is better than the linear fit (Fig. 6c). The polynomial trend suggests an acceleration of warming for extreme low temperatures in recent decades. The polynomial order k (for t = 1, 2, …, 6) is 1.2, 1.4, and 1.3 for the solid red curves in Figs 6c, 6d, and 6e, respectively, and 1.5 for the dashed red curve in Fig. 6d. This suggests that the rate of warming has been accelerating more rapidly in the North (k = 1.4 for T Nn ) than in the South (k = 1.2 for T Nn ; Figs 6c and 6d). Within the same season, warming is stronger in the extremes of daily minimum temperatures than in the extremes of daily maximum temperatures (i.e., the extreme T Nn has warmed up more than the extreme T Xn , and the extreme T Nx , more than the extreme T Xx ). This implies that warming is stronger during nighttime than daytime of the same season because daily minimum temperatures usually occur during nighttime and daily maximum temperatures during daytime.
For comparison, we also derive the 10 yr means of annual mean temperature (T mean ) and changes therein in the same way (i.e., first find the most suitable polynomial trend for each annual mean series, then derive the 10 yr means from the fitted trend line). Our estimates of annual mean temperature warming rate over the 1951-2010 period are 1.5°C and 3.3°C per century for southern and northern Canada, respectively, and 1.9°C per century when averaged over the 146 stations across Canada (Figs 6c to 6e). These are similar to previously reported values, although different numbers of stations were analyzed with different methods used to estimate the changes. Using 338 stations in Canada, Vincent et al. (2012) reported that the annual mean temperature had increased by 1.5°C over the 60-year period from 1950 to 2010, which is equivalent to 2.5°C per century and is a little higher than our estimate of 1.9°C per century for the country-wide average. They also reported that the annual mean of T max and T min has increased by 1.4°C and 1.7°C, respectively, over the 1950-2010 period. That is, the annual mean of T min has warmed up more than that of T max . As shown in Fig. 6, the annual mean temperatures have warmed up more than the summer extremes but less than the winter extremes.

Summary and conclusions
We have homogenized the time series of daily minimum and daily maximum surface air temperatures recorded at 146 stations across Canada and derived the annual extremes. We have chosen the most suitable GEV distribution for each time series of annual extremes from a GEV tree that consists of both stationary and non-stationary distributions. The most suitable fitted GEV is then used to estimate historical changes in one-in-20 year temperature extremes, including the extreme low daily minimum (T Nn ), extreme high daily minimum (T Nx ), extreme low daily maximum (T Xn ), extreme high daily maximum (T Xx ). Our results show that warming is strongest for the extreme low minimum temperature and weakest for the extreme high maximum temperature and is much stronger in the Canadian Arctic than in southern Canada. Also, warming is stronger in winter (i.e., for extreme low temperatures, T Nn and T Xn ) than in summer (i.e., for extreme high temperatures, T Nx and T Xx ), and stronger during nighttime than during daytime of the same season (T Nn versus T Xn or T Nx versus T Xx ). Over the 1951-2010 period, the rates of warming are estimated to be 3.9°C, 3.1°C, 1.5°C, and 0.8°C per century for the national average of one-in-20 year extreme T Nn , T Xn , T Nx , and T Xx , respectively. For the average over northern Canada (north of 56°N), the warming rates are as high as 6.8°C, 6.2°C, 1.7°C, and 1.7°C per century, respectively. The warming rate for the annual mean temperature is lower than that for the extreme low temperature but higher than that for the extreme high temperature. On a national average, it is 1.9°C per century over the 1951-2010 period.
Overall, our results are consistent with those reported in previous studies, particularly in the sense that Canada has become much less cold but not much hotter.