Evaluating Biosphere Model Estimates of the Start of the Vegetation Active Season in Boreal Forests by Satellite Observations

The objective of this study was to assess the performance of the simulated start of the photosynthetically active season by a large-scale biosphere model in boreal forests in Finland with remote sensing observations. The start of season for two forest types, evergreen needleand deciduous broad-leaf, was obtained for the period 2003–2011 from regional JSBACH (Jena Scheme for Biosphere–Atmosphere Hamburg) runs, driven with climate variables from a regional climate model. The satellite-derived start of season was determined from daily Moderate Resolution Imaging Spectrometer (MODIS) time series of Fractional Snow Cover and the Normalized Difference Water Index by applying methods that were targeted to the two forest types. The accuracy of the satellite-derived start of season in deciduous forest was assessed with bud break observations of birch and a root mean square error of seven days was obtained. The evaluation of JSBACH modelled start of season dates with satellite observations revealed high spatial correspondence. The bias was less than five days for both forest types but showed regional differences that need further consideration. The agreement with satellite observations was slightly better for the evergreen than for the deciduous forest. Nonetheless, comparison with gross primary production (GPP) determined from CO2 flux measurements at two eddy covariance sites in evergreen forest revealed that the JSBACH-simulated GPP was higher in early spring and led to too-early simulated start of season dates. Photosynthetic activity recovers differently in evergreen and deciduous forests. While for the deciduous forest calibration of phenology alone could improve the performance of JSBACH, for the evergreen forest, changes such as seasonality of temperature response, would need to be introduced to the photosynthetic capacity to improve the temporal development of gross primary production.


Introduction
The boreal zone is characterized by strong seasonality of vegetation.Climate change is expected to influence this seasonal cycle and thus annual carbon balances of terrestrial vegetation are also subject to change [1].In Finland, the growing season is predicted to lengthen in a warming climate with an earlier start in spring and a later cessation in autumn [2].Rising temperature enhances both carbon uptake from the atmosphere by photosynthesis (Gross Primary Production, GPP) and release of carbon to the atmosphere from biomass and soil carbon by autotrophic and heterotrophic respiration.Higher temperatures will enhance the forest carbon uptake in spring [3,4], while in autumn the ecosystem respiration may continue longer at a high level, reducing the annual net uptake of carbon by vegetation [5].Winter temperatures have been projected to rise even more than the mean annual temperatures [6], which may disturb the dormancy development of trees [7].The decrease in the number of frost days is most pronounced in the southern boreal zone, while, especially in the northerly regions, alternating freezing and melting periods may become more frequent [8], increasing the risk of frost damage [9].Correspondingly, frequent severe springtime night frosts that constrain photosynthetic recovery may counteract the positive response of photosynthesis and phenological development to warmer spring temperatures and may have a long-term effect on the forest carbon balances [10].Snow depth and snow-covered area will decrease in the warming climate [8,11].
In process-based modelling of the carbon cycle of vegetation, photosynthesis is constrained by the Leaf Area Index (LAI), whose development follows a modelled phenological cycle.Thus, any shortcomings in the modelled phenological cycle have marked impact on the predicted GPP [12][13][14].Currently, the description of the seasonal development of LAI needs improvement in global earth system models [15] and more validation and calibration for the boreal zone, especially for deciduous broad-leaf forests (DBF) [14,16].While local comparisons can be made with micrometeorological measurements of the CO 2 exchange at the ecosystems scale [14,[16][17][18]], a regional approach is needed to capture the spatial variability, and remote sensing is a useful tool for studying vegetation phenology over large areas [19].Benchmarking biosphere models e.g., [20], with remotely sensed phenology would provide steps to a more accurate presentation of the annual carbon balance and its response to climate variability [15,21].
Remote sensing of the vegetation phenology in the boreal region is hindered by frequent prolonged overcast periods, low sun elevation in autumn and winter and a short vegetation period with rapid phenological changes [22].Furthermore, a remote sensing method suitable for extracting the start of season of canopy photosynthesis capacity and GPP in one ecosystem does not necessarily work for all vegetation types [23].In DBF, the leaf development in spring is a prerequisite for the start of photosynthetic activity and strong correlations are observed between early season phenology (leaf out date, start of LAI development) and the onset of daily net carbon uptake [16,24].In contrast to DBF, photosynthesis in evergreen needle-leaf forests (ENF) increases gradually already before the production of new foliage and thus the seasonal recovery of photosynthesis has been difficult to observe in remote sensing observations [25,26].
Several remote sensing methods have been developed for detecting the start of the growing season in boreal areas [22,[27][28][29].These methods track the greening of the vegetation within the footprint of the satellite pixel and have been evaluated against in situ phenological observations of deciduous trees [27,28,30] or phenological indicators determined from micrometeorological CO 2 flux measurements [29].Because in ENF the onset of photosynthesis occurs before the production of new foliage, remotely sensed start of season, tracking the development of new leaves and shoots at the land surface scale, showed a delay of several weeks when compared to the onset of photosynthesis [16,23,29,31].Large shifts, of up to three months, between the greening phase from three global remotely sensed data sets of fraction of absorbed photosynthetically active radiation (fAPAR) were obtained for ENF in the boreal region [32].Instead, a proxy indicator of start of season, determined from satellite time series of Fractional Snow Cover (FSC), was found suitable for this forest type in boreal areas [33].
In this study, we use simulations by the biosphere model JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg) [34] to determine the start of season in DBF and ENF during the past decade in Finland.The start of season is here defined as the beginning of the photosynthetically active period of vegetation.JSBACH is the biosphere component of an Earth system model with a global scope and has, therefore, not been validated for regional applications.Evaluating such a model regionally is challenging, as there are hardly any suitable datasets available at this scale.Even though these large-scale models have been targeted to global scale, it is important to assess their performance at a regional scale to gain understanding about their applicability and characteristics that require further development.The carbon and water cycle of vegetated areas are modelled in JSBACH by considering the physiological response of vegetation to the climatic variables.In this work, JSBACH is run with climatic variables from a regional climate model.Moreover, to assess the sensitivity of the model to the climate drivers, JSBACH is run also with observed meteorological data from micrometeorological measurement sites in Finland.
We use remote sensing observations for the evaluation of JSBACH-modelled start of season of boreal forests.Because of structural and functional differences in the photosynthetic recovery in DBF and ENF, specific methods will be employed to extract the start of season date from remote sensing data in the respective forest type.We apply the method developed by Delbart et al. [27] and the method based on snow cover [33] for DBF and ENF, respectively.Both methods have been used previously for the mapping of the start of season in boreal areas.For DBF, we also assess here the accuracy of the remote sensing observations against ground observations of bud break of birch trees in Finland.We assess the spatial and temporal correspondence of JSBACH results at the regional scale with these satellite observations.In addition, we evaluate the simulated start of season dates for ENF locally with micrometeorological CO 2 flux measurements.

Study Area
The study area covers Finland and ranges from 19 ˝E to 32 ˝E, and from 59 ˝N to 71 ˝N (Figure 1).According to Köppen's climate classification, the area belongs to the temperate coniferous-mixed forest zone with cold and wet winters.The average temperature of the warmest month exceeds 10 ˝C and the average temperature of the coldest month is below ´2 ˝C [35].Finland's vegetation can be divided into four phytogeographical zones: the hemiboreal, southern, middle and northern boreal zones [36] (Figure 1).The main land cover type is forest (with 44% coverage of evergreen, deciduous and mixed forest), followed by shrubs and herbaceous vegetation, inland waters and wetlands (Supplementary Materials, Table S1).ENF is the dominant forest type, except for areas north of 68 ˝N where DBF is dominating (Figure 1b).The main tree species of forests are Scots pine (Pinus sylvestris) (67%), Norway spruce (Picea abies) (22%), and birch (10%), including silver birch (Betula pendula) and downy birch (Betula pubescens).More than half of the forests are pure stands (55%); stands with some mixing cover 31% of the forested area, and 13% of forest stands are mixed [37].

Phenological Observations of Bud Break of Birch
We used phenological observations of the bud break of birch (Betula pubescens Ehrh.) from 8 sites for the evaluation of satellite observations of start of season in DBF (DBF_SOSsat) in Finland.The sites are located in the southern (3), middle (1) and northern boreal zone (4) (Figure 1, Table 1)

Carbon Dioxide Flux Measurements
The net ecosystem CO 2 exchange (NEE) between atmosphere and vegetation can be directly measured by the micrometeorological eddy covariance (EC) method [38].The 30-min GPP value was here obtained by subtracting the modelled respiration (R) value from the NEE observation (GPP = NEE-R) utilizing standard flux partitioning procedures [39,40].By the same parameterisation, the NEE and GPP time series were gap-filled for comparison with the model results.In this study, we used EC data from two boreal Scots pine forests, Sodankylä in northern Finland (67 ˝21'N, 26 ˝38'E) [41] and Hyytiälä in southern Finland (61 ˝51'N, 24 ˝17'E) [42] (Figure 1a), from years 2001-2010.Supporting meteorological data from both flux sites were used in the site level JSBACH model simulations (Section 2.4.4).The starting date of the vegetation active season (ENF_SOS EC ) was determined for the two sites from the continuous EC data.ENF_SOS EC was obtained as the day on which the CO 2 uptake permanently exceeds 15% of the growing season maximum uptake level.The method has been described in more detail in a previous study by Böttcher et al. [33].

Phenological Observations of Bud Break of Birch
We used phenological observations of the bud break of birch (Betula pubescens Ehrh.) from 8 sites for the evaluation of satellite observations of start of season in DBF (DBF_SOS sat ) in Finland.The sites are located in the southern (3), middle (1) and northern boreal zone (4) (Figure 1, Table 1).Bud break observations were obtained from the Finnish Forest Research Institute (now the Natural Resource Institute Finland) for the period from 2001 to 2008.The phenological observations at Aulanko and Saariselkä had already stopped in 2006.According to the observer's instructions [43], bud break is observed when at least 50% of buds have broken throughout the tree crown.Phenological observations were carried out twice a week at the respective sites.The monitored trees (five) are of local provenance and were preferably naturally regenerated.The microclimate of the site must be representative for the local climate conditions.Monitored trees may have been changed from year to year until 2007.Further details on observations are given in Kubin et al. [43]. 1 Pixels with the sum of deciduous and mixed forest fractions >50% at the original spatial resolution of 0.005 were used for the averaging of the Normalized Difference Water Index time series to the coarser grid; 2 Altitude was calculated from a digital elevation model provided by the National Land Survey of Finland with a grid size of 25 m by 25 m and a height precision of 2 m; 3 SD, standard deviation.

The JSBACH Model
The biosphere model JSBACH [34,44] is used for the modelling of GPP.It is part of the Max Planck Institute's (MPI) Earth System Model and implemented as a subroutine of its atmospheric component ECHAM [45].JSBACH is used in global simulations, but it can also be applied at regional [46][47][48] and at site level scales [49], as is done in this work.JSBACH is a process-based model and calculates the exchanges of carbon, water and energy between the land surface and the atmosphere.
The phenology model (Logistic Growth Phenology; LoGro-P) of JSBACH provides leaf area index (LAI) at each time step to compute photosynthetic production.The model describes the phenological dynamics based on three subsequent phases; a growth (spring), a vegetative (summer) and a rest phase (autumn and winter).For evergreen and deciduous forest, the phenological cycle is dependent on air temperature.For the shift from dormancy to the growth phase in spring, the "alternating model" of Murray et al. [50] is employed in JSBACH.It consists of a heat sum approach (also often referred to as growing degree day approach) combined with a chilling condition.The counting of chilling days starts from the last autumn event (shift from vegetative phase to dormancy phase).The critical temperature, above which temperature contributes to heat sum and below which chilling days are counted, was set to 4 ˝C in both forest types.A detailed description of the JSBACH LoGro-P phenological model and the parameter values used in this work is given in Appendix A. The shift to the growth phase corresponds to the start of development of new shoots in ENF and to bud break in DBF.Because ENF have needles present throughout the year (and non-zero LAI in the model), photosynthesis can start before the shift to the growth phase and defines the start of season in ENF.The photosynthetic uptake increases with temperature and is enhanced by the increase of LAI after the beginning of the growth phase.
The photosynthesis model in JSBACH is based on the biochemical model developed by Farquhar et al. [51], as are many of the large-scale models.The temperature responses of photosynthesis model parameters describing the potential electron transport rate and maximum carboxylation rate stay the same throughout the year.The photosynthesis is coupled to the stomatal conductance, which is constrained at low soil water content [52].As a part of the water and energy balances, the model simulates snow cover.Fractional Snow Cover is calculated according to Roesch et al. [53].
Different vegetation types are described as Plant Functional Types (PFTs).Each PFT has specific physiological and physical properties, such as photosynthesis parameters, aerodynamic roughness and canopy albedo.Four different PFTs were described for each grid cell in our simulations and based on the Coordinated Information on the European Environment (CORINE) Land Cover data.In this work, only the tiles containing ENF and DBF are considered.As ENF is the most abundant PFT in Finland, for the vast majority of grid cells there was representation of ENF, and thus an estimate was available for the whole Finnish territory except for the very northern areas.Nevertheless, the cover of grid cells with DBF representation was even more complete than that of ENF because it is everywhere among the four most abundant PFT's and the distribution of DBF extends further in northern Finland.

Preparation of Meteorological Driving Data for Regional Runs
In the set-up used, JSBACH requires seven mutually consistent meteorological variables (2-m air temperature and specific humidity, 10-m wind velocity, short-wave and long-wave radiation, potential short-wave radiation and precipitation) in hourly time resolution.In order to produce such a data set for a regional domain it is necessary to use a climate model of suitable scale.For the JSBACH regional driving data production we used the regional climate model REMO [54,55] with an implementation of surface properties by Gao et al. [56].REMO resolves three-dimensional (3D) atmospheric variables such as air temperature, specific humidity and wind speed in a hybrid vertical coordinate system [57].The model also produces surface layer values of the atmospheric variables (2-m air temperature and 10-m wind velocity) as well as surface and soil variables such as snow cover and soil moisture.In this study, REMO was run for a domain covering Fennoscandia, with a grid resolution of 0.167 ˝, corresponding approximately to 18 km.ERA-Interim, a global atmospheric re-analysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) [58], was used as lateral boundary data for the 3D atmospheric variables.The forward run covered the period from 1979 to 2011.The forward run was preceded by a 10-year spin-up to equilibrate soil moisture and temperature.
Model-specific biases for distinct regions are inherent for regional climate models [59,60].REMO exhibits colder temperatures than observed in winter and spring in Eastern Europe [61], and the area of cold bias covers most of Finland [56].In summer and autumn the temperatures are too warm compared to observations.REMO overestimates precipitation in northern Europe throughout the year [56].In order to remove the biases, we corrected air temperature and precipitation against gridded homogenized weather data for 1980-2011 provided by the Finnish Meteorological Institute (FMI) [62].A quantile-quantile type bias correction algorithm was applied for daily mean temperature [63], while daily cumulative precipitation was corrected using parametric quantile mapping [64].Finally, the daily corrections were applied to the hourly modelled air temperature and precipitation values.To differentiate between the original REMO data and the bias-corrected REMO data, we will use the term regional climate model data when referring to the bias-corrected meteorological drivers in the following.
2.4.3.Regional JSBACH Runs JSBACH was driven using hourly regional climate model data obtained from REMO for the years 1980-2011.A spin-up period of 30 years was applied before the forward run to equilibrate soil water and soil temperature.The spatial domain of JSBACH is the same as that of the driving model, REMO, covering Fennoscandia.For the present study only Finnish territory is investigated.

Site-Level JSBACH Runs
Additionally, JSBACH was run at site level for the two Scots pine EC flux sites, Sodankylä and Hyytiälä, using both regional climate model data and local site observations as meteorological input data.Regional climate model data, described in the previous section, was taken from the closest grid point to the site in the regional domain for the site-level simulations.The local meteorological data set, including the same variables, were gap-filled from the in situ meteorological measurements, except for the theoretically derived potential short-wave radiation.The site-level simulations for ENF were done with the PFT coniferous evergreen trees prescribed at both EC flux sites.The PFT for DBF was not considered in the site-level runs, because there are no continuous EC measurements in deciduous forest in Finland.

Model-Derived Indicators for the Start of Season
The start of season in ENF for the regional JSBACH run (ENF_SOS mod ) as well as for the local site-level runs (ENF_SOS mod,loc ) was calculated from the canopy GPP of the respective PFT fraction.The method of calculating the ENF_SOS mod was similar to that for ENF_SOS EC (see Section 2.2).The JSBACH modelled start of season of DBF (DBF_SOS mod ) was the first day of the year when the LAI exceeded the tiny threshold of 0.01 m 2 ¨m´2 .This date typically coincides with the transition to the growth phase in LoGro-P.
Since we used a proxy indicator from satellite time series of FSC for ENF_SOS sat , we also included JSBACH modelled daily FSC in our assessment and calculated the proxy indicator from it using the method described in Section 2.5.2.The JSBACH FSC of fully snow-covered ground is by definition smaller or equal to 0.95 [53] and depends on terrain variation within the grid cell, whereas the FSC from satellite observation equals one in the case of fully snow-covered ground.Thus, to apply the remote sensing method, we scaled the modelled FSC based on the maximum simulated value for each grid cell during the period from 2003 to 2010, which is considered to correspond to full snow cover.
In addition to ENF_SOS mod and DBF_SOS mod described above, the date of the start of the JSBACH "growth phase" (Section 2.4.1) of DBF and ENF was used in this study.We compared the JSBACH growth phase date of DBF simulated with a regional climate model and FMI gridded temperature for the analysis of remaining bias in the start of season due to the driving temperature.For this comparison, the grid cells at the edge of the area with FMI temperature observations were excluded.

Remote Sensing Data
We used daily Moderate Resolution Imaging Spectrometer (MODIS) time series for the detection of the start of season in ENF and DBF.The pre-processing of the MODIS time series and the extraction of start of season dates will be described in the following sections.

Pre-Processing of MODIS Time Series
Terra MODIS Level-1B data were obtained from the National Aeronautics and Space Administration (NASA)'s Level 1 and Atmosphere Archive and Distribution System (LAADS) and from the receiving station of the FMI in Sodankylä for the period from February to October in 2003 to 2010.The pre-processing consisted of calculation of Top-of-Atmosphere (TOA) reflectance and geometric projection of images into a latitude/longitude grid (datum WGS-84).MODIS observations with solar zenith angles greater than 73 degrees and sensor zenith angles greater than 50 degrees were excluded.
For the detection of DBF_SOS sat , we calculated the Normalized Difference Water Index (NDWI) from MODIS TOA reflectance data at 0.005 ˝spatial resolution (corresponding approximately to a pixel size of 280-190 m (at 60 ˝N and 70 ˝N) by 550 m) as: where R NIR and R MIR refer to near-infrared and mid-infrared reflectance in MODIS band 2 (841-876 nm) and MODIS band 6 (1628-1652 nm), respectively.For the detection of ENF_SOS sat , we calculated the FSC based on the SCAmod algorithm by Metsämäki et al. [65,66] at the same spatial resolution as the NDWI.Clouds in daily FSC and NDWI products were masked based on an automatic cloud masking algorithm developed at the Finnish Environment Institute.

Satellite-Observed Start of Season for JSBACH Model Evaluation
For the later comparison of satellite observations with modelled start of season, we spatially aggregated the daily time series of FSC and NDWI to a grid resolution of 0.25 ˝by 0.25 ˝.Because the heterogeneity of land cover is one main factor for discrepancies between satellite-observed phenology and ground observations [67,68], we focussed on satellite pixels with dominant cover of ENF and DBF according to CORINE Land Cover data of the year 2000 [69] (see also Figure 1b).In the aggregation, the FSC time series were averaged for MODIS pixels with a dominant fraction (ě70%) of coniferous forest.NDWI time series were averaged for pixels with dominance of deciduous and mixed forest (sum of both vegetation types was ě70%).
A sigmoid function was fitted to the FSC time series and ENF_SOS sat was detected when the fraction of snow cover decreased from full snow cover, FSC = 1, to FSC = 0.99 as follows: where x 1 is the inflection point (day when FSC reaches 0.5) and x 2 controls the rate of change of the sigmoid function.The full details of the method are given in Böttcher et al. [33].The accuracy of ENF_SOS sat has been evaluated with ENF_SOS EC at three CO 2 flux measurement sites in Finland.The root mean square error (RMSE) was 5.3 days and the bias 1 day (late) [33].The quality of the estimate was similar for the southern and the northern boreal sites, but it may be less reliable for the hemiboreal zone in winters with low snow depth [33].
For the detection of the DBF_SOS sat , we followed the method by Delbart et al. [27], which is based on the NDWI.This index is useful for the detection of the start of season in boreal areas, because snow melt and greening-up of vegetation can be differentiated from the time series [70], whereas the observed spring increase in the widely applied Normalized Difference Vegetation Index (NDVI) can be largely influenced by reflectance change in the visible range due to the melting of snow, and thus the increase in NDVI corresponds actually to the time of snow melt in boreal areas [71,72].The method for determination of DBF_SOS sat from NDWI was first applied in central Siberia [27] and later to boreal Eurasia [73].DBF_SOS sat is the date when the vegetation within the footprint of the satellite pixel is greening up, which corresponds to the time of NDWI increase after snow melt [27].Comparison with in situ observations of leaf appearance of birch, larch and aspen in central Siberia with estimates from SPOT-VEGETATION NDWI time series showed an RMSE of 6.7 days and a negligible bias [70].
In this study, we smoothed averaged daily NDWI time series (spatial resolution 0.25 ˝) with the adaptive Savitzky-Golay filter [74,75] that is implemented in the software TimeSat [76].DBF_SOS sat was then determined from the smoothed time series as the last day when NDWI is lower than the minimum annual NDWI increased by 5% of the NDWI spring amplitude: where NDWI min is the minimum NDWI and ε equals 5%.The original threshold value ε of 20% that accounts for the noise level in the time series [27] was reduced to 5% because of the lower noise level in the averaged and smoothed NDWI time series here.Satellite-observed start of season dates were not determined for areas with scarce coverage of the respective vegetation type within a MODIS pixel.Furthermore, long temporal gaps in satellite time series in spring increase the uncertainty in the detection of SOS from satellite observations.Therefore, we excluded grid cells with a temporal gap larger than two weeks around the estimated start of season date for the comparison with JSBACH model results.

Evaluation of DBF_SOS sat against Bud Break Observations
As described in the previous section, the method for the detection of DBF_SOS sat was validated earlier in Siberia [70] for boreal forests with a different species composition than in Finland.As described in Section 2.5.2, we made small changes to the original method.Therefore, we used bud break observations of birch from phenological sites (Section 2.3) for evaluation of DBF_SOS sat in Finland.The DBF_SOS sat at the coarse grid for the evaluation of JSBACH model results is not appropriate for comparison with ground observation.Thus, the original MODIS NDWI time series with a pixel size of 0.005 ˝by 0.005 ˝were aggregated to 0.05 ˝by 0.05 ˝taking into account pixels with a dominance of deciduous and mixed forest (>50% coverage, according to National CORINE Land Cover) and excluding pixels that were partially (ě10%) covered by water bodies.DBF_SOS sat was extracted from the nearest grid cell (at 0.05 ˝resolution) to the phenological observation sites.The number of MODIS pixels that were used in the averaging of the NDWI time series is given in Table 1.The difference between the median of elevation within the grid cell and the elevation of the phenological site was less than 50 m (Table 1).For only two sites, Kevo and Koli, the variation in elevation within the pixel was larger than 50 m.
We analysed the correspondence between satellite and ground observations for all site-years, for each site and for the southern, middle and northern boreal zone, separately.Furthermore, we investigated the correspondence of the year-to-year anomalies between DBF_SOS sat and the bud break observations.The statistical metrics used in the comparison were the coefficient of determination (R 2 ), the root mean square error (RMSE) and the bias.The equations are given in Appendix B (Table B1).We excluded year 2002, because of missing data in the MODIS archive during the spring period.

Evaluation of JSBACH-Modelled Start of Season
For the comparison of model-and satellite-derived start of season dates, JSBACH results were projected to the same grid as the satellite SOS (0.25 ˝, see Section 2.5.2).We investigated correspondence between the dates for the whole of Finland and for the boreal sub-regions: the southern, middle and northern boreal zone (Figure 1a), and calculated Pearson correlations, RMSE and bias for the period from 2003 to 2010 separately for the two forest types (ENF, DBF).In order to assess correspondence regarding spatial patterns, we used the multi-year mean dates for the whole period for each grid cell.In addition, agreement in temporal patterns was investigated based on the year-to-year anomalies.
For site-level comparisons at Sodankylä and Hyytiälä, the ENF_SOS mod was extracted from the closest JSBACH grid cell.We calculated Pearson correlations, the RMSE and bias between ENF_SOS EC and ENF_SOS mod,loc , and between ENF_SOS EC and ENF_SOS mod .We also assessed the agreement between the respective year-to-year anomalies for the period 2001 to 2010.

Evaluation of DBF_SOS sat with Phenological Observations of Bud Break
Comparisons of DBF_SOS sat with bud break observations including all site-years resulted in an RMSE of one week and an early bias of 0.13 days (R 2 = 0.75, probability p < 0.001, Number of observations N = 48) (Figure 2a).The site-wise correlation between DBF_SOS sat and in situ observations were not significant (Table 2), but the number of observations (4-7 per site) was also very small.RMSE and bias ranged between 4 to 10 days and 0 to 5 days, respectively.The smallest RMSE and bias were obtained for northern sites, Saariselkä and Kevo, in Lapland.The largest RMSE (10 days) and bias (5 days early) were both found for Kannus in western Finland.Large deviations (>two weeks) between satellite and in situ observations were observed for Vesijako in 2005 (15 days) and Kannus in 2007 (19 days).
For the southern boreal zone, including sites Aulanko, Vesijako and Koli, the correlation of DBF_SOS sat with bud break observations was weak (R 2 = 0.32, p = 0.012, N = 19).The RMSE was 6.4 days and the bias 1.8 days (early).For the four sites in northern Finland: Värriö, Saariselkä, Kolari and Kevo, the correlation was higher than for the sites in the southern boreal zone (R 2 = 0.62, p < 0.001, N = 22), but the RMSE (6.9 days) and the bias (2.8 days late) were larger.The correlation between interannual anomalies of DBF_SOS sat and bud break was significant (Figure 2b), but could explain only less than 30% of the interannual variation in the field observations.

Evaluation of JSBACH-Modelled Start of Season
For the comparison of model-and satellite-derived start of season dates, JSBACH results were projected to the same grid as the satellite SOS (0.25°, see Section 2.5.2).We investigated correspondence between the dates for the whole of Finland and for the boreal sub-regions: the southern, middle and northern boreal zone (Figure 1a), and calculated Pearson correlations, RMSE and bias for the period from 2003 to 2010 separately for the two forest types (ENF, DBF).In order to assess correspondence regarding spatial patterns, we used the multi-year mean dates for the whole period for each grid cell.In addition, agreement in temporal patterns was investigated based on the year-to-year anomalies.
For site-level comparisons at Sodankylä and Hyytiälä, the ENF_SOSmod was extracted from the closest JSBACH grid cell.We calculated Pearson correlations, the RMSE and bias between ENF_SOSEC and ENF_SOSmod,loc, and between ENF_SOSEC and ENF_SOSmod.We also assessed the agreement between the respective year-to-year anomalies for the period 2001 to 2010.

Evaluation of DBF_SOSsat with Phenological Observations of Bud Break
Comparisons of DBF_SOSsat with bud break observations including all site-years resulted in an RMSE of one week and an early bias of 0.13 days (R 2 = 0.75, probability p < 0.001, Number of observations N = 48) (Figure 2a).The site-wise correlation between DBF_SOSsat and in situ observations were not significant (Table 2), but the number of observations (4-7 per site) was also very small.RMSE and bias ranged between 4 to 10 days and 0 to 5 days, respectively.The smallest RMSE and bias were obtained for northern sites, Saariselkä and Kevo, in Lapland.The largest RMSE (10 days) and bias (5 days early) were both found for Kannus in western   We compared ENF_SOS EC in Sodankylä and Hyytiälä with JSBACH-simulated dates based on the two different sets of meteorological drivers: ENF_SOS mod,loc and ENF_SOS mod (Figure 3).ENF_SOS EC ranged from mid-March to mid-April (33 days) at Hyytiälä, and from the end of April to the beginning of May (14 days) at Sodankylä.The ENF_SOS mod,loc was too early in comparison with observations (Table 3).The bias was twice as large for Sodankylä as compared to Hyytiälä and correlations between ENF_SOS mod,loc and ENF_SOS EC were not significant.For the southern boreal zone, including sites Aulanko, Vesijako and Koli, the correlation of DBF_SOSsat with bud break observations was weak (R 2 = 0.32, p = 0.012, N = 19).The RMSE was 6.4 days and the bias 1.8 days (early).For the four sites in northern Finland: Värriö, Saariselkä, Kolari and Kevo, the correlation was higher than for the sites in the southern boreal zone (R 2 = 0.62, p < 0.001, N = 22), but the RMSE (6.9 days) and the bias (2.8 days late) were larger.The correlation between interannual anomalies of DBF_SOSsat and bud break was significant (Figure 2b), but could explain only less than 30% of the interannual variation in the field observations.

Site-Level Start of Season in Evergreen Needle-Leaf Forest
We compared ENF_SOSEC in Sodankylä and Hyytiälä with JSBACH-simulated dates based on the two different sets of meteorological drivers: ENF_SOSmod,loc and ENF_SOSmod (Figure 3).ENF_SOSEC ranged from mid-March to mid-April (33 days) at Hyytiälä, and from the end of April to the beginning of May (14 days) at Sodankylä.The ENF_SOSmod,loc was too early in comparison with observations (Table 3).The bias was twice as large for Sodankylä as compared to Hyytiälä and correlations between ENF_SOSmod,loc and ENF_SOSEC were not significant.The bias between ENF_SOS EC and the start of season from the regional run, ENF_SOS mod , was much smaller than the bias between ENF_SOS EC and ENF_SOS mod,loc .Again, the early bias was larger for Sodankylä than for Hyytiälä.Correlations between ENF_SOS mod -modelled dates with ENF_SOS EC were significant only for Hyytiälä (Table 3), whereas the year-to-year anomalies correlated significantly for both sites (R 2 = 0.51, p = 0.001, N = 19).
An investigation of the reasons for differences in ENF_SOS mod and ENF_SOS mod,loc in Sodankylä showed that GPP simulated with JSBACH followed closely the development of air temperature during the springs, 2004 and 2008 (Figures 4 and 5).However, the simulated GPP showed a higher increase during warm spells in spring than the GPP from CO 2 flux measurements.This was even more pronounced for the JSBACH GPP simulations with the local driving data than for simulations with regional climate data (Figure 4a,b).The local temperature observations indicated warm spells in 2004 around the day of year (doy) 96, 108 and 119.Therefore, the simulated GPP exceeded the threshold for the start of season earlier than the observed GPP.Consequently, the JSBACH-simulated ENF_SOS mod,loc and ENF_SOS mod occurred earlier than ENF_SOS EC , on doy 96 and 117, respectively.Instead, ENF_SOS EC was delayed by cold spells that occurred around doy 110 and 115, and it was only observed on doy 121.Modelled snow depth from JSBACH simulations with both local and regional meteorological drivers was lower than the observed snow depth and the snow melted earlier in the model (Figures 4 and 5c).
The difference between JSBACH ENF_SOS mod and ENF_SOS mod,loc was smaller in 2008 than in 2004 (only 5 days) (Figures 3 and 5).ENF_SOS mod occurred on the same day as ENF_SOS EC observations (doy 123).Again, the local mean daily air temperatures were higher around the time of the start of season than the regional climate model temperatures, thus resulting in a higher GPP increase in JSBACH simulations with local drivers.Similarly to 2004, the snow melted earlier in JSBACH model simulations than indicated by snow depth measurements at Sodankylä (Figure 5c).The simulated GPP seemed to follow more closely the observed GPP at the site after the observed start of season date than during the time of spring recovery of photosynthesis.ENF_SOSEC was not available for 2002 at Sodankylä, therefore only 9 years were used in the analysis.
The bias between ENF_SOSEC and the start of season from the regional run, ENF_SOSmod, was much smaller than the bias between ENF_SOSEC and ENF_SOSmod,loc.Again, the early bias was larger for Sodankylä than for Hyytiälä.Correlations between ENF_SOSmod-modelled dates with ENF_SOSEC were significant only for Hyytiälä (Table 3), whereas the year-to-year anomalies correlated significantly for both sites (R 2 = 0.51, p = 0.001, N = 19).
An investigation of the reasons for differences in ENF_SOSmod and ENF_SOSmod,loc in Sodankylä showed that GPP simulated with JSBACH followed closely the development of air temperature during the springs, 2004 and 2008 (Figures 4 and 5).However, the simulated GPP showed a higher increase during warm spells in spring than the GPP from CO2 flux measurements.This was even more pronounced for the JSBACH GPP simulations with the local driving data than for simulations with regional climate data (Figure 4a,b).The local temperature observations indicated warm spells in 2004 around the day of year (doy) 96, 108 and 119.Therefore, the simulated GPP exceeded the threshold for the start of season earlier than the observed GPP.Consequently, the JSBACH-simulated ENF_SOSmod,loc and ENF_SOSmod occurred earlier than ENF_SOSEC, on doy 96 and 117, respectively.Instead, ENF_SOSEC was delayed by cold spells that occurred around doy 110 and 115, and it was only observed on doy 121.Modelled snow depth from JSBACH simulations with both local and regional meteorological drivers was lower than the observed snow depth and the snow melted earlier in the model (Figures 4 and 5c).The difference between JSBACH ENF_SOSmod and ENF_SOSmod,loc was smaller in 2008 than in 2004 (only 5 days) (Figures 3 and 5).ENF_SOSmod occurred on the same day as ENF_SOSEC observations (doy 123).Again, the local mean daily air temperatures were higher around the time of the start of season than the regional climate model temperatures, thus resulting in a higher GPP increase in JSBACH simulations with local drivers.Similarly to 2004, the snow melted earlier in JSBACH model simulations than indicated by snow depth measurements at Sodankylä (Figure 5c).The simulated GPP seemed to follow more closely the observed GPP at the site after the observed start of season date than during the time of spring recovery of photosynthesis.

Regional Assessment of JSBACH-Modelled Start of Season against Satellite Observations
The earliest dates of ENF_SOSsat occurred in the south-western coastal areas and the latest dates were found in the north-east, ranging from the end of March to the beginning of May (Figure 6).DBF_SOSsat occurred about one month later than ENF_SOSsat and ranged from the beginning of May to mid-June.For the regional assessment of JSBACH SOS, we first investigated the spatial agreement between the multi-year averages of ENF_SOSmod and DBF_SOSmod for the period 2003-2010 with satellite observations.The modelled mean dates were highly correlated with the respective remote sensing observations (Figure 7, Supplementary Materials: Table S2) showing, however, an early bias for ENF_SOSmod and a late bias for DBF_SOSmod.The RMSEs were five and six days for ENF_SOSmod and DBF_SOSmod, respectively.For DBF_SOSmod, this is a bit larger than the accuracy of the remote sensing estimates when considering the site-means (Table 2).

Regional Assessment of JSBACH-Modelled Start of Season against Satellite Observations
The earliest dates of ENF_SOS sat occurred in the south-western coastal areas and the latest dates were found in the north-east, ranging from the end of March to the beginning of May (Figure 6).DBF_SOS sat occurred about one month later than ENF_SOS sat and ranged from the beginning of May to mid-June.For the regional assessment of JSBACH SOS, we first investigated the spatial agreement between the multi-year averages of ENF_SOS mod and DBF_SOS mod for the period 2003-2010 with satellite observations.The modelled mean dates were highly correlated with the respective remote sensing observations (Figure 7, Supplementary Materials: Table S2) showing, however, an early bias for ENF_SOS mod and a late bias for DBF_SOS mod .The RMSEs were five and six days for ENF_SOS mod and DBF_SOS mod , respectively.For DBF_SOS mod , this is a bit larger than the accuracy of the remote sensing estimates when considering the site-means (Table 2).
For the boreal sub-regions, the highest spatial correlations between the mean start of season dates for the period 2003-2010 was observed in the middle boreal zone (R 2 values ě 0.79) for both, ENF_SOS mod and DBF_SOS mod , while R 2 was low in the northern boreal zone (<0.5).The bias varied regionally; it increased for ENF_SOS mod and decreased for DBF_SOS mod from the southern to the northern boreal zone (Supplementary Materials: Table S2).For ENF, it was largest (ě10 days) in the south-western coastal areas (hemiboreal zone, late bias of ENF_SOS mod ) and south-western Lapland (early bias of ENF_SOS mod ) (Figure 8a).It has to be noted that the satellite estimate based on FSC may be less reliable for the hemiboreal zone in winters with low snow depth [33].
Only a few grid cells in the hemiboreal zone and the most northern grid cells showed a larger bias than 10 days for DBF_SOS mod (Figure 8b) and the bias was smallest for the transient region between the middle boreal and northern boreal zone.
Throughout Finland, the multi-year average ENF_SOS mod took place 20-40 days earlier than the simulated date of the JSBACH start of "growth phase" of ENF.The time lag between these two phenological events decreased from south-west to north-east and increased again for northern areas of Lapland.
In the second step, we looked into year-wise comparisons of the start of season dates (Figure 9, Supplementary Materials: Table S3).The correspondence between the annual dates from JSBACH and from satellite observations showed large variations in R 2 , RMSE and bias.As shown in Figure 9a,c, for many adjacent grid cells ENF_SOS mod occurred on the same day (clumping), whereas we found a more gradual change from south to north in ENF_SOS sat .In order to assess the reasons for clumping, we analysed the change in 5-day mean temperature, which is a good predictor of spring recovery in ENF, according to Suni et al. [18], and the 5-day mean temperature often exceeded zero degrees concurrent with ENF_SOS mod .Thus, a change in weather pattern may simultaneously trigger the start of season in a large area.Clumping was not observed for DBF_SOS mod (Figure 9; lower panel).For the boreal sub-regions, the highest spatial correlations between the mean start of season dates for the period 2003-2010 was observed in the middle boreal zone (R 2 values ≥ 0.79) for both, ENF_SOSmod and DBF_SOSmod, while R 2 was low in the northern boreal zone (<0.5).The bias varied regionally; it increased for ENF_SOSmod and decreased for DBF_SOSmod from the southern to the northern boreal zone (Supplementary Materials: Table S2).For ENF, it was largest (≥10 days) in the south-western coastal areas (hemiboreal zone, late bias of ENF_SOSmod) and south-western Lapland (early bias of ENF_SOSmod) (Figure 8a).It has to be noted that the satellite estimate based on FSC may be less reliable for the hemiboreal zone in winters with low snow depth [33].
Only a few grid cells in the hemiboreal zone and the most northern grid cells showed a larger bias than 10 days for DBF_SOSmod (Figure 8b) and the bias was smallest for the transient region between the middle boreal and northern boreal zone.
Throughout Finland, the multi-year average ENF_SOSmod took place 20-40 days earlier than the simulated date of the JSBACH start of "growth phase" of ENF.The time lag between these two phenological events decreased from south-west to north-east and increased again for northern areas of Lapland.
In the second step, we looked into year-wise comparisons of the start of season dates (Figure 9, Supplementary Materials: Table S3).The correspondence between the annual dates from JSBACH and from satellite observations showed large variations in R 2 , RMSE and bias.As shown in Figure 9a,c, for many adjacent grid cells ENF_SOSmod occurred on the same day (clumping), whereas we found a more gradual change from south to north in ENF_SOSsat.In order to assess the reasons for clumping, we analysed the change in 5-day mean temperature, which is a good predictor of spring recovery in ENF, according to Suni et al. [18], and the 5-day mean temperature often exceeded zero degrees concurrent with ENF_SOSmod.Thus, a change in weather pattern may simultaneously trigger the start of season in a large area.Clumping was not observed for DBF_SOSmod (Figure 9; lower panel).The application of the method used for detection of ENF_SOS sat (Section 2.5.2) relating to the JSBACH-simulated snow cover was not possible for many years (2005, 2007, 2008 and 2009) in the hemiboreal, southern boreal zone and partly also in the middle boreal zone, because the simulated scaled FSC was below 1 during the whole spring period.For the remaining areas with sufficient snow cover, the modelled decrease of FSC occurred in mean earlier than ENF_SOS sat (8 days) for the period 2003-2010, which is also confirmed by local comparison to snow depth observations in Sodankylä (Figures 4c and 5c).The application of the method used for detection of ENF_SOSsat (Section 2.5.2) relating to the JSBACH-simulated snow cover was not possible for many years (2005, 2007, 2008 and 2009) in the We verified whether the driving temperature contributed to the differences between DBF_SOS sat and DBF_SOS mod .Correlations between the multi-year average dates of the start of the "growth phase" of DBF for the period 2003-2010 from LoGro-P simulations with two temperature drivers (regional climate model and FMI gridded temperature) were very high and the remaining bias for the whole of Finland was very small (R 2 = 0.98, RMSE = 1.66,Bias = ´0.05days).Nevertheless, the bias varied spatially (Figure 8c); late biases (1-2 days) occurred in central Finland and early biases were observed in north-eastern Finland (´2-´4 days).
The last part of this Section focuses on the inter-annual variations in the start of season dates.On the whole, inter-annual variations of JSBACH model-derived dates were in agreement with satellite observations (Figure 10), except in 2007, when DBF_SOS mod was 8 days delayed compared to DBF_SOS sat .The deviation from the respective mean date (period 2003-2010) at individual locations (grid cell) ranged in average between ´14 and 10, and ´12 and 8 days for ENF_SOS sat and ENF_SOS mod throughout Finland.Regarding DBF_SOS, inter-annual variations at one location were slightly smaller than for ENF_SOS, ranging in mean from ´8 to 9, and ´8 to 7 days for DBF_SOS sat and DBF_SOS mod , respectively.The inter-annual variations of mean air temperature for Finland in May correlated with the inter-annual variations in DBF_SOS mod ; that is, higher May temperatures lead to an earlier start of season, (R 2 = 0.71).Instead, for the evergreen forest, April temperatures are more important: the mean air temperature in April was correlated with the ENF_SOS mod time series (R 2 = 0.56).The inter-annual differences between DBF_SOS mod and ENF_SOS mod can be largely explained by the differences in the mean air temperature in April and May (R 2 = 0.70).
We observed an early anomaly in the start of season in ENF in 2007, both from satellite observations and JSBACH simulations (Figures 10 and 11c,f).The yearly anomaly for ENF_SOS sat was larger than ´20 days in south-western coastal areas and in central Finland (Figure 11c).The ENF_SOS mod anomaly agreed in many parts of Finland with the satellite-observed anomaly, but not for western areas at Bothnian Bay and areas in the northern boreal zone (Figure 11f).
ENF_SOSmod throughout Finland.Regarding DBF_SOS, inter-annual variations at one location were slightly smaller than for ENF_SOS, ranging in mean from −8 to 9, and −8 to 7 days for DBF_SOSsat and DBF_SOSmod, respectively.The inter-annual variations of mean air temperature for Finland in May correlated with the inter-annual variations in DBF_SOSmod; that is, higher May temperatures lead to an earlier start of season, (R 2 = 0.71).Instead, for the evergreen forest, April temperatures are more important: the mean air temperature in April was correlated with the ENF_SOSmod time series (R 2 = 0.56).The inter-annual differences between DBF_SOSmod and ENF_SOSmod can be largely explained by the differences in the mean air temperature in April and May (R 2 = 0.70).We observed an early anomaly in the start of season in ENF in 2007, both from satellite observations and JSBACH simulations (Figures 10 and 11c,f).The yearly anomaly for ENF_SOSsat was larger than −20 days in south-western coastal areas and in central Finland (Figure 11c).The ENF_SOSmod anomaly agreed in many parts of Finland with the satellite-observed anomaly, but not for western areas at Bothnian Bay and areas in the northern boreal zone (Figure 11f).Furthermore, agreement was found for parts of central and northern Finland regarding the early anomaly of DBF_SOS in 2006 (Figure 12b,e); although DBF_SOSmod indicated a stronger early anomaly than shown by DBF_SOSsat.Some disagreement between the anomalies was observed in 2007; an early anomaly (up to 15 days) occurred in the south-eastern parts of Finland, according to DBF_SOSsat, but this was not the case for DBF_SOSmod (Figure 12c,f).Furthermore, agreement was found for parts of central and northern Finland regarding the early anomaly of DBF_SOS in 2006 (Figure 12b,e); although DBF_SOS mod indicated a stronger early anomaly than shown by DBF_SOS sat .Some disagreement between the anomalies was observed in 2007; an early anomaly (up to 15 days) occurred in the south-eastern parts of Finland, according to DBF_SOS sat , but this was not the case for DBF_SOS mod (Figure 12c,f).

Quality of Remote Sensing Observations of the Start of Season
The accuracy of DBF_SOSsat was about one week for Finnish sites, which is in line with results obtained in Siberia [70].Moreover, our observations are in agreement with a map of satellite-derived greening-up of vegetation in Fennoscandia for the period 1982-2002 by Karlsen et al. [77]; Figure2a, but DBF_SOSsat occurred slightly later in the northern boreal zone than in the map for a relatively warmer period (2000)(2001)(2002)(2003)(2004)(2005)(2006) [28]; Figure 5.The sites Aulanko, Kevo, Saariselkä and Värriö were also used for evaluation of the NDVI-based method in studies by Karlsen et al. [28,78] that showed higher correlation between field observations and the NDVI-based start of season in Aulanko and Kevo, but larger RMSE and bias for all four sites than in our results (Table 2).
For the evaluation of DBF_SOSsat, we compared visual observation of bud break of birch on a few trees with the satellite signal on green-up for many pixels covering an area of several kilometres.As described by Badeck et al. [67], several factors could cause deviations between the in situ observation and DBF_SOSsat: (a) the inter-individual variability in bud break, (b) systematic time lags between bud break and the date derived from the satellite signal and (c) the heterogeneity in

Quality of Remote Sensing Observations of the Start of Season
The accuracy of DBF_SOS sat was about one week for Finnish sites, which is in line with results obtained in Siberia [70].Moreover, our observations are in agreement with a map of satellite-derived greening-up of vegetation in Fennoscandia for the period 1982-2002 by Karlsen et al. [77] Figure 2a, but DBF_SOS sat occurred slightly later in the northern boreal zone than in the map for a relatively warmer period (2000)(2001)(2002)(2003)(2004)(2005)(2006) [28] Figure 5.The sites Aulanko, Kevo, Saariselkä and Värriö were also used for evaluation of the NDVI-based method in studies by Karlsen et al. [28,78] that showed higher correlation between field observations and the NDVI-based start of season in Aulanko and Kevo, but larger RMSE and bias for all four sites than in our results (Table 2).
For the evaluation of DBF_SOS sat , we compared visual observation of bud break of birch on a few trees with the satellite signal on green-up for many pixels covering an area of several kilometres.As described by Badeck et al. [67], several factors could cause deviations between the in situ observation and DBF_SOS sat : (a) the inter-individual variability in bud break, (b) systematic time lags between bud break and the date derived from the satellite signal and (c) the heterogeneity in terrain and land cover.Concerning (a), no information on the observer error was given with the ground observations, but the visual observations at Finnish sites were carried out twice a week and we therefore assumed an observers error of about three days.The within-population variation in the date of bud break can range between a few days to up to four weeks and varies from year to year depending on weather patterns [79].The variation of microclimate for areas covering several kilometres and elevation differences lower than 50 m was estimated with 2.5 days in Badeck et al. [67].Regarding (b), satellite-observed start of season of DBF did not show a systematic time lag to bud break observations in this study for the whole of Finland.According to Delbart et al. [27,70], the timing of the NDWI minimum in spring after snowmelt corresponds to the time of greening-up of vegetation.Because of noise in the satellite time series, due to cloud cover and atmospheric influences, a threshold value was applied to avoid false detection of the start date (see Equation ( 3)).We minimized the effect of land cover heterogeneity (c) by selecting MODIS pixels with a dominance of deciduous and mixed forest cover.However, differences in forest structure and the contribution of the understorey to the satellite signal may be one reason for the different biases in DBF_SOS sat when looking at the southern and northern boreal sites separately (Section 3.1).More investigations are needed on a possible inherent regional bias in the remote sensing estimate for the start of season date in deciduous broad-leaf forest.Information on the forest structure and the understorey phenology, e.g., from web-camera observations, could be helpful in this respect.
In this study, we use the timing of the beginning of snow melt as a proxy indicator for the start of season in ENF that showed good correspondence with ENF_SOS EC for flux measurement sites in Finland.The quality of the estimate has been discussed in Böttcher et al. [33].Recent studies suggest that the photochemical reflectance index (PRI) can be used as a direct indicator of spring photosynthetic activation in boreal ENF [25] and can potentially be applied to satellite observations, e g., MODIS [80].However, satellite-based monitoring with the PRI will need to take into account, among other things, the varying background [25], and this is especially relevant for the spring recovery in evergreen boreal forest when the ground is still snow-covered.Alternatively, the Plant Phenological Index (PPI) seems insensitive to noise and snow and follows well the seasonal development of GPP in boreal ENF [81].In addition, new space-borne instruments, such as those for the Global Ozone Monitoring Experiment (GOME-2), the Greenhouse Gases Observing Satellite (GOSAT) and the Orbiting Carbon Observatory-2 (OCO-2), provide the opportunity of measuring solar-induced fluorescence and, therefore, ways to directly assess photosynthetic activity [26,82,83], and could be considered in future work.
The presence of long and frequent periods of cloud cover increases the uncertainty of the phenological dates extracted from optical satellite observations.In this study, we applied temporal interpolation of the daily FSC and NDWI time series, but excluded cases of continuous cloud cover longer than two weeks.Including the spatial dimension in the gap-filling based on geostatistical techniques [84][85][86] could improve the completeness of the time series.The quality of the gap-filled time series may still depend on the pattern of missing data in time and space [85].Another possibility would be the use of the all-weather capabilities of satellite radar and microwave observations alone [87] or in combination with optical observations.For ENF, our method relies on observations of snow cover.Microwave retrievals of the Snow Water Equivalent (SWE), such as provided by the GlobSnow project [88], could be used in combination with the optical observations to gain information on the snow during periods of cloud cover.Currently, the SWE time series covering the period from 1979-2012 is only available at a grid resolution of 25 km ˆ25 km, but a higher resolution SWE product (augmented using the optical FSC data) is under development at the FMI.

Modelling of Springtime Development and the Start of Season in JSBACH
The inspection of the modelled GPP time series (Figures 4 and 5) showed that the JSBACH model instantaneously responds to the increase of air temperature above zero ˝C resulting in a higher GPP than obtained from the CO 2 flux measurements and, consequently, earlier ENF_SOS mod than ENF_SOS EC .The photosynthesis module in JSBACH has the same temperature response for parameters describing the potential electron transport rate and maximum carboxylation rate for the whole year, thus indicating that the simulated vegetation is ready to immediately photosynthesize in spring with increasing air temperatures without any recovery period required, providing there is foliage present as in the case of ENF.However, the boreal coniferous forests have a recovery period in spring after the winter dormancy before they reach their full summertime photosynthetic capacity [89,90].This may explain the differences between observed GPP at Sodankylä and simulated GPP using local weather observations (Figures 4 and 5).Overestimation of the spring carbon uptake (a too large GPP) by JSBACH was also found by Dalmonech and Zaehle [21] in a study that evaluated, among others, the phenology of the JSBACH model.Furthermore, because of generally lower air temperatures from the regional climate model than from the site measurements, the ENF_SOS mod,loc was earlier than ENF_SOS mod .
The observed early bias of three days in ENF_SOS mod from JSBACH is relatively small in comparison to results by Richardson et al. [14] for the simulated start of photosynthetic uptake from 14 biosphere models (mean bias of ´11 ˘15 days) and is in the same range as for the best-performing models in the evaluation.The increase in the negative bias of ENF_SOS mod from south to north (Figure 8, Supplementary Materials: Table S2) found here, agrees with higher temperatures observed at time of ENF_SOS at northern sites compared to southern boreal sites [4,18,33] and suggests that the response of air temperature for ENF follows a climate gradient in the boreal zone.
For the performance of the JSBACH model regarding the estimation of ENF_SOS, the photosynthesis model plays the dominant role and deviations from ENF_SOS sat are not due to the LoGro-P model.Improvements to the model performance during the start of season of ENF could be achieved by implementing seasonality to temperature response of JSBACH photosynthesis parameters, and possibly also accounting for photoinhibition in coniferous forests, a phenomenon observed during spring at high light levels and freezing temperatures [91]: the light receiving apparatus in chloroplast gets damaged at high light levels, and for low temperature conditions the photosynthetic capacity may be notably reduced.The description of effects of soil water freeze on carbon uptake could also be improved.According to observations [18,92], the coniferous trees are able to photosynthesize before the soil thaw, but lack of liquid water from soil soon reduces the carbon uptake.In the model, the actual simulated proportions of liquid and frozen water in soil layers could be used to limit the carbon exchange.There are advances related to expressing frozen soil in JSBACH northern permafrost regions [93].
The modelled FSC by JSBACH does not serve as a reliable predictor of ENF_SOS (Section 3.3), as too little snow and too early snow melt is simulated for our study area.Having bias-corrected temperature and precipitation data does not guarantee precise snow cover, as the formulation of snow accumulation relies on a universal temperature threshold and parameterizations of snow processes are not adjusted for region and drivers.
In contrast to results for ENF, we found a late bias of DBF_SOS mod (Figure 7, Supplementary Materials: Table S2).This agrees with findings by Dalmonech and Zaehle for southern Finland [21], but due to the higher temporal resolution of the remote sensing data set used here, we can quantify the bias better and also give results for the northern parts of Finland that were masked out in the study above.As regards other biosphere models from the study by Richardson et al. [14], the bias for JSBACH DBF_SOS mod turned out to be relatively small, although direct comparisons are difficult due to differences in the environmental drivers for the model runs.
It has to be noted that the definition of DBF_SOS mod does not correspond exactly to DBF_SOS sat .While DBF_SOS mod is defined based on a minimal change in LAI, the visual observation of bud break may indicate a later stage of development.An optimized parameterization of the LoGro-P model for the boreal zone would allow an improvement of the overall bias for DBF_SOS mod ; variation in parameter values, such as the alternating temperature (T alt ) and the minimum value of the critical heat sum (S min crit ) can reduce the mean bias (Appendix A, Tables A1 and A2).A proper tuning of the LoGro-P parameters for our area of interest was outside the scope of this study and needs further investigation, but could be performed using satellite observations and Bayesian optimization techniques.
Similar to ENF, the bias in DBF_SOS mod showed a regional gradient; that is, the bias decreased from south to north (Figure 7, Supplementary Materials: Table S2).According to field observations, the heat sum requirement for bud break of birch decreases from south to north in the boreal region [94,95] and also the temperature threshold for the NDVI-based greening-up decreases [28,77].This seems to be in contradiction to our results for JSBACH DBF_SOS mod .The same critical temperature and heat sum requirement was used in the LoGro-P model for the whole study area.According to the above, one would expect a gradient with earlier bias in southern areas compared to northern areas.We observed the opposite, except for the northernmost parts of the study area (Figure 8b).Reasons for this could be an insufficient chilling predicted by the LoGro-P model, a regional gradient in bias in the bias-corrected driving temperature, or in the remote sensing estimate of DBF_SOS sat .Under the current climate, the chilling requirement is typically fulfilled in the boreal region [96] and this was also the case in our study area.Therefore, the chilling requirement cannot be used for a better calibration of the south-north gradient, and since the start date of the chilling decay follows a climate gradient, an increase in the decay parameter (C decay ) would further pronounce the delay in the DBF_SOS mod in southern regions.The S min crit parameter plays a role as well, but while the overall bias could be reduced with a lower parameter value, the regional gradient would still remain and this may suggest that optimal parameters for the boreal region are zone-specific.
Despite a very small temperature bias in the start date of the "growth phase" of DBF for the whole of Finland (Figure 8c), the slope of the regression line between DBF_SOS sat and the modelled date from LoGro-P was slightly higher (0.9 compared to 0.79) when FMI gridded temperatures were used instead of temperature data from the regional climate model.This indicates a small influence of the driving data to the regional gradient.Likewise, the late bias in DBF_SOS sat for phenological sites in northern, and early bias in southern Finland (Section 3.1) need further investigations and may partly explain differences in bias in JSBACH DBF_SOS mod in the different boreal sub-zones.

Conclusions
In this study, we evaluated the timing of the start of the photosynthetically active season in evergreen needle-and deciduous broad-leaf forest in Finland in the biosphere model JSBACH against satellite observations.The differences between the two forest types with respect to photosynthetic recovery in spring were taken into account.The accuracy of the satellite-derived dates for the deciduous forest was assessed with bud break observations.The NDWI-based method is applicable for deciduous forest in Finland and the estimated dates can be used for the evaluation of model results at a regional scale.The inaccuracy of the satellite-observed date was one week and the bias was small (Figure 2), but regional differences in bias need further evaluation.
The comparisons between JSBACH-simulated and satellite-observed start of season at the regional scale showed that the general spatial and temporal dynamics were rather well captured by JSBACH, for both ENF and DBF (Figures 7 and 10).The performance was slightly better for evergreen forest.We found an early bias of three days and a late bias of five days in the start of season of evergreen needle-and deciduous broad-leaf forest, respectively.The bias showed a south-north gradient for both forest types.While in the deciduous forest the start of the vegetation active season is determined by the phenological model LoGro-P, the photosynthesis model is more important for the evergreen forest type.The date of the predicted "growth phase" in evergreen forest occurred about one month after the start of the vegetation active season.Therefore, an improved calibration of the phenology model would not be sufficient for this plant functional type.Instead, JSBACH photosynthesis parameters would need to be adjusted to limit too strong early photosynthesis and hence to improve the behaviour of gross primary production in early spring (Figures 4 and 5).In the deciduous broad-leaf case, results could be slightly improved with different JSBACH phenology parameters, but the model formulation seems not optimal with respect to the south-north gradient in Finland.In further work, one could optimize the phenological parameters with a Bayesian method using the remote sensing observations to better understand the effectiveness of individual parameters and the relationships between them.Furthermore, using in situ measurements such as chlorophyll fluorescence observations and greenness change from phenology cameras could be used to improve the photosynthesis and phenology descriptions of the model.The logistic dynamics (Equation (A1)) is such that whenever Λ fell to zero it could never recover.Plants solve this problem by growing seeds or leaf buds containing biochemical resources sufficient to grow a leaf until they can grow further by carrying out photosynthesis.In LoGro-P, this is modelled by assuming a minimum value of Λ called "seed value"-when conditions are favourable for starting growth.
Environmental conditions enter the dynamics only via the growth and shedding rates k and p.Hence, the more complicated part of LoGro-P is the modelling of plant's dependence on temperature, soil moisture etc., which is different for the various phenology types, as explained below.

A3. Evergreen Phenology
Basic to the phenology model for evergreens is the assumption that the phenological dynamics can be divided into two subsequent phases, a "growth phase" (spring), characterized by non-zero k and absence of leaf shedding (p " 0); and a "rest phase" (all other seasons), where growth is zero (k " 0) and leaf shedding p is non-zero.If not zero anyway, growth and shedding rates are prescribed parameters.

A4. Summer-Green Phenology
The summer-green phenology assumes three phases: "growth phase" (spring), characterized by non-zero k and absence of leaf shedding (p " 0); a "vegetative phase" (summer), where growth is zero (k " 0) and leaf shedding p is small; and a "rest phase" (all other seasons), where growth is zero (k " 0) and leaf shedding p is non-zero.
Vegetative phase and rest phase differ only by the magnitude of the shed rate.The introduction of three phases instead of only those two of the evergreens is suggested by fAPAR data of summer-green vegetation, where typically, after a short time of rapid growth, a reduction of fAPAR is seen, much before the leaves are shed in autumn-also, here, if not anyway zero, growth and shedding rates are prescribed parameters.

A5. Spring and Autumn Events
A main issue in modelling the phenological phases is to determine when the vegetation shifts from one phase to the next.For the transition from the growth to the subsequent phase it is assumed that the growth phase lasts only for a fixed time, so that after that time the next phase automatically starts.However, for the shift from the rest phase to the growth phase (called spring event) and the shift from the vegetative phase to the rest phase (autumn event), the dependence on environmental conditions cannot be neglected.This is the topic of the following sections.

A6. The Spring Event
To determine the spring event, the alternating model of Murray et al. [50] is used (see also Botta et al. [97]): Let S pdq denote the value of the heat sum at day d: S pdq " where T(d) is the mean day temperature at day d, T alt is a parameter called alternating temperature (which has the function of a cut-off temperature in the heat sum), and d 0 is the starting date for temperature summation.In the model, it is taken care that this starting date falls in times after the buds have already experienced chilling (see below).
Another key quantity of the alternating model is the number of chill days, C pdq: this is the number of days with a mean day temperature below the alternating temperature T alt , where counting starts at the day d a of the last autumn event: C pdq " is fulfilled-if at the spring event Λ is less than the above mentioned seed value, Λ is set to this seed value to allow for sufficiently fast growth.

A7. The Autumn Event
At the autumn event the vegetative phase of summer-greens switches to the rest phase with rapid leaf shedding.The autumn event is calculated from a smoothed air temperature T s pnq at every time step n (see below).The smoothing is performed to prevent "false" autumn events upon short periods of coldness.The autumn event happens when during the vegetative phase T s falls below a critical temperature T crit s .To prevent that this event is detected in spring, the condition is added that the mean day air temperature T pdq should be smaller than T s , i.e., the autumn event happens at the first day d when T pdq ă T s pnq ă T crit s (A8) Heat summation starts in the northern hemisphere always on 1 January (day 1), and in the southern hemisphere on day 183.

Figure 1 .
Figure 1.(a) CO2 flux measurement sites, phenological sites and phytogeographical zones in Finland according to Ahti et al. (1968) (b) Red-Green-Blue (RGB) image of forest cover fractions determined from CORINE land cover 2000.

Figure 1 .
Figure 1.(a) CO 2 flux measurement sites, phenological sites and phytogeographical zones in Finland according to Ahti et al. (1968) (b) Red-Green-Blue (RGB) image of forest cover fractions determined from CORINE land cover 2000.
In addition, due to long periods of cloud cover, the DBF_SOS sat observations were not available for the following site-years: Aulanko 2001, Kolari 2005, Saariselkä 2001 and Kevo 2007 and 2008.

Figure 2 .
Figure 2. Relationship of the start of season in deciduous broad-leaf forests (DBF_SOSsat) with phenological observations of birch bud break for the period from 2001 to 2008: (a) for all site-years and (b) for the year-to-year anomalies.

Figure 2 .
Figure 2. Relationship of the start of season in deciduous broad-leaf forests (DBF_SOS sat ) with phenological observations of birch bud break for the period from 2001 to 2008: (a) for all site-years and (b) for the year-to-year anomalies.

Figure 3 .
Figure 3.Comparison of the start of season in evergreen needle-leaf forest determined from eddy covariance measurements (ENF_SOSEC) and satellite observations (ENF_SOSsat), together with JSBACH model-derived dates for 2001 to 2010 in (a) Hyytiälä, southern Finland and (b) Sodankylä, northern Finland.The start of season from JSBACH was determined from model runs using both REMO-generated bias-corrected meteorological data (ENF_SOSmod) and local site observations (ENF_SOSmod,loc).

Figure 3 .
Figure 3.Comparison of the start of season in evergreen needle-leaf forest determined from eddy covariance measurements (ENF_SOS EC ) and satellite observations (ENF_SOS sat ), together with JSBACH model-derived dates for 2001 to 2010 in (a) Hyytiälä, southern Finland and (b) Sodankylä, northern Finland.The start of season from JSBACH was determined from model runs using both REMO-generated bias-corrected meteorological data (ENF_SOS mod ) and local site observations (ENF_SOS mod,loc ).

Figure 4 .
Figure 4. (a-c) Daily GPP (Gross Primary Production), mean and minimum air temperature and snow depth for Sodankylä in spring 2004.Vertical lines indicate the start of season (ENF_SOS): JSBACH local (magenta) refers to ENF_SOSmod,loc; JSBACH regional (cyan) is ENF_SOSmod, and observation refers to ENF_SOSEC (black).

Figure 4 .
Figure 4. (a-c) Daily GPP (Gross Primary Production), mean and minimum air temperature and snow depth for Sodankylä in spring 2004.Vertical lines indicate the start of season (ENF_SOS): JSBACH local (magenta) refers to ENF_SOS mod,loc ; JSBACH regional (cyan) is ENF_SOS mod , and observation refers to ENF_SOS EC (black).

Figure 5 .
Figure 5. (a-c) Daily GPP, mean and minimum air temperature and snow depth for Sodankylä in spring 2008.Vertical lines indicate the start of season date: JSBACH local (magenta) refers to ENF_SOSmod,loc; JSBACH regional (cyan) is ENF_SOSmod.Note that ENF_SOSEC occurred on the same day as ENF_SOSmod shown with cyan line.

Figure 5 .
Figure 5. (a-c) Daily GPP, mean and minimum air temperature and snow depth for Sodankylä in spring 2008.Vertical lines indicate the start of season date: JSBACH local (magenta) refers to ENF_SOS mod,loc ; JSBACH regional (cyan) is ENF_SOS mod .Note that ENF_SOS EC occurred on the same day as ENF_SOS mod shown with cyan line.

Figure 6 .
Figure 6.Spatial distribution of mean start of season dates in Finland for the period 2003-2010.The upper panel shows estimates for evergreen needle-leaf forest: (a) determined from satellite observation, and (b) determined from JSBACH simulations.The lower panel shows estimates for deciduous broad-leaf forest: (c) determined from satellite observations and (d) from JSBACH simulations.

Figure 6 .
Figure 6.Spatial distribution of mean start of season dates in Finland for the period 2003-2010.The upper panel shows estimates for evergreen needle-leaf forest: (a) determined from satellite observation, and (b) determined from JSBACH simulations.The lower panel shows estimates for deciduous broad-leaf forest: (c) determined from satellite observations and (d) from JSBACH simulations.

Figure 7 .
Figure 7. Scatterplots of mean satellite-and JSBACH model-derived start of season for forested areas in Finland for the period 2003-2010 for: (a) evergreen needle-leaf forest (ENF_SOS); and (b) deciduous broad-leaf forest (DBF_SOS).The dotted line gives the 1:1 line and the black line the linear fit between JSBACH-modelled and satellite-observed dates.

Figure 7 .Figure 8 .Figure 8 .
Figure 7. Scatterplots of mean satellite-and JSBACH model-derived start of season for forested areas in Finland for the period 2003-2010 for: (a) evergreen needle-leaf forest (ENF_SOS); and (b) deciduous broad-leaf forest (DBF_SOS).The dotted line gives the 1:1 line and the black line the linear fit between JSBACH-modelled and satellite-observed dates.Remote Sens. 2016, 8, 580 15 of 31

Figure 8 .Figure 9 .
Figure 8. Grid cell-wise bias between start of season derived from the JSBACH model and remote sensing in Finland for the period 2003-2010: (a) in evergreen needle-leaf forest; and (b) in deciduous broad-leaf forest; (c) Bias in JSBACH-modelled starting date of the "growth phase" of DBF because of the driving temperature data.The bias was calculated as the difference in the start date of the "growth phase" simulated by running the LoGro-P-model with regional climate model data and with gridded temperature observations.Red colours show areas with a late bias and blue areas indicate an early bias of the modelled date.

Figure 9 .
Figure 9. (a-f) Scatterplots of satellite-observed and JSBACH-simulated start of season for individual years 2005-2007.Upper panel: evergreen needle-leaf forest (ENF_SOS); Lower panel: deciduous broad-leaf forest (DBF_SOS).The dotted line gives the 1:1 line and the black line the linear fit between satellite-and JSBACH-modelled estimates.

Figure 10 .
Figure 10.Interannual variations of the start of season in evergreen coniferous forest and deciduous broad-leaf forest in Finland derived from the JSBACH model and from remote sensing observations.

Figure 10 . 31 Figure 11 .
Figure 10.Interannual variations of the start of season in evergreen coniferous forest and deciduous broad-leaf forest in Finland derived from the JSBACH model and from remote sensing observations.Remote Sens. 2016, 8, 580 17 of 31

Figure 11 .
Figure 11.(a-f) Spatial distribution of yearly anomalies, the deviations (in days) from grid cell mean for 2003-2010, of start of season in evergreen needle-leaf forest.Upper panel: remote sensing estimate (ENF_SOS sat ).Lower panel: model-derived estimate (ENF_SOS mod ).

Figure 12 .
Figure 12. (a-f) Spatial distribution of yearly anomalies, the deviations (in days) from grid cell mean 2003-2010, of the start of season of deciduous broad-leaf forest.Upper panel: remote sensing estimate (DBF_SOSsat).Lower panel: model-derived estimate (DBF_SOSmod).

Figure 12 .
Figure 12. (a-f) Spatial distribution of yearly anomalies, the deviations (in days) from grid cell mean 2003-2010, of the start of season of deciduous broad-leaf forest.Upper panel: remote sensing estimate (DBF_SOS sat ).Lower panel: model-derived estimate (DBF_SOS mod ).

d ÿ d 1
"d a θ pT alt ´T pdqq (A5)where θ is the Heaviside step function.To prevent that under bad climate conditions C pdq influences the phenological behaviour of the next growing season, the number of chill days is limited to 365 days.From C pdq a critical heatsum S crit pdq is computed:S crit pdq " and C decay are parameters of the alternating model.Finally, the spring event happens at the first day d for which during the rest phase S pdq ě S crit pdq (A7)

Table 1 .
Phenological observation sites for the evaluation of DBF_SOS sat and characteristics of MODIS pixels at 0.05 ˝by 0.05 ˝resolution that were used in the comparison.

Table 2 .
Comparison between DBF_SOS sat and phenological observations of bud break of birch for sites in Finland.A positive bias means that the detected SOS is late compared to observations.

Table 2 .
Comparison between DBF_SOSsat and phenological observations of bud break of birch for sites in Finland.A positive bias means that the detected SOS is late compared to observations.

Table 3 .
Start of season in evergreen needle-leaf forest for the period 2001-2010 at Hyytiälä and Sodankylä as defined from CO 2 flux measurements (ENF_SOS EC ) and based on JSBACH model simulations with local (ENF_SOS mod,loc ) and regional driving data (ENF_SOS mod ).A negative bias means that the simulated date is early compared to observations.ENF_SOS EC was not available for 2002 at Sodankylä, therefore only 9 years were used in the analysis. 1

Table 3 .
Start of season in evergreen needle-leaf forest for the period 2001-2010 at Hyytiälä and Sodankylä as defined from CO2 flux measurements (ENF_SOSEC) and based on JSBACH model simulations with local (ENF_SOSmod,loc) and regional driving data (ENF_SOSmod).A negative bias means that the simulated date is early compared to observations.
the vegetation active season of DBF DBF_SOS mod Starting date of the vegetation active season of DBF determined from JSBACH regional runs DBF_SOS sat Starting date of the vegetation active season of DBF determined from satellite observations doy the vegetation active season in ENF ENF_SOS EC Starting date of the vegetation active season determined from continuous EC measurements ENF_SOS mod,loc Starting date of the vegetation active season of ENF determined from JSBACH local runs at CO 2 flux measurement sites ENF_SOS mod Starting date of the vegetation active season of ENF determined from JSBACH regional runs ENF_SOS sat Starting date of the vegetation active season of ENF determined from satellite observationsFor small Λ, this equation describes exponential growth at rate k reduced by leaf shedding at rate p.The term in brackets reduces growth when Λ approaches the maximum LAI Λ max of the particular PFT.Equation (A1) has the solution Φ pt `τq " pk ´pq Φ ptq kΦ ptq `pk ´p ´kΦ ptqq exp p´pk ´pq τq is called the "phenological state".To simulate the development of the LAI Equation (A2) is iterated at a daily time step τ.