IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS IN ROCKY MOUNTAIN NATIONAL PARK, COLORADO, U.S.A. GUIMING WANG 1⋆ , N. THOMPSON HOBBS 1 , FRANCIS J. SINGER 1 , DENNIS S. OJIMA 1 and BRUCE C. LUBOW 2 1 Natural Resource Ecology Laboratory, Colorado State University, Fort Collins, CO 80523, U.S.A. 2 Department of Fisheries and Wildlife, Colorado State University, Fort Collins, CO 80523, U.S.A. Abstract. Changing climate may impact wildlife populations in national parks and conservation areas. We used logistic and non-linear matrix population models and 35 years of historic weather and population data to investigate the effects of climate on the population dynamics of elk in Rocky Mountain National Park (RMNP), Colorado, U.S.A. We then used climate scenarios derived from Hadley and Canadian Climate Center (CCC) global climate models to project the potential impact of future climate on the elk population. All models revealed density-dependent effects of population size on growth rates. The best approximating logistic population model suggested that high levels of summer precipitation accelerated elk population growth, but higher summer minimum temperatures slowed growth. The best approximating non-linear matrix model indicated that high mean winter minimum temperatures enhanced recruitment of juveniles, while high summer precipitation enhanced the survival of calves. Warmer winters and wetter summers predicted by the Hadley Model could increase the equilibrium population size of elk by about 100%. Warmer winters and drier summers predicted by the CCC Model could raise the equilibrium population size of elk by about 50%. Managers of national parks have relied on effects of weather, particularly severe winters, to regulate populations of native ungulates and prevent harmful effects of overabundance. Our results suggest that these regulating effects of severe winter weather may weaken if climate changes occur as those that are widely predicted in most climate change scenarios. 1. Introduction Understanding the consequences of shifts in global climate for the operation of biotic processes has emerged as a pressing challenge confronting contemporary ecologists (Watson et al., 1996; Shaver et al., 2000). The effects of climate are often manifested as part of a composite of forces that may interact to shape ecosystem dynamics (NcNaughton, 1983). It follows that understanding the effects of climate requires understanding the interplay between climate and other agents of change (Shaver et al., 2000). Populations of large herbivores can exert strong effects on plant communities (McNaughton, 1985), and these effects can amplify or attenuate a variety of ecosystem processes (Hobbs, 1996). As a result, the impacts of a changing climate on herbivores could extend to many other components of natural and human ⋆ Address for correspondence: Department of Biological Sciences, Arkansas Tech. University, Russellville, Arkansas 72801, U.S.A., E-mail: guiming.wang@mail.atu.edu Climatic Change 54: 205–223, 2002. © 2002 Kluwer Academic Publishers. Printed in the Netherlands. 206 GUIMING WANG ET AL. dominated systems. In temperate environments, growth rates of populations of large herbivores are often limited by effects of winter weather on survival and recruitment (Mech et al., 1987; Post and Stenseth, 1999). Energy costs of thermoregulation and locomotion can increase dramatically during winter (Hobbs, 1989), and these elevated requirements must be offset by energy gained from food supplies that seasonally shrink in quantity and diminish in nutritional value. Imbalance between energy supplies and demands elevates starvation mortality during winter (Bartmann and Bowden, 1984; Hobbs, 1989). Given that severe winter weather imposes a strong limit on herbivore populations, climate warming may cause problems of overabundance, particularly if warming predominates during winter, as is widely predicted (National Assessment Synthesis Team, 2001; Dickinson, 1986). These effects may be particularly important at higher latitudes (Houghton et al., 1995, 1996). Excessive abundance could become particularly problematic in national parks and other protected areas. Many of these areas lack large capable predators. Culling by humans is proscribed by law or prevented by cultural preference (Berger, 1991). Thus, because many herbivore populations are not limited by the top-down effects of predation, the expression of bottom-up feedbacks from plants may result in enduring impacts on plant communities (Noy-Meir, 1975; Gill, 1992; Stromayer and Warren, 1997). This is simply because in some plant communities, negative feedback to herbivore population growth operates through a reduction in the abundance and distribution of palatable plants. However, these same plants can contribute in important ways to biological diversity and other ecosystem values (Berger, 1991; Garrott et al., 1993; Decalesta, 1994). Perceived impacts of elk (Cervus elaphus) populations on vegetation of national parks in the Rocky Mountains have stimulated a lively debate about how those populations should be managed (Cole, 1971; Houston, 1971; Chase, 1986; Hess, 1993; Boyce, 1998). In 1968, RMNP adopted a natural-regulation policy for management of ungulates, which allowed populations to fluctuate without human intervention within the Park (Boyce, 1991; Singer et al., 1998). As part of this policy, culling of elk in RMNP ended, and the population increased steadily thereafter (Lubow et al., 2002). Here, we report studies of potential effects of climate warming on dynamics of an increasing elk population in RMNP. In particular, we used model selection techniques to separate the role of climate from density dependent mechanisms of population regulation. We examine how plausible scenarios for future climate might amplify problems of overabundance. IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 207 2. Methods 2.1. MODELING APPROACH Weather in the Rocky Mountains displays marked interannual variability. We sought to use observed responses of the elk population to historic variability in weather to understand how changes in future weather might affect their population dynamics. Our overall modeling approach was to select the best approximating models of population growth incorporating the effects of past weather and population densities. We then used scenarios describing likely future climates based on the Hadley and CCC general circulation models to drive our population models. In the subsequent section, we describe the development of alternative models. We then describe our model selection procedures. Finally, we discuss the development of future climate scenarios and their application in population model projections. 2.2. MODEL DEVELOPMENT We constructed two models and fit them to the observed data on elk population sizes using weather factors as independent variables. The first model was a discrete logistic model, Nt +1 = Nt + Nt Rt (1 − Nt /K) , (1) where Nt +1 was the size of elk population at year t + 1, Nt the size of population at year t, K the equilibrium population size of elk, and Rt annual population growth rate. Rt was a function of weather variables, Rt = R0 + ai Wit , (2) where R0 was the intercept, Wit the weather covariate at year t, ai the coefficient for each weather covariate, and i = 1, 2, . . . , 12. R0 represented the population growth rate without the influences of the weather or when Wit = 0, but not the intrinsic population growth rate. R0 , K, and ai were unknown parameters and were estimated with a nonlinear optimization algorithm, described below. The second model was a nonlinear matrix model with five sex-age/stage groups: calf (NJ ), male yearlings (NYM ), female yearlings (NYF ), bull (NAM , >2 years), and cow (NAF , >2 years) (Figure 1). We assumed the survey occurred at mid winter and that young were born in a single pulse in June of each year. We defined recruitment as the proportion of adult females that produce young that subsequently live to the census, and calf survival as the proportion of young of the year that survived from their first census to their first birthday (Figure 1). We assumed that the litter size was one calf for the majority of reproducing cow. We modeled elk population size (N) with equations, NJ t +1 = NAF t Rr , (3) 208 GUIMING WANG ET AL. Figure 1. Schematic representation of the two-sex nonlinear matrix population model. Symbol mf is the fecundity rate of female, f 2 the survival rate of calves to yearling, s2 the survival rate of yearling to adults, s3 adult survival rate, and t, t + 1, t + 3 years. NYF t +1 = NJ t SJ rc , (4) NYM t +1 = NYM t SJ (1 − rc) , (5) NAM t +1 = 0.93 NAM t + NYM t SYM , (6) NAF t +1 = 0.93 NAF t + 0.93 NY F t , Nt +1 = NJ t +1 + NYF t +1 + NYM t +1 + NAM t +1 + NAF t +1 , (7) (8) where Rr was the recruitment, SJ calf survival rate, rc the proportion of females in calves, and SYM male yearling survival rate. Rr , SJ , rc, and SYM were unknown IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 209 parameters to be estimated. Calf recruitment rate, Rr , and the calf survival rate (SJ ) were modeled as p in the equation: logit(p) = a0 + a1 Nt + ai Wit , i = 2, 3, . . . , 13 , (9) where a0 was the intercept, p either survival or recruitment rate, and a1 , . . . , a13 coefficients; therefore, p = exp(a0 + a1 Nt + ai Wit )/(1 + exp(a0 + a1 Nt + ai Wit )) . We assumed that rc and SYM were constant over years and that male and female adults and female yearlings had the same and constant survival rate of 0.93 (Lubow et al., 2002). We used 12 weather variables in our discrete logistic and matrix models. The weather variables included the winter mean monthly minimum temperature, winter mean monthly maximum temperature, winter mean monthly temperature, winter mean monthly precipitation, winter maximum temperature, winter minimum temperature, summer mean monthly minimum temperature, summer mean monthly maximum temperature, summer mean monthly temperature, summer mean monthly precipitation, the number of days minimum temperature was below 0 ◦ C, and the number of days minimum temperature was below –17.8 ◦ C (0 ◦ F). The summer months were from April to September and winter months from October to March. We used data on the elk population size and composition of the RMNP elk herd from 1965 to 1999 to estimate model parameters and to select the best approximating models. This elk herd winters in the eastern portion of RMNP and is a demographically distinct elk herd (Lubow et al., 2002). Animals spend the summer at high elevation. A portion of these animals return seasonally to low elevation ranges within the Park and remain there throughout the winter (Larkins et al., 1995). Survey data on the wintering elk population were obtained as follows. Animals were counted from 1965 to 1999, either on the ground or from the air, using both a helicopter and a fixed-wing airplane. Data on the composition of the elk herd were collected from 1965 to 1978 and from 1984 to 1999. Surveys of elk counts and composition of the herd were conducted from 1965 to 1991 by National Park Service employees. Elk counts were adjusted for sightability following the procedures of Samuel et al. (1987) and Steinhorst and Samuel (1989). The sightability method takes into account elk groups present in the surveyed area, but not sighted because of different sighting characteristics, e.g., group size, group activity, tree and snow cover (Lubow et al., 2002). Data on the daily and monthly mean, minimum and maximum air temperatures, and precipitation were collected at the Estes Park Weather Station, National Weather Data Service. The number of days that minimum temperature was below 0 ◦ C or –17.8 ◦ C (or 0 ◦ F) was counted using daily temperature data. Although data 210 GUIMING WANG ET AL. on snow accumulation would have been valuable to our studies, these data were not consistently available and were not used in the analysis. 2.3. MODEL SELECTION We used information theoretic methods (Buckland et al., 1997; Burnham and Anderson, 1998) to quantify strength of evidence for alternative models and to estimate their parameters. The relative support for models in data can be assessed using likelihood theory and the Kullback–Leibler information discrepancy, I (f, g), which measures the information that is lost when model g is used to approximate truth f (Burnham and Anderson, 1998). Thus, I (f, g) estimates the relative ‘distance’ between a given candidate model and an unknown truth. This measure has a deep basis in information theory (reviewed by Burnham and Anderson, 1998), but rather than develop that basis here, we simply note that the Kullback–Leibler discrepancy, combined with maximum likelihood estimation (MLE), provides Akaike’s Information Criterion (AIC). AIC estimates the expected relative Kullback–Leibler distance. Thus, a single AIC value has no interpretation by itself, but comparing AIC values among alternatives allows us to assess the relative support in the data for two or more models. Using AIC allows us to select the model (from a set of alternatives) that achieves the best possible compromise between accuracy and parsimony and that provides the optimum approximation of a complex truth. This selection process proceeds as follows (Burnham and Anderson, 1998). Presume we have a set of R (parametric) models, denoted here simply as g1, . . . , gR. We also have a data set that can be used to estimate the parameters of the R models. We will refer to the set of parameters defining each model as θr . The likelihood function L(θr |data, gr)r is maximized, and the maximized likelihood, L(θ̂r ), is at the MLE of θr . Given L(θ̂r ), we calculate AIC as: AIC = −2 log(L(θ̂)) + 2K , (10) where K is the number of parameters in the model. When sample sizes are small, as they were in our study, AIC is adjusted as: AICc = −2 log(L(θ̂)) + 2K(n/(n − K − 1) , where n is sample size or the number of observations, K is the number of parameters, and L(θ̂) is the maximum value of the likelihood function. We are interested in the relative support in the data for alternative models. We assess that support using the difference between AIC values: r = AICr − min(AIC) , (11) where min(AICr ) is the model with the minimum AIC value. This model corresponds to the single best Kullback–Leibler model selected from the alternatives given the data at hand. IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 211 Because model parameters are estimated based on data, there is some uncertainty about which is the “best” model. This uncertainty can be estimated with Akaike weights, wr . Akaike weights are based on the likelihood of a model given the data as: L(θ̂r |data) = e
− 21 r  , (12) We can normalize these model likelihoods so that they are relative weights that sum to one: wr = e
− 21 r 
 R  − 12 r (13) . e i=1 The relative likelihood of model r versus model j is wr /wj . The wr may be thought of as probabilities that the estimated model r is the best approximating model for the data at hand, given the set of models considered. Thus, the model selection uncertainty (given the set of R models) can be compared from these Akaike weights. We fit models (1) and (2) to the estimated elk population sizes. For our matrix model, we fit Equations (3) to (8) to observations on population sizes, sex and age composition of the herd. The unknown parameters of the two types of models were estimated by maximizing likelihood. Assuming normally distributed errors, the likelihood function (loglf ) was computed as: loglf = − i j [(Xij − Yij )2 /SE2ij ] , (14) where Xij was the j th model-estimated value of the ith quantity, Yij the j th observed value of the ith quantity, and SEij the standard error at time j associated with the observed value of quantity i. In order to use information from the survey data as much as we could and obtain precise estimates of the unknown parameters, we included the observed quantities of population size and age-sex composition ratios in the likelihood function (14) for our matrix model (White and Lubow, 2002). We only used the estimated population sizes and associated variances in the likelihood function for the discrete logistic model. We programmed the above model Equations (1) to (9) using Visual Basic for Application (VBA, Sanna, 1997). Equations (1) to (2) and (3) to (9) were embedded in the Newton nonlinear optimization algorithm (Press et al., 1992, p. 362). To estimate the unknown parameters, the value of objective function (14) was maximized with the Premium Sovler Plus 3.5 (Frontline Systems, 1999) add-in in Excel. We used a forward selection procedure to examine effects of weather variables. We first fit the single-weather factor model to the observed data and obtained the value of AICc for each single-factor model. We selected the model of the smallest AICc as the best single-factor model (BSM), and then added one of the remaining 212 GUIMING WANG ET AL. candidate weather factors to the BSM to build the two-factor model. If the AICc value of the two-factor model was not improved, the model selection was finished; otherwise, we obtained the best two-factor model (BTM). We iteratively added one of the remaining candidate factors to the selected best model of the previous step to build the three-factor model, four-factor model, and so on, until the AICc value was not improved by adding any remaining candidate factors. In each step, each previously added-in factor was removed to test if the removal improved or reduced the AICc value. If the reduction occurred, the factor was removed from the selected best model at that step. For the matrix model, we built three different models for each weather factor at each step. Model one had a weather factor in both calf survival and recruitment regression equations; model two had the weather factor only in the calf survival rate equation; and model three had the factor only in the calf recruitment equation. We computed the AICc values for the three models, and selected the model of the lowest AICc value as the selected best matrix model for that selection step. 2.4. DEVELOPMENT OF CLIMATE SCENARIOS AND MODEL PROJECTION We developed future climate scenarios from data and climate projections developed for the U.S. National Assessment: The Potential Consequences of Climate Variability and Change (http://www.nacc.usgcrp.gov/National Climate) by the VEMAP Phase 2 Data Development Group at the National Center for Atmospheric Research (VEMAP Website: http://www.cgd.ucar.edu/vemap/ve298.html). We used the two general circulation model scenarios used in the U.S. National Assessment based on the CCC and Hadley Centre simulations. We chose the CCC and Hadley Models to maintain the consistency with the national and regional assessments of the effects of global change (National Assessment Synthesis Team, 2001). For the temperature data, the monthly differences between the baseline and the scenario years of the GCM data were combined with the baseline climatology of the historical period to create the scenario climates for a transient greenhouse gas and aerosol simulation out to 2100. The precipitation scenarios were computed slightly differently. The ratio of future/historical was first computed and was used to multiply the baseline climatological precipitation data to derive the scenario climate. We derived the Estes Park climate data based on the VEMAP data. For selected decades of the scenario for the CCC and Hadley Models, we generated decadal monthly means of temperatures and precipitation. These were then used to compute the average decadal changes in climate variables. These changes were added (multiplied for precipitation) to the historical baseline period to generate a modified thirty-year climate adjusted for the climate change period in a potential future decade. So no changes in the interannual variability of the baseline historical period were created, only those changes in the modal amplitude of the monthly climate for a thirty-year period were used for a base case analysis. 213 IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS Table I Parameter estimates and Akaike weights for competing discrete logistic models of elk population at Rocky Mountain National Park. Tmin is the minimum temperature Models Twofactor Threefactor Intercept Summer precipitation Summer mean Tmin 8.63 0.05 –0.29 9.25 0.05 –0.30 Winter precipitation 0 –0.06 Equilibrium size Akaike weights (wr ) 1049 0.77 1039 0.23 We projected the sizes of elk population using the best approximating matrix population models driven by the CCC or Hadley based scenarios of the weather variables. This projection allows us to evaluate how the elk population would have responded if the weather conditions in the past 35 years had resembled what the CCC or Hadley Model predicted. We also ran projections in an additive manner by substituting one CCC Model-predicted weather variable at a time into the model. The order of the CCC Model-predicted weather variables into the model in the additive projection was determined by the sensitivity of the elk population to the change in a weather variable in a single weather factor matrix model. 3. Results 3.1. DISCRETE LOGISTIC MODEL Different combinations of weather variables provided 30 alternatives of the discrete logistic model. Only two models had appreciable support in the historical data. The best approximating logistic model included effects of summer mean monthly minimum temperature and summer mean monthly precipitation (Table I). Summer mean monthly minimum temperature was negatively related to elk population size, whereas summer mean monthly precipitation was positively related to elk population size. The second best model had three weather variables, i.e., summer mean monthly minimum temperature, summer mean monthly precipitation, and winter mean monthly precipitation. Winter monthly precipitation was inversely related to elk population size. The other two weather factors in the second best model had the same effects on elk population size as those in the two-factor model in terms of the signs of coefficients (Table I). The best approximating model offered a reasonable fit to the observed elk population size. The Akaike weight of the best approximating model was 0.77, and 214 GUIMING WANG ET AL. Figure 2. Observed and the best discrete logistic model – fitted population sizes of elk at Rocky Mountain National Park from 1965 to 1999. 30 out of the total 35 annual elk population sizes predicted by the best model fell within the 95% confidence intervals of the observed elk population size (Figure 2). The value of the Akaike weights of the second best model was 0.23. Therefore, the evidence supporting the two-factor model was three times stronger than the evidence in the data supporting the three-factor model. The two-factor logistic model estimated equilibrium population size for RMNP was about 1,040. Examining the two logistic models that were supported by the data suggests that the growth of the elk population at RMNP was retarded by increasing density, warmer summers, and wetter winters, but was accelerated by wet summers. 3.2. NONLINEAR MATRIX POPULATION MODELS Including different terms for weather effects and effects of density provided 62 candidate matrix models. Of these candidates, four models emerged with approximately equal support in the data, of which the Akaike weights were greater than 0.05 (Table II). All models included coefficients representing negative feedback between population density and both calf survival and recruitment. This indicates that the growth rate of the elk population at RMNP was attenuated by increasing density as a result of diminished calf survival rate and recruitment (Table II). Positive effects of summer precipitation on calf survival emerged in all models. However, two of the four models included negative terms for effects of summer mean monthly minimum temperature on recruitment, while two models included positive terms for effects of winter mean monthly minimum temperature. These results show some evidence supporting increased recruitment during warm winters IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 215 Figure 3. Observed and the best nonlinear matrix model – fitted population sizes of elk at Rocky Mountain National Park from 1965 to 1999. and reduced recruitment during warm summers. These results were consistent with the finding of the two-factor discrete logistic model in terms of the effect direction on the elk population size. We used model four (Table II) for projection of climate change because: (1) 26 out of 35 predicted points by model four fell within the 95% confidence intervals of observed elk population size (Figure 3); (2) the model predicted that the equilibrium elk population size was about 1,000, which was close to observations and the prediction of equilibrium size by the discrete logistic model; and (3) model four included all the three identified weather variables by the other three matrix models. 3.3. PROJECTION OF ELK POPULATION SIZES UNDER CLIMATE CHANGES The CCC and Hadley Models projected warming summers in the region surrounding RMNP during 2000–2035. The CCC Model projected that the summer monthly mean temperature will increase by 2.7 ◦ C and the winter mean monthly temperature by 2.6 ◦ C on average by 2035. The Hadley Model projected an average increase of 1.6 ◦ C in the summer mean monthly temperature over the next 35 years, and 1.5 ◦ C in the winter mean monthly temperature. However, the CCC and Hadley Models had opposite predictions on the change in the future precipitation. The CCC Model predicted dry years and that the mean summer monthly precipitation will decrease by 0.26 cm, and the winter mean monthly precipitation by 0.86 cm. The Hadley Model predicted wet years at Estes Park for the future 35 years, the summer mean monthly precipitation will increase by 3.5 cm, and the winter mean monthly precipitation by 2.03 cm on average. 216 Table II Parameter estimates and Akaike weights for competing nonlinear matrix models of elk population at Rocky Mountain National Park. Tmin is the minimum temperature, rc the proportion of females in calves, and SYM the male yearling survival rate. Model Intercept Summer mean Tmin Winter mean Tmin Density Summer precipitation Model 1 2.6 0 0 –0.004 0.09 Recruitment 0.07 0 0.05 –0.0009 0 Model 2 Calf survival 2.57 Recruitment 0.07 0 –0.03 0 –0.004 0.09 0 –0.0009 0 Model 3 Calf survival 2.10 0 0 –0.003 0.11 Recruitment 0.08 0 0 –0.0009 0 Model 4 Calf survival 2.50 Recruitment 0.07 0 –0.02 0 –0.004 0.09 0.05 –0.0009 0 SYM Akaike weights (wr ) 0.26 0.30 0.41 0.25 0.31 0.29 0.19 0.56 0.23 0.25 0.31 0.07 GUIMING WANG ET AL. Calf survival rc IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 217 Figure 4. Projection of elk population dynamics using nonlinear matrix model and Canadian Climate Centre model-based scenarios for the future climate change. Baseline is the prediction of elk population sizes using historic weather data. Tmin is the mean minimum temperature. The matrix model projections indicated that the equilibrium population sizes could reach 1,600 animals under the CCC scenario, and 2,000 animals under the Hadley scenario. The no-change scenario (based on historic, unadjusted weather) projected a steady state of 1,000 animals (Figure 4), which is close to current population estimates. Our projection, based on future climate scenarios, indicated that the increase in winter mean minimum temperature exerted effects on population dynamics that were more pronounced than effects of increased summer temperatures (Figure 5). 4. Discussion and Conclusion Our models portrayed the response of the elk population in RMNP to historic variation in weather. Coupling this portrayal with projections of a future climate suggests that the elk in RMNP could increase markedly as a result of enhanced survival and recruitment of juvenile animals. However, we emphasize that there are profound uncertainties in applying scenarios based on global models to ecological models operating at local scales. Thus, we do not offer these results as quantitative forecasts of what is likely to happen to elk numbers. Rather, we interpret our results as qualitative evidence that ecosystems that are stressed by overabundance of ungulates may be vulnerable to future climate changes. These vulnerabilities are supported by work of others. Post and Stenseth (1999) found that fecundity of female red deer in Norway was positively related to previous warm and wet winters. However, Post and Stenseth (1999) revealed complicated or contrasting 218 GUIMING WANG ET AL. Figure 5. Projection of elk population dynamic using nonlinear matrix model and Canadian Climate Centre model-based scenarios for the future climate change in an additive manner. Baseline is the prediction of elk population sizes using historic weather data. demographic responses to warm winters among 16 populations of seven ungulate species in the east and west of the Atlantic Ocean using multiple linear or autoregression analyses. Both our matrix and logistic models suggested that increases in summer mean monthly minimum temperatures lowered elk population growth rates or juvenile recruitment. Miller (1979) found that high temperatures and low rainfall during summer can enhance the growth of heather, and consequently result in negative effects on the nutritional state and fecundity of red deer hinds (Putman et al., 1996). Climate could indirectly affect elk populations through altering vegetation in their habitats. Moreover, Loison and Langvatn (1998) concluded that winter harshness had no consistent effect on calf winter survival over all ungulate species and all populations of a given species. The complicated relationships between ungulate demographics and climates need a process-oriented analysis to elucidate the underlying mechanism. We believe our models support the general conclusion that population sizes of elk are likely to increase in response to plausible scenarios for future climate. However, this conclusion requires caveats. First, we observed that warmer, drier conditions during the growing season could retard population growth rates while warmer, drier conditions during winter could accelerate them. Our models projected a net increase in growth rates and a subsequent elevation in steady state density because the accelerating effects of milder winter weather offset the attenuating effects of summer weather. This means that it is particularly important to understand how seasonal differences in future climate may be manifested. The disagreement among climate predictions on future patterns of summer precipitation causes large differences in population model projections. Second, our model is IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 219 based solely on empirical relationships among weather variables, population size, and population dynamics. We do not model explicit climate effects on vegetation, which could feed back to the elk population. Finally, projecting a future population exceeding current densities requires assuming that the slopes of the past relationships between density and vital rates remain constant. This assumption is probably unrealistic because larger populations might change vegetation conditions in such a way that density-dependent feedbacks would intensify. Such effects would cause our projections to overestimate steady states under altered climates. Climate changes could alter the vegetation succession of ungulate habitats (Miller, 1979; Houston, 1982). Under the CCC and Hadley scenarios, coniferous forests would expand substantially and move upward to elevations currently dominated by tundra. In high elevations (>3450 m) of RMNP area, the proportion of tundra in total area could be expected to decline from current 80% to 2% under the CCC scenario, and 13% under the Hadley scenario at climax stage (Hobbs et al., unpublished data). These vegetation changes would alter the summer and winter habitats enormously in RMNP. The elk used a high proportion of grassland and open pine habitats (Larkins et al., 1995). The limitation of foods in winter is often believed as the mechanism for the density dependent effects in ungulate populations (Houston, 1982). Future warming winters could expand elk winter ranges to high elevation areas and lift or reduce the limitation of winter ranges on the population sizes of elk in RMNP in a relative short term. However, in the long run, substantial decreases in the areas of tundra and grasslands, as suggested by the above climax stage, would deteriorate the food quality in both summer and winter ranges of elk and force elk to use more low-quality shrubs or trees. Therefore, climate change could affect the populations of elk not only directly, but also indirectly through altering the relative proportion of grassland to forest vegetation. Our results are consistent with findings of others on the importance of density dependence and weather in regulating ungulate population dynamics. Increases in densities of ungulate populations have been associated with reduced juvenile survival and recruitment in elk (Merrill and Boyce, 1991; Coughenour and Singer, 1996), white-tailed deer (Odocoileus hemionus) (McCullough, 1979), Norwegian red deer (Cervus elaphus) (Forchhammer et al., 1998), soapy sheep (Ovis aries), and Saiga antelope (Saiga tatarica tatarica) (Coulson et al., 2000). On the other hand, climate or weather variables have also been identified as density-independent factors affecting ungulate populations (Putman et al., 1996; Coulson et al., 2000). Merrill and Boyce (1991) showed that severe winters could depress the yearling recruitment and population growth of elk in Yellowstone National Park (YNP). Coughenour and Singer (1996) found that winter and summer rainfall had a significant positive effect on calf recruitment rates of elk in YNP. Singer et al. (1997) demonstrated that the calf survival rate of elk during the winter was positively related to the winter severity index, meaning that a mild winter was favorable of high winter survival of elk calves in YNP. Lubow et al. (2002) also showed that 220 GUIMING WANG ET AL. winter temperature had significant positive effects on recruitment of elk at RMNP. Smith and Anderson (1998) found that the duration of winters was inversely related to the calf survival rates of the Jackson Elk Herd. Negative effects of cold and harsh winters were also found on the demography of other northern ungulate species. Bartmann and Bowden (1984) reported an inverse relationship of the winter mortality of mule deer (Odocoileus hemionus) and winter temperature in Middle Park, Colorado. The survival rates of red deer calves of both sexes in Norway were positively correlated with winter temperatures (Forchhammer et al., 1998; Loison and Langvatn, 1998). Loison et al. (1999) thought that winter climate played a major role in determining body size and calf survival of red deer in Norway. Similarly, spring and winter temperatures had a positive effect on neonatal survival of bighorn sheep (Ovis Canadensis) in Canada when population density was high (Portier et al., 1998). In general, in large mammals, young individuals are more sensitive both to stochastic environmental variation and to density changes. Annual variability in survival rates was found mainly in calf survival rates (Gaillard et al., 1998). However, the above studies used either a comprehensive climate index (i.e., the North Atlantic Oscillation (NAO) or winter severity index) or a few weather variables (i.e., mean summer and winter precipitations and temperatures) as candidate predictors. We screened a full range of weather variables that could be affected under plausible scenarios for climate change. Although our results resemble other findings, our interpretations may diverge from earlier work. Weather effects are often treated as if they operated independently of density (Putman et al., 1996). In contrast, we view these regulating forces as highly contingent – the effects of density depend in a meaningful way on the effects of weather. If recruitment is elevated over the long term because of effects of favorable climate, then density feedbacks will need to intensify if equilibrium densities are to remain unchanged. This is to say that in the absence of changes in the relationship between density and vital rates, future climate change will weaken density dependent feedback. Our models suggest that negative feedback from population density on vital rates may be attenuated by amplifying effects of milder winters anticipated under most scenarios of future climates. Managers of national parks have relied on the effects of severe winter weather to regulate populations of native ungulates (Cole, 1971; Houston, 1971; Boyce, 1998). Our results suggest that these regulating effects may weaken if climate changes occur as widely predicted. This means that managers of national parks and conservation areas may need to reevaluate policies for sustaining natural processes in ecosystems where native ungulates have the potential to become excessively abundant. IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 221 Acknowledgements This study was supported by a grant from the Science to Achieve Results (STAR) program of United States Environmental Protection Agency (Agreement Number R 827449-01-0). Dr. H. Galbraith and two anonymous reviewers made helpful comments on the manuscript. References Bartmann, R. M. and Bowden, D. C.: 1984, ‘Predicting Mule Deer Mortality from Weather Data in Colorado’, Wildl. Soc. Bull. 12, 246–248. Berger, J.: 1991, ‘Greater Yellowstone’s U.S.A. Native Ungulates Myths and Realities’, Conserv. Biol. 5, 353–363. Boyce, M. S.: 1991, ‘Natural Regulation or the Control of Nature?’, in Keither, R. B. and Boyce, M. S. (eds), The Greater Yellowstone Ecosystem: Redefining America’s Wildness Heritage, Yale University Press, New Haven, pp. 193–208. Boyce, M. S.: 1998, ‘Ecological-Process Management and Ungulates: Yellowstone’s Conservation Paradigm’, Wildl. Soc. Bull. 26, 391–398. Buckland, S. T., Burnham, K. P., and Augustin, N. H.: 1997, ‘Model Selection: an Integral Part of Inference’, Biometrics 53, 603–618. Burnham, K. P. and Anderson, D. R.: 1998, Model Selection and Inference: A Practical Information – Theoretic Approach, Springer, New York. Chase, A.: 1986, Playing God in Yellowstone, Atlantic Monthly Press, Boston. Cole, G. F.: 1971, ‘An Ecological Rationale for the Nature or Artificial Regulation of Native Ungulates in Parks’, Tran. N. A. Wildl. Nat. Res. Conf. 36, 417–426. Coughenour, M. B. and Singer, F. J.: 1996, ‘Elk Population Processes in Yellowstone National Park under the Policy of Natural Regulation’, Ecol. Appl. 6 (2), 573–593. Coulson, T., Milner-Gulland, E. J., and Clutton-Brock, T.: 2000, ‘The Relative Role of Density and Climatic Variation on Population Dynamics and Fecundity Rates in Three Contrasting Ungulate Species’, Proc. R. Soc. Lond. B 267, 1771–1779. Decalesta, D. S.: 1994, ‘Effects of White-Tailed Deer on Songbirds within Managed Forests in Pennsylvania’, J. Wildl. Manage. 58, 711–718. Dickinson, R.: 1986, ‘The Climate System and Modeling of Future Time’, in Bolin, B., Doos, B., Jager, J., and Warrick, R. A. (eds), The Greenhouse Effects, Climate Change, and Ecosystems, Wiley, Chichester, U.K., pp. 207–270. Forchhammer, M. C., Stenseth, N. C., Post, E., and Langvatn, R.: 1998, ‘Population Dynamics of Norwegian Red Deer: Density-Dependence and Climatic Variation’, Proc. R. Soc. Lond. B 265, 341–350. Frontline Systems: 1999, Solver User’s Guide, Incline Village, Nevada. Gaillard, J. M., Festa-Bianchet, M., and Yoccoz, N. G.: 1998, ‘Population Dynamics of Large Herbivores: Constant Adult Survival and Variable Recruitment’, Trends Ecol. Evol. 13, 58–63. Garrott, R. A., White, P. J., and White, C. A. V.: 1993, ‘Overabundance: An Issue for Conservation Biologists?’, Conserv. Biol. 7, 946–949. Gill, R. M. A.: 1992, ‘A Review of Damage by Mammals in North Temperate Forests. 1. Deer’, Forestry 65, 145–169. Hess, K.: 1993, Rocky Times in Rocky Mountain National Park: An Unnatural History, University of Colorado Press, Boulder. Hobbs, N. T.: 1989, ‘Linking Energy Balance to Survival in Mule Deer: Development and Test of a Simulation Model’, Wildl. Monogr. 101, 39. 222 GUIMING WANG ET AL. Hobbs, N. T.: 1996, ‘Modification of Ecosystems by Ungulates’, J. Wildl. Manage. 60 (4), 695–713. Houghton, J. T., Meira Filho, L. G., Bruce, J., Lee, H., Callander, B. A., Haites, E., Harris, N., and Maskell, K.: 1995, Climate Change 1994 Radiative Forcing of Climate Change and An Evaluation of the IPCC IS 92 Emission Scenarios, Reports of Working Groups I and III of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge. Houghton, J. T., Meira Filho, L. G., Callander, B. A., Harris, N., Kattenberg, A., and Maskell, K.: 1996, Climate Change 1995 The Science of Climate Change, Contribution of WG I to the Second Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge. Houston, D. B.: 1971, ‘Ecosystems of National Parks’, Science 172, 648–651. Houston, D. B.: 1982, The Northern Yellowstone Elk, Macmillan Publishing Co., Inc., New York. Larkins, K., Clarke, J. A., and Olmsted, C. E.: 1995, ‘Patterns of Habitat Use and Distribution of Elk Adjacent to Rocky Mountain National Park’, J. Colorado-Wyoming Acad. of Sciences 27, 27. Loison, A. and Langvatn, R.: 1998, ‘Short- and Long-Term Effects of Winter and Spring Weather on Growth and Survival of Red Deer in Norway’, Oecologia 116, 489–500. Loison, A., Langvatn, R., and Solberg, E. J.: 1999, ‘Body Mass and Winter Mortality in Red Deer Calves: Disentangling Sex and Climate Effects’, Ecography 22, 20–30. Lubow, B. C., Singer, F. J., Johnson, T. L., and Bowden, D.: 2002, ‘Dynamics of Interacting Elk Populations within and Adjacent to Rocky Mountain National Park’, J. Wildl. Manage. 66 (3). McCullough, D. R.: 1979, The George Reserve Deer Herd, The University of Michigan Press, Ann Arbor. McNaughton, S. J.: 1983, ‘Serengeti Grassland Ecology: The Role of Composite Environmental Factors and Contingency in Community Organization’, Ecol. Monogr. 53, 291–320. McNaughton, S. J.: 1985, ‘Ecology of a Grazing Ecosystem: The Serengeti’, Ecol. Monogr. 55, 259–294. Mech, L. D., McRoberts, R. E., Peterson, R. O., and Page, R. E.: 1987, ‘Relationship of Deer and Moose Populations to Previous Winters’ Snow’, J. Animal Ecol. 56, 615–627. Merrill, E. H. and Boyce, M. S.: 1991, ‘Summer Range and Elk Population Dynamics in Yellowstone National Park’, in Keiter, R. B. and Boyce, M. S. (eds.), The Greater Yellowstone Ecosystem: Redefining America’s Wildness Heritage, Yale University Press, New Haven, pp. 262–273. Miller, G. R.: 1979, ‘Quantity and Quality of the Annual Production of Shoots and Flowers by Calluna-Vulgaris in Northeast Scotland’, J. Ecol. 67, 109-129. National Assessment Synthesis Team: 2001, Climate Change Impacts on the United States: The Potential Consequences of Climate Variability and Change, Report for the U.S. Global Change Research Program, Cambridge University Press, Cambridge. Noy-Meir, I: 1975, ‘Stability of Grazing System: An Application of Predator-Prey Graphs’, J. Ecol. 63, 459–481. Portier, C, Festa-Bianchet, M., Gaillard, J. M., Jorgenson, J. T., and Yoccoz, N. G.: 1998, ‘Effects of Density and Weather on Survival of Bighorn Sheep Lambs (Ovis canadensis)’, J. Zool. 245, 271–278. Post, E. and Stenseth, N. Chr: 1999, ‘Climatic Variability, Plant Phenology, and Northern Ungulates’, Ecology 80, 1322–1339. Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: 1992, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, Cambridge. Putman, R. J., Langbein, J., Hewison, A. J. M., and Sharma, S. K.: 1996, ‘Relative Roles of DensityDependent and Density-Independent Factors in Population Dynamics of British Deer’, Mammal Rev. 26 (2–3), 81–101. Sanna, P.: 1997, Using Visual Basic for Applications 5, Special Edition, QUE, Indianapolis. Samuel, M. D., Garton, E. O., Schlegel, M. W., and Carson, R. G.: 1987, ‘Visibility Bias during Aerial Surveys in Northcentral Idaho’, J. Wildl. Manage. 51, 622–630. IMPACTS OF CLIMATE CHANGES ON ELK POPULATION DYNAMICS 223 Shaver, G. R., Canadell, J., Chapin III, F. S., Gurevitch, J., Harte, J., Henry, G., Ineson, P., Jonasson, S., Melillo, J., Pitelka, L., and Rustad, L.: 2000, ‘Global Warming and Terrestrial Ecosystems: A Conceotual Framework for Analysis’, BioSciences 50, 871–882. Singer, F. J., Harting, A., Symonds, K. K., and Coughenour, M. B.: 1997, ‘Density Dependence, Compensation, and Environmental Effects on Elk Calf Mortality in Yellowstone National Park’, J. Wildl. Manage. 61 (1), 12–25. Singer, F. J., Swift, D. M., Coughenour, M. B., and Varley, J. D.: 1998, ‘Thunder on the Yellowstone Revisited: An Assessment of Management of Native Ungulates by Natural Regulation, 1968– 1993’, Wildl. Soc. Bull. 26 (3), 375–390. Smith, B. L. and Anderson, S. H.: 1998, ‘Juvenile Survival and Population Regulation of the Jackson Elk Herd’, J. Wildl. Bull. 62, 1036–1045. Steinhorst, R. K. and Samuel, M. D.: 1989, ‘Sightability Adjustment Methods for Aerial Surveys of Wildlife Populations’, Biometrics 45, 415–425. Stromayer, K. A. K and Warren R. J.: 1997, ‘Are Overabundant Deer Herds in the Eastern United States Creating Alternate Stable States in Forest Plant Communities?’, Wildl. Soc. Bull. 25, 227– 234. Watson, R. T., Zinyowera, M. C., and Moss, R. H.: 1996, Climate Change 1995 Impacts, Adaptations and Mitigation of Climate Change: Scientific-Technical Analyses, Contribution of Working Group II to the Second Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge. White, G. C. and Lubow, B. C.: 2002, ‘Fitting Population Models to Multiple Sources of Observed Data’, J. Wildl. Manage 66, 300–309. (Received 10 July 2000; in revised form 26 October 2001)