Forests on the move: Tracking climate-related treeline changes in mountains of the northeastern United States

Remote sensing analysis Physical copies of true color high resolution historical aerial imagery (sub-meter resolution) were acquired from the Appalachian Mountain Club (AMC) and the USFS White Mountain National Forest Headquarters. Imagery for the Presidential Range was taken in 1978 and Katahdin imagery was taken in 1991. Hard copy images were scanned and converted to TIFF format at 300 dpi (resulting in 0.5 m resolution images). Spatial analyses of change in treeline positions over time were enabled by acquiring high resolution 2018 false-color near-infrared imagery from the National Agriculture Inventory Program (NAIP 2021). Both sets of imagery were taken during summer months (1:40,000 scale). Using ArcGIS 10.8 (ESRI 2011, Redlands, CA, USA), historic imagery was ortho- and georectified to newer imagery via a spline function along 60 ground control points, and then converted into one orthomosaic image (RMSE < 1m). Exact error was always below 5 m for each individual image. All areas above treeline were manually digitized based on observed tree cover for both sets of images, and the resulting polygons were converted to raster format at 2 m resolution (all raster pixels within each polygon had a value of 1). We identified forest cover only as areas with overlapping crowns and seen as green reflectance in historic imagery and red reflectance in contemporary false-color near-infrared imagery (no visible bare earth or easily identified alpine vegetation). Isolated tree island edges were also digitized and included as treeline if they were >20 m in diameter in any direction (determined in ArcGIS) and included an individual >2 m in height as validated in the field. Alpine rasters were aligned to and multiplied by Lidar-derived digital elevation models (DEMs; 2 m resolution) acquired from New Hampshire and Maine state GIS repositories in order to determine treeline elevations. A total of 400 random sample points (200 for each range, using the ArcGIS random sample point tool) were placed along the outer boundary of the alpine rasters derived from our contemporary imagery, and for each of them we established a paired point at the nearest location along the alpine raster boundary derived from our historic imagery. Field surveys Field sampling was carried out in the summer of 2021 to characterize tree demography and demographic variation among different treeline forms identified from the current imagery. A subset of contemporary points from our GIS-based sample point pairs (n = 54, 33 in the Presidential Range, 21 in the Katahdin Range, see above) were selected using a random number generator to serve as sites for establishing belt transects. Each belt transect was 100 m in length and 4 m wide (2 m on either side of transect for a total area of 400 m2) and perpendicular to elevation contours, spanning the ecotone between closed forest interior and open alpine habitat. The start of each transect (the lowest elevation on the transect, set as 0 m) was located 50 m downslope (straight-line distance) of contemporary sample points. The start and end of each belt transect were recorded using a Garmin GPSMAP 64 (Garmin, Olathe, Kansas, USA). Each tree > 0.1 m in height with a stem rooted within the transect was recorded noting species, basal diameter (10 cm from the ground), height, horizontal distance from the transect, and distance along the transect (to estimate stem density of trees). Slope, aspect, elevation, and soil depth to bedrock (using a metal soil probe) were recorded at 20 m intervals along the belt transect centerline (0 m, 20 m, 40 m, 60 m, 80 m, 100 m). For all belt transects, treeline form was assigned based on visual assessments (based on changes in tree height and density across the ecotone). Additionally, we visited a majority of our other accessible contemporary random sample points (~80%) in order to assign treeline form and ground-truth remote sensed treeline classifications. For all visited sample points we took a new GPS point at the field-verified treeline location (continuous canopy cover and at least one individual >2 m in height) nearest to our random sample points (assigned from our treeline delineation procedure). The new points were compared to the original sample point locations and assessed for accuracy (measuring linear distance between points). Eye-level photos of treelines were taken at all sample points to keep a permanent record of treeline appearance. We stress that because tree height could not be extracted or field validated from our historic imagery, some krummholz individuals (<2 m) may have been present above our treeline delineation using our classification scheme. Out of all 400 sample point pairs across both the Presidentials and Katahdin, 88 were classified as abrupt (22%), 70 as diffuse (17.5%), 84 as island (21%), and 162 as krummholz (40.5%). Spatial data processing To examine the factors potentially influencing the spatial dynamics of treeline advance, both climatological and topographical variables were extracted for the Presidential Range. We could not conduct a similar analysis for Katahdin given the lack of fine-scale climatological data in that area. Elevation was extracted from 2 m state produced DEMs. Using the Spatial Analyst toolbox in ArcGIS, topographical variables such as slope, aspect, and curvature (measure of convex or concave shape of the terrain ranging between -4 and 4) were extracted from our DEMs. Circular aspect data (measured in degrees, 0-360⁰) were converted to radians and linearized (east and west = 1, north and south = 0). Before linearization, aspect values were used to calculate degree difference from prevailing wind (DDPW - 290˚) and degree difference from south (DDS - 180˚) variables. DDPW is a proxy for exposure to strong winds that can cause both direct physical damage and damage from icing, as well as a proxy for the potential for snow accumulation. The prevailing wind direction for the Presidential range (290˚) was based on wind measurements from the Mount Washington Observatory. DDS is a proxy for the amount of direct solar radiation (in the northern hemisphere). Average monthly mean, maximum, and minimum temperatures as well as annual accumulated growing degree days (AGDD) were calculated from an array of 34 HOBO dataloggers (Onset Computer Corporation, Bourne, MA, USA) placed at various elevations and adjacent to Appalachian Mountain Club buildings in the White Mountains of New Hampshire. HOBO loggers have recorded hourly air temperature at ground level (0 m height) continuously since 2007. Air temperature means and AGDD were calculated from HOBO logger data; for AGDD calculations we used a base temperature of 4˚C, consistent with other studies examining growth patterns of balsam fir, the dominant species within studied treelines. AGDD was calculated as the accumulated maximum value of growing degree days (GDD) in a year. Gridded maps (90 m spatial resolution) of mean annual temperature (Tmean, between 2007 and 2020) and AGDD for the Presidential Range region were produced using a cokriging interpolation method. To do this, temperatures and AGDD response variables were first checked for normality using qq-plots. Next, correlation between response variables and potential covariates was assessed; both elevation and aspect were highly correlated with HOBO derived temperature and AGDD. We used normal-score simple cokriging with a stable semi-variogram model to interpolate (prediction map) climate variables over the entire spatial extent of the Presidential Range (RMSE ~ 1 for both Tmean and AGDD). Mean annual precipitation was estimated from 30-year normal PRISM climate data (1991-2020; PRISM Climate Group, Oregon State University, https://prism.oregonstate.edu).遥感分析 本数据集获取了来自阿巴拉契亚山脉俱乐部（Appalachian Mountain Club, AMC）与美国林务局（United States Forest Service, USFS）怀特山国家森林总部的历史真彩色高分辨率航空影像硬拷贝（亚米级分辨率，sub-meter resolution）。其中，总统山脉（Presidential Range）的影像拍摄于1978年，卡塔丁山（Katahdin）的影像拍摄于1991年。将硬拷贝影像扫描后以300dpi分辨率转换为TIFF（Tagged Image File Format）格式，最终得到0.5米分辨率的影像。为开展林线位置随时间变化的空间分析，本数据集还获取了来自国家农业调查计划（National Agriculture Inventory Program, NAIP 2021）的2018年假彩色近红外影像。两组影像均拍摄于夏季，比例尺为1:40000。 使用ArcGIS 10.8（ESRI 2011，美国加利福尼亚州雷德兰兹），通过沿60个地面控制点的样条函数，将历史影像进行正射校正与地理配准，使其对齐至当代影像，随后拼接为一幅正射镶嵌图（均方根误差Root Mean Square Error, RMSE < 1米），单幅影像的绝对误差始终低于5米。 基于两组影像中观测到的树木覆盖情况，手动数字化所有林线以上区域，将生成的多边形以2米分辨率转换为栅格格式（每个多边形内的所有栅格像素值均为1）。森林覆盖仅定义为冠层重叠、在历史影像中呈绿色反射、在当代假彩色近红外影像中呈红色反射的区域（不含裸露地表或易识别的高山植被）。孤立的树岛边缘也被数字化并纳入林线范畴，前提是其任意方向直径大于20米（在ArcGIS中测定）且包含至少一株经野外验证高度大于2米的个体。 将高山栅格与从新罕布什尔州和缅因州地理信息系统仓库获取的激光雷达（Light Detection and Ranging, Lidar）衍生数字高程模型（Digital Elevation Model, DEM，2米分辨率）对齐并相乘，以计算林线海拔。基于当代影像生成的高山栅格外边界，共放置400个随机采样点（每个山脉200个，使用ArcGIS随机采样点工具），并为每个点在历史影像生成的高山栅格外边界的最近位置设置配对点。 野外调查 2021年夏季开展野外采样，以表征树木种群动态以及从当前影像中识别的不同林线类型间的种群动态差异。使用随机数生成器选取GIS采样点对中的一个子集（n=54，总统山脉33个，卡塔丁山脉21个，见前文）作为样带设置点。每条样带长100米、宽4米（样带两侧各2米，总面积400平方米），垂直于等高线，跨越闭合森林内部与开放高山生境的交错带。每条样带的起点（样带最低海拔，设为0米）位于当代采样点下方50米直线距离的下坡处。使用Garmin GPSMAP 64（美国堪萨斯州奥拉西市Garmin公司）记录每条样带的起点和终点。记录样带内所有株高>0.1米、茎干扎根于样带内的树木，记录其物种、基径（距地面10厘米处）、树高、距样带的水平距离以及沿样带的距离（用于估算树木茎干密度）。沿样带中心线以20米间隔（0米、20米、40米、60米、80米、100米）记录坡度、坡向、海拔以及基岩以上土壤深度（使用金属土壤探针）。所有样带的林线类型基于目视评估（基于交错带内树高和密度的变化）进行划分。此外，我们访问了大部分其他可到达的当代随机采样点（约80%），以划分林线类型并实地验证遥感林线分类结果。对于所有访问过的采样点，在距离随机采样点最近的、经野外验证的林线位置（连续冠层覆盖且至少有一株高度>2米的个体）采集新的GPS点。将新点位与原始采样点位置进行比较并评估精度（测量两点间的直线距离）。在所有采样点拍摄林线的平视照片，以留存林线外观的永久记录。 需要说明的是，由于无法从历史影像中提取或实地验证树高，部分矮曲林（krummholz）个体（<2米）可能存在于本分类方案划定的林线以上区域。在总统山脉和卡塔丁山脉的全部400个采样点对中，88个被划分为突然型林线（abrupt，22%）、70个为扩散型林线（diffuse，17.5%）、84个为岛状型林线（island，21%）、162个为矮曲林型林线（krummholz，40.5%）。 空间数据处理 为探究可能影响林线推进空间动态的因素，本数据集提取了总统山脉的气候与地形变量。由于卡塔丁地区缺乏精细尺度的气候数据，无法开展类似分析。海拔从2米分辨率的州级DEM中提取。使用ArcGIS空间分析工具箱，从DEM中提取地形变量，包括坡度、坡向以及曲率（地形凹凸形状的度量，取值范围为-4至4）。将以度为单位的圆形坡向数据（0-360°）转换为弧度并进行线性化处理（东和西=1，北和南=0）。线性化之前，使用坡向值计算与主导风向（DDPW，290°）的角度差以及与南向（DDS，180°）的角度差变量。DDPW是暴露于强风的代理变量，强风可造成直接物理损伤和结冰损伤，同时也是积雪潜在积累的代理变量。总统山脉的主导风向（290°）基于华盛顿山天文台的风速测量数据确定。DDS是直接太阳辐射量的代理变量（北半球）。 从2007年起持续在地面高度（0米）记录小时气温的34台HOBO数据记录仪（Onset Computer Corporation，美国马萨诸塞州伯恩市），放置于新罕布什尔州怀特山脉的不同海拔高度且邻近阿巴拉契亚山脉俱乐部建筑。从HOBO记录仪数据中计算气温平均值和年度积温（Accumulated Growing Degree Days, AGDD）；AGDD的计算采用4℃的基准温度，这与其他研究中研究林线优势种香脂冷杉（balsam fir）生长模式的基准一致。AGDD计算为年度内生长度日（Growing Degree Day, GDD）的累积最大值。使用协同克里金插值法生成总统山脉区域2007-2020年的年平均气温（Tmean）和AGDD的90米空间分辨率栅格图。具体步骤为：首先通过QQ图检验气温和AGDD响应变量的正态性；随后评估响应变量与潜在协变量间的相关性，发现海拔和坡向与HOBO记录仪获取的气温和AGDD高度相关。我们使用带稳定半变异函数（semi-variogram）模型的正态得分简单协同克里金法，在总统山脉的整个空间范围内插值得到气候变量（Tmean和AGDD的RMSE约为1）。年平均降水量从30年平均PRISM气候数据（1991-2020；PRISM气候组，俄勒冈州立大学，https://prism.oregonstate.edu）估算得到。