The effects of climate change on the timing of peak fall foliage in Acadia National Park

In recent decades, most United States’ National Parks have experienced extreme temperature and precipitation regimes outside of their historical ranges of variability with unknown effects on fall phenology. We determine 1) how seasonal climate in Acadia National Park, Maine and the timing of peak fall foliage have changed between 1950 and 2021; 2) how changes in seasonal climate have affected fall foliage; and 3) how we might expect the timing of fall foliage to change given future climate projections. We use ERA5-Land data to analyze changes in climate. We couple remotely sensed data with archival research to determine changes in the timing of peak fall foliage. We use multivariate regressions to understand the relationship between climate and fall foliage. We use CESM2 data to predict the timing of peak fall foliage coloration through 2060. Minimum temperatures, maximum temperatures, precipitation, and the number of warm nights, hot nights, warm days, hot days, and downpour days have all significantly increased (p ≤ 0.05). The timing of peak fall foliage is now occurring almost two weeks later (p ≤ 0.05). September temperature and precipitation and May precipitation were positively correlated with delayed peak fall foliage. Early October precipitation was negatively correlated. Future climate projections predict the timing of peak fall foliage to occur between October 30 and November 2 by 2060. Understanding how climate is affecting leaf senescence both is crucial in a national park where fall tourism brings large gains to the local economy and provides key information to park managers planning for a resilient, sustainable future.

in visitation, particularly during the spring and fall shoulder seasons. Importantly, though, many National Parks and the businesses in the surrounding areas hire seasonal staff focused on the busiest summer months, and close facilities during the late fall, winter, and early spring. Interviews with national parks experts have highlighted climate change and tourism as the two biggest forthcoming challenges in park management (Feretti-Gallon et al. 2021). Thus, specifically understanding the effects of climate change on seasonality in the National Parks is essential for park managers, visitors, and those living and working in the surrounding towns. Here, we will specifically focus on the effects of climate change on Acadia National Park, Maine ( Fig. 1) in the fall.
The National Parks of the United States disproportionately experience the effects of climate change compared to the rest of the United States (Hansen et al. 2014;Gonzalez et al. 2018). Almost all National Parks but those in the southeast warming hole (Partridge et al. 2018) are experiencing more days with temperatures above their historical 95th percentile (Monahan and Fisichelli 2014). Spring is occurring earlier in over 75% of the National Parks, with over half experiencing springs that have advanced outside of the 95th percentile of historical  conditions (Monahan et al. 2016). And, although a precipitation signal across all National Parks is less clear, regional trends exist where parks in the northeast are getting wetter while those in the southwest are getting dryer (Monahan and Fisichelli 2014;Gonzalez et al. 2018).
Importantly, these changes in climate affect the numbers and timing of park visitation. A study on over 340 parks determined that warmer temperatures are associated with more park visitors, up to an inflection point: visitation decreases when temperatures rise above 25 °C (Fisichelli et al. 2015). Precipitation is, again, more complex. In the United States, precipitation was not correlated with increases or decreases in park visitation, but instead, associated with significant changes in where park visitors go once there: on rainier days, visitors stay closer to infrastructure once in the park (Wilkins et al. 2021). However, recent recreation-climate-change and public landsvisitorship studies have demonstrated that the spatial Fig. 1 Acadia National Park boundary is in green and Mount Desert Island, the specific focus of this work, is outlined in dark blue. The Schoodic Peninsula and Isle au Haut sections of Acadia National Park are also highlighted in green in the northeast and southwest, respectively. The McFarland Hill Global Historical Climatology Network daily (GHCNd) temperature station is highlighted with the orange circle, and two Southwest Harbor GHCNd precipitation stations are highlighted with purple squares. Inset. Red box showing extent of map, and location of Acadia National Park, relative to Maine, USA and temporal scale of interest matters. Changes in visitation vary seasonally, regionally (Wilkins et al. 2021), and even between parks within the same state. A Utah-focused National Parks study found that in three of the five parks studies (Arches, Canyonlands, and Capital Reef), visitation increased with temperature until maximum daily temperatures exceeded (25 °C), yet this visitation-temperature threshold was not evident in Bryce Canyon and Zion (Smith et al. 2018). And, work has shown that winter and spring visitation to parks in Canada has not changed, while autumn visitation is increased by warmer temps and decreased by precipitation (Hewer and Gough 2019).
Autumn has been referred to as "the neglected season in climate change research" . However, understanding the effects of climate change on autumn phenology is critical for regions like New England, where fall tourism ("leaf peeping") is a billion-dollar industry (Archetti et al. 2013). In 2021 alone, Maine's National Parks and Monuments, which include Acadia National Park, the Katahdin Woods and Waters National Monument, and the Saint Croix Island International Historic Site, added $770 million in economic output to the Maine Economy (Thomas et al. 2022). Acadia National Park, which draws a majority of Maine's National Parks and Monuments visitors, has experienced the greatest increases in visitation since 2000 in the fall months of September and October (IRMA 2022), Fig. 2). Moreover, fall foliage tourism specifically comprises 20-25% of total annual tourism visitation to Maine itself (Pratt-Kelley 2022). Thus, understanding how climate change affects fall foliage in Acadia National Park is an essential first step in understanding the relationship between climate change and park visitation.
Whether because of its complex drivers, the protracted nature of events, a lack of standardized evaluation methods, or simply being considered a less 'exciting' season, less than a third of phenology research has focused on autumn, while a majority of phenology research has focused on the spring season . Temperature and photoperiod are generally considered the two primary controls to determine the timing of plant senescence (Jolly et al. 2005;Richardson et al. 2013;Gill et al. 2015;Fu et al. 2018), but fall temperatures are not as correlated with senescence as spring temperatures are with leaf-out and bud burst (Menzel 2003). There is no consensus in the fall phenology literature on whether to expect advanced or delayed leaf senescence and coloration with changes in climate (Delpierre et al. 2009;Gordo and Sanz 2010;Keenan and Richardson 2015;Gill et al. 2015;Liu et al. 2016;Chen et al. 2020;Wu et al. 2022), likely because fall phenological responses to climate stressors tend to be species and location specific (Archetti et al. 2013; Xie et al. 2018a, b). However, remote-sensing and field-based studies correlate delayed onset of leaf senescence with warmer spring, summer, and autumn temperatures, and decreased autumn extreme precipitation events (Delpierre et al. 2009;Ibenez et al. 2010;Gunderson et al. 2011;Archetti et al. 2013;Fu et al. 2018, Chen et al 2020, and correlate earlier onset of fall colors with heat stress, increased autumn precipitation, drought, early frosts, and elevated nitrogen concentrations (Norby et al. 2003;Chmielewski et al. 2005;Leuzinger et al. 2005;Fu et al. 2014;Gallinat et al. 2015;Xie et al. 2015;Signarbieux et al. 2017;Zani et al. 2020;Chen et al. 2020).
Here, we seek to answer four questions: (1) To what extent has seasonal climate changed in Acadia National Park between 1950 and 2021; (2) To what extent has the timing of peak fall foliage changed in Acadia National Park during the same time period; (3) To what extent have changes in seasonal climate affected the timing of peak fall foliage; and (4) How might we expect the timing of peak fall foliage to change given future climate projections. We use ERA5-Land reanalysis data to investigate Acadia's changing climate. We couple remotely-sensed data with archival research to determine how the timing of peak fall foliage has changed. We use multivariate regressions to understand the relationship between changes in Acadia's climate and the timing of peak fall foliage. And, we predict the timing of peak fall foliage coloration through 2060 using Community Earth System Model (CESM2) shared socioeconomic pathway (SSP) projection data (O'Neill et al. 2017). Our results provide insight into understudied yet important seasonal dynamics in one of the most visited National Parks in the United States and will directly assist park visitors and managers adapt to ongoing change.

Acadia National Park, Mount Desert Island
Mount Desert Island is home to over 85% of Acadia National Park. We focus our analysis on Mount Desert Island (MDI) for two reasons. Firstly, a majority of the deciduous trees in Acadia are found on MDI (Fig. 3), where The Great Fire of 1947 burned over 40 km 2 of spruce and fir trees-including nearly one-third of the total park area on Mount Desert Island, and maples, birch trees, and aspens have grown in their place (NPS 2020). Secondly, in 2021, over 90% of the Acadia visitors remained on MDI; less than 1% visited Isle au Haut and the rest visited the Schoodic Peninsula (Broom 2021). We also extend our time period of analysis to the pre-satellite data era from 1950 to 2021 to facilitate robust statistical analyses.

Climate data and analysis
Because there is no consistent long-term observational record on MDI (see SI for MDI GHCNd station information), we downloaded 9 km resolution global land-based reanalysis ERA5-Land hourly precipitation and 2 m air temperature data from January 1, 1950-December 31, 2021(Muñoz-Sabater et al. 2021, and derived daily total precipitation, and minimum and maximum temperatures. ERA5-Land data is improved from the European Centre for Medium-Range Weather Forecasts' (ECMWF) fifth version of European ReAnalysis (ERA5) as it focuses solely on land, includes an elevation correction, and offers higher-resolution data (Muñoz-Sabater et al. 2021). Comparisons with GHCNd data from stations on MDI demonstrate that ERA5-Land data provide robust representations of weather and climate dynamics across the study area (SFigs. 1, 2, 3).
We used these daily data both to analyze overall trends in seasonal climate in Acadia National Park between 1950 and 2021 and to derive independent variables for our statistical modeling of the timing of peak fall foliage, as described in the "Statistical Analyses" section below. We first calculated seasonal averages for each variable across Mount Desert Island, i.e., average winter (December, January, February), spring (March, April, May), summer (June, July, August), and fall (September, October, November) minimum temperature, maximum temperature, and total precipitation.
We also determined 95th and 99th percentile temperature and precipitation thresholds over our entire time period. For each smaller time period of interest, we calculated meteorological variables including "warm days"-defined as the number of days warmer than the 95th percentile maximum temperature threshold; "hot days"-defined as the number of days warmer than the 99th percentile maximum temperature threshold; "warm nights"-number of nights warmer than the 95th minimum temperature threshold; "hot nights"-number of nights warmer than the 99th minimum temperature threshold; "cool nights"-number of nights between the 54% and 73% percentile (~ 4.5-10 °C); "rain days"-number of days with more precipitation than the 50th percentile of daily precipitation; "pour days"-number of days with precipitation greater than the 95th percentile; "downpour days"-number of days with precipitation greater than the 99th percentile; and "droughts"number of periods with more than seven consecutive days of no precipitation (See Supplementary Table 1 for a list of all variables). We note that with our "cool/ warm/hot nights" variables, we did not subset the hourly data to nightly hours and thus assume that the daily minimum temperature represents the overnight low. Occasionally, the daily minimum temperature occurs during the day, like when a strong cold front moves in, but on average, daily minimum temperatures provide reasonable approximation of overnight temperatures.
All statistical tests were run using R v4.2.1. We ran Mann-Kendall (McLeod 2022) tests to determine if there were statistically significant trends in the aforementioned values over each time period of interest. The output of this test also provided Kendall's tau correlation coefficient, which determines the Deciduous forest land is highlighted in orange strength of the trend (values close to − 1 suggest a strong decreasing trend and values close to 1 suggest a strong increasing trend). We calculated the median rate of change of each significant trend using the Sen's slope estimator (Pohlert 2020).
We also used the coefficients from our final ERA5-Land based model to predict 2022 fall foliage timing, first using January 1, 2022-October 31, 2022 ERA5-Land data as the input independent variables, and then using Global Historical Climatology Network daily (GHCNd) data. Because there is no station on Mount Desert Island that collects both temperature and precipitation, we downloaded daily temperature data from the McFarland Station (USR0000MMCF) and downloaded and averaged daily precipitation data from the two precipitation stations in Southwest Harbor (Southwest Harbor 2.6 SE-US1MEHN0003, Southwest Harbor 0.9 NW-US1MEHN0064).

Remotely-sensed data and analysis
The timing of prolonged events, such as the duration and stages of fall leaf senescence and coloration, can be difficult to define, and is particularly tricky across Mount Desert Island where over half of the vegetation is classified as evergreen forest land (Fig. 3). A majority of institutions and states use percentage of leaves that have changed color to categorize various stages of fall foliage: the Maine Department of Agriculture, Conservation, and Forestry, which releases weekly reports in the fall, has defined five periods of fall foliage: very low-0-10% of leaves have changed color; low-10-30% of leaves have changed color, moderate-30-50%, high-50-70%, and peak-70-100%; weather.com has four defined categories of fall foliage: no change-0-10%; patchy-10-50%, near peak-50-75%, and peak, 75-100%; and the long-term ecological research (LTER) site at Harvard Forest records the percentage of leaf fall and leaf color every three to seven days. The United States Forest Service, however, uses five descriptive categorical values (0-4) to monitor fall foliage in Hubbard Brook, New Hampshire, e.g., "0: Most leaves fallen" and "2: Most leaves becoming yellow or red, and a few fallen leaves." And, when discussing with park visitors and local residents about what is 'peak' fall foliage, one can sum up their answers as "I know it when I see it." Thus, not only is monitoring fall foliage using remotely-sensed data difficult, but so too is validating the remotely sensed results in areas without long-term monitoring stations.
We adjusted the normalized brownness index (NBI) methodology developed by Zhang and Goldberg (2011). To calculate NBI, Moderate Resolution Imaging Spectroradiometer (MODIS) data is used to calculate annual normalized difference vegetation index (NDVI) timeseries. Then, a sigmoidal model is applied to the senescent phase of the smoothed NDVI timeseries (Zhang and Goldberg 2011). Then, using principles of spectral mixture modeling, fractions of 'green material' and 'brown material' are isolated, and the fraction of brown material parameter, i.e., the normalized brownness index, can then be used to determine fall foliage phases corresponding to percentage of leaf coloration and leaf fall (Zhang and Goldberg 2011). See Zhang and Goldberg (2011) for a complete theoretical background and methodology.
The two alterations we made to the algorithm were related to the input data and the annual time series data. Following Zhang and Goldberg (2011), we downloaded projected, mosaicked, and subsetted 500 m resolution MODIS MC43A4 Version 6 Nadir bidirectional reflectance distribution function (BDRF)-adjust reflectance (NBAR) data between February 5, 2000 (first data available) and December 31, 2021 over our study area using the USGS application for extracting and exploring analysis ready samples (AρρEEARS) tool. However, Zhang and Goldberg (2011) used a previous version of the dataset with a temporal resolution of 8 days, whereas we used Version 6 which has a daily temporal resolution. Also, instead of calculating annual NDVI time series to back out the NBI, we calculated, smoothed, and fit a sigmoid function to annual enhanced vegetation index (EVI) time series, because EVI better corrects for atmospheric conditions, reduces noise and is more sensitive to changes in vegetation cover (Huete et al. 1999).
We compare the results of the remotely sensed data to an archival historical analysis of peak fall foliage described in the "Historical Analysis" section below as another means of validating the archival work, but we note that we do not include these remotely sensed data in the statistical analyses described in "Statistical analyses" section. Moreover, calculating NBI works best when a pixel is predominantly deciduous forest. As highlighted by Fig. 3, although we ran our algorithm over each of the 2719 MODIS pixels that cover Mount Desert Island, only 104 pixels were comprised of more than 50% of deciduous or deciduous and mixed forest cover, only 24 pixels were comprised of more than 50% of deciduous forest cover, and only three pixels were comprised of more than 75% of deciduous forest cover. Thus, we focused our remotely sensed data analysis on the 104 pixels that contain at least 50% deciduous or deciduous and mixed forest cover (approximately 26 km 2 ).

Historical analysis
Because we both wanted to extend our time series of fall foliage to the pre-MODIS data era (pre-2001), and we required data to validate our remotely sensed analysis, we performed an archival analysis to determine the timing of peak fall foliage. Much effort was made to contact individuals, businesses, and scientific organizations in and around Mount Desert Island to collect historical information on fall foliage. Moreover, although the Maine Department of Agricultural, Conservation and Forestry (DCAF) has been keeping records of the timing of fall foliage across Maine for decades, only recently have they reported these weekly records online, and a majority of the paper records were thrown out (M. Macaluso, personal communication). Thus, we used National Park Service Acadia Monthly Narrative Reports, newspaper archives, fall foliage reports, and citizen science submitted photos to reconstruct a time-series of peak fall foliage dating back to 1951 (Table 1). Screenshots of a subset of the reports, newspaper clips, and photos can be found in the Supplementary Information (SI). We recognize that 'peak' fall foliage can be a multiday event, and thus note that that choosing one day as 'peak' for each year during which we have data is not ideal. However, we do believe we have derived the best representation of the timing of peak fall foliage over our study period.
Similar to the climate data analysis above, we conducted a Mann-Kendall test and calculated the Sen's slope estimator to determine if there were statistically significant trends in the timing of peak fall foliage over our time period.

Statistical analyses
We used multivariate linear regressions to determine the broad climatic drivers of the timing of peak fall foliage across Mount Desert Island.

Correlation analysis
We derived over 130 climatic variables between 1950 and 2021 using the same metrics as highlighted in the 'Climate Data and Analysis Section,' but broken down into even smaller time periods including, seasonal, monthly, and weekly time scales (Supplementary Table 2). Precipitation variables included our aforementioned total precipitation, number of 'pour days,' number of 'downpour days,' and number of 'drought periods' within each timescale of interest. We categorized temperature variables in three ways: (1) we calculated mean temperatures within each time period of interest (week, month, season); (2) we calculated the average minimum and maximum temperature within each time period of interest; and (3) we calculated the number of warm nights, hot nights, cool nights, warm days, and hot days within each time period of interest. Following Archetti et al. (2013), we then determined the correlation coefficient between each derived variable and peak fall foliage, and we identified those that were significant at a  1951, 1952, 1955, 1962, 1964, 1966NPS Acadia National Park monthly summary reports 1970, 1971, 1972, 1975-1977, 1980, 1985, 1989 Newspaper archives Maine fall foliage report, CSS photos p-value ≤ 0.10 to then use as independent variables in our regression modeling.

Regression modeling and model selection
We used a variety of variable selection techniques to derive models of the timing of peak fall foliage in Acadia National Park. We performed variable selection using best subsets (Miller 2020), best subsets with a validation set, k-fold cross-validation (Davison and Hinkley 1997; Canty and Ripley 2021), and lasso selection approaches. The best subsets technique identifies the set of variables that maximizes an objective indicator of model fit, including the residual sum of squares (RSS), adjusted R 2 , Bayesian Information Criterion (BIC) and Mallows C p statistic values (Beale et al. 1967;Hocking and Leslie 1967). For our purposes, we used adjusted R 2 ("Best Subsets A" below) and BIC ("Best Subsets B" below). The validation set approach splits the dataset into a training and validation set and computes the MSE of each best subsets selection model to determine which has the lowest MSE. The full dataset is then used to compute the coefficients in the final model. The k-fold cross validation takes the validation set approach further, splitting the data into k parts (here we used k = 10), or folds, and using k−1 parts to train the sample and the remaining partition to validate it. The prediction errors are then averaged, and the model with the lowest error is selected. Lastly, lasso (least absolute shrinkage and selection operation) performs variable selection by penalizing complex models while still minimizing the sum of squares to determine the best model (Tibshirani 1996). Variables from each temperature category described above were never combined in our variable selection process, thus, each variable selection process was run three times using each category of temperature variables. We ensured each model's selected variables did not exhibit multicollinearity by calculating variance inflation factor (VIF) values and removing selected variables if the VIF was greater than 5. Then, to choose the best-fit model of those derived from our variable selection processes highlighted above, we compared Akaike information criterion (AIC) and BIC values calculated using the full dataset, and the root-mean-squared-error (RMSE) and adjusted R 2 values, which were calculated using leave-one-out-cross validation.
We used our best-fit model to predict the timing of peak fall foliage in 2022 using both GHCN-d and ERA5-Land 2022 data as independent input datasets and with the future climate projection analysis described below.

Future climate scenarios
To determine how the timing of peak fall foliage might be expected to change in the future, we compared climate projection data from three Shared Socioeconomic Pathways: SSP1-2.6, SSP2-4.5, and SSP5-8.5. These SSPs consider socioeconomic factors, and how the world might address both climate change adaptation and mitigation. SSP1-2.6 is the "Sustainability-taking the green road" pathway, where global temperatures stay below an excess 2 °C of warming and there are low challenges to mitigation and adaptation; SSP2-4.5-is considered the "Middle of the road" pathway with medium challenges to mitigation and adaptation, and SSP5 is the "Fossil-fuel development" pathway with high challenges to climate adaptation and mitigation (O'Neill et al. 2017). We downloaded coupled model intercomparison project (CMIP6) (CESM2) community earth system model daily maximum near-surface air temperature, daily minimum near-surface air temperature, and total daily precipitation data from January 1, 2022-December 31, 2060 (O'Neill 2016; Danabasolgu et al. 2020).

Acadia's changing climate
Minimum temperatures, maximum temperatures, precipitation, and the number of warm nights, hot nights, warm days, hot days, and downpour days have all significantly (p ≤ 0.05) increased since 1950 across Mount Desert Island ( Table 2). The number of warm nights has increased by an average of 26 days/ year over the 71 year time period, and the number of warm days has increased by an average of 17 days/ year. Both hot nights and hot days have increased an average of 5 days/year since 1950. And the number of days exceeding the 95th percentile of precipitation ("pour days"), has increased by an average of 4 days per year between 1950 and 2021. Annual minimum temperatures have increased by an average of 1.2 °C, while annual maximum temperatures have increased by 1.12 °C.
There were no significant changes in spring temperatures or precipitation. Summer, fall and winter minimum temperatures significantly increased (Table 2). Summer and fall maximum temperatures significantly increased, winter maximum temperatures increased (Table2). Sen's slope tests indicated that minimum temperatures (i.e. night temperatures) have increased at a faster rate than maximum temperatures (i.e. day temperatures) ( Table 2). We include significant results in Table 2, but the results from all Mann-Kendall tests can be found in Supplementary Table3.

Peak fall foliage
Both the remotely sensed time series and archival time series highlight the variability in the timing of peak fall foliage. The remotely sensed time series shows no trend in the timing of peak fall foliage since 2001 (Fig. 4). Similarly, the archival time series shows no trend over this shorter, twenty-year time period (Fig. 5). Linear regressions of both the remotely sensed and archival time series between 2001 and 2021 have R 2 values of 0.0003 and 0.0005 respectively, and both linear regressions have a slope of − 0.02. However the correlation between the two datasets is 0.29 (SFigure 4).
Extending the time series to 1950, the archival time series demonstrates that the timing of peak fall foliage has been significantly delayed (p ≤ 0.05) since the start of time period (Fig. 5), with a Sen's Slope value of 0.21, indicating that, on average, peak fall foliage is now occurring almost two weeks later than it did in 1950 (Supplementary Table 3). When a linear regression was fit to the data, the slope was 0.13, which would indicate peak fall foliage is on average 9 days later over the 70 year time period. The archival data show that between 1950 and 1959, the average date of peak fall foliage was October 7, but between 2011 and 2021, on average, peak fall foliage occurred on October 15, happening as late as October 20 in 2021.

The effects of climate change on fall foliage
Of the 20 models that resulted from the variable selection process, September temperature-based variables were present in all models, with warmer temperatures being consistently correlated with delayed timing of peak fall foliage. May precipitation variables were present in 11 models, September precipitation variables were present in nine models, and both variables were positively correlated with delayed peak fall foliage. Variables representing rainfall during the first week of October were present in nine models, with increased early-October precipitation correlated with advanced peak fall foliage. Late summer (July and August) periods of drought were present in five models: in three models as linear terms correlated with delays in fall foliage, and twiceincluding in the best-fit model-as squared terms, indicating a complex relationship between late summer drought and timing of fall foliage. Temperatures during the first week in October were present in six models, but with conflicting correlation coefficients.
In Table 3, we show the results of the models using average temperature as they had slightly lower AIC, BIC, and RMSE values and higher R 2 than models using minimum, maximum, and  (Table 3). All model results are presented in the Supplemental Material (Supplementary Tables 5, 6, 7), and, with the exception of the temperatures-during-the-first-week-of-October variables, all other variables had comparable coefficients across models in which they were present. The average RMSE across all models was approximately 4 days.
When 2022 ERA5-Land data and GHNCd data were put into our best model (Table 3), which explained 47% of the variance and had a RMSE of 3.74 days, the model predicted that fall foliage would occur on day of year (DOY) 287 (October 14) and DOY 288 (October 15) respectively. As highlighted by the Maine DCAF reports and photos from ANP (see SI), peak fall foliage began between October 12 and 13, 2022 and a large storm with heavy rainfall caused the majority of trees across MDI to lose their leaves on October 14, 2022. Figure 6 shows the date of peak fall foliage between 2022 and 2060 using CESM2 data. During 2050-2060, SSP1's average date of peak fall foliage was October 30, while both SSP2 and SSP5's average date was only three days later, November 2. This small difference appears to be driven by differences in precipitation across the scenarios (Table 4): between 2022 and 2060, average September temperatures ranged from 21.6 °C (SSP1) to 21.7 °C (SSP5). However, SSP1 is wetter than SSP5 (Table 4).

Discussion
Climatic changes in Acadia National Park are more pronounced in some seasons than others. Agreeing with other New England based studies (Young and Young 2021), our results showed that all seasons except spring have seen significant increases in temperature and that nighttime temperature in these seasons have increased more rapidly than daytime temperature. Although winter temperatures have been increasing at almost double the rate of summer temperatures when averaged across the state of Maine (Runkel et al. 2022), this seasonal difference is not evident across MDI, likely due to its coastal geography.
Studies on the timing of various events related to leaf senescence have demonstrated that fall foliage is occurring both later and earlier, and it appears that directionality of these results are time-period and study-location dependent. Both our 20 year, remotely-sensed dataset and 70 year, archival dataset of the timing of peak fall foliage highlight that there exists much variability in the timing of leaf senescence and that the length of the time period of study is integral to understanding overall trends. Similar to our results, studies over longer (20 + year) time periods find fall foliage generally delayed (Ibanez et al. 2010;Gill et al. 2015)-in Acadia, the timing of peak fall foliage occurs almost two weeks later now than it did in 1950; while studies focusing on the MODISsatellite-data era (2001-present) find the timing of fall senescence events shifting earlier (Xi et al. 2015;Wu et al. 2018;Zhou et al. 2020)-in Acadia, our archival and remote-sensing results highlighted no discernable trend in the timing of fall foliage between 2001 and 2021. Interestingly, a study on the timing of peak fall foliage on the East Coast of the United States observed delays throughout the 1980s and 1990s but advances more recently (Liu et al. 2017), but the mechanisms behind this potential shift remain unclear.
Our regression and future scenarios modeling provides more insight into these time-series trends. As highlighted throughout the senescence phenology literature, changes in daylight and air temperature provide the primary cues to begin shutting down chlorophyll production. Our results suggest that September temperatures seem to be the most important variable in predicting of timing of peak fall foliage, with cooler temperatures resulting in earlier timing of peak fall foliage and warmer Septembers resulting in delays in the timing of peak fall foliage (Table 3, Supplementary Tables 2, 5, 6, 7).
As also highlighted by both the regression and SSP modeling results, the effects of precipitation on the timing of fall foliage, particularly over a large spatial scale, are difficult to disentangle. Increased May precipitation is correlated with delayed senescence and early October rainfall is correlated with earlier senescence. And, both September precipitation and July-August drought explained more variance when a squared term was used, highlighting the nonlinear impact of precipitation and drought on the timing of fall senescence events. The presence of May precipitation in over half of our final models is surprising, but not completely unfounded. Multiple studies have highlighted complex relationships between spring and fall phenology (Richardson et al. 2010;Misson et al. 2010;Fu et al. 2014;Keenan and Richardson 2015;Xie et al. 2015Xie et al. , 2018bLu et al. 2023). Spring drought stress has been correlated with earlier dormancy dates (Xie et al. 2015); excessive drought and heat in the spring has been correlated with earlier leaf coloration (Bigler and Vitasse 2021); increases in  spring productivity have been correlated with earlier senescence (Zani et al. 2020); and spring weather and climate variables-particularly spring precipitationwas found to be more strongly correlated to leaf coloration than to leaf drop (Xie et al. 2018b). Moreover, a Harvard Forest-focused study highlighted that although increased precipitation in May and June was generally correlated with earlier leaf coloration and leaf drop, the strength of this relationship was species specific: less drought tolerant species, particularly paper birch, were the most sensitive to spring precipitation; and species with higher drought tolerances, like oaks, beech and maples were the least sensitive to variation in spring precipitation (Xie et al. 2018b).
Our results then, both agree with other studies which have demonstrated that seasonal dynamics across the entirety of the growing season affect fall phenology, and emphasize that the specific mechanisms by which various stressors and events affect fall phenology need to be further explored.
Our study is not without limitations. Our remotelysensed fall foliage record and archival data record are only moderately correlated (R 2 = 0.29). Because no clear definition of 'peak fall foliage' exists, the archival record is not a perfect representation of the past. Moreover, leaf coloration is a prolonged event. Even the most current records-those published by the Maine DACF-are only published weekly on Wednesdays, making the direct comparison between weekly data and daily remotely-sensed data difficult. However, we believe our archival record is robust, particularly as both the archival and remotely sensed datasets demonstrate the same trends in the timing of fall foliage between 2001 and 2021. And ideally, our remotely sensed NBI record can serve as a more 'objective' baseline of the timing of peak fall foliage going forward.
We also rely on reanalysis data for our regression analyses because no consistent, long-term meteorological dataset of MDI exists. The R 2 s between daily and monthly minimum and maximum temperature data range between 0.90 and 0.99 (SFigs. 1 and 2). ERA5-Land minimum temperature data show a slight warm-bias (SFig. 1) and maximum temperature data show a slight cold bias (SFig. 2). Precipitation data are more weakly correlated than the temperature data but do not appear to be biased (SFig. 3). While temperature trends are generally well represented in climate models and reanalysis data, precipitation, extreme precipitation, and drought events are much more difficult to parameterize, model, and predict (Ukklola et al. 2020). Error in reanalysis data is an inherent limitation for this analysis. The inclusion of variables based on the identification of days based on percentile exceedances rather than absolute values (e.g., "very hot days") mitigates some-but not all-of this error due to the strong collinearity between the reanalysis and observational datasets. However, our finding that predicted peak fall foliage for 2022 differed by only one day between the reanalysis and observational data suggests that differences in these datasets have only a marginal impact on the substantive outcome.
Using both ERA5-Land and GHCNd data as inputs, our final model successfully predicted the timing of peak fall foliage in 2022. The Maine DACF October 12, 2022 report designated MDI as 'High: 50-70%' foliage coloration, and the October 19, 2022 report designated the region as "Past Peak", skipping "Peak: 70-100%" altogether, due to a large storm that swept across MDI on October 14. Our model predicted that peak fall foliage would occur on October 14 (ERA5-Land input data) and October 15 (GHNCd input data). Generally, our results show that ANP visitors who plan their trips months in advance would likely have a better chance of catching the fall colors in mid-October, rather than early October, which was the case seventy years ago.
Understanding how climate is affecting seasonal events like leaf senescence, and the timing of peak fall foliage is crucial for future park management in a National Park where fall tourism brings large economic gains to the local economy and fall visitor ship is increasing. Many local businesses in and around MDI begin to shut down for the season after Indigenous People's Weekend (second Monday of October), and in 2022, ANP's free park shuttle ended its service on Indigenous People's day, much to the chagrin of the influx of those mid-October visitors. Disentangling the relationship between climate change, the timing of fall events, and visitorship-the latter of which is the subject of future work-is essential. Case studies such as this can assist both visitors and park managers plan for a resilient, sustainable future across the National Parks.
Author contributions SAS designed the study and performed a majority of the data synthesis and analysis. KRS provided expertise and guidance on the statistical analyses and which climate datasets to use. MOH performed the archival fall foliage analysis. MGG collected, organized and synthesized GIS data related to the study. The first draft of the manuscript was written by SAS and all authors commented on subsequent versions and approved the final draft.
Funding This work was supported by a Second Century Stewardship Research Fellowship, public-private partnership led by Schoodic Institute at Acadia National Park and the National Park Service.
Data availability All open data used in this study has been referenced and hyperlinked appropriately and the historical fall foliage dataset will be published on the Harvard Dataverse upon manuscript acceptance.

Competing interests
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.