Modeling the forest dynamics of the Sierra Nevada under climate change using SORTIE-ND

Model simulation results suggest that forests in the Sierra Nevada mountains of California will tend to increase in density and basal area in the absence of fire over the next century, and that climate change will favor increases in drought-tolerant species. Climate change is projected to intensify the natural summer drought period for Mediterranean-climate forests. Such changes may increase tree mortality, change species interactions and composition, and impact ecosystem services. To parameterize SORTIE-ND, an individual-based, spatially explicit forest model, for forests in the Sierra Nevada, and to model forest responses to climate change. We use 3 downscaled GCM projections (RCP 8.5) to project forest dynamics for 7 sites at different elevations. Basal area and stem density tended to increase in the absence of fire. Climate change effects differed by species, with more drought-tolerant species such as Jeffrey pine (Pinus jeffreyi A.Murray bis) and black oak (Quercus kelloggii Newb.) exhibiting increases in basal area and/or density. Increasing forest density may favor carbon sequestration but could increase the risk of high-severity fires. Future analyses should include improved parameterization of reproduction and interactions of disturbance with climate effects.


Introduction
Climate change impacts on forests have been a major research focus due to their importance for carbon storage (Oren et al. 2001;Bonan 2008;Earles et al. 2014) and other crucial functions (Flint et al. 2013;Goulden and Bales 2014;Grossiord et al. 2014), as well as the concern that long generation times could limit tree responses (Rice and Emery 2003; Aitken and Whitlock 2013). The forests of the Sierra Nevada mountains in California are adapted to a dry summer/wet winter Mediterranean climate. However, climate change is expected to intensify the water cycle, resulting in more extreme drought and precipitation (IPCC 2013;Reidmiller et al. 2018). Recently, the "hot drought" of 2012-2016 (Swain et al. 2014) was followed by extremely wet years in 2017 (Murray and Lohman 2018) and 2019 (CDWR 2019).
Such changes are expected to impact forest demography. "Background" mortality rates of adult trees in the western USA have increased due to temperature-associated increases in aridity (van Mantgem et al. 2009). Extreme droughts can produce mass mortality events (Allen et al. 2010); The 2012-2016 California drought killed over 100 million trees, with some local mortality rates exceeding 60% (Fettig et al. 2019). Growth of all age classes and seedling recruitment are affected by climate too. Our recent analysis found that hotter summer temperatures reduced the survival and growth of most species of Sierra Nevada tree seedlings (Moran et al. 2019). Shifts in forest composition have already been observed. A combination of fire suppression and climate change led to increasing forest density, decreasing numbers of large trees, and a greater abundance of shade-tolerant conifers during the twentieth century (Dolanc et al. 2014b;McIntyre et al. 2015). At lower elevations, oaks have become more abundant relative to pines (Dolanc et al. 2014b;McIntyre et al. 2015), and models predict further shifts in from needle-leaf to broadleaf species (Lenihan et al. 2003;Liang et al. 2017).
Forest responses depend on the demographic responses of different species and life stages, and on species interactions. Because of the long timescales involvedboth for tree lifecycles and the timelines of projected climate changesimulation models are crucial for better understanding forest dynamics in a changing environment. Individual-based forest simulators account for the interaction of each tree with the environment and the other individuals and species on the landscape, bridging the gap between small-scale individual-based studies and landscape studies. No individual-based forest simulator has been previously parameterized for the Sierra Nevada.
SORTIE-ND (http://sortie-nd.org/index.html) was originally developed for broadleaf forests in the Eastern USA (Pacala et al. 1996;Martin et al. 2010), but has since been used for multiple forest types (Uriarte et al. 2009;Bose et al. 2015;Ameztegui et al. 2017). Competitive interactions, based on the sizes of neighboring trees, can be accounted for. SORTIE-ND has accessible source code and a website (http://www.sortie-nd.org) where users can share newly developed functions. This accessibility was the main reason we chose this model for our study; comparisons to other models are discussed in Appendix 1. This paper documents the species-level parameterization of SORTIE-ND; we plan to incorporate individual-level variation and heritability of climate responses in future iterations.
We test the performance of the model by hindcasting the dynamics of three plots at different elevations. Then we project the forest responses to changing climate over the next century for seven plots at different elevations, using downscaled forecasts from the most recent Coupled Model Intercomparison Project (CMIP5) (IPCC 2013). We predicted that: 1. Total basal area will increase over time in most of these scenarios, as suggested by past densification and the positive growth responses of many species to warmer winters (Aubry-Kientz and Moran 2017). 2. Changes in tree density will depend on initial stand structure (many small trees that undergo self-thinning vs. a few larger trees). 3. Species that are found in hotter/drier environments or that show more positive growth or survival response to higher temperatures will tend to increase in abundance and/or basal area in the climate change scenarios relative to the control scenario.

The simulator
In SORTIE-ND each tree has a DBH (diameter at 1.4 m), species, xy coordinates, and status: seedling, sapling (<1.4 m tall, non-reproductive), or adult. Minimum adult DBH was computed with a linear model predicting DBH from age calibrated with Forest Inventory and Analysis (FIA) data (https:// www.fia.fs.fed.us/), and the species' minimum reproductive age (Burns and Honkala 1990). At each time step, trees reproduce, grow, and/or die; behaviors can depend on environmental variables. Values for all parameters and names of SORTIE-ND behaviors are given in Appendix 1.

Choice of climate variables
We limited the number of climate variables per behavior to two for simplicity and computational efficiency. For adult growth, the two best climate predictors were January minimum temperature and precipitation (Aubry-Kientz and Moran 2017) and for adult mortality July maximum temperature and precipitation (Appendix 2). For the survival of seedlings >10 cm in height, the best-fit model included July maximum temperature and total precipitation anomalies averaged over the current and past year, while for growth the current year July maximum temperature and total snow performed better (Moran et al. 2019). However, projections for snow are highly uncertain because they combine already uncertain precipitation projections with calculations of how much falls as snow and the melt rate. As current year precipitation also predicted seedling growth well, we used this and July maximum temperature for both seedling growth and survival. Table 1 shows the direction of climate effects on seedlings versus trees/sapling growth and survival.

Growth
Annual diameter growth (ΔD j,t ) for adults and saplings is computed as: where i indicates individual, s species, and t year; N i,t the total neighbors of i in year t; D j,t the DBH of neighbor j; dist i,j the distance between i and j; JMn t the minimum January temperature; and P t the precipitation. The β's are estimated parameters, with σ g 2 representing the error term.
Seedling behaviors in SORTIE-ND are usually based on diameter at 10 cm height. As our data were recorded by height bin without diameter, we created four size classes: 0 (new seedlings); 1 (> 2 years old but < 10 cm); 2 (10-50 cm); and 3 (50-140 cm). Transitioning to the next size class is a Bernoulli process: where θ i,t is the probability of transitioning; BA t is the total adult tree basal area in a 10-m radius; and JMx t is the maximum July temperature. β 0,s represents a size-specific intercept for size class s. Our model was fitted to seedlings that were > 10 cm tall. The size-specific parameter for class 1 was therefore based on seedlings that died back. The same parameter was applied to class 0.

Mortality
Adult and sapling mortality is modeled as a Bernoulli process: where θ i,t is the probability of dying (Appendix 2). For seedlings, survival is also a Bernoulli process: where θ i,t is the probability of surviving. Class 0 seedlings were initially given the size effect estimated for seedlings <10 cm (Moran et al. 2019), but this resulted in an unrealistic increase in tree density (A4) because the established seedlings ( >10 cm) used to fit the model tend to have much higher survival than younger seedlings. Therefore, the class 0 and 1 size effects were "hand-tuned" (see Appendix 4 for an explanation). The parameters had to be reduced (from −2.57) most for C. decurrens and A. concolor, which have very small seeds and first-year seedlings, and the least for P. jeffreyi and P. lambertiana, which have large seeds and robust seedlings. All species had a maximum local adult BA beyond which we never observed seedlings. Therefore, if local adult BA exceeds the threshold (Appendix 1), the probability of seedling survival is set to 0.

Reproduction
In order to reduce computation time, we implicitly combine seed dispersal and germination, such that the "seeds" actually represent first-year seedlings. The fecundity and dispersal functions are parameterized based on the distribution of firstyear seedlings relative to adult trees (Appendix 3), and germination probability is set to 1. Fecundity is simulated with a zero-inflated Poisson (ZIP) distribution: some individuals produce zero offspring (with probability p), while others produce offspring according to a Poisson distribution with parameter λ. Where z j,t is the number of new seedlings from tree j in year t: Many of our tree species exhibit mastinglocally synchronized mass seed production (Koenig et al. 2015;Gallego Zamorano et al. 2018). We therefore used two ZIP distributions for masting and non-masting years respectively (Appendix 3). The probability of masting is based on time since the last mast: where y is the probability of masting, X is years since last mast, and a and b are fitted parameters. We modeled dispersal using a two-dimensional Student's (2Dt) distribution (Clark et al. 1999a): first-year seedlings found in quadrat i of area A in year t; d i, j the distance between i and tree j; J the total potential parent trees; z j, t the number of first-year seedlings produced by j; and u the shape parameter of the 2D-t function. Because first-year seedlings are difficult to identify to species, the dispersal and fecundity model was fit at the genus level (Appendix 3).

Hindcasting
We ran 20 simulations using historical climate sequences for three 1 ha USGS plots -POFLABMA (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), and BBBPIPO andUPLOG (1997-2013) obtained from the California Basin Characterization Model (CA BCM) model, which downscales 800-km-scale global climate data to 270 m (Flint et al. 2013a). Initial plot maps used coordinates and 1997 or 1999 DBH of adult trees in these plots replicated 9 times in a 3 × 3 array (Appendix 1, Fig. 7). The buffer "ring" of replicated plots compensates for edge effects. To create a seedling map, we calculated the 1999 density within seedling subplots of classes 2 and 3, which are less affected by year-to-year variation, and replicated this density across the whole plot. We then compared the tree density and basal area for each species in the central square to measured values. We also compared species-level averages for individual growth and mortality. For this, we used only the trees initially present in the plot, so that the youngest trees would not skew the comparison.

Simulating climate change responses
Once the model appeared to be capturing the past behavior, we projected forest dynamics for seven plots: the three hindcasting plots plus four new 0.8 ha plots established in Sequoia National Park in 2015 (Table 2). Annual projected climate sequences for 2006-2099 were downloaded from CA BCM for three global climate models (GCMs) and one emissions pathway: RCP 8.5, a "worst-case scenario" that would be expected to produce the strongest effects. While this pre-Paris Climate Accord pathway may be avoided with the initial Nationally Determined Contributions, warming of 2.5-4°C by 2100 (between RCP 4.5 and 6) would still occur (IPCC 2013; Rogelj et al. 2016); several countries, including the USA, had not met those NDCs as of 2020 (Roelfsema et al. 2020). Different GCMs produce different outcomes from RCP 8.5: MIROC represents a relatively hot-dry scenario; CCSM a moderate-warm scenario; and CNRM a relatively warm-wet scenario (Appendix 1; Fig. 8). For each plot, a control climate sequence ("current") was created from historical 1974-2014 sequences randomly sampled for the same number of time steps. We ran ten replicates for each site-scenario combination.

Hindcasting
In nearly all cases, the observed mean species-level adult growth or mortality was within the 95% CI for the simulations and vice versa (Table 3). Where there were exceptions, the simulated mean was within the observed range, but the observed mean was higher than the simulated range. Measured total BA and density fell within or close to the 95% CI of 20 model runs for most species in most plots. However, A. concolor BA and density were somewhat overestimated in UPLOG, and C. decurrens BA in BBBPIPO and A. magnifica density in POFLABMA were underestimated (Fig. 1).

Climate change responses
Initial tree densities across plots ranged from 387.5 to 1189 trees/ha and total BA from 52.46 to 108.28 m 2 /ha. There was no correlation between initial tree density and initial BA, and while the lowest elevation plot (BBBPIPO) had an unusually high tree density, there was otherwise no relationship between tree density and elevation. There was, however, a trend toward higher initial tree BA at higher elevations (p = 0.17, adjusted R 2 = 0.20; Appendix 5). Total tree BA and stem density increased over the 93year simulations in all sites except the very dense lowest elevation site, where stem density decreased by about 20% in all scenarios (Figs. 2 & 3 and Appendix 5). Density increases in the other sites ranged from 41.8 to 527.6% and tended to be related negatively to initial density and positively to elevation, though this was not statistically significant (Appendix 5). Basal area increases ranged from 13.2 to 62.7% and also tended to be positively related to elevation (Appendix 5). The "warm-wet" scenario (CNRM) resulted in the greatest BA increase and the lowest drop in density at BBBPIPO, as well as the highest increase in BA and density at SP. Response to climate scenarios varied substantially by species. Pinus monticola, P. concolor, and Q. chrysolepis are omitted from the following species-level discussion due to low abundance. Of the two fir species, A. magnifica always increased in basal area and density, while A. concolor lost BA under climate change in sites SJP and SJM (Figs. 2, 3, 4). A. concolor also exhibited less BA growth under climate change relative to the control at SP and SJ, and reduced density increases at multiple sites. Abies magnifica, on the other hand, exhibited similar density changes across scenarios and a greater increase in BA under climate change in higher elevation plots.
Increases in P. lambertiana BA and density were reduced by warming in the two lowest elevation plots but increased under moister scenarios at SJ. At SJP and UPLOG, P. lambertiana BA declined but density increased substantially; with warming, the decline in BA was more pronounced and the increase in density smaller. P. ponderosa declined in BA and increased in density for all site-scenario combinations except for SP-Control, SJP-Control, and SJP-MIROC ("hot-dry") where BA increased (Figs. 2,3,4,5,6). Pinus jeffreyi increased in both density and BA at SJP and UPLOG, more so in the climate change scenarios, while at SJ they only increased substantially with climate change (Figs. 2, 3, 4).
Calocedrus decurrens decreased in density in all scenarios while usually increasing in BA (Figs. 2 & 3). The only site-scenario combination that decreased in cedar BA was SP-CCSM ("moderate-warm") but, as there are very few cedars at this site, that could have been due to stochasticity. Warming favored increased BA of Q. kelloggii (Figs. 2 & 5), while stem density always decreased, but less so with warming at SP and UPLOG (Figs. 3 & 6).

Discussion
Our results suggest that, in the absence of disturbance, forests in the southern Sierra Nevada would likely increase in stem density and basal area over the coming decades (Figs. 2 & 3), particularly at sites with low current density or high BA. The direction of change in total density and BA was the same for the three RCP 8.5 scenarios and 1974-2014 "control" at all sites, as well as for most site-species combinations. Therefore, we chose not to run scenarios RCP 6.0 or 4.5, as these would almost certainly fall within the same range. For instance, total BA would still likely increase a similar amount across all sites, while C. decurrens BA at UPLOG would increase more than the control but less than the RCP 8.5 scenarios. Site-level variation was much higher than between scenarios. The stem density decrease at the lowest site, BBBPIPO, was likely due to higher mortality among the smallest stems, as average pertree BA increased at this site across all scenarios despite modest total BA increases, a pattern that was also seen to a lesser extent in site SJ; increased recruitment decreased average tree BA at the other sites (Fig. 14, Appendix 5). Species-level responses varied more by climate scenario, though site-level variation was still considerable; changes were often consistent in direction between scenarios within a site.
While firs continued to dominate, A. concolor performed worse under warmer conditions at 4 out of 6 sites, while A. magnifica performed as well or better. This is consistent with the individual demographic behaviors: Higher temperatures favored more positive demographic rates in A. magnifica  (Barbour et al. 1991). Firs also tolerate shading, which was likely an advantage in the increasingly dense simulated stands. Mean A. concolor BA declined across all scenarios in 5 out of 6 sites due to some combination of reduced growth and increased mortality relative to recruitment, while treelevel A. magnifica BA increased under climate change scenarios more than the control in 3 out of 4 sites despite little change in density, reflecting increased growth (Fig. 14,  Appendix 5). The biggest winner in the climate change scenarios was P. jeffreyi, which exhibited both BA and density increases, though this species also increased at 2 out of 3 sites in the control scenario. The increases in BA were driven primarily by recruitment, leading to steady or declining mean tree BA (Fig. 14, Appendix 5). The positive effects of warm dry conditions on P. jeffreyi are consistent with previously observed  (Table 1) as well as its abundance on the dry eastern side of the Sierra Nevada and broader elevational range in the south (Haller 1959). The responses of P. lambertiana were complex, showing increases in density at all sites across climate scenarios but site-specific differences in BA changes. Increases in BA and density at low elevations were reduced with climate change at 1609-1806 m but increased in 2 of the 3 warming scenarios at 2197.6 m. At SJP and UPLOG, however, BA decreased by up to 28.9%, and density increased by as much as 785%, indicating large tree mortality across climate scenarios, which led to a steep decline in mean tree basal area (Appendix 5). Increased winter temperatures boost the growth of adult P. lambertiana, but higher summer temperatures negatively affect its survival (Table 1). Pinus lambertiana appears to have unusually high shade-tolerance for a pine (Moran et al. 2019), which may have contributed to its density increases.
Both C. decurrens and Q. kelloggii are predicted to decrease in density but increase in BA. This was driven by higher mortality and/or lower recruitment of small individuals, leading to increased average tree size (Fig. 14, Appendix 5). Warmer temperatures favor adult growth and survival in Q. kelloggii and adult survival in C. decurrens but negatively impact their seedlings (Table 1). Warmer conditions particularly favor increases in total Quercus BA; for C. decurrens, this was only consistently true at the highest site of occurrence (Fig. 2). Both species increased in density during the twentieth century (Dolanc et al. 2014a;McIntyre et al. 2015) and exhibited low mortality in the 2012-2016 "hot drought" (Fettig et al. 2019). Conversely, P. ponderosa usually decreased in BA but increased in density; higher temperatures negatively impact P. ponderosa survival but boost seedling growth (Table 1). The increase in P. ponderosa density was surprising given that this species is fairly shade-intolerant, but the sites containing it  Increases in forest density and declines in large trees occurred across California during the twentieth century (McIntyre et al. 2015), but this trend was most pronounced in private timberlands, and densification was greater on nonwilderness National Forest lands than in National Parks (Easterday et al. 2018). Both our USGS and newly established plots have had no management interventions other than, in some cases, prescribed fire, since at least the 1980s (Appendix 2). Trajectories for more disturbed sites might differ, so we would encourage those who are interested in National Forest or private timberlands to apply the parameterized model to their specific sites. Climate projections for various GCM-RCP combinations can be downloaded from the California Climate Commons (http://climate.calcommons.org/bcm).

Fig. 6
Tree density projections, site SP. Left -Q. kelloggii (QUKE, top) and P. ponderosa (PIPO, bottom) in all 4 climate scenarios: Control, MIROC ("hot-dry"), CCSM4 ("moderatewarm"), and CNRM ("warmwet"). Right -All species in control (top) and "hot-dry" (bottom) scenarios. Solid lines -median over 10 replicates. Shading -95% confidence interval. ABCO, A. concolor; CADE, C. decurrens; PILA, P. lambertiana; PIPO, P. ponderosa; QUCH, Q. chrysolepis; QUKE, Q. kelloggii Some of the patterns we observed, such as the increased basal area of Q. kelloggii under climate change (Fig. 5), appeared after the first 30 years. An analysis based on tree growth increment data was more consistent with our findings, projecting increased growth in California forests under climate change both with and without CO 2 -induced increases in water use efficiency (Charney et al. 2016). Prior landscape-scale vegetation models also predict that climate change favors the recruitment of more drought-tolerant species such as oaks and gray pines (Lenihan et al. 2003;Liang et al. 2017).
Our projections of increased overall forest density would likely change if fire effects were included, but the tendency toward density increases could make fuel management more challenging. Fire suppression since the 1880s has already led to increasing dominance of firs, higher forest density, and increased risk of severe fires (Beaty and Taylor 2007;Collins et al. 2011;Earles et al. 2014). In the Sierra Nevada, forest area burned per year increases with spring and summer temperatures (Keeley and Syphard 2016). A recent simulation study suggested that thinning and understory burning, while having little effect under contemporary fire regimes, could significantly decrease fire severity, increase carbon sequestration stability, and maintain net ecosystem exchange under projected extreme fire weather (Krofcheck et al. 2017). The Liang et al. model (2017) included climate-associated increases in wildfire and predicted a reduction in landscape C sequestration potential. The impacts of both low and high intensity fire should be included in future analyses if projections are to be used to guide management. SORTIE-ND includes several fire behaviors; however, our aim here was to examine direct climate effects.
Another source of mortality that can be influenced by climate is disease and pest outbreaks. In the 2012-2016 drought, various species of bark beetles were the proximate cause of death for many trees, particularly large pines (Fettig et al. 2019). The transmission and effects of fungal pathogens, such as the blister rust that affects sugar and western white pines, on the other hand, are favored by wet conditions (Tomback and Achuff 2010). However, these effects are challenging to include in forest simulators, as host and herbivore/pathogen dynamics interact.
Early seedling survival proved to have a strong influence on adult and sapling density. Parameterizing these dynamics can be challenging, as very young seedlings are difficult to detect and individually tag. We plan to address this in future iterations of the model. Including soil properties might also be important for local dynamics, as the model does not currently account for rock outcrops, soil depth, or low nutrients that might limit tree establishment and growth. However, such data are not always recorded for long-term forest monitoring plots.
In interpreting both total and species-level responses, it is important to remember that models are parameterized based on observed variation. For the 26 parameterization plots, maximum July temperatures increased 1-3°C since the 1970s, but changes in January minimums were smaller and average precipitation did not change (Moran et al. 2019). Growth sensitivity to climate has been observed to shift over time (Wilmking et al. 2020); this may also be true of mortality and fecundity responses. The same is true of emissions scenarios, which are subject to change due to human policies and behaviors. Therefore, projections beyond 50-100 years should be considered indicative of direction and magnitude of possible changes but not precise predictions.

Conclusion
This is the first time an individual-based simulator has been applied to California Mediterranean-climate forests. It is encouraging that, despite increasing temperaturedriven aridity, the model does not project a collapse of tree diversity or basal area due to climate change alone over the next century. Indeed, while shifts in forest composition favoring more heat-and drought-tolerant species are likely, tree densities and basal areas are projected to increase in the absence of disturbance, particularly above 2000-m elevation. However, these results should not be taken as reason to be complacent, as increasing density plus the drying effects of higher summer temperature could contribute to increased risk of severe wildfire. Our results also suggest a need to better understand tree reproduction, as changes in tree density were quite sensitive to these processes and seed production can be both directly and indirectly influenced by climate (Clark et al. 2021).

Appendix 1. Model choice and setup
Why SORTIE-ND?
As cited in the main text, several forest models have been parameterized for California forests. However, most are not individual-level models, which is what was required for our future plans of investigating the impact of individual-level variation and heritability in climate responses on forest dynamics under climate change. Liang et al. (2017) used the landscape-scale model LANDIS II, in which species are represented by biomass in age cohorts. Lenihan et al. (2003) used MAPSS-CENTURY 1, a dynamic vegetation model that simulates vegetation types in grid cells of > 900 m 2 over very large landscapes. CACTOS (Battles et al. 2008) does track individuals but was designed to capture short-term stand growth and includes only simple mortality and in-growth functions. Charney et al. (2016) used an entirely different approach, correlating growth rings with climate, grouping sites according to their climate responses, calculating vulnerability indexes based on projected climate change and local growth responsiveness, and using this to forecast changes in growth under different scenarios.
We considered various other individual-based modeling approaches (Bugmann 2001), including FORCLIM (Bugmann 1996), FORSKA (Prentice et al. 1993), JABOWA (Kienast 1991), and PPA (Purves 2008). However, FORCLIM uses cohorts to describe trees of the same species and age rather than fully treating them as individuals; the source code was also difficult to find. FORSKA does not track individuals under 1-cm DBH. The most recent JABOWA code was not available without purchase, and the model does not track seedlings and treats each canopy as covering an entire patch. PPA is a SORTIE descendant that includes physiological responses but does not track demography on as fine a scale. SORTIE-ND offered high accessibility with the ability to download old model behaviors and upload novel behaviors for other users to access; a history of success with applications to multiple forest types; and pre-existing behaviors that were relatively simple to modify for the type of forest demography data we had available. Moreover, the Battle lab at UC Berkeley is developing SORTIE-NDcompatible tree-shrub competition behaviors (personal communication) that could be combined with our future studies.

Climate sequences
As mentioned in the main text, we focused on January minimum temperature, July maximum temperature, and annual precipitation. These variables were chosen for consistency with the prior published statistical analyses used to parameterize the seedling growth and survival and adult growth behaviors. While the absolute minimum/ maximum temperatures do not always occur in January or July, these monthly minimums are closely correlated with overall winter and summer minima and maxima. Examples of changes in these values under different climate change scenarios are shown in the figure below. Fig. 7 Landscape setup. Darker green square is the focal area, lighter green area the buffer region. Each square is initiated with a replicate of the real tree map Fig. 8 Projected future climate sequences from 2006 to 2099 for three hindcast plots based on MIROC ("hot-dry"), CCSM ("moderate-warm"), and CNRM ("warm-wet") GCMs for RCP 8.5. Figure also includes the "current" or control sequence drawn from the observed historical sequence. July maximum temperature (°C), precipitation (mm), and January minimum temperature (°C)  We used a Bernoulli process to model the mortality process, and estimated the probability of dying with a linear model and a logit link function:

Plot setup
where Y i,t is the state of the tree i at time t ; p i,t is the associated probability of dying; the β's are a vector of parameters to estimate; and ε is an error term. Both DBH and DBH 2 are included as predictors because both very small and very large trees tend to exhibit elevated mortality. Competition is also known as an important predictor of mortality because it decreases the access to resources. We used a simple competition index: where N i is the number of neighbors of tree i, D j,t is the diameter of neighbor j in year t, and dist i,j is the distance between trees i and j.
To estimate the parameters of the model, we used a Bayesian calibration implemented in R (R core team 2017) using Gibbs sampling. The vector of β's is assigned a multivariate normal prior, and parameter σ 2 an inverse gamma prior. The model was fit for each species individually. To select the climate variables included in the model, we ran the model with pairs of variables: one related to temperature (JulMax or JanMin) and one related to water availability (precipitation, actual evapotranspiration (AET), or April snowpack). The combination of variables used in the final model is selected based on DIC values. Years are considered masting years if the number of seedlings is higher than the mean number of seedlings. The vector of zeros and ones corresponding to masting or non-masting years is called mast thereafter. We modeled dispersal using a 2Dt distribution function (Clark et al., 1999a). The number of seedlings of genus g found in quadrat i (area 5 × 5 m, or 25 m 2 ) in year t is: where d i,j is the distance between the center of quadrat i and the tree j; J g is the number of trees of genus g in the plot; z j,t is the number of seedlings produced by tree j in year t, and u is the parameter of the 2Dt function that determines its shape.
We estimated u across all plots, but assumed that trees in one plot do not disperse their seeds far enough to add seedlings in the quadrats of another. The latent variable z j,t must also be estimated. For u, the parameter of the 2Dt function for dispersion, we choose the same weak prior as in Clark et al., (1999b): a Gamma distribution with parameters (1, 0.01). For p m and p n , the first parameters of the zero-inflated Poisson (ZIP) models for masting and non-masting years respectively; the priors were uniform distributions between 0 and 1, as it is a probability. For λ m and λ n the second parameters of the ZIP for masting and non-masting years respectively, a Gamma distribution with parameters (2, 0.01).
We used a Gibbs sampler to estimate the parameters and the latent variables (masting years and fecundity for each tree and each year), and a Metropolis-Hastings step at each step of the Gibbs sampler to sample from full conditional distribution. At each time step, the algorithm estimates first u; then p m and p n ; then λ m and λ n ; then z j,t , the true number of seedlings produced by each tree. Results are presented in Table 6.
The results suggest that pines have the longest effective dispersal distances (median u equivalent to mean dispersal distance of 126.12 m), while cedars disperse the shortest distances (median u equivalent to mean dispersal distance of 36.35 m), and firs and oaks are intermediate. P in both masting and nonmasting years is bigger for pines and oaks, indicating a greater variation between individuals in seed output than for firs and cedars. λ m is the largest for firs and cedars, which on average produce over 45 new seedlings per tree in a mast year, while the larger-seeded oaks and pines produce 1-6 new seedlings. All species produce on average less than one seedling per tree in a non-mast year, but λ n is slightly larger for pines and oaks.  Appendix 5. Additional graphs and tables of initial conditions and changes in density and basal area Fig. 11 Initial conditions for the seven sites. Initial tree density vs. initial basal area, initial basal area vs. elevation, and initial tree density vs. elevation. In the latter case, the relationship is not quite significant (p = 0.14) but the adjusted R 2 = 0.25