Pathways of tundra encroachment by trees and tall shrubs in the western Brooks Range of Alaska

1 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Eric Post Editor-in-Chief: Miguel Araújo Accepted 7 January 2020 43: 1–10, 2020 doi: 10.1111/ecog.05015 43 1–10


Introduction
The Earth's surface warmed by an average of 0.85°C between 1880 and 2012 (IPCC 2014). Spatial variation in the warming trend is considerable, with high latitudes and high elevations warming at above average rates -phenomena known as Arctic Amplification (Serreze and Barry 2011) and Elevational Dependent Warming (Pepin et al. 2015) respectively. Arctic Amplification is largely attributable to reduced albedo associated with dramatic declines in Arctic sea ice extent (Screen and Simmonds 2010) and high latitude spring snow cover (Matsumura et al. 2014) and has led to rates of Arctic warming that are more than twice the global mean. Amplified warming has had profound effects on the species composition and functioning of Arctic ecosystems, with further dramatic changes anticipated for the future (Post et al. 2009). Likewise, Elevational Dependent Warming has substantial implications for glaciers, phenology of alpine ecosystems and hydrology (Pepin et al. 2015). Taken together, the two phenomena point to the potential for very rapid warming in arctic mountain ranges. Unfortunately, the paucity of temperature records from these regions leaves that hypothesis untested (Wang et al. 2016).
A shift in dominant vegetation type is one of the more dramatic changes that has occurred in many areas of the Arctic. Numerous studies have examined the pattern and process of deciduous shrub encroachment into arctic tundra (Sturm et al. 2001, Tape et al. 2006, Myers-Smith et al. 2011, Tape et al. 2012, Lantz et al. 2013). Meanwhile, relatively few studies have examined Arctic treeline advance (Lescop-Sinclair and Payette 1995, Suarez et al. 1999, Lloyd et al. 2003, Harsch et al. 2009, Lantz et al. 2019, despite the comparative ease of differentiating trees from tundra in historical imagery. Changes in tree abundance in the Arctic have the potential to dramatically alter surface energy exchange (Beringer et al. 2005), carbon (C) cycling (LaFleur et al. 2001), wildlife habitat (Fullman et al. 2017) and the availability of subsistence resources to local communities (Kofinas et al. 2010). Improved understanding of the pattern, process and magnitude of increasing tree abundance in the Arctic is therefore an important goal.
Temperature is widely thought to control global treeline position (Körner 1998, Körner 2012, with the transition from forest to tundra generally occurring where growing season mean air and soil temperatures are between 6 and 7°C (Körner and Paulsen 2004). The global correlation between treeline elevation and temperature implies that climate warming should lead to treeline advance. However, numerous studies have pointed out that moisture limitation (Moyes et al. 2015), nutrient limitation (Sullivan et al. 2015, Ellison et al. 2019, geomorphology (Macias-Fauria and Johnson 2013), competition with tall shrubs (Dial et al. 2016) and herbivory (Speed et al. 2010) may be more important constraints than temperature in some regions. A global synthesis revealed that only 50% of the treelines examined show evidence of recent advance (Harsch et al. 2009), lending support to the notion that treeline positions are controlled by more than temperature alone.
In this study, we aimed to improve understanding of the pattern, process and magnitude of changes in tall shrub and tree abundance in the well-studied Agashashok River watershed in the western Brooks Range of Alaska. To this end, we took advantage of orthorectified aerial and satellite images from 1952, 1979 and 2015. Our specific objectives were as follows: 1) identify the primary pathways of vegetation change between 1952 and 1979 and between 1979 and 2015, including transitional and apparent 'climax' states. 2) Quantify elevational shifts of vegetation classes over time and the key physiographic controls on tree abundance in 2015 (e.g. slope, aspect, elevation). 3) Examine the correspondence or lack thereof between the current (2015) treeline elevation and the 6 and 7°C growing season mean air temperature isotherm of Körner and Paulsen (2004).

Study area and climate data
We focused our study on a 5457 ha mountainous area roughly centered on the main stem of the Agashahok River ( Fig. 1), approximately 30 km upstream from its confluence with the Noatak River in the Baird Mountains of the western Brooks Range (67.5°, −162.2°). The area is accessible by bush plane and has been the subject of considerable biogeochemical, tree physiological and ecosystem ecology research. This area is unique in the sense that aerial imagery is available from the 1950s and climate data are available from multiple elevations within a small area. For this study, we took advantage of longterm air temperature data (1950-present) from the Kotzebue airport (66.89°, −162.61°, 3 m a.s.l., ~65 km south of the study area). We compared the Kotzebue air temperature record with shorter-term data from a station that we maintain in the floodplain of the study area (95 m a.s.l., 2006-present, Sullivan et al. 2015) and a ridge-top station that is maintained by the United States National Park Service (NPS, 405 m a.s.l., 2012-present). The Kotzebue air temperature data were obtained from the National Center for Environmental Information at the National Oceanic and Atmospheric Administration (<www. ncdc.noaa.gov/>), while the data from the NPS Asik station were obtained from the Western Regional Climate Center at the Desert Research Institute (<https://wrcc.dri.edu/>).

Remote sensing data
We used a time series of remote sensing imagery for three dates: an aerial black-and-white image acquired on 1952/07/20 (1:40 000), an aerial color image from 1979/07/14 (1:66 000) and a pan-sharpened Worldview-2 (WV-2) image from 2015/06/21 (2 m pixels). Spatial resolutions of the imagery were 1-2 m for aerial and 0.5 m for WV-2 image. The 1952 and 1979 aerial images were downloaded from the United States Geological Survey (USGS) EarthExplorer webpage (<https://earthexplorer.usgs.gov>), while the 2015 WV-2 image was acquired from the Polar Geospatial Data Center (PGC) at the Univ. of Minnesota (<www.pgc.umn. edu>). During August of 2017, we measured more than 100 ground control points on exposed rocks that were visible in all three images using a single-frequency GPS receiver with an overall accuracy up to 30 cm (Trimble, model GeoXH 2008; TerraSync software). All remote sensing images were orthorectified in ERDAS Imagine 2017 (Hexagon Geospatial, Madison, AL) using the Arctic digital elevation model (DEM, 5 m resolution, also acquired from the PGC). Treelines were manually digitized, such that the lines connect trees that are no more than 15 m apart. Digitization of the treeline was performed for illustrative purposes only and was not used for quantitative analyses. Red stars show the locations of the two weather stations (at 95 and 405 m a.s.l.) used to calculate lapse rates in recent years. The inset table shows the effect of elevation on air temperature when averaged across June, July and August and for July alone during the period of overlap between the two weather stations (2013)(2014)(2015)(2016)(2017). A small void in the Arctic DEM necessitated elimination of a small area in the northwest-central portion of the study area from our analyses.
For orthorectification of the 2015 WV-2 image, we used a rational polynomial function model with polynomial coefficients provided by Digital Globe and five control points measured in the field. For the aerial images of 1952 and 1979 we used a frame camera model, with camera parameters provided by the USGS and 42-50 control points measured in the field or from the reference WV-2 image. The root mean square error (RMSE) of orthorectification was 3.23, 2.01, 0.04 m for the 1952, 1979 and 2015 images, respectively.

Classification of random points
To examine changes in vegetation over time, we overlaid 3053 random points on all three orthorectified images. Each point was classified visually in each image as unvegetated, alpine tundra, arctic tundra, tall shrubs (> 1.5 m tall), trees with tall shrubs, open woodland, forest, floodplain tundra and alluvium based upon the cover type within a 3 m radius of the point. The unvegetated category was primarily unvegetated alpine terrain and rock outcrops at lower elevations and was distinguished from unvegetated terrain within the floodplain, which was classified as alluvium. The distinction between arctic and alpine tundra was made based upon vegetation type, taking advantage of our observations over 13 yr of fieldwork in the area. Drier heath vegetation on exposed ridges and on steeper slopes at higher elevation was classified as alpine tundra, while moist and wet tundra dominated by sedges and dwarf deciduous shrubs were classified as arctic tundra. Points classified as arctic tundra were generally lower in elevation and on shallower slopes than those classified as alpine tundra. If there was at least one tall shrub within the 3 m radius, the point was classified as tall shrubs. If there was only one tree within the 3 m radius and surrounding trees were consistently separated by > 3 m, the point was classified as open woodland. If the 3 m radius contained both a tree and a tall shrub, the point was classified as trees with shrubs. Finally, if there were multiple trees within the 3 m radius and surrounding trees were generally spaced at < 3 m, the point was classified as forest.
Relative changes in the abundance of vegetation classes were calculated as the difference between the number of points in the more recent observation (e.g. n 2015 ) and the number of points the previous observation (e.g. n 1952 ), divided by the number of points in the previous observation (e.g. relative change = (n 2015 − n 1952 )/n 1952 ). Areal changes in each vegetation class were estimated by assuming that our 3053 random points were representative of the 5457 ha study area. The effect of slope on area was addressed by increasing the area of each 3 m radius point as the secant (1/cos) of slope (Rinas et al. 2017). The study area was similarly increased to 5540 ha as the secant of the mean slope (10°).

Modeling elevational shifts and physiographic determinants of tree abundance
We sought to describe the spatial distribution of vegetation types using topographic variables derived from the DEM, such as elevation, slope angle, aspect and others, like 'northness' (cosine of aspect). For this purpose, we eliminated points classified as alluvium from the dataset and aggregated classes that included trees to create three vegetation classes: tree, arctic tundra and alpine tundra. Unvegetated areas at high elevation were lumped with alpine tundra and tall shrubs were split between alpine and arctic tundra, depending upon the adjacent tundra type. We attempted to model tall shrubs as a fourth vegetation class, but found that doing so decreased model performance, likely because tall shrubs represent a small proportion of the landscape. We then used multinomial modeling with the nnet package (Venables and Ripley 2002) in R 3.4.4 (R Core Team) to determine the probability of each sample point being one of the three vegetation types as a function of the topographic variables in 1952, 1979, following Rinas et al. (2017. Multinomial modeling with a GLM-style approach was preferred over machine learning (e.g. random forest), because the structure of algebraic models can be informative. For instance, models including quadratic (i.e. squared) terms indicate that vegetation classes may not respond monotonically to topographic variables. The most ecologically reasonable and most parsimonious multinomial model was identified with the aid of AIC. The best multinomial model was cross validated by fitting to a random subset (50% of the dataset) and then testing model predictions against the remainder of the dataset. Confusion matrices were produced using the caret package in R (Kuhn 2019) to assess model performance.
The best multinomial model was fit separately to the simplified point classifications (tree, arctic tundra, alpine tundra) for 1952, 1979 and 2015. While holding other variables in the model at their mean values, the elevation coefficients in each model were used to examine changes over time in the elevations at which there was a 0.5 probability of a point being classified as 'tree', alpine tundra or arctic tundra. In the case of quadratic terms, the elevation of the peak probability was calculated as in Dial et al. (2016).

Air temperature trends
The Kotzebue air temperature record shows statistically significant warming trends in June, July and September, with limited differences in trends between the earlier and later periods (Fig. 2). Air temperature outside of the growing season has also increased strongly in Kotzebue, such that mean annual air temperature increased by 0.4°C/decade between 1950 and 2018 (Supplementary material Appendix 1 Fig. A1). Between 2006 and 2018, air temperature at our weather station (95 m a.s.l.) in the Agashashok waterhshed was generally well correlated with Kotzebue air temperature, with a few exceptions (Fig. 2, Supplementary material Appendix 1 Fig. A2). The most notable exception is that June air temperature in the Agashashok watershed was generally warmer and not well-correlated with the Kotzebue data, likely reflecting the effects of sea ice and associated fog on June air temperature in Kotzebue.

Changes in vegetation classes of random points
During the earlier time period , 16% of the 3053 random points changed vegetation class (Supplementary  material Appendix 1 Table A1). During the later time period  18% of the points changed vegetation type (Supplementary material Appendix 1 Table A2). Over the full 63-yr period 32% of the random points changed vegetation classification at least once (Supplementary material Appendix 1 Table A3) and 3% of points changed vegetation classification twice. In general, the pathways and rates of vegetation change were similar between the earlier and the later time periods (Fig. 3, 4). Unvegetated areas, which generally occurred at the highest elevations, declined by 18% over the 63-yr period and were most likely to transition to alpine tundra. Despite the encroachment of alpine tundra into unvegetated areas, there was a 15% net loss of alpine tundra over the 63-yr period, owing to initial transitions to both tall shrubs and open woodland. Arctic tundra declined by 31% over the full period of the study, giving way initially to tall shrubs and open woodland, with the former being the more common transition. Tall shrubs increased strongly over the full study period (+86%) and more than one-third of the tall shrub communities present in 1952 gained sparse trees by 2015. The trees with tall  (−2%). Meanwhile, the closed canopy forest class increased strongly (+84%) and showed the greatest change in areal extent of all vegetation classes, with a gain of 73 ha over the 63-yr period. Notably, there was very little evidence of transitions from vegetation classes containing tall shrubs (tall shrubs and trees with shrubs) to forest. There was also little evidence of transitions between the trees with shrubs and forest classes, suggesting that they may be alternative stable states in the absence of disturbance. Alluvium was considered separately from the unvegetated class. While there was considerable bi-directional flux between alluvium and floodplain tundra, there was a net decline in alluvium of 16% over the 63-yr period.

Elevational shifts of vegetation classes
The elevation at which there was a 0.5 probability of a point being classified as 'tree' increased 18 m over 27 yr between 1952 and 1979 and an additional 20 m over 36 yr between 1979 and 2015 (Fig. 5). Meanwhile, the elevations associated with a 0.5 probability of alpine tundra also shifted up, though by lesser amounts (1952-1979: +6 m, 1979-2015: +9 m), in part because trees invaded both arctic and alpine tundra. The elevation associated with the maximum probability of encountering arctic tundra also increased by a few meters in each of the time periods. The small elevational shift of arctic tundra reflects its occurrence in areas with lower slope. On an areal basis, more arctic tundra was invaded by trees and tall shrubs than alpine tundra.

Physiographic controls on tree abundance
The best multinomial model of 2015 vegetation type (trees, arctic tundra and alpine tundra) was a model that was quadratic for elevation and slope, linear for northness and included an interaction of slope and northness (Supplementary material Appendix 1 Table A4, Fig. 6). The overall accuracy of the multinomial model fit was 0.71 in all three years (Supplementary material Appendix 1 Table  A5-A7). Trees were most abundant at lower elevations (< 200 m) on steep (~20°) south-facing slopes. Alpine tundra was most abundant at high elevations (> 200 m) on steep (> 20°) north-facing slopes, while arctic tundra was most abundant at mid-elevations (~200 m) in areas with shallow slope (< 20°), especially on north-facing aspects. Spatial plotting of the multinomial model predictions for 1952, 1979 and 2015 revealed a clear pattern of treeline advance out of the drainages and into upland arctic and alpine tundra types (Supplementary material Appendix 1 Fig. A3).

Discussion
Researchers have long related global treeline elevations to isotherms, suggesting that temperature or variables related to temperature, may explain treeline position. Daubenmire (1954) and references therein suggested that global timberline elevation is associated with the 10°C isotherm for the warmest month. More recently, Körner and Paulsen (2004) concluded that global treelines occur within a narrow range of growing season air temperature (6-7°C), which they argued equals growing season shallow soil temperature, despite acknowledged differences in seasonality between air and soil temperatures. Based upon an observed lapse rate of −4.6°C km −1 in the Agashashok watershed, the Daubenmire (1954) isotherm would occur at 800 m, while the Körner and Paulsen (2004) isotherm would occur at 1100-1200 m in our study area -elevations that are well above the highest peaks (maximum 600 m a.s.l.). Our multinomial model fit to the 2015 random point data revealed an inflection point at an elevation of 202 m, where there is an equal probability of finding a tree or tundra (a quantitative method of defining treeline elevation, Dial et al. 2016). Between 2012 and 2017, mean July temperature at 200 m was 12.7°C, while mean June-August air temperature was 11.2°C. If a mean July air temperature of 10°C or a mean growing season air temperature of 6-7°C were good predictors of treeline elevation in this region, then the study area would likely be almost completely forested, with the exception of areas that are too steep and where parent material movement is too great to allow for While air temperature in our study area is relatively warm by treeline standards, the soils are quite cold, owing to the presence of discontinuous permafrost and/or impeded drainage. For instance, at an arctic treeline on the west side of the Agashashok River, where the forest grades into tussock tundra, Sullivan et al. (2015) reported growing season mean temperatures of 11.6°C for the air at 2 m and 4.9°C for the soil at 10 cm depth. This large difference between air and soil temperatures is inconsistent with the assertion of equality between growing season air and soil temperatures and suggests the lower than expected treeline elevation, when based on growing season air temperature, could be related to cold soils. The notion that treeline elevation in our study area is related to soil rather than air temperature is supported by our observation that arctic treelines, where the forest grades into the colder soils of mesic and wet tundra, occur at much lower elevations than alpine treelines, where the forest interfaces with the comparatively warm soils of dry tundra (Supplementary material Appendix 1 Fig. A4). The occurrence of arctic treelines in areas with modest slope, where growing season air temperature likely varies little across the treeline, also points to the greater importance of soil, rather than atmospheric conditions. Recent work in the Agashahok watershed and throughout the Brooks Range has revealed that growth of treeline white spruce is strongly nutrient-limited (Sullivan et al. 2015, Ellison et al. 2019. We hypothesize that tree abundance may be similarly limited by the effects of cold soils on tree access to nutrients. Cold soils could limit tree access to nutrients through reduced fine root productivity, limited mycorrhizal colonization/productivity, reduced fine root lifespan and/or limited microbial processing of soil organic matter. It may be that climate warming has not yet fully lifted these constraints, as much of the additional energy is expended in the evaporation of water and thawing of ground ice. Thus, one component of the time lag between atmospheric warming and treeline rise may be the disparity between air and soil temperatures. Our analysis of vegetation change for > 3000 random points revealed two parallel and nearly mutually exclusive pathways of vegetation change. Arctic and alpine tundra were almost equally likely to transition to tall shrubs or open woodland. However, once established, open woodland either remained open woodland or transitioned to forest, while tall shrubs either remained tall shrubs or transitioned to trees with shrubs. Open woodland appeared to be a transitional state, with little change in abundance over time. There were very few instances of tall shrubs or trees with tall shrubs transitioning to forest. This suggests that, in the absence of disturbance, the establishment of tall shrubs may lead to an alternative 'climax' community, composed of tall shrubs with scattered trees, that was almost completely absent from the landscape in 1952 (1% of random points).
Between 1952 and 2015, our 5457 ha study area in the Agashashok watershed gained 42 ha of tall shrubs, 29 ha of trees with tall shrubs and 73 ha of forest. Meanwhile, 70 ha of arctic and 23 ha of alpine tundra were lost. These vegetation changes have probably affected ecosystem-atmosphere fluxes of C, H 2 O and energy, while changing the quality and quantity of wildlife habitat. The additional leaf area maintained by tall shrub and forest ecosystems has likely led to an increase in growing season C uptake (LaFleur et al. 2001), while the accumulation and preservation of deeper winter snow in the taller vegetation has probably led to warmer soils and greater overwinter C efflux (Sullivan 2010). The lower albedo associated with more forested areas has likely increased absorption of shortwave radiation and led to local heating that has accentuated the climate warming trend (Beringer et al. 2005). Encroachment of tall shrubs and trees into alpine tundra is probably not favorable for Dall's sheep, which depend upon low-statured dry heath tundra (Rachlow and Bowyer 1998). A woodier landscape is also unfavorable for caribou, which avoid forest and tall shrubs during migration (Fullman et al. 2017) and use habitat similar to Dall's sheep during times of insect harassment (Walsh et al. 1992). While many of the tall shrubs that have increased in abundance in our study area are alders, some are willows and an increase in their abundance is favorable for moose, which are new to the Agashashok watershed, having moved in around 1950 (Tape et al. 2016). Comparison of the current treeline elevation with the global mean isotherm of Körner and Paulsen (2004) suggests that there is considerable potential for future increases in woody vegetation within the Agashashok watershed. Even if the climate warming trend were to cease, we would still expect to see further increases in tree abundance, as many seedlings and saplings that have recently established, but are not yet visible in satellite imagery, will likely grow into mature trees and have increasing effects on ecosystem function.

Data availability statement
Data presented in this manuscript will be archived in the Arctic Data Center (<https://arcticdata.io/>).
Permits -The authors thank the United States National Park Service for permission to work in Noatak National Preserve under permit #NOAT-2017-SCI-0001.