Analyzing Daily Estimation of Forest Gross Primary Production Based on Harmonized Landsat-8 and Sentinel-2 Product Using SCOPE Process-Based Model

Vegetation top-of-canopy reflectance contains valuable information for estimating vegetation biochemical and structural properties, and canopy photosynthesis (gross primary production (GPP)). Satellite images allow studying temporal variations in vegetation properties and photosynthesis. The National Aeronautics and Space Administration (NASA) has produced a harmonized Landsat-8 and Sentinel-2 (HLS) data set to improve temporal coverage. In this study, we aimed to explore the potential and investigate the information content of the HLS data set using the Soil Canopy Observation of Photosynthesis and Energy fluxes (SCOPE) model to retrieve the temporal variations in vegetation properties, followed by the GPP simulations during the 2016 growing season of an evergreen Norway spruce dominated forest stand. We optimized the optical radiative transfer routine of the SCOPE model to retrieve vegetation properties such as leaf area index and leaf chlorophyll, water, and dry matter contents. The results indicated percentage differences less than 30% between the retrieved and measured vegetation properties. Additionally, we compared the retrievals from HLS data with those from hyperspectral airborne data for the same site, showing that HLS data preserve a considerable amount of information about the vegetation properties. Time series of vegetation properties, retrieved from HLS data, served as the SCOPE inputs for the time series of GPP simulations. The SCOPE model reproduced the temporal cycle of local flux tower measurements of GPP, as indicated by the high Nash–Sutcliffe efficiency value (>0.5). However, GPP simulations did not significantly change when we ran the SCOPE model with constant vegetation properties during the growing season. This might be attributed to the low variability in the vegetation properties of the evergreen forest stand within a vegetation season. We further observed that the temporal variation in maximum carboxylation capacity had a pronounced effect on GPP simulations. We focused on an evergreen forest stand. Further studies should investigate the potential of HLS data across different forest types, such as deciduous stand.


Introduction
Vegetation is an essential component of the terrestrial ecosystems that interacts with the atmosphere through the carbon and water cycles. The leaf stomata regulate the release of water by properties (i.e., vegetation phenology) are also needed to reflect the vegetation responses to climate variability and were found to improve the temporal GPP estimated by the SCOPE model [20,22].
Accurate description of the temporal variations in vegetation properties depends on the availability of rich time series of high quality remote sensing data. A single sensor may not provide the needed temporal coverage due to cloud cover and sensor temporal revisit limitations. More than 300 Earth observation satellites for optical imaging are orbiting the Earth, allowing us to combine observations from multiple sensors for improved temporal coverage. However, a simple combination of observations does not result in smooth spectral time series, e.g., due to the differences in radiometric characteristics and mismatch in spatial resolution between the sensors. It is not always suitable to assess the vegetation phenology and retrieve consistent time series of vegetation characteristics from the combined datasets. Recently, the National Aeronautics and Space Administration (NASA) produced the harmonized Landsat and Sentinel (HLS) data set by combining the surface reflectance data from the Operational Land Imager (OLI) and Multi-Spectral Instrument (MSI) sensors aboard Landsat-8 and Sentinel-2 satellites, respectively [30], to further improve the temporal resolution of the combined product. In HLS data, a smooth spectral time series is produced by accounting for the differences in spatial resolution, sensor-to-sensor differences in atmospheric correction approaches, and view geometry and radiometric characteristics of spectral bands. The potential of HLS time series data has not yet been fully explored with respect to GPP estimations. Recently, Lin et al. [31] applied a data-driven approach using HLS data to evaluate GPP variations in natural vegetation types. However, the full potential of HLS time series data has not yet been explored with respect to GPP estimation using data-assimilation approaches. The objective of this study was to fill this gap by assimilating multi-sensor spectral information obtained from the HLS data within the SCOPE model to retrieve the temporal variations in biochemical and structural vegetation properties to evaluate the improvement in GPP estimation. This research addresses the following two questions: (1) What information is preserved or lost in the retrieval of vegetation properties using multispectral HLS data? (2) What improvement in the GPP estimation can be expected by providing the SCOPE model with the temporal variations in the vegetation properties? We assumed that HLS multispectral data can provide sufficient information for retrieval of vegetation properties. However, whether such information can improve GPP estimation needs more investigation. We aimed to gain insights regarding these two aspects.

Study Area
We used measurements recorded at the Bílý Kříž ecosystem station in the Czech Republic, which is part of the Czech Carbon Observation System (CzeCOS; http://www.czecos.cz/) network. The Bílý Kříž site is situated in the northeast Czech Republic (49 • 30 07.474 N, 18 • 32 12.777 E) at an altitude of 875 m above sea level on a 13 • southwest-oriented planar slope that is approximately 100 m downslope of a mountain crest extending from west to east [32]. This site is a Class 2 candidate site of the Integrated Carbon Observation System (ICOS) (www.icos-cp.eu) and it regularly contributes to the FLUXNET dataset labeled as CZ-BK1 [33]. The site is predominated by a Norway spruce (Picea abies (L.) H. Karst) forest stand that represents the evergreen needle leaf forest vegetation class. The stand had a mean height of 17 m and the age of 35 years at the end of the 2016 vegetation growing season. This site is located in a moderately cold, humid, and precipitation-rich temperate climatic zone [34]. The soil at the Bílý Kříž site is an Entic Podsol type with a texture ranging from sandy loam to sandy clay with a 15-35% clay fraction [35]. The site is equipped with an eddy covariance flux tower system that has been measuring the energy and trace gas fluxes between the ecosystem and the atmosphere since 2004. We investigated the retrieval of vegetation properties and the estimation of GPP from remote sensing data for the 2016 vegetation growing season (April to September) when airborne hyperspectral data and ground measurements of vegetation biochemical and structural properties were available for this site.

Remote Sensing Data
We obtained a time series of 39 images at 30 m resolution from the HLS data product that covered the study area from 1 April to 30 September 2016. The HLS data are atmospherically corrected with spatial coregistration, bidirectional reflectance distribution function normalization, and spectral bandpass adjustments of Landsat-8 MSI and Sentinel-2 OLI sensors. The MSI surface reflectance bands are also adjusted to Landsat-8 characteristics by applying the spectral response function of the OLI sensor. The HLS surface reflectance data product is distributed globally within the Sentinel-2 tiling system and is abbreviated as S30 (obtained from Sentinel-2 MSI sensor) and L30 (obtained from Landsat-8 OLI sensor). For more information about the processing of MSI and OLI observations within the HLS data product, the reader is referred to Skakun et al. [36] and Claverie et al. [30]. We included the time series of HLS data available between 11:38 to 11:54 a.m. Central European Time during the 2016 vegetation growing season. Further filtering of the HLS data was performed based on the quality indicators using NDVI and near infrared (NIR) reflectance (see Section 3.2 for more details).
We also obtained the hyperspectral airborne data over the study area. We compared the retrieval of vegetation properties from both airborne and HLS data to examine if multispectral HLS bands preserve the required information compared to hyperspectral bands. The airborne data were acquired around local solar noon between 11 a.m. and 2 p.m. on 31 August 2016 with two push broom spectroradiometers (Compact Airborne Spectrographic Imager (CASI) and shortwave Infrared (SWIR) Airborne Spectrographic Imager (SASI) from Itres Ltd., Calgary, AB, Canada) onboard the Flying Laboratory of Imaging Systems operated by the Global Change Research Institute [37,38]. The CASI operates in the visible and NIR regions between 372 and 1044 nm (72 spectral bands with a sampling distance of 9.4 nm) and a nominal spatial resolution of 1.0 m. SASI operates in the shortwave infrared region between 957 and 2442 nm (100 spectral bands with a sampling distance of 15 nm) and nominal spatial resolution of 2.5 m. The airborne data were corrected radiometrically (using the factory calibration coefficients in RCX software (Itres Ltd., Calgary, AB, Canada)), geometrically (using GeoCorr software (Itres Ltd., Calgary, AB, Canada)), and atmospherically (using ATCOR-4 software [39]), and combined into a single hypercube according to the processing chain established at Global Change Research Institute [37]. The quality of the atmospheric corrections was evaluated by comparing the airborne spectra extracted for homogeneous artificial surfaces with spectral signatures measured in the field during the overflight, which produced maximum differences of 3% in NIR and SWIR wavelengths for some targets ( Figure S1). Both HLS and airborne data were used to retrieve vegetation properties at the Bílý Kříž site (Section 3.2). For this purpose, we extracted the surface reflectance of the vegetation canopy, i.e., top-of-canopy (TOC) reflectance, from the target area delineated by the ecosystem boundary (Section 3.1.2).

Ground Measurements
We obtained half-hourly eddy covariance (EC) measurements of net ecosystem exchange (NEE) to predict half-hourly GPP and standardized meteorological variables for the spruce forest stand during the 2016 growing season. The EC technique [40] was applied to measure turbulent fluxes of CO 2 at the Bílý Kříž forest site. The EC system consisted of Gill ultrasonic anemometer (HS-50, Gill Instruments, Hampshire, UK) and a LI-COR infrared gas analyzer (LI-7200, LI-COR, Lincoln, NE, USA), and was installed on a meteorological tower at 20.5 m as of 1 October 2013 and moved to 25 m as of 7 June 2016 above ground level. The EC raw data measured at 20 Hz were processed using spike detection and removal [41], time lag compensation, sonic temperature correction [42], and high- [43,44] and low [45]-frequency spectral corrections. Coordinates were rotated using the planar-fit method [46]. Fluxes were computed at half-hourly time intervals using the block averaging method. Fetch filtering was applied to exclude half hours for which more than 30% of the signal originated from outside of the target area. This was achieved by comparing fetch distance at a given upwind direction with estimated distance providing 70% contribution to the EC fluxes using footprint modeling [47,48]. Fetch distance for all directions is represented as an ecosystem boundary, i.e., the target area (Figure 1), that delimits the area that is sufficiently homogeneous considering the canopy height and the species represented (Norway spruce at the Bílý Kříž site).
All EC processing and calculations were performed using EddyPro software (v6.2.0, LI-COR, Lincoln, NE, USA). A thorough data quality checking procedure was applied to EC measurements in this study using R package openeddy (R Core Team, 2020; https://github.com/lsigut/openeddy; see McGloin et al. [32] for details). Measured CO 2 fluxes were filtered based on the friction velocity threshold computed using the moving point method [49]. After this filtering, half-hourly CO 2 fluxes were assumed to represent the NEE (µmol CO 2 m −2 s −1 ). Gaps were filled using marginal distribution sampling [50] and NEE was partitioned into GPP (hereafter, GPP EC ) and ecosystem respiration according to Lasslop et al. [6] using the R package REddyProc [51]. Half-hourly GPP EC were then aggregated to daily sums of carbon, expressed as g C m −2 s −1 .
Meteorological variables included half-hourly data of air temperature (T a , in°C), TOC incoming shortwave radiation (R in , in W m −2 ), TOC incoming longwave radiation (R li , in W m −2 ), vapor pressure deficit (e a , in hPa), air pressure (p, in hPa), and wind speed (u, in m s −1 ) for each day in the vegetation growing season of 2016. We additionally obtained half-hourly soil water content (SWC, in %) from the nine sensors installed at a 5 cm depth at different locations within the target area. For each half-hour, SWC values from all nine sensors were averaged to obtain the representative value of soil moisture of the target area. The half-hourly meteorological variables, mean SWC, along with HLS-, CASI-, and SASI-derived information served as the inputs to the SCOPE model to estimate (or simulate) half-hourly GPP, which was aggregated to daily sums of carbon (Section 3.3). We refer to simulated GPP as GPP SIM hereafter. We validated daily GPP SIM against daily GPP EC . We also compared daily GPP SIM with daily means of half-hourly measurements of meteorological variables to show the relationship between them. Table 1 lists the summary statistics of daily meteorological variables during the 2016 vegetation growing season.  Additionally, we obtained ground measurements of vegetation properties characterizing spruce forest conditions at the site, which were collected during a joint ground/airborne campaign on 31 August 2016. The location of the field sampling plots was pre-determined by previously existing forest inventory plots established by the Institute of Forest Ecosystem Research (IFER), Czech Republic, in the vicinity of the ecosystem station. The coordinates of the central point of each plot were measured using a differential global positioning system. The measured biochemical vegetation properties included leaf chlorophyll content (C ab ), leaf carotenoids content (C ca ), leaf water content (C w ), and leaf dry matter content (C dm ). The biochemical properties were measured in a laboratory (leaf pigments were extracted according to Porra et al. [52]) for needles sampled from the top and bottom part of the crown of three representative trees selected at each sample plot. The measured structural property was the leaf area index LAI, which was determined at the plot level using a Plant Canopy Analyser LAI-2200 instrument (LI-COR, Lincoln, NE, USA). More details about the ground measurements and data pre-processing are reported in Homolová et al. [38]. These ground measurements were compatible with airborne hyperspectral data acquired on 31 August 2016. For this study, we used the minimum, mean, and maximum measurements of three plots that were located inside and right at the edge of the target area ( Figure 1). These three plots represented a Norway-spruce-dominated forest stand. Table 2 provides the means of the measured vegetation properties. We validated the retrieved vegetation properties against the measured means.

Retrieval of Vegetation Properties Time Series Using HLS and Airborne Data and the SCOPE Model
The SCOPE model is an integration of four modules [19] that interact with each other to simulate the optical properties and physiological state of plants. The modules include: (a) the radiative transfer module in the optical domain (RTMo) for simulating TOC reflectance by tracking the propagation of incident solar and sky radiation; (b) the radiative transfer module for thermal radiation (RTMt), to simulate TOC outgoing thermal radiation ; (c) energy balance and biochemical routines to simulate heat fluxes and photosynthesis; and (d) a radiative transfer module that simulates the TOC spectrum for chlorophyll fluorescence (RTMf). In this study, we inverted the RTMo module as a separate model to retrieve vegetation properties and finally simulated the photosynthesis (i.e., GPP SIM ) using forward SCOPE modeling.
The leaf level model PROSPECT5 [53] and the canopy-level reflectance model 4SAIL (scattering by arbitrary inclined leaves) [54] are combined in the RTMo model of SCOPE. The PROSPECT5 model simulates the leaf-level reflectance and transmittance from the input biochemical leaf properties. Leaf-level optical properties are further upscaled on canopy scale by accounting for the canopy architecture using the 4SAIL model. We also incorporated the effect of soil background (soil brightness and SWC) on the simulated reflectance using the brightness-shape-moisture (BSM) submodel as previously suggested [20,21]. Table 2 provides the input parameters of the RTMo model and the BSM submodel. Table 2. Input parameters for the radiative transfer module (RTMo) of the Soil Canopy Observation of Photosynthesis and Energy fluxes (SCOPE) model together with the lower bound (LB), upper bound (UP), prior mean (µ), prior standard deviation (SD), and mean of the measurements acquired in August 2016 from three field sampling plots ( Figure 1). The range of the measured soil water content (SWC) in the brightness-shape-moisture (BSM) submodel is provided from Table 3.

Parameter
Symbol We evaluated the quality of HLS data before the retrieval of vegetation properties. For the 2016 vegetation growing season, we downloaded HLS data for 39 days. The HLS data are distributed with a per-pixel quality layer to mask the pixels of poor quality (e.g., due to the presence of clouds or other processing artifacts). However, the accuracy of the internal cloud mask is insufficient [30], and its use may lead to the omission of cloud detection. This was also the case for our study area for some days. Therefore, we did not use the internal cloud mask of the HLS data, but propose our method for filtering the time series based on the empirical threshold of NDVI and NIR reflectance. We only considered pixels to be of high quality when the NDVI varied between 0.5 and 0.99 (vegetation pixels) and the NIR reflectance between 0.1 and 0.4 (typical spruce canopy reflectance values as, for instance, reported by Rautiainen et al. [55]). Thresholds were set based on our previous experience with the reflectance characteristics of Norway spruce in the study area during the growing season. For each day, we calculated the mean of TOC reflectance in each HLS band within the target area if more than 50% of the pixels were of high quality (fulfilling both the NDVI and NIR thresholds). Of the original 39 days, 14 days were retained after the filtering. For the retrieval, we selected six HLS bands relevant for studying the vegetation properties that were common in both Sentinel-2 and Landsat-8 data. These selected bands were blue (482 nm), green (561.4 nm), red (654.6 nm), NIR narrow (864.7 nm), and two shortwave infrared bands (SWIR 1 (1608.9 nm) and SWIR 2 (2200.7 nm)). Table 3 shows the 14 days with the percentage of high quality pixels within the target area.
For each of the 14 days with HLS data, we used an optimization method to invert the RTMo model against mean TOC reflectance to retrieve the dynamics of the vegetation properties, as described in van der Tol et al. [24]. To tune the parameters, we used nonlinear least squares optimization, implemented in the MATLAB (the MathWorks Inc., Natick, MA, USA) built-in function "lsqnonlin", to minimize a cost function that calculated the sum of squared differences between the simulated and measured mean TOC reflectance in the selected six HLS optical bands. In this study, we chose to retrieve five biochemical parameters (C ab , C w , C dm , C s , and C ca ; Table 2) and three structural parameters (LAI, leaf inclination distribution function (LIDF) and leaf structural parameter (N); Table 2). To run the optimization method, prior information, i.e., range, mean, and standard deviation of each vegetation property to be retrieved ( Table 2) was needed. We used non-informative prior information following Verhoef et al. [21] and Bayat et al. [22] to ensure the significant impact of the measured HLS reflectance on the retrieved parameters. We used the middle of the total range (i.e., (lower bound + upper bound)/2) as a prior mean of each vegetation property. Assuming uniform prior distribution over the range between the lower and upper bound produced values of standard deviation equal to 1/ √ 12 ≈ 0.3 times the range of each vegetation property. We fixed the SWC based on the measurements of mean SWC within the target area measured at the observation time of HLS data (Section 3.1.2 and Table 3). We used the BSM submodel to describe soil background reflectance. The BSM submodel is based on the library of global soil vectors [56], which is extended via a brightness-shape component and SWC effect. This has two main advantages: (1) the model can simulate soil reflectance, providing an effective alternative to soil spectrum measurements; and (2) the model can scale the dry soil spectrum for any given SWC. Further, the BSM submodel can be inverted to retrieve BSM input parameters (soil brightness (B), Lat, and Lon) from the given dry soil spectrum. We took advantage of this possibility and retrieved these three input parameters by inverting the BSM submodel against a dry soil reflectance extracted from the airborne data.
We also calculated a mean TOC reflectance within the target area from the airborne hyperspectral data acquired on 31 August 2016. It was used to retrieve the vegetation properties using the RTMo inversion, where we used the same non-informative prior information of the vegetation properties and the BSM input parameters as we used for HLS data. We further used the ground measurements (Table 2) to compare the retrieved vegetation properties obtained from HLS and airborne data. Table 3. Overview of the harmonized Landsat-8 and Sentinel-2 (HLS) data used in this study, with the percentage of high-quality pixels and mean NDVI within the target area defined by flux footprint (delineated in Figure 1). CET: Central European time, NDVI: Normalized difference vegetation index, L30: Landsat-8 data, S30: Sentinel-2 data, SWC: soil water content.

Simulating Gross Primary Production Using SCOPE Model
After retrieving the time series of the forest's biochemical and structural properties, we set up the SCOPE model to simulate the half-hourly GPP SIM during the 2016 vegetation growing season. These half-hourly values were summed to daily GPP SIM values calculated in g C m −2 s −1 . The SCOPE model inputs for GPP SIM included: (1) half-hourly meteorological variables, i.e., R in , R li , T a , p, e a , and u (introduced in Section 3.1.2); (2) vegetation properties, i.e., C ab , C w , C dm , C s , C ca , LAI, LIDF, and N (Table 2); (3) maximum rate of carboxylation (V cmax , in µmol m −2 s −1 ); and (4) the Bell-Berry stomatal parameter m. We ran the SCOPE model with three scenarios to investigate the impact of information retrieved from the HLS data on the accuracy of GPP SIM : "Fixed" scenario: We assumed no variations in the retrieved vegetation properties during the growing season. We fixed them either at the available measured means ( Table 2) or at the prior means. LIDF a and LIDF b were fixed at -0.35 and -0.15, respectively, showing the spherical characteristics of Norway spruce [54]. We further fixed V cmax = 80 µmol m −2 s −1 and m = 9, which are the typical values for evergreen needle leaf species as suggested in the SCOPE model. This scenario allowed us to monitor the temporal variation in GPP SIM determined by changes in radiation and atmospheric demand.
"HLS info" scenario: We used the time series of vegetation properties retrieved from the HLS data as the inputs to the SCOPE model. In this scenario, we also fixed V cmax and m at the typical values for evergreen needle leaf species, since these two cannot be directly retrieved from the optical reflectance [20].
"HLS info & var Vcmax" scenario: We used the same time series of retrieved vegetation properties as in the HLS info scenario. Here, instead of using a fixed typical value of V cmax , we estimated its time series from the retrieved time series of C ab from the HLS data. For this, we used the empirical relationship between V cmax [µmol m −2 s −1 ] and C ab [µg cm −2 ] proposed by Houborg et al. [57] for C 3 plants: where a = 2.529 and b = -27.34. This scenario allowed us to monitor the temporal variation in GPP SIM when we relied on optical domain information both directly from the retrieved vegetation properties and indirectly from the empirical relationships.

Statistical Evaluation of Model Performance
We determined the performance of the SCOPE simulations using statistical criteria that evaluated the efficiency with which the SCOPE reproduced the flux tower GPP EC . These criteria provided a measure of the SCOPE model efficiency in simulating daily GPP SIM over the 2016 vegetation growing season. We used the following criteria: where n is the number of daily simulated GPP SIM (y i ) and daily measured GPP EC (z i ), and z and y represent the mean over the vegetation growing period. The RMSE criterion has the unit of GPP. A low RMSE indicates high accuracy. It is appropriate to compare the simulations under different scenarios (Section 3.3). The NSE [58] can range from −∞ to 1. An NSE value close to 1 indicates a perfect match of simulations to the observations. Following Dumont et al. [59], we assumed that an NSE ≥0.5 indicates adequate accuracy in GPP SIM . The COR represents the linear relationship between simulations and observations, and can vary from -1 to 1. We also used Equations (3) and (4) to evaluate the efficiency with which the RTMo produced TOC reflectance on each day. In this case, z i and y i represent the observed and the RTMo-simulated HLS/airborne TOC reflectance at the wavelength i, respectively. We compared the vegetation properties retrieved from the HLS data with those retrieved from the airborne data, and the measurements using the percentage difference: where A i and B i represent the ith vegetation property retrieved from HLS data and airborne data, respectively. For comparison with the measurements, B i represents the ith measured vegetation property. Figure 2a shows the temporal variation in the mean TOC reflectance extracted from the HLS data within the target area. The highest variation was observed in the NIR band with increasing NIR reflectance from the start of the vegetation growing seasons, which then stabilized mid-season, and decreased at the end of vegetation growing season. Other bands exhibited the same pattern, but with smaller variations in the magnitude of reflectance compared to the NIR band. The observed temporal variation in the TOC reflectance likely did not occur due to the variation in LAI because the evergreen Norway spruce species can be expected to have stable LAI (Section 4.3) throughout the year. Instead, change in leaf age distribution caused temporal variation in the reflectance. Both young and mature leaves contribute to the TOC reflectance when captured from the top of the crown by the sensors of HLS data. Young leaves, in particular, start becoming greener from the beginning of the vegetation growing season, and then they stay sufficiently green until September and then start browning. This phenological development of the young leaves was likely the cause of variation in the TOC reflectance observed in this study. Figure 2b compares the mean and standard deviation, showing the variability of TOC reflectance within the target area, extracted from the HLS and airborne data acquired on 31 August 2016. To facilitate comparison, we also plotted the airborne reflectance resampled at the band centers of the HLS data. Both the HLS and airborne data exhibited similar magnitudes of mean TOC reflectance, with slightly higher NIR reflectance (0.19) in the HLS compared to airborne data (0.16). We found that both the HLS and airborne data exhibited similar variability within the target area even though the stand variability was poorly captured in the HLS data with coarse spatial resolution compared to the airborne data.

Simulation of TOC Reflectance
We assessed the performance of the RTMo model in retrieving vegetation properties by: (1) using the goodness-of-fit measured between the simulated and observed TOC reflectance; and (2) comparing the retrieved with the measured vegetation properties (Section 4.3). Figure 3 illustrates the goodness of fit using COR and NSE, and provides a spectral comparison between the observed mean TOC reflectance (by HLS and airborne) and the simulated TOC reflectance (by the RTMo model) for all 14 days during the vegetation growing season. We found that the NSE was close to one on each day, indicating a perfect match between the simulated and observed TOC reflectance. Strong positive correlations between simulations and observations were indicated by a COR value close to one. These showed that the RTMo model was able to successfully reproduce the observed HLS and airborne TOC reflectance. We presented the observed and simulated TOC reflectance curves for the selected Julian days in Figure S2.  Figure 1). The correlation coefficient (COR) and Nash-Sutchliffe Efficeince (NSE) between observed and simulated TOC reflectance are shown in each plot.  (Table 2). An optimization process constrained the RTMo input by the TOC reflectance (either HLS or airborne). Each vegetation property was retrieved in a narrow range during the growing season, indicating the reduction in prior uncertainty after the optimization.

Retrieval of Vegetation Properties
From HLS data, we were able to retrieve the vegetation properties at least once per month. Only in June was it not possible to find a reliable observation for use in the retrieval. These retrievals provide insight into the temporal patterns of vegetation properties. Except for C ab and C dm , we did not observe any apparent temporal variation in any of the studied vegetation properties. We retrieved a low estimate of C ab at the start of the growing season (23.4 µg cm −2 on 3 April), which increased to 50 µg cm −2 in mid-season (5 July), and then showed an overall decreasing trend toward the end. C dm followed a similar pattern as C ab during the growing season, i.e., both decreased or increased together between any two successive days. The retrieval of LAI showed consistency during the growing season, and we retrieved the LAI between 7 and 8 m 2 m −2 , except for 23 May, where it dropped to 6 m 2 m −2 . This drop is, however, not an unrealistic estimate of LAI for Norway spruce.
We obtained the percentage difference between the vegetation properties retrieved from the HLS and airborne data using Equation (5) to show the similarity between them. We observed low percentage differences in most of the vegetation properties. The percent differences for C ab , C dm , C ca , and LAI were less than 10%, with the lowest being 0.3% for C ab . C w and LIDF a showed a difference of 25%. The maximum differences in C s and N were observed at 63% and 42%, respectively. (a-i) Retrieved vegetation properties (C ab , C w , C dm , C s , C ca , LAI, LIDF, and N); (j) maximum rate of carboxylation (V cmax ) from the available time series of the harmonized Landsat-8 and Sentinel-2 data during the 2016 growing season. Retrievals from the airborne hyperspectral data for one specific acquisition on 31 August are indicated by the plus sign. Wherever field measurements were available, mean values (diamond sign) and minimum-maximum (Min-Max) ranges are indicated. Information about the vegetation properties is given in Table 2.
Based on the percentage difference, we compared the vegetation properties retrieved from the HLS data on 31 August with the mean of the available ground measurements (Table 2 and Figure 4). We found that the percentage differences for C ab , C w , C dm , C ca , and LAI were less than 30%. We found the lowest difference for LAI, with 7%, whereas C ab showed a difference of 17%. These percentage differences appear to be acceptable given the aggregated effect of the target area compared with the three separate plots of measurements ( Figure 1). Moreover, all the retrieved vegetation properties were within the minimum and maximum, i.e., the range of the measurements (Figure 4), indicating the retrievals were reasonable, at least on 31 August 2016.

Simulations of GPP SIM with SCOPE Model
Finally, we evaluated the performance of the SCOPE model by comparing the simulated daily GPP SIM with the daily GPP EC measured during the 2016 vegetation growing season. Figure 5 illustrates the temporal variation in daily GPP SIM for three designed scenarios (Section 3.3), together with the daily GPP EC . We also present the scatterplots between daily GPP SIM and GPP EC to show their correlation between them. Under both the Fixed and HLS info scenarios, daily GPP SIM closely followed the daily GPP EC (Figure 5a,c), visually indicating the accuracy of GPP SIM . The quantitative description of the GPP SIM accuracy was confirmed with the estimated NSE value (Equation (3)), which was found to be greater than 0.5 for both scenarios. The high NSE value indicated the sufficient quality of the predictive power of the SCOPE model for simulating gross primary production. We, however, observed some under-and over-estimation in GPP SIM compared to the GPP EC . The scatterplots (Figure 5b,d) further explain the positive linear relationship between GPP SIM and GPP EC , with the value of COR being high (close to 0.9 for both scenarios). This high value indicated that the variation of GPP SIM between consecutive days were synchronized with that of GPP EC . We also concluded that both Fixed and HLS info scenarios led to similar magnitude and variation in daily GPP SIM during the vegetation growing season, and both scenarios quantitatively produced similar NSE values, together with similar RMSE values (Equation (2)).
Under the HLS info & var Vcmax scenario, we estimated V cmax from the time series of C ab (Equation (1)) as ranging from 31.8 to 108.5 µmol m −2 s −1 during the vegetation growing season (Figure 4j). The variation in V cmax resulted in the synchronization between daily GPP SIM and GPP EC , but substantial differences in their magnitudes (Figure 5e) compared to the other two scenarios. These indicated the decrease in the quality of the predictive power of the SCOPE model for GPP SIM . The drop in NSE below 0.5 and the increase in RMSE provided quantitative evidence of this decrease. The linearity in the relationship between GPP SIM and GPP EC also reduced, which resulted in a drop in COR (Figure 5f).
Most of the variability in daily GPP SIM can be explained by its response to meteorological variables. We showed that R in was the main driver of GPP SIM , with strong dependence especially up to daily mean R in = 200 W m −2 (Figure 6a). The GPP SIM response to R in was further modified by T a (Figure 6b). The GPP SIM response to T a can be roughly separated to two regimes by a threshold of ≈13°C. The colder regime (<13°C) was associated with generally lower GPP SIM than during the warmer regime (>13°C). Most of the colder regime days occurred during the beginning of the growing season, whereas the warmer regime occurred during the core of the growing season. The GPP SIM values for the end of the growing season fell during both regimes but typically represented the lower end of the GPP SIM range within a given regime (Figure 6b). Although T a was the lowest mostly in the beginning of the growing season and could be expected to form the lower bound of the light response curve, GPP SIM values for the end of the growing season showed even lower values when accounting for R in .

Discussion
This study provided insight into the potential of using newly produced HLS data for the retrieval of temporal variations in vegetation properties using the RTMo model, and we investigated if these temporal variations could improve GPP SIM . In addition, we investigated the information contained about vegetation properties in the multispectral HLS data compared to hyperspectral airborne data.
The amount of HLS data contributing to the time-series was reduced by cloud coverage condition over the site.The cloud mask is not properly defined in the HLS data, and this was also observed in this study. The non-cloud pixels were defined based on the NDVI range (0.4 to 0.5) and NIR range (0.1 to 0.4). The NDVI range adopted in this study was based on the broad range of NDVI for vegetation, and the NIR range was based on our previous experience on the characteristics of Norway spruce during the growing period in the study area. The variations in these ranges can be observed for other species and the study period [60]. For example, the NIR reflectance can drop to below 0.1 outside the vegetation growing period and can rise to above 0.4 for grey alder tree species. Therefore, NDVI-and NIR-based filtering of cloud-free HLS data may be suitable at the local level and for shorter periods, like in this study. However, this filtering solution might be unsuitable if the study involves several tree species in a broad area and a longer study period. Therefore, the existing cloud mask of the HLS data, i.e., irrespective of tree species/land cover types, should be made more robust.
The RTMo model successfully reproduced the observed HLS and airborne TOC reflectance ( Figure 3).
We retrieved the time series of vegetation properties from the HLS data. We had a limited number of ground measurements for validation of the retrievals. However, the following points should be noted: 1. A temporal increase in C dm is positively correlated with C ab . A previous study, however, did not focus on the temporal aspects of evergreen as well as deciduous tree species, but also found a positive correlation between C dm and C ab [68]. 2. We observed a good fit between the retrieved vegetation properties and the ground measurements (Section 4.3) on 31 August. At least for one day in the time series, we could therefore validate the retrievals from the HLS TOC reflectance. The retrieved C ab , however, did not exhibit substantial temporal variation. Nevertheless, these variations followed the characteristics of the evergreen tree species, which shows young leaf development phenology during the growing season. These findings agree with those of a previous study on the global spatio-temporal distribution of leaf chlorophyll [69], which highlighted the consistency in the temporal leaf chlorophyll profile of evergreen tree species across the year with an increasing concentration within new needles in spring.
Considering the points described above, we inferred a successful retrieval of vegetation properties from the HLS data. In a future study, these two points could be considered together for validation of retrievals, at least for evergreen tree species, if limited measurements are available.
We compared the retrievals from HLS and airborne data. The low percentage differences in most of the vegetation properties indicated that switching from hyperspectral to multispectral observations did not prominently influence the retrieval of vegetation properties in our case. In other words, the multispectral HLS data preserve a considerable amount of information about vegetation properties. This finding is in agreement with Croft et al. [70] who reported similarity between the C ab values over the range of tree species belonging to different biome types, retrieved from hyperspectral measurements and Landsat observations. Further, our finding is in line with Darvishzadeh et al. [71], who exploited the potential of Sentinel-2 data for the mapping of C ab of Norway spruce. Our study highlighted the potential of the newly produced HLS data set to monitor the vegetation properties by taking advantage of its denser time series.
We further examined if the time series of vegetation properties, retrieved from the HLS data, could improve the simulation of GPP SIM by the SCOPE model. We observed that both constant (Fixed scenario) and time-varying vegetation properties (HLS info scenario) led to almost identical GPP SIM during the vegetation growing season (Section 4.4). These results confirmed that the GPP SIM in Norway spruce is mainly driven by meteorological conditions (Figure 6). Therefore, the time-varying vegetation properties' information could be significant for GPP SIM only when the tree is either under stress condition (e.g., due to drought) or exhibits strong temporal phenology (such as deciduous tree species). In a previous study, Bayat et al. [20] reported an improvement in temporal GPP SIM of the plant under stress using the time-varying vegetation properties input to the SCOPE model. In our study, the SCOPE model was mainly driven by the meteorological conditions because the Norway spruce did not show pronounced temporal variations in vegetation properties (Section 4.3) and were under no significant stress condition during the selected episode [32]. Therefore, constant vegetation properties were sufficient for simulating temporal GPP SIM with reasonable accuracy. We further observed the influence of time-varying V cmax (HLS info & var Vcmax scenario) on GPP SIM . We expected that some fluctuations in GPP SIM (Figure 5a,c) would be further improved with time-varying V cmax instead of using a constant value throughout the vegetation growing season. The time-varying V cmax , however, reduced the accuracy of GPP SIM (Figure 5f), which could be attributed to the empirical relationship of V cmax and C ab defined for C 3 plants (Equation (1)). A more robust empirical relationship between the C ab and V cmax must thus be defined for this site. A biome specific definition [72] could improve the estimation of V cmax , but this should be tested at the local site in a future study.

Conclusions
The present analysis showed the potential of the newly produced multispectral harmonized Landsat-8 and Sentinel-2 (HLS) data to retrieve time series of vegetation properties (i.e., biochemical and structural), which were further used as inputs to the Soil Canopy Observation, Photochemistry and Energy fluxes (SCOPE) model to simulate the time series of gross primary production GPP SIM . The study led to the following conclusions: 1. HLS data can provide the dense time series of surface reflectance at the desired locations because it combines data from two existing satellites: Landsat-8 and Sentinel-2. The results demonstrated that the HLS data are capable of preserving the needed information about vegetation properties. We investigated the retrieval only for one evergreen tree species. Nevertheless, our analysis has an important implication for the future use of HLS data for the dense time series retrieval of vegetation properties, which is needed for other tree species, such as deciduous species, showing strong temporal phenology. 2. We did not observe any improvement in GPP SIM using the time-varying vegetation properties retrieved from the HLS data. We observed the influence of the time-varying maximum rate of carboxylation (V cmax ) on GPP SIM . However, the empirical relationship used to estimate V cmax from leaf chlorophyll content (C ab ) decreased the accuracy of GPP SIM . Future studies need to redefine this empirical relationship.

Conflicts of Interest:
The authors declare no conflict of interest.