Journal of Earth Science & Climatic Change
Open Access

Climate trends; Index; Temperature; Precipitation; Great lakes region; Trend detection; Extreme climate events

Introduction

As global warming continues the upward trends in annual mean temperature across North America since 1960 are widespread but spatially non-uniform. Pronounced polar amplification of warming is observed in high latitudes [1]. Extreme high temperature records across North America (NA) are being set more frequently than extreme cold records and the probability of cold extreme events is reduced [2]. Trends in daily maximum and minimum temperature are significant in high latitudes [3-5]. Annual precipitation has increased in recent decades over the eastern NA, and decreased over the western NA, with considerable spatial variability. The intensity and frequency of one-day heavy precipitation events have very likely increased since the mid-20th Century across most of the US but no detectable trend is reported in Canada [6]. In Canada during the 20th century, the annual mean temperature has increased between 0.5 and 1.5°C in the south. The warming is greater in minimum temperature than in maximum temperature in the first half of the century, resulting in a decrease of diurnal temperature range (DTR). Annual precipitation has also increased from 5% to 35% in southern Canada over the same period (Zhang et al. 2000). Spatially, for the precipitation indices, no significant trends are observed for annual total precipitation and extremely wet days in the southwest and the central Ontario [7-8]. For temperature indices, cool days and warm nights have significant trends in more than 90% of southern Ontario [9]. The trend is projected to continue in this century [10-11]. In general, the spatial variability of precipitation indices is much higher than that of temperature indices [12-13].

Climate change has caused widespread adverse impacts and related losses and damages to nature and people, beyond natural climate variability. Accelerating climate change hazards pose significant risks to the wellbeing of populations, food production, freshwater resources, supply security and the natural, managed and human systems on which they depend. In addition, these climate extremes have resulted in many casualties and caused enormous damage to property, infrastructure, and human health [14-15]. Such as scorching heat waves in Canada and the United States of America (US and Canada heat wave 2021).

The Great Lakes region (GLR), an engine of growth in North America, is at risk from climate change. Furthermore, any adverse impacts on the GLR affect its immediate surrounding areas as well. An investigation of climate change at regional scale is very important to identify the impacts of climate change on all aspects of society. Although there are many studies about climate change in the GLR based on historical observation and future projection datasets; most of these studies used station data or low-resolution grid data. Extreme weather/climate events often happen at a local and short time scales. High spatial-temporal resolution data is critical for such studies. Accurate and timely observation will allow a deeper understanding of climate change and associated impacts. The recently released ERA5 data have a high spatial and temporal resolution, which makes them suitable for further investigation of climate extremes in these regions [16]. The climate change research community has defined a number of indices as indicators of these extreme events [17-18]. These indices have been used in many studies for the detection and attribution of extreme events [19-22]. The applicability of these indices depends on the geographic location and specific purpose of the research interest. Using the recently released high resolution ERA5 data, this study will focus on trend detection of 38 climate indices over the GLR that are widely used in Canada and the USA.

Many parametric and non-parametric statistical tests have been developed to assess trends in climate change. Such as linear regression [23-24], Mann–Kendall (MK) trend test with Sen Slope estimation [25], Spearman’s rho [26-27], innovative trend analysis [28], partial trend analysis [29] and crossing trend analysis. Each trend identification technique has a set of restrictive assumptions and limitations. Among these assumptions, the most significant ones are the normal (Gaussian) probability distribution function (PDF) and serially independent structure of a given time series. [30] Reviewed these methodologies (2020). Furthermore, for monotonic trend detection the most widely used technique is the non-parametric MK test [31], which has no normal distribution restriction. The MK trend test is frequently used for significant trend identification by numerous researchers [32-38]. The classical Mann–Kendall (MK) trend test has been employed frequently in climate trend detection [39] with Sen Trend slope calculation.

In this paper, we will examine trends in climate change indices in the GLR and surrounding area over the past 40 years from January 1981 to December 2020. We will not include data for the period 1950- 1980 due to the difference in data production process [40]. Because most of temperature and precipitation related indices generally do not satisfy the Gaussian assumption, their trends will be detected with the MK method; however, the MK has an assumption that the time series should be serial independent. To deal with this issue, it is essential to convert available time series to serially independent versions to get reliable results. For this purpose, pre-whitening [41] and more effective over-whitening (Şen 2017) procedures were suggested. To investigate the impacts of the presence of serial correlation in a time series and to deal with the associated uncertainty, four versions of Mann-Kendall test with Sen’s slope estimator are used to detect trends in climate change indices in the GLR and surrounding areas.

The structure of the paper is as follows: Data is described in section 2, Methods follow in section 3, Results of trend detection with different methods are described in section 4, and section 5 will summarize the conclusions of this study.

Data

The purpose of this study is to examine trends in climate change indices in the Great GLR and surrounding areas over the past 40 years from January 1981 to December 2020. The reanalysis product ERA5-Land recently released by the European Centre for Medium- Range Weather Forecasts (ECMWF) was suitable for this study. A major advantage of ERA5-Land is the higher temporal and spatial resolution than other global reanalysis products. The ERA5-Land has finer horizontal resolution (9 km) than the ERA5 (30 km) and provides a more accurate description of the water and energy cycles. The assimilation of a much larger number of reprocessed datasets has also improved this reanalysis product [42] compared air temperature data from ERA5-Land, Climatic Research Unit gridded Time Series (CRU_TS), and Global Forecast System (GFS) and found that ERA5- Land and GFS datasets with higher spatial and temporal resolution have comparable performance to CRU_TS data, thus having a greater value in applications.

The study domain is within a rectangle of 100–70°W and 35–65°N (Figure 1) that surround the five interconnected freshwater bodies known as the Great Lakes, and parts of the neighboring Canadian provinces and US states. It includes eight US states: Minnesota, Wisconsin, Illinois, Indiana, Michigan, New York, Ohio, Pennsylvania, and two Canadian provinces: Ontario and Quebec. In this study, the Great Lakes region (GLR) is defined by the red rectangle of 94- 74°W and 40-51°N (Figure 1). The ERA5-Land 2-m air temperature data in the study region is downloaded using the Climate Data Store (CDS) application-programming interface. The data has a spatial resolution of 9 km, and a temporal resolution of 1 h. ERA5-Land does not have minimum and maximum temperature data, so we downloaded the hourly minimum and maximum temperature from ERA5 and interpolated them to the ERA5-Land grid points using bi-linear interpolation; however, at some grids, the interpolated minimum and maximum temperatures are not consistent with hourly average temperature, i.e. the minimum temperature is higher than the average temperature or the maximum temperature is lower than average temperature. To deal with this inconsistence, we take the daily maximum temperature from the 24 hourly 2-m air temperature and the 24 interpolated maximum, and take the daily minimum temperature from the 24 hourly 2-m air temperatures and the 24 interpolated minimum values total precipitation, indicating the total daily precipitation for the previous day.

Based on the daily data, 38 climate change indices (Table A1) are calculated for each year. These indexes include the Expert Team on Climate Change Detection and Indices (ETCCDI) core extreme temperature and precipitation indices [43] and extra indices frequently used in Canada, including the heating degree days (HDD18C), cooling degree days (CDD18C), etc. This study focuses on the long-term trend in annual index time series. Therefore, the period for calculating each index value is one-year from January 1 to December 31 of each year. For the percentile-based indices, we use the most recent 30-year (1991-2020) as reference period rather than 1971-2000 used by many previous studies. We will detect the trends of these indices at each grid point with the four versions of MK methods, and analyze their spatial patterns.

Table A1: Definition of the 38 temperature and precipitation indices analyzed in this study.

The normals defined as the 30-year (1991-2020) averages are calculated using annual averaged temperature (Tm) and total precipitation (Pr). In Figure 1, the shaded colors represent the spatial distribution of the normal of precipitation, and black contours represent the spatial distribution of the temperature normal. In general, normal are determined by physiographic factors, such as elevation, landscape, latitude and longitude, etc. Elevation and latitude can explain much of the spatial variability of temperature and precipitation normal. Other factors that have significant influence on the variability include the position of the grid with respect to lakes (e.g. the Great Lakes) and mountains (e.g. Appalachian Mountains) and, at small/local scale, terrain attributes (aspect and morphology or the relief).Therefore, the isotherms are parallel to the latitude lines and subject to local corrections by other factors. Precipitation is more significantly influenced by atmospheric circulation and large mountains. In the study area, average precipitation decreases from southeast to northwest with maximum values over the Appalachian Mountains.

The normals of these indices are defined as their 30-year (1991- 2020) averages. Generally, heat events indexes TR (Figure A1a), SU (Figure A1e), TX30GE (Figure A1f), TX35GE (Figure A1g) and CDD18C (Figure A2e) decrease with increasing latitude. Temperature indices TXx (Figure A1h), TXn (Figure A1i), TNx (Figure A2a) and TNn (Figure A2b), decrease with increasing latitude either. GSL (FigureA3c) decrease with increasing latitude and cold events FD (Figure A1b), ID (Figure A1c), TN10LT (Figure A1d), HDD18C (Figure A2d) increase with increasing latitude. Over large swaths of the north, some absolute-threshold-based hot events decreased to zero, while in the south, some absolute threshold-based cold events decreased to zero. As a result, these indices are zero for the entire period in most of the study area, therefore trend detection for these indices in these areas is not useful. The percentile-threshold-based indices TX90p (Figure A2f), TX10p (Figure A2g), TN90p (Figure A2h), TN10p (Figure A2i), CSDI (Figure A3a), WSDI (Figure A3b) have almost uniform spatial distribution. Precipitation indices PRTOT (Figure A3f), R1mm (Figure A3g), R10mm (Figure A3h), R20mm (Figure A3i), SDII (Figure A4a), Rx1day (Figure A4b), Rx5day (Figure A4c), R95pTOT (Figure A4g) and R99pTOT (Figure A4i) decrease with increasing latitude.

Methods

Usually, the trend of a time series X (t) is extracted by the ordinary least squares (OLS) model (Solomon et al. 2007; Stocker and coauthors 2014) as

Where β is the linear trend γ, is the intercept and ε is the independent and identically distributed (i.i.d.) noise with zero mean [44-46]. It has been recognized that, if is auto correlated, the effect of autocorrelation in ε lead to an increase of the uncertainty in the trend. The MK test (Mann, 1945) is a non-parametric test that does not require the data to be normally distributed (e.g. most extreme indices are not normally distributed); therefore it has been the most widely used non-parametric test for trend detection in climatological studies. The non-linear trend and the turning point can be derived from Kendall test statistics (Kendall, 1975). Many studies have employed the MK test [47-48]. To deal with the autocorrelation problem, many variants of the MK test have been developed. In this study, we examine trends in annual climate time series using the four most commonlyused MK trend detection methods to address the uncertainties in the detected trends. We define a simple metric to quantify uncertainties. Although the methods have been described in detail in many papers, for completeness, we briefly describe them as follows.

Method 1: The original MK test

The original MK test does not consider serial correlation. It is a ranking-based method that sequentially compares each value of a time series with the rest of the values [49]. The test statistic S is given by:

where

x_j and x_k the sequential data values, and n is the length of the dataset. A positive value of S indicates an uptrend, a negative value of S indicates a downtrend, and zero indicates no trend.

When sample n > 10 (n =40 in this study), the test statistic S follows a normal distribution [5050], with Expectation (E) and Variance (Var) as follows:

Where t_pis the number of data points in the p^th tied group and q is the number of tied groups in the dataset. The standardized test statistic (Z) is calculated as:

Here, the value of Z is the MK test statistic that follows a standard normal distribution with mean being 0 and variance being 1. Confidence intervals of 90, 95 and 99% (p<0.10, p<0.05 and p<0.01, respectively) are taken to classify the significance of positive and negative temperature, precipitation and index trends. In this study, we apply the three intervals to annual temperature and precipitation data to examine the effects of the interval on spatial distribution of significant trends, and take the 95% confidence interval for other indices trend tests.

The slope of the linear trend is estimated using the Theil–Sen estimator, also known as Sen's slope estimator (Sen, 1968). The slope estimates of N pairs of data are first computed by

Where x_jand x_k are data values at times j and k (j>k), respectively. The median of these N (=n*(n-1)/2) values of is Q_i the Sen’s estimator of the slope.

Method 2: The modified MK test using pre-whitening

This test was suggested by Yue and Wang (2002) to use prewhitening (PW) on the time series before the application of trend test. To match the theoretical background of the MK test to actual records, serial correlation elimination is suggested from the original series through the PW treatment, which offsets the serial correlation structure. The trend is less than prior to PW implementation. [51] Have shown that PW procedure is successful only in the case of first-order Auto Regressive (AR) processes with an increasing absolute trend and positive serial correlation coefficient only.

Method 3: The modified MK test using trend-free pre-whitening

This test was also proposed by Yue and Wang (2002) to remove trend components and then PW the time series before application of trend test. This procedure suggests that the removal of trend is necessary prior to the conversion of the time series into independently structured form, and hence, the damage on the trend is reduced significantly,but still there are problems although this approach provides better application opportunity for MK trend test application [52].

Method 4: The hamed and rao modified MK test

This modified MK test was proposed by Hamed and Rao (1998) to address serial autocorrelation issues. They suggested a variance correction approach to improve trend analysis. It considers the first n significant lag auto-correlations; in this study, we used the default settings which considered all significant lags.

Uncertainty in trend detection

The four methods mentioned above can be used to detect trends in time series. To detect trends in indices, we choose the commonly used significant level α=0.05, and use the four methods to detect significant trends for each of the grids. Each method has advantages and disadvantages, and we do not know which method is the best for detecting trends in a specific index. We define a measure to indicate the uncertainty of a detected trend. At each grid point, we count the number of MK test methods (NM) that have detected a significant trend at the point for the index. NM is in the range from0 (no significant trend detected by any of the four methods) to 4 (significant trend detected by all four methods) and indicates how many methods have detected a significant trend at this point. The higher the NM, the less uncertainty in the detected trend.

Result

Autocorrelation analysis results

Autocorrelation in time series could significantly affect trends in time series. When there is a large positive autocorrelation coefficient in time series, the trend will be significant. However, there are two types of auto correlations: as shown in equation (1), one is from the term βt and another is from the noise εt generally, the higher the first term, the larger the autocorrelation and stronger the trend in the time series. In addition, the second term is noise and the autocorrelation due to this term will cause uncertainty of the estimated trends. Before detecting trends in climate indices, we first investigate the autocorrelation in the data. Since our data are annual values, the autocorrelation is weak. Only the lag-1 autocorrelation may be significant, while the long-term lag correlation is negligible in most indices. Therefore, next, we only provide results for lag-1 autocorrelation results.

In the study region, there are 71,633 land grid points that have data for autocorrelation calculation. Grid points on the Atlantic Ocean and Hudson Buy have no data and are outside our study area. We calculate the lag-1 autocorrelation coefficient for each grid point for each index based on 40 years of data, and then count the number of grid points with significant correlations. The threshold used in this study is the 0.05 significance level (95% confidence interval). The corresponding threshold of correlation coefficient is 0.36. Therefore, if the correlation coefficient is greater than 0.36, we consider the lag-1 autocorrelation to be significant, otherwise the correlation is not significant.

To understand the overall autocorrelation situation in the study region we count grid points with significant autocorrelation and express the result as a percentage of grid points. The results are shown in the column “Autocorr” in Table A2, which indicates that most indices have no significant autocorrelation at most grid points (Figure 2). Only six temperature-related indices have significant autocorrelation at 10%~32% of the grid points. They are TN10p, which has 31.1% of grid points with significant lag-1 autocorrelation, followed by TX10p with 18.8%, DTR with 15.8%, T2M with 14.2%, TN90p with 12.3% and HDD18C with 10%.

Table A2: Percentage of grids (PG) with significant trend detected by different methods (Original, Pre_Whiten, Free_pre_whitening and Hamed_Rao), PG with significant lag-1 auto correlation (AutoCorr) and averged trend (TrendAve).

Spatial distribution of autocorrelation coefficients

Figure 3 shows the spatial distribution of lag-1 autocorrelation coefficients of the six indices, which have the most grid points with significant autocorrelation. Grid points with significant autocorrelations at the 0.05 significance level are shaded in the maps. Significant positive autocorrelations are mainly in the area around the Hudson Bay and in south of 40°N for TN90p and east of 80°W for TN10p. The areas with significant autocorrelation are much smaller for the other four indices.

Comparison of trends detected by the four methods

As described in section 3, the main difference between the four methods is in handling autocorrelation in the noise term. To compare the effects of the methods on the significant trend detection, we applied these four methods to annual average temperature and annual total precipitation time series under three significant levels (α=0.01 , 0.05, and 0.1). Figure 4 show that there is a significant increasing trend in temperature in the area around Hudson Buy, Quebec, and the region south of the 40N. Among the four trend detection methods, the Hamed and Rao Modified MK test (method-4) always detect largest areas with significant trend and the Modified MK test using Pre-Whitening method (method-2) always detect smallest areas with significant trend. The reason is that method-4 considered all significant lags and prewhitening processes in method-2 only considered lag-1 auto correlation and consequently underestimated the trends at some locations. In the region around the Great Lakes, the increasing trend does not pass the significant test at the three significant levels in three of the four testing methods. In the region around the Hudson Bay, the increasing trend passes the significant test at significant levels α=0.01 and α=0.05 . Therefore, significance of the detected trends are sensitive to trend

detection methods and significant levels. Figure 5 shows the spatial distribution of trend of annual total precipitation amount. The four methods are not significantly different in precipitation trend detection due to the very weak autocorrelation of the precipitation time series. The common spatial pattern of trends in precipitation is the increasing trend in Quebec with the trend greater than 100mm/40-year and the decreasing trend in Kansas of the USA. In other areas, there is no significant trend for annual precipitation. The significant level affects extension of the area with significant trends. The higher level (e.g.α=0.01 ), the smaller areas and the lower level (e.g.α=0.1 ), the larger areas with significant trends.

Percentage of grid point with significant trend in climate indices

Significance level α=0.05 is usually used in such tests. Therefore, from now on, we use this setting for significant tests. To understand the overall situation in trend detection of the indices, count the number of grid points with significant trend and express the result as a percentage of total grids in the study area. Comparing the results of the four trend detection methods, method-4 detected the most grid points with significant trends, while method-2 detected the least grid points with significant trends. Figure 2 shows that the patterns of percentage variation of the indices are similar to each other among the four methods. We can see that method-1 and method-3 resulted in similar results but with slightly lower percentages for CSDI and CCD18C from method-3. Method-2 produced the lowest percentage for most indices due to the Pre-Whitening process, which weakens the trend because it does not consider signals in higher lags regression terms. Method-4 resulted in the highest percentage for almost all indices because it considers signals in all significant lags. The results show that most of the grid points (approximately 85%) do not have significant trends in 14 out of the 38 indexes. Among the 38 variables, only six temperature variables, which have high autocorrelation coefficients, have more than 50% grids with significant trends (e.g. TN10p,TX10p and CDSI). There are fewer grid points with significant trends for indices that represent ultra-extreme climate events (e.g. TXx, WSDI and TX35GE). Generally, there are more grid points with significant trends for temperature indices than for precipitation indices. In addition, there are more grid points with significant trends for cold event related indices than warm event related indices. Figure 6: Significant trends detected by the four methods. Plotted is NM defined in section 3.5. The shaded areas represent at least one of the four detection methods that have detected a significant trend. Over the lightest red area NM=1, over the darkest red area NM=4, and the white area NM=0.

Spatial pattern of trend in climate indices

The spatial pattern of significant trends in annual temperature and precipitation detected with the four methods at the three significant levels are similar, but the areas with significant trends change with methods and levels. We use NM defined in section 3.5 to measure the uncertainty.

Trend in annual mean temperature and related indices

Figure 6 shows the spatial distribution of NM for 24 temperature related indices. In the region around the Hudson Bay and coastal region of the Atlantic Ocean, the significant increasing trend of Tm (Figure 6a) and decreasing trend of HHD18C (Figure 6n) are significant with high confidence. The uncertainty increases (i.e. NM decreases) in the western areas. The upward trend in the Great Lakes and its immediate vicinity is highly uncertain. Other daily average temperature related indices include CDD18C (Fig 6o), GSL (Figure 6v), GSL-Start (Figure 6w) and GSL-End (Figure 6x). As temperature continues to increase, growing season length increases over the entire study areas due to the earlier growing season start (decreasing GSL_start) and the later growing season end (GSL_end). GSL_end increases much in the northern areas. The GSL trend has large uncertainty over most of the study area. CDD18C increases significantly in the region south of 42.5°N. In the northern areas, because CDD18C are close to zero, there is no trend.

Trend in extreme cold events related indices

All four methods have detected significant decreasing trends in extreme cold events. Figures 6j and 6l show significant increasing trends in the lowest daily maximum temperature (TXn) and the lowest daily minimum temperature (TNn) respectively over most of the study area with high confidence. Correspondingly, Figures 6q and 6s present significant decrease trends in cold days TX10p) and cold nights (TN10p) respectively over almost the entire study area with great certainty. Similarly, the consecutive cold spell index (CSDI, Figure 6t) decreased significantly over most of the study area as well. As temperature increases, frost days (FD, Figure 6c) and icing days (ID, Figure 6d) decrease significantly with great certainty in area north of 60°N. The index TN10LT (Figure 6e) in the far north, southwest and the over the Great Lakes.

Trend in extreme hot events related indices

There are significant increasing trends in extreme hot events over some regions. For example, tropical nights (TR, Figure 6b) and summer days (SU, Figure 6f) have significantly increasing trends in the south of 42°N with great certainty. Hot nights (TN90p, Figure 6r) and hot days (TX90p, Figure 6p) increased significantly over areas north of 60°N, south of 42°N, and eastern regions with high confidence. The increasing trends in TX30GE (Figure 6h) are significant mainly over the southwest areas. Moreover, even the extreme indices TXx (Figure 6i) and TNx (Figure 6k) have significantly increasing trends at some locations. The trend in warm spell duration index (WSDI, Figure6u) is significant over some scattered area with large uncertainty. However, the decreasing trends in extreme cold events are more significant than that the increasing trend in extreme hot events.

Trend in extreme wet events related indices

Figure 7 shows the distribution of the NM for 12 precipitation related indices. Precipitation and wet days have increased in the past 40 years. Among the precipitation indices, Pr (Figure 7a), PRCPTOT (Figure 7b), R10mm (Figure 7d), SDII (Figure 7f), R95pDays (Figure 7k) and R95pTOT (Figure 7l) have significantly increased over Quebec and some other spotty areas. For other precipitation indices, such as R1mm (Figure 7c), R20mm (Figure 7e), Rx1day (Figure 7g), Rx5day (Figure 7h), CDD (Figure 7i) and CWD (Figure 7j), their increasing trends are significant at scattered locations with large uncertainty due to the characteristics of extreme precipitation events, which often happen at local scale with short duration.

Summary and Discussion

In this study, trends in temperature and precipitation related indices over the Great Lakes region and surrounding area from 100°W to 70°W and 35°N to 65°N have been detected and analyzed using the newly released 40-year ERA5-land data. Four versions of the Mann Kendall trend detection methods are used to deal with uncertainty caused by the time series’ autocorrelation. A simple indicator is defined to measure the uncertainty due to different noise processing in the trend detection methods. Significant warming trends are detected by many indices; for example, the extreme cold events significant decreased over most of the study area in the past 40-years with high confidence; specifically, cold nights (TN10p) and cold days (TX10p) have significant decreasing trends. The hot events increased with less significant trends over most areas. Precipitation and some extreme wet event related indices increased significantly over Quebec and some spotty areas; the rest of the wet event related indices do not have significant trends over most of the study area.

Six out of the 38 indices analyzed have more than 50% of the grid points over the study area with significant trends; they are TN10p, TX10p, HDD18C, Tm, TNn and TXn. No significant trend is detected by any methods in some indices may be due to the length of data record not long enough for trend detection as trends are usually weaker for shorter temporal variations. The scattered location with significant trends in precipitation events is due to their high spatial variation; for example, storms often happen at small/local scales.

At large spatial scale, our results are consistent with previous studies using station observation data or gridded reanalysis data at coarser resolutions. This study used higher resolution data and multiple trend detection methods to more accurately locate where the most significant trends are and how confident we are about these detected trends. The warming trend may continue based on future climate projections in Ontario and other regions within this study domain [53]. Such trends will significantly affect all aspects of society. For example, lakes are experiencing less ice cover, with more than 100,000 lakes at risk of having ice-free winters if air temperatures increase by 4C. Ontario’s three principal wine regions are represented by the presence of cool fresh nights during the ripening period. As cool nights continue decrease, wine grape varieties best suited for cool climates may not be suitable anymore. Significant increases in temperature would cause risks of heat stress, therefore restrict productivity and degrade health and wellbeing with significant economic consequences. As both temperature and precipitation increase in the study region, the precipitation increase trend is less significant with significant temporal and spatial variations, which may lead to more flood and/or drought over the study area and cause significant damage to the society. The warming has already caused many extensive wildfires which may happen more frequently in the future as the significant warming trend continues.

While trends of some extreme heat events, such as TR, SU, TX30GE, and TX35GE, are more significant over the southern region, the northern boundaries of these events have moved northward, and would continue to move north ward in the future due to the rising temperature; and their occurrence frequencies will increase at their current locations.

ERA5 data has been widely used in weather and climate change studies. It has been compared with conventional observations and other reanalysis datasets; although it is still not perfectly consistent with conventional observations, it is recognized as one of the best datasets that suitable for climate change studies. In some regions, some absolute threshold-based indices (e.g. TX30GE, TX35GE) remain zero for a long time and then continue to increase, while other absolute threshold-based indices (FD, ID, TN10LT) decrease to zero and then remain zero for a long time. The trend detection methods used in this study may not be applicable to such transitional regions and our future research will address such issue.

1. Austria, Polioptro F, Erick RB (2016) Maximum temperatures and heat waves in Mexicali, Mexico: trends and threshold analysis.Air Soil Water Res9:21-28.

Google Scholar

2. Bloomfield P, Douglas N (1992) Climate spectra and detecting climate change. Clim Change21(3):275-287.

Indexed at, Google Scholar

3. Burn, Donald H. and Mohamed A. Hag Elnur (2002) Detection of hydrologic trends and variability. J Hydrol 255(1-4):107-122.

Indexed at, Google Scholar, Crossref

4. Chattopadhyay S, Dwayne RE (2016) Long-term trend analysis of precipitation and air temperature for Kentucky, United States.Climate4(1):10.

Indexed at, Google Scholar, Crossref

5. Chen F, Philipp GM, Holger K, Tung F, et al. (2022) Trends in auto-correlated temperature series. Theor Appl Climatol147(3):1577-1588.

Indexed at, Google Scholar, Crossref

6. Citakoglu H, Necmiye M (2021) Trend analysis and change point determination for hydro-meteorological and groundwater data of Kizilirmak basin.Theor Appl Climatol145(3):1275-1292.

Indexed at, Google Scholar, Crossref

7. Deng Z, Xin Q, Jinliang L, Neal M, Xiaogang W, et al. (2016) Trend in frequency of extreme precipitation events over Ontario from ensembles of multiple GCMs.Clim Dyn 46(9):2909-2921.

Indexed at, Google Scholar, Crossref

8. Deng Z, Jinliang L, Xin Q, Xiaolan Z, et al. (2018) Downscaling RCP8. 5 daily temperatures and precipitation in Ontario using localized ensemble optimal interpolation (EnOI) and bias correction.Clim Dyn 51(1):411-431.

Indexed at, Google Scholar, Crossref

9. Deng Z, Jinliang L, Xin Q, Xiaolan Z, Hamed B, et al. (2018) Projection of Temperature and Precipitation Related Climatic Design Data Using CMIP5 Multi-Model Ensemble: A case study for Ontario, Canada under RCP 6.0.J Build Sus1(1):39-54.

Indexed at, Google Scholar, Crossref

10. Dergunov AV, Yakubailik OE (2020) Comparative analysis of data on air temperature based on current weather data sets for 2007-2019. Earth Environ Sci 548(3):032034.

Indexed at, Google Scholar, Crossref

11. Douglas EM, Chelsea AF (2011) Is precipitation in northern New England becoming more extreme? Statistical analysis of extreme rainfall in Massachusetts, New Hampshire, and Maine and updated estimates of the 100-year storm. J Hydrol Eng 16(3):203-217.

Indexed at, Google Scholar, Crossref

12. Franzke C (2013) A novel method to test for significant trends in extreme values in serially dependent time series.Geophys Res Lett40(7):1391-1395.

Indexed at, Google Scholar, Crossref

13. Güçlü YS (2018) Multiple Şen-innovative trend analyses and partial Mann-Kendall test.J Hydrol566:685-704.

Google Scholar

14. Gan TY (1998) Hydroclimatic trends and possible climatic warming in the Canadian Prairies. Water Resour Res 34(11):3009-3015.

Indexed at, Google Scholar

15. Hamed KH, Ramachandra R (1998) A modified Mann-Kendall trend test for auto correlated data. J Hydrol 204, no. 1-4: 182-196.

Indexed at, Google Scholar, Crossref

16. Helsel DR, Robert MH (1992)Statistical methods in water resources. Technometrics 49:323-324.

Indexed at, Google Scholar, Crossref

17. Hirsch RM, James RS, Richard AS (1982) Techniques of trend analysis for monthly water quality data.Water resources research18, no. 1: 107-121.

Indexed at, Google Scholar, Crossref

18. Helsel DR, Lonna MF (2006) Regional Kendall test for trend.Environ Sci Technol40(13):4066-4073.

Indexed at, Google Scholar, Crossref

19. Hersbach H, Bill B, Paul B, Shoji H, et al. (2020) The ERA5 global reanalysis.Q J R Meteorol Soc146(730):1999-2049.

Indexed at, Google Scholar, Crossref

20. Hussain Md, Ishtiak M (2019) pyMannKendall: a python package for non parametric Mann Kendall family of trend tests.J Open Source Softw4(39):1556.

Indexed at, Google Scholar, Crossref

21. Kendall M (1975) Rank correlation measures. Charles Griffin London 202:15.

Indexed at, Google Scholar, Crossref

22. Karl TR, David R (1999) Easterling. Climate extremes: selected review and future research directions.Clim Change42(1):309-325.

Indexed at, Google Scholar, Crossref

23. Kim YH, Min SK, Zhang X, Zwiers F, et al. (2016). Attribution of extreme temperaturechanges during 1951–2010. Clim Dyn 46(5-6):1769-1782.

Indexed at, Google Scholar, Crossref

24. Karpouzos DK, Kavalieratou S, Babajimopoulos C (2008) Trend analysis in hydro-meteorological data. Technical Report pg: 5-4.

Indexed at, Google Scholar, Crossref

25. Luo Y, Shen L, Sheng LF, Jingshi L, et al. (2008) Trends of precipitation in Beijiang River basin, Guangdong province, China. Hydrol Int J 22(13):2377-2386.

Google Scholar

26. Mann HB (1945) Nonparametric tests against trend. J Econom 245-259.

Indexed at, Google Scholar

27. Martinez A, Polioptro F, Erick RB (2017) Temperature and heat-related mortality trends in the Sonoran and Mojave Desert region.Atm8(3):53.

Indexed at, Google Scholar, Crossref

28. Meshram SG, Ercan K, Chandrashekhar M, Mohammad AG, et al. (2020) Long-term temperature trend analysis associated with agriculture crops.Theor Appl Climatol 140(3):1139-1159.

Indexed at, Google Scholar, Crossref

29. Min SK, Zhang X, Zwiers FW, Shiogama H, et al. (2013) Multimodel detection and attribution of extreme temperature changes. J Clim 26(19):7430-7451.

Indexed at, Google Scholar, Crossref

30. Muñoz SJ, Emanuel D, Anna AP, Clément A, et al. (2021) ERA5-Land: A state-of-the-art global reanalysis dataset for land applications.Earth Syst Sci Data13(9): 4349-4383.

Indexed at, Google Scholar, Crossref

31. Navarro EJ, Agustín RM, Enrique RV, Jorge EZ, et al. (2018) Observed trends and future projections of extreme heat events in Sonora, Mexico.Int J Climatol38(14): 5168-5181.

Indexed at, Google Scholar, Crossref

32. US and Canada heatwave (2021) US and Canada heatwave: Pacific Northwest sees record temperatures. 5.

Indexed at, Google Scholar, Crossref

33. Vincent LA, Zhang X, Mekis É, Wan H, et al. (2018) Changes in Canada's climate: Trends in indices based on daily temperature and precipitation data. Atmos Ocean 56(5): 332-349.

Indexed at, Google Scholar, Crossref

34. Von SH, Geesthacht H, Navarra A (1995) Misuses of statistical analysis in climate. Springer 11-26.

Indexed at, Google Scholar, Crossref

35. Vose RS, David R, Easterling KE, Kunkel AN, et al. (2017) Temperature changes in the United States.NCA4. GSFC-E-DAA-TN49028.

Indexed at, Google Scholar, Crossref

36. Wang F, Wei S, Haijun Y, Guangyuan K, et al. (2020) Re-evaluation of the power of the mann-kendall test for detecting monotonic trends in hydrometeorological time series.Front Earth Sci8:14.

Indexed at, Google Scholar, Crossref

37. Dos S, Carlos A, Christopher MN, Tantravahi VR, et al. (2011) Trends in indices for extremes in daily temperature and precipitation over Utah, USA. Int J Climatol 31(12):1813-1822.

Google Scholar

38. Rivard C, Harold V (2009) Trend detection in hydrological series: when series are negatively correlated.Hydrol Process An Int J23(19):2737-2743.

Indexed at, Google Scholar

39. Rivard C, Harold V, Andrew RP, Marie L, et al. (2009) Groundwater recharge trends in Canada.Can J Earth Sci46(11):841-854.

Indexed at, Google Scholar, Crossref

40. Sen PK. Estimates of the regression coefficient based on Kendall's tau. (1968) J Am Stat Assoc 63(324):1379-1389.

Google Scholar

41. Şen Z (2017) Hydrological trend analysis with innovative and over-whitening procedures. Hydrol Sci J 62(2):294-305.

Indexed at, Google Scholar, Crossref

42. Solomon S, Martin M, Melinda M, Dahe Q, et al. (2007)Climate change 2007-the physical science basis: Working group I contribution to the fourth assessment report of the IPCC. Cambridge University press 4.

Indexed at, Google Scholar, Crossref

43. Stocker Thomas ed (2014)Climate change 2013: the physical science basis: Working Group I contribution to the Fifth assessment report of the Intergovernmental Panel on Climate Change. Cambridge university press.

Indexed at, Google Scholar, Crossref

44. Sun Y, Hu T, Zhang X, Li C, et al. (2019) Contribution of global warming and urbanization to changes intemperature extremes in Eastern China. Geophys Res Lett 46(11):426-434.

Indexed at, Google Scholar, Crossref

45. Tosunoglu F, Kisi O (2017) Trend analysis of maximum hydrologic drought variables using Mann–Kendall and Şen's innovative trend method.River Res Appl 33(4):597-610.

Google Scholar

46. Yin H, Sun Y (2018) Detection of anthropogenic influence onfixed threshold indices of extreme temperature. J Clim 3(16):6341-6352.

Indexed at, Google Scholar, Crossref

47. Yue S, Chun YW (2004) The Mann-Kendall test modified by effective sample size to detect trend in serially correlated hydrological series. Water Resour Manag 18(3):201-218.

Indexed at, Google Scholar, Crossref

48. Yue S, Paul P, Bob P, George C, et al. (2002) The influence of autocorrelation on the ability to detect trend in hydrological series.Hydrol Process16(9):1807-1829.

Indexed at, Google Scholar, Crossref

49. Zhang XA, Hegerl L, Jones GC, Tank P, et al. (2011) Indices for monitoring changes in extremes based on daily temperature and precipitation data. Wiley Interdiscip Rev Clim Change 2(6):851-870.

Google Scholar

50. Zhang J, Benjamin SF, Tara J (2016) Extreme precipitation drives groundwater recharge: the northern high plains aquifer, central United States, 1950-2010. Hydrol Process 30(14):2533-2545.

Indexed at, Google Scholar, Crossref

51. Zhang Liang (2020) Climate change projections of temperature and precipitation for the great lakes basin using the PRECIS regional climate model.J Great Lakes Res 46(2):255-266.

Indexed at, Google Scholar, Crossref

52. Zhang X, Greg F, Megan KY, Lucie V, et al. (2019) Changes in temperature and precipitation across Canada.CCCR 112-193.

Google Scholar

53. Zhu H, Jinliang L, Xiaolan Z, Xiaoyu C, et al. (2020) The Ontario Climate Data Portal, a user-friendly portal of Ontario-specific climate projections.Sci Data 7(1):1-10.

Indexed at, Google Scholar, Crossref