Online Resource 2: Inflow projections

Choosing Rainfall Scenarios

Analysis of precipitation and temperature gradients across the Prairie Pothole region (PPR) of North America suggest that the climate has been getting warmer and wetter during the 20th century(Millett et al., 2009). Across the region, average annual precipitation increased by 9%, however this increase was greatest in the eastern part of the PPR. The eastern PPR has become wetter while the western PPR has been drying (Millett et al., 2009). The results of previous analyses extend over an area of ~750 000km2 and use data from weather stations that are situated at a large distance from the refuges of interest in this study. Because the refuges in our study are roughly mid-way between the eastern and western extents of the PPR defined by Millett, it is unclear what rainfall trend might be expected from the Millett analysis alone.

To investigate likely precipitation trends in the region near the refuges, we analyzed precipitation data from the nearest weather station to each refuge. Monthly precipitation data was obtained from the US Historical Climatology Network (Menne et al., 2011). At each of the 5 stations, the data record spanned the years 1895-2010.

Methods/Models fitted

The 5 data stations were analyzed using a mixed model to distinguish a fixed trend from the random effects (Pinheiro and Bates, 2000). The mixed model framework enables synthesis of all 5 weather stations to separate the trend from the random fluctuation in each model. In each model we allowed for fixed effects due to year and site (b) and a random effect for site-to-site variability (bi). These factors were the same in each model that we considered. Different candidate models arise by considering whether or not the precipitation can be predicted from recent years—i.e., if the rainfall for the current year is similar to previous years or if it changes randomly. We had reason to suspect that rainfall is autocorrelated from personal communication with refuge staff, who noted that the refuges tended to go through sustained periods of drought and flood, suggesting correlation between years.

We tested models with a range of autocorrelation structures to find the most parsimonious structure. The candidate models differ by how many years of data affect the prediction of the current precipitation. The linear model assumes that residuals are normally distributed around the fixed trend, assuming no relation between the residuals between years. The linear AR(n) models predict the new precipitation using the last n years of residual data. We were unable to investigate lag times longer than 10 years due to computational limitations. The formulae for each of the models, the specification of their error structures, and the criteria for choosing the best model are contained in Table 1, Table 2 and Table 3.

Table 1: Description of candidate models tested in the analysis

Model / Formula / Error
Linear / /
Linear with AR(1) autocorrelation / /

Linear with AR(n) autocorrelation / /

Table 2: Goodness of fit statistics for different models 1895-2010. The linear AR(n) models predict the new precipitation using the last n years of residual data.

Model / AIC / Log-Likelihood / p-value
Linear / 3483.19 / -1737.59 / -
Linear AR(1) / 3482.83 / -1736.42 / 0.125
Linear AR(2) / 3474.55 / -1731.27 / 0.0018*
Linear AR(3) / 3455.54 / -1720.77 / <0.0001*
Linear AR(5) / 3452.94 / -1717.47 / <0.0001*
Linear AR(10) / 3439.10 / -1705.55 / <0.0001*

From the likelihood ratio test, all linear AR(n) models other than n=1 were significantly different (p<0.05) to the linear model. The maximum likelihood ratio was between the linear model and linear AR(10), with L.ratio of 64.08.Model AR(10) was also significantly different to model AR(5), with likelihood ratio 23.84 and p-value 0.0002.In fact all models were significantly different to each other except for the linear and linear AR(1). This suggests that the addition of the autocorrelation terms explains a significant amount of additional variance, so the correct model to select for the period 1895-2010 is the linear AR(10) model.

Table 3: Goodness of fit statistics for different models 1960-2010. The linear AR(n) models predict the new precipitation using the last n years of residual data

Model / AIC / Log-Likelihood / p-value
Linear / 1488.73 / -740.37 / -
Linear AR(1) / 1489.71 / -739.85 / 0.3113
Linear AR(2) / 1491.01 / -739.52 / 0.4284
Linear AR(5) / 1493.68 / -737.84 / 0.4088
Linear AR(10) / 1491.62 / -731.81 / 0.0720

From the likelihood ratio test, none of the linear AR(n) models were significantly different (p-value<0.05) to the linear model (p=0. 072 for AR(10)). The maximum likelihood ratio was between the linear model and the linear AR(10) model with a likelihood ratio of 17.11. Adding the autocorrelation terms did not explain a significant amount of additional variance, so the most parsimonious model is the linear model for the period 1960-2010.

Predictions from 1895-2010 record

For the period 1895-2010, the most parsimonious model was an autoregressive AR(10) model (Figure 1), suggesting that there were long phases of correlated precipitation events (e.g., sustained periods of flood or drought). Estimated model parameters are included in Table 4. We were unable to investigate lag times longer than 10 years due to computational limitations, but the data suggest that there were correlations that lasted at least 10 years in the full rainfall record. There was no significant linear trend in this model (i.e., year was not a significant predictor of precipitation; p=0.627), meaning that the precipitation can be explained by a location variable (random intercept, p<0.00001) and the previous 10 years of precipitation data.

Table 4: Parameter estimates for the Linear AR(10) mixed model fit to the precipitation data 1895-2010.

Parameter / Estimate / SE / p-value
Fixed Effects
Intercept / 20.469 / 1.2264 / <0.00001
Year / -6.476x10-3 / 0.0133 / 0.6272
Random Effects
Intercept / 1.803 / 4.9216 (residual) / -
Autocorrelation parameters
(ρ1, ρ2,…, ρ10) / 0.0378
0.1014
0.1832
0.0865
0.0154
-0.0031
0.1262
-0.0429
0.1502
0.0070 / - / -

From personal communication with refuge staff, we knew that the area went through a significant dry period during the 1930s, and further dry periods in the late 1950s and early 1960s, before becoming increasingly wet in the last 50 years. This trend is evident by examining the top left plot of Figure 1, where we can see a general decrease in the early years of the century (c1900-1940), followed by an increase since the 1960s. This is not captured by the linear trend as the decrease and increase tend to cancel each other out so that there is no evidence of a net trend over the whole century. Because climate effects are expected to be most evident in more recent data, we also modeled the past 50 years of data to determine the rate of the recent wet period.

Figure 1: Linear mixed model fits to the 5 rainfall gages closest to the 5 study refuges for the period 1895-2010. Data are annual precipitation in inches. For these data the most parsimonious model was a linear model with AR(10) autocorrelation. This means that the past 10 years of data influence the precipitation in the current year, or that there are rainfall ‘cycles’ of at least 10 years in this data set. The top-left figure shows the data from all 5 stations and the mixed model mean fit to the data. The other figures show the data from one station as well as the mixed model mean fit (solid line), and a realization of the data using the random effects parameters for this station (dashed line). There has been no significant linear change in mean annual precipitation at these stations over the measurement period.

Predictions from 1960-2010 record

For the period 1960-2010, the most parsimonious model was the linear model with no autocorrelation (Figure 2). Random fluctuations about the mean trend line are explained with a normal random model and do not depend on past values. The trend lines in the figure are significant (p=0.0001). Estimated model parameters are included in Table 5. The period was characterized by a trend of increasing precipitation at a rate of 0.07 inches/year (p=0.0001). From an initial average precipitation of 17.79 inches in 1960 (p<0.0001), the modeled annual precipitation increased to 21.36 inches in 2010, an increase of 20% over 50 years (0.4% per year).

Table 5:Parameter estimates for the best-fit Linear mixed model fit to the precipitation data 1960-2010.

Parameter / Estimate / SE / p-value
Fixed Effects
Intercept / 17.794 / 0.7489 / <0.0001
Year / 0.0714 / 0.0184 / 0.0001
Random Effects
Intercept / 1.1767 / 4.3167(residual) / -

Summary of predictions

The most parsimonious mixed model fits for the periods 1895-2010 (no significant trend) and 1960-2010 (trend increasing at 0.4%/year) are plotted in figure 2.

Figure 2: Mixed model fits to the annual precipitation data for four National Wildlife Refuges (Arrowwood (North Dakota), Lacreek (South Dakota), Sand Lake (South Dakota) and Tewaukon (North Dakota)) for 1895-2010. (a) AR(10) (Autoregressive, 10 step) linear fit to all data 1895-2010; (b) linear fit to 1960-2010. There was no significant linear trend in mean annual precipitation, 1895-2010 (p = 0.62), but precipitation increased by 0.4%/year, 1960-2010 (p < 0.0001).

Variability Predictions

As well as changes in the mean, the variability of precipitation may change over time. The mixed model framework assumes a stationary variance around the trend, so this modeling approach cannot be used to determine if the variability is changing over time. To investigate whether variability has changed in the rainfall record we placed the rainfall data for each decade into interval bins and found the standard deviation of each bin. The resulting standard deviations give a measure of the variability of the rainfall during the relevant decade for each station (Figure 3). If there is a systematic trend across the region then the best fit linear regression line should be significant and the same trend should appear in at least the majority of the stations. We found the regression line was a poor fit to the data, with low R2 values and high p-values at each station. There does not appear to be a strong trend in variability over the rainfall record. Even if we restrict the data to the wetting period from 1960-2010, we do not obtain a significant trend in the variability (Figure 4).

Figure 3: Linear mixed model fits to the 5 rainfall gages closest to the 5 study refuges for the period 1960-2010. Data are annual precipitation in inches. For these data the most parsimonious model was a linear model with no autocorrelation. The top-left figure shows the data from all 5 stations and the mixed model mean fit to the data. The other figures show the data from one station as well as the mixed model mean fit (solid line), and a realization of the data using the random effects parameters for this station (dashed line). There has been an average mean increase in precipitation of 0.4% per year at these stations over the past 50 years (this corresponds to a total increase of 19% since 1960).

Figure 4: Standard deviation of annual precipitation for each decade from 1900-2010 at each of the 5 rainfall stations. Solid lines give the best fit linear regression line to the data. None of the trend lines were significant at the alpha=0.95 level, suggesting no change in the variability of the rainfall over the period of the precipitation record.

Figure 5:Standard deviation of precipitation for each decade from 1960-2010 at each of the 5 rainfall stations. Solid lines give the best fit linear regression line to the data. None of the trend lines are significant at the alpha=0.95 level, suggesting no change in the variability of the rainfall over the period of the precipitation record.

Choosing inflow scenarios for the adaptive management modeling

Modeling the rainfall recorded at the gage nearest each of the study refuges provides an indication of how precipitation entering the system has changed over the past century. From our modeling, there has been no net trend in precipitation over the past 115 years. However there are temporal lags in the data that represent wet and dry periods in the rainfall record. These correlated periods are at least 10 years long. Anecdotal evidence of the rainfall history of the area suggests that the precipitation decreased in the early part of the 20th century before increasing again after the early 1960s, and this is confirmed by inspecting the rainfall record.

There are insufficiently long inflow records to carry out a similar analysis on inflow data at the refuges, however we postulate that the general trends in rainfall will carry through the system and be reflected in the inflow to the refuges. This assumption ignores the effect of any changes to the albedo or evaporation as a result of changing temperatures, as well as changes to the groundwater systems in the area. Despite the apparent limitations of these assumptions, we note that we are not seeking an exact figure for inflows to the refuges, but instead trying to propose inflow scenarios that capture the magnitude of the changes we might see in the future.

From the precipitation data, it is unclear what the future holds for this area. Although the precipitation has been increasing steadily for the past 50 years, we have evidence that the system evolves slowly on a long time scale. We also know that there was a sustained dry period that immediately preceded the current wetting period. We cannot say with certainty if the current wetting period is indicative of a regime shift that will continue, or if the system will return to a drying period in the near future.

What we do know about the precipitation is the rate that it is currently increasing. Although the average precipitation has not changed significantly over the past 115 years, the rate of change over the past 50 years was 0.4% per year. The trend could continue at the current rate but it could also return to its long term average by becoming negative. We propose two possible inflow scenarios to capture this uncertainty—one in which the current trend of +0.4%/year continues, and one in which the direction of change reverses and becomes a drying trend, so that we see a change of -0.4% per year. Because the past rainfall record suggests that the variability is not changing systematically over time we leave the standard deviation of the inflow the same as the present value in both scenarios.

In both of the inflow scenarios, we assume that the system has no memory (i.e., no autocorrelation). This is consistent with the best fit model for 1960-2010. It is also necessary for the Markov model that we will use to specify the adaptive management plan for each refuge. This means that our projections are unlikely to hold for long periods of time, as the rainfall record suggests that over long time periods that autocorrelation is important in the system. However making predictions into the distant future is always error-prone, and we argue that this approach is reasonable in the short-term. The best approach will be to manage for the short-term future and monitor how the climate evolves, and then revise the models as necessary. Our model projects inflows for 20 years into the future based on the linearly increasing trend observed since 1960.

References

Menne M, Williams C, Vose R (2011) United States Historical Climatology Network. Accessed: 2nd November 2011.

Millett B, Johnson WC, Guntenspergen G (2009) Climate trends of the North American prairie pothole region 1906-2000. Climatic Change 93:243-267.

Pinheiro J, Bates D (2000) Mixed-Effects Models in S and S-PLUS. Springer-Verlag, New York, NY.