## Abstract

Assessment of future ecosystem risks should account for the relevant uncertainty sources. This means accounting for the joint effects of climate variables and using modelling techniques that allow proper treatment of uncertainties. We investigate the influence of three of the IPCC's scenarios of greenhouse gas emissions (special report on emission scenarios (SRES)) on projections of the future abundance of a bryophyte model species. We also compare the relative importance of uncertainty sources on the population projections. The whole chain global climate model (GCM)—regional climate model—population dynamics model is addressed. The uncertainty depends on both natural- and model-related sources, in particular on GCM uncertainty. Ignoring the uncertainties gives an unwarranted impression of confidence in the results. The most likely population development of the bryophyte *Buxbaumia viridis* towards the end of this century is negative: even with a low-emission scenario, there is more than a 65 per cent risk for the population to be halved. The conclusion of a population decline is valid for all SRES scenarios investigated. Uncertainties are no longer an obstacle, but a mandatory aspect to include in the viability analysis of populations.

## 1. Introduction

Our ability to anticipate future ecosystem changes is of great concern, particularly under the scope of climate change [1]. In projecting ecosystem changes into the future, it is important to account for the relevant sources of uncertainty to provide a robust picture for evaluating risks [1–3]. This means accounting for the joint effects of climate and ecosystem variables, and using up-to-date modelling techniques that allow proper treatment of uncertainties [1]. In this way, we can quantify the overall projection uncertainty that arises from interactions among model formulations, parameter estimates and their inherent uncertainties [4,5]. In this paper, we show the influence of sources of variability and uncertainty in climate models and in a population dynamics model on population projections using an epixylic bryophyte as a model species.

Global climate models (GCMs) are designed to represent the climate system on a global scale, and to provide a realistic description of the large-scale climatic features. However, they are incapable of describing the regional or local climate characteristics owing to their coarse spatial resolution. It is therefore necessary to downscale the GCMs to provide projections at regional or local scales. The required improvement in spatial resolution can be achieved by dynamical downscaling using regional climate models (RCMs). By design, RCMs maintain the physical consistency among the variables produced by the GCMs [6]. Thus, using RCMs for dynamical downscaling of global climate scenarios, which are produced by GCMs, is often a necessary step for enabling analyses of climate change impacts.

The variation in projections by climate models with different approaches for modelling physical processes (model formulation) or different parameter sets are often larger than the variation resulting from the climate scenarios investigated by a specific model [2,7–9]. Averaging over ensembles of climate models accounts for this uncertainty [10]. For example, the impact of climate change on the distribution of different taxa have been averaged among climate models and also among different types of model for the distribution of the taxa [11,12]. Fewer studies account for the uncertainty in the biological process models [2,13–15]. However, none of these studies increases our understanding of how uncertainty sources in the climate and ecosystem modelling influence the projections of the biological systems.

The influence of climate model uncertainties is well studied by climatologists [16–20], but our understanding of their impacts on projections of biological systems is scarce. We use the whole chain GCM—RCM—population dynamics model to, first, investigate the influence of three greenhouse gas emission scenarios described by the IPCC (SRES scenarios henceforth) [21], and in addition, to investigate the relative importance of five uncertainty sources on population projections: (i) the probability distribution of the population dynamics parameters, (ii) the GCM formulation, (iii) the climate sensitivity, essentially the change in equilibrium global mean temperature resulting from a doubling of the greenhouse gas emissions [22], (iv) the models' initial conditions resulting in variability at decadal and longer time-scales, and (v) the spatial resolution of the regional downscaling. We used an ensemble [18] based on the Rossby Centre RCM RCA3 [20] to downscale a relatively large set of GCM outputs used. Using a single RCM limits the information regarding the performance of different RCMs in future climates, but makes it easier to disentangle other aspects of the uncertainty [13]. Moreover, since the weather variables used in this study are aggregated to monthly and seasonal values, the difference between different RCMs can be assumed to be less pronounced [16,23].

Although projections of species or ecosystems into the future have increased in the last years [2,14,24–26], there are large biogeographic and taxonomic gaps in the scientific literature [27,28]. An example of a large poorly studied group are non-vascular plants, such as bryophytes, which have a key role in ecosystem processes [29,30], and constitute a large proportion of the biodiversity of many biomes [31]. Often, they are not even mentioned in reviews [28,32]. For a common ground floor species investigated (*Hylocomium splendens*), a positive response to climate change has been suggested [30]. However, epixylic species, such as our model species, may be more sensitive as the microclimate of their substrate, dead wood, which is less buffered against climatic variation [29]. Moreover, many epixylic species are short-lived which is assumed to be associated with sensitivity to increments in environmental variability [33].

The aim of this study was to investigate the influence of different SRES scenarios, and to compare the variability among them with the variability introduced by uncertainty sources inherent in the climate projection and in the population dynamics model on estimates of the future abundance of a population. We illustrate this by projecting the population dynamics of the model species *Buxbaumia viridis*, a short-lived epixylic moss of conservation concern. We make use of an ensemble of 13 local climate scenarios produced by one RCM forced by five GCMs to drive a weighted average of two hierarchical Bayesian models for the inter-annual variation of the species [3]. We ask: (i) what is the relative importance of different uncertainty sources on the estimate of future population viability, and (ii) what is the probability of decreased population abundance in the future given the different uncertainty sources?

## 2. Material and methods

### (a) Model species

*Buxbaumia viridis* (DC.) Moung. & Nestl. is a dioicous moss that grows on strongly decayed dead wood patches in coniferous or deciduous forest [34]. It has a circumboreal distribution, but it appears rather uncommon wherever it occurs [35]. The species is included in the Annex II of the EU Habitats Directive [36]. The model applied [3] is based on data on abundance of mature sporophytes of *B. viridis* collected between 1996–2003 and 2008, from six local sampling plots of 25 × 25 m in a boreo-nemoral forest situated in eastern central Sweden (Vipängen: 59°49′ N, 17°39′ E).

### (b) Modelling the sporophyte population dynamics

We use an evaluated, weighted average of two Bayesian hierarchical generalized linear models for simulating the inter-annual variation of sporophyte abundance of *B. viridis*. See the study of Ruete *et al*. [3] for further model details*.* The model has a linearized form of an exponential population growth model as the deterministic skeleton of population regulation, which is reasonable for bryophytes [37]. Specifically, we assume that the local plot-scale, *p*, abundance of sporophytes in year *t* followed a Poisson distribution, *N*_{y,p} ∼ Poisson(*μ*_{y,p}), with mean sporophyte abundance *μ*_{y,p}, where,
2.1where *N _{y,p}* is the observed abundance of sporophytes, and

*N*

_{y}_{−1,p}is the abundance of sporophytes in the preceding year. The ‘effect size’ parameter

*β*

_{1}determines the proportion of spores produced by the sporophytes in the preceding year that gives rise to sporophytes in the focal year. The parameter

*ɛ*

_{y,p}is a normally distributed overdispersion. Moreover, the logarithm of the regional, year-specific growth rate

*λ*

_{y}, is a linear function of the number of preceding autumn frost days (

*F*), the sum of spring precipitation (

*P*), and the mean spring temperature (

*T*). Specifically, , and 2.2where

*σ*

_{λ}quantifies the local-scale variation of log

*λ*

_{y}, unknown regional processes such as unexpected interaction effects between weather variables, or measurement error,

*φ*

_{0}is an ‘intercept’ parameter and

*φ*

_{i}|

*i*= 1, 3 are ‘effect size’ parameters.

We conducted simulations of sporophyte abundance at the local scale for the period 1961–2098, driven by the 13 climate projections (table 1). Each climate projection provides yearly values of the climate variables (*F _{y}*,

*P*and

_{y}*T*; equation (2.2)) that, together with the preceding years population size (

_{y}*N*

_{y}_{−1,p}; equation (2.1)) drive the population projection (

*N*

_{1961–2098,p}). Moreover, we accounted for regional and local-scale variation and data uncertainty in the simulations by, for each of the 2000 replicates, drawing a new set of values of the model parameters (e.g.

*φ*

_{i}|

*i*= 1, 3) from the joint full posterior distribution of all model parameters, which are standard output from Bayesian estimation [3]. Simulations were started off in 1960 with the plot-specific median levels observed in years 1996–2003,

*N*0 = {4, 3, 5, 2, 7, 6}.

In order to disentangle the impact of parameter variability and process uncertainty on the population projections, and to isolate the uncertainty inherent in the population dynamics parameters from the uncertainty present in the different sets of climate projections, we also conducted a second set of simulations (2000 replicates) using only the mode (i.e. point estimates) of the estimated parameters of the population dynamics model. These simulations only include deterministic estimates (the point estimates) of the effect of the climate variables on the population dynamics (equations (2.1) and (2.2); e.g. *φ*_{i} | *i* = 1, 3), a deterministic effect of the abundance of the preceding year (*β*_{1}; equation (2.1)), a deterministic estimate of the expectancy for the normally distributed residual variation (*ɛ*_{y}_{,p}; equation (2.1)), but the stochastic variability inherent in the Poisson nature of the sporophyte abundance (demographic stochasticity modelled by the Poisson distribution; equation (2.1)). These simulations thus exclude the variation or uncertainty accounted for when each replicate is a draw from the joint posterior distribution of all model parameters (see above). All simulations were run using the software R v. 2.13.0 [38].

### (c) Climate datasets

We obtained observed local weather data for the period 1961–2008 from Ultuna weather station for calibrating the RCM, and for evaluating the accuracy of this calibration (see §2*e*). The station is located 100 m north of the study forest. The weather data were summarized as the number of days with a temperature below 0°C for the period October–November of the preceding year, and as the sum of precipitation (mm) and the mean of monthly temperature means (°C) for the period March–June.

We studied five sources of variation and uncertainty in climate modelling (table 1) using an ensemble of regional climate change scenarios produced by the Rossby Centre regional atmospheric climate model RCA3.0 [20,39]. See Kjellström *et al.* [18] and Nikulin *et al.* [19] for a presentation of the ensemble members along with an ensemble analysis of the projected climate change. The main target of the ensemble is the uncertainty due to different driving coupled atmosphere–ocean GCMs under the SRES A1B scenario [21,40] that was greatly used in the most recent IPCC fourth assessment [41]. For further details of GCM characteristics, see overviews in Randall *et al*. [40] and Meehl *et al*. [42] and, specifically for HadCM3, Collins *et al.* [10]. Additional emission scenarios included in the ensemble are the SRES A2 (comparatively high emissions) and B1 (optimistic low emissions) used to force ECHAM5.

For a given SRES scenario, different GCMs will not project the same global warming trends owing to different approaches for modelling specific physical processes or due to choices of different model parameter settings. A collective effect of model differences is that different GCMs will have different climate sensitivity, which is defined as the change in equilibrium global mean temperature that results from a doubling of the greenhouse gas emissions [22]. One can explore the effects of different parameter settings on the climate sensitivity of a single GCM by varying the settings and comparing the resulting climate sensitivities. We account for this source of variation based on the ensemble [18], specifically using the three runs with the Hadley Centre GCM, where HadCM3-Q0 is the reference run using optimal parameter settings, HadCM3-Q3 is the simulation with parameter settings representing lower climate sensitivity, and HadCM3-Q16 is the simulation representing increased high climate sensitivity [10].

Even for a given greenhouse gas emission scenario and a specific parameter setting of a GCM, runs may give different answers depending on differences in the specification of the initial conditions. Essentially, differences in initial conditions of model runs result in variation of the slow components of the systems, mainly the oceans. They will not be in phase with each other, or with the observed variability at (multi-)decadal time-scales. When summarizing projections for some period, say 2059–2098, some variation between different runs will be due to this (multi-)decadal scale variability. This variation reflects natural variability which should be accounted for [19]. We do this by using three different runs with the exactly same set-up of ECHAM5 (−r1, −r2 and –r3) forced by the A1B scenario and differing only in the initial conditions.

Most regional climate change scenarios of the ensemble have a spatial resolution of about 50 km covering the same European-wide domain. However, for the ECHAM5-r1 run forced by emission scenario A1B, the ensemble includes regional downscalings at 50, 25 and 12 km, and this allowed us to investigate the effect of the spatial resolution of the RCM on the projected population abundance. When using monthly and seasonally aggregated data, as we do in the current study, large effects are less likely but we wanted this confirmed and tested it. The representation of extreme situations, which show large spatial variability, depends on the spatial resolution of the ensemble members [23].

While the RCM provides higher resolution data than the GCM, there is still a scale mismatch between the point station data used as a reference for calibration of the ecological model and climate model output. That is, the station, Ultuna, may not be representative for the whole grid cell area represented by the RCM output. Thus there is a need for a further downscaling and calibration step using statistical methods. This procedure draws on the same general methodology as that commonly employed for direct empirical statistical downscaling of GCM data, which includes corrections of biases introduced by both the GCM and the RCM. Here, we use the distribution-based scaling (DBS) method to correct this bias [43]. The correction factors are derived by comparing the RCM output to the observations in the reference period (here the station data from 1960 to 2008). The factors are used to adjust RCM outputs to make them statistically comparable with the observed weather data in terms of mean and standard deviation. For precipitation, rainfall probability is also concerned. Derived correction factors are then also applied to the rest of the time series for the future in the same climate projection. The method assumes that the biases are systematic and constant for the entire climate projection. To implement DBS adjustment, it is important to use long enough local weather observation data. These requirements are met by the Ultuna station data.

### (d) Simulation summary

We summarize the simulated sporophyte abundance as the summed abundance among the six sample plots for three 40-years periods: 1961–2000 (reference period), 2019–2058 (near future) and 2059–2098 (far future). By having sampled the joint probability distribution of all parameters of the ecological model in the simulations, we obtain the probability distribution of the abundance for each year of simulation. The 50 and 90 per cent highest posterior density intervals (HPDIs) for the posterior distributions of the abundance were calculated around the median.

We calculated the probability of change in population abundance for the two 40-years periods, and for the different sources of uncertainty in climate projections (table 1). These probabilities were calculated as the percentage change in simulated sporophyte abundance by comparing the abundance for each year of the early and far future for each replicate (40 yr × 2000 replicates) with the abundance in a random year and replicate (without replacement) for the reference period. We then calculated the cumulative probability of decline (from 0 to −100%). As a reference probability distribution, we finally calculated the same probability for the reference period against itself. The most likely abundance change is where the difference between the reference probability distribution curve and a focal probability distribution curve is the greatest.

### (e) Model evaluation

We evaluated the performance of the population simulations and how the observational variability in weather data for the reference period was reflected in the downscaled data by comparing the probability of change in population abundance between simulations using observed weather and using simulated weather based on downscaled projections. As expected for well-calibrated RCM projections, the simulations driven by projected climatic data for the reference period gave similar population projections as the simulations driven by observed data for the same period (electronic supplementary material, figure S1). We could therefore safely investigate the difference between the projected population abundance for the reference period with that of the early and far futures.

## 3. Results

We show a future decline of the *B. viridis* abundance under the three IPCC SRES scenarios forcing ECHAM5 (figure 1). The median abundance was markedly reduced towards the end of this century under all SRES scenarios (figure 1) and other uncertainty sources (electronic supplementary material, figure S2). The decline in population abundance is also illustrated by an increase in probability of population decline to low levels (figure 2). The most likely decline in the far future (2059–2098) is large: a reduction of 59 per cent (±s.d. 7.9) compared with the reference period when averaging over SRES scenarios and uncertainty sources (figure 2). Moreover, there is a more than 65 per cent risk for the population to be halved (i.e. decline by 50% or more) in the far future under SRES B1 scenario (figure 2). The simulations using point estimates of the parameters show the same patterns of decline, but underestimate the uncertainty in the population projections. Therefore, the probabilities of decline were higher, i.e. seemingly more certain, than those using simulations with the full probability distribution of the parameters (figure 2).

The order of severity among the SRES scenarios follows the expectations (B1 least severe, A2 most severe; figure 2*a,b*), but this is not noticeable until the far future. In the near future (2019–2058), the population abundance can be expected to decrease by 49 per cent or more assuming SRES A2 scenario, by 33 per cent or more assuming A1B and by 37 per cent or more assuming B1. In the far future (2059–2098), the corresponding decreases are 66, 66 and 49 per cent, respectively.

The greatest source of uncertainty is GCMs formulation. The variation in the probability of population decline in the near future among different SRES scenarios is smaller than the variation among simulations of the A1B scenario driven by different GCMs formulations (figure 2*c*,*d*). However, this difference in variation among emission scenarios and GCM formulation tends to disappear in the far future: there is similar variation among GCM formulations as among emission scenarios. The projections based on the GCM CNRM deviate from the others, as this model predicts considerably lower probabilities of decline and higher probabilities of increase in the near future than the other GCM formulations (figure 2*c*,*d*).

The model with the reference level in climate sensitivity (HadCM3-Q0) presented the highest increase in probability of decline in the near future. All probabilities of decline increased in the far future, but the difference between HadCM3-Q0 and HadCM3-Q16 disappeared (figure 2*e*,*f*).

The simulations of population abundance using different initial conditions in the GCMs, representing natural variability, show larger variation in the near future than in the far future (figure 2*g*,*h*). The variation introduced by initial conditions in the near future is even larger than the variation introduced by emission scenarios in this period.

As expected for variables that do not represent extreme weather events, there was very little variation in probability of decline among simulations driven by a GCM (ECHAM5-r3) with different spatial resolutions (figure 2*i*,*j*), and this variation was small in comparison with the variation introduced by other uncertainty sources.

## 4. Discussion

Our results show that the strength in statements about population changes in the future may vary significantly depending on the number of uncertainty sources that are accounted for in the projections. We should be aware of this for providing robust pictures for evaluating risks. Population viability analyses acknowledging the uncertainty in climate change are increasing [2,11–14,24,44–46], but the influence of different uncertainty sources on the statements about the population development has, to our knowledge, never been disentangled. Uncertainties are no longer an obstacle, and should be incorporated when assessing the robustness of projected trends. This additional information helps to prioritize among competing demands, and should be a mandatory aspect in the assessment of populations. Increased understanding of the relative importance of different uncertainty sources also helps choosing which sources to include when conducting studies on the impact of climate change.

We have expressed the uncertainty in the form of probabilities of change in population abundance to different levels. These probabilities of reaching different levels can be faced against other interests, and are hence straightforward to apply in conservation planning and decision-making. We could report probabilities of declines because we used a population model that was fitted using the Bayesian statistical modelling approach [47]. Thereby, the natural variability and uncertainty in the data used for model development were included in the parameter estimates, and transferred into uncertainties in the projections of population abundance [1–3]. This approach can be applied to any species whose dynamics or distribution pattern can be generalized using a probabilistic model, e.g. a population dynamics model that is fitted to observed data on dynamics. Ignoring these uncertainty sources generally gives an unwarranted impression of confidence in results. Here, the interpretation would be biased towards certainty in population decline (cf. e.g. figure 2*a*,*b*).

Three out of four of the studied sources of uncertainty related to climate modelling contribute to the uncertainty in the projected levels of *B. viridis* abundance. First, results are mainly sensitive to the GCM formulation. The variation in probability of population decline in the near future among GCMs is larger than the variation in probability of decline in the near future given by different SRES scenarios (figure 2). A limitation with this comparison is that a single GCM is used to compare scenarios. However, the relatively large variability introduced by different GCMs in the near future agrees with earlier population projections and analyses of variation among climate models [2,18]. The deviating projection provided by the downscaled CNRM GCM has been previously reported [18,20].

Second, the GCM parametrization (i.e. climate sensitivity) had influence on the population projections (figure 2*e*,*f*). This is illustrated by the population simulations driven by the three HadCM3 runs, and the same has been found in pure climate simulations [18]. The model with the reference level in climate sensitivity (HadCM3-Q0) presented the highest increase in probability of decline in the near future. All probabilities of decline increased in the far future, but the difference between HadCM3-Q0 and HadCM3-Q16 disappeared (figure 2*e*,*f*).

Third, the natural variability in the slow components of the climate system, mainly the oceans phase, is an important source of uncertainty (figure 2*g*,*h*). We represent this source by variability in initial conditions. However, its importance decreases the longer the simulations are run, in accordance with pure climate simulations [18]. Moreover, its importance seems to be lower than GCM formulation or climate sensitivity (figure 2).

The fourth source of uncertainty, the spatial resolution of the RCM had a negligible influence on the variation in probability of bryophyte population decline (figure 2*i*,*j*). The spatial resolution is known to be important when focusing on extreme climate events showing large spatial variation [19], but not when using climate data that have been aggregated over longer time periods, as we do.

The runs conducted with ECHAM5 allowed us to explicitly quantify the relative importance of two sources of uncertainty, and show that the initial conditions assumed for GCMs may have greater influence on the projected abundance for the near future than the emission scenarios (figure 2). In contrast, the spatial resolution may be arbitrarily chosen.

### (a) Implication for *Buxbaumia viridis* dynamics and viability

There is a clear risk for a large decline of the *B. viridis* population under all the greenhouse gas emission (SRES) scenarios investigated. The difference in decline among emission scenarios will not be evident until the far future. The order of probability of population decline among the SRES scenarios follows the expectations (B1) is least severe. The most likely event (the greatest increment in probability) concerns great population declines.

We could not have anticipated the species' response to the changing climate without making simulations which explicitly accounts for the relative importance of the climate variables during different times of the year. The assumed main climate variables driving the *B. viridis* dynamics are namely frost days in the autumn and moisture balance in spring, as reflected by a positive effect of precipitation and a negative effect of temperature [3]. In the coming century, both the mean and the variability in all seasons' temperature, and spring precipitation are projected to increase [7]. Even more, in Scandinavia there is, and will be, a positive correlation between the change in temperature and in precipitation [18]. We now show that the negative effect of increased spring temperature will be greater than the positive effect of increased precipitation.

The projected decline in *B. viridis* abundance may reflect a future decline in additional similar species. In particular, we have projected future changes in the moisture balance, and it is likely that other species which are sensitive to changes in moisture conditions in a way similar to *B. viridis* will also respond in the same way. It is known that the RCM used in this study (RCA3.0) is biased by overestimating the summer precipitation in northeastern Europe [18]. This means that the future of this and similar species may be even worse. Moreover, with this change in the moisture balance, we may expect lower future levels of dead wood as a result of an increased rate of decomposition [48]. A way to alleviate this projected decline is to increase the amount of dead wood in managed forest ecosystems [29]. Moreover, even with this future prospect, a change in greenhouse gas emission policies and the social awareness contemplated by the SRES B1 scenario will help to ameliorate the negative effects of anthropogenically caused climate change on this and many other species.

## Acknowledgements

We are thankful to Per Nyman for providing us with the weather data from Ultuna weather station, and to three anonymous reviewers for providing constructive comments on the manuscript. The work was funded by grants 2005–933 and 2006–2104 from FORMAS to T.S. L.B. and W.Y. contributed to this study under the Mistra-SWECIA research programme funded by the Swedish strategic research foundation Mistra.

- Received February 24, 2012.
- Accepted March 8, 2012.

- This journal is © 2012 The Royal Society