Proxy Indicators for Mapping the End of the Vegetation Active Period in Boreal Forests Inferred from Satellite-Observed Soil Freeze and ERA-Interim Reanalysis Air Temperature

Triggered by decreases in photoperiod and temperature, evergreen needle-leaved trees in the boreal region downregulate photosynthetic activity and enter dormancy in autumn. Accompanying changes in canopy structure and chlorophyll content are small and precede the cessation of photosynthetic activity. Low solar elevation and cloud cover during this period pose additional challenges for the use of optical satellite instruments. Alternatively, environmental variables that correlate with photosynthesis, such as soil freeze, can be detected from satellite microwave observations independent of weather and illumination conditions. We tested for the first time the usability of satellite-observed soil freeze from the Soil Moisture and Ocean Salinity (SMOS) instrument as a proxy indicator for the end of vegetation active period (VAPend) at six eddy covariance sites in Finland and Canada. The time when soil freeze commenced over the large SMOS pixel can be employed to estimate VAPend (R2 = 0.84, RMSE = 7.5 days), defined as the time when the photosynthetic capacity of the forest drops below 10% of the growing season maximum. In comparison to satellite-based soil freeze timing, an air temperature-based proxy from ERA-Interim reanalysis data showed better performance (R2 = 0.92, RMSE = 5.2 days). VAPend was mapped in the boreal forest zone in Finland and Canada from both indicators based on linear regression models.


Introduction
Ecosystem processes related to carbon, water and nutrient cycling are directly linked to vegetation phenology. Climate change is affecting the timing of phenological events, and vegetation phenology has, therefore, been proposed as an indicator for long-term biological impacts of climate change on terrestrial ecosystems (Schwartz 1998;Menzel and Fabian 1999;Richardson et al. 2013). Lengthening of the growing season is expected as the consequence of increasing temperatures caused by climate change (Fronzek et al. 2012;Laapas et al. 2012). An earlier start of the growing season in the boreal forest increases carbon sequestration (Black et al. 2000;Pulliainen et al. 2017). Warmer autumns resulted in a greater increase in respiration compared to gains in photosynthesis (Piao et al. 2008;Vesala et al. 2010). Still, the impacts of predicted warming with rising temperatures in autumn and the foreseen lowered temperature seasonality (Xu et al. 2013) are uncertain and photoperiod may constrain potential future lengthening of the growing season in autumn (Mäkelä et al. 2006). Despite the importance of autumn phenology for the variability of productivity, it has received less consideration than spring phenology Gonsamo and Chen 2016).
Trends in autumn phenological events, such as leaf colouring, were analysed from in situ observations of deciduous species and from satellite-observed phenology based on vegetation indices. Results vary regionally depending on the considered time period. Satellite observations show delayed leaf senescence in the boreal forest in Eurasia, parts of Fennoscan-dia and in Northern America for the period 1982-2008 (Jeong et al. 2011). A later end of the growing season was observed in southern Fennoscandia for the period 1982-2006 (Karlsen et al. 2009). Instead, from ground observations of leaf colouring of birch in Finland, no significant changes were found for the period 1997-2006 (Pudas et al. 2008). Autumn Normalized Difference Vegetation Index (NDVI) increased during the period 1982-2006 in large parts of Eurasia (Piao et al. 2011). However, trends in autumn senescence from the Phenological Index were not significant for most parts of the boreal region for the period 1999-2013 (Gonsamo and Chen 2016).
To protect their photosynthetic apparatus against frost damage during winter, coniferous trees in the boreal zone develop frost hardiness and enter dormancy. The decrease in photosynthetic capacity and the cessation of growth are triggered by shortening of the photoperiod and decreasing temperature (Vogg et al. 1998;Öquist and Huner 2003;Suni et al. 2003;Sevanto et al. 2006;Vesala et al. 2010;Bauerle et al. 2012). The oldest leaves of evergreen needle-leaved trees are shed in autumn with the lifespan and turnover rates depending on the tree species. The lifespan increases and turnover rates decrease from south to north in the boreal evergreen forests (Reich et al. 2014). After the oldest foliage drops, evergreen needle-leaved trees continue to photosynthesize until the occurrence of severe frosts (Vogg et al. 1998;Suni et al. 2003;Wu et al. 2012). The development of frost hardiness leads to a loss of chlorophyll and an increasing capacity to transform absorbed light into heat (Öquist and Huner 2003). The process of downregulation of photosyn-thesis in autumn can be divided into two stages. First, an initial decline is observed as a response to low temperatures. During this stage, photosynthesis recovers fast when conditions become favourable. Second, a more persistent stage of dormancy is reached with harsher winter conditions and photosynthesis recovers more slowly even when conditions become favourable (Kolari et al. 2014). We refer to the latter stage as the end of the vegetation active period (VAPend), where the vegetation active period (VAP) is defined as the period of the year when vegetation is assimilating carbon through photosynthesis.
Satellite methods that track changes in the leaf area index (Jin and Eklundh 2014), the pigment content (Gamon et al. 2016;Ulsig et al. 2017), and the end of the growing season (Gonsamo and Chen 2016) in boreal evergreen needle-leaved forest from vegetation indices were developed, but the detection of VAPend from satellite observations remains challenging. Observations with optical instruments are hindered by frequent cloud cover and low solar elevation in autumn at high latitudes. Recently, the re-emission of excess energy, not used in the photosynthetic process, has been observed as sun-induced chlorophyll fluorescence (SIF) from space (Joiner et al. 2014;Walther et al. 2016). Satellite retrievals of surface fluorescence were found to be less affected by atmospheric scattering than those of surface reflectance (Frankenberg et al. 2012;Frankenberg et al. 2014;Köhler et al. 2015); however, at low solar elevation, the noise in the SIF signal increases and is possibly biased by rotational-Raman scattering (Joiner et al. 2014). Moreover, currently available satellite measurements of chlorophyll fluorescence have a large footprint which complicates comparison with CO 2 flux measurements (Joiner et al. 2014).
Therefore, the development of robust satellite indicators for VAPend in boreal evergreen needle-leaved forest is still lacking. In addition to observations of changes in canopy structure, chlorophyll content and, more directly, chlorophyll fluorescence, environmental variables that correlate with photosynthesis, such as temperature, soil freeze (Jarvis and Linder 2000;Barr et al. 2009;Barichivich et al. 2013), snow cover, and incoming radiation could help to track autumn seasonal changes in boreal forests from space. In this work, we will focus on satellite indicators of soil freeze.
Despite the fact that soil freeze is not a primary cause for the start of photosynthetic downregulation in autumn, we hypothesize that it can be used as a proxy indicator for VAPend in evergreen boreal trees. Photosynthesis is always associated with water loss due to transpiration that needs either be balanced by soil water uptake or change in stem water storage (Sevanto et al. 2006). As soil freeze will block the water transport from the soil to trees, it may be more directly related to photosynthesis (Hölttä et al. 2017), even though photosynthesis has been observed when the soil and tree stem were partially frozen in boreal coniferous forest (Suni et al. 2003;Sevanto et al. 2006).
Soil freeze can be detected from active and passive microwave observations (Kimball et al. 2004;Smith et al. 2004;Bartsch et al. 2007;Kim et al. 2011;Kim et al. 2012;Roy et al. 2015;Rautiainen et al. 2016;Derksen et al. 2017) with the advantage of independence from weather conditions and low solar elevation during autumn and winter in the high latitudes. New observations of the soil freeze state from L-band (1-2 GHz) active radar and passive microwave radiometers (Roy et al. 2015;Rautiainen et al. 2016;Derksen et al. 2017) bring the advantage of higher emission depth and lower sensitivity to vegetation and snow cover compared to observations from higher frequencies (Kimball et al. 2004;Kim et al. 2011) that have been used earlier to derive information on soil freeze (Rautiainen et al. 2012;Roy et al. 2015).
The objective of this study was to investigate the possibility of applying satellite-observed soil freeze as a proxy indicator for VAPend in boreal evergreen needle-leaved forests. For this, we utilized observations of soil freeze state detected from L-Band passive microwave observations from the European Space Agency's Soil Moisture and Ocean Salinity (SMOS) mission ). In addition to satellite-derived soil freeze, we also included an autumnal freeze indicator from ERA-Interim reanalysis air temperature data in our analysis. We described the late stages of photosynthetic downregulation with three dates for VAPend (see Sect. 2.4) that were determined from CO 2 flux measurements for six sites in the boreal forest zone in Finland and Canada. The relationship between VAPend and indicators for air and soil freeze was analysed, and based on that, proxy indicators for mapping of VAPend at continental to global scales were developed.

Measurement Sites
We selected six CO 2 flux measurement sites in evergreen needle-leaved boreal forest. Three sites are located in Finland and three sites in Canada (Fig. 1, Table 1). The Canadian sites (CA-Obs, CA-Ojp and CA-Qfo), and one Finnish site (Hyytiälä, FI-Hyy) belong to the southern boreal phytogeographical zone, whereas two Finnish sites, Kenttärova (FI-Ken) and Sodankylä (FI-Sod), belong to the northern boreal zone. All Finnish sites are located further north (62-68°N) than the Canadian sites (50-54°N). The solar elevation at solar noon decreases below 30°at the end of October at CA-Qfo but at the beginning of September at FI-Ken. The dominant tree and understory species vary by site.
CA-Obs and CA-Ojp are located at the southern edge of the boreal forest about 100 km north-east of Prince Albert in central Saskatchewan, Canada. The distance between the two sites is approximately 30 km. Both sites are characterized by old growth forests that regenerated after fire events in 1879 and 1929 at CA-Obs and CA-Ojp, respectively (Kljun et al. 2006). CA-Obs is dominated by black spruce (Picea mariana Mill. BSP) with sparsely distributed tamarack (Larix laricina Du Roi) and jack pine (Pinus banksiana Lamb), and the ground cover consists of mosses (sphagnum and feather mosses) and sparse understory. Further site characteristics are given in Kljun et al. (2006), Bergeron et al. (2007) and Krishnan et al. (2008). CA-Ojp is a homogenous jack pine forest with a uniform age structure (Griffis et al. 2003;Kljun et al. 2006). The understory is dominated by lichen and green alder (Alnuscrispa Ait.). Due to its location in a topographic depression, CA-Obs is poorly drained and has a high water table. Soil at CA-Ojp is coarse textured and well drained (Kljun et al. 2006). CA-Qfo, located 30 km south of Chibougamau (Quebec), is dominated by black spruce and the forest floor consists mainly of feather mosses, sphagnum and lichen. The forest stand covers the tower flux footprint, except in the south-east direction where a peatland is situated. Further details of the site are given in Bergeron et al. (2007) and Coursolle et al. (2012).
FI-Hyy in southern Finland belongs to the SMEAR II station (Station for Measuring Ecosystem-Atmosphere Relations) of the University of Helsinki. The stand is dominated by Scots pine (Pinus sylvestris L.) sown in 1962 with some naturally regenerated Norway spruce (Picea abies) and scattered deciduous trees. It is a typically managed forest (Vesala et al. 2010). The forest floor is composed of dwarf shrubs and mosses (Kulmala et al. 2008).
Sites FI-Ken and FI-Sod are located above the Arctic Circle in Finnish Lapland. The Scots pine forest in FI-Sod is naturally regenerated and is representative of pine-dominated forested areas in central Lapland (Thum et al. 2009). The site belongs to the Arctic Research Centre of the Finnish Meteorological Institute. FI-Ken, a homogenous Norway spruce (Picea abies [L.] H. Karst.) forest, is situated in the Pallas-Yllästunturi National Park (Aurela et al. 2015).
According to the GlobCover map (Arino et al. 2012), open (15-40%) needle-leaved deciduous or evergreen forests is the dominant land cover type for the surrounding areas of all CO 2 flux measurement sites corresponding to the area covered by the SMOS pixels. The coverage differs by site showing the highest proportions for CA-Ojp (84%) and lowest for CA-Qfo (42%) (Online Resource Table S1). The second highest contributions come from closed to open (> 15%) mixed broadleaved and needle-leaved forest for FI-Sod and FI-Hyy. For CA-Obs and CA-Qfo, the proportion of open water bodies within the SMOS pixel is relatively high (14% and 24%, respectively).

CO 2 Flux Measurements
The in situ fluxes of CO 2 were measured by the micrometeorological eddy covariance method which provides direct measurements of CO 2 exchange between atmosphere and biosphere averaged on an ecosystem scale. In the eddy covariance method, the vertical CO 2 flux is obtained as the covariance of the high-frequency (10 Hz) observations of vertical wind speed and the CO 2 concentration (Baldocchi 2003). The measurements were made 5-10 m above the mean forest height.
The eddy covariance fluxes were calculated as half-hourly averages taking into account the appropriate corrections following Aubinet et al. (2012). The measurement systems and the post-processing procedures have been presented in more detail for FI-Hyy by Kolari et al. (2004), Vesala et al. (2005) and Mammarella et al. (2009), for FI-Sod and FI-Ken by Aurela et al. (2009) and Aurela et al. (2015), for CA-Obs and CA-Ojp by Barr et al. (2004), and for CA-Qfo by Bergeron et al. 2007. The AmeriFlux Base data products (version 1-1) were used for the Canadian sites (Black 1999a, b;Margolis 2003).

Soil Temperature Measurements at Sodankylä
To compare satellite-based soil freeze with in situ observations, we utilized measurements from 8 measurement sites of the automatic network of in situ soil moisture and soil temperature observations at the Arctic Research Centre at Sodankylä (Ikonen et al. 2016). The sites cover typical soil types and land cover classes in the area and allow better representation of soil freeze conditions for the large areas covered by the SMOS footprint compared to the local data at the CO 2 flux sites. Sites are located in 170 m to 23 km distance to FI-Sod. At each soil measurement site, the vertical profile of soil properties (dielectric constant, electric conductivity and soil temperature) is measured at different depths with 5TE (Decagon Devices Inc., Pullman, WA, USA) or CS655 (Campbell Scientific Inc., Logan, UT, USA) soil moisture sensors. Two additional sensors are located 10 m away from each vertical profile in opposing directions and measure at 10 and 5 cm depths. Further details of the measurement setup are described in Ikonen et al. (2016). In this study, we utilized daily average soil temperature observations at the 5 cm depth for the period 2011-2016. Measurements from the available sensors were averaged for each site. To account for inaccuracy in the temperature measurements (± 1°C for 5TE and ± 0.5°C for CS655) and dependence of the freezing point on soil water content (Kozlowski 2004(Kozlowski , 2016, we defined the freezing point for each site individually during the melting period in spring as the temperature of the temporary plateau that is observed until all ice has melted. The first date when soil temperature dropped below the freezing point and the first date of the period with continuous soil freeze were determined. The number of available sites for each year ranged from 5 to 8.

Indicators of the End of the Vegetation Active Period
VAPend was determined for the 6 CO 2 flux measurement sites (Table 1) from the continuous CO 2 flux data measured by the micrometeorological eddy covariance method. The downregulation of photosynthesis in autumn is a gradual process and it is difficult to calculate definite dates for VAPend from CO 2 flux measurements. Different approaches, based on curve-fitting and absolute or relative threshold values, were utilized to describe the seasonal cycle of photosynthesis of plant communities from time series of net or gross ecosystem photosynthesis (e.g. Gu et al. 2009;Thum et al. 2009;Richardson et al. 2010). Here, the day on which the CO 2 uptake falls below a selected threshold level is defined as VAPend (Fig. 2). This is similar to the method for the extraction of the start of VAP in Böttcher et al. (2014) and Pulliainen et al. (2017). In practice, the dates were obtained from the annual cycle of the gross photosynthesis index (GPI), which indicates the photosynthetic capacity of the ecosystem on a daily time scale. It was calculated as a difference between daily averages of nighttime (photosynthetic photon flux density (PPFD) < 20 µmol m −2 s −1 ) and day- The obtained GPI value does not strictly represent gross ecosystem photosynthesis since the nighttime respiration differs from the daytime respiration, but it provides a reasonable estimate when used as a relative measure. We tested three different GPI threshold (GPIth) values: 5%, 10% and 15% of the growing season maximum, which is the 90th percentile of the daily GPI values during the month of the highest uptake for  the multi-year measurement period. Different threshold values were included because we do not have prior knowledge as to which one would be most useful indicator for VAPend. As an additional condition, the GPI has to stay below the GPIth for three consecutive days, and the number of following days with GPI > GPIth may not exceed the number of days with GPI < GPIth. The GPI time series was linearly gap-filled in order to estimate VAPend. The length of the gap and the timing of VAPend within the gap are recorded and taken into account when estimating the uncertainties of the data. To obtain information about multi-year mean and interannual variability of VAPend for each site, we calculated the dates of VAPend based on the three threshold (VAPend th ) values for all available observation years (see Table 2).

ERA-Interim Air Temperature Data
ERA-Interim global atmospheric reanalysis air temperature (2 m) data produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Dee et al. 2011) was utilized to derive an indicator on autumn freeze. ERA-Interim reanalysis data are not direct observations. The data set is produced with a forecast model that includes coupled components of the atmosphere, land surface and the ocean waves. Available observational data sets from satellites and conventional observing systems are assimilated into the forecast model to constrain the state of the atmosphere and the Earth surface. Interpolated temperature observations from the synoptic station network over land are combined with the latest state of the atmosphere from the forecast model and utilized to update temperature estimates of the land surface model, thus indirectly influencing the background forecast model (Dee et al. 2011). The accuracy of ERA-Interim surface temperature is higher for areas with a dense station network (Simmons et al. 2010). The air temperature data set in the original spatial resolution of 0.75°× 0.75°was resampled to the Equal Area Scalable Earth (EASE) grid 2.0 (Brodzik et al. 2012) corresponding to the SMOS soil state data. We calculated the 10-day moving average air temperature (T 10 ) and defined the first date of autumnal freeze as the day when T 10 ≤ − 1°C. This threshold corresponds to a considerable reduction of ecosystem carbon uptake in boreal coniferous forest according to Hollinger et al. (1999).

Soil Freeze from Microwave Observations
Soil freeze/thaw state was determined from SMOS satellite data using the method by Rautiainen et al. (2016). We used global SMOS daily gridded level 3 brightness temperature from the Centre Aval de Traitement des Données SMOS (CATDS) (Al Bitar et al. 2017) as input data. The soil freeze/thaw detection algorithm is based on variations in soil permittivity that occur during soil freeze and thaw (permittivity of liquid water~80 and of ice~3 in L-Band) (Mätzler et al. 2006). This change in permittivity leads to higher L-Band brightness temperature of frozen soil compared to thawed moist soil. The algorithm uses a pixel-based relative frost factor calculated from the normalized ratio of brightness temperatures at vertical and horizontal polarizations which is compared against pixel-based thaw and freeze references to determine the soil freeze and thaw state. To reduce noise and the effect of short-term variations on the frost factor, a moving average with a window size of 20 days was applied. Ancillary information on air temperature from the ERA-Interim reanalysis data (see Sect. 2.5) is used to remove false summer freeze and winter thaw detections. Further details of the method are described in Rautiainen et al. (2016). The data set gives daily information on global soil freeze and thaw state in Northern Hemisphere at a spatial resolution of 25 km × 25 km in the EASE grid 2.0 for the period 2010-2016. The soil state is categorized into three classes: frozen, partially frozen and thawed. The SMOS instrument retrieves information from the top-soil layer (down to depths of 10-20 cm, depending on soil type and vegetation). Thus, the frozen soil state can be understood as a soil frost depth of at least 10 cm over the whole pixel. The class partially frozen corresponds to the transition period between first soil frost occurrences to completely frozen soil within one pixel. The soil data set was compared with (1) global atmospheric reanalysis (ERA-Interim) soil temperature data and (2) local, point scale in situ measurements of soil temperature and the volumetric liquid water content ). For comparison with global (1) and local (2) data, the unbiased root mean square error was 20 and 17 days, respectively. The SMOS product indicates systematically later (about 2 weeks) first soil freeze than that using the ERA-Interim Level 2 soil temperature (7-28 cm layer) data set. The bias was smaller for the comparison with in situ data sets (7 days for soil temperature and 1 day for volumetric liquid moisture content), but the correlations were moderate which is probably due to the fact that in situ observations do not represent soil state within the large SMOS pixel adequately . Time series of freeze and thaw state were extracted for the six CO 2 flux sites in Finland and Canada using the ascending (6 am) and descending SMOS orbits (6 pm), respectively. Different orbits were used to avoid signal contamination from man-made radio frequency interference (Soldo et al. 2015). We extracted the first day of partially frozen and frozen soil from time series of the soil freeze and thaw state for the comparison with VAPend. For years with a sudden freeze of the soil, when partially frozen soil was not observed, the first date of partially frozen soil was set to the first date of soil freeze.

Proxy Indicators for Mapping the End of the Vegetation Active Period
We analysed the correspondence of the VAPend th dates with the satellite-derived first date of partially frozen soil, the first date of frozen soil, and the first date of autumnal freeze based on air temperature. To assess the relationship between the freeze indicators and VAPend th , the coefficient of determination (R 2 ), the root mean square error (RMSE) and the bias were calculated. We selected the satellite soil freeze indicator and the VAPend th date which showed the highest correspondence to define a proxy indicator for VAPend based on linear regression. The performance of the proxy indicator from satellite soil freeze was compared with the air temperature-based indicator.
To produce spatial maps of VAPend for the average period 2010-2016, we applied the respective linear regression equations to the gridded soil freeze and air freeze indicator data set covering the boreal vegetation zone in Canada, Finland and surrounding areas. Areas outside the boreal forest zone (see Fig. 1) were masked out.

End of Vegetation Active Period at Measurement Sites
We utilized three dates for VAPend based on threshold values 15%, 10% and 5% of the maximum GPI (Fig. 2) to describe stages of the gradual change in photosynthetic downregulation in boreal evergreen needle-leaved forest (Sect. 2.4). Among all six sites, FI-Hyy in the southern boreal zone showed the latest VAPend dates and the highest interannual variability in the observed dates (

Comparison of the End of the Vegetation Active Period with Autumnal Freeze and Satellite-Derived Soil Freeze
The VAPend dates were significantly correlated with the satellite-derived first date of partially frozen soil and the first date of frozen soil; R 2 ranged between 0.3 and 0.84 (p < 0.05, see Online Resource Table S2, Fig. S1). The highest correspondence was found between the first date of partially frozen soil and VAPend (Fig. 3a). The first date of partially frozen soil explained 84% of the variation in VAPend 10 , and the transition date occurred on average, 15 days before the soil was partially frozen within the SMOS pixel. The first date of frozen soil to a depth of more than 10 cm over the whole SMOS pixel was detected much later than VAPend; it lagged behind VAPend 5 and VAPend 15 by 24 to 41 days, respectively. There was a smaller delay between VAPend 5 and the first date of partially frozen soil (4 d) than for VAPend 10 . The slope of their relationship was close to unity. However, there was higher scatter in the dates (R 2 0.74, Online Resource Fig S1a, Table S2) which could be partly due to higher uncertainties in VAPend 5 . Due to the lower radiation levels during late autumn, there were less data available for calculating the daytime averages. If there were not enough data with PPFD > 600 µmol m −2 s −1 , data at lower radiation Doy is day of year starting from 1 January. When the transition date occurred in the following year the doy is added to the last day of the previous year (31 December) levels were used with a scaling factor. If there were still not enough daytime data, the GPI values were interpolated between nearest observations. For some years, the determination of VAPend 5 relied on interpolated values of the GPI (e.g. Fig. 2). Particularly large discrepancies between VAPend 5 and the date of partially frozen soil were observed in FI-Hyy in 2013 and 2016 when the VAPend 5 lagged behind the first date of partially frozen soil by 20 and 27 days, respectively. In both years, there was a large gap in the GPI time series (1 month in 2013 and 10 days in 2016) leading to high uncertainty in the estimated date. Therefore, we consider VAPend 5 a less reliable indicator than VAPend 10 .
The progress of the GPI and satellite soil freeze is illustrated for 2 years with different temperature development in autumn at FI-Sod (Figs. 4 and 5). In 2013, temperature declined fast and GPI 15 to GPI 5 thresholds occurred within a few days, starting on 18 October and ending on 20 October (doy 291-293, Fig. 4). Due to the absence of snow, soil temperature changes followed closely the air temperature, and soil temperatures at the CO 2 flux site reached below-zero a few days days before GPI 15 was attained. Dry soils (Haplic Arenosols) in open and forested areas froze first. The latest soil freeze date was detected in forested areas with wetter soil type (Umbric Gleysol), and in wetlands with organic soil (Fig. 4). Soil freeze started within the SMOS pixel, i.e. partially frozen soil was detected, on doy 302. The mean soil frost date from soil temperature observation was observed 12 days before, on doy 290 (Online Resource Table S3). Soil freeze for the whole SMOS pixel was detected for the first time on 3 November (doy 307). Except for one wetland site, in situ soil temperature measurements remained below the freezing point after 21 November (doy 325). This coincided well with the time of the continuous period of soil freeze detected from the SMOS pixel.
In contrast to 2013, autumn 2015 was characterized by alternating air temperature values (Fig. 5). A short frost period around doy 283 resulted in a reduction of photosynthetic carbon uptake that exceeded GPI 15 and GPI 10 for some days. Photosynthetic capacity recovered during the following warmer period until the subsequent frost period during which VAPend 15 was determined on 27 October (doy 300). This was again followed by above zero temperature and slightly increased photosynthetic capacity. Similarly, soil temperature measurements indicated soil freeze-thaw cycles and large variations of more than a month in the date when soil temperature dropped below the freezing point for the first time in areas with different vegetation cover and soil type (Online Resource Table S3). Again, soil freeze occurred first in dry mineral soils in open areas and latest in organic soils in wetland areas (Histosol). Thus, there were areas within the SMOS pixel for which the soil was frozen over the depth of 10-15 cm whereas other areas remained non-frozen for longer time. Satellite-observed soil freeze over the whole pixel was detected only after mid-December. Difference in soil wetness may also explain later (16 d) satellite soil freeze at CA-Ojp in 2010 compared to neighbouring site CA-Obs, while VAPend 10 took place only one day later than at CA-Obs. We compared the performance of the satellite-derived soil freeze state proxy for VAPend with the date of autumnal freeze from ERA-Interim reanalysis air temperature data. For all three VAPend dates, the agreement was better than that for the satellite soil freeze (Online Resource Table S2, Fig.  S1 c, f, i). Again, the date of autumnal freeze showed the highest correspondence with VAPend 10 (Fig. 3b). VAPend 10 , on average, occurred 5 days before the autumnal freeze.
The slope of the linear regression lines of the two interrelated freeze indicators with VAPend 10 was almost identical (~0.8). Autumn freeze reduced the carbon uptake. During warm winters, the time difference between VAPend and air and soil freeze increased. The mean time lag between VAPend 10 and first date of air freeze and satellite-observed partially frozen soil was 3 days longer for the warmest site FI-Hyy than for sites in northern Finland (2010-2016, Online Resource Table S4). For example, in 2011, the autumnal freeze was observed 24 days after VAPend 10 in FI-Hyy. During the warm winter, T 10 did not drop below − 1°C until 7 January 2012. Instead, the photosynthetic capacity of the forest decreased below GPI 10 by mid-December. Low CO 2 uptake above GPI 5 was observed until the end of 2011, and VAPend 5 occurred 2 days after the first date of the autumnal freeze. This suggests that other factors, such as short photoperiod, also exert control on VAPend.
However, there was considerable year-to-year variation in VAPend 10 (see Table 2). For the Finnish sites, air temperature-derived autumn freeze and the satellite-derived first date of partially frozen soil explained well the deviations from the multi-year mean (2010-2016) with R 2 0.82 and 0.60 (N 20, p < 0.001), respectively. Autumnal freeze from T 10 took place considerably earlier than VAPend 10 once during the period 2010-2016. This was the case in 2015 at FI-Sod when it was observed early on 15 October (doy 288) during the first cold period with sub-zero temperature (Fig. 5).
ERA-Interim global atmospheric reanalysis air temperature data used to derive the air temperature proxy (Fig. 3b)  Fig. 5 Time series of the gross photosynthesis index (GPI), air and soil temperature (5-cm depth), snow depth and satellite-derived soil freeze/thaw at Sodankylä in 2015. Red vertical lines indicate the end of the vegetation active period based on the GPI thresholds at the CO 2 flux measurement site FI-Sod. The black vertical line indicates first date of autumnal freeze from ERA-Interim reanalysis air temperature data. HA haplic arenosol, UG umbric gleysol, His histosol were also used as ancillary information to avoid false detections of soil freeze in summer from SMOS. During the retrieval of soil state, partially frozen soil and frozen soil are only assigned when a summer temperature flag has been released . In consequence, it is possible that partially frozen or even frozen soil state detections are not independent of air temperature data, and soil freeze is detected only one day after the summer flag has been released. This phenomenon was observed in our study during 2010-2011 for FI-Ken and during 2015 for FI-Hyy. However, when removing these three points from the data the quality of the linear regression equation between satellite soil freeze with VAPend 10 remained the same (N 22, R 2 0.83, p < 0.001).

Mapping of the End of Vegetation Active Period Using Proxy Indicators
The spatial distribution of VAPend 10 in the boreal forest zone in Finland and Canada was determined from the satellite-observed first date of partially frozen soil and the first date of autumnal freeze. For that, the linear regression equations in Fig. 3 were used to model VAPend 10 . The RMSE of the linear regression model for the calibration was 7.5 days for the satellite-based and 5.2 days for the air temperature-based proxy (Online resource Fig. S2). For the average period 2010-2016, VAP ended before 10 October in the northwestern part of Canada and shifted to 20 November to 9 December in Southern Quebec and the island of Newfoundland (Fig. 6a). In northwest Finland, VAPend 10 occurred in mid-October and shifted to 10 December in the most southern parts (Fig. 6b). In general, the estimations based on the autumnal freeze produced a similar but smoother spatiotemporal gradient in the timing of VAPend 10 as satellite soil freeze. Despite the fact that the proxy indicator from autumnal freeze produced more accurate estimates for the CO 2 flux measurement sites, the satellite-derived soil freeze allows the mapping of VAPend with higher spatial detail at a global scale.

Discussion
In our analysis of CO 2 flux measurement sites in the boreal zone, we found a strong connection between satellite-derived soil freeze and autumnal freeze from air temperature with VAPend. Our results are supported by Barr et al. (2009) who found a strong relationship of the first date of freeze with the end of photosynthetic activity at CA-Obs and CA-Ojp for an earlier period. In previous comparisons with satellite observations, only low correspondence of VAPend with estimates of landscape freeze from the SeaWinds scatterometer was observed for evergreen forests in northern America, including CA-Obs. This might be attributed to the low sensitivity of the K u -band radar to freeze-thaw processes in autumn (Kimball et al. 2004). In contrast to the SMOS L-Band, the SeaWinds radar signal is not directly sensitive to the soil processes, thus observing landscape freeze from vegetation and the very skin of the soil surface or snow cover during winter (Kimball et al. 2004). In the study by Barichivich et al. (2013) also, only a weak link between the end of growing season with the timing of autumn freeze from satellite observations and the end of thermal growing season was found on the continental and global scale. The correlation between daily landscape freeze from Special Sensor Microwave Imager (SSM/I) and the end of growing season from the NDVI was not significant. The authors suggested, therefore, a decoupling between autumn phenology and thermal conditions for the northern high latitudes. This difference, relative to our study, may be due to the fact that the end of season indicator from NDVI describes changes related to leaf senescence of deciduous vegetation that precedes VAPend 10 in boreal evergreen forest ecosystems.
In almost all observations, the GPI 10 was reached days or weeks before T 10 dropped below − 1°C and before satellite soil freeze, suggesting that the relationship between autumnal soil freezing and decline of photosynthetic capacity is correlative rather than causal. The strong correlation found between autumn freeze and VAPend 10 is probably due to the fact that both processes, soil freezing and transition of plants towards a dormant state, are causally driven by a cumulative effect of declining temperature, similar to the effects of rising temperature in spring (Mäkelä et al. 2004;Pulliainen et al. 2017). Furthermore, the slope of the regression is < 1 (Fig. 3), such that delay time between VAPend 10 and the autumnal freeze is longer, the warmer the autumn. This suggests that other factors, such as photoperiod, also affect the change in photosynthetic capacity and that they become relatively more important during warmer autumns. The importance of the freeze-based proxy indicators could, therefore, be reduced in the future as autumn frost and soil freeze will likely occur later in a warming climate.
In our comparison of GPI threshold values, time of year (or photoperiod) rather than time of autumn freeze was more important when the threshold percentage used was higher, as shown by the slopes of the regression of VAPend th against the autumnal freezing date (Online resource Figure S1). However, for high percentage GPI thresholds, this could simply be linked with the strongly reduced light levels after the equinox, rather than marking a transition towards a dormant state. Consistent with this, Suni et al. (2003) found for the FI-Hyy site that day length was a better predictor than temperature for the date when net ecosystem exchange of CO 2 fell below 20% of the maximum summer uptake. In contrast, we found that VAPend 5 usually coincided with or followed the date of autumnal freezing and that the slope of the regression was close to 1, suggesting that autumnal freezing may describe an actual temperature-related limitation to photosynthesis. This is further supported by the finding that the latest VAPend 10 date in our investigation was observed after the winter solstice on 25 December 2015 at FI-Hyy, showing that the GPI can remain even above the 10% level during days with the shortest period of light at this site. However, more research is needed to ascertain the relative significance of the different environmental drivers in the autumnal downregulation of photosynthesis.
Air temperature-based indicators were suggested earlier by Thum et al. (2009) as a proxy for VAPend in boreal forest ecosystems, including FI-Hyy and FI-Sod. In contrast to the uniquely defined temperature threshold that was used here, Thum et al. (2009) adapted the temperature-based indicators for each site, which hindered the application to wider areas. In agreement with other studies (e.g. Hollinger et al. 1999;Kimball et al. 2004), we observed a reduction of the photosynthetic carbon uptake with the freezing temperatures. The drop of T 10 below − 1°C does not necessarily lead to a sustained decrease in GPI. Photosynthetic activity can recover to a certain level when temperatures and other factors, such as photoperiod, are favourable (Fig. 5). This is not accounted for with the simple air temperature-based threshold. The method could be adjusted in a way that it indicates a more persistent air temperature decrease as, for example, done for the GPIbased VAPend (Fig. 2) or a cumulative frost sum.
The VAPend 10 transition date always occurred before the detection of the first date of partially frozen soil from SMOS (Fig. 3a). Soil temperature observations at FI-Sod and from sites with a similar soil type in the surrounding areas (Figs. 4 and 5) indicated an earlier occurrence of frozen soil than that determined for the SMOS pixel, but the GPI 10 was in most years exceeded before the occurrence of in situ soil freeze at FI-Sod (four out of six, e.g. Figs. 2 and 5), implying that soil freeze is not a necessary condition for VAPend 10 . There is a spatial mismatch between the footprint of the CO 2 flux measurement site, extending only several hundred meters from the site, and the SMOS footprint of 40-50 km in diameter. Soil freezing is strongly affected by the soil type (see Online Resources Table S3), snow conditions and the openness of the vegetation layer. As shown for the FI-Sod site, soils with finer texture and higher water content froze at a slower pace than dry coarse-textured soils with lower heat capacity (e.g. Fig. 5) and the delay time depended on air temperature development. At FI-Sod, the forested area on dry mineral soil, representing best the CO 2 flux footprint, covered only a small fraction (7%) of the SMOS pixel. Moreover, other land cover classes, such as large water bodies (e.g. CA-Obs 14% and CA-Qfo 23% of the pixel, see Online resources Table S1) and wetland areas (FI-Sod 13% of the pixel), affect the soil state retrieval. Thus, scaling errors lead to inevitable discrepancies between satellite soil freeze and site measurements and possibly higher scatter in results compared to autumnal freeze (Fig. 3). L-band active radar measurements could provide soil freeze and thaw state at higher spatial resolution (Derksen et al. 2017), but the NASA Soil Moisture Active Passive (SMAP) radar failed in 2015. The currently operating Advanced Land Observing Satellite-2 Phased Array L-Band Synthetic Aperture Radar (ALOS-2 PALSAR) has a narrow swath and, thus, a too long revisit time (14 days) for the soil freeze and thaw monitoring. For site-specific more detailed studies, tower-based L-Band radiometer instruments, such as the ELBARRA-II (ETH L-Band Radiometer) (Schwank et al. 2010) installed at Sodankylä, could be utilized.
In comparison to satellite-derived soil freeze, the global ERA-Interim air temperature data set with a coarser resolution appears to represent the condition at the CO 2 flux measurement sites rather well (Figs. 4 and 5). The spatial variation of air temperature is smaller than for soil temperature, and, additionally, in situ air temperature observations that are ingested in the data assimilation system to produce the ERA-Interim global reanalysis data (Dee et al. 2011) probably stem from nearby synoptic weather stations, for example at Sodankylä. However, the accuracy of the reanalysis temperature is expected to be lower for areas with lack of near-surface observations and especially at high latitudes, because of difficulties in modelling surface temperatures over snow cover (Simmons et al. 2010).
We presented maps of the timing of VAPend in evergreen needle-leaved forest in the boreal region based on proxy indicators from air temperature and satellite-derived soil freeze (Fig. 6). We selected VAPend 10 as a reference value for the model calibration, although VAPend 5 may be closest to the dormant state of vegetation from the three tested GPI threshold values. This was done since, the calculation of VAPend 5 from CO 2 fluxes has higher uncertainty than the determination of VAPend 10 because of lower light levels and resulting gaps in the GPI time series. The model parameterization for the calculation of VAPend 10 from the two freeze indicators is based on a limited number of observation sites and years. Future work is, therefore, needed to validate the results throughout the boreal forest region, especially for areas that are not well represented by the sites in this study, such as the northern boreal zone in North America and the southernmost part of the boreal zone in Europe (Fig. 6). However, publicly available eddy covariance data are scarce and sites are not well distributed in this region (Falge et al. 2017).
Furthermore, the possibilities for comparison with other independent data sets are limited. Recently, maps of the end of growing season for the northern hemisphere (1999-2013) based on a vegetation index (Phenological Index) were provided by Gonsamo and Chen (2016). For the evergreen needle-leaved boreal forest sites, the correspondence with gross primary production (GPP)-based autumn phenology indicators was only moderate (r 0.34, RMSE 16 days) due to the time delay between greenness decrease and cessation of photosynthesis, and noise in the satellite observation during the transition period . In our study, the multi-year mean VAPend dates (Table 2) occurred when solar elevation at solar noon was below 30°at all sites; hence, uncertainties in reflectance retrievals from optical satellite observations increase due to a longer atmospheric path length and shadowing effects of the canopy (Gutman 1991;Privette et al. 1995). A cutoff threshold of the solar zenith angle > 70 degrees is applied for satellite retrievals of SIF from GOME-2 (Joiner et al. 2014) and also for operational snow cover products (Metsämäki et al. 2015) that could otherwise provide complementary information.

Conclusion
We investigated whether a new satellite-derived soil freeze prototype product from SMOS can be used as a proxy indicator for the mapping of the end of the vegetation active period in boreal evergreen needle-leaved forest. Our results showed that the "partially frozen" soil state derived from satellite observations can be utilized to infer the date when the photosynthetic capacity of the forest decreases below 10% of the growing season maximum with good accuracy (R 2 0.84, RMSE 7.5 days) at six CO 2 flux measurement sites in the boreal forest zone in Finland and Canada. Results for an air temperature-based indicator, defined as the day when 10day average air temperature decreased below − 1°C, were consistent with satellite-derived soil freeze. Nevertheless, air freeze showed a stronger correlation with the end of vegetation active period than satellite-derived soil freeze (R 2 0.92, RMSE 5.2 days). Lower accuracy for the satellite-based proxy may be due to contributions to the SMOS signal from land cover and soil types, other than those dominant at the CO 2 flux measurement sites. Although further studies covering more sites and observation years are recommended, the results suggest that autumnal air temperature and satellitebased soil freeze can be utilized for the mapping of the start of the dormant season in boreal evergreen needle-leaved forest at continental to global scale. The maps complement existing satellite-based data on the end of the growing season in deciduous forests in the boreal region from vegetation indices and provide spatially explicit information that could help in monitoring the end of the vegetation active period and in the calibration and validation of biosphere models.
Acknowledgements Open access funding provided by Finnish Environment Institute (SYKE). The work was supported by the European Commission through the Life + project MONIMET (K.B, A.N.A, M.A., A.M., P.K.); Grant agreement LIFE12 ENV/FI000409) and the Envibase project funded by Ministry of Finance, Finland. The work of K.R has been financially supported by the Academy of Finland, decision 309125, Academy of Finland Centre of Excellence, decision 118780, and "SMOS + Frost2 Study" funded by ESA (contract 4000110973/14/NL/FF/lf). We thank Jaakko Ikonen (Finnish Meteorological Institute) for providing information on the distribution of soil types in the Sodankylä area. We acknowledge Fluxnet-Canada and the Canadian Carbon program for the AmeriFlux data records: CA-Obs, CA-Ojp, CA-Qfo. In addition, funding for the AmeriFlux data resources was provided by the U.S. Department of Energy's Office of Science. We are grateful to the site data provider Hank Margolis (Université Laval, Ca-Qfo) for making the flux measurements available. Flux measurements and ancillary data for the Hyytiälä site at SMEAR (Station for Measuring Ecosystem-Atmosphere Relations) research stations of the University of Helsinki were downloaded from the open data service of the University of Helsinki (https://avaa.tdata.fi). GlobCover (Arino et al. 2012) (Global Land Cover Map, © ESA 2010 and UCLouvain) was downloaded from the ESA DUE GlobCover website: http:// due.esrin.esa.int/page_globcover.php. SMOS data were obtained from the "Centre Aval de Traitement des Données SMOS" (CATDS), operated for the "Centre National d'Etudes Spatiales" (CNES, France) by IFREMER (Brest, France). We thank the ECMWF for making ERA-Interim data available (http://apps.ecmwf.int/datasets/). We also thank two anonymous reviewers for their helpful comments to the manuscript.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.