Fine-scale modeling of bristlecone pine treeline position in the Great Basin, USA

Great Basin bristlecone pine (Pinus longaeva) and foxtail pine (Pinus balfouriana) are valuable paleoclimate resources due to their longevity and climatic sensitivity of their annually-resolved rings. Treeline research has shown that growing season temperatures limit tree growth at and just below the upper treeline. In the Great Basin, the presence of precisely dated remnant wood above modern treeline shows that the treeline ecotone shifts at centennial timescales tracking long-term changes in climate; in some areas during the Holocene climatic optimum treeline was 100 meters higher than at present. Regional treeline position models built exclusively from climate data may identify characteristics specific to Great Basin treelines and inform future physiological studies, providing a measure of climate sensitivity specific to bristlecone and foxtail pine treelines. This study implements a topoclimatic analysis—using topographic variables to explain patterns in surface temperatures across diverse mountainous terrain—to model the treeline position of three semi-arid bristlecone and/or foxtail pine treelines in the Great Basin as a function of growing season length and mean temperature calculated from in situ measurements. Results indicate: (1) the treeline sites used in this study are similar to other treelines globally, and require a growing season length of between 147–153 days and average temperature ranging from 5.5°C–7.2°C, (2) site-specific treeline position models may be improved through topoclimatic analysis and (3) treeline position in the Great Basin is likely out of equilibrium with the current climate, indicating a possible future upslope shift in treeline position.


Introduction
The treeline ecotone on a mountain is the transition zone between closed montane forest and treeless alpine landscape, encompassing the highest locations where mature trees are found (Wardle 1971, Scuderi 1987, Jobbagy and Jackson 2000, Körner 2012). In the absence of disturbance-related conditions and substrate prohibiting tree growth, the treeline position represents a boundary between areas in which climatic conditions allow for physiological activity in mature trees, and areas where tree growth is not possible. Research suggests this life-form boundary is climatelimited; regardless of species, elevation, or latitude, treeline positions globally share common climatological characteristics (Wardle 1971, Jobbagy and Jackson 2000, Körner 2012, Weiss et al 2015. Two independent studies Paulsen 2004, Paulsen andKörner 2014) provide evidence of a common growing-season isotherm around 5°C-6°C present at many different treeline sites globally.
Accordingly, climate-limited treelines are valued as paleoclimatic indicators of environmental change as regional treeline positions have been shown to track centennial-scale changes in climatic conditions (Scuderi 1987, Lloyd and Graumlich 1997, Salzer et al 2013. In the American southwest, Great Basin bristlecone pine (Pinus longaeva, D. K. Bailey) forms climate-limited treelines throughout Nevada and California. This species is a valuable climate proxy due to its extremely long-lived nature (e.g. Currey 1965) and the tendency of its annual rings to correlate with the most growthlimiting environmental factor. Ring-width chronologies from the upper-forest-border (at and just below Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. treeline) have been widely used as a proxy for temperature (e.g. LaMarche Jr 1974), while ring-width chronologies from the more arid lower-forest-border have been used a proxy for summer precipitation (e.g. Hughes and Funkhouser 2003). These findings indicate the primary growth-limiting factor operates on a gradient, changing from moisture limitation at the lower-forest-border to temperature limitation at the upper-forest-border (Kipfmueller and Salzer 2010).
Past research has shown topography influences climate-and subsequently biological systems-on the scale of tens to hundreds of meters (Weiss et al 1988, Lookingbill and Urban 2003, Dobrowski et al 2009, Geiger et al 2009, Adams et al 2014. This phenomena is referred to as topoclimate, and has been the subject of our recent research regarding the climate response of near-treeline bristlecone pine (Bunn et al 2011, Salzer et al 2013. Bunn et al (2011) discovered that topographic position affects the growth response of trees; individual trees growing well below the upperforest-border in areas of cold air pooling displayed distinctly different ring-width patterns from nearby trees (within tens of meters) outside areas of cold air pooling. Further, the climate signal of low-elevation trees in areas of cold air pooling was very similar to the classic temperature-limited signal characteristic of the upper-forest-border. Salzer et al (2014) built on Bunn et al (2011) by constructing treeline and below-treeline chronologies from north and south-facing aspects. The authors identified a divergence in growth patterns between north and south facing aspects, as well as a climate-response-threshold between moisture and temperature limitation approximately 60-80 vertical meters below treeline.
This study models treeline positions from a topoclimatic perspective. Combining evidence of climate-driven treeline formation with in situ temperature measurements, we present three site-specific models in the Great Basin predicting bristlecone pine treeline position as a function of topoclimate.

Study areas
We chose three Great Basin treeline sites for this analysis (figure 1); (1) Mount Washington, Snake Mountain Range, NV (MWA, 38.91°N. lat., 114.31°W. long., treeline position approximately 3400 m.a.s.l.), (2) Chicken Spring Lake, Sierra Nevada, CA (CSL, 36.46°N. lat., 118.23°W. long., treeline position approximately. 3600 m.a.s.l.), and (3) Sheep Mountain, White Mountains, CA (SHP, 37.52°N. lat., 118.20°W. long., treeline position approximately. 3500 m.a.s.l.). Sites MWA and SHP support Great Basin bristlecone pine treelines, while the CSL treeline is formed by mostly foxtail pine (Pinus balfournaia, Grev. & Balf.), a closely related species to bristlecone pine with a slightly shorter life-span and similar climate growth-response (Lloyd and Graumlich 1997). Environ. Res. Lett. 12 (2017) 014008 meter in living trees, dispersed across varying topographic features within a 1 km 2 -2 km 2 area. Our primary goal was to capture differences in temperature between different topographic positions, so the relative differences in temperature between all the sensors at a given site were equally as important to the raw recorded temperatures. Because we recorded temperatures for only one calendar year at any given site, our data reflect the weather conditions at that site specific to the period of deployment, rather than a long-term climatic average (figure 2). To more accurately represent the average climate at each treeline location, we calculated monthly anomalies between the temperature during deployment and the climate normal for each location, and applied these corrections to our sensor data (PRISM Climate Group 2004). This provided a data set that captured relative differences in temperature due to topography, which the raw values representative of the average climate, rather than anomalous weather during the period of deployment (figure 2). We then applied a warming correction to our data set to more accurately represent the climate when Great Basin treelines stabilized their current positions. Salzer et al (2013) report treeline positions in the Great Basin moved downslope up to 100 meters below their highest positions during the Holocene climatic optimum, and established their current positions (well below the maximum positions during the climatic optimum) in the early 1300s A.D. (also see Carrara and McGeehin 2015). The authors present a multi-millennial Great Basin climate reconstruction from bristlecone pine chronologies of previous September-August temperature anomalies relative to a baseline period A.D. , which shows an approximate warming of 1.5°C between the period when treeline positions in the Great Basin stabilized in the early 1300's and present day. Therefore, we subtracted 1.5°C from our observed temperatures so that the topoclimate dataset would most accurately represent the climate that influenced the current treeline positions when they established in the early 1300's, rather than today's climate that has no influence on treeline positions formed in the past.

Topoclimate analysis
From the observed hourly temperatures, we calculated values of two climate variables unique to each sensor: average monthly temperatures were calculated by averaging all hourly values within a given month, yielding twelve values per sensor; annual sum of degree hours above 5°C was also calculated, yielding one value per sensor. We used lasso regression models (Kuhn 2015) (10-fold cross-validation, ten repeats per fold) to model each climate variable as a function of topographic variables at ten meter resolution. The topographic variables used for prediction are: elevation, slope, aspect-derived Eastness and Southness indices, topographic position and convergence indices, and solar radiation loads. The models were used to predict the variables across areas above 3000 m.a.s.l. at each site, yielding thirteen topoclimate raster surfaces for each study location representing values of average monthly temperature and degree hours above 5°C. Model skill was relatively high but fluctuated between variables, and relied most      Environ. Res. Lett. 12 (2017) 014008 on elevation and solar radiation values as predictors (see appendix A for more on this process and measures of model skill).

Treeline position models
Paulsen and Körner (2014) present a model that predicts treeline positions globally as a function of three parameters; a threshold temperature (DTMIN, measured in°C) above which physiological activity is possible, a growing season length (LGS, measured in days) that includes all days with an average daily temperature above DTMIN, and a seasonal mean temperature (SMT, measured in°C) that is the average of all days within the growing season. Using the authors' best fit value of DTMIN (0.9°C), we adopted their methods to calculate LGS and SMT raster surfaces at each site from our predicted monthly topoclimate surfaces. We used cubic splines to interpolate daily temperatures from the modeled monthly topoclimate rasters, and summed the number of days with average temperatures above 0.9°C for the growing season length, and averaged the daily temperatures of all days within the growing season to find the seasonal mean temperature. We built classification models using the LGS and SMT raster variables to predict treeline position as the boundary between two mountainous biomes; a subalpine region of closed montane forest, and a treeless alpine region above the upper-forest-border ( figure 3 panel (a)). We defined the boundaries of each biome around the treeline position through multi-step process: (1) Using Google Earth we digitized treeline position at the landscape scale (the red line in figure 3 panel (a)). Conventions set by Körner (2007Körner ( , 2012 define treeline position at a larger scale by connecting straight lines between the upper reaches of mature trees. We altered this method because our 10 meter resolution topoclimate variables allowed for a more resolved definition of treeline position. We were very deliberate in the areas of treeline used to build the models, selecting only stretches of treeline that were obviously climate-limited, and not influenced by disturbances such as slope, rockfall, lack of substrate, etc. (2) We then set a 25 meter upslope and downslope buffer for the boundary of each biome nearest to treeline, to ensure a conservative separation between the upper boundary of the subalpine and the lower boundary of the alpine, and set the width of each biome to 250 meters.
With the biome regions delineated, we obtained training data for the classification models by extracting values of LGS and SMT specific to each biome from randomly spaced points with a density of 500 points per square kilometer. Classification models were then developed through an iterative process at each site; we generated three models with maximum branch lengths of one, two, and three splits, and compared the accuracy, complexity, and cost of adding additional splits between each model. The simplest, most accurate model was chosen by balancing the prediction accuracy and the complexity of each model, with the fewest number of splits and terminal nodes representing the simplest model. For example, if the prediction accuracy was similar between models of different complexities (one split vs two or three splits), preference was given to the model with the fewest number of splits.

Results and discussion
3.1. Treeline prediction The classification trees (figure 4) at all sites suggest seasonal mean temperature is the best predictor of Environ. Res. Lett. 12 (2017) 014008 treeline position. SMT thresholds of 7.2°C at MWA (89% overall accuracy), 5.4°C at CSL (81% overall accuracy), 6.2°C at SHP (72% overall accuracy) separate the alpine and subalpine regions. When secondary and tertiary splits were allowed the overall accuracy increased by 0.4%, 0.6%, and 0.1% respectively, thus we present the simplest model at each site with a single threshold in SMT separating the biomes. The Snake Range (MWA) model was the most accurate (table 1 and figure 5). Rates of misclassification were different between the predicted biomes (producer accuracy: alpine 85%, subalpine 92%) indicating a model slightly biased towards subalpine prediction. Cohen's kappa statistic-a measure of how different a prediction is from a randomized classification-is 0.78, indicating this classification is different than random with substantial agreement. The Sierra Nevada Range (CSL) model was slightly less accurate than the MWA model (table 1 and figure 6), with a slightly larger bias towards subalpine prediction (producer accuracy: alpine 66%, subalpine 93%). Cohen's kappa indicates substantial agreement the prediction is different than random. The White Mountains (SHP) model provided the least accurate prediction (producer accuracy: alpine 63%, subalpine 82%), yet was consistent with other sites in a bias toward subalpine prediction, and Cohen's kappa indicates moderate agreement the prediction is different than random. It is important to note that at all sites higher rates of producer accuracy for the subalpine region come at the cost of overpredicting this region. This results in a tendency to predict the treeline position slightly higher than a model without this bias, and is most pronounced at CSL and SHP.
A drawback to this method lies in the discontinuity of easily identifiable climate limited treelines. In the absence of other unrelated factors that influence treeline position, the entire treeline at a site would be climate-limited. However this is rarely the case, at each site there are many stretches of treeline that are clearly driven by factors other than climate-limitation, such as substrate availability, slope, or snowpack. The challenge lies in identifying enough continuous area that display obvious climate limitation to allow a treeline position model. As a result, appropriate stretches of treeline for modeling were usually limited to a single aspect due to the geographies of each site and locations of climate limited treelines: at MWA the climate-limited treeline faces west, at CSL the climate limited foxtail pine treeline faces southwest; and at SHP the climate limited treeline faces east. Topoclimate modeling requires a dense network of sensors in a single area, and thus we needed to locate our sensors on a single side of each mountain range near the climate limited treelines to maximize out predictive power (Bruening 2016). While this method allows for accurate climate prediction near the sensors, the models lose predictive power on the other sides of the mountain ranges. Consequently, these models are inherently site-specific and are in no way intended as a comprehensive treeline position model for the Great Basin. This analysis is different from other treeline models conducted at larger scales, yet our results regarding the physiological constraints of treeline position at these sites are consistent with the global models.

Global model comparisons
This analysis provides comparisons to previous treeline studies while accounting for topocli-matic effects on treeline position (Körner and Paulsen 2004, Körner 2012, Paulsen and Koorner 2014. We model the length of the growing season and its average temperature at a scale previously unavailable to treeline researchers, which have provided insights into the physiology of near-treeline bristlecone pine growth in the Great Basin. In a related analysis, Tran et al (2017) examined the climate-growth response of bristlecone pine at MWA, CSL, and SHP, and used cluster analysis to identify the primary growth patterns in the ring-widths of trees sampled at these sites. They found that Bristlecone pine growth nearest to treeline was controlled by SMT (calculated via the same methods as described in this analysis), while Bristlecone pine growth farther downslope and away from the upper-forest-border was more influenced by    2. Note: the authors report that the LGS value is at least 94 days, and may be longer depending on the climate region). These values are not specific to any one species or ecozone, and are mostly consistent with our topoclimate calculations. The similarity in seasonal mean temperatures suggests that Great Basin treelines are likely subject to a similar SMT isotherm as other sites globally, despite differences in growing season length.
According to Körner (2012) the length of the growing season at treeline varies significantly between climate regions, shown by plots of daily mean root zone temperature from 32 different treeline sites from various different climate regions (pages 40-47). Each site has an estimated length of the growing season calculated from root-zone and air temperatures, with values spanning from around 100 days in the subarctic and boreal zones up to 365 days in the equatorial tropics. Growing season lengths from warm-temperate and cool-temperate climate zones (the Great Basin falling somewhere in between these climate zones) range between 122-190 days, consistent with our modeled values of LGS. For further validation of our modeled treeline variables, we obtained observed root zone temperatures on Mount Washington (MWA) which allowed for independent calculations of SMT and LGS (Scotty Strachan, unpublished data). We calculated LGS to be 152 days (modeled value is 147 days, table 2) and SMT to be 7.7°C (modeled value is 7.2°C, table 2).
While the length of the growing season and seasonal mean temperature may be useful predictors of treeline (as defined and calculated by Körner (2012), Paulsen and Körner (2014) and in this analysis), these variables explain little about tree physiology. Two treeline sites may have comparable season lengths and average temperatures with contrasting annual mean temperature profiles (Körner 2012). Paulsen and Körner (2014) guard against interpreting these variables too literally-TREELIM calculates a best fit of 6.4°C as 'a common isotherm of low temperature for forest limits' , and they stress the absence of physiological relevance represented by this value. While it is close to the physiological limit for woody biomass accumulation, they argue the seasonal mean temperature 'reflects an arithmetic mean that is subsuming the combined action of low temperature, integrated over time, on a suite of processes associated with tissue formation, from root tips to apical meristems' .
3.3. Treeline position models, treeline potential, and growth-limiting factors Modeled treeline positions were obtained for each mountain range by applying the SMT threshold calculated in each classification model to the SMT raster surface at each site (figures 5-7). We developed modeled treeline positions under two different climate scenarios.
The first ((a) in figures 5-7) represents predictions of the current treeline position using SMT representative of the climate from the early 1300s (the period of current treeline stabilization, see section 2.2), the SMT raster used to build the models. The second ((b) in figures 5-7) represents the same models applied to today's climate, which has warmed 1.5°C since the treeline positions in the Great Basin stabilized their current position according to Salzer et al (2013).
The differences between the treeline position model projections (figure 4) are the result of a relatively sensitive threshold between the subalpine and alpine regions. By projecting our model using climate data representative of the current climate, we conclude treeline positions at our sites are likely out of equilibrium with the Great Basin climate; at all sites the treeline position modeled from today's climate ((b) in figures 5-7) moved upslope from its current position ((a) in figures 5-7). It is important to note that b is a projection of treeline position potential based on warming at these sites since the current treelines established their positions. We speculate that current treeline positions have yet to 'catch up' with the the current climate, a result of the slow nature of treeline dynamics (Körner 2012). The demographic processes that cause an upslope shift in treeline position are slow and lag behind changes in climate (Scuderi 1987, 1994, Lloyd and Graumlich 1997, however an in-depth discussion of how and why this lag exists is outside the scope of this analysis. It is unknown exactly how long behind changes in climate treeline position lags, but it is speculated to be as long as hundreds of years depending on the ecozone, species, slope, etc. As the Great Basin climate gradually warms, uninhabitable areas above treeline start to become Table 2. Values of LGS and SMT from from the treeline positions at each site, compared to a global model's best fit presented by Paulsen and Korner (2014 (2007). The authors provide Bristlecone pine range projections indicating a single degree of warming could be enough to initiate an upslope migration of Bristlecone pine treeline by tens to hundreds of meters. This is supported by our topoclimate analysis and the sensitivity of the treeline position models ( figure 4). An increase in seasonal mean temperature would force the position of the SMT isotherm upslope and make conditions more favorable for Bristlecone pine growth above the current treeline position ((b) in figures 5 to 7). The extent to which the SMT isotherm will move upslope is dependent on the topoclimatology of each treeline position; for example a gradual slope would enable farther upslope migration of treeline position than a steep slope. An in depth analysis of microrefugia and topoclimatology in the Great Basin would benefit the study of bristlecone pine treeline dynamics.

Conclusion
We predicted the position of three Great Basin treelines formed by bristlecone and foxtail pine exclusively from climate and topography data. Through a topoclimatic analysis we captured landscape-scale effects of topography on climateand consequently on treeline position-that are independent of changes in elevation. Our results indicate the average temperature throughout the growing season (SMT) is the most dominant factor influencing treeline position on the landscape, regardless of species or elevation, in agreement with previous research Jobbagy and Jackson (2000), Körner (2012), Paulsen and Koorner (2014). At the sites in this analysis, treelines form in areas on a mountain slope where average growing season temperatures range between 5.5°C and 7.2°C, and the growing season length is approximately 150 days (defined by all days with an average temperature above 0.9°C). We also provide an estimate of the climate sensitivity of Great Basin treelines, and demonstrate that the treelines in this analysis are likely out of equilibrium with the current climate. Dendroclimatological studies would benefit from a comprehensive investigation of these findings. In conjunction with the recommendations of Salzer et al (2014), our results suggest that over time the spatial window of near-treeline temperature-sensitivity shifts on a centennial-millennial timescale. Treeline position modeling throughout the Holoscene would enable a more accurate isolation of a purely temperaturelimited signal. For any given year, chronologies reconstructing temperatures should only use trees that are within the temperature-sensitive window. As treeline position shifts up and down, individual trees Environ. Res. Lett. 12 (2017) 014008 should be added and removed from the reconstructive chronology depending on their location relative to the temperature-sensitive window. This would enable the cleanest temperature signal and potentially more accurate calibration for global climate models.