Historical drought patterns over Canada and their relation to teleconnections 1 2 3

Zilefac Elvis Asong, Howard Simon Wheater, Barrie Bonsal, Saman Razavi, Sopan Kurkute 4 5 6 7 8 Global Institute for Water Security and School of Environment and Sustainability, University of 9 Saskatchewan, 11 Innovation Blvd, Saskatoon, SK, Canada S7N 3H5 10 11 Environment and Climate Change Canada, 11 Innovation Blvd, Saskatoon, SK, Canada S7N 3H5 12 13 14 *Corresponding author: 15 Phone: +1 306 491 9565 16 Email: elvis.asong@usask.ca 17 18


Introduction
Drought is a naturally occurring environmental phenomenon and a major natural hazard that can have devastating impacts on regional water resources, agriculture, industry and other social-ecological systems, with far-reaching impacts in an increasingly globalized and uncertain world (IPCC, 2013;Sternberg, 2011).
Although still among the least understood extreme weather events affecting larger areas, droughts have proved to be the costliest and most widespread of natural disasters (Bryant, 2005;Wilhite, 2000a).This is primarily due to their usually lengthy duration, severity and large spatial extent, sometimes reaching continental scales and lasting for many years (Sheffield et al., 2009).Generally, droughts can affect all components of the hydrological cycle, from its origin as a deficit in precipitation-P (Dai, 2011;Palmer, 1965;McKee et al., 1993), to its combination with high evapotranspiration losses that can lead to a deficit in soil moisture and subsequent manifestation into a hydrological drought (Tallaksen and Stahl, 2014).
A review of key drought concepts (e.g.classification and indices) and the relation between droughts and large-scale climate indices has been carried out by Mishra and Singh (2010).However, due to the wide variety of sectors affected by droughts, their diverse spatial and temporal structures, the inter-dependence across climatic, hydrologic, geomorphic, ecological and societal variables, and the demand placed on water supply by different users, there is no universal definition of droughts and associated impacts.The most used drought classification is that initially proposed by Dracup et al. (1980) and later integrated by Wilhite and Glantz (1985) and Wilhite (2000b); Wilhite and Glantz (1985).Based on the degree of water deficit, droughts are often classified into three types including (1) meteorological, (2) agricultural, and (3) hydrological.Further details on drought classification and definitions are found in (Mishra and Singh, 2010;Dai, 2011;Van Loon et al., 2016).
Studies on regional drought characteristics are important and should be incorporated in water resources management efforts (Mishra and Singh, 2011;Wheater and Gober, 2013).Of particular interest is the analysis of drought occurrence over Canada, a country in which drought is among the most costly natural hazards, particularly in the interior Prairie region, e.g., Bonsal et al. (2011).During the period 1950-2010, nationwide annual mean surface air temperature-T increased by 1.5 o C (Vincent et al., 2012).Being the second largest country in the world and with a large continental interior, this rapid warming has been accompanied by significant changes in many other hydroclimatic elements in different parts of the country, including increases in P (Mekis and Vincent, 2011), decreases in the duration of snow cover (Brown and Braaten, 1998), and decreases in annual streamflow (Zhang et al., 2001).Climate projections also indicate that many regions of Canada will likely experience increasing drought risk by the end of the 21 st century (Masud et al., 2017;Bonsal et al., 2013;Dibike et al., 2017).
Historically, most areas of Canada have experienced periodic droughts with different durations, severities, and marked spatial extent, but the agricultural belt of the Canadian Prairies has tended to be highly susceptible to droughts due in part to its location in the lee of the Rocky Mountains and its strong dependence on rain-fed agriculture (Shabbar and Skinner, 2004).In particular, devastating drought events over western Canada during the 1890s, 1910s, 1930s, 1960s, 1980s, 1999 -2005, and most recently in 2015 have been identified by using a variety of drought indicators at various scales (Bonsal et al., 2011b;Bonsal and Regier, 2007;Szeto et al., 2016).The 1961 drought (the worst single year drought on the Prairies, with about 50% of normal growing season precipitation) led to a total net farm income drop of 48% ($300 million) compared with the previous year (Bonsal et al., 1999).The drought of 1988 had many impacts on the agricultural sectors of Canada, including wind erosion, livestock, incomes, farm management, crop production, and prices.In addition, the sparse snow cover and high spring Ts resulted in little or no spring runoff from Prairie watersheds in 1988, such that the mean runoff volume was 60% to 70% of normal (Wheaton et al., 1992).Furthermore, the 1999 -2005 drought which was at its most severe between 2001 -2002 was felt across Canada but concentrated on the Prairies, and cost the regional economy an estimated $3.6 billion in lost agricultural output (Council of Canadian Academies, 2013).
The uncertainty of drought characterization in Canada in an era of changing climate and increasing pressure from competing water users poses a major challenge to sustainable water management.A better understanding of the spatial distribution of drought, and its frequency, intensity and duration is thus required.Increased knowledge of these drought characteristics and their relationship to large-scale oceanatmosphere forcing is necessary for predicting seasonal drought severity, as well as for planning for impacts Oscillation (PDO) and El Niño-Southern Oscillation (ENSO) have been associated with warm winter T in western and central Canada (Bonsal et al., 2001;Shabbar and Yu, 2012;Shabbar and Khandekar, 1996) and a reduction of snow cover in western Canada (Brown and Braaten, 1998).A review of the association of large-scale variability and low streamflows over Canada was made by Bonsal and Shabbar (2008).They found a higher frequency of low-flow events to coincide with warmer/drier conditions during El Niño events and positive phases of the PDO and the Pacific North American (PNA) pattern.Nazemi et al. (2017) investigated the major drivers of annual streamflow variability in the headwaters of the Canadian Prairies during the 20 th century and found the PDO to significantly determine monotonic trends and shifts in the central tendency of annual mean streamflow.Fleming and Quilty (2006)  inaccuracies associated with cold climate processes (Wang and Lin, 2015;Wong et al., 2016;Asong et al., 2017;Asong et al., 2016a).For this purpose, it is worthwhile to study the long-term time series of P and T regarding their nonhomogeneous climatic and hydrological properties.It is also important to identify homogeneous regions within Canada with distinct drought features for improved drought risk assessment and for a more efficient water resources management at the regional level.So far, most studies on droughts have been limited to Canada south of 60 o N (Masud et al., 2015;Bonsal et al., 2013;Dibike et al., 2017).
Nevertheless, we have not come across studies that attempt to establish the link between nation-wide drought characteristics (e.g.spatial structure, temporal patterns, periodicities) with the large-scale oceanatmospheric modes of variability in a comprehensive manner.
This study aims to fill these gaps by providing a comprehensive analysis of historical droughts over the whole of Canada.Drought events are characterized by the Standardized Precipitation-  (Zhang et al., 2000), this study utilizes and compares two common Canada-wide gridded data sets (monthly P and T) for the period 1950 -2013.
The paper is organized as follows.Section 2 provides a description of the data and analysis methods.
Section 3 discusses the detailed characteristics of drought over different regions of Canada by applying the different aforementioned statistical analyses to the drought index.This section also discusses the physical

Study area
The study area comprises the entire Canadian landmass (Fig. 1).The region includes several major river systems including the Great Lakes -St.Lawrence River system, which is one of the largest freshwater resources globally.Topography plays an important role in shaping regional climates ranging from wet maritime on the coasts, to dry continental across the Prairies and Boreal Plain.Snowfall is restricted to winter months (approximately October to April depending on region).The occurrence, intensity, and timing of seasonal P greatly influence the functioning of ecosystems in various terrestrial ecozones in this region (Hogg et al., 2000).Based on the period 1950 -2013, mean annual P varied from more than 2460 mm on the west and 1260 mm on the east coasts regions, to less than 360 mm in the interior Prairie (southern central) and northern (above 60 o N) regions (Fig. 2).The long-term monthly (January -December) minimum and maximum Ts ranged from 30  to 15  o C (see Section 2.2 below for data sources).
Characterized by a highly variable hydroclimate and diminishing water resources, southern parts of Canada are home to cities with the highest population densities and support a vibrant agro-based economy that was hard-hit by the most severe and prolonged droughts of 1988, 1999 -2005 and 2015, as well as severe floods of 2011, 2013 and 2014 (Wheater and Gober, 2015;Pomeroy et al., 2016).

Gridded observations-ANUSPLIN
The Australian National University Spline (ANUSPLIN) implementation of trivariate thin plate smoothing splines (Hutchinson et al., 2009) 46cc-8990-01862ae15c5f) with over 330 locations for T and 460 for total P (note that these numbers are not constant over time).These data have been quality controlled and adjusted to account for known changes in measurement practices.In particular, records from stations separated by less than 10-km were merged so that correlations between stations would be small.See Mekis and Vincent (2011) for a detailed discussion on merging techniques and trends in the mean climatologies of these data.For CANGRD, these data were interpolated to evenly spaced (50-km) grids using Gandin's optimal interpolation (Gandin, 1966) technique.
As in the case of ANUSPLIN, monthly P, Tmax and Tmin values were extracted from 1950 -2013 and used in the analyses.

Drought index calculation and drought identification
Many quantitative metrics have been developed and used for identification and monitoring of droughts (Mishra and Singh, 2011;Raible et al., 2017).A variety of these indices measure, in most cases, how much P and T for a given period deviate from historical averages.Examples include the Palmer drought severity index (Palmer, 1965), Palmer hydrological drought index (Palmer, 1965), self-calibrated Palmer drought severity index (Wells et al., 2004), standardized precipitation index -SPI (McKee et al., 1993), the SPEI (Vicente-Serrano et al., 2010), and multivariate standardized drought index (Hao and AghaKouchak, 2013).In Canada, most drought analyses (Dibike et al., 2016;Masud et al., 2017) have utilized climate-based indices since the T and P variables are readily available for longer periods and span larger areas, compared to hydrologic variables such as soil moisture and streamflow.In this study, we make use of the SPEI as a meteorological proxy for drought quantification.
As detailed in Vicente-Serrano et al. ( 2010), the SPEI is a multi-scalar drought index based on a water balance approach that uses the monthly difference between P and potential evapotranspiration (PET) to analyze wet/dry spells over multiple time scales.SPEI involves the calculation of monthly PET, and then subtracting this from the corresponding monthly P to obtain the climatic water balance.Several derivations have been put forth for calculating PET, including the widely used Penman (Penman, 1948), Thornthwaite (Thornthwaite, 1948), Priestley-Taylor (Priestley and Taylor, 1972) and Hargreaves (Hargreaves et al., 1985) methods.However, most of these approaches require long records for solar radiation, Ts, wind speed, and air pressure which are not readily available in Canada and many regions of the world.The Hargreaves method (Hargreaves, 1994) which simply uses Tmin and Tmax for estimating PET is employed in this study.Once PET is calculated, the difference between P and PET for the month j is calculated as in Eq. (1): where j Q values represent monthly water surplus or deficit.months are used to represent hydrological drought, and are useful for monitoring surface water resources (Beguería et al., 2014;Hayes et al., 2011).To ascertain the variability of spatiotemporal patterns for different types of droughts, the SPEI was used at different time scales, namely, at 1 (SPEI1), 3 (SPEI3), 6 (SPEI6), and 12 (SPEI12) consecutive months from January -December; at 6 months during the warm season (April -September, SPEI6Apr-Sept); and at 12 months during the hydrologic year (October -September, SPEI12Oct-Sept).SPEI1, SPEI3, SPEI6, and SPEI12 account for the sub-annual variability of droughts while SPEI6Apr-Sept and SPEI12Oct-Sept account for the inter-annual variability.A drought event occurs any time the SPEI is continuously negative and reaches an intensity of -1.0 or less, and ends when SPEI becomes positive.Whenever a drought event has been detected with a start and termination month, drought properties such as duration, magnitude and frequency can be determined.SPEI categories are shown in Table 1.Unless indicated, we focus mostly on other SPEI time scales except SPEI1 which at one month is mainly a meteorological drought index.

Trend Analysis
Long-term trends in drought intensity and variability (inter-annual) are analyzed using the SPEI time series at each grid point.This enables investigation of the percentage of grid points with increasing/decreasing trends during the period 1950-2013 based on the ANUSPLIN-and CANGRDderived SPEI values.Pre-whitening as described in Yue et al. (2002) is first applied to the monthly SPEI anomalies to remove lag-1 auto-correlation, since serial correlation is generally recognized to influence trends in auto-correlated time series, which may distort the power of Mann-Kendal test.Pre-whitening is limited to a low order autoregressive model, i.e., AR(1), since it can falsify the structure of variability in time series across time scales (particularly with higher order models) (Razavi and Vogel, 2018).The pre- trends on a pixel basis.Further details on Mann-Kendall trends can also be found in Hamed and Rao (1998).

Empirical Orthogonal Function Analysis
Hydro-climatological data are often characterized by non-linearity and high dimensionality.Thus, the challenging task is to find ways to reduce the dimensionality of the system to as few modes as possible by expressing the data in such a way as to highlight their similarities and differences (Hannachi et al., 2007).
Empirical Orthogonal Function (EOF) analysis (also known as principal component analysis-PCA) is among the most widely and extensively used method in hydroclimate sciences for accomplishing such tasks (Richman, 1986;Wilks, 2011;Preisendorfer and Mobley, 1988).EOF is employed to find hydroclimate subregions that experienced similar drought features during the period 1950-2013.The EOF consists of computing the covariance matrix of the SPEI series with the corresponding eigenvalues and eigenvectors (Uvo, 2003).
Following the rule by North et al. (1982), the sampling errors at 95% confidence level of the eigenvalues associated with the leading components were estimated in order to establish how many modes to retain for rotation.To achieve more stable and physically explainable patterns, a rotation of the retained components with the varimax procedure (Richman, 1986)  correlation between each REOF and the original SPEI series).Finally, the time variability of the selected RPCs of SPEI were examined for possible trends in the identified sub-regions using linear regression.The slopes of the trends were computed by applying the method of least squares linear fitting to the time series.

Wavelet Analysis
After delimitating hydroclimate sub-regions that experienced similar drought features, the RPC time series corresponding to the REOFs was then analyzed in a time-frequency domain (to reveal dominant and their duration can be distinguished (Fugal, 2009;Grinsted et al., 2004).
The wavelet power spectrum (WPS) is defined as the squared modulus of the CWT (Jiang et al., 2014;Hao et al., 2012), and represents the signal energy at a specific scale (period) and time.From the WPS, the various periodicities and the time intervals of their occurrence can be determined.For a given time series  n y , with 1, 2,3,..., , n N  the CWT is given by: where   n Ws are the wavelet coefficients, n is the time index describing the location of the wavelet in time, s is the wavelet scale, and t  is the sampling interval.The function  is the mother wavelet, and the asterisk (*) denotes its complex conjugate.The Morlet wavelet was used since it provides a good balance between time and frequency domains and is suitable when the purpose is to extract dominant signals (Grinsted et al., 2004).The statistical significance of the CWT was assessed against a red noise background at a 95% confidence interval.The CWT function creates a cone of influence (COI) that delimitates a region of the WPS beyond which edge effects become important and the power could be suppressed (Torrence and Compo, 1998).
Wavelet analysis can also be used to identify the covariance between two time series.This can be done using the concept of wavelet coherence (Grinsted et al., 2004).Wavelet coherence (WCO) reveals local similarities between two time series and may be considered a local correlation coefficient in the timefrequency (time-period) plane.For climatological time series, WCO can be used to identify their possible Ws represent the WPS of the time series   i x and  i y , respectively.The statistical significance (p < 0.05) of the WCO is estimated using Monte Carlo methods with respect to a red-noise spectrum resulting in significant periodicities of coherence delineated by significance contours.As in the CWT, regions outside of the COI should be interpreted with caution (Torrence and Compo, 1998).

Spatial structure of long-term climatic water balance
Fig. 3 shows the average monthly (January -December) water surplus/deficit (mm) as defined in Eq.
(1) for ANUSPLIN and CANGRD.It is evident that, apart from the Pacific/Atlantic maritime, most of the continental Canadian interior experienced moisture deficits, with the Prairie being the most affected region (this region also has the highest climatological Tmin and Tmax, and some of the least mean annual P during the study period-Fig.2).Other than the coastal areas, there is a general east to west moisture deficit pattern, mostly determined by the P and T pattern.For comparison, 5.8% (33.3%) of CANGRD points experienced moisture surplus (deficit), and 4.3% (41.6%) for ANUSPLIN.This implies that ANUSPLIN showed a relatively higher tendency towards dryness relative to CANGRD.Table 2 shows the percentage of grids with decreasing and increasing trends for ANUSPLIN and CANGRD.For both data sets, the percentage generally decreases with increasing SPEI time scale.ANUSPLIN has a higher tendency towards dryness (decreasing trends) unlike CANGRD.Using SPEI12

Long-term trends
as an example, 10% of ANUSPLIN pixels experienced decreasing trends compared to 4% in the case of CANGRD.Conversely, apart from SPEI1, CANGRD showed a strong tendency towards wetness (increasing trends) relative to ANUSPLIN.These differences can be attributed partly to the inputs (e.g. the number of stations considered for gridding ANUPLIN is larger than the number of stations of the CANGRD dataset), estimation methods and spatial resolution which are different for both data sets (see Section 2.2).

Spatial patterns of drought
Concerning the EOF analysis, and taking into account the percentage of variance explained by each rotated component (REOF), two main patterns were identified for subsequent discussion and analysis.Table 3 presents the explained variances of varimax rotated components relative to the considered SPEI time scales and data sets.The first two components (for both data sets) explain about 28% to 33.9% of the total variance depending on the time scale, with the minimum and maximum variances observed for SPEI1 and SPEI12, respectively, in the case of ANUSPLIN.For CANGRD, SPEI6Apr-Sept and SPEI12 explained the minimum and maximum variances, respectively.
Fig. 5 shows that between the two main components (REOF), the regions with higher correlations (>0.7) do not overlap, with clearly spatially disjunctive structure.For all SPEI time scales, the loading patterns of the first component (REOF1) with a maximum explained temporal variance of 20.6% (SPEI12, ANUSPLIN) and 19.5% (SPEI1, CANGRD) highlights mainly the drought evolution on the interior Prairie ecozone of Canada.This semi-arid region is relatively low-lying and characterized by high hydroclimate variability and some of the least annual mean P. The second component (REOF2) explains a maximum variance of around 13.8% (SPEI3) in the case of ANUSPLIN and 13.1% (SPEI12Oct-Sept) for CANGRD.
REOF2 is mainly representative of the Northern central part of Canada (with the least annual mean P), including the Taiga Shield, Taiga Plains, Southern Artic, and Taiga Cordillera ecozones, with positive correlations across all time scales and data sets.In summary, the foregoing analyses suggest that the Prairies and Northern central regions are the leading droughts modes for Canada, in that they are well captured by the two data sets with high positive loadings at all investigated time scales.

Temporal drought characteristics
The RPCs of REOF1 and REOF2 for ANUSPLIN and CANGRD, as well as correlation coefficients between the data sets, are shown in Fig. 6.It is evident that the RPCs from the two data sets are strongly correlated (

cor 
) for all SPEI time scales, except for RPC2 of SPEI3 where 0.45 cor  (the centroids of SPEI3 REOF2 in Fig. 5 are slightly different for both data sets).Primary focus is on SPEI6Apr-Sept and SPEI12Oct-Sept which reflect the agricultural season (when droughts are most critical for rain-fed and irrigated agriculture), and hydrological year in Canada, respectively.Fig. 6 reveals that the Prairies (REOF1 and RPC1) experienced moderate to extreme (Table 1) drought episodes starting in the mid 1950s, early and late 1960s, the 1970s and 1980s, the period from mid 1997 -2005, and the summer of 2009.Some of the drought episodes were extremely dry and severe as they were prolonged in time such as the one in the late 1970s, mid 1980s, and 1997 -2005.These identified drought events correspond well with the findings of Bonsal et al. (2011) who identified large-scale Prairie droughts in 1961, 1967, 1988, and 2001 using ANUSPLIN and CANGRD data.In the Northern central region (REOF2 and RPC2), although fewer droughts were detected for the hydrological year (SPEI12Oct-Sept) with the most intense occurring in 2000 (based on the ANUSPLIN data set), the seasonal SPEI series (SPEI6Apr-Sept) detected several drought events in this region.Extremely dry episodes in this region occurred in the early 1960s, 1980s, early 1990s and 2000s.
Linear trends of the RPCs for the two sub-regions are depicted in Fig. 7 and Table 4.For the Prairie sub-region, both ANUSPLIN and CANGRD showed insignificant decreasing trends for SPEI6, SPEI6Apr-Sept and SPEI12Oct-Sept.However, in the Northern central region, the trend direction is not the same for both data sets at all SPEI time scales.For ANUSPLIN (CANGRD), decreasing (increasing) trends are found in SPEI6Apr-Sept.Conversely, CANGRD (ANUSPLIN) shows decreasing (increasing) trends in SPEI3.
Regionally, SPEI at different time scales tend to display more significant trends in the Northern central relative to the Prairies.Most of the trends are significant in the case of CANGRD compared to ANUSPLIN.The apparent differences in trends between the two data sets and sub-regions may be attributed to the low station density in areas above 60 o N which can introduce higher uncertainty in the gridded P and T fields (Vincent et al., 2015).

Frequency estimation/periodicity of drought
To detect the dominant frequencies in the RPCs of the SPEI series in Fig. 6, the time series were further analyzed using a CWT.The wavelet power spectrum (WPS) of the RPCs is shown in Fig. 8.
Periodicities of drought will be identified for SPEI6Apr-Sept and SPEI12Oct-Sept and their relationship to teleconnection indices examined.These two time scales correspond to the warm season and water year, respectively, and represent periods when moisture shortages are most critical for various sectors in Canada.
Other SPEI time scales which explain sub-annual variability are included in Fig. 8 for the interested reader.
For SPEI6Apr-Sept, from the WPS of RPC1 in Fig. 8, significant interannual variability of between 8 and 32 months is evident throughout much of the entire lengths of the SPEI time series, coinciding with 2-6 year frequencies usually associated to ENSO, e.g., (Shabbar and Skinner, 2004;Gobena and Gan, 2006).However, it is concentrated most heavily during the mid 1950s, early and late 1960s, the 1970s and 1980s, the period from mid 1997 -2010.For SPEI12Oct-Sept, a strong frequency band is centered mostly around 16 -66 months (~4 years), and concentrated during the mid 1950s and late 1960s.Moreover, 32 -64 months frequencies are dominant in the mid 1980s, and the years between 1990 -2010.
In the Northern central region (RPC2), for SPEI6Apr-Sept, significant periodicity of between 8 -40 months cycle as a dominant period of variability is evident.It is concentrated in the early 1960s, 1980s, late 1990s, and 2000s.For SPEI12Oct-Sept, the WPS indicates a significant high power for relatively lowfrequency (16 -100 months, i.e. ~7 years) signals, concentrated over the period 1956 -1975 and 1995 -2013.The foregoing analysis reveals that the Prairie region (REOF1) of Canada is dominated by highfrequency power signals with high cycles of oscillation for both SPEI6Apr-Sept and SPEI12Oct-Sept, while the northern central region (REOF2) is dominated by relatively low-frequency power signals and low cycles of oscillation.The analysis further indicates that significant interannual periodicities (<10 years) dominate drought variability over the two identified sub-regions across Canada.The dominant periodicities and the intervals during which they occurred are summarized in Table 5.

Coherence between drought and large-scale climate drivers
The WCO technique was used to identify both frequency bands and time intervals at which SPEI6Apr-Sept and SPEI12Oct-Sept, and large-scale climate indices are covarying (Torrence and Webster, 1999).The results of WCO coefficients between the RPCs of SPEI6Apr-Sept and SPEI12Oct-Sept and teleconnection indices are shown in Figs. 9 -11.In these plots, the coloured shading represents the magnitude in the coherence as shown in the colour bar, which varies from 0 to 1 and indicates the time scale variability in the correlation between the two time series.As in Fig. 8, the black contours represent the significant regions.The relative phase relationships are shown as arrows, with in-phase pointing right (positive correlation), and anti-phase pointing left (negative correlation), whereas a vertically up (down) arrow indicates that the second time series lags (leads) the first in phase by 90 o .If association exists between two time series, a slowly varying phase lag would be expected, and the phenomena would be phase locked, i.e., phase arrows point only in one direction for a given wavelength (Grinsted et al., 2004;Gobena and Gan, 2006).For this study, phase angle associations were noted strictly as either being in-phase or antiphase locked.To simplify and limit the length of the paper, only results for large-scale climate indices having the strongest correlations with frequencies of drought are shown.Other results are included as supplementary material.1978 -2000 (in the case of ANUSPLIN).As in Fig. 8, it can be seen that the Prairie region has more shortterm periodicities, compared to the Northern central region.It is also evident that although the MEI covaried with drought events, the frequency is higher but short-lived in the Prairies relative to the Northern central region.
Fig. 10 shows that for both the seasonal (SPEI6Apr-Sept) and annual (SPEI12Oct-Sept) drought time series, sporadic but significant coherence is observed intermittently from year-to-year with the PDO.For SPEI6Apr-Sept, dominant strong coherence occurred between 16 -32 months over the period 1988 -2003 for the Prairies region, whereas fluctuation was intermittently observed between 1955 -2000 for SPEI12Oct-Sept over the Prairie and Northern central regions ranging from 16 -64 months.Unlike the PDO, the PNA showed strong in-phase relationship with drought over Canada (Fig. 11).It is apparent that the PNA covaried significantly with SPEI6Apr-SeptPC1 over most of the years during the study period.Oscillations in the PNA are manifested in the SPEI6Apr-SeptPC1 on wavelengths varying from 2-108 months (~9 years), suggesting that the PNA actively mirrors SPEI6Apr-SeptPC1.Particularly, the PNA was phase-locked with drought during the period 1980 -2001 at 16 -32 months scale.Apart from the late 1990s and early 2000s at the 8 -16 months scale, no significant coherence was found between the PNA and SPEI6Apr-SeptPC2.
Although with less coherence relative to SPEI6Apr-Sept, the PNA also co-varied with SPEI12Oct-Sept, more so for the Northern Central region compared to the Prairie region.In the Prairie region, the strongest coherence between SPEI12Oct-SeptPC1 and the PNA occurred over the 16-32 months scale, while for the Northern central region, significantly cross power between the PNA and SPEI12Oct-SeptPC2 occurred over the 32 -64 months scale, spanning the periods 1960 -1973 and 1985 -2000.The relationship between drought and the AMO is shown in Fig. S1.It is clear that the AMO did not co-vary with SPEI6Apr-Sept over both homogenous drought regions.However, it is important to note that Shabbar and Skinner (2004) found significant correlation pattern between the winter AMO index and following summer PDSI in the north of the Prairies provinces.Here, we only made use of the AMO index for the months April -September of each year.For SPEI12Oct-Sept, significant and high-energy existed for the Prairie region only (SPEI12Oct-SeptPC1) and mainly distributed in the 8 -64 months (~5 years) band and spanning the period 1970 -2005.For the AO (Fig. S2), in-phase significant coherence existed with SPEI6Apr-Sept PC2 and SPEI12Oct-SeptPC1.In the Northern central region (SPEI6Apr-Sept PC2), 16 -32 months high-energy regions in this band were found over the period 1963 -1978, and 8 -16 months significant coherence also existed from 1994 -2002.In the Prairie region (SPEI12Oct-SeptPC1), drought was in-phase with the AO especially from 1981 -2003 where significant cross power and coherence mainly concentrated in the 64 -128 months band (based on CANGRD).The results further indicate that the AO was in antiphase with SPEI6Apr-Sept PC1 and SPEI12Oct-SeptPC2 during the study period from 1950 -2013.In terms of the NAO (Fig. S3), sporadic significant coherence is noticed with the seasonal SPEI6Apr-Sept at higher frequency ranging between 4 -20 months and mainly concentrated over the period 1975 -1990.There are also statistically significant in-phase coherence between SPEI12Oct-SeptPC1 and the NAO at around 70 -128 months band in the late 1970slate 1990s, based on CANGRD.
The foregoing analysis have shown significant covariance between drought variability over Canada and large scale climate indices, especially the MEI, PNA and PDO.This is in line with previous studies e.g.(Bonsal and Shabbar, 2008;Fleming and Quilty, 2006;Tremblay et al., 2011;Perez-Valdivia et al., 2012) that have established links between Canadian hydroclimate and teleconnection indices.Furthermore, one should expect most of the climate indices to yield similar findings given that they appear to be interrelated at several time scales (Sheridan, 2003;Gan et al., 2007;Ng and Chan, 2012).We recommend that future studies examine the degree to which such interrelations can affect the findings reported here by eliminating, for example, the influence of the PDO on the MEI via partial wavelet coherence analysis (Ng and Chan, 2012).Also, a lead/lag response of the identified drought frequencies as well their correlations to positive and negative phases of various teleconnections constitute an area for future research.

Summary and conclusions
This study performs a comprehensive analysis of historical droughts over the whole of Canada, of these sub-regions with similar drought variability and characteristics can be useful for drought risk management at a regional scale in Canada.Two of the most important river basins are both in the Prairies region (Saskatchewan River Basin) and Northern region (MacKenzie River Basin).Lastly, this study is the first of its kind to identify dominant periodicities in drought variability over the whole of Canada in terms of when the drought events occur, the duration, and how often they do so over the Prairies and Northern central regions. Hydrol Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.due to future climate change.Previous studies have documented significant links between low-frequency internal climate variability and Canadian hydroclimate.For example, positive phases of the Pacific Decadal quantified the effects of organized modes of climate variability upon groundwater resources, by examining the influence of ENSO on water levels in shallow aquifers in British Columbia.They found water levels to be above average during La Niña years and below average during El Niño years, an indication of variability in winter and spring P that recharges the aquifer systems.Similarly,Tremblay et al. (2011) analyzed the variability of groundwater systems for three different regions across Canada and their linkage to the North Atlantic Oscillation (NAO), the Arctic Oscillation (AO), the PNA, and ENSO.Their findings indicated that groundwater variability in the Prince Edward Island region is mostly modulated by the NAO and AO, in Manitoba it is influenced by the PNA, while for Vancouver Island NAO, AO, and ENSO showed the highest influence.Perez-Valdivia et al. (2012) found variability in groundwater levels in the Canadian Prairies in the 2-7 and 7-10 year bands to be highly influenced by ENSO while oscillation modes in the 18-22 year band reflected a negative correlation with the PDO index.The development of a comprehensive drought monitoring system capable of providing early warning of a drought's onset, severity, persistence, and spatial extent in a timely manner is a critical component in establishing a national drought policy or strategy.However, such nation-wide drought assessments in Canada are hampered partly by observational uncertainties.The paucity and heterogeneous distribution of Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.P and T estimates are an important limitation for drought characterization.Ground-based measurements (e.g.gauges) are limited especially over the Rocky Mountains and north of the 60 th parallel, and suffer from Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.and dynamical mechanisms driving the observed dry episodes in the country.Finally, a summary and conclusions are given in Section 4.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.To compute SPEI, the monthly j Q values are first standardized with respect to the long-term monthly mean values.One-month SPEI is generally representative of meteorological drought, while time scales between 3 and 6 months are considered as an agricultural drought index.Longer scales such as 6 and 12 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.whitened SPEI values on different time scales are then used for detecting statistically significant (p<0.05) was applied.The patterns defined in this way are referred to as rotated empirical orthogonal functions (REOF).In summary, the REOF (spatial patterns) values indicate the spatial representativeness of each rotated temporal component-RPC (temporal patterns).Subsequently, we obtained the most revealing patterns of drought evolution across Canada and determined the spatial extent of each component series by mapping the factorial matrix values (i.e. Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.oscillations) by means of the continuous wavelet transform (CWT).Subsequently, by utilizing the wavelet coherence (WCO) technique (Grinsted et al., 2004;Torrence and Compo, 1998;Addison, 2016), the relationships between the dominant oscillations and large-scale climate indices that have possibly modulated historical drought patterns across Canada are investigated.The CWT projects the spectraltemporal characteristics of a time series onto a time-frequency plane from which the dominant periodicities 2) Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.teleconnection with large-scale atmospheric drivers.The WCO between two time series can be computed by normalizing and smoothing their cross wavelet spectrum:

Fig. 4
Fig. 4 depicts the spatial structure of long-term (1950 -2013) SPEI trends at various time scales.Only grid points with significant trends are shown.It is noteworthy that significant trends largely occur in spatially isolated blocks.Decreasing trends are observable in the southern parts of the Prairies, the foothills of the Rocky Mountains and Pacific maritime regions, whereas increasing trends are limited to a small area in the north, and parts of Atlantic coast to the east, similar to the water balance shown in Fig.3.

Fig. 9
Fig. 9 illustrates the squared WCO between the temporal patterns of drought and MEI.It is apparent that, for SPEI6Apr-Sept, the MEI had significant coherence with PC1 mainly concentrated in the 8 -16 months band and from 1960 -1980 in the Prairie region, indicating that the two time series are phase-locked over this time interval.Also, the strongest coherence between PC1 and MEI occurred over the 32-50 month scale, spanning the period 1986 -2002.In the case of SPEI6Apr-Sept PC2 (Northern central region), discontinuous in-phase coherence patterns were detected in the 2 -32 months bands.The dominant high energy coherence occurred in the 16 -32 month scale, spanning the period 1985 -2005.Similarly, for the hydrological year time scale (SPEI12Oct-Sept), drought in the Prairies (SPEI12Oct-SeptPC1) experienced significant high power with the MEI at the 16 -32 month scale from 1980 -2005 based on ANUSPLIN (1986 -2005 based on CANGRD).Whereas, in the Northern central region (SPEI12Oct-SeptPC2), significant in-phase cross power between drought and the MEI occurred at the 32 -64 month period over the years Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.
considering the role teleconnections by analyzing different monthly P and T products for the period 1950 -2013.SPEI, a climatological drought index, is applied over various temporal scales to evaluate various drought characteristics such as trends, spatiotemporal patterns of long-term change, inter/intra-annual Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.variability, and periodicity/frequency.In addition, potential prominent modes of low-frequency variability such as the PNA, AO, MEI, PDO, AMO and NAO which are well established to influence the hydroclimate of Canada and North America are investigated as precursors to historical drought occurrence.The main conclusions from the analyses are:  Apart from the Pacific/Atlantic maritime regions, most of the continental Canadian interior experienced moisture deficits, with the Prairie region being the most affected region between 1950 -2013.Based on a trend analysis derived from the water balance, significant decreasing trends in SPEI values are observed in the southern Prairies, the foothills of the Rocky Mountains and Pacific maritime regions, whereas increasing trends are limited to a small area in the north, and parts of Atlantic coast to the east.Therefore, southern parts of the country showed a trend towards drier conditions and vice versa for the northern regions.Note that the northern region (above 60 o N) has lower station density and as such higher uncertainty in the gridded P and T fields. EOF identifies two main spatially disjunctive sub-regions of drought variability over Canada -the Prairie and Northern central regions.Based on both seasonal (SPEI6Apr-Sept) and annual (SPEI12Oct-Sept) SPEI values, the Prairie sub-region experienced moderate to extreme droughts episodes starting in the mid 1950s, early and late 1960s, the 1970s and 1980s, the period from mid 1997 -2005, and the summer of 2009.Some of the drought episodes were extremely dry and severe as they were prolonged in time such as the ones in the late 1970s, mid 1980s, and 1997 -2005.In the Northern central region, although fewer (likely due to below 0 temperatures for most of the cold season) droughts were detected for the hydrological year (SPEI12Oct-Sept) with the most intense occurring in 2000, the seasonal SPEI series (SPEI6Apr-Sept) detected extremely dry episodes in this region in the early 1960s, 1980s, early 1990s and 2000s.However, drought variability in the Prairies in particular experienced largely insignificant trends, a finding similar to numerous previous studies. Wavelet analysis was particularly useful to detect periodical signals in the SPEI time series patterns for each sub-region.For SPEI6Apr-Sept, the analysis reveals clearly the presence of a dominant periodicity of between 8 and 32 months, persisting approximately from 1955 to 2001 in the Prairie Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.region, while in the North central region, significant periodicity of between 8 -40 months cycle as a dominant period of variability spanning the years 1955 -2000 is apparent.For SPEI12Oct-Sept, a strong power frequency band over the Prairie region, centered mostly around 16 -66 months (~4 years) and spanning the period 1955 -1968 is found.Moreover, 32 -64 months periodic high-power signals are dominant during the years 1970 -2002.In the Northern central region, a significant high power for relatively low-frequency (16 -100 months, i.e. ~7 years), spanning the period 1956 -1975 and 1995 -2013 is detected.Therefore, the Prairie region of Canada is dominated by high-frequency (i.e.more frequent and shorter cycles of dry events) power signals for both SPEI6Apr-Sept and SPEI12Oct-Sept, while the Northern central region is dominated by low-frequency (i.e. less frequent cycles of dry events) power signals.The analysis further indicates that significant interannual periodicities (with a period of <10 years) dominate drought variability over the two identified major regions of drought variability across Canada.The identified drought short-time (long-time) interannual periodicities in the Prairie (Northern central region) are likely associated with the immediate and significant influence of the MEI and PNA in particular as these large scale climate indices have maximum regions with a 5% significance level in the WCO plots.For the MEI index, 8 -16 and 32 -50 months was the most predominant and effective period on the drought occurrence in Canada, while for the PNA, 2-108 months period was the most predominant over the Prairie region and for SPEI6Apr-Sept compared with the Northern central region, while for SPEI12Oct-Sept, the PNA showed more coherence with drought in the 32 -64 months scale.The foregoing analysis has indicated the need to consider various observational data sets in drought characterization, given the uncertainty in data.In terms of trends, the ANUSPLIN data set indicated a higher tendency for drought over the study period relative to CANGRD.Furthermore, irrespective of the time scale of accumulation, ANUSPLIN tends to reveal more drought severity compared to CANGRD although the correlation between the time series of the two data sets from each of the homogenous drought subregions is very strong.Therefore, further applications using other gridded data sets to verify the role played by the spatial resolution of the input data on regional drought patterns are recommended.The identification Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-122Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 14 March 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 1 .
Figure captions Figure 1.Study area showing topographic and hydrographic features of Canada.

Figure 4 .
Figure 4. Trends in SPEI at each grid point for ANUSPLIN and CANGRD.Only significant values are shown on the map.Brown indicates decreasing and green is increasing.a=SPEI at 1 month, b=SPEI at 3 months, c=SPEI at 6 months, d=SPEI at 12 months, e=SPEI at 6 months from April -Sept, f=SPEI at 12months of the water/hydrologic year i.e.Oct -Sept.Trends are significant at 95% confidence level.

Figure 5 .
Figure 5. Spatial patterns of the first two REOFs of SPEI at various time scales.The spatial extent of the first two REOFs was characterized by mapping the values of the factorial matrix.See Table3for information on variances explained by each REOF pattern.

Figure 6 .
Figure 6.Temporal patterns (RPCs) of the first two rotated principal components (REOFs) of SPEI at various time scales.Indicated within each box is the pattern correlation between ANUSPLIN and CANGRD.

Figure 8 .
Figure 8. Wavelet power spectrum of the time series (RPCs) shown in Fig. 6.The black contour designates the 95% confidence level against red noise, and the cone of the influence (COI) where edge effects might distort the picture is shown as a lighter grey shade.

Figure 9 .
Figure 9. Squared wavelet coherence between the MEI and the temporal patterns of drought (SPEI6Apr_Sept and SPEI12Oct_Sept).Phase arrows pointing right indicate signals are in phase, whereas a left-pointing arrows indicate an antiphase relationship.Arrows deviating from the horizontal are indicative of lead-lag relationships between the two signals.The black contour designates the 95% confidence level against red noise, and the cone of the influence (COI) where edge effects might distort the picture is shown as a lighter grey shade.

Figure 10 .
Figure 10.Squared wavelet coherence between the PDO and the temporal patterns of drought (SPEI6Apr_Sept and SPEI12Oct_Sept).Phase arrows pointing right indicate signals are in phase, whereas a left-pointing arrows indicate an antiphase relationship.Arrows deviating from the horizontal are indicative of lead-lag relationships between the two signals.The black contour designates the 95% confidence level against red noise, and the cone of the influence (COI) where edge effects might distort the picture is shown as a lighter grey shade.

Figure 11 .
Figure 11.Squared wavelet coherence between the PNA and the temporal patterns of drought (SPEI6Apr_Sept and SPEI12Oct_Sept).Phase arrows pointing right indicate signals are in phase, whereas a left-pointing arrows indicate an antiphase relationship.Arrows deviating from the horizontal are indicative of lead-lag relationships between the two signals.The black contour designates the 95% confidence level against red noise, and the cone of the influence (COI) where edge effects might distort the picture is shown as a lighter grey shade.

Figure 5 .
Figure 5. Spatial patterns of the first two REOFs of SPEI at various time scales.The spatial extent of the 885 first two REOFs was characterized by mapping the values of the factorial matrix.See Table3for 886 information on variances explained by each REOF pattern.887

Figure 6 .
Figure 6.Temporal patterns (RPCs) of the first two rotated principal components (REOFs) of SPEI at various time scales.Indicated within each box is the pattern correlation between ANUSPLIN and CANGRD.

Figure 9 .
Figure 9. Squared wavelet coherence between the MEI and the temporal patterns of drought (SPEI6Apr_Sept and SPEI12Oct_Sept).Phase arrows pointing right indicate signals are in phase, whereas a left-pointing arrows indicate an antiphase relationship.Arrows deviating from the horizontal are indicative of lead-lag relationships between the two signals.The black contour designates the 95% confidence level against red noise, and the cone of the influence (COI) where edge effects might distort the picture is shown as a lighter grey shade.

Figure 10 .Figure 11 .
Figure 10.Squared wavelet coherence between the PDO and the temporal patterns of drought (SPEI6Apr_Sept and SPEI12Oct_Sept).Phase arrows pointing right indicate signals are in phase, whereas a left-pointing arrows indicate an antiphase relationship.Arrows deviating from the horizontal are indicative of lead-lag relationships between the two signals.The black contour designates the 95% confidence level against red noise, and the cone of the influence (COI) where edge effects might distort the picture is shown as a lighter grey shade.

Table 2 :
Percentage of grids with decreasing (Dec.) and increasing (Inc.)trends for the SPEI 800

Table 3 :
Percentage of variance explained by the first two varimax rotated loadings (REOFs) of the SPEI 802 at various time scales computed using observations from ANUSPLIN and CANGRD data sets 803