Assessment of Spatio-Temporal Patterns of Black Spruce Bud Phenology across Quebec Based on MODIS-NDVI Time Series and Field Observations

Satellite remote sensing is a widely accessible tool to investigate the spatiotemporal variations in the bud phenology of evergreen species, which show limited seasonal changes in canopy greenness. However, there is a need for precise and compatible data to compare remote sensing time series with field observations. In this study, fortnightly MODIS-NDVI was fitted using double-logistic functions and calibrated using ordinal logit models with the sequential phases of bud phenology collected during 2015, 2017 and 2018 in a black spruce stand. Bud break and bud set were spatialized for the period 2009–2018 across 5000 stands in Quebec, Canada. The first phase of bud break and the last phase of bud set were observed in the field in mid-May and at the beginning of September, when NDVI was 80.5% and 92.2% of its maximum amplitude, respectively. The NDVI rate of change was estimated at 0.07 in spring and 0.04 in autumn. When spatialized on the black spruce stands, bud break was detected earlier in the southwestern regions (April–May), and later in the northeastern regions (mid to end of June). No clear trend was observed for bud set, with different patterns being detected among the years. Overall, the process bud break and bud set lasted 51 and 87 days, respectively. Our results demonstrate the potential of satellite remote sensing for providing reliable timings of bud phenological events using calibrated NDVI time series on wide regions that are remote or with limited access.


Introduction
Plant phenology is the study of recurring life cycle events, such as growth reactivation and dormancy, leaf emergence and senescence, and flowering.Phenology is considered as a sensitive bio-indicator of climate and its changes [1][2][3].Over the past three decades, phenology studies have been conducted using field-based observations, bioclimatic models, near-surface remote sensing and satellite remote sensing based techniques [4][5][6][7].Researchers have applied different approaches depending on the question addressed and the spatial scale of investigation.Satellite data provide wide coverage with varying temporal, spectral and spatial resolutions.Remote sensing is therefore a flexible, reliable and widely recognized tool to study phenology in a number of ecosystems, from forests to grasslands, which may be too difficult using other data collection methods [8][9][10].Current sensors on board satellite platforms record spectral signals that can be used to monitor seasonal and interannual variations in vegetation cover and determine the timings of phenological events and the growing season across the landscape [11,12].
Imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS) on board the Terra and Aqua satellites is widely used to study vegetation phenology [4,13].Among the available tools, the Normalized Difference Vegetation Index (NDVI) [14,15] is a reliable spectral index for the reconstruction of phenological transitions of different vegetation types, including croplands, forests and grasslands, because it is related to the amount of green-leaf biomass [13,16,17].NDVI is considered as a potential remote sensing proxy to investigate the effect of climate on vegetation phenology at regional to continental scales due to its close relationship with plant activity [18,19].
In the last two decades, substantial advances have been made in predicting length and start/end dates of the growing season from both field observations and satellite remote sensing data [20][21][22].However, the scarcity of ground observations in some locations is a challenge for studies of growth patterns over large or remote regions of the world [23,24].Another constraint exists in stands of evergreen species, which show only small seasonal variations in canopy greenness compared to deciduous species [21,25,26].Phenology and growing season in North America have been studied extensively using the NDVI at medium (500 m) and coarse (1 and 8 km) spatial resolutions [22,27,28].However, direct observations in the field confirming the results from remote sensing in remote areas remain scarce.
To address these issues, we propose an innovative statistical approach to calibrate MODIS NDVI time series and demonstrate how it can be used to describe annual changes in black spruce [Picea mariana (Mill.)bud phenology across a large portion of the coniferous forest in a study area in Quebec, Canada.This approach provides a robust alternative to field-based monitoring, especially in areas where field observations are difficult due to the remoteness or extreme weather conditions.

Study Area and Field Observations Site
The overall study area is located in Quebec, Canada, and covers nearly 800,000 km 2 between 45.5 • -55 • N and 64 • -80 • W. Using the Quebec 1:20k forest map from the Quebec Government [29], we selected 5000 polygons representing stands dominated by black spruce trees (>75%) and that had not been disturbed over the past 30 years (Figure 1a).
Field measurements in this study were taken in a black spruce stand of the Forêt d'Enseignement et de Recherche Simoncouche (48 • 22 N, 71 • 25 W, 338 m a.s.l.) in Quebec, Canada.The study site is located within the balsam fir (Abies balsamea L. Mill.)-white birch (Betula papyrifera Marsh.)bioclimatic domain, at the southern border of the boreal forest (Figure 1b).The site has a typical boreal climate characterized by short cool summers and cold harsh winters.Mean the annual temperature is 0.9 • C, with absolute minimum and maximum temperatures of −36.7 • C in January-February, and 31.1 • C in July, respectively.Mean annual precipitation is 1162 mm [30].

Identification of the Phenological Phases
We collected observations of bud phenology on 105 black spruce saplings each week from May to October during 2015, 2017 and 2018.We collected data on saplings due to the need to perform reliable observations directly on the buds, which would have been difficult to perform on tall mature trees.All phenological phases of bud break and bud set were recorded according to [31].Six phases were defined for bud break: open bud (B1), with a pale spot at the tip of the bud; elongated bud (B2), with lengthening brown scales; swollen bud (B3), with smooth and pale-coloured scales but no visible needles; translucent bud (B4), with needles visible through the scales; split bud (B5), with open scales but needles still clustered; and exposed shoot (B6), with needles completely emerged from the surrounding scales and spreading outwards.Five phases of bud set were defined: white bud (S1), with the presence of a white bud; beige bud (S2), with beige scales around the bud; brownish bud (S3), with a significant increase in volume; brown bud (S4), with needles starting to spread outwards; and spread needles (S5), with the needles in the whorl spreading outwards.(S3), with a significant increase in volume; brown bud (S4), with needles starting to spread outwards; and spread needles (S5), with the needles in the whorl spreading outwards.

MODIS Data
We acquired time series of 250 m Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1, version 6) from the NASA Earthdata website (https://earthdata.nasa.gov).MOD13Q1 is derived from atmospherically corrected bi-directional surface reflectance imagery and contains vegetation index data as well as the pixel reliability layer needed for quality checking [32].We used 16 day NDVI maximum composite values in ISIN tiles h13v03 and h13v04 from 2009 to 2018, thus covering a large portion of the eastern Quebec boreal and temperate forests.In total, this represented 460 NDVI granules (230 for each ISIN tile for the whole period).All the images were filtered to exclude unreliable NDVI pixels using a pixel reliability layer supplied by MODIS vegetation index production team.MODIS scenes were re-projected to NAD83/Quebec Lambert projection and stacked by band (i.e., vegetation index) and by date of the composite period for each tile.

Assessing and Calibrating NDVI
We selected 40 polygons with an area:perimeter ratio greater than 13 for pure black spruce stands adjacent or close to the Simoncouche study site [29].The forest polygons selected for all analyses included stands composed by >75% of black spruce.We used these polygons to derive standlevel NDVI time series by calculating, for each NDVI image, the mean of all pixels located within the

MODIS Data
We acquired time series of 250 m Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1, version 6) from the NASA Earthdata website (https://earthdata.nasa.gov).MOD13Q1 is derived from atmospherically corrected bi-directional surface reflectance imagery and contains vegetation index data as well as the pixel reliability layer needed for quality checking [32].We used 16 day NDVI maximum composite values in ISIN tiles h13v03 and h13v04 from 2009 to 2018, thus covering a large portion of the eastern Quebec boreal and temperate forests.In total, this represented 460 NDVI granules (230 for each ISIN tile for the whole period).All the images were filtered to exclude unreliable NDVI pixels using a pixel reliability layer supplied by MODIS vegetation index production team.MODIS scenes were re-projected to NAD83/Quebec Lambert projection and stacked by band (i.e., vegetation index) and by date of the composite period for each tile.

Assessing and Calibrating NDVI
We selected 40 polygons with an area:perimeter ratio greater than 13 for pure black spruce stands adjacent or close to the Simoncouche study site [29].The forest polygons selected for all analyses included stands composed by >75% of black spruce.We used these polygons to derive stand-level NDVI time series by calculating, for each NDVI image, the mean of all pixels located within the limits of the polygons.NDVI phenology curves were then fitted for each year using a double-logistic model, where time (t) was represented by Day Of the Year (DOY): Three pairs of coefficients were included: minimum and maximum NDVI (min and max), timings of the inflection points (S and A), and rates of change at the inflection points (mS and mA) [20,33].The estimated double-logistic functions were standardized (std NDVI ) to the range 0-1 and compared to field observations of phenological phases.
We applied logit models for bud break (six phases, i.e., B1-B6) and bud set (five phases, i.e., S1-S5) to estimate the probability of observing each sequential phenological phase along the std NDVI measurements.This procedure fits the proportional odds model to ordinal response data, here represented by the sequential phases of bud phenology.The ordinal logit models were applied using the LOGISTIC procedure in SAS [34].We assessed the performance of the model three times, each time using two years of bud phenology observations for model training and the remaining year for validation.Observations and predictions were compared by using linear regression and the goodness of fit (R 2 ) and the root mean square error (RMSE) in prediction as performance metrics [35].

Spatializing Black Spruce Phenology
We again extracted NDVI values but this time for the period from 2009 to 2018 and for 5000 pure black spruce sites [29] distributed over two MODIS tiles covering our study area.Polygons with a ratio area:perimeter > 13 were selected and for each polygon, NDVI pixels were extracted and used to calculate mean polygon NDVI for each MODIS image.The area:perimeter threshold value was used to ensure that we used mostly clumps of pixels and excluded long and narrow polygons, this decreasing the risk of including mixed pixels in our dataset.We applied the thresholds estimated by the above-mentioned models to estimate the beginning (B1) and ending (S5) of bud phenology.We performed a hotspot and coldspot analysis using Getis-Ord Gi* [36] in ArcGIS 10.2 (Environmental Systems Research Institute, Inc., Redlands, CA) to identify spatial clustering in phenological events [37].We used the inverse distance weighted interpolation method to interpolate the resulting values over the entire area [38].

MODIS NDVI Data
Annually, NDVI showed a bell-shaped pattern with a slow increase in spring, followed by a rapid increase and culmination in July and a decrease in autumn until reaching a minimum value in winter (Figure 2).On average, NDVI ranged from 0.32 in winter and spring to 0.9 in summer.The double-logistic function represented NDVI well during the whole season, including the asymmetry in the transitions between spring, summer and autumn.The satisfactory fitting was confirmed by the analysis of the residuals, which were uniformly distributed around zero for the three studied years.Only 3% of the residuals exceeded the range between −0.2 and 0.2, confirming the model reliability (Figure 2).
Compared to 2017 and 2018, the year 2015 showed the highest NDVI (0.917) in summer (Table 1).In 2015 and 2018, the inflection points (coefficients S and A) occurred at approximately the same time, at the beginning of May (DOY 121 and 122) and the end of October (DOY 284 and 290), respectively.In 2017, NDVI started to increase earlier (DOY 118) and decreased later (DOY 311) compared to the other two years, which resulted in a 20% increase of the growing season (DOY 192) compared to 2017 (DOY 162) and 2018 (DOY 167).On average, coefficients mS and mA were estimated at 0.07 and 0.04, respectively, demonstrating that autumnal rates of reduction in NDVI were slower than increases in the spring (Table 1).
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 17 at 0.07 and 0.04, respectively, demonstrating that autumnal rates of reduction in NDVI were slower than increases in the spring (Table 1).

Comparing NDVI with Bud Phenology
Bud break (B1-B6) and part of bud set (S1-S2) occurred under increasing (Figure 3).varied from 0.81 to 0.97 during bud break (i.e., B1-B6), culminating at phase S2.Phases S3 to S5 occurred under decreasing , at values ranging between 0.92 and 1. Larger standard deviations were observed for bud break phases than for bud set, except for phase S5.The probability of occurrence of each bud phenological phase according to was estimated with ordinal logit models.Model validation produced satisfying estimations of the timings of the phenological phases (Figure 4).The results of the model were significant, as indicated by the high R 2 (0.87) and a RMSE of 12.85.The regression estimated a slope of 0.97, which was not statistically different from one (p > 0.05).
With < 0.70, the probability of observing a dormant bud (phase B0) was 0.94, which reduced at increasing (Figure 5).The probability of observing B1 exceeded dormant bud at of 0.75, meaning that at this value, the buds were more likely to be active than dormant.B1 culminated at of 0.79, reaching the highest probability of occurring.B2 exceeded B1 at of 0.84, and culminated with of 0.87.B3, B4 and B5 culminated at values ranging between 0.91 and 0.96.Above 0.96, corresponded with the highest probability of observing B6, which represented the end of the bud burst process (Figure 5).S1-S3 showed the highest probability of being observed at values between 0.97 and 1.In autumn, S4 culminated at of 0.95.At < 0.93, S5 exceeded S4, indicating that the buds

Comparing NDVI with Bud Phenology
Bud break (B1-B6) and part of bud set (S1-S2) occurred under increasing std NDVI (Figure 3).std NDVI varied from 0.81 to 0.97 during bud break (i.e., B1-B6), culminating at phase S2.Phases S3 to S5 occurred under decreasing std NDVI , at values ranging between 0.92 and 1. Larger standard deviations were observed for bud break phases than for bud set, except for phase S5.The probability of occurrence of each bud phenological phase according to std NDVI was estimated with ordinal logit models.Model validation produced satisfying estimations of the timings of the phenological phases (Figure 4).The results of the model were significant, as indicated by the high R 2 (0.87) and a RMSE of 12.85.The regression estimated a slope of 0.97, which was not statistically different from one (p > 0.05).
With std NDVI < 0.70, the probability of observing a dormant bud (phase B0) was 0.94, which reduced at increasing std NDVI (Figure 5).The probability of observing B1 exceeded dormant bud at std NDVI of 0.75, meaning that at this std NDVI value, the buds were more likely to be active than dormant.B1 culminated at std NDVI of 0.79, reaching the highest probability of occurring.B2 exceeded B1 at std NDVI of 0.84, and culminated with std NDVI of 0.87.B3, B4 and B5 culminated at std NDVI values ranging between 0.91 and 0.96.Above 0.96, std NDVI corresponded with the highest probability of observing B6, which represented the end of the bud burst process (Figure 5).S1-S3 showed the highest probability of being observed at std NDVI values between 0.97 and 1.In autumn, S4 culminated at std NDVI of 0.95.At std NDVI < 0.93, S5 exceeded S4, indicating that the buds were more likely to be dormant.Consequently, the bud set process was considered complete when std NDVI was below 0.9 (Figure 5).
were more likely to be dormant.Consequently, the bud set process was considered complete when was below 0.9 (Figure 5).were more likely to be dormant.Consequently, the bud set process was considered complete when was below 0.9 (Figure 5).were more likely to be dormant.Consequently, the bud set process was considered complete when was below 0.9 (Figure 5).

Spatial Pattern of Black Spruce Phenology
The dates of bud break (B1) were not homogenous across the study area (Figure 6).It generally occurred earlier in the southwestern region (end of April to May) and later in the northeastern region (mid to end of June), with the larger variations observed along the longitudinal gradient.Overall, the earliest and latest bud break dates were observed in 2010 (DOY 165) and 2012 (DOY 180), respectively.The highest spatial heterogeneity in bud break dates across the study area was observed in 2012, as shown by the large variances in Table 2. Bud break was more homogenous in 2013.The earliest bud break occurred in the last week of April (DOY 115-120) during 2010, 2012 and 2013.On average, bud set (S5) started at the beginning of August (DOY 215) and ended at the end of October (DOY 303) (Figure 7).The spatial pattern of bud set was not constant and no clear trend could be found among years.Late bud set occurred more frequently in western regions, mainly in 2009 and 2014.When comparing the years, the bud set date was more homogeneous across the study area in 2017 and more heterogeneous in 2009 (Table 2).The earliest bud set were estimated in 2013, 2016, and 2018.The range between the first and last bud break was 51 days, which is almost half the range observed in bud set (87 days).Therefore, the spatial pattern of the growing season length seems to be dominated by bud set.The longest growing seasons were observed in 2012 and 2014, while the shortest ones were observed in 2013 and 2016.The southwestern and central-eastern regions showed the longest growing seasons, while the shortest were detected in the north-central region of the study area.

Discussion
Time series of satellite-derived vegetation indices are a reliable tool to describe the patterns of bud break and bud set dates across large and remote regions.In this study, we used field observations of all black spruce phenology phases to calibrate time series of MODIS NDVI, which we used to spatialize bud break and bud set dates across the black spruce stands in Quebec, Canada.The first phase of bud break (B1), occurring around mid-May, corresponded to a standardized NDVI of 80.5%.NDVI culminated at the end of July when the process of bud set was occurring.Winter bud (i.e., the last phase of bud set) was reached when standardized NDVI reached 92.2% of its maximum amplitude at the beginning of September.Our study provided a new statistical approach based on ordinal logit modelling to statistically connect the sequential events of phenological phases with fortnightly NDVI data.
Remotely sensed phenology of the boreal forest has been performed across North America or in combination with the entire northern hemisphere [9,19,22,39].Wide area coverage provided by satellite remote sensing offers a unique point of view at synoptic scale but lacks in fine spatial details.To our knowledge, species-specific calibrations of the major parameters of boreal forest bud phenology based on satellite imagery are still missing.Fu et al. [40] calibrated MODIS vegetation phenology at global scales (i.e., Global Vegetation Phenology) using only remote sensing phenological observations from 2001 to 2010, observing the summer culmination in NDVI and late beginning of the growing season (in mid-April to mid-May) in the boreal and cool regions of North American evergreen coniferous forests.White et al. [27] calibrated NDVI at 8-km resolution using ground phenology records from 1982 to 2006 for North American forests but they explicitly removed evergreen species from their analysis.Our results improve the phenological description by calibrating MODIS NDVI time series on field observations of bud phenology at higher spatial resolution (250 m), enhancing the reliability of the estimations at the stand level for two important phenological phases, the beginning and end of bud activity.
Other studies have used satellite-based retrievals of solar-induced fluorescence (SIF) to study phenological changes in evergreen canopies [41,42].However, apart from the spatial resolution of the current satellites used for SIF (from ~3 km 2 (OCO-2) to >1500 km 2 (GOME-2)), the main difference between the approach presented here and SIF-based phenology resides is their respective definitions of phenological events.For example, SIF, which is mostly driven by plant fluorescence yield and by the amount of incoming photosynthetically active radiation (PAR), is an indicator of plant carbon uptake (gross primary production, photosynthetic light-use efficiency).Consequently, SIF-and other phenological indicators based on plant physiology [43,44] define, for example, the start of the growing season (SOS) as the beginning of photosynthetic activity (e.g., date of spring recovery in evergreen trees).In contrast, with our approach, SOS is defined as the bud break date which, from a spectral perspective, will result in visible changes in canopy color, namely greenness.Seasonally, photosynthetic activity and changes in canopy greenness and structure are not well synchronized, either in the spring [45][46][47] or in the fall and as a result, this may lead to differences in estimates of the growing season length between physiology-based and structure-based indicators of plant phenology [42].
Compared with autumn, we observed higher slopes of the double-logistic function in spring.The estimated winter NDVI (coefficient min of the function) was able to reduce the influence of snowmelt during winter, thus representing an effective greening up of evergreen trees during spring [20,48].In addition, bud break in black spruce starts at the end of May when snowmelt is completed [30].These results are in accordance with a previous study that used NDVI time series from satellite remote sensing to investigate xylem growth and timing of plant phenology across the boreal forests of Quebec [49].In this study, we demonstrated the existence of a relationship between black spruce bud phenology phases at a wide geographical scale and a remote sensing-derived vegetation index (NDVI).Our results show that establishing a link between bud phenology observations and NDVI time series provides a comprehensive view of boreal forest activity patterns and trends in remote areas where direct observations or recurrent samplings are unachievable.
Methods to detect the beginning and end of the growing season from NDVI are already available in the literature [13,50].In optical remote sensing, the variability in spectral signatures is considered as a noise disturbance.In general, these noises are related to varying atmospheric conditions, the presence of snow on coniferous canopies and the spatial and spectral variability of understorey vegetation [51][52][53].Phenological events may therefore be harder to estimate for boreal evergreen coniferous forests, which also exhibit lower seasonal changes in foliage biomass and thus, in vegetation indices such as the NDVI [4,20].
In this study, the double-logistic model was able to suitably describe the seasonal pattern of NDVI, despite the wide spectral variability within and between forest stands.It is well known that such a function can effectively reduce the noise in remote sensing time series, thus adequately describing the timings and duration of the growing season [53,54].In previous studies, the inflection points of double logistic models fitted to NDVI data more accurately estimated spring than autumnal events in boreal species [20,22].This was attributed to the longer and slower decrease in canopy greenness occurring in autumn [4,55,56].In this study, we were able to adequately describe and detect the inflection point in autumn using the double-logistic function that represents the ending phase of bud set.Over the last two decades, NDVI has been the most commonly used vegetation index for remote sensing-based modelling of phenology [13,17,27,56].Our method has great potential to be applied to a wider range of species worldwide, after species-specific calibration.By providing a novel and reliable statistical approach, this study could allow the improvement of the performance of phenological models for evergreen forests to the same levels as has been achieved for grasslands and deciduous forests [57,58].

Relationship between NDVI and Field-Based Phenology
Our study assesses the explicit link between satellite imagery and phenological events recorded in the field [13,59].We observed a quick increase in NDVI during the bud break process from mid-May to the end of June and a slow decrease in NDVI during the bud set phases from July to October.The punctual phases of bud break were identified along the continuous process of gradual variation in greenness.The variability in NDVI at bud set was lower than that at bud break because NDVI was close to its maximum when bud set was already occurring.By comparing bud phenological phases with NDVI curves, the peak standardized NDVI occurred at S2.It was observed four weeks and one after the completion of bud burst and the start of bud set, respectively.This increase in NDVI could be related to changes in needle biomass resulting from a gradual increase during bud burst phases and shoot elongation during spring and summer.Furthermore, leaf expansion, i.e., increment in leaf length and width during spring, is relatively rapid due to the transition from dormancy to active photosynthesis [60][61][62].
The decrease in NDVI during late summer and autumn may be associated with needle aging and declining pigmentation [63].Moreover, in evergreen species, some of the old needles, commonly those with ages varying from five to seven years, change color at the end of summer and fall in autumn.Given that MODIS vegetation indices (NDVI, EVI) are primarily sensitive to changes in leaf chlorophyll content and structure [18,19,64], NDVI may be more representative of phenological changes in evergreen species during the first part of the growing season (from mid-spring to summer) than in the autumn [4,22].This could be attributed to the reduction in photosynthetic efficiency and the lower carbohydrate demand during dormancy [45], which results in changes in light absorbance and reflectance.For evergreen trees, detection of these changes is difficult as they usually require specific and very narrow wavebands, still lacking in most current orbiting sensors [48,65].However, in this study, the calibration of NDVI time series with phenological observations and the resulting threshold values alleviates some of the difficulties encountered for evergreen species when using other approaches.This calibration is specific to black spruce stands because different thresholds of NDVI could be estimated for other tree species of the boreal forest, mainly for deciduous species.We expect this model to be reliable for closed stands strongly dominated by black spruce.Overall, our results suggest that satellite NDVI may be used as a large-scale complement to field observations of bud phenology, and therefore, support the general applicability of satellite-based vegetation indices to estimate both the beginning and end of the growing season in boreal evergreen forests, as often defined by bud break and bud set dates.

Spatial and Temporal Changes of Bud Phenological Phases
The beginning of bud break (phase B1) appeared later in the northeastern regions, which are also closer to sea, resulting in a shorter duration of the overall bud break process compared with bud set.On average, bud break lasted 51 days across the study area, which covered approximately 5 • in latitude.This trend confirms results from other studies; however, we extended our analysis to 5000 sites compared to a previous study that analyzed altitudinal and latitudinal gradients of xylem phenology at only six sites black spruce sites [49].The timings of bud break were observed from the beginning of May to the end of June and are related to a number of factors, including the fulfilment of the requirement in winter chilling, the lengthening of photoperiod and warming in spring temperature [66,67].In contrast, bud set lasted longer, 87 days, and lacked a clear and constant spatial pattern.Bud set is an important phenological event and explains most of the variation in tree growth [68,69], while the bud set in black spruce is a typical photoperiod-dependent process [70].In conifers, the cessation of shoot elongation and development of terminal buds indicates vegetative maturity.It therefore seems to result from exposure to the shortening days of late summer.During this time, the temperature remains warm so it could be a response to some endogenous signals triggered by the shortened photoperiod [66,71].Overall, satellite remote sensing offers the possibility of modeling phenology by recording spectral information at regular time intervals due to their wider coverage area and high temporal resolution [22,72].Our model utilized the potential of time series satellite data to provide the spatio-temporal patterns of bud phenology by calibrating NDVI across the Quebec region of Canada from 2009 to 2018.We expect that the presented calibration approach could be tested on a wider basis for other sites and tree species.

Conclusions
With the availability of limited field-based observations due to extreme weather, remote locations and tedious sampling, validation with the help of ground observations data is a key issue in satellite remote sensing-based phenology analyses over large areas [73].Uncertainties increase with the magnitude in scale mismatch between collected field data and coarse resolution remote sensing-based observations [13].Calibration of bud phenological events using NDVI time series on the one hand indicates the limitations of this approach for spatio-temporal pattern analyses, but on the other hand, also demonstrates the ability of remote sensing to upscale bud phenology over larger forested areas.The most confounding issues are related to the selection of vegetation index suitable for the estimation of phenological metrics during the start and end of the growing season in evergreen forests, which exhibit smaller seasonal variations in canopy optical properties compared to deciduous forests [22,74].In this study, we defined a new methodology that aims to retrieve the bud phenological phases for black spruce using MODIS-NDVI time series.We aggregated NDVI time series at stand level and used a double-logistic curve to filter and describe NDVI-time series over time.We also used a unique set of field-based observations to model bud phenology using ordinal logit models and proposed statistical-based thresholds to identify the bud phenological phases for a boreal coniferous species in Quebec, Canada.These thresholds were used to generate a spatial and temporal portrait of bud phenology on a yearly basis over large areas.We calibrated NDVI for black spruce, but different threshold values are expected for other tree species.Thus, our study can be reliably applied to mostly pure black spruce stands, while applications to other species should be considered carefully.Overall, our results suggest that satellite-derived NDVI can be used as a proxy for monitoring phenological events when suitable statistical methods are applied to calibrate remote sensing data with field observations.

Figure 1 .
Figure 1.Location of the 5000 black spruce stands in Quebec, Canada (a) selected from two ISIN tiles (h13 v03 and h13 v04, not shown) of MODIS, and the stand where phenological phases were measured (b).The zoomed-in area shows the 40 black spruce polygons (in red color) used for calibration within the buffer of 25 km 2 around the field observation site.

Figure 1 .
Figure 1.Location of the 5000 black spruce stands in Quebec, Canada (a) selected from two ISIN tiles (h13 v03 and h13 v04, not shown) of MODIS, and the stand where phenological phases were measured (b).The zoomed-in area shows the 40 black spruce polygons (in red color) used for calibration within the buffer of 25 km 2 around the field observation site.

Figure 2 .
Figure 2. NDVI time series (grey circles) and fitted double-logistic curves (solid black line) for three years in a black spruce boreal stand of Quebec, Canada.Each day of the year shows the mean NDVI values for 40 black spruce stands used for model calibration.

Figure 2 .
Figure 2. NDVI time series (grey circles) and fitted double-logistic curves (solid black line) for three years in a black spruce boreal stand of Quebec, Canada.Each day of the year shows the mean NDVI values for 40 black spruce stands used for model calibration.

Figure 3 .
Figure 3. Correspondence between phenological phases measured in the field and mean standardized MODIS NDVI values of black spruce stands used for model calibration (n = 40).Dots and error bars represent mean and standard deviation, respectively, for the three years used for calibration.

Figure 4 .
Figure 4. Relationship between observed DOY and simulated DOY in three years.The black line is the regression line between the observed DOY and the simulated DOY.The small panel shows the distribution of the residuals.

Figure 5 .
Figure 5. Probability of occurrence of phenological phases derived from field observations and standardized MODIS NDVI values of boreal black spruce stands in Quebec, Canada.

Figure 3 .
Figure 3. Correspondence between phenological phases measured in the field and mean standardized MODIS NDVI values of black spruce stands used for model calibration (n = 40).Dots and error bars represent mean and standard deviation, respectively, for the three years used for calibration.

Figure 3 .
Figure 3. Correspondence between phenological phases measured in the field and mean standardized MODIS NDVI values of black spruce stands used for model calibration (n = 40).Dots and error bars represent mean and standard deviation, respectively, for the three years used for calibration.

Figure 4 .
Figure 4. Relationship between observed DOY and simulated DOY in three years.The black line is the regression line between the observed DOY and the simulated DOY.The small panel shows the distribution of the residuals.

Figure 5 .
Figure 5. Probability of occurrence of phenological phases derived from field observations and standardized MODIS NDVI values of boreal black spruce stands in Quebec, Canada.

Figure 4 .
Figure 4. Relationship between observed DOY and simulated DOY in three years.The black line is the regression line between the observed DOY and the simulated DOY.The small panel shows the distribution of the residuals.

Figure 3 .
Figure 3. Correspondence between phenological phases measured in the field and mean standardized MODIS NDVI values of black spruce stands used for model calibration (n = 40).Dots and error bars represent mean and standard deviation, respectively, for the three years used for calibration.

Figure 4 .
Figure 4. Relationship between observed DOY and simulated DOY in three years.The black line is the regression line between the observed DOY and the simulated DOY.The small panel shows the distribution of the residuals.

Figure 5 .
Figure 5. Probability of occurrence of phenological phases derived from field observations and standardized MODIS NDVI values of boreal black spruce stands in Quebec, Canada.

Figure 5 .
Figure 5. Probability of occurrence of phenological phases derived from field observations and standardized MODIS NDVI values of boreal black spruce stands in Quebec, Canada.

Figure 6 .
Figure 6.Hotspot analysis of bud break (phase B1) dates (DOY) across boreal black spruce stands in the study area in Quebec, Canada.

Figure 6 .
Figure 6.Hotspot analysis of bud break (phase B1) dates (DOY) across boreal black spruce stands in the study area in Quebec, Canada.

Figure 7 .
Figure 7. Hotspot analysis of bud set (phase S5) dates (DOY) across boreal black spruce stands in the study area in Quebec, Canada.

Figure 7 .
Figure 7. Hotspot analysis of bud set (phase S5) dates (DOY) across boreal black spruce stands in the study area in Quebec, Canada.

Table 1 .
Coefficients of the double-logistic function fitted to three years of MODIS NDVI values extracted for 40 pure black spruce boreal stands in Quebec, Canada.min and max represent minimum and maximum NDVI, S and A represent timings of the inflection points and rates of change at the inflection points represented by mS and mA.

Table 1 .
Coefficients of the double-logistic function fitted to three years of MODIS NDVI values extracted for 40 pure black spruce boreal stands in Quebec, Canada.min and max represent minimum and maximum NDVI, S and A represent timings of the inflection points and rates of change at the inflection points represented by mS and mA.

Table 2 .
Variability in bud phenology across the study area.