Using remote sensing to monitor the spring phenology of Acadia National Park across elevational gradients

. Greenup dates and their responses to elevation and temperature variations across the mountains of Acadia National Park are monitored using remote sensing data, including Landsat 8 surface re ﬂ ectances (at a 30-m spatial resolution) and VIIRS re ﬂ ectances adjusted to a nadir view (gridded at a 500-m spatial resolution), during the 2013 – 2016 growing seasons. The 30-m resolution provides a better scale for studying the phenology variation across elevational gradients than the 500-m resolution, as greenup dates monitored at 30-m scale have better agreement with leaf-out dates recorded in the ﬁ eld alongside the north – south-oriented hiking trails on three of the park ’ s tallest mountains (466 m, 418 m, and 380 m), and can provide landcover-speci ﬁ c analysis. The spring phenology responses to temperature and elevation vary among different spatial scales. Greenup dates of Acadia National Park monitored at 30-m scale show a weak advancing trend with higher spring temperature, while greenup dates monitored at 500 m show a weak delaying trend. The species mix within landcover at 30-m scale could weaken the advancing trend detected at ﬁ eld observation level. The landcover mix and elevation variation within 500-m scale could alter the spring phenology response to spring temperature variation. Greenup dates monitored at both 30-m and 500-m scales vary among different elevational zones, aspects, landcovers, and years. However, the relationship between greenup dates and elevation is rather weak.


INTRODUCTION
Spring vegetation phenology is the timing of physiological and morphological changes in vegetation that herald the start of new growth and are characterized by the timing of leaf out, budburst, and flowering. Environmental factors that influence spring phenology, such as temperature, sunlight, and precipitation, tend to vary widely in mountainous regions due to the influence of complex topography, leading to a variety of diverse niches within the broader ecosystem (Haslett 1997, Bales et al. 2006). Thus, the study of vegetation phenology in these topographically complex mountainous areas may provide early indications of changes in ecosystem processes resulting from climate change using field observation (Kimball et al. 2014, McDonough MacKenzie et al. 2019) and satellite monitoring (Fontana et al. 2008, Liu et al. 2016. Studies have documented Hopkins' bioclimatic law, where spring phenology shifts across elevation and latitude (Hudson Dunn and de Beurs 2011, Richardson et al. 2019). However, Vitasse et al. (2018) showed that spring phenology across elevations is becoming increasingly uniform. Even in a certain region, such as in Acadia National Park, though the field observation showed advancing spring phenology in warmer microclimates, there is no clear pattern in phenology progressing from lower elevation to higher elevation (McDonough MacKenzie et al. 2019).
Vegetation phenology has been historically monitored at the ground level through field surveys (Hameed and Gong 1994, Estrella et al. 2009, Vitasse et al. 2009). This approach is still widely used today to record phenological transitions and to study its responses to environmental factors, such as temperature, precipitation, and elevation, in various ecosystems (Dickinson et al. 2010, Polgar et al. 2014, Ma et al. 2016, McDonough MacKenzie et al. 2019. While ground-level phenological observations are straightforward to measure, manual monitoring can be quite laborious and time-consuming, especially in complex ecosystems or in mountainous regions where sampling is often constrained to small areas. Satellite imagery has been increasingly used to capture land surface phenology (LSP) over various ecosystems worldwide , Ganguly et al. 2010, Zhang et al. 2012. Land surface phenology describes the overall phenological dynamics of vegetation within a satellite sensor field of view (de Beurs and Henebry 2004). Generally, LSP is monitored using relatively coarse resolution imagery (250-1000 m), from sensors including the Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), the Visible Infrared Imaging Radiometer Suite (VIIRS), and Satellite Pour l'Observation de la Terre (SPOT)-VEGETATION, as these sensors can provide daily or near-daily observations worldwide (Fontana et al. 2008, Zhang et al. 2012, 2018a, Zhao et al. 2015. These sensors have been used in studies focused on phenology monitoring of mountainous areas in North America (Hudson Dunn and de Beurs 2011, Hwang et al. 2014), in Asia (Qiu et al. 2013, Liu et al. 2016, and in Europe (Fontana et al. 2008, Guyon et al. 2011. However, in mountainous environments, where elevation varies significantly within one pixel, these sensors may not provide the best spatial scale to describe the dynamics of vegetation phenology. Imagery from satellites such as Landsat, HuanJing, Sentinel 2, and GaoFen provides higher spatial resolutions and is more appropriate to measure vegetation phenology in heterogeneous regions (Walker et al. 2012, Melaas et al. 2013). However, due to the relatively low temporal frequency of these sensors, these higher resolution sensors often only acquire a few cloud-free images over an entire growing season. Efforts have been made to improve upon the temporal resolution of Landsat and other higher resolution imagers by using multiple-year aggregation methods (Avitabile et al. 2012, Melaas et al. 2013 and data fusion algorithms (Gao et al. 2006, Hilker et al. 2009).
The satellite-monitored phenology appears to vary at different spatial scales, especially in heterogeneous areas (Liu et al. 2017a, Peng et al. 2017. Also, single point measured field observation may not be representative for the whole satellite pixel monitoring , Liu et al. 2017a). This could possibly change the responses of monitored phenology to climate drivers among different spatial scales (Polgar et al. 2014). However, there is a lack of studies addressing this difference in scale with respect to phenology across elevational gradients. As mountainous region plays important role in v www.esajournals.org phenology studies, how the phenology responses to elevation and temperature vary across spatial scales needs to be further explored.
Acadia National Park was created in 1916 to conserve the natural environment of coastal Maine in the northeastern United States and minimize direct human impacts (Shafer 1999). Three field transects across elevational gradients have been established at Acadia, to specifically record the spring vegetation phenology and to study the phenology-temperature relationship from 2013 to 2016 (McDonough MacKenzie et al. 2019). These meticulously collected field data provide an ideal set of ground observations to assess the ability of satellite-derived products to capture the greenup dates of this topographically complex area of Acadia National Park.
In this study, we investigate variations in the spring phenology of vegetation along elevational gradients in Acadia National Park at both a finer spatial resolution by augmenting Landsat 8 data with Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) simulation, and a coarser resolution by using 500-m VIIRS daily Nadir Bidirectional Reflectance Distribution Function-Adjusted Reflectance. This study aims to: 1. Analyze the variation of the satellite-derived greenup dates among different elevational zones; 2. Study the spring phenology responses to elevation and temperature variations at different spatial scales.
In addition, the differences between satellitedetected greenup dates and field-observed leaf-out dates, and the spatial heterogeneities of elevation and landcover at 500-m scale are also examined.

Research area
Acadia National Park is made up of areas of Mount Desert Island, the Schoodic Peninsula, and other nearby offshore islands along the coast of Maine. The research here is focused on the Mount Desert Island unit of the park, which is mostly covered by evergreen forest, mixed forest, wetland, and shrub/scrub (Fig. 1). This topographically complex area contains the tallest coastal mountain on the North Atlantic seaboard, Cadillac Mountain, at an elevation of 466 m. This area experiences relatively long and cold winters, and abundant precipitations.

Satellite data
The satellite data used for LSP detection in this study include the 30-m atmospherically corrected surface reflectances from Landsat 8 Operational Land Imager (OLI) (Vermote et al. 2016), and the gridded 500-m daily Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (NBAR) from the NASA operational VIIRS products (Liu et al. 2017b).
The surface reflectance products from Landsat 8 OLI have a higher spatial resolution (30 m), but a lower temporal resolution (a 16-d repeat cycle). Fortunately, Acadia National Park falls in the overlap zone between Landsat path 11/row 29 and path 10/row 29; therefore, the number of Landsat observations is doubled. Landsat 8 OLI imagery is acquired with near nadir viewing geometry (only up to 7.5°off-nadir); thus, atmospherically corrected surface reflectances, downloaded from EarthExplorer, are used in this study directly, without any additional angular correction. The 16-d interval of Landsat hinders phenological extraction; thus, we used a fusion model to improve the temporal resolution (Gao et al. 2006). This fusion model pairs the daily coarse resolution MODIS NBAR imagery with the 16-d fine-resolution Landsat 8 OLI imagery, providing 30-m imagery interpolated at a daily timescale.
The 500-m VIIRS NBAR products are part of the VIIRS BRDF, Albedo, and NBAR product packages (VNP43I), which are available on a daily time step (Liu et al. 2017b). The VNP43I products, starting with the first launch of VIIRS in 2012, on board the Suomi-NPP satellite, have been developed to produce a continuous record equivalent to that from the MODIS record, the MCD43 products (Schaaf et al. 2002, Wang et al. 2018. The continuous NBAR products (VNP43I4) provide reflectance values with consistent nadir viewing geometry that remove the variation in the directional surface reflectances caused by the viewing geometry variation of VIIRS (AE56°across-track) (Liu et al. 2017b). In order to provide solar illumination geometry consistent with the Landsat acquisitions, the solar zenith angle of the Landsat overpass time is v www.esajournals.org used to calculate the VIIRS NBARs from the BRDF model parameter products (VNP43I1) for this study.

Greenup detection
Hybrid piecewise logistic model functions (Zhang et al. 2003, Zhang 2015 are used to detect the greenup from the satellite-derived Enhanced Vegetation Index (EVI) time series. Each parameter of the model can be assigned to a biophysical meaning related to vegetation growth or senescence. This method has been used to derive phenology from various satellite sensors, such as MODIS (Zhang et al. 2003), VIIRS (Zhang et al. 2018b), and Landsat (Zhang et al. 2020).
Annual EVI time series of vegetation pixels are composed of increase and decrease sections, representing growth and senescence cycles, respectively. Pixels that do not have sustained EVI increase and decrease are considered as Missing values due to clouds in the time series are filled by linear interpolation of two nearest dates. Also, Savitzky-Golay filters are applied to further remove irregular variations. Then, the EVI time series are separated into single growth or senescence cycle, and each cycle is modeled using the following logistic function: where t is time in days, EVI(t) is the EVI value of day t, a and b are fitting parameters, c + d is the maximum EVI value, and d is the initial background EVI value. The two points with maximal curvatures in each cycle represent the days when phenology changes from one phase to another. The greenup, one of the parameters of the model, which is defined as the onset of the EVI increases, is the day of the first maximal curvature of growth cycle.
Since the 500-m VIIRS products do not include blue band spectral information, the EVI2 is substituted for EVI (Jiang et al. 2008) to derive greenup dates at 500 m. The Landsat 8 EVI is calculated using Eq. 2, and the VIIRS EVI2 is calculated using Eq. 3.
where ρ R , ρ NIR , and ρ B represent the reflectance of the red, near-infrared (NIR), and blue bands, respectively.
With the establishment of satellite-derived EVI time series, piecewise logistic functions are fitted to the various time series to compute the greenup dates of Acadia National Park.

Field-observed phenology along the hiking trails
Three north-south field transects have been established along hiking trails on Cadillac, Pemetic, and Sargent mountains in Acadia to record the leaf-out and flowering phenology from 2013 to 2016 (Fig. 2). This "trail-astransects" approach was designed as a transition from a monitoring scheme run by a researcher to a citizen science program; therefore, it represents an easy approach for hikers to collect data while they climb the mountain trails (McDonough MacKenzie et al. 2019). This approach is similar to the method of the USA National Phenology Network (Denny et al. 2014).
Each hiking trail is separated into different elevational zones based on aspect (north and south) and elevation (low, medium, high, and summit). The altitude range of each elevational zone is standardized across these Acadia mountains and aspects: "low" is below 183 m, "medium" represents 183-274 m, "high" represents 274-366 m, and "summit" is above 366 m. The hiking trails on Cadillac, Pemetic, and Sargent mountains have sections in each elevational zone. A fourth hiking trail, Giant Slide, connects with Sargent Ridge (Fig. 2). Therefore, across all trails, there are three distinct transects in the low, medium, high, and summit elevational zones on both south and north aspects, respectively. In total, there are 24 elevational zones.
In each elevational zone, phenophases of understory species and of overstory deciduous trees along the hiking trail are recorded twice a week from 1 April through 30 June. Note that the evergreen trees in these regions are not monitored. The leaf-out date of each species is defined as the first observed leaf-out date of any plant for that species in each elevational zone along the hiking trails. In addition, within each monitoring zone, species are categorized as "common," "occasional," and "uncommon" based on the abundance occurring in that monitoring zone (Appendix S1: Tables S3, S4).
HOBO pendant temperature loggers were deployed in all elevational zones along the hiking trails to record temperature from July 2013 through July 2016. Mean spring temperature (mean temperature of March and April) is used to describe microclimate variations along the elevational gradients and to model phenological response to temperature variations in Acadia National Park (McDonough MacKenzie et al. 2019).

Analysis of greenup dates
The greenup dates as detected by satellites are analyzed through the following steps: 1. Satellite-estimated greenup dates are compared with field-monitored leaf-out dates by v www.esajournals.org considering each elevational zone on each ridge and aspect as a separate community. Field-monitored community leaf-out date of each elevational zone is derived by averaging field-monitored leaf-out dates of "common" species of that zone. Satellitemonitored community greenup dates of the whole elevational zone and of each landcover in the zone are derived by averaging the satellite-estimated greenup dates of all vegetation pixels and of each vegetation landcover pixels, respectively; 2. The variation of satellite-monitored greenup dates with elevation is analyzed using ANOVA and boxplots at both 30-m and 500m scales. Also, mixed linear regression models (Harrison et al. 2018) between elevation and greenup dates by taking year, aspect, and landcover variations into consideration are performed at both 30-m and 500-m v www.esajournals.org 6 December 2021 v Volume 12(12) v Article e03888 scales to study phenological response of the greenup dates to elevation at different scales. We used a fit linear mixed-effects model (fitlme) function in MATLAB on greenup dates with uncorrelated random effects for elevation and intercept grouped by landcover, year, and aspect in the following form.
greenup ∼ elevation þ ð1jlandcoverÞ þð1jyearÞ þ ð1jaspectÞ þðelevation À 1jlandcoverÞ þðelevation À 1jyearÞ þðelevation À 1jaspectÞ (4) 1. The phenological response of the greenup dates to the temperature variation is analyzed using satellite-estimated community greenup dates and the mean spring temperature recorded by the HOBO temperature loggers in the field.

Satellite-detected greenup dates of Acadia National Park
The greenup dates as detected by the VIIRS EVI2 and by the Landsat 8 EVI augmented with the simulated STARFM (Landsat 8 + STARFM) EVI are displayed in Fig. 3. At 30-m scale, obviously more spatial details are captured than at the 500-m scale in this topographically complicated area, as the 500-m pixel edges can be seen from these figures.
In addition, the greenup dates detected using these two spatial scales are not always close, as VIIRS-detected greenup dates are earlier than Landsat 8 + STARFM-detected dates of years 2015 and 2016.
The spatial heterogeneities of landcover and elevation within 500-m pixels of Acadia National Park are displayed in Appendix S1: Fig. S1, which show that less than 10% 500-m pixels in this area contain only one landcover, more than v www.esajournals.org 90% 500-m pixels contain at least two landcovers, indicating a majority of 500-m pixels are mixed with different landcovers. The difference between the highest elevation and the lowest elevation in a 500-m pixel could reach to more than 200 m (Appendix S1: Fig. S1b). Only around 10% 500-m pixels have elevational difference smaller than 15 m. Therefore, Acadia National Park is landcover-mixed and topographically heterogeneous at 500-m scale.
Comparison of satellite-detected greenup dates with field-recorded leaf-out dates The differences between satellite-monitored greenup dates and field-recorded leaf-out dates over the hiking trails are displayed in Fig. 4, where Fig. 4a shows the comparison for all vegetation pixels at 500-m scale, and Fig. 4b shows the comparison for all vegetation pixels at 30-m scale. In addition, as the field-monitored common species are herbaceous and shrubs, the 30-m greenup dates of pixels, whose landcovers are herbaceous and shrubs, are compared with corresponding field-monitored community leaf-out dates (Fig. 4c). This landcover-specific comparison cannot be performed at 500-m scale as most 500-m pixels are landcovermixed (Appendix S1: Fig. S1) in Acadia Nation Park.
Over the hiking trails of these three mountains, 30-m satellite-derived community greenup dates of all vegetation pixels have lower bias (−3.34 d) with the field-observed community leaf-out dates than 500-m satellite-derived community greenup dates of all vegetation pixels (−6.37 d). The agreement between satellite monitoring and field observation improved when only herbaceous and shrub pixels at 30-m resolution are used for comparison as displayed in Fig. 4c (bias = 1.76 d, RMSE = 8.27 d).
Therefore, the greenup dates derived at 30-m scale agree better with the field-monitored leafout dates than the greenup dates derived at 500m scale do. Also, taking landcover and species function type into consideration improves the agreement between satellite-derived and fieldmonitored spring phenology. These indicate that a finer scale is more suitable for phenology monitoring in this area, and landcover mix could drive the disparity between field-measured and remotely sensed phenology.

Phenology response to spring temperature over the hiking trails
The linear regressions between satellitedetected greenup dates and field-recorded temperatures have different trends at 500-m and 30-m scales (Fig. 5), where the 30-m satellite-derived greenup dates have slightly advancing trend with increasing temperature, and the 500-m satellite-derived greenup dates have slightly delaying trend with increasing temperature. Note that neither model is significant (Fig. 5).
Landcover-specific analyses are performed at 30-m scale (Fig. 6). The greenup dates of most landcovers display advancing trend with higher temperature as well, except for mixed forest (Fig. 6), which is similar to the trend at 500-m scale. However, the correlation is not statistically significant for all landcovers. Different landcovers display various advancing trends, as the greenup dates of deciduous forest advance 1.63 d, and the greenup dates of evergreen forest advance 0.44 d with temperature increase of 1°C. The greenup dates of shrub/scrub and herbaceous pixels have a similar trend that advance around 0.8 d with 1°C increase of mean spring temperature.

Phenology response to elevational gradients in Acadia National Park
The ANOVA results (Table 1) indicate that the greenup dates of the whole park varied significantly among different elevational zones at both 30-m resolution (F = 517.293, P < 0.001) and 500m resolution (F = 240.30, P < 0.001). Also, the greenup dates of the whole park varied significantly among different aspects and years. The ANOVA results for each landcovers at 30-m resolution show significant variations in greenup dates among the different elevational zones except for wetland (F = 0.703, P = 0.550) as displayed in Table 1.
The boxplot further illustrates the satellite greenup date variation along elevational gradients. The satellite greenup dates monitored at different spatial scales display different responses to the elevation in different years (Fig. 7). At 30-m scale, greenup dates and elevations are well correlated in 2013, as the low elevational zone had the earliest greenup dates and the summit zone had the latest. The median value of the greenup dates of each elevational zone (the red central mark of the boxplot) advances 1.3 d on v www.esajournals.org average as one moves to lower elevations (Fig. 7). However, in 2014-2016, the relationship between greenup dates and elevational zones is less clear. In these years, the lower elevational zones do not consistently show the earlier greenup dates, and the phenology appears to be decoupled from the elevational gradient (Fig. 7). At 500-m scale, greenup dates and elevations are well correlated in 2014, as the low elevational zone had the earliest greenup dates and the summit zone had the latest (Fig. 7). However, for other three years, the relationship between greenup dates and elevational zones is less clear. Though the summit elevational zone has latest greenup dates in 2013 and 2015, the low elevational zone does not have the earliest greenup dates (Fig. 7).
The landcover-specific boxplots of each landcover at 30-m resolution show that the greenup dates of each landcover have different responses to elevation (Fig. 8). In mixed forests, the greenup dates are correlated with elevational zones (earlier greenup at lower elevations) for all four years. Similar trends are also found in evergreen forest in 2014, 2015, and 2016, and in deciduous forest in 2013 and 2014 (Fig. 8). Herbaceous pixels generally reflects later greenup dates at higher elevations. However, herbaceous pixels in the summit zone deviates from this pattern; here, the median greenup date occurs 2 d before the high elevational zone.
The mixed linear regressions between elevation and greenup dates at both 30-m and 500-m scales are displayed in Table 2, where the intercept indicates the greenup dates at elevation of 0, and the coefficient indicates the rate of greenup dates that advance with higher elevation (day per meter; Fig. 9). It shows that there is a weak advancing trend of greenup dates when elevation increases at both 30-m (coefficient = −0.0064, P = 0.44) and 500-m (coefficient = −0.0049, P = 0.49) ( Fig. 4. Continued) v www.esajournals.org scales. However, the intercepts and coefficients are different at two scales, the intercept at 500-m scale is 7 d earlier than at 30-m scale. At 30-m scale, the yearly variation has the biggest random effect to the intercept, and aspect has the biggest random effect to the coefficient, and so does at 500-m scale. Landcover also affects the relationship between elevation and greenup dates, especially to the intercept. This indicates in this study area that there is no significant correlation between elevation and greenup dates, due to the high year-to-year variation and diversity of landcover and aspect.

Satellite-monitored greenup dates detected at different spatial scales
The differences in the satellite-detected greenup dates between 30-m and 500-m observations appear to be mainly caused by spatial differences, as more spatial details can be captured by the 30-m data than by the 500-m data. Also, the greenup dates detected at 30-m resolution appear to agree better with the ground-measured leaf-out dates than at 500-m resolution. A majority of 500-m pixels of Acadia National Park are mixed with various landcovers and also are topographically complicated. Given this spatial variation in the spring phenology, 30 m is a more appropriate scale for phenology monitoring in this topographically complicated Acadia National Park. Other studies also supports that in heterogeneous area, 30 m can provide more spatial details in phenology monitoring (Walker et al. 2012, Liu et al. 2017a. However, it must be recognized that the STARFM results, which are used for phenology detection at 30-m resolution, are not values truly observed by satellite sensors, and this process can introduce its own uncertainties. Other approaches to improve the temporal resolution may be used through combining images from a number of sensors that have similar higher spatial resolution designs, such as Landsat, HuanJing, and Sentinel 2 (Wu et al. 2015). However, the differences between all of these sensors may bring additional uncertainties to the time series needed for phenology monitoring.

Satellite-monitored greenup dates responses to elevation and temperature at different scales
Elevation and spring temperature have been shown to affect the greenup dates at both 30-m and 500-m scales in Acadia National Park. However, the greenup date responses to elevation and spring temperature variations are more complicated in Acadia National Park than the findings in other regions where the spring phenology generally advances with higher spring temperature and at lower elevational zones (Guyon et al. 2011, Hudson Dunn and de Beurs 2011, Qiu et al. 2013. The correlation between satellite-derived v www.esajournals.org greenup dates and spring temperature is not statistically significant. In addition, the satellitedetected greenup dates showed different responses to spring temperature at different scales in Acadia National Park, with an advancing trend with higher spring temperature at 30-m scale and a delaying trend at 500-m scale. Elevational zones with higher altitude do not always correspond to later greenup dates, and the elevational trend in greenup dates varies across different years at both 30-m and 500-m scales. Also at 30-m resolution, the elevational trend of spring phenology also varies with landcover. In the field, a significant linear relationship between leaf out and spring temperature was v www.esajournals.org found at the species level in this region (McDonough MacKenzie et al. 2019). Possible reasons for the weak relationship between greenup dates and spring temperature when using satellite imagery for phenology monitoring are the variation in phenological responses to temperature among difference species and the heterogeneous microclimates, which is represented by temperature here, within each elevational zone. McDonough MacKenzie et al. (2019) found that magnitude of the relationship between spring phenology and temperature varies across species in this region. Diez et al. (2012) also illustrated significant variation among species in sensitivities to climate. Temperature recorded by one HOBO pendant temperature logger in the field is used to represent the temperature of each elevational zone; therefore, the temperature variation in each  elevational zone is not considered in the analysis here. A pixel-by-pixel analysis at 30-m or even finer spatial scales may help us to understand the greenup date response to spring temperature instead of using the elevational zone average. However, the thermal sensors of satellites generally have lower spatial resolution; it is difficult to directly derive 30-m or finer spatial scale temperature using satellite imagery (Roy et al. 2014).
In contrast to this study, earlier researchers found significant effects of elevation on spring phenology in mountainous region at continent scale (Vitasse et al. 2018), at national scale (Richardson et al. 2019), and in region with larger elevation variation (Norman et al. 2017). A roughly 20-d advance trend per 1000 m was found (Vitasse et al. 2018, Richardson et al. 2019. While the elevation difference in Acadia National Park is less than 500 m, the spring phenology advance trend with elevation was expected to be less. Also, the park's heterogeneous species composition across elevational zones could complicate elevational greenup trends. In addition, elevation and local temperature are often decoupled in the mountains of Acadia National Park; higher elevational zones do not always correspond to colder spring temperature (McDonough MacKenzie et al. 2019). These factors could weaken the spring phenology response to elevation that detected by satellites here. Therefore, to study elevational effects on phenology in mountainous region with low elevation range, such as the Acadia National Park, imagery with high temporal resolution is just as important as high spatial resolution. Notes: L represents random-effects covariance parameters of landcover, Y represents random-effects covariance parameters of year, and A represents random-effects covariance parameters of aspect; ". . ." indicates the parameter is not considered in the regression.
The regression equation is: greenup dates = intercept + elevation × coefficients. Factors causing the differences between the fieldobserved leaf-out dates and satellite-detected greenup dates Comparison of the 30-m greenup dates of herbaceous and shrub pixels with the fieldobserved leaf-out dates with corresponding species function type illustrates a good agreement. However, the agreement of field observations and satellite detections at both 30-m and 500-m scales over all vegetation pixels is less strong.
The field observations focus on collecting data along hiking trails, which represent a smaller area than that captured by a 30-m Landsat pixel, and a way much smaller area than that captured by a 500-m VIIRS pixel. The species composition varies significantly within the park; thus, the species composition along the hiking trails may be quite different than the species composition of an entire Landsat pixel that covers the trails. In addition, landcover type affects the greenup dates in the park. Therefore, the sampling bias between field observation and satellite, and the landcover mix along the hiking trails could contribute to the differences between the fieldobserved leaf-out dates and the satellite-detected greenup dates.
PhenoCam has been established in Acadia National Park, which could also be an approach to validate the satellite-detected phenology (https://phenocam.sr.unh.edu/webcam/sites/acadia/). However, overlaying the PhenoCam field of view with satellite pixels may be challenging in such a topographically complicated area.

CONCLUSION
This study documented elevational variation of greenup dates in Acadia National Park and its responses to temperature across different spatial scales, and while leveraging field-observed phenology to validate remote sensing monitoring over this topographically complex area.
We found variation in greenup dates across the elevational zones, aspects, and landcovers in Acadia, although the relationship between greenup dates and elevation changes among years, landcovers, and spatial scales. The spring phenology appears to be decoupled from the elevational gradient; there is no clear trend as monitored by satellite at both scales. Meanwhile, the responses of spring phenology to spring temperature showed different trends at different spatial scales. The greenup dates advance with higher spring temperature at 30-m scale, but delay with higher spring temperature at 500-m scale. Therefore, the relationships of spring phenology to temperature and elevation do not hold across spatial scales. The heterogeneous species composition and their varied responses to elevation and temperature affect the response of satellitemonitored greenup dates to elevation and temperature.
Acadia National Park is both landcover-mixed and topographically heterogeneous at 500-m scale. The 30-m scale is more suitable for phenology monitoring in this area, as more phenology details can be detected at this scale than at the more commonly used coarse resolutions of satellite data. Also, the greenup dates derived from satellite at 30-m scale agreed well with the fieldmonitored first leaf-out dates and can perform landcover-specific analysis. The low biases and RMSEs between field-observed and satellitemonitored phenology assure that, although it is difficult to monitor individual species phenology using satellite, the overall phenological dynamics detected by satellite can be used to fill the gap and expand the geographic scope of field-based research in topographically complex areas that are hard to cover on the ground.
For this study, data fusion (of Landsat and MODIS) is utilized to improve the temporal resolution of the time series required for phenology monitoring. In the future, higher resolution images from other platforms, such as Sentinel 2 A/B, HuanJing, and GaoFen, may also be incorporated to produce higher spatial and temporal resolution time series for phenological monitoring. With the availability of these higher resolution satellite data, monitoring and analyzing phenology using satellite imagery in highly complex topography will likely become more viable and informative.

ACKNOWLEDGMENTS
The processing of MODIS and VIIRS has been supported by NASA grants NNX14AI73G and NNX14AQ18A, and the Landsat processing has been supported by USGS Award G12PC00072. The first author is also supported by NSFC 41901367 during the