Drivers of contemporary and future changes in Arctic seasonal transition dates for a tundra site in coastal Greenland

Climate change has had a significant impact on the seasonal transition dates of Arctic tundra ecosystems, causing diverse variations between distinct land surface classes. However, the combined effect of multiple controls as well as their individual effects on these dates remains unclear at various scales and across diverse land surface classes. Here we quantified spatiotemporal variations of three seasonal transition dates (start of spring, maximum normalized difference vegetation index (NDVImax) day, end of fall) for five dominating land surface classes in the ice‐free Greenland. Using a distributed snow model, structural equation modeling, and a random forest model, based on ground observations and remote sensing data, we assessed the indirect and direct effects of climate, snow, and terrain on seasonal transition dates. We then presented new projections of likely changes in seasonal transition dates under six future climate scenarios. The coupled climate, snow cover, and terrain conditions explained up to 61% of seasonal transition dates across different land surface classes. Snow ending day played a crucial role in the start of spring and timing of NDVImax. A warmer June and a decline in wind could advance the NDVImax day. Increased precipitation and temperature during July–August are the most important for delaying the end of fall. We projected that a 1–4.5°C increase in temperature and a 5%–20% increase in precipitation would lengthen the spring‐to‐fall period for all five land surface classes by 2050, thus the current order of spring‐to‐fall lengths for the five land surface classes could undergo notable changes. Tall shrubs and fens would have a longer spring‐to‐fall period under the warmest and wettest scenario, suggesting a competitive advantage for these vegetation communities. This study's results illustrate controls on seasonal transition dates and portend potential changes in vegetation composition in the Arctic under climate change.


| INTRODUC TI ON
Arctic vegetation has undergone intricate changes since the 1980s as a consequence of Arctic amplification of climate warming (Post et al., 2019;You et al., 2021) and associated increased precipitation (Bintanja, 2018), varying between regions (Loranty et al., 2016), latitudes (Miles & Esau, 2016), and vegetation types (Frost & Epstein, 2014).On average, vegetation greenness has increased (Jia et al., 2009;Pearson et al., 2013) in the Arctic during the past four decades, but decreased in the boreal zone and southern West Siberia from 1980 to 2014 (Ju & Masek, 2016;Miles & Esau, 2016;Verbyla, 2008).Variations in vegetation greenness are associated with changes in seasonal transition patterns (Gao et al., 2023), which can have a strong correlation with gross primary productivity (Piao et al., 2019) and affect the total amount of carbon absorbed by vegetation (La Puma et al., 2007;Street et al., 2013), thus comprehending the effects of climate change on the seasonal transition shifts is important for assessing the Arctic ecosystem carbon budget.
To reflect vegetated-land surface phenology, seasonal transition dates are widely extracted using the normalized difference vegetation index (NDVI) from satellite data (Fawcett et al., 2021;Wu et al., 2017).Here we adopt the terms "start of spring (SOS)" and "end of fall (EOF)" because we focused on land surface phenology, which is related but not identical to plant phenology (White et al., 2009).Moreover, the SOS and EOF are strongly linked to the timing of vegetation leaf-out and leaf senescence, respectively (Klosterman et al., 2014), indicating the changes in growing season length can be strongly portended by the SOS and EOF.In the Arctic, longer growing seasons were observed due to earlier onset and delayed ending during 1982-1999 but shorter growing seasons were found due to delayed onset and earlier ending during 2000-2014 (Park et al., 2016).26.8% of the Arctic region showed an earlier timing of peak photosynthetic activity, with a rate of −1.09 ± 0.29 days/ decade during 2000-2016 (Park et al., 2019).The shifts in seasonal transition dates are linked to interactions between changes in snow cover (Wu et al., 2023;Xie et al., 2018), air temperature (Berner et al., 2020;Krab et al., 2018;Reichle et al., 2018), precipitation patterns (Huang et al., 2022;Shen et al., 2011) and some other local factors, suggesting that the effect of individual factors as well as the combined effect of multi-factors on seasonal transition dates warrants further investigation.
Previous studies have documented that the overall advancing trend of SOS was attributed to increasing temperatures since the 1980s (Piao et al., 2019;Prevéy et al., 2017) but spatial heterogeneous, with positive, negative, or no changes observed (Zeng et al., 2011).Different trends in the SOS can be partly attributed to snow cover (Johansson et al., 2013;Pedersen et al., 2018).Snow cover affects the SOS by warming up soil temperature during winter (Ge & Gong, 2010), and altering soil moisture (Potopová et al., 2016) and nutrient cycles (Brooks et al., 2011) during spring.Satellite and experimental data have verified that earlier snow ending day (SED) primarily accounts for an earlier growing season in the Arctic (Kelsey et al., 2021;Slatyer et al., 2022;Zheng et al., 2022), with a 0.56-day advancement in the growing season per 1-day advancement in SED across the Pan-Arctic region during 1982-2015 (Wu et al., 2023).

Furthermore, snow metrics play a crucial role in vegetation vigor
during the growing season (Buus-Hinkler et al., 2006;Pedersen et al., 2018;Westergaard-Nielsen et al., 2017), subsequently impacting the maximum NDVI (NDVI max ) day (Kelsey et al., 2021).The EOF is driven by multiple factors, including photoperiod, temperature, and nutrients (Gill et al., 2015), and later snow onset day (SOD) can postpone the EOF (Cooper, 2014).In addition, the decline in wind is found to significantly extend the date of autumn foliar senescence (Wu et al., 2021).The effects of climate and snow cover on seasonal transition dates vary between different vegetation communities (Krab et al., 2018;Wipf & Rixen, 2010) and terrain (An et al., 2018;Tomaszewska et al., 2020).In-situ studies have observed that shrub phenology is more responsive to snowmelt time than graminoids, which are rather driven by temperature (Legault & Cusa, 2015;Wipf, 2010).The SOS occurs earlier in shrub-rich tundra compared to wet tundra (Jia et al., 2009).Long-term satellite observations show that the SOS exhibited a more rapid advancement in higher-elevation mountain areas in Canada's Arctic from 1985 to 2013 (Chen et al., 2021), but summer NDVI was negatively related to latitudinal gradients (Wu et al., 2020).Terrain affects vegetation phenology indirectly in winter, with slopes being more influential than aspects in establishing links between snow and vegetation phenology in highland regions (Tomaszewska et al., 2020).The effects of diverse controls on different seasonal transition dates are thus complex and hard to quantify across different land surface classes.
Significant temperature and precipitation increases are projected for the Arctic, with air temperature expected to rise by 3.5-10.4°C(Cai et al., 2021) and precipitation to increase by over 50% (Bintanja & Andry, 2017) by 2100relative to the start of the century.
More intense warming and wetting are expected in autumn and winter, accompanied by more winter rainfall (Bintanja, 2018;Overland et al., 2019).Correspondingly, snow mass is projected to increase by up to 15% over some parts of the Arctic by 2050 relative to the 20th century (Callaghan et al., 2011) although snow cover extent is projected to shrink (Mudryk et al., 2020).In addition, it is reported more than 50% of the Arctic area will shift to a different vegetation class with the expansion of woody shrubs and trees under various scenarios by the 2050s (Pearson et al., 2013).However, the response of snow cover and vegetation to climate change is varied and complex between different scales, combined with projected vegetation redistribution and its impacts on ecosystem feedback, intensifying concerns about the implications of climate change on future land surface seasonal transition dates.
Several studies have investigated relationships between vegetation growing seasons and climatic factors from the plot (Elmendorf et al., 2012;Niittynen et al., 2020) to large scale (Reichle et al., 2018) in the Arctic.However, seasonal transition dates and their controls vary with different scales (Niittynen et al., 2020); the scarcity of ground-based or high-resolution satellite snow data in the Arctic has limited our understanding of the impact of snow cover on the seasonal transition dates at a small or medium scale.In addition,

| Study site
The study region is located on the south part of Disko Island, Western Greenland (69°16′ N, 53°27′ W), covering an area of approx.21 km 2 (Figure 1).The area is characterized as a transition region between the low and high Arctic with an annual mean air temperature of −3°C and precipitation of 418 mm (of which approx.one-third falls as snow) during 1991-2017 (Zhang et al., 2019).Most Greenlandic land surface classes are found in this region (Karami et al., 2018), including abrasion (2.05 km 2 , 10%), fen (0.70km 2 , 3%), dry heath (6.43 km 2 , 31%), wet heath (5.05 km 2 , 24%), and tall shrub (0.74 km 2 , 4%), and well-drained mesic tundra heath is the most dominating ecosystem type (Figure 1b).In our study area, abrasion is characterized as dry with very low vegetation activity during the growing seasons and includes very sparse Dryas or graminoids; fen is covered with grasses and mosses with mostly saturated soil; dry heath is mainly composed of Betula and Vaccinum; wet heath is a mix of Betula, Vaccinum, Salix, and Empetrum but the height of the plant doesn't exceed 40cm; tall shrub includes Salix and Betula with taller than 40cm height (Karami et al., 2018).Ten meteorological stations from GEM (https:// data. g-e-m.dk/ ) and DMI (https:// www.dmi.dk/ ) are distributed in this region and provide hourly air temperature, precipitation, wind speed and direction, and relative humidity measurements from August 2010 to December 2020.We used the ArcticDEM (https:// www.pgc.umn.edu/ data/ arcti cdem/ ) data with 32 m spatial resolution to acquire the elevation, slope, and aspect in the study area.Three snow fences are located in the mesic tundra heath.These are photographed once per day from October through March and three times per day during the rest of the year.The snow depth (Liu et al., 2023) and vegetation coverage reflected from the photos can serve as validation data for simulated snow cover (see Figure S1) and extracted seasonal transition dates (see Figure 2).

| Precipitation correction
Precipitation gauge records have long been realized as underestimations of true precipitation amounts due to under-catch in windy conditions, and solid precipitation is more affected by wind speed than rainfall (Legates, 1995;Wolff et al., 2015;Yang et al., 2001) methods for the Geonor automated weighing gauge.The Pluvio gauge has the same physical characteristics relevant for precipitation catchment as the Geonor and Wolff's correction method has been used to adjust precipitation amounts from Pluvio in Canada (Milewska et al., 2019).Therefore, we corrected the precipitation based on Wolff's method in our study area: where p T is the corrected precipitation, p M is measured precipitation, T is the hourly air temperature, and v is the wind speed measured at gauge height.

| SnowModel
SnowModel is a distributed snow-evolution model (Liston & Elder, 2006a) and is composed of four submodels.MicroMet is a data assimilation and interpolation model to produce high-resolution meteorological forcing data spatially (Liston & Elder, 2006b).EnBal is responsible for surface energy calculations based on near-surface atmospheric conditions provided by Micromet (Liston et al., 1999).
SnowPack is used to define snowpack changes with changing meteorological conditions (Liston & Hall, 1995).SnowTran-3D simulates the snow depth by considering wind-induced redistribution (Liston & Sturm, 1998).The model inputs are corrected precipitation, wind speed and direction, air temperature, humidity, topography (DEM data), and vegetation type (Karami et al., 2018).All input meteorological data are hourly and from ground observations.The spatial resolution of outputs from SnowModel is determined by the input DEM.We took one-year data as an example to simulate the snow cover distribution by inputting 10 and 32 m spatial resolution DEM separately and found no difference in simulation accuracy between the two spatial resolutions.Thus, we used the DEM data with 32 m spatial resolution to simulate snow for improving computing efficiency.We made full use of all input data and simulated snow cover evolution from August 2010 to December 2020.In addition to simulated daily snow depths from SnowModel, we used daily precipitation, wind speed, and air temperature interpolated in space by MicroMet based on ground observations.We found the intraand inter-annual variations of simulated snow depth were consistent with the observed values (R = 0.94, RMSE = 0.12 m, Figure S1).
(1) | 5 of 17 There is no difference between the simulated and observed snow onset day, and the mean simulated snow ending day was delayed 1.6 ± 1.5 days compared to observations during 2012-2020.

| Seasonal transition dates
To quantify the seasonal transition dates, we used NDVI derived from Sentinel-2 MultiSpectral Instrument (Level-1C) images during 2016-2020 based on Google Earth Engine (https:// devel opers. google.com/ earth -engine/ datas ets/ catal og/ COPER NICUS_ S2).We performed an atmospheric correction (Yin et al., 2022) on the images before calculating NDVI.The months from May to October were set as the study period each year.The quality control process includes three steps: (i) the cloud was masked according to the QA60 band; (ii) images were removed if the number of pixels with NDVI values outside the range of −1 to 1 exceeds 30% of the total pixels while extracting the median value of each date; (iii) NDVI outliers resulting from cloud mask errors (Coluzzi et al., 2018) and sporadic snow were deleted pixel by pixel.NDVI outliers mentioned here appear as a sudden drop to almost zero in the growing season and do not form a sequence in this study (Komisarenko et al., 2022).To identify outliers, we iterated through every two consecutive NDVI values in the time series and calculated the difference between the second and first values for each pixel every year.We defined anomalous NDVI differences as points outside of the percentiles threshold [10 90], and if the NDVI difference is positive, then the first NDVI value used to calculate the difference will be the outlier, otherwise, the second one will be the outlier.Finally, 215 images were used to reflect seasonal transition dates in all five study periods of 2016-2020 after the quality control.Each image was resampled with 32 m spatial resolution to match the resolution of the ArcticDEM data and SnowModel outputs.To detect seasonal transition dates, we used a double sigmoid model to fit the NDVI changes on time series, and points where the curvature changes most rapidly on the fitted curve, appear at the beginning, middle, and end of each season (Klosterman et al., 2014).The applicability of this phenology method in the Arctic has been demonstrated (Ma et al., 2022;Westergaard-Nielsen et al., 2013, 2017).We focused on three seasonal transition dates, that is, SOS, NDVI max day, and EOF.The NDVI values for some pixels are still below zero in spring and summer due to topographical shadow.We, therefore, set a quality control rule before calculating seasonal transition dates for each pixel, that is, if the number of days with positive NDVI values from June to September is less than 60% of the total number of observed days, the pixel will not be considered for subsequent calculations.As verification of fitted dates, the seasonal transition dates in dry heaths and corresponding time-lapse photos acquired from the snow fence area are shown in Figure 2.
Snow cover extent is greatly reduced and vegetation is exposed with lower NDVI values on the SOS.All visible vegetation is green on the NDVI max day.On EOF, snow cover distributes partly, and NDVI decreases to a value close to zero.

| Structural equation modeling
Structural equation modeling can be used to analyze the complex networks of causal relationships in ecosystems (Fan et al., 2016).In this study, we defined climatic conditions (Table 1) and terrain conditions (elevation, slope, and aspect) as environmental factors to assess their direct and indirect effects on SOS, NDVI max day, and EOF.The original theoretical model is shown in Figure S8.The freezing index and growing degree days were calculated based on daily air temperature (0°C as threshold) and were used to determine the freezing severity of the winter season (Bilotta et al., 2015;Frauenfeld et al., 2007) and assess the impact of thermal conditions on plant growth (Wypych et al., 2017), respectively.Mean snow depth, SED, freezing index, and monthly precipitation from March to May are used to explain SOS changes.The NDVI max day appears in late July and early August in this study area, so we used the SED, growing degree days and total precipitation in May and June, and averaged wind speed in June to explain it.The climatic factors used to explain EOF, which appeared from mid-September (Figure S3), include the growing degree days and total precipitation as well as averaged monthly

TA B L E 1
The list of climatic conditions for the structural equation modeling and random forest model.

Snow onset day (SOD)
The first day of continuous accumulation of snow cover (continuous approx.7 months for most areas of the valley)

Snow ending day (SED)
The last day of continuous accumulation of snow cover (i.e., 0% snow-covered ground after this day)

Mean snow depth
The mean value of snow depth during the snow cover duration (m) Averaged monthly wind speed in June, July, August, and September (Wspd Jun., Wspd Jul. , Wspd Aug. , and Wspd Sep. ) The mean value of daily wind speed for each month (m/s) Freezing index from October to February of next year, in March, and April (FI Oct.-Feb., FI Mar. , FI Apr. ) The absolute value of the sum of air temperature below 0°C during a certain period (°C) Growing degree days in May, June, July, and August (GDD May , GDD Jun. , GDD Jul. , and GDD Aug. ) The sum of air temperature above 0°C for each month (°C) Total precipitation in March, April, May, June, July, August, and September (Pre Mar., Pre Apr. , Pre May , Pre Jun. , Pre Jul. , Pre Aug. and Pre Sep. ) The sum of daily precipitation for each month (mm) wind speed in June, July, and August, and SOD of the following snow cover season.For those areas where EOF appeared after September, we also considered the temperature, precipitation, and wind speed in September to explain it (Figure S5).The terrain factors are combined with different climatic conditions to explain each seasonal transition date.We built the direct and indirect pathways (Figure 4; S5) between all variables based on the piecewiseSEM function in R (version 4.2.3).PiecewiseSEM is more advanced in coping with complex models without the requirement of normal data than other structural equation methods (Lefcheck, 2016).The model performance was assessed based on the Fisher's C statistic, with a model being accepted if the associated p > .05(Shipley, 2000).The path model with all hypothesized pathways for each seasonal transition date was depicted in Figure S8, and here we report the final accepted model with only significant pathways.The number of samples is 8607 for the abrasion, 3102 for the fen, 26,938 for the dry heath, 21,936 for the wet heath, and 3094 for the tall shrub (corresponding to the number of analyzed pixels), which meets the minimum requirement of a sample size to analyze complex models with 100-200 observations (Fan et al., 2016;Grace et al., 2015).

| Future scenarios
Under various RCPs (Representative Concentration Pathways), the Arctic is projected to experience amplified warming and increased precipitation, particularly during autumn and winter with more intense magnitude (Bintanja & Andry, 2017;Cai et al., 2021;Overland et al., 2019).We set the study period (2016-2020) as the presentday period, and the annual mean-and autumn/winter mean-air temperature will increase by varying degrees under RCP 2.6, 4.5, and 8.5 by 2050 relative to the present-day period in the Arctic (Overland et al., 2019).Precipitation is projected to increase at about 4.5% per degree of warming (Bintanja & Selten, 2014).Thus, we assume the annual mean-and autumn/winter mean-precipitation will also in-  2).Then we input the manipulated data into SnowModel to simulate the snow cover situations after 2.5 decades.

| Random forest
Random forest regression can model complex and non-linear relationships between predictor and response variables (Hutengs & Vohland, 2016).In this study, we randomly split the present-day (2016-2020) dataset (predictor and response variables) of each land surface class into two parts, that is, 80% of the data was used to train models, and 20% was used as test data.Then we input all climatic conditions (Table 1) under six future climate scenarios, combined with terrain conditions, into the trained models to predict We observed an inter-annual variation in the three seasonal transition dates (Figure 3), that is, the standard deviation ranged from  The NDVI max day appeared in late July or early August, and the later NDVI max day corresponded to a later SOS.EOF appeared in October for most parts of the study area and showed less spatial heterogeneity than SOS.

| Controls on seasonal transition dates
The coupled climatic and terrain factors accounted for 46%-61%, 44%-54%, and 42%-58% of the changes in SOS, NDVI max day, and EOF, respectively, across the five land surface classes (Figure 4; Figure S4).Second to abrasion areas, dry heath is considered the most widespread vegetated surface class in Greenland (Karami et al., 2018).We, therefore, highlight dry heath as an example to show the intuitively direct and indirect effects of all controls on SOS, NDVI max day, and EOF separately (Figure 4).Terrain conditions affected seasonal transition dates by modulating snow distribution, air temperature, and precipitation, and the direct effects of them were weak with coefficients below |±0.1|, except for the effect of elevation on the NDVI max day and EOF with a coefficient of −0.38 and −0.19, respectively.The direct effects of precipitation from March to May on SOS were stronger among all controls with −1.03 to 1.95 (Figure 4a).Precipitation in May and June as well as SED exerted more directed effects on NDVI max day ranging from 0.3 to 0.57 (Figure 4b).EOF was more affected by summer air temperature and precipitation than SOD of the following snow cover season (Figure 4c).The effects of all controls on seasonal transition dates showed generally similar patterns in other land surface classes as in dry heaths (Figure S4), and differences in the effect of each control on seasonal transition dates between five land surface classes are illustrated in Figure 5.
In terms of total effects, climatic conditions determined changes in seasonal transition dates more than terrain conditions (Figure 5).

| Future projection
Increasing air temperature and precipitation can decrease snow depth by 28-40 cm, advance SED by 10-26 days, and delay SOD by 5-20 days between different future scenarios and land surface classes, compared to the present-day period (Figure S6).The prediction accuracy of the SOS, NDVI max day, and EOF is 0.69-0.75,0.48-0.56,and 0.69-0.80across the five land surface classes, separately (Table S1).SOS appears earlier for all land surface classes, the number of days in advance increases from scenarios of RCP 2.5 to RCP 8.5, ranging from 2 to 14 days between different land

| Effects of multi-factors on seasonal transition dates
In this study, the inter-annual variations in SOS, NDVI max day, and EOF did not exhibit similar trends, attributed to the distinct climatic controls on different seasonal transition dates.The surface reflectivity is significantly altered (Thevenard & Haddad, 2006) after the snow begins to melt, making SED a key positively associated indicator for SOS observed by remote sensing.Lower freezing indexes mean higher air temperatures and larger snow depths mean stronger snow insulation effects (Park et al., 2015), thus the combination of both could warm soil temperatures for vegetation and further advance SOS.Conversely, higher freezing indexes and smaller snow depths could delay SOS, explaining the positive direct effect of FI Oct.-Feb.(Figure 4; Figure S4) and the negative direct effect of snow depth on SOS in dry heaths.In addition, the positive total effects of snow depth and FI Oct.-Feb.on SOS stem from the positive direct effects of both on SED.Interestingly, we found that the effects of freezing index and precipitation on SOS in March are opposite to those in April (Figure 5), especially in dry heaths, wet heaths, and tall shrubs, which could be explained by considering snow conditions.
Snow began to melt after March in this study region, lower air temperatures (higher freezing indexes) in March could increase the ratio of snowfall to precipitation, preventing snow from melting thereby maintaining a low snow thermal conductivity and warmer soil (Calonne et al., 2011), contributing to the advanced SOS; higher March precipitation means more snowfall and could delay SOS by postponing SED.In April, air temperature started to rise (Figure S6), and more meltwater in the snowpack could weaken the snow insulation effect and make soil temperatures more susceptible to the fluctuations of air temperatures, thus higher freezing indexes in April delayed the SOS.In addition, more precipitation in the form of rainfall in April could speed up snow melting (Sui & Koehler, 2001) and provide water resources for the early stage of vegetation growth, partly explaining the negative effect of April precipitation on SOS.
Although the pre-spring precipitation and air temperatures showed more effects on SOS, they modulated SOS based on influencing snow properties, verifying the importance of snow conditions on

SOS.
The positive effect of SED on NDVI max day stems from the advanced SOS caused by earlier SED, which means advanced plant nutrient uptake and vegetation growth (Ernakovich et al., 2014).We found that higher GDD Jun.corresponded to advanced NDVI max day, which is consistent with previous conclusions (Chen et et al., 2016).Correspondingly, NDVI max is found to be significantly positively correlated with summer temperature in the Arctic (Vickers et al., 2016).The snow melted totally in June across most parts of the study area, and higher temperatures in the early growing season could reduce the risk of freezing damage for plants after the snow melts (Pau et al., 2011), then further promote vegetation growth.
However, more precipitation in May and June led to a later NDVI max day.That is because the precipitation amount in both months was lower than the precipitation in July and August in this study, especially during 2016-2018 with only 2-5 mm of June precipitation (Figure S6).Moreover, the precipitation in May could delay the SED due to some spring snowfall events and further delay the NDVI max day.Wind speed in June showed positive effects on NDVI max day for each land surface class.This may be because higher wind speeds result in greater vapor pressure deficits (He et al., 2022), which is negatively related to net ecosystem production and slows down vegetation growth rates (Yuan et al., 2019).In addition, an increase in wind speeds can raise evapotranspiration rates, which induces more soil water losses and inhibits plant vegetation (Wu et al., 2021).Thus we hypothesize that the negative effect of precipitation was limited and NDVI max day was dominated by SED, air temperature (Kelsey et al., 2021), and wind speed in this study.
We found the EOF was disentangled from previous winter conditions and more related to spring and summer climate.Warmer spring (higher GDD Jun. ) resulted in advanced EOF by making the NDVI max day earlier, suggesting that a warmer spring and earlier snowmelt may not extend the growing season length but shift each seasonal transition date earlier.Gamon et al. (2013) also indicated that vegetation productivity is not related to earlier snowmelt and found no trend of lengthening growing seasons in herbaceous plants of the Arctic under climate warming (Gamon et al., 2013;Zhao et al., 2015).
However, the warmer and wetter summer could delay EOF (Feng & Liu, 2015;Walker et al., 2012) because vegetation productivity is associated primarily with precipitation and soil moisture, and secondarily with air temperature in Arctic tundra ecosystems (Gamon et al., 2013).Therefore, the higher GDD Jul. and GDD Aug. combined more precipitation in July and August caused delayed EOF by promoting vegetation productivity.It was reported that the decline in wind speed significantly delayed the autumn foliar senescence over latitudes above 50° (Wu et al., 2021).However, the effects of wind speed on EOF varied between different months and land surface classes (Figure 5).The monthly (Jun-September) wind speed ranged from 1.5 to 2.6 m/s (Figure S7), thus we hypothesized such a small range of wind speed could not determine the EOF variations in this study.Sporadic snowfall events before SOD of the following snow cover season could change land surface characteristics and soil thermal environment, but the vegetation photosynthesis is resistant to cool soil temperatures resulting from autumn snow melt (Barrere et al., 2018;Stoy et al., 2022).Moreover, senescence is more synchronous than the onset of the growing season across large areas because it is cued by not only air temperature and precipitation but also photoperiod and radiation (Borner et al., 2008;Ernakovich et al., 2014).Thus, the effect of SOD for the following snow cover season on EOF was less than the effect of air temperature and precipitation during summer in this study.
The spatial heterogeneity of terrain conditions could partly result in uneven seasonal transition date distributions.Moderate slopes can be beneficial for vegetation growth, but flat and exposed snow-free areas are prone to damage to vegetation due to wind erosion and extreme temperatures during winter (Larcher et al., 2010), and further result in low soil temperatures and inhibited vegetation growth.On the contrary, more snow cover accumulates in low-lying areas, which impacts vegetation growth because of later snow ending day (Riihimäki et al., 2017), so the effect of SED on SOS was strongest in fens (Figure 5).In this study, slope and aspect had limited effects on seasonal transition dates (Figure 4; Figure S4), while elevation emerged as the predominant terrain factor for all land surface classes, except for fens in low-lying areas.We found higher elevations corresponded to delayed SOS and advanced EOF, aligning with previous studies (Piao et al., 2011;Tranquillini, 1964).It is indicated that SOS is strongly negatively related to spring growth at all elevations, and increased summer growth can result from delayed SOS at higher elevations (Mei et al., 2021), supporting our finding of earlier NDVI max days at higher elevations.However, the negative effect of elevation  & Nilsson, 2009).Essentially, the effects of elevation on seasonal transition dates were mostly from the regulation of elevation on air temperature, precipitation, and changes in snow cover (Figure 4; Figure S4).

| Differences in seasonal transition dates between land surface classes
SOS appeared earliest in tall shrubs and dry heaths (average of DOY 143), while it was latest in abrasions (average of DOY 153).
In contrast, EOF occurred the latest in tall shrubs (average of DOY 292) and earliest in abrasions (average of DOY 285).In our study area, species such as Dryas or graminoids found in abrasions responded to air temperature and snowmelt slower than species in other land surface classes, due to their slow-growing strategy.
This was reflected by their low NDVI values and can also partly explain the earliest EOF in abrasions (Figure S2).Species in fens and wet heaths are sensitive to soil moisture, and they may respond to changes in snow regimes and hydrology rather than air temperature alone (Björk & Molau, 2007), which can be reflected by the higher total effects of SED and snow depth on SOS in fens and wet heaths (Figure 5).In addition, these surface classes are likely to have higher snow depths resulting from snow accumulation in the associated depressions.Our findings highlight that tall shrubs had the earliest onset of greenness with the longest growing season among all tundra vegetation types at the site, which is in line with other studies in the Arctic (Jia et al., 2004;Kelsey et al., 2021).
The interaction between snow cover and shrubs can be one of the reasons for earlier SOS in tall shrubs (Bewley et al., 2010;Marsh et al., 2010), that is, tall shrubs buried by snow undergo 'spring-up' during melt and the snow melt rate generally increase with shrub abundance (Pomeroy et al., 2006).Although NDVI reached its maximum simultaneously in all land surface classes, their associated NDVI max values were significantly different (Figure S2), which is related to specific species growth forms in each land surface class (Raynolds et al., 2008).Generally, lower, wetter regions dominated by graminoids are more productive than higher, drier locations with a higher percentage of lichens and mosses (Gamon et al., 2013).The light-colored barren abrasion type had the lowest NDVI max values, due to its sparse vegetation.

| Future changes in seasonal transition dates
Future warming (Tokarska et al., 2020)  Fens, as one of the most endangered ecosystems, has been found that higher temperatures could result in increased gross ecosystem photosynthesis (Sullivan et al., 2008).The projected longest spring-to-fall period in fens under RCP 8.5 results from being less limited by Arctic low temperatures and beneficial from adequate soil moisture.The projected delay in EOF is not as dramatic as the projected advance in SOS for all land surface classes, especially in tall shrubs without delay in EOF under all future climate scenarios, which means EOF may be limited by other factors, such as the above-mentioned photoperiod.Although the extended shrub growing season length has been reported (Mekonnen et al., 2021), our results show that more than 3°C increase in temperature may no longer be sufficient to support the extension unless other favorable conditions are also enhanced in the Arctic.In contrast, the abrasion has the shortest spring-to-fall period although its SOS is advanced and EOF is delayed under future scenarios, which is limited by its slow vegetation growth.Nevertheless, our projections suggest an extension of the growing season length in all land surface classes, which aligns with long-term monitoring and experimental studies that demonstrated the positive response of shrub functional groups to warming, indicating their potential extension in the future (Bjorkman et al., 2020;Martin et al., 2017;Vowles & Björk, 2019).Projected shifts in vegetation distribution can alter ecosystem carbon balances by influencing intricate interactions between soil, plants, and the atmosphere (Pearson et al., 2013;Stewart et al., 2018;Yu et al., 2017).

| Uncertainty analyses
The studied climate, snow cover, and terrain conditions accounted for 44%-61% of the variations in seasonal transition dates across different land surface classes, meaning that other factors must be considered to fully understand the variations.Arctic tundra is often moisture-and nutrient-limited (Gamon et al., 2013;Huemmrich et al., 2010;Westergaard-Nielsen et al., 2021), particularly in the early growing season before peak production (Humphreys & Lafleur, 2011).Our findings support this pattern, with higher NDVI max values observed in moisture-dependent graminoids in low and moist areas compared to drought-tolerant vegetation types (prostrate shrubs, mosses, and lichens) in elevated regions.The effects of soil temperature differ between vegetation communities and plant functional types (Sloan et al., 2016).Lower soil temperatures limit the translocation of nutrients and carbohydrates from below-ground stores in shorter vegetation, and thus leaf production occurs later in the growing season (Chapin III et al., 1980;Rasmussen et al., 2022).In contrast, above-ground nutrient reserves enable rapid leaf production in tall shrubs in response to rising air temperatures (Sprugel, 2002).Because our study area is near the coast, seasonal transition dates may be influenced by sea ice dynamics on long-time series (Kerby & Post, 2013).However, variations in spring phenology are best explained by snow melt date rather than regional spring drop in sea ice extent in some Arctic sites (Assmann et al., 2019), which is consistent with our results in this study.Nevertheless, sea ice decline could affect snow cover changes in coastal areas (Callaghan et al., 2011) on a long time series, which may result in more winter precipitation than assumed in this study.These local boundary conditions are not easily integrated into the SnowModel.
As a result, the projected snow cover patterns and corresponding changes in seasonal transition dates might be more representative of areas further from the coast.
We acknowledge that the study period (

CO N FLI C T O F I NTE R E S T S TATE M E NT
The corresponding author confirms on behalf of all authors that there are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the Dryad Digital Repository at https:// doi.org/ 10. 5061/dryad.jsxks n0hp.

| 3 of 17 LIU
et al. the interactions among multi-environmental factors and their collective impacts on seasonal transition dates remain unclear.In this study, we simulated snow evolution at 32 m spatial resolution and extracted seasonal transition dates (SOS, NDVI max day, EOF) for five land surface classes from Sentinel-2 across an Arctic tundra ecosystem in Greenland.We explored the direct and indirect links between seasonal transition dates and snow cover, air temperature, precipitation, wind speed, and terrain, using structural equation modeling.Moreover, we manipulated meteorological data to mimic future climate change patterns and predicted future snow cover.Ultimately, we projected the response of seasonal transition dates to future climate changes based on random forest models trained by presentperiod data.We aim to address the following questions: (1) What are the differences in seasonal transition dates among the five tundra land surface classes?(2) How does the combination of multiple conditions affect seasonal transition dates directly or indirectly?(3) How do seasonal transition dates of different land surface classes respond to future climate change?
. As one of the most important inputs of the used snow model, underestimated precipitation needs to be corrected.Precipitation is measured by Pluvio in our study area, which automatically records precipitation amounts per hour at a gauge height of 2 m.Wolff et al. (2015) developed a nonlinear relationship between under-catch and wind speed during winter precipitation events and associated correction F I G U R E 1 (a) Map of the study area and ground-based observations distribution.The region circled by black lines is the study area.Blue dots represent the meteorological stations that provide input data for the used models and green dots are three snow fences used to validate snow depth simulations.(b) The distribution of land surface classes in the study area.
a) The change of median normalized difference vegetation index (NDVI) across the dry heath in the study area in 2020.(b)-(d) show the corresponding photos of the start of spring (SOS), maximum NDVI (NDVI max ) day, and end of fall (EOF) of 2020 in the well-drained mesic tundra heath where snow fences are set up (69°18′ N, 53°30′ W, green dots in Figure1a).
crease by varying degrees under different RCPs by 2050.We set up two future scenarios based on the present-day climate to mimic the possible climate conditions by 2050 for each RCP: (1) increase in air temperature and precipitation across the entire year (scenario I); (2) increase more air temperature and more precipitation from September to February, and stay same to the scenario I during the rest of months (scenario II).Finally, six future climate change scenarios were produced (Table

|
the seasonal transition dates by 2050.All predictor variables in the random forest analysis are consistent with environmental factors in the structural equation modeling when predicting the SOS, NDVI max day, and EOF.Seasonal transition dates variations Seasonal transition dates varied between different land surface classes (Figure 3).During the five growing seasons, the annual mean SOS was the latest (day of the year [DOY] 153) in abrasions and ranged from DOY 144 to 146 in the other four land surface classes.The annual mean NDVI max day occurred on DOY 209 for all five land surface classes although the NDVI max values ranged from 0.35 in abrasions to 0.54 in tall shrubs (Figure S2).Abrasions had the earliest EOF at DOY 285, whereas tall shrubs had the latest at DOY 292, and the remaining three land surface classes experienced it at DOY 291.

| 7 of 17 LIU
et al. in fens to 13 days in abrasions for EOF, respectively.Thus the SOS exhibited the most drastic inter-annual variation in comparison to the NDVI max day and EOF.The spatial distribution of seasonal transition dates is illustrated in Figure S3.Spring started in May across the valley and in June across the northeast and northwest hillsides.
SED, snow depth, and FI Oct.-Feb.positively affected SOS in all land surface classes.Increased FI Mar. with reduced March precipitation, and decreased FI Apr. with increased April precipitation, and greaterMay precipitation could advance SOS in all land surface classes except for the fen, where these conditions had adverse or no effects.On the other hand, GDD May showed a positive total effect on SOS with 0.35 in fens, but few effects in other land surface classes.Higher elevation caused delayed SOS, with total effects ranging from 0.09 to 0.22, while slope negatively impacted SOS (−0.18 to −0.03) in all land surface classes except for the fen.The NDVI max day was primarily influenced by SED, May and June precipitation, and averaged June wind speed, all with positive total effects, and by GDD Jun. with negative effects in all land surface classes.EOF was mainly affected by air temperature and precipitation during spring and summer (Figure5).Higher GDD Jun.corresponded to earlier EOF, but warmer July and August with more August precipitation could postpone EOF for all land surface classes, despite varying effects of June and July precipitation.Monthly wind speeds from June to August exerted varying effects on EOF between different land surface classes, but these effects were lower compared to air temperatures and precipitation.The SOD of the following snow cover season has a certain positive effect on EOF with 0.17-0.44across all land surface classes except for the fen.Increased elevation corresponded to earlier NDVI max day and EOF in all land surface classes except for the fen.

F
Changes of the mean (a) start of spring (SOS), (b) maximum normalized difference vegetation index (NDVI max ) day, and (c) end of fall (EOF) across the study area during 2016-2020.The error bars are one standard deviation of the annual means.surfaceclasses (Figure6a).The NDVI max day remains the same in fens but advances by 2-9 days in abrasions, and it is delayed by 2 days in dry heaths, by 2-3 days in wet heaths, and by 4-8 days in tall shrubs under six future scenarios (Figure6b).Under all future scenarios except Scenario I RCP2.6 , EOF is postponed by 1-6 days in abrasions and 1-5 days in dry heaths.Under Scenario I RCP2.6 , Scenario II RCP2.6 , and Scenario I RCP4.5 , EOF advances by 1-3 days or remains unchanged in fens and wet heaths, but it delays by 3-7 days in fens and 1-2 days in wet heaths under the remaining three scenarios.Tall shrubs are expected to experience one day advance in EOF under Scenario I RCP2.6 , Scenario II RCP2.6 , and Scenario I RCP4.5 but no changes under the other three scenarios (Figure6c).The spring-to-fall period is projected to extend by 5-20 days in abrasions, 1-21 days in fens, 4-17 days in dry heaths, 1-13 days in wet heaths, and 2-12 days in tall shrubs compared to 2016-2020, with greater extensions corresponding to warmer and wetter climates.The spring-to-fall period in abrasions is projected to be the shortest (132-147 days) among all land surface classes, aligning with the duration observed in 2016-2020.Under Scenario I RCP2.6 , Scenario II RCP2.6 , and Scenario I RCP4.5 , tall shrubs have the longest spring-tofall period, lasting 142-145 days, but which becomes the second longest with 146-151 days under the remaining three scenarios.The dry heath has a spring-to-fall period of 139-145 days, longer than wet heaths and fens under Scenario I RCP2.6 , Scenario II RCP2.6 , and Scenario I RCP4.5 , but the fen has the longest spring-to fall period with 148-156 days under the remaining three scenarios, although fens, dry heaths, and wet heaths had a similar seasonal transition pattern during 2016-2020.

F
I G U R E 6 Comparison of the mean (a) start of spring (SOS), (b) maximum normalized difference vegetation index (NDVI max ) day, and (c) end of fall (EOF) between the present-day climate, and six future scenarios for five land surface classes.Numbers attached to each subfigure indicate the differences in SOS, NDVI max day, and EOF between future scenarios and the present-day period.
and increasing winter precipitation(McCrystall et al., 2021) have been forecasted in highlatitude areas, which cause changes in winter soil temperature and snow melting time.In this study, the projected seasonal transition dates varied between five land surface classes under different future climate scenarios.Although all climate scenarios suggested decreased snow depths, advanced SED, and delayed SOD compared with the 2016-2020 period (FigureS6), the snow depth are expected to be higher, with later SED and earlier SOD in scenario I than in scenario II (warmer and wetter autumns and winters), especially under RCP 4.5 and 8.5, highlighting the complex responses of snow cover to various magnitudes of future warming and wetting, which has been verified by large-scale snow cover projections(Callaghan et al., 2011;Mudryk et al., 2020) and further contributes to the changes in seasonal transition dates.The projected spring-to-fall period is the longest in tall shrubs (142-145 days) under Scenario I RCP2.6 , Scenario II RCP2.6 , and Scenario I RCP4.5 , but in fens (148-156 days) under Scenario II RCP4.5 , Scenario I RCP8.5 , and Scenario II RCP8.5 , although the abrasion is projected to have the most advance in SOS and delay in EOF, attributing to differences in response mechanisms to climate change and vegetation activities between various vegetation communities.
5 years) is insufficient to capture long-term trends of seasonal transition dates.However, the climatic conditions had significant inter-annual variations, thus we could observe the responses of seasonal transition dates to climate changes.Meanwhile, we analyzed interactions between multiple factors and projected future changes on a pixel scale, thus we have enough samples to support our conclusions.The need for high-resolution NDVI data from Sentinel-2 and validation of snow simulations limits the study period to start from 2016 at the earliest and end in 2020 at the latest.5 | CON CLUS IONS The spatiotemporal variations in seasonal transition dates (SOS, NDVI max day, and EOF) were revealed in this study for five land surface classes in the Arctic during 2016-2020.The climatic and terrain conditions can explain 44%-61% of the changes in seasonal transition dates across different land surface classes.Our results suggest that the combination of earlier SED, warmer air temperature, and more precipitation in summer could potentially contribute to a longer growing season while increasing elevation had the opposite effect.Under the wetter and warmer Arctic climate, we projected that the spring-to-fall period would extend for all land surface classes due to advanced SOS or delayed EOF.The length of spring-to-fall in fens will be the longest under RCP 4.5 and 8.5, which is the opposite of patterns during 2016-2020.This research demonstrates the necessity of analyzing seasonal transition dates from more perspectives, such as selecting different predictor variables for various transition dates, to understand the responses of surface phenology to climate change.Conceptualization; data curation; formal analysis; methodology; software; validation; writing -original draft; writing -review and editing.Peiyan Wang: Methodology; software; writing -review and editing.Bo Elberling: Funding acquisition; methodology; resources; writing -review and editing.Andreas Westergaard-Nielsen: Conceptualization; formal analysis; funding acquisition; methodology; resources; software; supervision; writing -review and editing.ACK N OWLED G M ENTS This work was supported by the Danish National Research Foundation (CENPERM DNRF100), the VILLUM FOUNDATION grant 42069, the National Natural Science Foundation of China (32201360), and the China Scholarship Council (CSC).We acknowledge the infrastructural support by Arctic Station (University of Copenhagen), Greenland Ecosystem Monitoring (GEM) for climate data, and SnowModel developed by Glen Liston.We thank the journal reviewers for their helpful comments.
days in fens to 17 days in abrasions for SOS, from 13 days in tall shrubs to 15 days in abrasions for NDVI max day, and from 11 days Increases in temperature and precipitation by 2050 under six climate change scenarios.
Note: Spring and summer range from March to August, and autumn and winter include the remaining months of the year.
Freezing index Oct.-Feb., Freezing index Mar., and Freezing index Apr.], growing degree days in May, June, July, and August [GDD May , GDD Jun. , GDD Jul. , and GDD Aug. ], total precipitation in March, April, May, June, July, and August [Pre Mar., Pre Apr. , Pre May , Pre Jun. , Pre Jul. , and Pre Aug. ], averaged wind speed in June, July, and August [Wind speed Jun., Wind speed Jul., and Wind speed Aug. ], elevation, slope, Standardized total effect of environmental factors (snow depth [SD], snow ending day [SED], snow onset day [SOD], freezing index from October to February and in March and April [FI Oct.-Feb., FI Mar. , and FI Apr.], growing degree days in May, June, July, and August [GDD May , GDD Jun. , GDD Jul. , and GDD Aug. ], total precipitation in March, April, May, June, July, and August [Pre Mar., Pre Apr. , Pre May , Pre Jun. , Pre Jul. , and Pre Aug. ], averaged wind speed in June, July, and August [Wspd Jun., Wspd Jul. , and Wspd Aug. ], al., 2018; Xu F I G U R E 5 on NDVI max day was slight in this study because NDVI reached the maximum value in late July and early August across different elevations.Moreover, the response of growing season length to climate change can be more sensitive in high-elevation areas with an overall faster-extending rate under a warming climate (Wilson