Vegetation Growth Analysis of UNESCO World Heritage Hyrcanian Forests Using Multi-Sensor Optical Remote Sensing Data

: Freely available satellite data at Google Earth Engine (GEE) cloud platform enables vegetation phenology analysis across different scales very efficiently. We evaluated seasonal and annual phenology of the old-growth Hyrcanian forests (HF) of northern Iran covering an area of ca. 1.9 million ha, and also focused on 15 UNESCO World Heritage Sites. We extracted bi-weekly MODIS-NDVI between 2017 and 2020 in GEE, which was used to identify the range of NDVI between two temporal stages. Then, changes in phenology and growth were analyzed by Sentinel 2-derived Temporal Normalized Phenology Index. We modelled between seasonal phenology and growth by additionally considering elevation, surface temperature, and monthly precipitation. Results indicated considerable difference in onset of forests along the longitudinal gradient of the HF. Faster growth was observed in low- and uplands of the western zone, whereas it was lower in both the mid-elevations and the western outskirts. Longitudinal range was a major driver of vegetation growth, to which environmental factors also differently but significantly contributed ( p < 0.0001) along the west-east gradient. Our study developed at GEE provides a benchmark to examine the effects of environmental parameters on the vegetation growth of HF, which cover mountainous areas with partly no or limited accessibility. and NIR at 10 m, red-edge and SWIR at 20 m, and atmospheric bands at 60 m spatial resolution. We selected Level-2A surface reﬂectance product available at GEE and computed NDVI. This corresponded to a processed dataset made available from European Space Agency at GEE cloud platform.


Introduction
Studying seasonal or inter-annual phenology is among the most significant topics in forest ecosystems, with multiple implications for both research and practice. As a result, phenological research aims to study ecosystem behavior in relation to climatic or environmental changes [1,2] and thus, the results play an important role in ecosystem management in response to its long-term, inter-annual, or seasonal variabilities. Remote sensing data and methods have been numerously employed to mimic phenological behavior and fluctuations, often as the most efficient and practical methods for gaining a better understanding of vegetation phenology dynamics at regional and global scales [3,4].
Optical satellite remote sensing provides information for studying large-scale forest changes in near real-time with a comparable spatial and higher temporal resolutions compared with field surveys, which are often constrained by logistics and difficult-toaccess terrain. Among the available methods, those based on broadband vegetation indices (VIs) derived from multispectral data have been widely used [5][6][7] for vegetation classification [8], phenological monitoring [9], change detection [10], and the retrieval of forest biophysical and structural attributes [11]. Since 2015, the Sentinel-2 archive has provided continuous imaging information on the Earth's surface, with which researchers have lately begun to investigate ecosystem changes, inter alia using approaches based on common broadband VIs. For vegetation-related studies, the normalized difference vegetation index (NDVI) is the most common and commonly used metric [12,13]. The reader is referred to [14] for a synthetic review of NDVI, including its origin, availability, advantages, and limitations.
Statistical approaches based on time series analysis of VIs from satellite data have been widely applied to track changes in forest and land cover dynamics in respect to their spectral reflectance from the vegetation overstorey. Other factors such as land surface temperature, rainfall, terrain, and soil quality were thought to be significant contributors to changes in vegetation health [15] and can be coupled with VIs for phenology monitoring over forest ecosystems. NDVI provides information about vegetation greenness, which records photosynthesis activity of a plant or tree leaves [16]. To understand change in photosynthesis activity between different seasons or vegetation growth cycle, change in NDVI values may not be sufficient for the entire temporal sequence of a vegetation growth period [9,17,18]. To solve this numerical limitation of NDVI, a new measure called the temporal normalized phenology index (TNPI) was proposed and recommended as a superior option for analyzing the temporal phenology cycle between two time steps of the maximum and minimum plant growth period [9]. Furthermore, with time series Landsat-8 data, the sensitivity of fluctuating NDVI to individual remote sensing-derived topography components and land surface temperature was effectively tested using TNPI, with the main benefit that it reduced the need for long-term monthly records to understand the forest phenology [9].
In this study, we aim to assess a recently developed remote sensing-only approach based on TNPI to understand the phenology-based vegetation growth and the effects of temperature, precipitation, and elevation variations on Hyrcanian forests (HF) along its entire longitudinal gradient. Apart from areas located in the Republic of Azerbaijan, the old-growth HF are mainly stretched as a thin belt (20-70 km wide and 800 km long) of mostly broadleaf deciduous forests located between the northern slopes of the Alborz Mountains and the southern coast of the Caspian Sea (CS) in Iran. They comprise an altitudinal distribution from the sea level to over 2000 m above sea level [19,20]. Forest vegetation above 2000 m is gradually replaced by forest/steppe or steppe vegetation in the form of ecotones, with their extent and altitudinal level controlled by climate, anthropology (human settlements/cattle grazing), or both. The three northern provinces of Guilan, Mazandaran, and Golestan (from west to east) occupy nearly 1.9 million ha of forests in various qualitative stages from fully degraded to richly stocked stands. In this realm, the history of HF can reveal remarkable details about past Quaternary vegetation of the northern hemisphere. The international attention to HF has been raised upon the inscription of multiple sites therein as UNESCO Natural Heritages in 2019 (https://whc.unesco.org/en/list/1584/; accessed on 24 July 2021). In addition, several long-term challenges during the recent decades, including extensive land-use conversions, cattle grazing by local inhabitants, illegal poaching, illegal harvest, and wood smuggling, are among the major threats to HF and thus call for their constant monitoring in terms of structure, composition, and function, also with respect to their phenology.
To this aim, we constrained our analysis to the 15 structurally and compositionally different sites distributed all over the HF latitudinal and longitudinal range. These sites correspond to the full list of UNESCO-inscribed sites in 2019. We followed a two-phase approach starting from Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data at 250 m spatial resolution to uncover large-scale phenological patterns and then continued with Sentinel-2 data on 10 m spatial resolution for local-scale analysis of vegetation growth. The hypothesis was that the previously proven abilities of the TNPI in highlighting the phenological patterns when applying Landsat-8 data [9] will be used to understand the vegetation growth between onset and ending of phenological phases when applying the methodology on Sentinel-2 data with higher spatial resolution. A secondary aim was also to concentrate a major part of the analysis on the cloud-computing capabilities of Google Earth Engine, which would lead to developing a full remote sensing, GEE-based workflow for future analyses with potential implications across various geographical scales and timespans.

Subsection
The study areas comprise all 15 UNESCO-inscribed sites distributed across the 8000 km 2 of temperate deciduous HF that run along the southern coast of the CS and through the provinces of Guilan, Mazandaran, and Golestan. ( It is distinguished by a varied plant environment and a variety of ecological conditions. The 15 inscribed sites ( Figure 1 and Table 1) reflect the majority of compositional, structural, and functional variations existing within the HF and are selected due to their unique characteristics as summarized and evaluated by the International Union of Conservation of Nature (IUCN) World Heritage Panel progress report [24], to which the reader is referred for further information on the selection criteria.

MODIS NDVI Data
We acquired time series of 250 m Terra MODIS Vegetation Indices (MOD13Q1, version 6) NDVI provided by the NASA LP DAAC (https://lpdaac.usgs.gov/products/mod13q1v006/; accessed on 10 July 2021), which is archived in the Google Earth Engine (GEE) platform (https://earthengine.google.com/, Mountain View, CA, USA; accessed on 10 July 2021). The NDVI follows the basic equation where ρ and ρ correspond to MODIS band 1 (620-670 nm) and band 2 (841-871 nm) spectral reflectance values. MOD13Q1 is derived from atmospherically corrected bidirectional surface reflectance imagery and contains VI data as well as the pixel reliability layer required for quality checking [25]. We computed 16-days NDVI maximum composite values for the selected sites using the JavaScript code editor in the GEE platform from 2017 to 2020. We extracted time series mean NDVI values for the selected 15 World Heritage forest polygons within the Hyrcanian forests. We selected a time period of 4 years to understand the general phenology pattern of each forest site.

Sentinel-2 Data
Sentinel-2 is a wide-swath, high-resolution, multispectral imaging mission with a global 5-day revisit frequency. The Sentinel-2 Multispectral Instrument (MSI) samples 13 spectral bands: visible and NIR at 10 m, red-edge and SWIR at 20 m, and atmospheric bands at 60 m spatial resolution. We selected Level-2A surface reflectance product available at GEE and computed NDVI. This corresponded to a processed dataset made available from European Space Agency at GEE cloud platform.

MODIS NDVI Data
We acquired time series of 250 m Terra MODIS Vegetation Indices (MOD13Q1, version 6) NDVI provided by the NASA LP DAAC (https://lpdaac.usgs.gov/products/ mod13q1v006/; accessed on 10 July 2021), which is archived in the Google Earth Engine (GEE) platform (https://earthengine.google.com/, Mountain View, CA, USA; accessed on 10 July 2021). The NDVI follows the basic equation where ρ Red and ρ NIR correspond to MODIS band 1 (620-670 nm) and band 2 (841-871 nm) spectral reflectance values. MOD13Q1 is derived from atmospherically corrected bidirectional surface reflectance imagery and contains VI data as well as the pixel reliability layer required for quality checking [25]. We computed 16-days NDVI maximum composite values for the selected sites using the JavaScript code editor in the GEE platform from 2017 to 2020. We extracted time series mean NDVI values for the selected 15 World Heritage forest polygons within the Hyrcanian forests. We selected a time period of 4 years to understand the general phenology pattern of each forest site.

Sentinel-2 Data
Sentinel-2 is a wide-swath, high-resolution, multispectral imaging mission with a global 5-day revisit frequency. The Sentinel-2 Multispectral Instrument (MSI) samples 13 spectral bands: visible and NIR at 10 m, red-edge and SWIR at 20 m, and atmospheric bands at 60 m spatial resolution. We selected Level-2A surface reflectance product available at GEE and computed NDVI. This corresponded to a processed dataset made available from European Space Agency at GEE cloud platform.

Elevation Data
A Shuttle Radar Topography Mission (SRTM)-based digital elevation model (DEM) was used to extract the elevation information of HF [26]. We used SRTM version 4 tagged data available on GEE, with the spatial resolution of 90 m, orthorectified and preprocessed by the data provider.

NDVI Curve Fitting
Time series MODIS NDVI was fitted with a double-logistic function where NDVI min and NDVI max are the minimum and maximum values measured in the winter and summer, respectively, start of season (SOS) and end of season (EOS) are the inflection points when the curve rises and falls, and slope1 and slope2 are the rates of increase and decrease of the curve at the inflection points ( Figure 2) [27,28]. This function describes asymmetrical patterns, leading to a reliable estimation of the trajectory in canopy greenness [27]. We extracted phenology parameters based on the 4 years of time series data. Length of growing season (LOS) was computed by calculating the difference between SOS and EOS [28].

Elevation Data
A Shuttle Radar Topography Mission (SRTM)-based digital elevation model (DEM) was used to extract the elevation information of HF [26]. We used SRTM version 4 tagged data available on GEE, with the spatial resolution of 90 m, orthorectified and preprocessed by the data provider.

NDVI Curve Fitting
Time series MODIS NDVI was fitted with a double-logistic function where NDVI and NDVI are the minimum and maximum values measured in the winter and summer, respectively, start of season (SOS) and end of season (EOS) are the inflection points when the curve rises and falls, and slope1 and slope2 are the rates of increase and decrease of the curve at the inflection points ( Figure 2) [27,28]. This function describes asymmetrical patterns, leading to a reliable estimation of the trajectory in canopy greenness [27]. We extracted phenology parameters based on the 4 years of time series data. Length of growing season (LOS) was computed by calculating the difference between SOS and EOS [28].

Computation of TNPI From Multi-Temporal NDVI
We selected Sentinel-2 datasets between 2017 and 2020 for two different months defining the maximum and minimum values from the MODIS NDVI-derived SOS and EOS phenology parameters of each forest site. Each monthly composite dataset included filtering and cloud-masking for the selected dates. Sentinel-2 image collections comprised NIR and Red spectral bands at 10 m spatial resolution that were used for the computation of the NDVI.
We computed the Temporal Normalized Phenology Index (TNPI) [9] to quantify the change in trajectories of NDVI during two-time steps of the vegetation growth cycle. TNPI Figure 2. Concept of curve-fitting [27] and identification of SOS and EOS for TNPI computation. The highlighted red area represents the vegetation growth period between SOS and EOS.

Computation of TNPI from Multi-Temporal NDVI
We selected Sentinel-2 datasets between 2017 and 2020 for two different months defining the maximum and minimum values from the MODIS NDVI-derived SOS and EOS phenology parameters of each forest site. Each monthly composite dataset included filtering and cloud-masking for the selected dates. Sentinel-2 image collections comprised NIR and Red spectral bands at 10 m spatial resolution that were used for the computation of the NDVI.
We computed the Temporal Normalized Phenology Index (TNPI) [9] to quantify the change in trajectories of NDVI during two-time steps of the vegetation growth cycle. TNPI can be computed between two NDVI images and represents changes that occurred between two phenological growth stages of each selected forest site. Khare et al. (2017) proposed the TNPI concept as follows: In this case, we considered the time between SOS and EOS, where NDVI corresponds to SOS and EOS date of each forest site identified previously using MODIS data. Therefore, in this case, TNPI will represent the vegetation growth period of individual forest sites. The modified TNPI equation will be as follows: We computed TNPI Growth for the entire HF both (1) using MODIS NDVI data to understand the general pattern of vegetation growth at coarser spatial resolution, and (2) using high spatial resolution Sentinel-2 NDVI data over the selected 15 forest polygons for spatially detailed investigation.

Precipitation Data
We used the GEE platform to extract precipitation chronologies from ERA-5 monthly averages (5th generation) and understand the effect of precipitation on vegetation growth. This is climate reanalysis data provided by ECMWF (European Centre for Medium-Range Weather Forecasts) with a spatial resolution of 31 km [29]. We extracted precipitation data by averaging the time period between 2017 and 2020 for SOS and EOS months separately (Equation (2)). The change in precipitation (∆PP) was computed using the following equation: where PP EOS is precipitation at the EOS and PP SOS is the precipitation at the SOS.

MODIS Land Surface Temperature Data
We used the GEE platform to extract MODIS average 8-day land surface temperature (LST) (MOD11A2 version 6, 1 km spatial resolution) data to understand the effect of temperature on vegetation growth. We extracted LST by averaging the time period between 2017 and 2020 for SOS and EOS months separately (Equation (2)). The change in LST (∆LST) was computed using the following formula:

Grouping Analysis
We initially studied the effect of environmental factors of elevation, LST, and precipitation on vegetation growth for MODIS-based TNPI Growth covering the entire range of HF along the west-east gradient. We randomly generated 1000 sample points and then extracted longitude, raster data values of MODIS based TNPI Growth , elevation, LST, and precipitation using the point sampling tool in QGIS version 3.16 [30]. To identify the longitudinal changes in the vegetation growth and the effects of the above-mentioned environmental factors, exploratory cluster analysis was conducted on 1000 locations using the Grouping Analysis tool within ArcMap version 10.7 [31]. This tool uses a K-means algorithm to find a solution where different parameters are grouped together such that attributes within one group are possibly similar, while groups are possibly different from each other [32].

Modelling Sentinel-2-Derived TNPI Growth with Environmental Factors
HF are mainly stretched in the longitudinal direction (Figure 1), thus, the 15 inscribed forest sites were divided into three forest zones ( Figure 1): (i) the west forest (11,12,13,14,15), (ii) the east forest (1,2,3,4,5), and (iii) the middle forest (6,7,8,9,10). We randomly generated 100 sampling points within each forest shapefile and then extracted raster data values of elevation, LST, precipitation, and Sentinel-2-based TNPI Growth using the procedures described in Section 2.2.8. The extracted values were combined and categorized for each forest zone. We used these extracted values to model TNPI Growth by multiple regressions as a function of elevation, LST, and precipitation for the west, east, and middle forest zones separately. A global ordinary least squares (OLS) model was first developed for each forest zone utilizing TNPI Growth as the response variable. Elevation, ∆PP, and ∆LST were all included as explanatory factors. A second-order polynomial form was used to express each global model: where Y represents the response variable (TNPI Growth ), α is the intercept, β 1 . . . β 6 are regression coefficients of the respective explanatory variables, and ε is the random error term. In addition to using each predictor (elevation, ∆PP, and ∆LST) as a factor, these predictors were allowed to interact with one another [9,33] to enable including interactions terms. Each global OLS model was screened in a thorough model selection procedure via the repeated model fitting of all possible submodels as described by [34]. The second-order Akaike Information Criterion (AICc) was calculated for each TNPI Growth model between the SOS and EOS time points to identify the submodel with the highest AICc [35,36]. The "MuMin" package in R was applied to calculate the subsets and AICc values and to select the submodels [34].

Phenology Patterns Derived from MODIS NDVI
The NDVI followed a bell-shaped pattern, with a gradual increase followed by a decline after a period of time. The NDVI was well-represented by the double-logistic function throughout the season.  Table 2). This resulted in average SOS length of two weeks in the eastern forest, 20 days in the western forest, and two weeks in middle forests.   The eastern forest's average EOS varied through the first week of November (DOY 306-310) along the longitudinal gradient. The average EOS in the western forest varied from mid-October (DOY 283) to the beginning of December (DOY 335). Similarly, the average EOS in the middle forest varied from the last week of October (DOY, 304) until the mid of November (DOY 318). The average EOS length was found to be less than a week for the eastern forest, approximately eight weeks in the western forest, and two weeks in the middle forest. Furthermore, the average NDVI for the eastern forest ranged from 0.42 to 0.84, whereas the average NDVI for the western forest ranged from 0.39 to 0.82. Similarly, the average NDVI range for the middle forest varied between 0.28 and 0.75. Furthermore, the standard deviation for Slope1 was the largest in the west forest (0.031), followed by that of the east forest (0.011) and the mid forest (0.005). Similarly, the standard deviation for Slope2 was highest in the west forest (0.014), followed by that of the mid forest (0.008) and the east forest (0.007) (Figure 3, Table 2).

MODIS-Based Vegetation Growth Analysis and Effects of Environmental Factors
MODIS-derived vegetation growth and effects of environmental factors along the longitudinal direction were observed by means of grouping analysis. Overall, longitude was a key attribute in grouping analysis, with the highest impact on all other parameters, Remote Sens. 2021, 13, 3965 9 of 17 including elevation, ∆PP, and ∆LST, with higher R 2 (0.7203) compared to that of other environmental factors (Figure 4, Table 3).

MODIS-Based Vegetation Growth Analysis and Effects of Environmental Factors
MODIS-derived vegetation growth and effects of environmental factors along the longitudinal direction were observed by means of grouping analysis. Overall, longitude was a key attribute in grouping analysis, with the highest impact on all other parameters, including elevation, ΔPP, and ΔLST, with higher R 2 (0.7203) compared to that of other environmental factors (Figure 4, Table 3).  Results of vegetation growth-environmental factor grouping indicated that the western forests included Groups 1-3 and 5, while the middle forests consisted of Groups 1, 4, and 5 ( Figure 5, Table 3). The eastern forests were mainly dominated by Group 4 along with a slight combination of Groups 1 and 2. According to the parallel box plot, areas of Group 4 were situated at higher altitudes with average elevation (1020.45 m), indicating that plant growth was slowed down due to lower ΔLST (−2.62 °C) and ΔPP (−10.12 mm). Group 5, on the other hand, which came at average height but elevation below the mean value (624.71 m) grew faster than Group 4 at a higher elevation. As a result of the lower ΔLST (−1.99 °C) and greater ΔPP (112.05 mm), plant growth was slightly higher than the mean value. Further, Group 1 was located at a slightly lower altitude than Group 5, but at a higher elevation (1590 m), thus indicating a faster growth than Group 5. This could be related to ΔLST (−2.31 °C) with a value lower than the mean as well as ΔPP (41mm) with a value higher than the mean. In case of Group 2, longitude was close to the mean value, the elevation (380.74 m) was minimal compared to that of other locations, and ΔLST (−1.18 °C) and ΔPP (59.51 mm) values were both higher, resulting in increased plant growth. However, Group 3, with longitude values below the mean at moderate elevation (1134.89 m) showed much less growth than other locations, since ΔLST (−1.93 °C) was close to the  Results of vegetation growth-environmental factor grouping indicated that the western forests included Groups 1-3 and 5, while the middle forests consisted of Groups 1, 4, and 5 ( Figure 5, Table 3). The eastern forests were mainly dominated by Group 4 along with a slight combination of Groups 1 and 2. According to the parallel box plot, areas of Group 4 were situated at higher altitudes with average elevation (1020.45 m), indicating that plant growth was slowed down due to lower ∆LST (−2.62 • C) and ∆PP (−10.12 mm). Group 5, on the other hand, which came at average height but elevation below the mean value (624.71 m) grew faster than Group 4 at a higher elevation. As a result of the lower ∆LST (−1.99 • C) and greater ∆PP (112.05 mm), plant growth was slightly higher than the mean value. Further, Group 1 was located at a slightly lower altitude than Group 5, but at a higher elevation (1590 m), thus indicating a faster growth than Group 5. This could be related to ∆LST (−2.31 • C) with a value lower than the mean as well as ∆PP (41mm) with a value higher than the mean. In case of Group 2, longitude was close to the mean value, the elevation (380.74 m) was minimal compared to that of other locations, and ∆LST (−1.18 • C) and ∆PP (59.51 mm) values were both higher, resulting in increased plant growth. However, Group 3, with longitude values below the mean at moderate elevation (1134.89 m) showed much less growth than other locations, since ∆LST (−1.93 • C) was close to the mean, and ∆PP (23.99 mm) was below the mean, resulting in its poor vegetation growth ( Figure 5 and Supplementary data).
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 17 mean, and ΔPP (23.99 mm) was below the mean, resulting in its poor vegetation growth ( Figure 5 and Supplementary data).

Multiple Regression Relationships between Sentinel-2-Derived TNPIGrowth and Environmental Variables
In the western forest, precipitation and elevation significantly contributed to vegetation growth (p < 0.0001), followed by a significant interaction between precipitation and elevation (Adj. R 2 0.5224, p < 0.0001), resulting in an increase in TNPIGrowth during the maximum precipitation period. Therefore, change in vegetation greenness was higher due to higher precipitation and increasing elevation. ΔLST showed no effect on the western forest, as the interaction between ΔLST and ΔPP was insignificant (Figures 6 and 7, Table 4).
In the eastern forest, ΔPP, ΔLST, and elevation all showed significant influence on TNPIGrowth (p < 0.0001). The observations also showed a significant interaction between ΔPP and ΔLST (Adj.R 2 = 0.6024, p < 0.0001), resulting in a major increase in TNPIGrowth during the minimum and medium precipitation periods (Figures 6 and 7, Table 4).
Similarly, all factors but ΔLST showed significant impact on TNPIGrowth in the middle forest, and a significant interaction was observed between elevation and ΔPP (Adj.R 2 =

Multiple Regression Relationships between Sentinel-2-Derived TNPI Growth and Environmental Variables
In the western forest, precipitation and elevation significantly contributed to vegetation growth (p < 0.0001), followed by a significant interaction between precipitation and elevation (Adj. R 2 0.5224, p < 0.0001), resulting in an increase in TNPI Growth during the maximum precipitation period. Therefore, change in vegetation greenness was higher due to higher precipitation and increasing elevation. ∆LST showed no effect on the western forest, as the interaction between ∆LST and ∆PP was insignificant (Figures 6 and 7, Table 4).
TNPIeast, and TNPImid represent the vegetation growth model for western, eastern, and middle forest zones, respectively, as described in Section 2.2.9.   Cross-sectional plots of three forest zones for Sentinel-2-derived TNPIGrowth models with an interaction between elevation and precipitation, and elevation and LST. The first panel depicts the western forests, the second panel depicts the eastern forests, and the last panel shows the middle forest zone as described in Section 2.2.9.

The Observed Phenology Patterns
The observed patterns for SOS and EOS along the longitudinal gradient of HF were in line with the observations throughout the entire region, as well as with those of the only previously published remote sensing-based study [37], which was carried out based on Figure 7. Cross-sectional plots of three forest zones for Sentinel-2-derived TNPI Growth models with an interaction between elevation and precipitation, and elevation and LST. The first panel depicts the western forests, the second panel depicts the eastern forests, and the last panel shows the middle forest zone as described in Section 2.2.9. Table 4. The final selected models and their performances for three forest zones using Sentinel-2 NDVI data. TNPI west , TNPI east , and TNPI mid represent the vegetation growth model for western, eastern, and middle forest zones, respectively, as described in Section 2.2.9. In the eastern forest, ∆PP, ∆LST, and elevation all showed significant influence on TNPI Growth (p < 0.0001). The observations also showed a significant interaction between ∆PP and ∆LST (Adj.R 2 = 0.6024, p < 0.0001), resulting in a major increase in TNPI Growth during the minimum and medium precipitation periods (Figures 6 and 7, Table 4).

Selected
Similarly, all factors but ∆LST showed significant impact on TNPI Growth in the middle forest, and a significant interaction was observed between elevation and ∆PP (Adj.R 2 = 0.5618, p < 0.0001), resulting in an increasing trend of TNPI Growth for all precipitation categories (Figures 6 and 7, Table 4).

The Observed Phenology Patterns
The observed patterns for SOS and EOS along the longitudinal gradient of HF were in line with the observations throughout the entire region, as well as with those of the only previously published remote sensing-based study [37], which was carried out based on coarser-resolution GIMMS NDVI3g time series data. The latter resulted in averaged DOY 77-158 (mean DOY 103) for SOS, as well as DOY 256-318 (mean DOY 299) for EOS across the HF, which was close to our TNPI-based estimates, but with two major differences: (1) the NDVI3g-based estimates did not imply any differentiation based on the longitudinal gradient, and (2) their average estimates for vegetation onset and senescence were generally a bit earlier than those from our analysis, which might be due to the reports that suggested that GIMMS data return earlier phenological estimates [38]. The study by [37] also confirmed an overall lengthening of the growing season across the entire Hyrcanian region, mainly as a result of climate change-induced factors. Whereas this study has been the only published example of analysis on spatiotemporal vegetation dynamics and phenological parameters across the HF, our study succeeded to expand this into a more spatially differentiated analysis, considering not only higher spatial resolution for the entire Hyrcanian region, but also by further focusing on the previously described World Heritage Sites on 10 m spatial resolution within a dense time series, which inherently contributes to reducing the multiple caveats of coarse resolution data like NDVI3g to assess vegetation phenology across heterogeneous landscapes, like mixed signal [39] and frequent missing values [40]. Among other biogeographic zones, [28] also reported highly heterogeneous patterns of forest phenology phases like bud break along the longitudinal gradients across black spruce (Picea mariana E.E. STERNS) dominated stands in boreal Quebec-Canada using MODIS time series, which yet cannot be directly compared with our results within a temperate zone. This also applies to other studies encompassing multiple biogeographic zones generally reported on a coarser spatial scale by a previous MERIS-based analysis of vegetation phenology across Europe [41].
A striking observation was in terms of differentiated EOS observed along the westeast gradient, where average EOS showed a reduction from east to west, from less than a week in the eastern zone to approximately eight weeks in the west. The climate over HF has been reported to range from warm Mediterranean in the east to Mediterranean in the west [20], implying a variance in duration of vegetation activities from the very eastern part (receiving ca. 700 mm precipitation) to the very western part (receiving ca. 1700 mm precipitation). Apart from that, the damming effect of the Alborz mountains for humidity (affecting plant activity and growth) reduces from west to east, as the distance between the mountain and the CS gradually increases (see, e.g., the longitudinal variations in environmental gradients in Figure 4). Previous field-based vegetation studies like [42] and [43] have previously shown that the most influential environmental factors that drive the floristic composition within the Hyrcanian forests include elevation above sea level, followed by annual precipitation and temperature, while topographic features are mainly affective within a given elevation zone or vegetation type. Other microclimatic factors, in particular soil physical and chemical properties, are of secondary importance [42] and were also out of the scope of our study.

MODIS-Based Vegetation Growth Analysis
Our vegetation growth-environmental factor grouping analysis on MODIS time series data suggested the longitude to be among the most important factors controlling the vegetation dynamics and growth across Hyrcanian forests. This was complementary to a recently reported benchmark study on floristic composition based on detrended correspondence analysis (DCA) and canonical correspondence analysis (CCA) of ca. 1600 vegetation plots of 802 vascular taxa covering the entire longitudinal and latitudinal range of HF [43]. Whereas proving the sole the effect of altitude on vegetation activities within the Hyrcanian zone is not new [44], the former recent study suggested the longitudinal range from the west to the east to be a major driver of species dynamics across the region. Our remote sensing analysis showed that it is presumably also a driver of vegetation activities and growth. However, we absolutely suggest this to be cautiously interpreted due to the absence of large-scale detailed field data on vegetation growth across the region.
Moreover, the absolute domination of clustered Groups 1-3 and 5 in the western longitudinal range generally suggested faster vegetation growth concentrated in low-and uplands of the western zone, whereas both mid-elevations and the western outskirts of the Hyrcanian forests (in all elevations) were associated with lower vegetation growth, i.e., Group 3. The higher growth observed in western longitudes can be attributed to the closer distance between the CS and the Alborz Mountains, which fosters a higher precipitation rate throughout the year [43]. A contrary observation was made almost homogenously across the eastern zone, which was dominated by low-growth Group 4 (see Figure 5), regardless of altitude. However, one may note the generally farther distance to the CS in the eastern zone compared with other longitudinal zones. CS has been considered to be the main source of precipitation across the region [45]. The highest level of precipitation that occurs in autumn is caused by the higher temperature of the sea water, leading to its more rapid evaporation, lowering air stability and thus creating a latitudinal gradient in precipitation between lowland, mountainous, and the upper-mountainous areas [43]. Somewhat similarly, the dominance of Groups 1,4 and 5 in the mid longitudinal range can be partially due to the rather wide sea-mountain distance and comparably smooth latitudinal gradient (resulting in lower ∆LST and ∆PP). This pattern showed a concentration of samples grouped in 1 in higher altitudes, that is presumably portraying the alpine habitats of oriental beech (Fagus orientalis Lipsky.) as the main tree community occurring in both pure and mixed forms in the mid-longitudinal range of HF. There are indeed a few national and international published field-based studies on the effect of elevation above sea level on tree growth indicators (growing stock, DBH growth, aboveground biomass), which suggested its main role on tree growth within beech communities [46,47]. However, it is rather uncertain to fully assign the observed patterns to specific forest communities at this spatial resolution, since the occurrence of communities and vegetation alliances in the Hyrcanian region, in particular in latitudinal gradient, is a function of multiple environmental and anthropogenic factors, which partly replace the tree communities with ecotonal and open-land species as a result of climate, but also factors like grazing, upland cultivation, and eradication by land-use conversion [47].

Multiple Regression Results on the Sentinel-2 Level
The observations made on the MODIS level across the entire Hyrcanian region were largely confirmed at higher spatial resolution based on spatial locations of the 15 inscribed World Heritage locations. Although we collectively treated the sites in each of the three longitudinal zones for our regression analysis, i.e., without differentiation among the single inscribed sites, the conducted causal analysis returned significant relationships among precipitation, elevation, and vegetation growth, and confirmed the previously reported influence of these two factors on TNPI Growth in both main and interaction terms. This was also the case for all sites located in the mid-zone, regardless of how much precipitation each site received. Finally, ∆LST was expectedly not an influential factor, mainly due to the higher latitudinal gradient (sea-mountain distance) and a significantly higher receipt of precipitation overall across the western zone. This was, however, shown to be different in the eastern zone, in which ∆LST was influential as both the main and interaction term due to the generally warmer climate, less precipitation, and lower regulatory effect introduced by CS [43].
All in all, the amount of precipitation did not prove to notably contribute to the direction of relationship between the growth rate and either elevation or ∆LST factors, confirming that HF stands are presumably still receiving enough humidity (both direct and indirect) to sustain their growth under current precipitation regimes. However, recent reports confirm emerging more frequent drought events within and at the margins of HF [48], which cause deficit of forest water content, leading to reduction in greenness and forest growth [49]. Temporal increase in intensity and duration of such events necessitates further remote sensing studies of their effects on forest growth indicators, in particular, using spatially higher resolution data than that of those applied in the mentioned former studies that made use of MODIS imagery. As such, our workflow can be further applied to test the hypothesis whether this effect can be shown by time series of high-resolution data, though the duration of currently available 10 m Sentinel-2 data and its limited availability due to cloud cover in many places is still a constraint for studying long-term drought events.

Conclusions
Our multi-resolution, multi-temporal analysis of MODIS and Sentinel-2 data was, to the best of our knowledge, the only study so far to combine freely available data, big data analysis, and modeling within a remote sensing-only analysis of forest growth across HF. Cautiously speaking, we also did not come across any other comparable study over the temperate forests of the northern hemisphere. As such, this analysis can be regarded as a benchmark for future analyses on comparably structured forest sites. The results suggested that SOS and EOS both change with environmental and topographic factors. They also showed that vegetation growth can be correlated with longitude (on larger spatial scale) and latitude (on finer spatial scale) that were generally in-line with former reports in terms of vegetation taxa and species associations. The results presented here can be further augmented as larger cloud-free composites of remote sensing. Climate data will also be available via GEE and open new horizons for both science and practice for rapid monitoring of vegetation dynamics across these constantly changing World Heritage ecosystems at no cost.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, Hooman Latifi (H.L.), upon reasonable request.