Searching for refuge: A framework for identifying site factors conferring resistance to climate‐driven vegetation change

Climate change is occurring at accelerated rates in high latitude regions such as Alaska, causing alterations in woody plant growth and associated ecosystem patterns and processes. Our aim is to assess the magnitude and speed that climate‐induced changes in woody plant distribution and volume may be reduced and/or slowed by relatively static landscape features like physical characteristics (e.g. depth to gravel, mineral cover percent and slope degree) and/or edaphic properties (e.g. soil organic matter, soil pH and site wetness rating) that resist climate‐vegetation responses.


| INTRODUC TI ON
Understanding woody plant response to climate change in Alaska is important, both globally and locally (Hinzman et al., 2005;Tape et al., 2016). From a global perspective, changes in the distribution of boreal forest have been shown to affect atmospheric processes, which, in turn, can cause amplification of global warming effects on woody plant distribution across these ecosystems (Helbig et al., 2016). Vegetation in Alaska is changing dramatically and at an accelerated pace resulting from complex climate-environment interactions (Beck et al., 2011;Brodie et al., 2019;Myers-Smith et al., 2020). Interior forests exhibit decreased vegetation caused by increased fire frequency and permafrost degradation caused by higher temperatures (Goetz et al., 2005;Lara et al., 2019;Pastick et al., 2019;Verbyla, 2008). These vegetation mortality trends may have many future consequences, including altering wildlife habitat, releasing increased amounts of carbon into the atmosphere and degrading socio-economic systems Lloyd & Bunn, 2007). By contrast, in alpine and Arctic ecosystems, there has been an increase in shrub distribution (Tape et al., 2006). Proposed consequences of increased vegetation in the Arctic include amplified warming (Betts & Ball, 1997), increased flammability (Frost & Epstein, 2014), changing snow deposition (Sturm et al., 2001) and reductions in species diversity in tundra plant communities (Roland et al., 2017). In addition to increased shrub volume, northward tree expansion has been widely observed over the 20th century (Harsch et al., 2009). These changes may initiate transitions in ecosystem trajectories, affecting biological and physical processes at multiple scales (Chapin et al., 2004;Hinzman et al., 2013;Pastick et al., 2019;Scheffer et al., 2012).
Understanding the compounding consequences of these interrelated components across a vast, heterogeneous spatial domain requires a detailed understanding of the ecological processes involved (Roland et al., 2013. Increased disturbance and disturbance severity has largely been driven by increases in fire in the boreal regions of Alaska (Barrett et al., 2010;Beck et al., 2011;Brown & Johnstone, 2012;Fang et al., 2015) and increases in permafrost degradation in Arctic regions of Alaska (Jorgenson et al., 2006). While it has been suggested that this increasing fire frequency may lead to a state shift or change in dominant vegetation type (Scheffer et al., 2012), it has also been found to be a variable spatial pattern that may be influenced by site-level species response that could lead to increased growth response from dominant species (Johnstone & Chapin, 2006). Increasing permafrost degradation may lead to both increased amount of suitable habitat for vegetation expansion (Roach et al., 2011(Roach et al., , 2013 and coastal erosion decreasing the amount of available habitat in some eco-regions of Alaska Jorgenson et al., 2006). Increased growing season length has led to increased drought stress amongst dominant spruce species (e.g. Picea glauca; Barber et al., 2000).
While much research in the past decade has focused on understanding the patterns and drivers of woody plant change, it is also important to assess which landscape features may inhibit vegetation expansion or offer plant refugia for particular plant species. Accelerating vegetation changes may be tempered by relatively 'static' landscape features in contrast to 'dynamic' climate Roland et al., 2019;Swanson, 2015).
Vegetation resistance to change may be facilitated by static terrain-mediated and ecosystem-protected processes, reflecting a gradient in relative degrees of 'inhibition' or 'refuge' (Stralberg et al., 2020). Static landscape features like site wetness in bogs and peatlands have been hypothesized to be ecosystem-protected features that allow for continued growth in the face of temperatureinduced drought or increased fire frequency . Ecosystem-protected landscape features are not usually included when forecasting woody plant migration at higher latitudes because plant volume response to these features has not been robustly determined across a variety of ecosystems. Alternatively, some static landscape features may inhibit vegetation expansion causing sites to be unsuitable for plant growth and provide terrainmediated resistance to climate changes. We seek to determine the extent that static variables may play a larger role in woody plant volume responses either inhibiting expansion or offering refugia in Alaska. It has been shown that, with increasing latitudes, the physical and edaphic environment becomes an increasingly dominant control on woody plant growth (Callaway et al., 2002;Hulshof et al., 2013;Klanderud et al., 2015;Pierce et al., 2017;Schemske & Mittelbach, 2017). Quantifying the effects of these static terrainmediated and ecosystem-protected landscape features on the potential rate of vegetation change in different ecological contexts in Alaska would greatly benefit management planning and reduce uncertainty in our understanding of plant response to dynamic climate .
We seek to improve our understanding of potential woody plant responses to climate change in Alaska by determining landscape features that may inhibit vegetation expansion or provide plant refugia.
In our analysis, we consider climate and climate-driven variables such as precipitation and burn status relatively 'dynamic' variables compared to the relatively 'static' physical landscape characteristics and edaphic conditions (hereafter we omit 'relatively'). In our definition, static landscape features are both terrain-mediated and ecosystemprotected processes that lead to a particular site being resistant to external forces. We refer to these sites as relatively 'robust' to climate changes or having higher 'robustness.' While our data set is comprehensive, we stress that robustness is not an absolute future state in our framework but a condition of higher resistance to future climate that can be estimated because of past independence from climate.
Due to its large extent and variety of ecosystems, Alaska encompasses substantial gradients that directly affect patterns in woody plant distribution and abundance, all of which need to be considered to effectively predict responses to climate change . We comprehensively assess woody plant volume response to dynamic climate and permafrost variables using field data collected across this wide spatial domain. In our study, we have three objectives: (1) to describe patterns in this data set with a uniquely large geographic extent and then leverage the data to develop spatially heterogeneous estimates of plant species robustness to climate (and climate-related disturbance) changes in woody plant volume across northern Alaska; (2) to learn where plant species may be either excluded or else harboured from climate changes, we used a large plant occurrence and volume data set and modelled species dependence on dynamic climate variables, assuming that sites where volume is accurately predicted without climate have higher 'robustness' to climate changes; and (3) to align our results with past work measuring temporal vegetation change to a lesser spatial extent with information from satellite imagery (Pastick et al., 2019) and oblique photograph pairs (Brodie et al., 2019) in northern Alaska.
We expected that the majority of sites would be classified as robust to climate changes because of the harsh-growing conditions and terrain-mediated, static barriers to vegetation expansion in Alaska (Roland et al., 2016Swanson, 2015). We also expected that landscape types with consistent access to surface-water such as peatlands or bogs provide static ecosystem-protected refugia to species that have the ability to occupy wet soils (Stralberg et al., 2020). Two specific species examples are Picea mariana (picmar), which are well known to frequently occupy wet, acidic soils (Mack et al., 2008;Viereck et al., 1983), and Salix pulchra (salpul), which have been found to respond differently to climate depending on soil wetness (Ackerman et al., 2017). We expected that robustness estimates for these species would depend largely on site static edaphic properties such as wetness because site wetness may allow these species to escape climate change pressures and static terrainmediated characteristics related to site wetness, like slope angle, may exclude them from expansion.

| Data
We analysed site-level woody plant volume data that were collected throughout a network of Alaskan National Parks and Preserves in- Bering Land Bridge National Preserve (BLNP) (Figure 1). Complete sampling design methodologies can be found in Roland et al. (2013) and Swanson (2015). These data consist of 2062 sites that are clustered at sample locations with approximately 20 sites per location.
Each site is 200 m 2 . These sites span wide latitudinal (62.11-68.48), longitudinal (−164.78 to −141.01) and elevational (2.26-1674 m) gradients. Along with the large spatial domain, these sites cover a broad climatic domain with total summer precipitation ranging from 11.44 cm to 77.66 cm and average July temperature ranging from 7.89 C to 16.45 C across sites ( Figure 2). These data were collected over 15 years (i.e., 2001-2015). The sites were typically measured once because of remote access barriers in the region. These data F I G U R E 1 Map of sample locations each containing approximately 20 sites and park boundaries in an Alaskan Network of National Parks represent our most current knowledge of plant volume and edaphic conditions for the sites, so we treat these data as a static snapshot (Roland et al., 2013).
The natural history of each species is distinct, which allows us to assess plant response variability throughout the woody plant community in terms of robustness to changes in climate on plant volume (Table 1).
Basal area (BA) for each tree species was derived using the diameters of all main stems at breast height (1.37 m, DBH, see Roland et al., 2013). For shrub species, canopy volume was calculated using vertically integrated point transects measured in the field and following the methods of Swanson (2015) where individual shrub height was measured, summed and then divided by total plot area (units = m 3 /m 2 ). We chose these measurements instead of using a common metric like biomass or percent cover for two reasons. First, biomass derivations for these shrub species are unreliable because allometric relationships are not well established and likely vary over this spatial domain. Second, while percent cover could be calculated for both trees and shrubs, there is less information in percent cover than either basal area or volume, which integrate the heights of the plants into the metric and are thus a better proxy for total volume than simple horizontal cover. Our basal area and volume data allow us to model the response of woody plant volume to changes in climate beyond what we could learn from occurrence-only or percent cover data sets. Throughout, we refer to both tree basal area and shrub volume as plant "volume." Comprehensive site-level covariate information, fully described in Roland et al. (2013), was also measured at each of the sites. We selected a subset of the covariate information that was found to most affected plant volume in previous studies using previous iterations of this data set (Brodie et al., 2019;Roland et al., 2013Roland et al., , 2019Swanson, 2015). Categorization of covariates into 'static' or 'dynamic' was based on expert understanding of temporal change of each covariate . We categorized covariates as 'static' if they were largely considered site ecosystem-type characterizations.
Some covariates categorized as static (i.e. flood disturbance, mineral cover percent, soil organic matter, soil pH, wetness rating, latitude and longitude) may change suddenly in rare circumstances such as a climate-related disturbance event, but typically do not change over long period (i.e. >100 years). We categorized covariates as 'dynamic' if they typically change and are recorded on an annual or semiannual basis. While there are many dynamic covariates to choose from, these variables were chosen based on previous research showing these variables most affected plant volume (Brodie et al., 2019;Roland et al., 2013Roland et al., , 2019. The static variables are wetness rating, depth to gravel (cm), equivalent latitude, flood disturbance (binary variable), mineral cover F I G U R E 2 Histograms of static and dynamic covariate values across 2062 vegetation sites in an Alaskan National Parks and Preserve Network. The top two rows are static variables while the bottom two rows are the dynamic variables (orange boxes). The time windows written in text in the upper right hand corners of each plot represent average temporal change based on ecological knowledge and previous work by Roland et al. (2013Roland et al. ( , 2019, Swanson (2015), and Brodie et al. (2019). The sites with an asterisk (*) indicate covariates that typically change slowly but may change suddenly in certain circumstances Fast growing.

PICMAR Black spruce
Coniferous tree Common in mesic sites.

Tolerant
High mortality. Limited postfire recruitment. Slower growing but tolerant of cold, wet soils and has shallow roots.

Tamarak
Coniferous tree Limited to Denali National Park.
Typically found in bogs.

POPBAL Balsam Poplar
Deciduous tree River valleys and floodplains.

Intolerant
Moderate mortality. Rapidly recruits or regenerates postfire by seed and resprouting.
Floodplain dominant, most northerly distributed tree species in Alaska.

Aspen
Deciduous tree Common on south facing slopes.
Occupies warm, dry end of forest gradient in Alaska.

Alaksa Birch
Deciduous tree Common amongst spruce species.
Somewhat tolerant High mortality. Rapidly recruits postfire by seed and resprouting. Occurs with both white and black spruce in a relatively wide variety of forest types.

Diamond leaf willow
Willow shrub Dominant tundra plant commonly grows with sedges and near water.

Intolerant
Often sprouts after wildfire.
Most widely distributed willow species, relatively tolerant to cold, organic soil conditions.

SALGLA Gray willow
Willow shrub Often found in riparian areas but can be found in forests and woodlands.
Important food for ungulates in winter.

Alaska willow
Willow shrub Ravine species.

Intolerant
Often resilient to fires.
Floodplain specialist, preferred moose browse.

Wooly willow
Willow shrub Forms closed thickets in open tundra often associated with Alnus viridis.
Low tolerance Sprouts rapidly from roots postfire.
Considered an indicator of better growing conditions at high elevation.

ALNVIR Green alder
Nonwillow shrub Requires moist soil. Commonly found amongst willows.

Tolerant
Fast growing colonist after disturbance.
Nitrogen fixing. Shallow root system.

Dwarf birch
Nonwillow shrub Well drained, nutrient poor, acidic sites. Often found at higher elevations.
Dominant species over large areas of shrub tundra.
Note: Divisions represent taxonomic groups conifer, broadleaf, willow and nonwillow. Notable characteristics are based on expert knowledge of the system. percent, soil organic matter depth (cm), slope degree and soil pH. We also used squared variables for pH, organic matter depth, percent mineral cover and depth to gravel because these variables are known to saturate in this study area (Roland et al., 2013). The dynamic variables are burn status (binary variables for old burn and recent burn derived from year of burn data), July mean temperature (C), summer total precipitation (cm), thaw depth (cm) and snow-free Julian date.
Burn status was recorded as true or false for both old burns (20-80 years) and recent burns (less than 20 years). July mean temperature, total summer precipitation and snow-free Julian date were taken from gridded data for Alaska by the PRISM Climate Group, which represent rolling 30-year climate normals of the most recent 30 years.
Thaw depth, depth to gravel and soil organic matter depth were collected at each site from a small (30-40 cm) soil pit. Wetness rating is based on the weighted averaging method from Federal Interagency Committee for Wetland Delineation (1989) (see Atkinson et al., 1993) where each species in a site are assigned a wetness rating. Ratings for plant species were developed by an interagency expert process (Lichvar et al., 2016) and downloaded from the USDA Plants database (USDA, NRCS, 2020), but a small number of missing taxa were assigned a wetness rating by the authors. The wetness rating of a site was the average of the constituent species' wetness ratings, weighted by their canopy cover. A higher rating indicates drier site conditions.
Throughout this manuscript, we refer to sites with higher wetness ratings as being "drier" or having more "dryness" to avoid confusion.
See Roland et al. (2013) for covariate definition information.

| Statistical analysis
We use a hierarchical Bayesian model to estimate where dynamic covariate information improves predictions of plant volume for a variety of species (adapted from Scharf et al., 2021). We consider sites 'robust' if plant volumes are largely predictable by site-level static variables, for instance, percent mineral coverage or soil pH ( Figure 3 right-hand side, p < .33 described below). In contrast, sites are clustered into the 'nonrobust' category if they depend on information from dynamic variables such as July mean temperature, total summer precipitation, depth to thaw or recent burn status after accounting for static sitelevel variables (Figure 3 left-hand side, p > .67 described below). In particular, our modelling framework identifies site characteristics that could facilitate climate refugia for particular species by determining sites that are currently predictably occupied without including dynamic climate or climate-related disturbance covariates. Similarly, sites that are predictably unoccupied by a particular species without information from dynamic covariates are considered 'excluding' landscapes for that species. Figure 4 provides a conceptual picture to the model specification described below. Furthermore, while this is a study using spatial information only, we also make brief comparison of our results in the discussion to photograph pair images in the same study area (Swanson, 2015) to validate our results.

| Model specification
Our model consists of three components in a Bayesian hierarchical framework: data, process and parameter models (Figure 4; Berliner, 1996). We utilized plant log-volume data y i as response variables collected across sites i = 1, …, n (Equation 1). Observations of logvolume are assumed to be unbiased and centred at the true logvolume ( i ). We modelled the detection process explicitly using a censored data model where observed log-volume is greater than the lower detection limit y * i (i.e. the minimum volume observed for each species). Our specification results in a mixture Tobit left-censored data model for observed log-volume y i that we express as where z i is a latent cluster indicator that equals 0 when a site is robust and 1 when nonrobust. In our case, the left-censored value is the minimum species volume detected in the data set. Using the Tobit model allows us to estimate unobserved plant volumes below the lower Our definition of robustness along climate (or dynamic) variable independence (horizontal axis) and plant volume response (vertical axis) gradients. The four categories within our diagram describe useful ecological concepts that are identifiable within our framework. The variable p represents the probability of a particular site being categorized as robust or nonrobust detection limit and also allows us to separate the data level uncertainty from the clustering process described below (i.e. Figure 4 middle from bottom level). We defined two true log-volumes ind i and dep i as functions of observed site-level covariates such that The covariates comprise static x (s) i and dynamic x (d) i environmental variables described above, and (s) and (d) represent their associated coefficients. Each also includes a coefficient for the intercept through the addition of a column of "1" in the covariate matrix. Our model coincides with the understanding that the dynamic variables should not affect the latent log-volume process for the sites considered robust. This model specification falls into the broader class of mixture models and allows us to discern sites that are robust (z i = 0 ) from nonrobust (z i = 1) based on whether we require information about dynamic variables (e.g. climate and climate-related disturbance) to predict log-volume. We show in Appendix S1 that the model formulation in Equation (1) represents a type of regularization where the dynamic variable effects shrink to zero for robust sites (Hooten & Hobbs, 2015).
To account for potential confounding amongst covariates, we restricted the dynamic covariates (x (d) i ), such that they are orthogonal to the static covariates (x (s) i ) for each species. This procedure is commonly used in spatial statistics to alleviate confounding (e.g. Hanks et al., 2015). Specifically, using the full static design matrix X (s) , we computed the restricted dynamic design matrix as where X (u) is the unrestricted dynamic design matrix containing climate-related variables. This transformation allows the component of the model based on dynamic variables to involve only climate-related information that is not already accounted for in the static variables (and reduces overall multicollinearity). We use ⊥ to denote that the dynamic variables are orthogonal to the static variables. For description of restricted regression and use of orthogonal variables in ecology, see Hanks et al., 2015. We specified the latent cluster membership model as i (p) characterizes the heterogeneity in robustness across the landscape. We define robust sites as sites with posterior mean cluster probabilities (p i ) less than 0.33 while nonrobust sites are sites with posterior mean cluster probability (p i ) greater than 0.67. We consider sites with posterior mean values between 0.33 and 0.67 as sites that cannot be categorized with confidence (i.e. uncategorized). We selected these thresholds because they indicate a clear majority. A variety of choices can be used to define x (p) i . For example, we could use the set of static covariates x (p) = x (s) . In this setting, our model allows us to understand which slowly changing landscape features may limit vegetation encroachment given dynamic changes in the environment. We can account for additional flexibility in the relationship between the static variables and robustness using a semiparametric relationship in which case x (p) is defined in terms of basis functions that represent the static covariate space (Hefley et al., 2017). These basis functions could be splines that span the static covariate space or, as we use in (4) X (d) = (I − X (s) (X �(s) X (s) ) −1 X �(s) )X (u) , F I G U R E 4 Directed acyclic graph (DAG) of our model showing the model levels: data model, process model, and model parameters and covariates. Circles indicate stochastic nodes estimated by the model while squares are deterministic nodes. Diamonds indicate model inputs. y i is the plant volume data by site i. x ′(s) i (i.e. wetness rating, depth to gravel (cm), equivalent latitude, flood disturbance (binary variable), mineral cover percent, soil organic matter depth (cm), slope degree and soil pH as well as squared depth to gravel, mineral cover percent, soil organic matter and soil pH), x ′(d) i (i.e., burn status, July mean temperature (C), summer total precipitation (cm), thaw depth (cm) and snowfree Julian date), and x ′ (p) i (i.e. same as (s)) are the covariate sets used along with coefficients (s) , (d) and (p) to obtain posterior estimates of cluster membership (z i ) and cluster membership probabilities (p i ) what follows and describe in more detail in Appendix S1, the basis functions could be derived from a first-order spectral representation of a Gaussian process model on logit(p i ). This semiparameteric representation allows us to account for highly nonlinear relationships between robustness and the complete set of site-level static variables. To complete the hierarchical model, we specified priors for the coefficients as multivariate normal such that, (s) ∼ Normal(0, (s) ), (d) ∼ Normal(0, (d) ) and (p) ∼ Normal(0, (p) ) where the coefficients are assumed independent here but are expressed generally for future implementations (see S2). Similarly, we expressed the prior for the variance generally for future implementation where 2 ∼ IG (q, r). We fit the model described above for each species individually.
Analysis on covariate predictive ability and model parsimony can be found in Appendix S1. More information about model fitting and model fit evaluation can be found in Appendix S2. We briefly evaluate our model using photograph pair images collected by Swanson (2013)  Information with specific code descriptions in the Data S1.

| Volume and covariate data
The most abundant tree species across our study domain was P.
glauca, which occurred in 36% of sites (basal area = 7.41 m 2 /ha on average across parks). The least abundant tree species was L. laricina, which was only present in DNPP (occupying 4% of all sites, basal area = 0.23 m 2 /ha on average in DNPP). Tree species occurred much more frequently in southeastern parks than northern parks (Appendix S1: Table S1). There were more shrubs than trees in Arctic  Table 2 and Figure S6.  Note: P represents "present-only" sites. O represents "overall" sites. Table describing summary statistics for volume measurements across species where basal area (m 2 /ha) was measured for trees and volume (m 3 /m 2 ) for shrubs relatively high elevation sites. Lastly, years since burned was skewed toward mostly unburned sites (>100 years since burned Figure 5). were not present at recently burned sites (<60 years since burned).

| Coefficient estimates
The relationship between wetness rating and plant volume was estimated with high precision and nonoverlapping zero for the majority of species where it was an included covariate. Wetness rating was estimated to have a slightly positive effect on volume across F I G U R E 5 Violin plots showing covariate distributions across all sites and across sites where species are present. White dots represent the median of each distribution, indicating the environment type where a particular species is most likely to be present in our data set. We include elevation to demonstrate the differing niche space of the species, but elevation was not included in the final model because it covaries very closely with July mean temperature (see Roland et al., 2019). Similarly, 'year since burn' was used to derived binary variables 'old burn' (20-80 years) and 'recent burn' (less than 20 years) that were used in the final model fitting. Years since burned > species (mean = 0.48, SD = 1.14), meaning drier than average sites had higher woody plant volumes. Also, across species, organic matter depth and depth to gravel were most positively associated with plant volume (organic matter depth mean = 1.34, SD = 2.81; depth to gravel mean = 1.11, SD = 2.97). Thaw depth ⊥ coefficients were estimated with the most precision across species followed by recent burn status ⊥ . Overall species thaw depth ⊥ was estimated to be positively related to volume (mean = 0.60, SD = 0.82) while recent burn status ⊥ was estimated to be negatively related to volume (mean = −0.26, SD = 1.00) with the exception of P. tremuloides for which recent burn status ⊥ was positively related to volume (mean = 2.27, SD = 0.37).

| Robustness categorization
The Lastly, we found that our robustness categorizations for P. mariana and P. glauca closely matched observations from photograph pair change detection.
We found that wetness rating and organic matter depth (SOM) were often indicative of robustness ( were more robust (mostly excluded) in sites with higher organic matter depths. Across species, the distribution of static covariates was not significantly different in most cases for robust sites relative to all sites (75% of cases had more than 0.8 overlap in distribution, i.e. semitransparent in Appendix S1: Figure S7). This was because most of the site-species combinations were estimated to be robust (mean = 72.2%) with B. nana having the most nonrobust sites (51.5%).

| DISCUSS ION
We accomplished three objectives with this work. First, we described the patterns of a novel data set from a network of Alaskan National Parks and used that data set to develop a model for estimating plant robustness to climate (and climate-driven) changes.
Second, we determined which landscape types may either inhibit vegetation expansion or offer climate refugia to specific plant species. In this discussion, we elaborate on these findings and elaborate on our third objective of aligning our results with past work focussing on temporal change.
The majority of site-species combinations were categorized as excluding landscape types (Figure 3 lower right; Figure 7). This trend correlated with species occurrence patterns and terrain-mediated patterns suggesting certain landscape types that may prevent woody plant encroachment. Soil and landscape characteristics in Alaskan ecosystems exert important influences on species distribution and volume patterns (Brodie et al., 2019;Callaway et al., 2002;Hulshof et al., 2013;Klanderud et al., 2015;Pierce et al., 2017;Roland et al., 2013Roland et al., , 2019Schemske & Mittelbach, 2017;Van Cleve et al., 1986;Viereck et al., 1983). Soil characteristics, such as soil pH or soil organic matter, can prevent woody plants from colonizing new areas because certain soil conditions inhibit establishment and/or growth of particular species (e.g. S. pulchra, Ackerman et al., 2017). For example, the nonconifer tree species in Alaska rarely, if ever, occur in high volumes on soils with a deep organic horizon (Alexander et al., 2012, Figure 5). Specifically, we found that site dryness, slope angle and organic matter depth were important indicators of robustness to dynamic variables across species (Table 3). Species P. mariana, S. pulchra and B. nana were excluded from growing in sites that were drier with steeper slopes and lower organic matter depths. Conversely, the tree species P. glauca, P. tremuloides and B. neoalaskana rarely occur in wetter sites, with shallower slopes and higher organic matter depths. Our results suggest that terrain-mediated characteristics, such as slope degree or mineral cover percent, in boreal and Arctic ecosystems may inhibit the expansion of certain species.
While exclusionary landscape types were dominant in our analysis, we also improved our understanding of conditions where species have the best chance of finding climate refuge (Figure 3 upper right).
Soil wetness rating was the dominant driving factor for determining species robustness ( Figure 6) and may offer spruce species, in particular, ecosystem-protected refuge (Stralberg et al., 2020) and attenuate ecosystem transition to deciduous forest (Johnstone & Chapin, 2006;Scheffer et al., 2012). Boreal peatlands are a notable example of an ecosystem-protected refuge where, in the face of environmental stressors such as severe drought, high water table depths can continue to promote plant growth for specific boreal species (Shur & Jorgenson, 2007;Stralberg et al., 2020). For instance, recent work illustrates that peatlands provide refuge at the southern range margins of P. mariana distribution (Langdon et al., 2020;Spei & Kashian, 2018). In our data set, boreal peatlands are represented via the 'wetness rating' and 'organic thickness' static covariates.
Over 50% of sites across the Alaskan network of National Parks can be classified as wetter (with a rating less than 3). We found that dryness was strongly related to decreased robustness across tree species (Table 3), indicating that peatland in this region may provide refugia for bog-dwelling boreal species as temperatures increase or fires remove top organic layers . Furthermore, peatlands currently make up 19% of Far North ecosystems (Tarnocai et al., 2009) and may provide substantial refugia from environmental stressors not yet accounted for forecasts of carbon or habitat dynamics. While past studies have hypothesized that wet ecosystems offer the potential to alleviate soil moisture growth limitations and therefore offer a potential long-term refugia for boreal species, our results are the first to confirm this finding using an expansive field data set collected across broad environmental gradients. These results bolster existing support for targeted protection of Far North wetlands (Kåresdotter et al., 2021).
The strength of our modeling framework is to highlight areas that are robust to climate changes. However, we also have the ability to identify which landscape types are particularly sensitive (i.e. nonrobust). Arctic and boreal vegetation may be able to expand into previously unoccupied sites as climate becomes more conducive for growth, unlocking suitable habitat (Euskirchen et al., 2009, Sites that were occupied and nonrobust ( Figure 3,  Note: Visualization of distributions and further description can be found in Appendix S1: Figure S7. differences in stand ages and therefore structural variability of spruce species amongst sites (Lloyd & Bunn, 2007;Nicklen et al., 2019). Variation in site moisture drove differences in relative sensitivity between the two spruce species. We found that P. mariana was present and robust ( Figure 3, upper right, refugia) because of the site dryness differences, while P. glauca was present and nonrobust ( Figure 3, upper left, structural variability) across a wider range of dryness (Table 3). Furthermore, we found that P. glauca volume responded more negatively to recent burn while P. mariana volume responded more positively to mean summer temperature, possibly initiating variability in growth response across our wide gradient of sites. Willow volume has been shown to respond more uniformly to increasing temperatures than tree volume (Ackerman et al., 2018).
Our results agree with this assessment where willow species had a more even response to dynamic variables than other species groups ( Figure 6, bottom).
To illustrate our results for particular species, we highlight P.
mariana and S. pulchra. Increased fire frequency and severity poses a threat to seedling establishment of P. mariana in Alaska. Previous studies of P. mariana have found that effective establishment of seedlings can occur on mineral soil, thin organic soil and sphagnum because they can be continually moist (Mack et al., 2008;Viereck et al., 1983), but if high severity fire reduces organic matter depth, P. mariana may be out-competed by faster growing P. glauca F I G U R E 6 Tile plots of posterior median of the coefficient for the static covariates ( (s) , top) and the dynamic covariates ( (d) , bottom). These signs can be generally interpreted as the marginal effects on plant volume. In the static coefficients, we interpret the squared term as concave up (+) or down (−) and the linear term as the inflection point in the effect of response to plant species volume. The dynamic covariates were used as part of a restricted regression (denoted with ⊥) therefore are not directly interpretable but highly correlated with named dynamic variables. Points on a square to indicate coefficients that were estimated with high probability for a nonzero effect, meaning the 95% credible intervals did not overlap zero  (Roland et al., 2013;Wirth et al., 2008). Two trajectories are possible for P. mariana to escape these pressures: Seeking refuge in bogs and peatlands that are less affected by fires and less occupied by competitors and expanding north into potential habitat . Our results show that both of these scenarios are possible but that more sites provided refuge (5.43%) than potential habitat (3.09%) where soil conditions such as site wetness, soil organic matter depth and depth to gravel were the best indicators of robustness. Furthermore, our results indicate that P. mariana volume may respond positively to increasing summer temperatures ( Figure 6) at sites that are suitable for P. mariana occupancy (i.e. potential habitat in lowland areas to the north). Shrub species are expanding in both distribution and volume because of increased growing season length in some areas of the Alaskan Arctic (Tape et al., 2006). S. pulchra has been shown to have higher growth rates in sites with wetter conditions (Ackerman et al., 2017 1%). We used a logistic regression to quantify the relationship between our robustness score (p) and the real change observed in the photograph pairs. We found that, for both spruce species, the posterior mean of robustness was significantly related to actual change in the photograph pairs in the Arctic National Parks from Swanson (2013). Specifically, in the more southern photograph pairs, there was evidence of increased tree growth. The photograph pair evidence matches our estimates of robustness where less robustness to climate changes at lower latitudes in the Arctic National Park occurs because site conditions are more conducive to vegetation growth. In Arctic Alaska, shrubs have expanded more than trees (Tape et al., 2006). Specifically, A. viridis has exhibited the most change across the photograph sample interval in Swanson (2013).
Our results identified increased thaw depth as predictive of increasing A. viridis volume, which may be a causal factor for the observed changes caused by amplifying temperatures that the Far North has experienced over the past 30 years Jorgenson et al., 2006; Figure 6).
Our model facilitated inference about landscape types in Alaska that may temper vegetation change caused by changing dynamic variables like climate or disturbance regimes. However, three limitations suggest important future directions for our approach are: (1) Our data do not contain information about temporal change at specific locations; thus, we compare volume across sites that are similar landscape types. As described above, we can compare our estimates of robustness to actual observed change between photograph pairs, but future work should seek to combine snapshot field data sets with time-series data sets. Our field data set is derived from a long-term monitoring program that, with repeated sampling, will allow future researchers to assess our projections directly over time.
(2) While the variation of the predictive variables in our data set is extensive, the climate variability in Alaska is changing beyond the historic normals (Hinzman et al., 2005), and therefore, our data set may not completely represent future climate regimes.
(3) We segregated static from dynamic covariates based on the relative temporal scales of change. However, some of the static covariates (i.e. soil organic matter) in our model may change more rapidly under climate change scenarios (Euskirchen et al., 2009) based on warming experiments (Natali et al., 2011) (Figure 2 top with *). Future work should consider the possibility that static covariates may become more dynamic in the future. Similarly, we have also included some static covariates that may be confounded with dynamic covariates (e.g. latitude and longitude; see Figure S1).
To confront this issue, we assessed the correlation between static and dynamic covariates ( Figure S1) and determined that the static covariates were not overly confounded with the dynamic covariates. However, in future research where static and dynamic covariates may be more confounded, researchers may be more explicit including spatial random effects, removing latitude and longitude from the restricted regression or utilizing the principle components instead of restricted regression. We specified our model generally to allow for future model implementations flexibility to adjust this model specification if needed.
Including woody plant volume data in our analyses allowed us to make inference beyond simple species distribution and pinpoint mechanisms that may exclude species volume changes affected by rapid warming. Our model can also be extended to accommodate temporal data, as they accumulate in the future. As more plant volume and microclimate data become available from remote sensing operations (Zellweger et al., 2019), it will be advisable to incorporate those data to improve our understanding of plant resistance to environmental change. Furthermore, understanding species differential responses to climate change will help with predictive mechanistic modeling. High latitudes are facing accelerated effects of global climate change. To predict and possibly prevent unwanted ecosystem change, it is critical to understand which factors may be accelerating or alleviating woody plant response to climate change.

ACK N OWLED G EM ENTS
This research was funded by the National Park Service and NSF DEB 1927177. We thank numerous technicians and pilots for assistance in the field. We thank A. Southwould, E. Debevec and D. Wilder for programming and database assistance. Any use of trade, firm or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13492. This submission uses novel code, which is provided in the Supporting Information and described in Data S1.