Tree-ring-based estimates of long-term seasonal precipitation in the Souris River Region of Saskatchewan, North Dakota and Manitoba

Abstract

Historically unprecedented flooding occurred in the Souris River Basin of Saskatchewan, North Dakota and Manitoba in 2011, during a longer term period of wet conditions in the basin. In order to develop a model of future flows, there is a need to evaluate effects of past multidecadal climate variability and/or possible climate change on precipitation. In this study, tree-ring chronologies and historical precipitation data in a four-degree buffer around the Souris River Basin were analyzed to develop regression models that can be used for predicting long-term variations of precipitation. To focus on longer term variability, 12-year moving average precipitation was modeled in five subregions (determined through cluster analysis of measures of precipitation) of the study area over three seasons (November–February, March–June and July–October). The models used multiresolution decomposition (an additive decomposition based on powers of two using a discrete wavelet transform) of tree-ring chronologies from Canada and the US and seasonal 12-year moving average precipitation based on Adjusted and Homogenized Canadian Climate Data and US Historical Climatology Network data. Results show that precipitation varies on long-term (multidecadal) time scales of 16, 32 and 64 years. Past extended pluvial and drought events, which can vary greatly with season and subregion, were highlighted by the models. Results suggest that the recent wet period may be a part of natural variability on a very long time scale.

En 2011, il y avait des inondations historiques sans précédente dans le bassin de la rivière Souris de Saskatchewan, du Dakota du Nord et du Manitoba. Elles sont arrivées pendant une période à long terme des conditions humides dans le bassin. Pour développer un exemple des futures tendances, il faut évaluer les effets pluridécennaux de la variabilité du climat au passé et/ou la possibilité de changements climatiques sur la précipitation. Dans cette étude, on a analysé des séries dendrochronologiques et des données historiques sur les précipitations dans la région mentionnée pour développer des modelés de régression qu’on peut utiliser pour prédire les variations à longue-terme sur les précipitations. Pour mettre l’accent sur la variabilité à longue terme, on a modélisé la moyenne glissante des précipitations de 12 ans en cinq sous-régions (déterminées par l’analyse de groupage de mesures de la précipitation) de la partie supérieure du bassin de la rivière Souris pendant trois saisons (novembre à février, mars à juin et juillet à octobre). Les modèles s’en sont servis de multirésolution de décomposition (une décomposition additive à base de puissances de deux en utilisant une transformée en ondelette discrète) des séries dendrochronologiques du Canada, et des États-Unis et la moyenne glissante des précipitations saisonnières basée sur les données du climat canadien ajustées et homogénéisées et les données climatiques du US Historical Climatology Network data. Les résultats montrent que la précipitation varie sur des échelles de longue durée (pluridécennal) de 16, 32 et 64 ans. Les évènements antérieurs pluviaux et des sècheresses prolongées, qui peuvent varier considérablement selon la saison et la sous-région, ont été mis en évidence par les modèles. Les résultats suggèrent que la période humide récente peut être une partie de la variabilité naturelle sur une échelle de très longue durée.

Introduction

Historically unprecedented flooding occurred in the Souris River Basin of Saskatchewan, North Dakota and Manitoba in 2011, and caused extensive damage to Minot, North Dakota, and numerous smaller communities in the Basin. Near-record snowfall occurred over parts of the region, and “record-setting rains” occurred in May and June (Service Assessment Team 2012). In addition, the stage had been set by wet antecedent conditions for a very large flood. There was an “anomalously wet fall season” (Service Assessment Team 2012) that continued a period of more than a decade of wet conditions, which filled available soil moisture storage, shallow groundwater aquifers, wetlands and lakes. Observations from the Gravity Recovery and Climate Experiment (GRACE) satellites show that from 2003 to 2012, parts of the basin experienced a fresh water storage rate of change of an additional 3 cm per year (National Aeronautics and Space Administration [NASA] Earth Observatory 2013). The severe flooding prompted the International Souris River Board to create a Souris River Flood Task Force, which prepared a plan of study for evaluating potential reservoir operation changes and flood control measures to manage future floods and droughts (International Souris River Board 2013). The task force plan indicated a need to evaluate the effects of multidecadal climate variability on future flood and drought risk.

As indicated by Razavi et al. (2015, p. 1813), there is a need to “broaden the understanding of hydrologic characteristics in any basin beyond the limited observational records.” The purpose of this work is to evaluate available records from meteorological stations and tree-ring climate proxy data from the study area to document low-frequency climate variability in the region surrounding the Souris River Basin. This work contributes to a better understanding of climate variability in the Souris River Basin and helps in developing a stochastic model for simulating future streamflows that are consistent with climatic uncertainty, cover the full range of possibilities from drought to pluvial periods, and provide estimates of future risk. The stochastic model under development can use the observational precipitation record to model high-frequency variability; however, this work was undertaken to better describe the low-frequency variability captured in much longer tree-ring climate proxy data.

Definition of study area

To understand regional long-term precipitation, the study area extends beyond the Souris River Basin boundary and includes parts of other basins in the region, including the Missouri River, Red River of the North (Red River), Qu’Appelle River, Assiniboine River and Saskatchewan River. The potential meteorological stations and the tree-ring sites were placed on the map, and a regional study boundary was determined by examining one- to four-degree buffers around the basin in relation to ecological regions (Commission for Environmental Cooperation 1997) in order to obtain a regional view of climate. The final study area was a four-degree buffer around the Souris River Basin (Figure 1). This places the southwestern boundary just north of the Black Hills of South Dakota and Wyoming – a mountainous area. The western boundary also results in the selection of meteorological stations east of the Cypress Hills region of Saskatchewan, an area in the Northwest Forested Mountains ecological region (Commission for Environmental Cooperation 1997). The four-degree boundary extends east into Minnesota, but excludes the forested Lake Superior region of Minnesota. The four-degree boundary provides maximal climate data and a larger set of potential tree-ring data sets in a fairly homogenous plains and prairies setting. The boundary, meteorological stations and potential tree-ring sites are shown in Figure 1.

Figure 1. Locations of the Souris River Basin, study area boundary (four-degree buffer around Souris River Basin), meteorological stations and potential tree-ring sites. The tree-ring sites ultimately used are labeled with the site names.

Data and methodology

Climate data were downloaded from two different national repositories because the Souris Basin is an international basin. Monthly precipitation, average temperature, maximum temperature and minimum temperature data were downloaded from the Adjusted and Homogenized Canadian Climate Data (AHCCD) website on 9 July 2014, from links at https://ec.gc.ca/dccha-ahccd/default.asp?lang=En&n=B1F8423A-1 (Environment Canada 2014). The AHCCD web site stated:

These data were created for use in climate research including climate change studies. They incorporate a number of adjustments applied to the original station data to address shifts due to changes in instruments and in observing procedures. Sometimes the observations from several stations were joined to generate a long time series. (Environment Canada 2014)

The data are further described by Mekis and Vincent (2011) and Vincent et al. (2012).

Monthly precipitation, average temperature, maximum temperature and minimum temperature were downloaded from the United States Historical Climatology Network (USHCN) on 21 May 2014, from the ftp site http://cdiac.ornl.gov/ftp/ushcn_v2.5_monthly/ (Menne et al. 2014). The USHCN web site provides further information about the data:

The United States Historical Climatology Network (USHCN) is a high-quality dataset of daily and monthly records of basic meteorological variables…. The USHCN has been developed over the years at the National Oceanic and Atmospheric Administration’s (NOAA) National Climatic Data Center (NCDC) to assist in the detection of regional climate change. Furthermore, it has been widely used in analyzing US climte [sic]. The period of record varies for each station. USHCN stations were chosen using a number of criteria including length of record, percent of missing data, number of station moves and other station changes that may affect data homogeneity, and resulting network spatial coverage. (Carbon Dioxide Information Analysis Center 2014)

Data for each of the AHCCD and USHCN meteorological stations were summarized as seasonal total precipitation and seasonal average maximum, minimum and mean temperatures. The seasons were defined starting with winter (season 1) as November–February because this is the period when much of the precipitation is stored in the basin in frozen soils or snowpack. The other two seasons followed from the winter definition. March–June (season 2) is generally the period when peak streamflow occurs, and July–October (season 3) represents summer and fall conditions.

Tree-ring data were downloaded from The International Tree-Ring Data Bank (ITRDB) website, http://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring (NOAA NCDC 2014) on 16 July 2014. Tree-ring sites within and near the study boundary were considered for inclusion in the analysis (Figure 1).

Hierarchical agglomerative cluster analysis

Since individual meteorological stations are spatially correlated, the goal was not to develop models for each site, but to develop models for climatically similar subregions. The subregions were determined using hierarchical agglomerative cluster analysis (HACA). In HACA, each observation (as described by the all the variables associated with it – in this case, the meteorological variables associated with a meteorological station) forms its own cluster, and then pairs of clusters are successively merged on the basis of a similarity measurement and a linkage method. There are N − 1 (where N is the total number of observations) merges in which the closest two clusters are merged into a single cluster, resulting in one less cluster at that merge. In this analysis, the similarity was computed by the Euclidean distance between observations. The HACA routine is performed by the function agnes in the cluster package (Maechler et al. 2014) of R (R Core Team 2014).

The linkage method used to merge clusters was Ward’s method, which uses an analysis of variance (ANOVA) approach to evaluate differences between clusters (Güler et al. 2002). The two clusters that are merged are the pair “that leads to the smallest increase in the sum of the within-group sums of squares” (Insightful Corporation 2001, p. 141). The within-group sum of squares is the sum of the squared Euclidean distances between observations at the center of its parent group.

The dendrogram reports the agglomerative coefficient (AC). The AC is a dimensionless quality index for measuring the clustering structure of the data set, and is between 0 and 1. An AC close to 1 indicates that there is a “very clear clustering structure” in the data; however, the structure needs to be analyzed to determine that it is reasonable (Kaufman and Rousseeuw 1990, p. 213), such as replicating the east–west gradient in precipitation in the study region.

The potential clustering variables were chosen to represent long-term differences between meteorological stations, not year-to-year differences, and included mean and median seasonal temperature, and mean and median seasonal total precipitation. The sites have varying numbers of observations and varying start and end dates, some missing years, and some Canadian sites were operated seasonally for parts of their period of record. Therefore, the cluster variables were based on data from the period 1900–2010 for sites with at least 60 years of seasonal values. This removed parts of the beginning and end of the record for some sites in order to have a more consistent period of record for the cluster analysis.

Tree-ring analysis

Potential data sets were initially selected for maximum record length. The ring-width series (the series of growth measurements from individual trees at a particular site) from potential data sets in the ITRDB were imported, plotted and, detrended, and mean-value chronologies (the site-level time series) were built using the Dendrochronology Program Library in R (dplR) – a package for dendrochronologists to handle data processing and analysis (Bunn 2008, 2010; Bunn et al. 2014) for the software R (R Core Team 2014).

The detrending, or standardization, process removes variability that is caused by non-climatic tree-stand or biological effects (such as differing growth rates as the tree ages) and is an important part of dendroclimatic studies (Cook and Kairiukstis 1990; Fritts 2001). The detrending method used was that of fitting a negative exponential curve to each series. This is “probably the most common method for detrending” (Bunn et al. 2014, p. 6) and can recover annual, decadal, to multicentennial signals in tree-ring chronologies (Bunn et al. 2004). The model is where is the ring-width growth at time is the initial year, the exponent is the slope of the decrease in growth of the ring widths and is the growth per year after the initial growth spurt (Fritts et al. 1969; Bunn et al. 2004). Building a mean-value chronology is a process of averaging across the detrended ring-width series at each site (Fritts 2001; Bunn 2008, 2010; Bunn et al. 2014).

It readily became apparent that some of the tree-ring chronologies seemed promising from the standpoint of a long period of record; however, some were based on tree-ring series with less overlap than others, creating small sample sizes for parts of the record. Therefore, the common interval for each set of tree-ring width series was examined. There are three ways to do this: (1) find the common interval that maximizes the number of series; (2) find the common interval that maximizes the number of years; or (3) find the common interval that is the best compromise between the two (Bunn et al. 2014). The compromise option was used to identify those series with a common interval of less than 50 years, which were dropped from further consideration.

Correlation of segments within potential data sets was examined using the corr.rwl.seg function in the dplR package (Bunn et al. 2014). This function calculates “correlation serially between each tree-ring series and a master chronology built from all the other series” in the raw ring-width series. “Correlations are done for each segment of the series where segments are lagged by half the segment length” (Bunn et al. 2014), with the segment length used being the default 50 years. “Correlations are calculated for the first segment, then the second segment and so on. Correlations are only calculated for segments with complete overlap with the master chronology” (Bunn et al. 2014, p. 19). This analysis highlighted some low correlations and led to an examination of the comments in the associated correlation statistics files in the ITRDB, where available. The result of this investigation was the removal of segments from one tree-ring width series (Burning Coal Vein), as suggested by the comments (Meko and Sieg 2014a), and then detrending the revised series.

A measure of reliability of the tree-ring signal is the express population signal (EPS) or subsample signal strength (SSS; EPS calculated over subsamples of the tree-ring period of record). These were defined by Wigley et al. (1984) and are available in dplR (Bunn et al. 2014). For the remaining possible tree-ring chronologies, EPS over 30-year moving windows that overlap by 15 years (Büntgen et al. 2012), the number of cores for each window, number of trees, and the signal to noise ratio (SNR) were examined. Examination of the EPS, SSS and SNR highlighted the shortcomings of some of the series in that their period of record seemed promising from a total number of years standpoint, but only had one core at the beginning of the record or periods within the record of low SSS or SNR. The sites considered potential candidates for modeling precipitation, based on SSS being above 0.85, were:

•	Boundary Bog – Tamarack (Larix laricina; MacDonald and Case 2014);
•	Burning Coal Vein – Ponderosa Pine (Pinus ponderosa Douglas ex C. Lawson; Meko and Sieg 2014a);
•	Theodore Roosevelt National Park – Rocky Mountain Juniper (Juniperus scopulorum Sarg.; Meko and Seig 2014b); and
•	Cedar Butte – Ponderosa Pine (Pinus ponderosa Douglas ex C. Lawson; Meko and Sieg 2014c).

These sites are labeled in Figure 1. SSS was also used to define truncation dates for reconstruction. The final period for reconstruction was determined to be 1700 to 1990.

This resulting set of potential tree-ring chronologies represents three different tree species that may respond to precipitation differently (Wettstein et al. 2011) and that might contribute differing explanatory information to the models. The study area is represented well by the inclusion of Boundary Bog at the northern edge of the study area, Cedar Butte at the southern edge, and the sites Theodore Roosevelt National Park and Burning Coal Vein closer to the center of the study area. Case and MacDonald (2003) used another tree-ring chronology (Picea mariana) from Boundary Bog to reconstruct Saskatchewan River annual streamflow at the Saskatchewan–Manitoba border, and this location is within cluster 5 in our study. Razavi et al. (2015) found that Larix laricina had a statistically significant correlation with Red Deer River and Bow River average water year flows in Alberta, directly west of cluster 4 in our study. Mamet and Kershaw (2011) compared the effects of climate on radial grow of Larix laricina, Picea glauca and Picea mariana in Manitoba and found that the species most sensitive to climate was Larix laricina. Pinus ponderosa has previously been used to identify periods of drought in North Dakota, South Dakota and Montana (Meko 1982; Stockton and Meko 1983). The three sites and two species representing chronologies from Meko and Sieg (2014a, 2014b, 2014c) were examined for their dendroclimatic potential by Sieg et al. (1996), along with many of the other sites initially considered in this study. Based on field research and statistical analyses (which did not include Larix laricina), Sieg et al. (1996, p. 301) stated that

the greatest potential for multi-century dendroclimatic reconstructions from trees in the northern Great Plains [of the United States] is limited to the western part of the Dakotas, and that efforts there are most likely to be successful if concentrated on ponderosa pine and Rocky Mountain juniper.

They also examined bur oak, indicating that it had “excellent” potential for dendroclimatology at some sites “if the time period of interest is restricted to post-1800” (Sieg et al. 1996, p. 301).

Modeling precipitation

A number of studies have identified a decadal-scale signal in precipitation in North America (10–20 years; Ault and St. George 2010), in western North America (greater than a 7-year period; Cayan et al. 1998), in the United States (in fall precipitation with a period of 12 years; Small and Islam 2008, 2009), in the Great Plains (Garbrecht and Rossel 2002) and in the central US (interdecadal variations of annual precipitation on 12- and 20-year time scales; Hu et al. 1998). Yang et al. (2007) found that September, October and November precipitation was strongly related to Niño-3.4 sea-surface temperatures on a 5.5–8.5-year timescale. Small and Islam (2008, 2009) identified a statistically significant signal in autumn precipitation in the central US with a periodicity of approximately 12 years, with the signal strongest in the Midwest and Great Plains. Ryberg et al. (2014) focused on long-term (multidecadal) variability and thus shorter term, quasi-periodic signals were treated as “nuisance” variability and smoothed out by using a 12-year moving average when analyzing precipitation, potential evapotranspiration and runoff. Therefore, in this study 12-year moving average precipitation was modeled using the tree-ring chronologies.

Annual tree-ring chronologies (commonly used to model annual or seasonal precipitation) may not adequately represent 12-year average precipitation and, as Razavi et al. (2016) indicate, annual tree-ring series contain information on longer than annual time scales; therefore, variables representing the multiresolution decomposition of the tree-ring chronologies were used. In multiresolution analysis of a time series, specific frequency components are extracted from a time series (Bunn et al. 2014). The example plot (Figure 2) shows a one-dimensional multiresolution analysis using a level-J additive decomposition of the time series using the pyramid algorithm (Mallat 1989) as implemented in the dplR package (Bunn et al. 2014). This is an application of wavelet analysis and the additive decomposition is for each power of two in the period of record, where the number of powers is trunc(log(nYrs)/log(2)) − 1, in which nYrs is the number of years in the tree-ring chronology. The R package waveslim (Whitcher 2013) was used to do the multiresolution decomposition. Each frequency component is scaled for plotting by dividing by the standard deviation. Supplemental figures in Appendix B of Ryberg (2015) show the multiresolution decomposition for the three other tree-ring chronologies.

Figure 2. Example multiresolution decomposition of tree-ring chronology at Boundary Bog, Saskatchewan, Canada.

In Figure 2, the line labeled D1 on the left and 2 years on the right represents the short-term variability over a 2-year time scale. D7 represents the variability over a 128-year time scale. This is an additive decomposition, so adding D1 through D7, plus a shift (not shown), produces the original tree-ring chronologies. The interest here is in the longer term variability, beyond 12 years; therefore, the D1, D2 and D3 wavelet voices were not considered, leaving the D4, D5, D6 and D7 wavelet voices as potential explanatory variables.

To further examine the frequency bands expressed in the tree-ring chronologies, a continuous wavelet analysis was done using a continuous Morlet wavelet transform (Bunn et al. 2014; Torrence and Compo 1998). Figures B.5–B.8 in Appendix B of Ryberg (2015) depict continuous wavelet analysis, and the colored portion of each figure shows a measure of wavelet power relative to white noise. Each continuous wavelet figure shows that the four sites have power in the 16-, 32-, 64- and 128-year time bands for at least part of the period of record. However, all or almost all of the 128-year time bands are within a cone of influence that occurs at the beginning and end of the time series, where edge effects are an issue. This indicates that the lengths of the tree-ring chronologies used are not long enough to have a reliable 128-year wavelet voice. Therefore, to model precipitation the D4 (16 years), D5 (32 years) and D6 (64 years), wavelet voices were used as potential explanatory variables.

Model selection

All subsets regression (Lumley 2009) was used to examine potential 12-year moving average precipitation models using the D4, D5 and D6 subset of the wavelet voices of the four chronologies. All-subsets regression performs an exhaustive search for the best subsets of the potential explanatory variables for predicting the response with a multiple regression model. The procedure examines all models with one explanatory variable and picks the best fit (or more than one depending on how many the user chooses to save), then examines all models with two explanatory variables and so on, up to the total number of possible explanatory variables or an upper limit set by the user. It was assumed that all wavelet voices of the four chronologies were potentially correlated with the seasonal precipitation in the five groups. The all-subsets regression modeling process dropped out the wavelet voices not significantly correlated with a particular group and season.

All-subsets regression provides measures of model quality: adjusted coefficient of multiple determination (), Schwartz’s information criterion (BIC), and Mallow’s . allows for the comparison of models that have differing numbers of explanatory variables by penalizing models that have additional coefficients (Helsel and Hirsch 1995). The closer is to 1 the better the model is considered, but this does not guarantee the model has good predictive capabilities. BIC takes into account the goodness of fit of the model and applies a penalty for increasing the number of parameters in a model. It approximates the Bayes factor, which is a “summary of the evidence provided by the data in favor of one scientific theory, represented by a statistical model, as opposed to another” (Kass and Raftery 1995, p. 777). In model selection using BIC, the goal is to minimize BIC while making a tradeoff between model fit and number of parameters in the model. The criterion is a measure of the total mean squared error, and an indicator of model bias (Neter et al. 1996). In comparing models, the models with the lowest values are considered those with the least bias.

Given the large number of potential explanatory variables and results of exploratory analysis, the following additional steps were taken to refine the model selection process:

•	Preference was given to models with fewer potential explanatory variables by using BIC as the selection criterion instead of . BIC penalizes additional explanatory variables more than does, although is still useful as a report of model quality.
•	The meteorological stations had varied start dates, depending on the national network and particular sites, so the precipitation data used to build the model were limited to a start date between 1888 and 1899, depending on the subregion, and ending in 1990 (to match the tree-ring chronologies).
•	Precipitation was predicted only for the period 1700–1990; this removed some of the edge effects attributed to the tree-ring chronologies where there were fewer cores in the early parts of their record.
•	Using seasonal 12-year moving averages would result in non-independent, serially correlated residuals. Independent residuals result in the best linear unbiased estimator of the response variable, and allow one to test hypotheses about the regression relation and estimate confidence and prediction intervals (Helsel and Hirsch 1995). Therefore, to reduce correlation in the residuals and leave some data to test the models, the seasonal 12-year moving average series were thinned to every fourth year before developing the precipitation models.
•	The Durbin–Watson test was performed to examine autocorrelation of the residuals. In the Durbin–Watson test (Durbin and Watson 1950, 1951, 1971; Zeileis and Hothorn 2002), the null hypothesis is that the autocorrelation of the residuals is zero and it is tested against the alternative hypothesis that the residuals are autocorrelated. A statistically significant result is undesirable and indicates that the residuals are autocorrelated.
•	The Breusch–Pagan test was performed to examine the variance of the residuals. In the Breusch–Pagan test (Breusch and Pagan 1979; Zeileis and Hothorn 2002), the null hypothesis is that the residuals are independent and normally distributed, with a constant variance. The alternative is that the variance is not constant. A statistically significant result is undesirable and indicates that the residuals are heteroskedastic, violating one of the underlying assumptions of regression necessary to do more than predict the response, such as obtaining the variance for the prediction, obtaining the best linear unbiased estimator, and estimating confidence and prediction intervals (Helsel and Hirsch 1995).

In addition to the quantitative measures, a literature review was done of historical accounts of weather in the study area and of past tree-ring, precipitation and streamflow studies that indicated periods of wet and dry conditions prior to 1900 (Carlyle 1984; Rannie 1998; Sauchyn and Beaudoin 1998; Thorleifson et al. 1998; Red River Basin Board 2000; St. George and Nielsen 2002, 2003; Brooks et al. 2003; Case and MacDonald 2003; St. George and Rannie 2003; Severson and Sieg 2006; Lapp et al. 2013). Severson and Sieg (2006) and Brooks et al. (2003) summarized the work of a number of other studies. In both the US and Canada, the eastern side of the study area is more information rich based largely on the work of Rannie (1998), who focused on the Red River Basin, and Severson and Sieg (2006), who focused on eastern North Dakota. When hindcasting precipitation outside the period of instrumental observation, models were examined in light of these past studies, keeping in mind that wet and dry periods would vary with respect to season, duration and spatial extent.

Results of cluster analysis

There was very little difference in clustering based on mean and median values, and clusters were well defined using precipitation only. Using the mean resulted in more geographically cohesive cluster groups. Therefore, the variable subset chosen for effectiveness and parsimony was mean season 1 total precipitation, mean season 2 total precipitation and mean season 3 total precipitation at each meteorological station with a minimum of 60 observations in each season over its period of record from 1900–2010.

The dendrogram (a complete graphical description of the hierarchical clustering) for the final cluster analysis is shown in Figure 3. In the dendrogram, the vertical lines at the bottom (arranged so that branches of the dendrogram do not cross) represent the individual meteorological stations. Horizontal lines connecting clusters represent the merges of similar clusters. The y-axis represents the distance (in three-dimensional space defined by mean season 1 total precipitation, mean season 2 total precipitation and mean season 3 total precipitation) between the two clusters being merged. The higher on the y-axis that merges occur, the more distant the clusters are from each other in terms of mean seasonal precipitation. The AC of the HACA in this study was high, 0.97, indicating distinct clustering structure.

Figure 3. Dendrogram for cluster analysis with horizontal line intersecting five groups branched below the line.

Determination of the number of cluster groups is subjective in that the analyst determines the number of groups by balancing interest in interpretable groups and avoidance of needless splitting of the data. The number of groups is determined by drawing a line across the dendrogram and examining the main clusters branching out beneath that line. Depictions of three, four and five groups were plotted on maps and examined for reasonable climate similarity, with each group having numerous meteorological stations (with the intent of having no micro-clusters). In addition, boxplots of mean seasonal precipitation for each number of groups were examined for significant difference between groups. By means of a dashed horizontal line, Figure 3 shows the data divided into five groups. The five branches below the dashed line represent the major cluster groups. They are major cluster groups because the linkage distance at which they combine with each other is relatively large, indicating that there are relatively large Euclidean distances between the sites in groups 1, 2, 3, 4 and 5. The clustering with five groups represented well the north–south gradient in winter precipitation and the east–west gradient in precipitation across the study area as shown in the boxplots of seasonal mean precipitation for each group (Figure 4).

Figure 4. Boxplots of mean seasonal total precipitation for each group.

The meteorological stations were mapped in a geographic information system (GIS) with their cluster group assignments, (1) southeast, (2) south central, (3) southwest, (4) northwest and (5) northeast. Voronoi polygons (Voronoi 1908) were drawn around each site used in the cluster analysis (Voronoi polygons are drawn so that every location within a polygon is closer to the meteorological station in that polygon than to any other meteorological station). Then, the polygons were extended as necessary to cover the entire study area and clipped to the study boundary. Each polygon was assigned to the cluster group of the climate site around which it was drawn. The numerous polygons were then dissolved in GIS so there were polygons representing each group. Those sites with shorter periods of record were assigned to the cluster group in which they fell; Figure 5 shows the final five clusters. Part of cluster 4 occurs within cluster 5. A higher density of meteorological stations in this area might result in a continuous cluster for cluster 4; however, the analysis indicated that the outlier station belonged in cluster 4, and it was subsequently grouped with the rest of cluster 4.

Figure 5. The five clusters in the study area. (The additional meteorological stations had shorter or seasonal periods of record and were not used for the cluster analysis.)

The nonparametric Kruskal–Wallis rank sum test (Higgins 2004) was used to test whether all groups have the same mean annual precipitation distribution function or whether at least one of the groups has a different location parameter. For season 1, the Kruskal–Wallis test was statistically significant at the α = 0.01 significance level, which means at least one of the groups has a mean significantly different from the other means. The test was repeated for season 2 and season 3 precipitation and the results were statistically significant. Boxplots (Figure 4) visually confirm the differences.

Results of precipitation modeling

The equations for reconstructing the seasonal 12-year moving average precipitation for each season and cluster group are shown below.(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15)

where Gi represents the ith group and Sj represents the jth season; BB, Boundary Bog tree-ring chronology; TRNP, Theodore Roosevelt National Park tree-ring chronology; BCV, Burning Coal Vein tree-ring chronology; CB, Cedar Butte tree-ring chronology; D4, the 16-year wavelet voice from the associated tree-ring chronology; D5, 32-year wavelet voice; and D6, 64-year wavelet voice.

Model information and quality criteria are shown in Table 1. The models generally show a high coefficient of multiple determination indicating that they explain most of the variability in the 12-year moving average seasonal precipitation. The quality of the models varies with group and season, with better results in the more eastern groups and more difficulty modeling the drier, western areas. Group 3 (the southwestern part of the study area) in particular had low coefficients of multiple determination in seasons 2 and 3. In these seasons the residuals failed the test for serial correlation in the residuals (Durbin–Watson test), as did group 4, season 1. These results can indicate the need for additional explanatory variables, but adding variables did not improve the models (not shown). The existence of serial correlation does not mean the models have no utility. When models have serial correlation in the residuals, the regression coefficients are not the most efficient estimate and the variance may be underestimated (confidence and prediction intervals for the estimated precipitation would be too narrow), but the coefficients are unbiased (Helsel and Hirsch 1995). The patterns provided by the three models with autocorrelation in the residuals are still indicative of long-term drought and pluvial periods. Group 3, season 2, also failed the Breusch–Pagan test for constant variance. Regression models that do not meet the assumptions of independent, normally distributed residuals can still be used to predict the response variable (Helsel and Hirsch 1995) and are still informative for the purposes of this study.

Table 1. Information and associated statistics for precipitation models.

Group	Season	Predictor variables	Calibration period, every fourth year in interval	R^2	R^2a	Durbin-Watson test p-value	Breusch-Pagan test p-value
1	1	BB.D6, TRNP.D5, BCV.D6, CB.D4, CB.D5	1888–1990	0.78	0.74	0.446	0.267
	2	BB.D5, TRNP.D6, CB.D4, CB.D6	1888–1990	0.66	0.62	0.102	0.107
	3	BB.D4, BB.D5, BB.D6, TRNP.D5, BCV.D5, BCV.D6, CB.D4, CB.D5	1888–1990	0.84	0.78	0.045	0.570
2	1	BB.D6, TRNP.D4, TRNP.D5, TRNP.D6, CB.D5	1888–1990	0.71	0.64	0.765	0.949
	2	BB.D5, TRNP.D6, BCV.D6, CB.D4, CB.D6	1888–1990	0.64	0.58	0.016	0.630
	3	BB.D4, BB.D5, BB.D6, TRNP.D6, CB.D5	1888–1990	0.81	0.77	0.172	0.891
3	1	BB.D4, BB.D6, TRNP.D6, BCV.D5, BCV.D6, CB.D6	1890–1990	0.93	0.92	0.053	0.435
	2	BB.D5, TRNP.D6, CB.D6	1890–1990	0.38	0.32	0.001	0.005
	3	BB.D5, CB.D5	1890–1990	0.48	0.45	< 0.001	0.235
4	1	BCV.D4, BCV.D6	1899–1990	0.66	0.63	0.001	0.290
	2	BB.D4, BCV.D5, BCV.D6, CB.D5	1899–1990	0.60	0.54	0.048	0.953
	3	BB.D6, TRNP.D5, BCV.D5, BCV.D6	1899–1990	0.84	0.81	0.014	0.459
5	1	BB.D5, BB.D6, TRNP.D5, BCV.D4, BCV.D6, CB.D5, CB.D6	1888–1990	0.91	0.88	0.268	0.203
	2	BB.D4, TRNP.D5, TRNP.D6, BCV.D6, CB.D5	1888–1990	0.74	0.69	0.419	0.401
	3	BB.D4, BB.D5, BB.D6, TRNP.D6, BCV.D6, CB.D5	1888–1990	0.84	0.81	0.651	0.301

Group 1, southeast cluster of the Souris River Region; group 2, southcentral cluster; group 3, southwest cluster; group 4, northwest cluster; group 5, northeast cluster; season 1, November–February; season 2, March–June; season 3, July–October; BB, boundary bog; TRNP, Theodore Roosevelt National Park; BCV, burning coal vein; CB, Cedar Butte; R², coefficient of multiple determination, R²_a, coefficient of multiple determination adjusted for degrees of freedom.

As an example, the model results for group 1 (southeastern group, mainly the Red River Basin), season 2 (March-June) 12-year moving average precipitation (Equation 2) are shown in Figure 6. Using this model, modeled precipitation matches well with observed precipitation and corresponds well with known wet and dry periods. The early 1700s (approximately 1703–1720) had periods of drought or dry years in the northern Great Plains and central North Dakota (Severson and Sieg 2006). Lapp et al. (2013) documented 24 sustained droughts over the past 600 years in the northwestern Great Plains and found droughts in 1701–1708 and 1717–1721, which matches well with the severe drought shown in Figure 6. Of the droughts they found, the most intense was the drought of 1717–1721, which also supports Figure 6 showing the early 1700s drought being worse than the 1930s drought. St. George and Nielsen (2003) documented floods on the upper Red River in 1726, 1727 and 1741; 1753–1762 was quite dry in parts of the northern Great Plains (Severson and Sieg 2006), and from the 1820s to about 1861, conditions seem to have been quite wet. There were five successive high or very high runoff years on the Red River from 1823 to 1828 and five successive high or very high runoff years again on the Red River in 1847–1852 (Thorleifson et al. 1998). The year 1849 stands out in records documented by Severson and Sieg (2006) as a very wet year, whereas 1852 saw one of the largest floods on the Red River at Winnipeg and extreme flooding on the Assiniboine River. The Red River again experienced a large flood at Winnipeg in 1861 (Red River Basin Board 2000). St. George and Nielsen (2002) described a pronounced wet interval in the 1850s in southern Manitoba.

Figure 6. Group 1 (southeast), season 2 (March–June) modeled and observed 12-year moving average precipitation.

While predicted 12-year moving average precipitation stayed above the mean from the 1820s into the 1850s, there was a decline in the middle of this period (Figure 6). Case and MacDonald (2003, p. 713) report that

There is ample historical documentation of meteorological drought during the mid-1800s across the northern Great Plains (e.g. Mock 1991; Blair and Rannie 1994). Tree ring reconstructions of precipitation have also indicated drought during the mid-19^th Century in the southern Canadian Prairies (Sauchyn and Beaudoin 1998), Rocky Mountain foothills (Case and MacDonald 1995), and Montane regions (Watson and Luckman 2001).

This drought may have caused the dip in precipitation just before 1850 and the sudden shift from wet to dry conditions after 1850 (Figure 6).

Dry, drought and/or low streamflow conditions were reported from about 1861 to about 1880 in the Red, Wild Rice (North Dakota), Sheyenne, Souris and Missouri River Basins; the Great Plains, southern North Dakota; eastern North Dakota; and central North Dakota (Thorleifson et. al 1998; St. George and Nielsen 2003; Severson and Sieg 2006). In the early 1880s, there was a shift from dry to wet, and historical photographs (Ryberg 2015) and instrumental records (US Geological Survey 2014) show floods in the Missouri, Red and Souris river basins. From 1900 on, the modeled and observed precipitation match well (Figure 6). The lowest recorded annual peak streamflow on the Souris River above Minot, North Dakota, occurred in 1931 during the widespread 1930s drought, and this model (for an area to the southeast of the Souris River Basin) estimated that dry period. The larger observed precipitation values at the end of the record are of approximately the same magnitude as those modeled over 150 years ago around 1850.

Season 3 (July–October) 12-year moving average precipitation (Equation 15) is shown for group 5 (northeast; lower Red River Basin, southern Manitoba and southeastern Saskatchewan) in Figure 7. The early part of the modeled precipitation shows a drought centered around 1750. St. George and Nielsen (2002) documented a severe drought in 1753 in southern Manitoba, with annual precipitation estimated at more than two standard deviations below the mean. There appears to have been a widespread regional drought in the 1750s to the early 1760s. Meko (1982) found a period of drought or dry years in the western Great Plains from 1753 to 1762 based on Ponderosa pine in North Dakota, South Dakota, Nebraska, Wyoming and Montana. Lapp et al. (2013) found a sustained drought in the northwestern Great Plains from 1755 to 1761, and Stockton and Meko (1983) found a major historical drought in the mid to late 1750s in the Great Plains. There is another severe drought around 1800 (Figure 7). Thorleifson et al. (1998, p. 191) described the period of 1793–1828 as one of “high variability” with numerous floods and “several drought episodes.” Rannie (1998) indicated that there was low runoff in the Assiniboine, Red and Clearwater Rivers in 1800. In group 5 (Figure 7), farther north than group 1 (Figure 6) and in a different season, the pronounced wet period during the 1800s started earlier, in this case in the 1820s. This suggests that the generally wetter period in the first half of the 1800s may have started with wetter falls. Wet falls provide antecedent conditions that contribute to large floods in the region (Ryberg et al. 2007) and may have contributed to the large floods documented on the Red and Assiniboine in the 1800s (Harrison and Bluemle 1980; Miller and Frink 1984; Rannie 1998; Severson and Sieg 2006). The modeled and observed precipitation match well, and the wet period of the end of the record is able to be modeled.

Figure 7. Group 5 (southern Manitoba and southeastern Saskatchewan), season 3 (July–October) modeled and observed 12-year moving average precipitation.

For season 1 (November–February), good historical or analytical accounts of snow amounts are difficult to obtain because blizzards can occur with highly variable snow amounts and season 1 is not a growing season; therefore, no graph is shown. The 1820s did see numerous large floods on the Red River (most likely driven by snowmelt) and there were accounts of snowstorms in this decade. There were reports of deep snowpack near the Red River Settlement and throughout the southern Red River Basin in the winter of 1825–1826 (St. George and Rannie 2003). On 20 December 1826, a snowstorm described as “fearful” drove away the bison from the Pembina region, and resulted in the loss of many horses and 33 people (Severson and Sieg 2006, p. 35).

As a check of the models developed for the five groups, modeled precipitation from the three seasons (November–February, March–June, July–October) was added to get the total 12-year moving average precipitation for each group. Despite being modeled separately, the seasonal precipitation summed to reasonable annual totals. The seasonal totals for groups 1 and 3 are shown in Figure 8. The summed seasonal precipitation matched annual totals well in groups 1 and 5 and did less well in groups 3 and 4 where it matched general patterns but did not match the highs and lows, and this is evidenced in the lower coefficients of multiple determination for groups 3 and 4 (Table 1).

Figure 8. Twelve-year moving average annual precipitation for group 1 (southeast) and group 3 (southwest).

Conclusions

Using the multiresolution decomposition of the tree-ring chronologies in the study area shows that the multidecadal precipitation variability can be predicted using tree-ring variability on 16-, 32- and 64-year time scales. The time scales vary with location in the study area and with season. The most frequently used tree-ring site was Boundary Bog (indicated by BB in the equations and in Table 1). In the Great Plains, cold, dry air masses from central Canada interact with warm, moist air from the Gulf Mexico to set the stage for the outbreaks of severe thunderstorms that can bring heavy rainfall to localized areas (Wishart 2004). While regression analysis is based on correlation and is not a causal analysis, use of the Boundary Bog site may capture patterns in the cold, dry air masses moving south in Canada. The most frequently used wavelet voice over all sites was the 64-year wavelet (indicated by D6 in the equations and in Table 1). While tree rings are more sensitive to drought conditions than to wet conditions, the cluster groups within the greater Souris River Basin that were best modeled with the multiresolution decomposition of the tree-ring chronologies were those in the wettest part of the study area, groups 1 (southeast group, mainly the Red River Basin in the US) and 5 (northeast group, southern Manitoba and southeastern Saskatchewan).

The modeled precipitation with the three seasons summed generally compares well with the observed data, showing that the individual seasonal models are good at predicting precipitation during the period of instrumental observation. The models also generally match low precipitation (high precipitation) with known drought (pluvial) periods. The model for group 3, the southwest, did the least well at matching observed precipitation (Figure 8). This group had some of the poorest seasonal models (Table 1). However, the resulting modeled precipitation pattern seems reasonable, although it misses the wet period in the 1940s.

One issue is that the tree-ring chronologies did not have tree rings recent enough to model the wet period at the end of the precipitation period of record. The tree-ring record ends at approximately the same time as an abrupt change from dry to wet conditions in the early 1990s. However, the recent wet period seems to have started earlier in season 3, particularly for the eastern side, groups 1 and 5. Season 3 for group 5 is an example in which the high precipitation in recent times started earlier in season 3 and was able to be modeled using the decomposed chronologies (Figure 7). The model also shows very high precipitation in the first half of the 1800s, and historical accounts corroborate this. The recent wet period may be similar to that of the 1800s and may be part of natural variability on a very long time scale, and in many cases there probably are not enough instrumental data to detect it.

The modeled precipitation and the accounts of past pluvial and drought periods are characterized by sudden shifts from wet to dry and dry to wet. These sudden shifts have been documented elsewhere in the region and appear to be characteristic of climate in the northern Great Plains (Vance et al 1992; Shapley et al. 2005; Vecchia 2008). By describing regional long-term variability in seasonal precipitation, these results will be used to inform models of future Souris River streamflow, thereby helping the International Souris River Board and the Souris River Flood Task Force evaluate future water-resource management options.

Acknowledgements

This project was funded by the North Dakota State Water Commission and the US Geological Survey Cooperative Water Program. The Souris River Basin boundary was provided as a GIS shapefile by Tara Gross (US Geological Survey, written communication, 11 July 2014). Thanks to Kevin Vining of the US Geological Survey North Dakota Water Science Center and Gregory Pederson of the US Geological Survey Northern Rocky Mountain Science Center for providing valuable early reviews, and to two anonymous reviewers for their constructive comments. Thanks as well to Valerie Kling for translating the abstract into French.