2. Materials and Methods
2.1. Sample Collection and Preparation
We collected tree-ring samples from living P. engelmannii trees at four sites (FLS, RPS, MWS, and SLS) across the GYE in northwest Wyoming, US (44.4–45.0°N, 109.5–110.5°W), two located in Yellowstone Na- tional Park and two on the Shoshone National Forest (Figure 1a and Table S1). At each site, we collected in- crement cores from living trees (2 per tree) (Stokes & Smiley, 1968), and at the FLS site, we used a chainsaw to remove cross sections from downed remnant, subfossil P. engelmannii logs. All cores and cross sections were processed using standard dendrochronological techniques (Stokes & Smiley, 1968). Because samples did not exhibit visible color changes between heartwood and sapwood or between living and subfossil sam- ples, no chemical processing was required prior to collecting LWB measurements.
2.2. Chronology Development
We obtained LWB measurements for all samples using an Epson Expression XL 12000 commercial flatbed scanner in combination with the software CooRecorder (Larsson, 2014), the absolute dates of which were then validated statistically with the software COFECHA (Holmes, 1983). We inverted the raw LWB values to allow for detrending by methods similar to those used for MXD data (Rydval et al., 2014). We removed age-related trends by detrending all LWB series using the age-dependent spline using the power trans- formed residuals in the SignalFree framework (Melvin et al., 2007).
2.3. Climate Response and Reconstruction Model
We compare each site LWB chronology with monthly CRU TS 4.03 (land) 0.5Tmax and Tmean data (Harris et al., 2014) averaged over the GYE study region (44.00–44.50° N, −109.00 to 109.50°W) for the period 1905–2015. We use Pearson correlation analysis to test temporal (e.g., target reconstruction season) and spatial relationships between temperature and LWB data using KNMI Climate Explorer (Trouet & Van Old- enborgh, 2013). We then employ a nested principal components regression (PCR) model, which accounts for the decrease in number of predictor chronologies back through time (E. R. Cook et al., 1999). LWB chronologies are retained for modeling if they are positively and significantly (p < 0.01) correlated to the target instrumental data during the common period. To validate/cross-validate the reconstruction model, we split the common period into two equal periods for calibration on the early period and verification on the late period. To assess model uncertainty, we use the maximum entropy bootstrapping method (MEBoot) (Vinod & López-de Lacalle, 2009) to produce 300 reconstruction replicates. We conduct a severity–duration analysis (González & Valdés, 2003) on reconstructed values of at least 2+ years to determine the magnitude (average cumulative departure from the long-term mean), duration (time span of each event), and intensity (duration divided by the magnitude) of warm and cool events. The top 10 warm and cool events are based on the overall rank (duration rank + magnitude rank).
3. Results and Discussion
3.1. Chronology Development
Of the four GYE LWB chronologies, FLS is the only record that spans >1,000 years (737–2016 CE) due to the successful inclusion of remnant material into the chronology (Table S1). RBAR values indicate that all LWB chronologies show adequate signal coherence, ranging from 0.38 (FLS) to 0.29 (MWS) over the period during which EPS > 0.80. RBAR values for these chronologies are similar to those found in other BI studies (Heeter et al., 2020; Wilson et al., 2014, 2019). Correlations between each site LWB chronology indicate a high level of similarity (Figure S1), and thus all four LWB chronologies advanced to the candidate pool of predictors for the PCR model.
3.2. Climate Response and Seasonality
We found the strongest significant and positive relationships between all four LWB chronologies and cur- rent-year August Tmax (r = 0.43–0.76, p < 0.001). This is consistent with other studies, which suggest that LWB typically shows stronger relationships with Tmax than with Tmean (e.g., Heeter et al., 2020; Wilson et al., 2014, 2019). Because the LWB/instrumental data relationships were consistently strong after ca. 1920, we use 1920 to be the starting year for our subsequent analysis. We find the climate response of all four sites to be representative of Tmax variability across the GYE and broader western US (Figure S2). Across all four sites, we note a positive relationship between site elevation and temperature signal strength; the highest elevation sites yielded the strongest spatial and temporal temperature responses (Table S1). Wil- son et al. (2019) and Salzer, Larson, et al. (2014) show a similar relationship between site elevation and the signal strength of LWB chronologies across the Yukon, Canada, and TRW chronologies in the White Mountains of California, USA, respectively. We limit our reconstruction target to August Tmax due to a combination of several factors. First, recent studies have shown the growing season for high-elevation P. engelmannii, particularly at the FLS site, is restricted to late July through early September, such that late- wood tracheids are most likely differentiating and maturing during the month of August (Rochner, 2019). Second, the LWB metric is shown to be most strongly representative of summer temperature variability (Rydval et al., 2014; Wilson et al., 2014; Wilson, D'Arrigo, et al., 2017), and at the 4 GYE sites, August temperatures appear to be the most limiting factor on densiometric growth. Heeter et al. (2020) document similar seasonal reconstructions (e.g., August–September Tmax) due to shortened growing seasons at high elevations in the Southwest US.
3.3. GYE Tmax Variability 770–2019 CE
We provide a 1,250-year reconstruction of August Tmax that accounts for 60% of the instrumental period variance (Figure 1). Beta values for the four LWB chronologies (Table S1) indicate that all chronologies are retained for the reconstruction model. Over the common period (1822–2015), our model strength (r2 = 0.60) is comparable to other temperature reconstructions for the western US (e.g., MXD-derived Briffa et al., 1992 r2 = 0.40 at grid point #14 and TRW-derived Salzer, Bunn, et al., 2014 r2 = 0.46) (Figure S3). Validation/ cross-validation statistics remain robust for each forward (n = 2) and backward (n = 3) nest spanning the period 770–2019 and ensure that the reconstruction is time-stable with high explanatory power (Table S2). Within the long-term context, the reconstruction exhibits multidecadal to centennial-scale temperature variability from 770 to the middle twelfth century, after which cool temperatures dominate until the middle fifteenth century (Figure 2a). Prolonged warm temperatures start ca. 1450 and persist through the middle seventeenth century. Temperature conditions over next ca. 250 years are characterized by low variability until the turn of the twentieth century.
3.3.1. Notable Warm and Cool Events
The intense warming trend since ca. 2000 CE far exceeds conditions during any other period over the past 1,250 years, particularly during the MCA. The severity–duration analysis indicates that 2016 was the warmest single-year anomaly since 770, followed by the year 1056, which was centered during the MCA (Figures 2b, 2c, and 3a; Table S3). Yet across the GYE, the MCA is characterized by substantial multidec- adal summer temperature variability rather than one prolonged homogeneous warming event (Figures 2a and 4a). Within the time frame of the MCA, however, is the top-ranked warm event (based on intensity, duration, and magnitude) lasting ca. 20 years from 1050 to 1070. Although the period 2015–2019 is the tenth-ranked event, the warming over the past 5 years is twice as intense as warming during the MCA; hence, the overall ranking of the 2015–2019 event is lower due to the short duration. As 2020 rivaled 2016 for the highest global surface air temperatures on record (Voosen, 2021), current and impending projections of increasing anthropogenic warming (IPCC, 2013), especially for the midlatitudes (Li et al., 2012), suggest that the modern warming will likely soon surpass the MCA (as defined by the severity–duration analysis).
The twentieth century contains two noteworthy warming anomalies: the year 1988 and the 1931–1940 pe- riod (Figure 3a and Table S3). The fourth-highest-ranking warm year since 770 CE, 1988, coincides with the timing of the 1988 Greater Yellowstone Area fires, which collectively formed the largest complex fire recorded for the GYE (Christensen et al., 1989; Romme & Despain, 1989). Multiple studies suggest that climate, rather than fuels, was the primary driver of the 1988 Complex (Romme & Turner, 2004; Turner et al., 2003) and that continued warming could push GYE fire regimes past a tipping point, to where the system is qualitatively altered, by the middle 21st century (Westerling et al., 2011). The GYE summer tem- perature reconstruction characterizes 1988 as highly abnormal, with a return interval in excess of 500 years (Figure 2c). Our reconstruction provides long-term contextualization of extreme warm year abnormalities like 1988, which alludes to the increased role of temperature as a driver of fire activity in recent decades. As the influence of warming temperatures on aridification and increased fire activity across many parts of the western US continues to be explored (Abatzoglou & Williams, 2016; Westerling et al., 2011), historical con- textualization of modern climatic trends (e.g., identifying historical analogs) is important for developing a more comprehensive understanding of consequences of future climatic change for the GYE. The 1931–1940 period ranks as the fifth-warmest event and coincides with the Dust Bowl, a well-documented period of prolonged drought across much of the western US (B. I. Cook et al., 2008, 2009). The spatial imprint of warming during the 1930s is extensive, as this event can be seen in both regional and broader temperature reconstructions for the NH (Figure 4d).
As we found with the MCA, our reconstruction provides further evidence that the Little Ice Age (LIA; ca. 1300–1850) was characterized by multidecadal variability, rather than by a solitary, sustained cool period (Briffa et al., 1992; Salzer & Kipfmueller, 2005; Wise, 2014) (Figure 2a). The Maunder Minimum (MM; ca. 1645–1715; Eddy, 1976) is particularly pronounced, with the second highest-ranked cool event occurring from 1697 to 1709 (Figure 2b and Table S3). This cool event is strongly and negatively correlated (p < 0.01) with the highest period of snowpack anomalies (ca. 1695–1735 CE) for the GYE presented by Pederson et al. (2011) (Figure 4c). We also show evidence of the Dalton Minimum (DM; 1790–1830, Wagner & Zori- ta, 2005), characterized by more moderately cool, late-season Tmax and high cool-season snowpack.
3.3.2. Broad-Scale Forcings on GYE Temperature Variability
The period ca. 1500–1640 is the longest-sustained, multidecadal-scale warming in the reconstruction, and it aligns with the widely documented sixteenth century Megadrought over much of North America (E. R. Cook et al., 2007; Stahle et al., 2000), reduced precipitation within the Upper Missouri, Columbia, and Colorado River basins (Gray et al., 2007), anomalously frequent wildfire activity across the GYE (Whit- lock et al., 2008), and a ridge pressure pattern over the region, as demonstrated in the geopotential height reconstructions by Wise & Dannenberg (2014). As we note a strong, positive relationship between Tmax and upper air pressure (Figure S4), we posit that this ca. 140-year warm event, initially brought about by increased upper air pressure, caused prolonged drought and increased fire across the GYE.
While various decadal to multidecadal cooling events align with well-documented periods of reduced solar activity (e.g., Spörer Minimum; Figures 4b and S5), volcanic forcing likely also enhanced the expression of prolonged cool periods such as the MM and DM. We find strong and significant (p < 0.01; n = 16) depar- tures in GYE summer temperatures following many of the largest volcanic eruptions since 1100 CE. The SEA shows significant reductions in August Tmax 1 (p < 0.01) to 2 (p < 0.05) years following major volcanic eruptions (Crowley & Unterman, 2013; Sigl et al., 2015) (Figure 4b). Based on the lag value denoted by SEA, the ninth-ranking coldest single-year anomaly, 1993 (Table S3), can be partially explained by the ef- fects of the June 1991 eruption of Mount Pinatubo, Philippines (Figure 4a). Although both periods show pronounced reductions in solar irradiance, most reconstructions show solar activity to be lower during the MM compared to the DM (Hoyt & Schatten, 1993; Lean et al., 1995). The MM coincides with the 1641 Parker and 1696 Serua eruptions, thereby demonstrating the influence of both solar and volcanic forcing.