Characterisation of Functional-Trait Dynamics at High Spatial Resolution in a Mediterranean Forest from Sentinel-2 and Ground-Truth Data

: The characterisation of functional-trait dynamics of vegetation from remotely sensed data complements the structural characterisation of ecosystems. In this study we characterised for the ﬁrst time the spatial heterogeneity of the intra-annual dynamics of the fraction of absorbed photosynthetically active radiation (FAPAR) as a functional trait of the vegetation in Prades Mediterranean forest in Catalonia, Spain. FAPAR was derived from the Multispectral Instrument (MSI) on the Sentinel-2 satellite and validated by comparison with the ground measurements acquired in June 2017 at the annual peak of vegetation activity. The validation results showed that most of points were distributed along the 1:1 line, with no bias nor scattering: R 2 = 0.93, p < 0.05; with a root mean square error of 0.03 FAPAR (4.3%). We classiﬁed the study area into nine vegetation groups with different dynamics of FAPAR using a methodology that is objective and repeatable over time. This functional classiﬁcation based on the annual magnitude (FAPAR-M) and the seasonality (FAPAR-CV) from the data on one year (2016–2017) complements structural classiﬁcations. The internal heterogeneity of the FAPAR dynamics in each land-cover type is attributed to the environmental and to the speciﬁc species composition variability. A spatial autoregressive (SAR) model for the main type of land cover, evergreen holm oak forest ( Quercus ilex ), indicated that topographic aspect, slope, height, and the topographic aspect x slope interaction accounted for most of the spatial heterogeneity of the functional trait FAPAR-M, thus improving our understanding of the explanatory factors of the annual absorption of photosynthetically active radiation by the vegetation canopy for this ecosystem.


Introduction
Characterisation of the dynamics of vegetation functional traits complements structural characterisation of ecosystems [1]. These characterisations are of special interest in ecological studies of climate change, because functional traits tend to respond faster than structural traits to environmental changes [2,3].
The dynamics of functional traits of vegetation can be assessed using various synthetic variables of the annual activity (e.g., annual magnitude and seasonality) derived from remote sensing data at different spatial and temporal scales [1,[3][4][5][6][7][8][9][10][11]. Time series of satellite data are particularly useful for studying the temporal dynamics of the fraction of absorbed photosynthetically active radiation (FAPAR). The FAPAR corresponds to the fraction of photosynthetically active radiation (radiation at wavelengths of 0.4-0.7 µm) absorbed by the canopy. It depends on canopy structure, vegetation element optical properties and illumination conditions. FAPAR is one of the main components of primary productivity models [12]. It is a key variable for monitoring vegetation, studying the carbon cycle, and modelling climate. FAPAR has been recognized as an essential climate variable by the Global Climate Observing System (GCOS) [13].
The Multispectral Instrument (MSI) on the Sentinel-2 (S2) satellite (https://earth.esa.int/web/ sentinel) provides data with a high frequency of acquisition (5 days after S2A and S2B are in tandem and 10 days when only S2A was on orbit) and a high spatial resolution (10 m, 20 m or 60 m depending on the spectral band), which are adequate for the characterisation of the functional-trait dynamics of vegetation cover at a high level of detail. FAPAR can be processed on demand using the Sentinel-2 data service platform (https://s2.boku.eodc.eu) [14], which provides atmospherically corrected images and biophysical variables, including FAPAR, for any land surface on Earth.
Validation exercises are essential in remote-sensing studies for assessing the uncertainty associated to the satellite products. We directly validated the Sentinel-2 satellite products at our study site by comparison with ground-truth measurements [15,16]. Ground measurements acquisition and sampling design follow the international protocols of best practices standards for the validation of satellite products of the Committee on Earth Observation Satellites (CEOS) Working Group on Calibration and Validation (WGCV), Land Product Validation (LPV) subgroup [17], which are regularly applied to validate vegetation products generated by the Copernicus Land Monitoring Service [18].

Objectives
The main objective of this study was to characterise the spatial heterogeneity of FAPAR dynamics as a functional trait of a Mediterranean forest in Prades, Catalonia, using remote-sensing techniques and in situ measurements.
The specific objectives were to: I validate a high-spatial-resolution FAPAR product generated from the MSI data for a representative area of Prades forest with ground-truth measurements; II characterise the intra-annual FAPAR activity with two synthetic variables derived from FAPAR seasonal curves: annual magnitude (FAPAR-M) and seasonality (FAPAR-CV); III analyse the relationship between the topographic variables and the spatial heterogeneity of FAPAR-M of the main type of land cover at the study site: evergreen oak forest (Quercus ilex).

Materials and Methods
The study site is first described (Section 2.1). Then the databases (Section 2.2) and methods (Sections 2.3-2.5) are presented. Figure 1 shows a flow chart of the methodology.
Remote Sens. 2018, 10, x FOR PEER REVIEW 2 of 23 (FAPAR). The FAPAR corresponds to the fraction of photosynthetically active radiation (radiation at wavelengths of 0.4-0.7 µ m) absorbed by the canopy. It depends on canopy structure, vegetation element optical properties and illumination conditions. FAPAR is one of the main components of primary productivity models [12]. It is a key variable for monitoring vegetation, studying the carbon cycle, and modelling climate. FAPAR has been recognized as an essential climate variable by the Global Climate Observing System (GCOS) [13]. The Multispectral Instrument (MSI) on the Sentinel-2 (S2) satellite (https://earth.esa.int/web/ sentinel) provides data with a high frequency of acquisition (5 days after S2A and S2B are in tandem and 10 days when only S2A was on orbit) and a high spatial resolution (10 m, 20 m or 60 m depending on the spectral band), which are adequate for the characterisation of the functional-trait dynamics of vegetation cover at a high level of detail. FAPAR can be processed on demand using the Sentinel-2 data service platform (https://s2.boku.eodc.eu) [14], which provides atmospherically corrected images and biophysical variables, including FAPAR, for any land surface on Earth.
Validation exercises are essential in remote-sensing studies for assessing the uncertainty associated to the satellite products. We directly validated the Sentinel-2 satellite products at our study site by comparison with ground-truth measurements [15,16]. Ground measurements acquisition and sampling design follow the international protocols of best practices standards for the validation of satellite products of the Committee on Earth Observation Satellites (CEOS) Working Group on Calibration and Validation (WGCV), Land Product Validation (LPV) subgroup [17], which are regularly applied to validate vegetation products generated by the Copernicus Land Monitoring Service [18].

Objectives
The main objective of this study was to characterise the spatial heterogeneity of FAPAR dynamics as a functional trait of a Mediterranean forest in Prades, Catalonia, using remote-sensing techniques and in situ measurements.
The specific objectives were to: I. validate a high-spatial-resolution FAPAR product generated from the MSI data for a representative area of Prades forest with ground-truth measurements; II. characterise the intra-annual FAPAR activity with two synthetic variables derived from FAPAR seasonal curves: annual magnitude (FAPAR-M) and seasonality (FAPAR-CV); III. analyse the relationship between the topographic variables and the spatial heterogeneity of FAPAR-M of the main type of land cover at the study site: evergreen oak forest (Quercus ilex).

Materials and Methods
The study site is first described (Section 2.1). Then the databases (Section 2.2) and methods (Sections 2.3-2.5) are presented. Figure 1 shows a flow chart of the methodology.

Study Site
The study was conducted at a site 3.5 × 3.5 km in a Mediterranean forest (coordinates of the centre in the WGS84 geographic coordinate system: 41 • 20 26.76"N, 1 • 2 16.77"E) in Prades, Tarragona, Catalonia, Spain, in a mountainous area with altitudes ranging between 622 and 1138 m a.s.l. We defined the study area taking into account the representativeness of natural conditions of the forest and minimising anthropic environments. The main land-cover types at the study site are evergreen oak forest (Q. ilex), Scots pine forest (Pinus sylvestris), Austrian pine forest (P. nigra), and Pubescent oak forest (Q. pubescens) (Figure 2). Climatic conditions are typical of montane Mediterranean climates. The mean annual temperature is 12 • C and average annual rainfall is 658 mm. Summer drought lasts from mid-June to mid-September [19]. The mean annual temperature for the study period (15 October 2016 to 15 October 2017) was 15 • C, and the total cumulative rainfall was 600 mm, calculated using data from the State Meteorological Agency (AEMET) [20].

Study Site
The study was conducted at a site 3.5 × 3.5 km in a Mediterranean forest (coordinates of the centre in the WGS84 geographic coordinate system: 41°20′26.76′′N, 1°2′16.77′′E) in Prades, Tarragona, Catalonia, Spain, in a mountainous area with altitudes ranging between 622 and 1138 m a.s.l. We defined the study area taking into account the representativeness of natural conditions of the forest and minimising anthropic environments. The main land-cover types at the study site are evergreen oak forest (Q. ilex), Scots pine forest (Pinus sylvestris), Austrian pine forest (P. nigra), and Pubescent oak forest (Q. pubescens) (Figure 2). Climatic conditions are typical of montane Mediterranean climates. The mean annual temperature is 12 °C and average annual rainfall is 658 mm. Summer drought lasts from mid-June to mid-September [19]. The mean annual temperature for the study period (15 October 2016 to 15 October 2017) was 15 °C, and the total cumulative rainfall was 600 mm, calculated using data from the State Meteorological Agency (AEMET) [20]. Types and canopy coverages (cc) of the land covers in the study site. The percentage of the area of each land cover type over the entire study site is shown in detail. Data obtained from the Land Cover Map of Catalonia [21]. Reference System ETRS89 Universal Transverse Mercator (UTM) zone 31 N.

Satellite Products
We used time series of FAPAR derived from MSI data on board Sentinel-2 (https://earth.esa.int/web/sentinel). The FAPAR was processed on demand by the Sentinel-2 data service platform on the Earth Observation Data Centre (EODC) (https://s2.boku.eodc.eu) [14]. The platform processes European Space Agency's (ESA) Level-1C top-of-atmosphere reflectance using the Sen2Cor algorithm (version 2.4) to obtain atmospherically-corrected bottom-of-atmosphere (BoA) reflectance (Level-2A). The FAPAR products were generated using a neural network algorithm [22] calibrated on simulations from the coupled PROSPECT and SAIL radiative transfer models. The instantaneous FAPAR value at the time of the satellite overpass is computed directly from the radiative transfer model in the green parts of the canopy. The FAPAR products have a spatial resolution of 10 m and a temporal resolution of 10 d (5 d after the launch of the second twin satellite Sentinel-2). We checked the spatial accuracy of Sentinel-2 products using ground control points (GPS coordinates of bare soil areas, roads, and firebreaks), and we did not identify geolocation problems for our study site. Types and canopy coverages (cc) of the land covers in the study site. The percentage of the area of each land cover type over the entire study site is shown in detail. Data obtained from the Land Cover Map of Catalonia [21]. Reference System ETRS89 Universal Transverse Mercator (UTM) zone 31 N.

Satellite Products
We used time series of FAPAR derived from MSI data on board Sentinel-2 (https://earth.esa. int/web/sentinel). The FAPAR was processed on demand by the Sentinel-2 data service platform on the Earth Observation Data Centre (EODC) (https://s2.boku.eodc.eu) [14]. The platform processes European Space Agency's (ESA) Level-1C top-of-atmosphere reflectance using the Sen2Cor algorithm (version 2.4) to obtain atmospherically-corrected bottom-of-atmosphere (BoA) reflectance (Level-2A). The FAPAR products were generated using a neural network algorithm [22] calibrated on simulations from the coupled PROSPECT and SAIL radiative transfer models. The instantaneous FAPAR value at the time of the satellite overpass is computed directly from the radiative transfer model in the green parts of the canopy. The FAPAR products have a spatial resolution of 10 m and a temporal resolution of 10 d (5 d after the launch of the second twin satellite Sentinel-2). We checked the spatial accuracy of Sentinel-2 products using ground control points (GPS coordinates of bare soil areas, roads, and firebreaks), and we did not identify geolocation problems for our study site.
We identified and distinguished between cloud-free satellite products for one year (15 October 2016 to 15 October 2017) by visual interpretation of Sentinel-2 atmospherically corrected BoA reflectance (Level-2A) data. We generated true-and false-colour RGB images for each date to identify the shape, size, pattern, tone, texture, shadowing, and association of the elements. We evaluated all images available for the study period, ultimately selecting a total of 19 products for the construction of the seasonal curves.

Topographic Variables
We used a digital elevation model (DEM) with high spatial resolution (pixel size of 10 m) provided by the Data Services Platform for Sentinel-2 (http://s2.boku.eodc.eu/). We obtained the variable height and derived the variables slope, soil water distribution, topographic position and aspect for the main type of land cover in the study site (evergreen oak forest, Q. ilex) (Figure 3). of the seasonal curves.

Topographic Variables
We used a digital elevation model (DEM) with high spatial resolution (pixel size of 10 m) provided by the Data Services Platform for Sentinel-2 (http://s2.boku.eodc.eu/). We obtained the variable height and derived the variables slope, soil water distribution, topographic position and aspect for the main type of land cover in the study site (evergreen oak forest, Q. ilex) (Figure 3).
The soil water distribution, which is affected by topography, was calculated using the Topographic Wetness Index (TWI) [23], determined as: where a is the area of a specific catchment, the local upslope area draining through a certain point per unit contour length, and β is the local slope. TWI was calculated using the Multiple Flow Direction algorithm based on maximum downslope gradient (MFD-md) [24,25], available in the System for Automated Geoscientific Analyses (SAGA 2.3.2) programme.
The topographic position of the terrain was calculated using the Topographic Position Index (TPI), defined as the difference in height (positive or negative) between a central pixel and its surrounding cells within a predetermined radius (landscape unit considered) [26,27]. We predefined a landscape unit radius of 150 m to highlight minor valleys and ridges by the visual interpretation of very-high-spatial-resolution (VHR) satellite imagery using Google Earth.  The soil water distribution, which is affected by topography, was calculated using the Topographic Wetness Index (TWI) [23], determined as: where a is the area of a specific catchment, the local upslope area draining through a certain point per unit contour length, and β is the local slope. TWI was calculated using the Multiple Flow Direction algorithm based on maximum downslope gradient (MFD-md) [24,25], available in the System for Automated Geoscientific Analyses (SAGA 2.3.2) programme.
The topographic position of the terrain was calculated using the Topographic Position Index (TPI), defined as the difference in height (positive or negative) between a central pixel and its surrounding cells within a predetermined radius (landscape unit considered) [26,27]. We predefined a landscape unit radius of 150 m to highlight minor valleys and ridges by the visual interpretation of very-high-spatial-resolution (VHR) satellite imagery using Google Earth.

Land-Cover Types
The main land-cover types at the study site ( Figure 2) were analysed by the high-resolution thematic cartography of the Land Cover Map of Catalonia and ground measurements. Appendix A provides further details and descriptions of the main land-cover types in the study area.
We characterised the intra-annual dynamics of FAPAR by type of land cover using the shapefile for the Land Cover Map of Catalonia (level 5 of the legend) [21] and calculated the areas using the ETRS89 UTM zone 31 N reference system.

Ground Measurements and Data Processing
The satellite products (see Section 2.2.1) were validated following the international protocols of best practices standards of the CEOS WGCV LPV subgroup [17].
Ground-truth data were acquired using digital hemispherical photography (DHP), which allowed the calculation of FAPAR, and other variables of canopy architecture (e.g., Leaf Area Index (LAI), Fraction of green vegetation cover (FCOVER)), based on directional measurements of the fraction of gaps in the vegetation cover [28][29][30][31][32]. We used a digital Nikon camera and an extreme wide-angle lens (180 • fisheye lens) for each measurement.
We used CAN-EYE v6.4 software (https://www6.paca.inra.fr/can-eye) for processing the RGB hemispherical photographs and extracting the FAPAR canopy values. We defined the characteristics of the camera + fisheye lens system (optical centre at the centre of the image, equidistant projection function) and all information required for the software [33] for the processing. We also masked the field of view of the lens to values within the range 0-60 • to ensure that the area captured was within the sampling unit. The software allowed us to extract the FAPAR value for each sampling point of the study site by a supervised-classification process using the green elements of the RGB digital images. The CAN-EYE FAPAR values [33,34], corresponded to the direct instantaneous FAPAR, "Black Sky FAPAR", obtained from the fraction of intercepted photosynthetically active radiation as 1-gap fraction at solar zenith corresponding to 10:30 (Sentinel-2 revisit time for the study site: 10:50) (https://sentinel.esa.int/web/sentinel/missions/sentinel-2/acquisition-plans).
We acquired two sets of hemispherical photographs for sampling units with over and understory: upward looking (overstory) and downward looking (understory). The two sets of acquisitions were processed separately to derive FAPAR [16,35]. FAPAR was calculated assuming the independence of the gaps inside the understory and inside the overstory. The biophysical variable was calculated as:

Sampling Scheme
We designed the sampling scheme following the recommendations for running a field campaign of the FP7 ImagineS project (http://fp7-imagines.eu), consistent with the Validation of Land European Remote sensing Instruments (VALERI) project (http://w3.avignon.inra.fr/valeri/) [15].
We first surveyed the study area to identify the types of vegetation cover, determine the accessibility of the terrain, and plan the logistics of collecting the field data. We then programmed the date for the ground measurements to 15 June 2017 taking into account: (a) the period of maximum vegetation greenness for the study site (analysed with time series of the Normalized Difference Vegetation Index (NDVI)) to minimise problems due to the presence of non-green vegetation elements [15], (b) the acquisition plans for Sentinel-2 (https://sentinel.esa.int/web/sentinel/missions/ sentinel-2/acquisition-plans), and (c) the weather forecast, increasing the probability for clear skies.
The elementary sampling unit (ESU) consisted of a 10 × 10 m plot to match the pixel size of the Sentinel-2 satellite biophysical product to be validated. We acquired 8-15 measurements (DHP) for each ESU to obtain a reasonable estimate of the gap fraction [29,33]. We followed a predefined sampling scheme to locate the photos within the ESU as independently as possible ( Figure 4). from the borders of roads and surrounded by pixels with the same type of vegetation, because the radiative transfer of a single pixel may be influenced by the neighbouring pixels by lateral fluxes [36]. The centre of each ESU was geolocated using a global positioning system receiver.
Although our sampling took into account the main land covers and the environmental heterogeneity of the site, only a fraction of the FAPAR range of the study site ( Figure 5) was finally represented by the ESUs due to the difficulty of accessing the terrain. The FAPAR values sampled in situ were mostly between 0.58 and 0.90 (range 0.32), with a mean of 0.69 and coefficient of variation of 13%. Most of the ESUs were evergreen oak forest (15 of the 28 sampling units), followed by Scots pine (5), Maritime pine (3), Pyrenean oak (2), Spanish chestnut (2), and Pubescent oak (1) forests ( Figure 5). Finally, due to the lack of ESUs with FAPAR values lower than 0.58, we indirectly geolocated by visual interpretation of VHR satellite imagery using Google Earth two additional units of bare soil where we assume zero values for FAPAR. Appendix B provides further details of the ESU geolocation and FAPAR values.  We installed 28 ESUs for six forest classes to represent the main vegetation cover types and FAPAR values at the study site ( Figure 4). We established the ESUs at distances higher than 50 m from the borders of roads and surrounded by pixels with the same type of vegetation, because the radiative transfer of a single pixel may be influenced by the neighbouring pixels by lateral fluxes [36]. The centre of each ESU was geolocated using a global positioning system receiver.
Although our sampling took into account the main land covers and the environmental heterogeneity of the site, only a fraction of the FAPAR range of the study site ( Figure 5

Validation
We performed linear regressions to estimate the strength of the relationship between the ground-truth FAPAR values and the Sentinel-2 FAPAR product values. We assume no temporal variation in FAPAR within the three days between satellite acquisition on 12 June and the ground measurements on 15 June. We compared the slope and intercept with the 1:1 line [37]. We also calculated the root mean square error (RMSE) to represent the standard deviation of the differences between the values of the Sentinel-2 FAPAR product and ground-truth FAPAR to complement the validation. RMSE was calculated as: where Xin situ is the ground-truth FAPAR value and Xsentinel2 is the Sentinel-2 FAPAR product value.

Analysis and Characterisation of Functional-Trait Dynamics
We used one-year of cloud-free time series of FAPAR (see Section 2.2.1) to describe, classify, and map the intra-annual functional activity of the vegetation cover considering existing methodologies for functional characterisations of ecosystems [1,[4][5][6][7][8][9]11].
We calculated eight semi-season composites (late autumn, early winter, late winter, early spring, late spring, early summer, late summer, and early autumn) as the average of the satellite data

Validation
We performed linear regressions to estimate the strength of the relationship between the ground-truth FAPAR values and the Sentinel-2 FAPAR product values. We assume no temporal variation in FAPAR within the three days between satellite acquisition on 12 June and the ground measurements on 15 June. We compared the slope and intercept with the 1:1 line [37]. We also calculated the root mean square error (RMSE) to represent the standard deviation of the differences between the values of the Sentinel-2 FAPAR product and ground-truth FAPAR to complement the validation. RMSE was calculated as: where X in situ is the ground-truth FAPAR value and X sentinel2 is the Sentinel-2 FAPAR product value.

Analysis and Characterisation of Functional-Trait Dynamics
We used one-year of cloud-free time series of FAPAR (see Section 2.2.1) to describe, classify, and map the intra-annual functional activity of the vegetation cover considering existing methodologies for functional characterisations of ecosystems [1,[4][5][6][7][8][9]11].
We calculated eight semi-season composites (late autumn, early winter, late winter, early spring, late spring, early summer, late summer, and early autumn) as the average of the satellite data available for each period (Table 1). We performed linear interpolations for late autumn and early winter because of: (i) the lack of cloud-free images for these semi-seasons; and (ii) the detection of contaminated pixels by terrain shadows, due to the lower solar angle in these periods. Linear interpolations were performed between early autumn and late winter values in the QGIS 2.18.14 environment. We generated one FAPAR seasonal curve with the semi-season composites for each pixel of the study area ( Figure 6) and calculated two synthetic variables from the seasonal curves of FAPAR that capture important traits of the FAPAR dynamics: FAPAR-M as the mean of the eight semi-season composites and FAPAR-CV as the ratio of the standard deviation to the mean. available for each period (Table 1). We performed linear interpolations for late autumn and early winter because of: (i) the lack of cloud-free images for these semi-seasons; and (ii) the detection of contaminated pixels by terrain shadows, due to the lower solar angle in these periods. Linear interpolations were performed between early autumn and late winter values in the QGIS 2.18.14 environment. We generated one FAPAR seasonal curve with the semi-season composites for each pixel of the study area ( Figure 6) and calculated two synthetic variables from the seasonal curves of FAPAR that capture important traits of the FAPAR dynamics: FAPAR-M as the mean of the eight semi-season composites and FAPAR-CV as the ratio of the standard deviation to the mean. We analysed the spatial heterogeneity of FAPAR-M and FAPAR-CV over the entire study area and per land-cover type (see Section 2.2.3). We identified vegetation groups with different dynamics of FAPAR. The groups were defined using fixed limits between classes to allow future inter-annual comparisons, capturing the year-to-year dynamics [11]. The range of variation of FAPAR-M and FAPAR-CV were divided into three classes (low, medium, and high), fixing intervals at 1/3 of the ranges within 2.5 and 97.5 percentiles of the distributions to avoid long tails in the histograms. The combination of the generated classes (3 × 3) allowed us to identify, in a single map, the pixels for the vegetation groups, based on different dynamics of radiation absorption. We assigned a two-number code to each vegetation group. The first number corresponded to the FAPAR-M class, and the second number corresponded to the FAPAR-CV class. The coding in both cases was 1 for the low class, 2 for the medium and 3 for the high. We analysed the spatial patterns of the vegetation groups generated and evaluated the correspondence between them and the main vegetation cover types of the study area. We analysed the spatial heterogeneity of FAPAR-M and FAPAR-CV over the entire study area and per land-cover type (see Section 2.2.3). We identified vegetation groups with different dynamics of FAPAR. The groups were defined using fixed limits between classes to allow future inter-annual comparisons, capturing the year-to-year dynamics [11]. The range of variation of FAPAR-M and FAPAR-CV were divided into three classes (low, medium, and high), fixing intervals at 1/3 of the ranges within 2.5 and 97.5 percentiles of the distributions to avoid long tails in the histograms. The combination of the generated classes (3 × 3) allowed us to identify, in a single map, the pixels for the vegetation groups, based on different dynamics of radiation absorption. We assigned a two-number code to each vegetation group. The first number corresponded to the FAPAR-M class, and the second number corresponded to the FAPAR-CV class. The coding in both cases was 1 for the low class, 2 for the medium and 3 for the high. We analysed the spatial patterns of the vegetation groups generated and evaluated the correspondence between them and the main vegetation cover types of the study area.

Explanatory Topographic Factors of FAPAR-M Spatial Variability
We analysed the relationship between the topographic variables (see Section 2.2.2) and the spatial heterogeneity of FAPAR-M for the main land cover class in the study area, that is, evergreen oak forest (Q. ilex) using spatial autoregressive (SAR) models.
We compiled a spatially explicit database integrating FAPAR-M as the dependent variable and height (DEM), slope, TWI, TPI, and aspect of the terrain (derived from the DEM) as independent variables. We thus integrated all variables in a geographic information system (GIS). We re-sampled the spatial resolution of the satellite imagery to 50-m cells to minimise the spatial autocorrelation of the data and considered for the models only cells with evergreen oak forest (Q. ilex) cover (cells number = 2246). We obtained the average FAPAR-M and the height, slope, TWI, TPI, and aspect of the terrain for each 50-m cell. The topographic aspect data, whose original units were azimuth angles, were transformed into a four-level categorical variable (north, south, east, west). The collinearity of the explanatory topographic factors were inspected using Pearson correlation coefficients and generalised variance-inflation factors (GVIFs) [38,39], using the R 'car' package [40]. Appendix C provides further details of the collinearity analysis.
All models were fitted using the R 'spdep' package and the lagsarlm function [41], with the queen criterion on the matrix of spatial weights [42]. The normality, linearity, and homoscedasticity of the residuals were inspected visually. We applied model selection to include those terms that minimised the value of the Akaike Information Criterion (AIC) [43]. Statistical analysis was performed in R 3.4.2.

Validation of Satellite Products
The validation results (Figure 7) showed that the Sentinel-2 FAPAR product at the study site agreed well with the ground-truth data (n = 28) with an RMSE of 0.03 (4.3% of the mean), and RMSE of 0.031 (4.8% of the mean) with the two additional bare soil sampling units (n = 30). The data points were distributed along the 1:1 line, with no bias or scattering. The coefficients of the regression were close to one: R 2 = 0.93 (p < 0.05) for the n = 28 ESUs, and R 2 = 0.98 (p < 0.05) including the two additional sampling units (n = 30).

Explanatory Topographic Factors of FAPAR-M Spatial Variability
We analysed the relationship between the topographic variables (see Section 2.2.2) and the spatial heterogeneity of FAPAR-M for the main land cover class in the study area, that is, evergreen oak forest (Q. ilex) using spatial autoregressive (SAR) models.
We compiled a spatially explicit database integrating FAPAR-M as the dependent variable and height (DEM), slope, TWI, TPI, and aspect of the terrain (derived from the DEM) as independent variables. We thus integrated all variables in a geographic information system (GIS). We re-sampled the spatial resolution of the satellite imagery to 50-m cells to minimise the spatial autocorrelation of the data and considered for the models only cells with evergreen oak forest (Q. ilex) cover (cells number = 2246). We obtained the average FAPAR-M and the height, slope, TWI, TPI, and aspect of the terrain for each 50-m cell. The topographic aspect data, whose original units were azimuth angles, were transformed into a four-level categorical variable (north, south, east, west). The collinearity of the explanatory topographic factors were inspected using Pearson correlation coefficients and generalised variance-inflation factors (GVIFs) [38,39], using the R 'car' package [40]. Appendix C provides further details of the collinearity analysis.
All models were fitted using the R 'spdep' package and the lagsarlm function [41], with the queen criterion on the matrix of spatial weights [42]. The normality, linearity, and homoscedasticity of the residuals were inspected visually. We applied model selection to include those terms that minimised the value of the Akaike Information Criterion (AIC) [43]. Statistical analysis was performed in R 3.4.2.

Validation of Satellite Products
The validation results (Figure 7) showed that the Sentinel-2 FAPAR product at the study site agreed well with the ground-truth data (n = 28) with an RMSE of 0.03 (4.3% of the mean), and RMSE of 0.031 (4.8% of the mean) with the two additional bare soil sampling units (n = 30). The data points were distributed along the 1:1 line, with no bias or scattering. The coefficients of the regression were close to one: R 2 = 0.93 (p < 0.05) for the n = 28 ESUs, and R 2 = 0.98 (p < 0.05) including the two additional sampling units (n = 30).

Analysis and Characterisation of Functional-Trait Dynamics
FAPAR-M had a mean of 0.61 and a spatial standard deviation of 0.09. FAPAR-M ranged between 0.05 and 0.86 (range 0.81), with a negatively skewed distribution (Figure 8), where the values between 0.63 and 0.67 were the most frequent. Values were extreme for rocky outcrops and deciduous trees other than those of our land-cover types (≥20% canopy coverage (cc)), respectively.

Analysis and Characterisation of Functional-Trait Dynamics
FAPAR-M had a mean of 0.61 and a spatial standard deviation of 0.09. FAPAR-M ranged between 0.05 and 0.86 (range 0.81), with a negatively skewed distribution (Figure 8), where the values between 0.63 and 0.67 were the most frequent. Values were extreme for rocky outcrops and deciduous trees other than those of our land-cover types (≥20% canopy coverage (cc)), respectively.
The geographically explicit data set for FAPAR-M allowed observation of the spatial distribution of the data range (Figure 8), which identified specific contrasting patterns in the pixel groupings across the study site. FAPAR-CV at the study site had a mean of 0.06 and a spatial standard deviation of 0.04. FAPAR-CV ranged between 0.01 and 0.72 (range 0.71), with a positively skewed distribution (Figure 9), where the values between 0.03 and 0.04 were the most frequent. Values were extreme for Maritime pine forest (≥20% cc) and Pyrenean oak forest (≥20% cc), respectively.
The geographically explicit data set for FAPAR-CV allowed observation of the spatial distribution of the data range (Figure 9), which identified specific contrasting patterns in the pixel groupings across the study site, as with FAPAR-M. The geographically explicit data set for FAPAR-M allowed observation of the spatial distribution of the data range (Figure 8), which identified specific contrasting patterns in the pixel groupings across the study site.
FAPAR-CV at the study site had a mean of 0.06 and a spatial standard deviation of 0.04. FAPAR-CV ranged between 0.01 and 0.72 (range 0.71), with a positively skewed distribution (Figure 9), where the values between 0.03 and 0.04 were the most frequent. Values were extreme for Maritime pine forest (≥20% cc) and Pyrenean oak forest (≥20% cc), respectively.
The geographically explicit data set for FAPAR-CV allowed observation of the spatial distribution of the data range (Figure 9), which identified specific contrasting patterns in the pixel groupings across the study site, as with FAPAR-M. The integration of the FAPAR seasonal curves from the geographically explicit data set with data from the Land Cover Map of Catalonia allowed us to analyse the FAPAR intra-annual dynamics of the main types of land cover ( Figure 10). The integration of the FAPAR seasonal curves from the geographically explicit data set with data from the Land Cover Map of Catalonia allowed us to analyse the FAPAR intra-annual dynamics of the main types of land cover ( Figure 10). The integration of the FAPAR seasonal curves from the geographically explicit data set with data from the Land Cover Map of Catalonia allowed us to analyse the FAPAR intra-annual dynamics of the main types of land cover ( Figure 10). The main type of land cover, evergreen oak forest (≥20% cc), had a moderate FAPAR-M (mean = 0.619), with a spatial heterogeneity (calculated using the spatial standard deviation) of 0.085 (13.7% The main type of land cover, evergreen oak forest (≥20% cc), had a moderate FAPAR-M (mean = 0.619), with a spatial heterogeneity (calculated using the spatial standard deviation) of 0.085 (13.7% of the mean). Other deciduous trees had the highest FAPAR-M (mean = 0.715), and evergreen oak forest (5-20% cc), shrubland, and scree had the lowest values (mean FAPAR-M = 0.47, 0.48, 0.489, respectively). The remaining types of land cover had intermediate values of FAPAR-M ( Figure 11). The spatial heterogeneity of FAPAR-M was highest for scree, shrubland, and evergreen oak forest (5-20% cc), with spatial standard deviations of 0.13, 0.09, and 0.08 and means of 26, 18, and 18%, respectively. The other types of land cover had spatial heterogeneities ranging from moderate to low (Figure 11).  Figure 11). The spatial heterogeneity of FAPAR-M was highest for scree, shrubland, and evergreen oak forest (5-20% cc), with spatial standard deviations of 0.13, 0.09, and 0.08 and means of 26, 18, and 18%, respectively. The other types of land cover had spatial heterogeneities ranging from moderate to low ( Figure 11). Evergreen oak forest (≥20% cc) had the second lowest intra-annual temporal variability (mean FAPAR-CV = 0.052), with a spatial heterogeneity (calculated using the spatial standard deviation) of 0.03 (mean = 56.5%). Pyrenean oak forest had the highest seasonal variability (mean FAPAR-CV = 0.154), and Maritime pine forest had the lowest (mean FAPAR-CV = 0.041). FAPAR-CV had higher values of spatial heterogeneity than FAPAR-M for all the types of land cover ( Figure 12).  Evergreen oak forest (≥20% cc) had the second lowest intra-annual temporal variability (mean FAPAR-CV = 0.052), with a spatial heterogeneity (calculated using the spatial standard deviation) of 0.03 (mean = 56.5%). Pyrenean oak forest had the highest seasonal variability (mean FAPAR-CV = 0.154), and Maritime pine forest had the lowest (mean FAPAR-CV = 0.041). FAPAR-CV had higher values of spatial heterogeneity than FAPAR-M for all the types of land cover ( Figure 12).  Figure 11). The spatial heterogeneity of FAPAR-M was highest for scree, shrubland, and evergreen oak forest (5-20% cc), with spatial standard deviations of 0.13, 0.09, and 0.08 and means of 26, 18, and 18%, respectively. The other types of land cover had spatial heterogeneities ranging from moderate to low ( Figure 11). Evergreen oak forest (≥20% cc) had the second lowest intra-annual temporal variability (mean FAPAR-CV = 0.052), with a spatial heterogeneity (calculated using the spatial standard deviation) of 0.03 (mean = 56.5%). Pyrenean oak forest had the highest seasonal variability (mean FAPAR-CV = 0.154), and Maritime pine forest had the lowest (mean FAPAR-CV = 0.041). FAPAR-CV had higher values of spatial heterogeneity than FAPAR-M for all the types of land cover ( Figure 12).  The categorisation of the gradients of FAPAR-M and FAPAR-CV and their subsequent combination allowed us to identify nine vegetation groups with different dynamics of absorption of photosynthetically active radiation ( Figure 13). The map shows the most frequent patterns of the FAPAR intra-annual dynamics for the study site (2016-2017): low seasonality with medium and high magnitudes of annual FAPAR (group codes 2-1 and 3-1, respectively). The map highlights the clear and contrasting patterns across the study site, although the remaining groups had lower percentages of occupation (Table 2). The categorisation of the gradients of FAPAR-M and FAPAR-CV and their subsequent combination allowed us to identify nine vegetation groups with different dynamics of absorption of photosynthetically active radiation ( Figure 13). The map shows the most frequent patterns of the FAPAR intra-annual dynamics for the study site (2016-2017): low seasonality with medium and high magnitudes of annual FAPAR (group codes 2-1 and 3-1, respectively). The map highlights the clear and contrasting patterns across the study site, although the remaining groups had lower percentages of occupation (Table 2).  Finally, we evaluated the correspondence between the vegetation groups with different dynamics of absorption of photosynthetically active radiation and the land-cover types across the study site (Table 3). This analysis did not evaluate the classification of functional-trait dynamics, but objectively evaluated the degree of correspondence with the structural classification for descriptive purposes.  Table 3. Correspondence between the vegetation groups with different dynamics of absorption of photosynthetically active radiation and the land-cover types across the study site. The following data are indicated for each group code: land-cover type, the total area in hectares occupied by each type (in parentheses, the percentage of the type area relative to the total surface type area of the study site), and the percentage of correspondence (proportion of all the area covered by a vegetation group that corresponds to a given land-cover type). Only the main types that cover the vegetation groups are included.

Explanatory Topographic Factors of FAPAR-M Spatial Variability
The SAR model that minimised the value of the AIC incorporated the topographic aspect, slope, height, and the topographic aspect: slope interaction as explanatory factors of FAPAR-M for the main vegetation cover in the study area: evergreen oak forest. The value of the spatial autoregressive parameter in the model Rho = 0.82 (p < 2.22 × 10 −16 ) indicated a high spatial autocorrelation of the data. The Lagrange Multiplier Test indicated that the autocorrelation of the residues was not significant (test value = 0.89; p = 0.34). Table 4 shows the parameters from the application of the SAR model, with the statistical significance for the model estimates. Figure 14 shows the density scatter-plots between FAPAR-M and the explanatory topographic variables, analysed for the topographic aspects. The following data are indicated for each group code: land-cover type, the total area in hectares occupied by each type (in parentheses, the percentage of the type area relative to the total surface type area of the study site), and the percentage of correspondence (proportion of all the area covered by a vegetation group that corresponds to a given land-cover type). Only the main types that cover the vegetation groups are included.

Explanatory Topographic Factors of FAPAR-M Spatial Variability
The SAR model that minimised the value of the AIC incorporated the topographic aspect, slope, height, and the topographic aspect: slope interaction as explanatory factors of FAPAR-M for the main vegetation cover in the study area: evergreen oak forest. The value of the spatial autoregressive parameter in the model Rho = 0.82 (p < 2.22 × 10 −16 ) indicated a high spatial autocorrelation of the data. The Lagrange Multiplier Test indicated that the autocorrelation of the residues was not significant (test value = 0.89; p = 0.34). Table 4 shows the parameters from the application of the SAR model, with the statistical significance for the model estimates. Figure 14 shows the density scatter-plots between FAPAR-M and the explanatory topographic variables, analysed for the topographic aspects.  Finally, complementary analyses with general linear models (GLMs) using data from random points obtained by spatial samplings of the images (n = 200, minimum distance between points of 100 m to reduce the spatial autocorrelation of the spatial data) also indicated the topographic aspect, slope, height, and the topographic aspect: slope interaction as explanatory factors of FAPAR-M spatial variability, with a coefficient of determination, R 2 , of 0.35 (p < 7.58 × 10 −15 ). The normality and homoscedasticity of the residuals for the complementary GLM model were also visually inspected.

Validation of Satellite Products
We defined the study area and geolocated specific coordinates, taking into account the representativeness of natural conditions of Mediterranean forests to minimise anthropic conditions, for facilitating future studies of the relationships between the inter-annual dynamics of the seasonal curves and environmental changes.
The ground measurements represented only a fraction of the range of FAPAR values at the study site due to the difficulty of accessing the terrain. The absence of measurements with low FAPAR values was indirectly countered by the geolocation of the two additional units of bare soil by the visual interpretation of VHR satellite imagery using Google Earth. FAPAR values lower than 0.58 could not be sampled and two bare soil points with FAPAR = 0 where added.
The validation results of the Sentinel-2 FAPAR product for the study site and sampling period agreed well with the ground-truth data. Multi-temporal campaigns may be required to assess the temporal accuracy of the satellite products for the study site, especially for late autumn and early winter, but are not strictly mandatory when validating biophysical products [15].
We highlight the possibility of using the ground-truth database, as derivations of the validation, to validate FAPAR global products, with lower spatial resolution but higher revisit frequency (e.g., PROBA-V and MODIS). Up-scaling ESU ground-measurement techniques (e.g., following the Validation of Land European Remote Sensing Instrument, VALERI guidelines, http://w3.avignon.inra. fr/valeri/) could increase the set of remote sensing products validated for the study site, strengthening the possibilities of monitoring these ecosystems at different spatial and temporal scales.

Analysis and Characterisation of Functional-Trait Dynamics
The characterisation of ecosystems based on functional-trait dynamics derived from remotely sensed data complements structural studies for describing the heterogeneity of ecosystems [1,3]. Several studies have been developed for characterising functional activity of vegetation cover at a regional scale [1,[4][5][6][7][8][9]11]. Most of these studies covered large regional scales, with remote sensors of medium and low spatial resolution (e.g., MODIS and NOAA/AVHRR), where seasonal dynamics were characterised using the NDVI or the Enhanced Vegetation Index (EVI) as proxies of the FAPAR. We applied these approaches and remote-sensing techniques to characterise the functional-traits dynamics of the vegetation cover for a representative area of the Prades Mediterranean ecosystem, at a smaller landscape scale, higher level of detail (10 × 10 m pixels) and using FAPAR seasonal curves directly. The construction of a temporal data series from the MSI data [14] allowed us to obtain highly detailed information for the FAPAR intra-annual activity over the study site and consequently for the dynamics of potential primary productivity [12], a key aspect of ecosystem functioning [44,45].
The methodology to characterise the functional-trait dynamics of the vegetation cover at high spatial resolution in Prades Mediterranean Forest is objective and repeatable over time using the same set of variables. Each pixel of the remotely sensed data was defined using variables with a clear biological meaning and analysed from their individual seasonal behaviours.
The remotely sensed data used for the construction of the seasonal curves corresponded to the period 2016/2017. The FAPAR dynamics therefore did not correspond to the mean dynamics of the study site but characterised the behaviour of the vegetation activity for that period. Future studies will be able to characterise the mean curves for each pixel of these ecosystems or will be able to design monitoring protocols to analyse their behaviours over time.
We performed linear interpolations for the late autumn and early winter FAPAR products for constructing the seasonal curves, due to the lack of satellite images, in addition to the detection of contaminated pixels by terrain shadows and the uncertainty of the accuracy of the images for these periods. The low solar angle at the revisit times for Sentinel-2 in these periods generated terrain shadows in areas of steep topographic slopes that could substantially decrease the FAPAR values.
We did not attempt to assimilate the classes of vegetation groups to particular land-cover types, because different types of vegetation covers may have similar FAPAR seasonal dynamics, and the same type may have different patterns of FAPAR annual activity [1]. This internal heterogeneity in each land-cover type can be attributed to the environmental variability of the study area (e.g., soil depth, soil fertility, and water availability) and to the specific species composition in each pixel. The land-cover types used in this study for characterising functional-trait dynamics corresponded to a classification of the main plant species at each site (Appendix A). Pixels with a main forest species mixed with one or more secondary forest species, and the characteristics and dynamics of the understory, can contribute substantially to the heterogeneity of FAPAR dynamics. The FAPAR seasonal curves for each land-cover type may therefore differ from the typical curves expected for a particular plant species, because they characterise ecosystem vegetation activity.

Explanatory Topographic Factors of FAPAR-M Spatial Variability
The spatial heterogeneity of the FAPAR-M of the main vegetation cover of the study area, evergreen oak forest, was due mainly to the topographic aspect, slope, height, and the topographic aspect: slope interaction. A complementary analysis with GLMs, using data from random points obtained by spatial samplings of the images, also identified topographic aspect, slope, height, and the topographic aspect: slope interaction as the most important explanatory topographic variables of FAPAR-M heterogeneity, explaining 35% of the spatial heterogeneity.
Topographic aspect had the lowest FAPAR-M values in the south-facing slopes and the highest values in the north-facing slopes, as expected due to the higher solar radiation received on the southern orientations at these latitudes, with effects on evapotranspiration and storage of soil moisture. The range of FAPAR-M values, however, was highly dispersed for each aspect. Height and topographic slope, despite being selected in the explanatory SAR models with a weak and negative slope, had highly scattered data points. FAPAR-M decreased greatly with the topographic slope for south-facing slopes, suggesting an increasing soil moisture limiting factor.
The percentage of FAPAR-M variability not accounted for may have been due to other environmental factors not included in the present study, such as soil depth. Soil depth is associated with soil fertility, water availability, and the development of root systems and so could be relevant to the heterogeneity of FAPAR-M in these ecosystems.
The spatial heterogeneity of FAPAR-M was not correlated with TWI or TPI at the study site, contrary to the expected relationship with water availability. The lack of correlation of FAPAR-M with TWI and TPI could be due to diverse factors. First, the ability of TWI and TPI to explain spatial patterns of water availability are limited, because water availability does not depend only on the topography (e.g., it is also influenced by factors such as radiation, physical and chemical properties of the soil, and microrelief conditions not captured by the spatial resolution of the DEM [46,47]). Second, we used a predefined TPI landscape unit of 150 m to highlight minor valleys and ridges and used the MFD-md algorithm, because it had the best significant, albeit low, correlation with the data for in situ soil moisture in a study conducted in a mountainous system in Poland [25], but it is highlighted that the TPI landscape units considered, or the TWI flow algorithms selected, could generate different results [25,[48][49][50].
The spatial autocorrelation of the data was not completely eliminated, despite resampling the dependent and independent variables to cells of lower spatial resolution. An increase in the size of the cells per unit area could help to reduce the spatial autocorrelation of the data, but with the potential loss of detail in the relationships between FAPAR-M and the topographic variables.
We aimed to characterise the relationships between the annual activity of FAPAR (though FAPAR-M) and topographic variables. The intra-annual temporal variation of the spatial heterogeneity and the specific relationships between FAPAR and the topographic factors for each semi-season were not analysed in this study. Incorporating specific analyses for each semi-season could improve the characterisation of functional-trait dynamics and our understanding of the explanatory factors of the FAPAR heterogeneity through the year.

Conclusions
This study is the first to characterise the spatial heterogeneity of the intra-annual dynamics of the fraction of the photosynthetically active radiation absorbed by the vegetation canopy for a portion of a Mediterranean forest in Prades, Tarragona, using high-spatial-resolution data from the Sentinel-2 satellite. The methodology used was objective and repeatable over time using the same set of variables to enable comparisons with future characterisations.
We validated the FAPAR satellite products with ground-truth data for the period of maximum greenness of the vegetation, with very good agreement (RMSE of 0.03 FAPAR). Conducting multi-temporal field campaigns, however, would be appropriate to ensure the seasonal stability of the satellite products given the characteristics of the topographic slopes on the study area, mainly for periods characterised by terrain shadows due to lower solar angles that could substantially decrease the FAPAR values. Alternative methodologies to the DHP should be explored for these periods, due to the potential higher interference of non-green elements in the digital images.
Slope, height, and the topographic aspect: slope interaction accounted for most of the heterogeneity of FAPAR-M for the main vegetation cover. Some of the heterogeneity not accounted for may have been due to other topographic factors not included in the study, such as soil depth. We recommend that incorporating soil depth into the models could be relevant for improving our understanding of the explanatory factors of the spatial heterogeneity of the annual magnitude of the absorption of photosynthetically active radiation by the vegetation canopy for these ecosystems.
Author Contributions: S.S. and A.V. conceived the research approach and designed the experiment with contributions from I.F. and J.P. The data was collected and processed by S.S. and A.V., S.S., A.V., I.F. and J.P. contributed to the analysis and discussion of results. The manuscript was written by S.S. with contributions from A.V., I.F. and J.P.  A   Table A1. Description of the categories of the land-cover map of Catalonia [21] (level 5 of the legend) for the study site.

Land-Cover Type (Canopy Coverage) Description
Evergreen oak forest (≥20%) Natural forests with a tree canopy consisting mainly of Quercus ilex (evergreen oak) with a tree density >20%.
Scots pine forest (≥20%) Natural forests with a tree canopy consisting mainly of Pinus sylvestris (Scots pine) with a tree density >20%.
Austrian pine forest (≥20%) Natural forests with a tree canopy consisting mainly of Pinus nigra (Austrian pine) with a tree density >20%.
Pubescent oak forest (≥20%) Natural forests with a tree canopy consisting mainly of Quercus pubescens (Pubescent oak) with a tree density >20%.
Evergreen oak forest (5-20%) Natural forests with a tree canopy consisting mainly of Quercus ilex (evergreen oak) with a tree density of 5-20%.

Scree
Surface where boulders and sharp pebbles have accumulated on the flanks of a mountain. The area is dry with sparse vegetation.

Shrubland
Formations with a substantial covering of shrubs or shrubby trees, when tree cover is <5%.
Maritime pine forest (≥20%) Natural forests with a tree canopy consisting mainly of Pinus pinaster (maritime pine) with a tree density >20%.
Other deciduous trees (≥20%) Natural forests with a tree canopy consisting mainly of deciduous trees (e.g., northern red oak, maple, nettle tree) with a tree density >20%.
Spanish chestnut forest (≥20%) Natural forests with a tree canopy consisting mainly of Castanea sativa (Spanish chestnut) with a tree density >20%. The following data are indicated for each elementary sampling unit (ESU): latitude and longitude coordinates in the WGS84 Geographic Coordinate System, the type of land cover (percent canopy coverage), the number of digital hemispherical photographs (DHP) acquired, and the ground-truth values of the fraction of absorbed photosynthetically active radiation (FAPAR) obtained. ESUs 29 and 30 correspond to the extra units of bare soil obtained from very-high-resolution satellite imagery.