Observed earlier start of the growing season from middle to high latitudes across the Northern Hemisphere snow-covered landmass for the period 2001–2014

Vegetation phenology in spring has received much attention for its importance to terrestrial ecosystem carbon exchange and climate–biosphere interactions studies. Through control on surface water and heat balance, snow cover largely impacts on spring vegetation phenology. However, under the background of global warming and rapid reduction of spring snow cover extent across the Northern Hemisphere (NH), the responses of spring vegetation phenology have not been well documented. Using two satellite-derived land cover dynamic datasets and 420 in situ vegetation phenology observations from five filed datasets, this study evaluated the accuracy of satellite-derived vegetation phenology datasets and explored the changes of start of the growing season (SOS) across the NH snow-covered landmass for the period 2001–2014. Compared with MEaSUREs VIPPHEN, the MODIS SOS maps displayed higher accuracy in capture the real SOS climatology by validating with in situ observations (R2 = 0.67, bias = −3.99 d). Moreover, evidences from MODIS SOS maps pointed out that the SOS advanced by approximately 2.36 d in NH middle to high latitudes (43.5°N–70.0°N), but delayed by about 0.53 d in lower latitudes (33.0°N–43.5°N) from 2001 to 2014. The contrast SOS anomalies across the NH snow-covered landmass were further proved by changes in spring NDVI derived from GIMMS in the corresponding period. In addition, the observed changes in SOS were consistent with the spatiotemporal pattern of spring snow end date found in previous studies, indicating vegetation phenology changes should be taken into account in estimating the impacts of snow in climate–biosphere interactions studies.


Introduction
The start of the growing season (SOS) in spring is an essential component of environmental science influencing terrestrial ecosystem atmosphere carbon exchange (Joiner et al 2014), snow-albedo feedback calculation (Loranty et al 2014), energy budget estimation (Peñuelas et al 2009), dry deposition (Martin et al 2014), fire disturbance (Goetz et al 2005), heat waves (Lorenz et al 2013), and climate-biosphere interactions studies (Richardson et al 2012) at hemispheric and regional scales. As the largest single component of the cryosphere, seasonal snow cover obviously impacts on spring SOS through its influences on surface energy and moisture fluxes, e.g. the dynamic of snow cover is considered as an essential factor in phenological changes in Arctic tundra (Zeng and Jia 2013), Nepal Trans Himalaya (Paudel and Andersen 2012) and Tibetan Plateau (Wang et al 2013). Recent studies have proven that the snow cover extent has experienced a well-documented rapid decrease since satellite observations began in later 1960s, and climate projections suggest that the snow cover extent will continue to shrink in the future accompanied with global warming Robinson 2011, Stocker et al 2013), especially in the spring season. However, due to highly spatial heterogeneity of snow cover and vegetation distribution, poor quality of snow cover and vegetation dynamic datasets, as well as difficult in collection of in situ phenology observations, the response of spring SOS to snow cover changes remain unknown.
Global climate change has emerged as a major driver of spring SOS in terrestrial ecosystems ecosystem (Julien and Sobrino 2009, Forkel et al 2015, Fu et al 2015, Wang et al 2019. For example, water availability control on global interannual variability and trends in land surface phenology and greenness (Forkel et al 2015). Compared with water availability, land surface temperature plays a greater role in the phenological changes of vegetation, e.g. Krishnaswamy et al (2014) reported that the spring greenness was strongly associated with the observed increase in temperature over the world's mountain ecosystems located in the pan-tropical belt (30°N-30°S) for the period 1982-2006. However, the global warming hiatus has weakened the relationship between spring vegetation and temperature. Xu et al (2013) found that temperature and vegetation seasonality diminishment over northern lands for the period 1982-2010. Fu et al (2015) demonstrated that the apparent response of leaf unfolding to climate warming has significantly decreased from 1980 to 2013 in seven dominant European tree species. Wang et al (2019) reported that no trends in spring and autumn phenology during the global warming hiatus between 1998 and 2012.
Among the global vegetation phenology changes, the spring SOS anomalies in Northern Hemisphere (NH) have drawn great attentions for its large spatial distribution. A number of studies have quantified spring SOS conditions across the NH and explored its response to climate change in recent decades. At continental scales, Schwartz et al (2006) found that the onset of spring started earlier across most temperate NH land regions during 1955. Jeong et al (2011 concluded that the SOS advanced by 5.2 d in the early period (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) but advanced by only 0.2 d in the later period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) over the NH. Wang and Fensholt (2017) reported that temporal changes in coupled vegetation phenology and productivity are biome-specific in the NH, which indicates complex relationships and interactions between vegetation phenology and productivity. At regional scales, vegetation phenology changes in Fennoscandia  (Menzel and Fabian 1999, Bogaert et al 2002, Fu et al 2014 have been separately studied. However, due to different spatial coverage and time span, a robust and coherent continental scale account on spring SOS changes both using satellite and ground observation has been lacking.
Knowledge of spring SOS and its response to climate change over the snow-covered landmass is essential to climate-biosphere interaction studies and ecosystem predictions. Although, snow cover largely impacts on spring vegetation phenology, few of abovementioned studies have quantified the response of spring vegetation phenology to snow cover changes. Therefore, it is necessary to explore the response of spring SOS to climate change and snow cover extent reduction under this background. For this purpose, this study estimated the accuracy of two satellitederived land cover dynamic datasets using in situ observations and qualified the changes in SOS across the NH. Section 2 gives a brief review of the datasets and methods used in this study. The results and analysis are presented in section 3. The discussion and conclusions are given in sections 4 and 5.
2. Data and methods 2.1. Study area Land cover changes may led to uncertainties in contributing SOS anomalies across the NH. To focus on changes in spring SOS, this study confined the study area to vegetation-covered landmass over the NH within the stable snow-covered regions. The stable snow-covered regions are described in Chen et al (2016). The grid cells with non-vegetated, urban, water, permanent ice and snow are excluded from this study. To eliminate the impacts of land cover change from vegetated to non-vegetated, we generated a multiyear stable vegetation-covered regions from the annual Moderate Resolution Imaging Spectroradiometer (MODIS) land cover data (MCD12C1) from 2001-2014 and used it as the boundary of study area. According the International Geosphere Biosphere Program (IGBP) land cover classification system from MCD12C1, the land surface is divided into 17 types, including 11 natural vegetation types, 3 land use and land mosaic types, and 3 vegetation-free land types (Friedl et al 2010). The stable vegetated land cover map over the NH snow-covered landmass derived from MCD12C1 for the period 2001-2014 is shown in figure 1.

MODIS land cover dynamic dataset
Remotely sensed vegetation phenology datasets have been much touted for their potential to provide spatially explicit, temporally rich information on vegetation dynamics and patterns at landscape at regional to global scales (Konrad et al 2011). To identify spring SOS changes from 2001 to 2014, the MODIS Global Vegetation Phenology product (MCD12Q2) Version 6 were employed in this study. The MCD12Q2 algorithm identifies phenophase transition dates based on logistic functions fit to time series of the 2-band Enhanced Vegetation Index (EVI2) (Huete et al 2002) and provides estimates of the onset of greenness, green-up midpoint, maturity, peak greenness, senescence, green-down midpoint, dormancy, as well as overall and phenology metricspecific quality information at global scales with 500 m spatial resolution. In this study, the onset of greenness derived from MCD12Q2 was used as MODIS SOS maps in the following context, in which the date when EVI2 first crossed 15% of the segment EVI2 amplitude was defined as onset of greenness in spring (Huete et al 2002). Comparison of vegetation phenology retrieved from MODIS with in situ measurements shows that these metrics provide realistic estimates of the SOS according Zhang et al (2006).

MEaSUREs VIPPHEN dataset
Compared with MODIS SOS maps, the newly developed National Aeronautics and Space Administration (NASA) Making Earth System Data Records for Use in Research Environments (MEaSUREs) Vegetation Index and Phenology (VIPPHEN) global datasets provides another choice for SOS studies. MEaSUREs VIPPHEN were created using surface reflectance data from the Advanced Very High Resolution Radiometer (AVHRR) N07, N09, N11, and N14 datasets (1981-1999) and MODIS surface reflectance data (2000-2014) (Didan and Barreto 2016). MEaSUREs VIPPHEN provides estimates of vegetation phenological metrics such as the start, peak, and end of season as well as the rate of greening and senescence globally at the spatial resolution of 0.05°. In addition, a reliability band is included to provide context on the quality of the input data. To compare with MODIS SOS maps, the VIPPHEN based on the time series of EVI2 was employed in this study. In addition, the start of season 1 derived from MEaSUREs VIPPHEN was used as MEaSUREs SOS map in the following context, in which the start of season 1 in VIPPHEN is defined as the point in time for which the EVI2 value has increased by 10% of the distance between the minimum and maximum values (Jönsson and Eklundh 2004). According Didan et al ( 2018), the date dependent phenology parameters error was 5-30 d with an average error of ±15 d. For its long term series and complete spatial coverage, MEaSUREs VIPPHEN has been widely used in Phenology, climate change and urbanization studies (Beresford et al 2018, Qiu et al 2020).

GIMMS NDVI dataset
To cross-compare with SOS maps derived from MODIS and MEaSUREs, the Normalized Difference Vegetation Index (NDVI) derived from Global Inventory Monitoring and Modeling System (GIMMS) was employed in this study. The GIMMS NDVI is widely in measuring global vegetation phenology for its long time span and large spatial coverage

Ground vegetation phenology observations
To verify the performance of MODIS and MEaSUREs SOS maps in capture the distribution of 'real' greenness onset date, five in situ vegetation phenology data sets are used in this study, including the Canada Nature Watch (CNW), the Chinese Ecosystem Research Network (CERN), the Pan European Phenology Project (PEP), the Lilac leafing observations (Lilac) and the United States of America (USA) National Phenology observations (NPN).

Canada nature watch
The CNW was launched in 2000 as a partnership between Environment Canada, the environmental Non-Governmental Organization Nature Canada, and several other organizations, with the aim of getting the Canadian public to help researchers track changes in the natural environment. Data collected through CNW has been used in numerous scientific reports and studies over the past few years. The Plant Watch module of CNW, which got its start at the University of Alberta in 1995, has generated multiple articles in scientific journals describing changing climatic conditions in Canada (Beaubien and Hamann 2011). In this study, 122 stations covering from 2001 to 2014 across the Canada were selected to validate the performance of satellite-retrieved SOS maps.

Chinese ecosystem research network
The plant phenological observation dataset of the CERN is the integration of plant phenological observation data from 20 CERN ecological stations of more than 660 plant species (Song et al 2017). In this study, the beginning of leaf unfolding phase of woody plant and turning green phase of herbaceous plant from 6 stations were used as SOS dates to validate the performance of satellite-retrieved SOS maps.

Pan European Phenology project
The PEP is a project funded by the Zentralanstalt für Meteorologie und Geodynamik, the Austrian ministry for science and research, and the network of European meteorological services (EUMETNET), with the goal to establish an open access database with plant phenology data sets for science, research and education (Templ et al 2018). In this study, 162 stations with first leafing records ranging from 2001 to 2014 over the Pan European were selected to validate the perforce of satellite-retrieved SOS maps.

Lilac leafing observations
The Lilac leafing observations collected across the continental United States from 1956 to 2014 for purple common lilac, a cloned lilac cultivar and two cloned honeysuckle cultivars (Rosemartin et al 2015). Applications of this observational dataset range from detecting regional weather patterns to understanding the impacts of global climate change on the onset of spring at the national scale. Lilac and honeysuckle phenology data have proven robust in both model development and climatic research. In this study, 36 stations ranging from 2001 to 2014 across the United States were employed to validate the perforce of satellite-retrieved SOS maps.

USA national phenology network
The USA-NPN seeks to engage a diverse range of citizen scientist volunteers, federal, state, and nongovernmental organizations, educators and professional research scientists to collect phenological observations of plants and animals using consistent standards and to contribute their observations to a national data repository (Schwartz et al 2012). The aim of USA-NPN are developing a diversity of materials, tools, techniques, and protocols to assist decision making and education related to ecology, wildlife, human health, ecosystem services, natural resource management, biological conservation, and climate change adaptation. In this study, 94 NPN ground phenology stations ranging from 2001 to 2014 across the United States were used to validate the perforce of satellite-retrieved SOS maps. Details of the above 5 in situ datasets are listed in table 1. Stations with incomplete records, missing values, low latitudes and altitudes without seasonal snow cover were excluded in this study. The distribution of selected 420 in situ vegetation phenology observations are displayed in figure 2. In consideration of the definitional differences between five in situ datasets, the leafing out date in CNW, the beginning of leaf unfolding in CERN, the first leafing date in PEP, Lilac, and NPN records were used as in situ SOS observations.

Data processing
Before calculation, all of the satellite-retrieved datasets used in this study were regridded at a spatial resolution of 0.05°in a geographic latitude/longitude projection with the help of MRT and gdalwarp (https://gdal. org/). To obtain accurate results of MODIS SOS maps and compare with other datasets used in this study, the MODIS SOS in 500 m spatial resolution were aggregated to product a 0.05°SOS maps in a geographic latitude/longitude projection. To generate robust MODIS SOS observations, only selected quality control flags of onset of greenness with '1', representing good quality, were used in the resample process. In addition, to cross-compare SOS maps derived from MODIS and MEaSUREs, the GIMMS NDVI in spring season (March-April-May, MAM) with a spatial resolution of 8 km for the period 2001-2014 were employed to produce a robust MAM-averaged spring NDVI with a 0.05°spatial resolution in a geographic latitude/longitude projection. In consideration of the differences in spatial coverage, the analysis was carried out in gridcells with valid values in MCD12C1, MCD12Q2, MEaSUREs VIPPHEN, and GIMMS NDVI data sets.

Comparison between satellite-retrieved SOS maps
Comparison with higher spatial-resolution images is widely used in the validation of moderate-resolution satellite images, such as by Hall et al (1995). In this study, the comparison between MODIS and MEa-SUREs SOS maps was carried out for the period 2001 to 2014. The root mean square error (RMSE) and bias were used as criteria to evaluate the differences between MODIS and MEaSUREs SOS maps. The RMSE and bias of the MEaSUREs to MODIS SOS maps are expressed as follows:

Evaluation by in situ SOS observations
In addition to compare between satellite-retrieved SOS maps, validating moderate-resolution satellite images by field measurements is another widely used approach in previous studies. However, validation satellite images by field measurements is difficult because a single grid cell from satellite measurements usually contain the information from an extremely large area, which may overestimate or underestimate the information from a field measurement. Even though, the field vegetation phenology records are still the most convincing data to test the reliability of satellite-retrieved products. To estimate the accuracy of MODIS and MEaSUREs SOS maps, the evaluation by ground SOS observations also carried out in this study. The coefficient of determination (R 2 ) and bias are employed as criteria to estimate the ability of MODIS and MEaSUREs in capture the distribution of in situ SOS records.       0.05°, which cannot catch the 'real' value of SOS in specific spot location. In addition, the different definitions of SOS in MODIS and MEaSUREs datasets may also contribute to the different performance in comparison with in situ observations. Since new vegetation phenology dataset with high spatial resolution and large spatial coverage is not available at this moment, the SOS uncertainties caused by different definition is beyond the scope of this study.

Comparison between MODIS and MEaSUREs SOS maps
Previous research has proved that the raw in situ observations would give results that are highly dependent on the particular locations and reporting periods of the actual weather stations. Such results would mostly reflect those accidental circumstances rather than yield meaningful information about the climate (Hansen et al 2010). Even though, as shown in figure 4, the MODIS SOS maps still skillfully capture in situ SOS distributions over the NH snow-covered landmass. Thus, in the following context, we will use MODIS in the SOS change analysis.
3.3. Observed changes in spring SOS over the NH Using aggregated MODIS SOS dataset, the SOS changes across the NH snow-covered landmass for the period 2001-2014 are shown in figure 5. Changes are expressed as linear trends multiplied by the time interval in this study. The SOS date over most parts of the NH high latitudes (60°N-70°N) was advanced by approximately 4.18 d in average, especially in high latitude Eurasia where SOS advanced by more than 15 d at the 95% significance level. Compared with North America, the most notable advance of SOS occurred over Eastern Europe and Western Asia ( figure 5(a)). Moreover, our results pointed out a contrast SOS changes across the northern middle to high latitudes during 2001-2014. The SOS has delayed by about 1.13 d in average across the mid-latitudes of Central Europe, Central North America and Eastern Asia.
In addition, there were large differences in zonally averaged SOS anomalies between regions located in northern middle and high latitudes ( figure 5(b)). As shown in figure 5(b), the sign of the zonally averaged changes in SOS turns at around 43.5°N. The average SOS changes below 43.5°N was 0.53 d, with the most obvious delay occurred around 41.2°N (2.56 d). In contrast, the average SOS changes between 43.5°N-70°N was −2.36 d, with maximum advances distributed around 70°N (−6.35 d).
The observed contrasting changes in SOS between the NH middle and high latitudes are highly similar to the spatiotemporal distribution of anomalies in spring NDVI derived from GIMMS over the corresponding period ranging from 2001 to 2014 (figure 6(a)). Moreover, there are generally negative correlation between the SOS and spring NDVI over the studied area ( figure 6(b)). In response to the earlier SOS observed in high latitudes, spring NDVI increased due to vegetation growth. On the other hand, delayed SOS observed in mid-latitudes, resulted in decrease of spring NDVI. (±39.08) in DOY, respectively. Quantifications of land surface phenology and greenness dynamics are impaired by differences between satellite data sets and phenology detection methods (Forkel et al 2015). Thus, the differences in SOS definition and low spatial resolution of MEaSUREs dataset may explained this discrepancies. Although both MCD12Q2 and MEa-SUREs VIPPHEN used MODIS surface reflectance data as primary input data for the period 2000-2014, they were developed by using Huete et al (2002) and White et al (2009), with different methodology and EVI2 threshold. However, since the new vegetation phenology dataset with high spatial resolution and large spatial coverage is not available at this moment, the SOS uncertainties caused by different definition is beyond the scope of this study.

Discussion
The validation of satellite-retrieved SOS maps by using in situ observations shown that MODIS can generally catch the 'real' SOS distributions. Compared with MODIS, there are greater dispersion between MEaSUREs SOS dates and in situ observations. However, the five in situ datasets are affiliated to different organizations. Since there are no overlap locations of the five in situ datasets, it is difficult to estimate the consistencies between each data sets. Even though, as the 'real' SOS distributions in validating of two satellite-retrieved SOS maps, the differences in SOS definitions of five in situ observations will not influence the results.
This study pointed out a contrast changes in the SOS distributed in the NH middle latitudes for the period 2001-2014, in which there are advanced SOS distributed in NH high latitudes and delayed SOS detected in NH middle-latitudes. The SOS dates over most parts of the NH high latitudes (60-70°N) was advanced by approximately 4.18 d in average, especially in high latitude Eurasia, which is consistent with the general greening trend occurred in Europe (Fu et al 2014), Fennoscandia (Høgda et al 2013) and Siberia (Forkel et al 2015). Meanwhile, the SOS has delayed by about 1.13 d in average across the midlatitudes of Central Europe, Central North America and Eastern Asia. As proved in previous studies, global climate change has emerged as a major driver of spring SOS in terrestrial ecosystems ecosystem. The Arctic amplification (Cohen et al 2014, Francis et al 2017, Pacific Ocean surface cooling (Kosaka and Xie 2013) and delayed snow end date in spring (Chen et al 2015) resulted in declined spring surface temperature in NH mid-latitudes, which is also contributing to the delayed SOS detected in the NH mid-latitudes.

Conclusions
Under the background of global warming and snow cover extent reduction over the NH, the changes of spring vegetation phenology have not been well documented. In this study, two satellite-retrieved SOS datasets were evaluated by using five in situ observations across the NH snow-covered landmass. The changes in spring SOS were explored and further cross-compared with an independent spring NDVI for the period 2001-2014.
Comparisons between MODIS and MEaSUREs SOS maps for the period 2001-2014 displayed large bias in NH high latitudes around 50°N. Compared with MEaSUREs, the MODIS SOS maps displayed higher accuracy by validating with in situ observations (R 2 =0.67, bias=−3.99 d). Moreover, evidences from MODIS SOS maps pointed out contrast spring SOS anomalies occurred in NH snow-covered landmass, in which the SOS advanced by approximately −2.36 d in NH high latitudes (43.5°N-70°N) but delayed by about 0.53 d in mid-latitudes (33.0°N-43.5°N) from 2001 to 2014. The contrast SOS anomalies across the NH snow cover landmass were further proved by changes in spring NDVI derived from GIMMS in the corresponding period. Moreover, the observed changes in SOS were consistent with the spatiotemporal pattern of snow end date during the spring season found in previous studies (Chen et al 2015), indicating vegetation changes should be taken into account in estimating the impacts of snow in climate-biosphere interactions studies. Figure 6. 14 year (a) changes in spring NDVI calculated from GIMMS and (b) correlation coefficient between spring NDVI and MODIS SOS over the NH from 2001 to 2014. Black dots in (a) indicate that changes are statistically significant at the 95% level. Black dots in (b) indicate that the linear correlation between SOS and spring NDVI is statistically significant at the 95% level.
Compared with previous research, the present study evaluated the accuracy of two representative satellite retrieved SOS datasets and explored the changes in SOS over the NH snow covered landmass, which are important for current climate change studies and future climate projections. However, several unresolved issues remain in our study. Compared with in situ SOS observations from ground meteorological records, satellitebased SOS maps exhibit advantages in terms of the large spatial coverage and high temporal resolution. However, the different definition in SOS may result in incomparable between different datasets, which largely restrict the application of remote sensing technology in vegetation phenology studies. Although this study estimate the accuracy of satellite retrieved SOS across the NH, the 14 year period is still too short to fully understand SOS changes under climate change. Thus, high spatial resolution, long time span snow cover products are expected in future vegetation phenology studies.