Improving the Capability of the SCOPE Model for Simulating Solar-Induced Fluorescence and Gross Primary Production Using Data from OCO-2 and Flux Towers

Solar-induced chlorophyll fluorescence (SIF) measured from space has shed light on the diagnosis of gross primary production (GPP) and has emerged as a promising way to quantify plant photosynthesis. The SCOPE model can explicitly simulate SIF and GPP, while the uncertainty in key model parameters can lead to significant uncertainty in simulations. Previous work has constrained uncertain parameters in the SCOPE model using coarse-resolution SIF observations from satellites, while few studies have used finer resolution SIF measured from the Orbiting Carbon Observatory-2 (OCO-2) to improve the model. Here, we identified the sensitive parameters to SIF and GPP estimation, and improved the performance of SCOPE in simulating SIF and GPP for temperate forests by constraining the physiological parameters relating to SIF and GPP by combining satellite-based SIF measurements (e.g., OCO-2) with flux tower GPP data. Our study showed that SIF had weak capability in constraining maximum carboxylation capacity (Vcmax), while GPP could constrain this parameter well. The OCO-2 SIF data constrained fluorescence quantum efficiency (fqe) well and improved the performance of SCOPE in SIF simulation. However, the use of the OCO-2 SIF alone cannot significantly improve the GPP simulation. The use of both satellite SIF and flux tower GPP data as constraints improved the performance of the model for simulating SIF and GPP simultaneously. This analysis is useful for improving the capability of the SCOPE model, understanding the relationships between GPP and SIF, and improving the estimation of both SIIF and GPP by incorporating satellite SIF products and flux tower data.


Introduction
Solar-induced chlorophyll fluorescence (SIF), emitted by vegetation in the pigment beds of photosystems, is an indicator of the efficiency by which photons are transmitted into photochemical reaction centers [1]. SIF has recently emerged as a novel tool for the estimation of terrestrial gross primary production (GPP), the amount of carbon fixed by terrestrial vegetation via photosynthesis [2,3]. SIF is seen as a more promising proxy for GPP than the traditional vegetation indices that are derived from reflectance data [4][5][6]. GPP is a key component of the global carbon cycle, and its accurate estimation at region and global scales is vital for understanding the interactions between the terrestrial carbon cycle and climate change [7,8]. Previous studies have shown strong relationships between ecosystem SIF and GPP using SIF data from ground measurements (e.g., [9,10]), airborne imaging spectrometers (e.g., [4,11]), and satellite sensors such as the Greenhouse Gases the most widely used sequential method is perhaps the Kalman filter [28]. Since the conventional optimization methods are computationally intensive, particularly for complex models with a large number of parameters, efforts have been made to reduce the computational burden of parameter optimization. The surrogate-based optimization, replacing the original expensive models with simplified statistical models during the optimization process, is one of the commonly used ways to efficiently optimize complex dynamic models [36][37][38]. Some previous studies have used surrogate models or emulators of SCOPE (e.g., [39,40]). The performance of the emulators was typically determined by the machine learning algorithm, dimensionality reduction method, and the number of components during the modeling [37,40]. Given the efficiency and effectiveness of the surrogate-based optimization, here we used the surrogate-based method in our study.
Meanwhile, previous studies have demonstrated the importance of optimizing the key parameters in simulating photosynthesis process using the SCOPE model (e.g., [19,41]). The maximum carboxylation capacity (V cmax , µmol CO 2 m −2 s −1 ) of RuBisCO is one of the most important parameters in the photosynthesis submodel [42]. SIF is related to the electron transport rate in the process of photosynthesis [3], and there are complex and intrinsic linkages among the absorbed photosynthetically active radiation (APAR), fluorescence radiance, and the rate of photochemistry [16,43]. Declining V cmax reduces the saturated rate of photosynthesis and switches the light-limiting condition to the light-saturating condition at lower PAR, leading to the maximum photochemistry rate at which photochemical yield declines [43]. With the increase in PAR and nonphotochemical quenching, fluorescence yield first increases and then declines, while the fluorescence radiance (expected SIF) rises first nonlinearly and then linearly [43]. At a midday satellite overpass with high PAR, SIF tends to be higher for leaves with high V cmax than for leaves with low V cmax , leading to a positive relationship between SIF and V cmax [43]. Due to its intrinsic relationship with vegetation photosynthetic activity, SIF has also been used to estimate V cmax to better simulate GPP and SIF using the SCOPE model [41]. Several studies have evaluated the sensitivity of SIF and GPP to this parameter [39,41,44], and have shown that the SIF-V cmax relationship was not stable and instead varied between different versions of the SCOPE model [19,39,45,46]. Two recent studies showed that SIF was less sensitive to V cmax in SCOPE v1.6 [19,44]. However, the sensitivity of SIF to V cmax in the SCOPE v1.7 model across different biomes has not been examined yet. In addition, no studies have used satellite SIF observations and flux tower GPP data to jointly constrain the parameters of the SCOPE model for improving SIF and GPP simulations.
Several previous studies have used satellite SIF products (e.g., GOSAT and GOME-2) to improve the GPP estimation in cropland and grassland (e.g., [20,41]), which has demonstrated the potential of SIF observations in constraining GPP [19,41]. However, the satellite-derived SIF products from GOSAT or GOME-2 have coarse spatial resolution (e.g., 10 km diameter for GOSAT, 40 × 80 km 2 for GOME-2), leading to large scale mismatch between satellite grid cells and eddy covariance (EC) flux tower footprints [2]. The OCO-2 satellite, launched in 2014, has enabled the retrieval of finer-resolution (1.3 × 2.25 km 2 ) SIF than the previous satellites/sensors [47]. The ground area of the OCO-2 SIF soundings is close to the typical footprint of EC flux towers, and therefore the OCO-2 SIF data can be better matched with EC observations to examine the relationship between SIF and GPP or to simultaneously constrain models [5]. To date, only a few studies have used the OCO-2 SIF products to evaluate or improve the SIF simulation of the SCOPE model (e.g., [44]). Here, we used SIF observations in the far red region from OCO-2 and GPP data from flux towers to improve the performance of the SCOPE model for simulating SIF and GPP by optimizing the key uncertain parameters in the model. We chose two temperate forests (Park Falls and Willow Creek in the USA) as our study sites. The objectives of this study were to (1) identify the key parameters of the SCOPE model in predicting SIF and GPP; and (2) examine the effects of parameter optimization on the simulation of SIF and GPP, and assess how well the SCOPE model simulations can be constrained by OCO-2 SIF and flux tower GPP, both separately and together. Our results are useful for improving the capability of the SCOPE model, understanding the relationships between GPP and SIF, and revealing the potential of improving the SIF and GPP estimation by incorporating the SIF products and flux tower data.

Site Description and Flux Tower Observations
Our study sites consisted of two EC flux sites in northern Wisconsin, USA: Park Falls (US-PFa, 45.945 • N, 90.273 • W) [48,49] and Willow Creek (US-WCr, 45.806 • N, 90.080 • W) [50]. Both sites are located in the Chequamegon-Nicolet National Forest, and are characterized by an interior continental climate [49]. Park Falls is dominated by a mixed landscape of upland forests (mainly deciduous but with a significant coniferous coverage) and open wetlands [49]. The Willow Creek site is located about twenty km to the southeast of Park Falls. Willow Creek is a 80-100-year-old deciduous broadleaf forest (DBF) stand consisting mainly of sugar maple, basswood, and green ash, with a closed canopy of about 24 m in height and a leaf area index (LAI) of 5.3 [50].
We used EC and meteorological measurements from these two forest sites for the period 2014-2017 to constrain the SCOPE model, and investigated the relationships between satellite-derived SIF and flux tower GPP for the same period. EC flux and meteorological measurements have been made continuously at Park Falls [48] and Willow Creek since 1996 and 1998, respectively. We used the hourly GPP and meteorological data at Park Falls and half-hourly data from Willow Creek for the period from 2014 to 2017 for constraining the SCOPE model, and also tested the model with the data since 2003 at Park Falls. The net ecosystem exchange (NEE) measurements obtained from the EC technique were partitioned to GPP and ecosystem respiration using a non-linear regression of nighttime NEE to surface soil temperature [51]. The details of the flux data processing procedures of these two sites have been described in previous studies [51,52]. We also used meteorological data such as wind speed, air temperature, incoming solar radiation, specific humidity, and soil moisture as the driving data of the SCOPE model. The canopy height was 19 m and 24 m at Park Falls and Willow Creek, respectively, while wind speed was measured at the height of 30 m and 29.6 m, respectively. The roughness length and zero plane displacement were calculated from vegetation height and LAI in the SCOPE model as described in [53].

Solar-Induced Chlorophyll Fluorescence (SIF) Observations from Orbiting Carbon Observatory-2 (OCO-2)
The OCO-2 has been collecting SIF measurements with a repeat frequency of approximately 16 days since September 2014. The OCO-2 SIF was retrieved using a simplified forward transition model exploiting the infilling of the Fraunhofer lines at 757 nm and 771 nm wavelengths from the O 2 A-band observation records [47]. Although the spatial resolution of OCO-2 is 1.3 km × 2.25 km in the nadir mode, the total swath width is only 10.3 km. With the sparse coverage of the instrument, temporally dense SIF data are not available for most locations over the globe. The Park Falls site is a very tall tower and a validation site for the OCO-2 mission [2], and OCO-2 has been collecting SIF data in the target mode, leading to temporally well-spaced retrievals of SIF measurements near the flux tower site (Figure 1).
We obtained the level 2 SIF Lite product (OCO2_L2_Lite_SIF, version 8r) from the Land Processes Distributed Active Archive Center (LP DAAC). This product contains biascorrected SIF along with other selected fields from the IMAP-DOAS algorithm [47]. For each site, we extracted all the OCO-2 SIF observations that were available within the 2.5 km × 2.5 km area surrounding the flux tower during the study period. The instantaneous SIF observations were calculated by averaging all the available samples within the surrounding area. Then, we converted the instantaneous SIF measurements to daily SIF using the daily correction factor contained in the product [47]. We obtained the level 2 SIF Lite product (OCO2_L2_Lite_SIF, version 8r) Land Processes Distributed Active Archive Center (LP DAAC). This product bias-corrected SIF along with other selected fields from the IMAP-DOAS algorit For each site, we extracted all the OCO-2 SIF observations that were available w 2.5 km × 2.5 km area surrounding the flux tower during the study period. The ins ous SIF observations were calculated by averaging all the available samples w surrounding area. Then, we converted the instantaneous SIF measurements to d using the daily correction factor contained in the product [47].

Moderate Resolution Imaging Spectroradiometer (MODIS) Data Products
We obtained the MODIS reflectance product (MCD43A4, V6) and leaf ar (LAI) product (MCD15A2H, V6) from the Land Processes Distributed Active Center (LP DAAC) for the study period. MCD43A4 offers 500 m-resolution spe flectance adjusted by the nadir bidirectional reflectance distribution function based on reflectance data from both Terra and Aqua at the daily timescale, and th correction can minimize the effects of varying illumination and observation cond The base map is the MODIS land cover product (MCD12Q1) at 500 m spatial resolution.

Moderate Resolution Imaging Spectroradiometer (MODIS) Data Products
We obtained the MODIS reflectance product (MCD43A4, V6) and leaf area index (LAI) product (MCD15A2H, V6) from the Land Processes Distributed Active Archive Center (LP DAAC) for the study period. MCD43A4 offers 500 m-resolution spectral reflectance adjusted by the nadir bidirectional reflectance distribution function (BRDF) based on reflectance data from both Terra and Aqua at the daily timescale, and the BRDF correction can minimize the effects of varying illumination and observation conditions on surface reflectance [54]. The MCD15A2H LAI product was retrieved from both Terra and Aqua with 8-day intervals. For each MODIS product, we extracted the pixel values for the area (2.5 km × 2.5 km) surrounding each flux tower. For each 8-day interval, we averaged the values with good quality based on the quality assurance (QA) flags. Since the vegetation indices are widely used to monitor vegetation greenness and have significant relationships with vegetation GPP [7,8,55], we also derived the widely-used vegetation index, enhanced vegetation index (EVI), and examined its relationship with GPP to illustrate the performance of SIF.

Model Description and Parameterization
To improve the capability of SCOPE in simulating SIF emission and photosynthesis at the ecosystem level, we assessed how well OCO-2 SIF and flux tower GPP could constrain key uncertain model parameters. The SCOPE model is a 1-D vertical schematization of integrated radiative transfer model (RTM) and energy balance model that enables the simulation of canopy reflectance and fluorescence as well as the heat, carbon, and water fluxes between soil, vegetation, and atmosphere [15,16]. The radiative transfer and chlorophyll fluorescence routines of the SCOPE model are calculated with the module Fluspect [15], which is based on the vastly used canopy RTM, PROSAIL [56], with the addition of backward and forward fluorescence spectra. The PROSAIL model is a RTM that couples a leaf-level model, PROSPECT [57], with a canopy RTM, SAIL [58]. PROSAIL is perhaps the most commonly used RTM for reflectance and fluorescence applications. The updated version of SCOPE (V1.7) uses a new version of PROSPECT, PROSPECT-D [59], which adds anthocyanins to carotenoids and chlorophylls, the two plant pigments in the current model. The biochemical model in SCOPE is based on [42] and [60] for C3 and C4 plants, respectively. The SCOPE model can simulate fluorescence emission and photosynthesis simultaneously, and more details are provided in previous studies [15,16].
SCOPE contains a series of vegetation physiological parameters and requires meteorological data such as incoming solar radiation, wind speed, air temperature, and relative humidity. The LAI and leaf chlorophyll content (Cab) are two important structural and physiological parameters. Given the consistence of SCOPE and PROSAIL in radiation transfer routines, we retrieved the seasonal variations of Cab and LAI using the MODIS reflectance products of band 1 to band 7 and the PROSAIL model. In this study, the genetic algorithm (GA) approach was used to retrieve the structural and physiological parameters. The GA algorithm is a global optimization approach based on an analogy that consists of natural selection and evolutionary genetics [26,61]. The GA algorithm employs various operations of evolutionary strategies: encoding, fitness evaluation, parent selection, genetic operation, and replacement [61]. Assuming that the best chromosome is a raw fitness, f best , the fitness of the i-th chromosome (f i ) in the ordered list is conducted using a linear function during the evaluation: where d is the decrement rate. The average objective value of the population is mapped into the average fitness by Equation (1). The roulette wheel selection is commonly used to implement the proportionate scheme in GA. Giving the chromosome x with a fitness value f (x,t), its growth rate p i can be defined as follows: where F(t) is the average fitness of the population. We minimized the following objective function (Equation (3)): where Ref i,sim and Ref i,obs are the ith simulated and observed reflectance band, respectively, and n is the number of reflectance bands. The prior ranges of the parameters (e.g., Cab, LAI) of the PROSAIL model are given in Table 1. The LAI and Cab inverted by the PROSAIL model at Park Falls showed the seasonal patterns that were expected for temperate deciduous forests ( Figure 2). The inverted LAI was generally consistent with the MODIS LAI, and the inverted Cab showed a similar seasonal cycle as the LAI.

Parameter Sensitivity Analysis and Inverse Estimation of Key Paramete
The physiological parameters related to SIF and GPP simulation Table 1. The modules rely on different empirical parameterizations, w limited set of leaf-level data under some specific conditions. To exam the SCOPE model to predict SIF at satellite overpass time and asses parameters such as the Vcmax of RuBisCO can be constrained by OCO-2

Parameter Sensitivity Analysis and Inverse Estimation of Key Parameters
The physiological parameters related to SIF and GPP simulation are summarized in Table 1. The modules rely on different empirical parameterizations, which are based on a limited set of leaf-level data under some specific conditions. To examine the capability of the SCOPE model to predict SIF at satellite overpass time and assess how well the key parameters such as the V cmax of RuBisCO can be constrained by OCO-2 SIF and flux tower GPP, we conducted a parameter sensitivity analysis and estimated the key parameters. Three different objective functions were used to assess the sensitivity of physiological parameters: optimization using SIF as the constraint (F SIF , Equation (4)), optimization using GPP as the constraint (F GPP , Equation (5)), and optimization using both SIF and GPP as constraints (F SIF_GPP , Equation (6)). The performance of the model was evaluated using the coefficient of determination (R 2 ) of the linear regression between the measured and estimated values of SIF and GPP, measuring how much variation in the observations was explained by the models. The root mean square error (RMSE) was also used to evaluate the performance of the model. We minimized the following objective functions, respectively.
where SIF i,sim and SIF i,obs are the simulated and observed SIF at each time step, respectively; and GPP i,sim and GPP i,obs are the simulated and observed GPP at each time step, respectively; m is the number of timesteps; and std GPP and std SIF are the standard deviations of the observed GPP and SIF, respectively. The cost function for the optimization based on both SIF and GPP as constraints (Equation (6)) used the summed normalized errors of SIF and GPP to represent the total model error in simulating both SIF and GPP.
To identify the key parameters that are responsible for most of the variability in GPP and SIF simulated by SCOPE, both GSA method [31] and local SA approach were performed. In this study, we applied the Saltelli's method [31] (i.e., extended Fourier amplitude sensitivity test, EFAST), a quantitative variance-based SA method that can identify the effects of both parameters and the interactions of parameters on model simulations [31]. The EFAST combined Sobol's method and the Fourier amplitude sensitivity test (FAST) together [32], which can effectively quantify the main sensitivity effects (i.e., the first-order sensitivity index, SI 1st, representing the contribution of each input variable to the variance of model output), and total sensitivity effects of input variables (i.e., the total sensitivity index, SI total, the contribution to the total variance by the interactions between parameters) [32]. The sensitivity indices, SI 1st and SI total, are expressed as follows: where V i is the partial variance of the i-th parameter on output Y; V(Y) is the total unconditional output variance; S ij is the contribution to the total variance by the interactions between parameters i and j; and X~i denotes variation in all input parameters. Following previous studies [32,64], the sensitivity indices of the EFAST can be derived at the cost of nk model evaluations for robust results (n is the sample size and k is the number of input parameters). Aside from the GSA approach that is based on the SIMLAB, we also used a local SA approach, the "one-factor-at-a-time" (OAT) approach to analyze the sensitivity of SIF and GPP simulations to V cmax and light conditions. Since GOME-2 and OCO-2 SIF products are based on different wavelengths, we conducted a sensitivity analysis with different V cmax values at different wavelengths to identify the feasibility and sensitivity of satellite-derived SIF in constraining GPP simulations.
With the parameter sensitivity analysis, we selected the key parameters for GPP and SIF simulations, and inverted these parameters by minimizing the differences between the observed and modeled values. The SCOPE model employs complicated radiative transfer processes [15], which has a high demand for computational resources and thus hinders the estimation of V cmax easily. To reduce the computational burden, we also used an adaptive surrogate modeling-based optimization method [38] to optimize the key parameters of the model in simulating GPP and SIF. Although the surrogate-based method might not be able to provide exactly the same optimal solutions as the original models, it has been proven to be effective and efficient in obtaining approximate optimal parameters for land surface and radiative transfer models [36,37,62]. The surrogate-based method can acquire acceptable optimized results with fewer simulations. Following a previous study [62], we applied an adaptive, nonlinear regression method, Gaussian processes regression (GPR), to approximate and replace the original computation model with the surrogate model [38]. Previous studies [37,62] showed that the GPR method outperformed other surrogate methods. According to the results of the parameter sensitivity, we selected the key physiological parameters (i.e., V cmax , Rd, and fqe) that are related to SIF and GPP estimation in the model for optimization. The SCE method [35] was used to find the optimal parameters. Similar to the parameter sensitivity analysis, the estimation of the parameters was also based on three different objective functions using SIF and GPP to optimize parameters separately and also using both SIF and GPP to constrain parameters at the same time. We used the data from 2015 to 2017 for calibration, and the remaining data for validation. We evaluated the performance of the SCOPE model for simulating SIF using both OCO-2 and GOME-2 SIF observations.

Parameter Sensitivity Analysis of the Soil Canopy Observation, Photochemistry and Energy Fluxes (SCOPE) Model
We performed a sensitivity analysis of all parameters related to photosynthesis and SIF simulations. The first sensitivity index (SI 1st) and the total sensitivity index (SI total) of model parameters for SCOPE-simulated SIF (in OCO-2 SIF wavelengths of 757 nm and 771 nm) and GPP are illustrated in Figure 3. The results of the sensitivity analysis showed that for the SIF estimation, fqe was the most sensitive parameter, while V cmax was one of the least sensitive parameters (Figure 3a). The simulated SIF was sensitive to fqe, Cab, and LAI with the total-order impact ratio of 48.0%, 18.1%, and 16.7%, respectively. The other parameters had relatively low impact ratios (less than 5%). However, for the GPP simulation, the most sensitive parameters were V cmax , LAI, Rd, and Cab with the totalorder impact ratio of 41.7%, 31.4%, 16.1%, and 5.4%, respectively, and other parameters had relatively low impact ratio (Figure 3b). For the objective function of both SIF and GPP, both fqe and V cmax were the most sensitive parameters (Figure 3c); both GPP and SIF were sensitive to V cmax , Cab, Rd, fqe, and LAI with the total-order impact ratio of 21.9%, 21.2%, 17.7%, 15.3%, and 13.2%, respectively, while other parameters had relatively low impact ratio. Both LAI and Cab were used as input variables for the simulation of GPP and SIF and were therefore not optimized in this study.
To determine the feasibility and sensitivity of satellite measured SIF derived from different sensors (i.e., GOME-2 and OCO-2) to V cmax , we conducted a sensitivity analysis with different V cmax values ranging from 10 to 200 µmol m −2 s −1 (Figure 4). The SIF spectrum ranged from 640 to 850 nm with two peaks centered at 685 nm and 740 nm. The sensitivity of SIF to V cmax was higher at 757 nm than at 771 nm. The GOME-2 wavelength (740 nm) is right at one of the SIF emission peaks, which is also the V cmax sensitivity band for SIF prediction. The sensitivity of SIF to V cmax at the GOME-2 wavelength (740 nm) was higher than that of the OCO-2 wavelengths (757 nm and 771 nm).
Remote Sens. 2021, 13, x FOR PEER REVIEW To determine the feasibility and sensitivity of satellite measured SIF deriv different sensors (i.e., GOME-2 and OCO-2) to Vcmax, we conducted a sensitivity with different Vcmax values ranging from 10 to 200 μmol m −2 s −1 (Figure 4). The SIF s ranged from 640 to 850 nm with two peaks centered at 685 nm and 740 nm. The se of SIF to Vcmax was higher at 757 nm than at 771 nm. The GOME-2 wavelength (7 right at one of the SIF emission peaks, which is also the Vcmax sensitivity band for diction. The sensitivity of SIF to Vcmax at the GOME-2 wavelength (740 nm) wa than that of the OCO-2 wavelengths (757 nm and 771 nm). To investigate how both GPP and SIF responded to Vcmax and light condit further performed a sensitivity analysis with varying Vcmax values and light co  To determine the feasibility and sensitivity of satellite measured SIF derived from different sensors (i.e., GOME-2 and OCO-2) to Vcmax, we conducted a sensitivity analysis with different Vcmax values ranging from 10 to 200 μmol m −2 s −1 (Figure 4). The SIF spectrum ranged from 640 to 850 nm with two peaks centered at 685 nm and 740 nm. The sensitivity of SIF to Vcmax was higher at 757 nm than at 771 nm. The GOME-2 wavelength (740 nm) is right at one of the SIF emission peaks, which is also the Vcmax sensitivity band for SIF prediction. The sensitivity of SIF to Vcmax at the GOME-2 wavelength (740 nm) was higher than that of the OCO-2 wavelengths (757 nm and 771 nm). To investigate how both GPP and SIF responded to Vcmax and light conditions, we further performed a sensitivity analysis with varying Vcmax values and light conditions ( Figure 5). With varying Vcmax at the constant irradiance, SIF did not significantly increase with increasing Vcmax; for a given Vcmax value, SIF increased with increasing solar irradiance (Figure 5a). In contrast, GPP increased with increasing Vcmax, particularly at higher light conditions, and light saturation of GPP occurred under lower light conditions with lower Vcmax values (Figure 5b). To investigate how both GPP and SIF responded to V cmax and light conditions, we further performed a sensitivity analysis with varying V cmax values and light conditions ( Figure 5). With varying V cmax at the constant irradiance, SIF did not significantly increase with increasing V cmax ; for a given V cmax value, SIF increased with increasing solar irradiance (Figure 5a). In contrast, GPP increased with increasing V cmax , particularly at higher light conditions, and light saturation of GPP occurred under lower light conditions with lower V cmax values (Figure 5b).

Parameter Estimation of the SCOPE Model
The optimized values of the three most important parameters (Vcmax, Rd, and fqe) SIF and GPP simulations are given in Table 2. These optimized values were used for SCOPE model simulation at Park Falls. Figure 6 shows the performance of SCOPE estimating GPP and SIF at Park Falls, which illustrates the seasonal trajectories of G and SIF simulated by SCOPE with the objective function of SIF and GPP, respectively. T simulated GPP and SIF were validated with the flux tower-based GPP and OCO-2 data, respectively. We found that the simulated GPP was well correlated with flux tow GPP with R 2 values of 0.45 and 0.88 (RMSE = 2.11 and 1.20 gC m 2 day −1 ) at the daily ti scale for the objective function of SIF and GPP, respectively. The simulated SIF w strongly correlated with OCO-2 SIF with R 2 values of 0.77 and 0.76 for the 757 nm wa length (the RMSE were 0.14 and 0.54 W m −2 μm −1 sr −1 ), and 0.75 and 0.74 for the 771 n wavelength (the RMSE were 0.06 and 0.27 W m −2 μm −1 sr −1 ) for the objective functions SIF and GPP, respectively. Our results showed that the model performance varied w the objective function: GPP was better simulated with GPP as the objective function, a similarly, SIF was better simulated with SIF as the objective function.

Parameter Estimation of the SCOPE Model
The optimized values of the three most important parameters (V cmax , Rd, and fqe) in SIF and GPP simulations are given in Table 2. These optimized values were used for the SCOPE model simulation at Park Falls. Figure 6 shows the performance of SCOPE for estimating GPP and SIF at Park Falls, which illustrates the seasonal trajectories of GPP and SIF simulated by SCOPE with the objective function of SIF and GPP, respectively. The simulated GPP and SIF were validated with the flux tower-based GPP and OCO-2 SIF data, respectively. We found that the simulated GPP was well correlated with flux tower GPP with R 2 values of 0.45 and 0.88 (RMSE = 2.11 and 1.20 gC m 2 day −1 ) at the daily time scale for the objective function of SIF and GPP, respectively. The simulated SIF was strongly correlated with OCO-2 SIF with R 2 values of 0.77 and 0.76 for the 757 nm wavelength (the RMSE were 0.14 and 0.54 W m −2 µm −1 sr −1 ), and 0.75 and 0.74 for the 771 nm wavelength (the RMSE were 0.06 and 0.27 W m −2 µm −1 sr −1 ) for the objective functions of SIF and GPP, respectively. Our results showed that the model performance varied with the objective function: GPP was better simulated with GPP as the objective function, and similarly, SIF was better simulated with SIF as the objective function.         To further evaluate the performance of the optimized model, we validated the model optimized with both SIF and GPP as constraints independently at the Willow Creek site. We compared GPP and SIF simulated by the SCOPE model based on the default and optimized parameters against flux tower GPP and OCO-2 SIF observations at Willow Creek ( Figure 8). The model based on the optimized parameters estimated both GPP and SIF fairly well (R 2 = 0.88 and RMSE = 1.47 gC m 2 day −1 for daily GPP, R 2 = 0.70 and RMSE = 0.15 W m −2 µm −1 sr −1 for SIF at 757 nm, R 2 = 0.61 and RMSE = 0.09 W m −2 µm −1 sr −1 for SIF at 771 nm). The simulated GPP and SIF based on the optimized parameters were more consistent with flux tower GPP and OCO-2 SIF in seasonal trajectories than the simulated GPP and SIF based on the default parameters (Figure 8).
To further evaluate the performance of the optimized model, we valida optimized with both SIF and GPP as constraints independently at the Will We compared GPP and SIF simulated by the SCOPE model based on the d timized parameters against flux tower GPP and OCO-2 SIF observations at (Figure 8). The model based on the optimized parameters estimated both fairly well (R 2 = 0.88 and RMSE = 1.47 gC m 2 day −1 for daily GPP, R 2 = 0.7 0.15 W m −2 μm −1 sr −1 for SIF at 757 nm, R 2 = 0.61 and RMSE = 0.09 W m −2 μm − 771 nm). The simulated GPP and SIF based on the optimized parameters w sistent with flux tower GPP and OCO-2 SIF in seasonal trajectories than GPP and SIF based on the default parameters (Figure 8).

Predicting Long-Term SIF and GPP Using the Optimized Model
We further predicted the long-term SIF and GPP from 2003 to 2017 u mized SCOPE model. The GPP and SIF series simulated by the SCOPE mod tower GPP and satellite-based SIF data at Park Falls are illustrated in Fig  The seasonal trajectories of GPP and SIF simulated based on the optimize were compared against those based on the default parameters. The simula based on the optimized parameters was more consistent with the flux tow pared with that based on the default parameters. Similarly, the simulated SI nm, and 740 nm) based on the optimized parameters were also more consi OCO-2 and GOME-2 SIF observations compared with that based on the de ters. Figure 10 shows the scatter plots of predicted GPP against flux tower G to 2017. The optimized model had improved performance in predicting model with the default parameters at both hourly and daily timescales. W that the model with the optimized parameters had higher performance (low

Predicting Long-Term SIF and GPP Using the Optimized Model
We further predicted the long-term SIF and GPP from 2003 to 2017 using the optimized SCOPE model. The GPP and SIF series simulated by the SCOPE model against flux tower GPP and satellite-based SIF data at Park Falls are illustrated in Figures 9 and 10. The seasonal trajectories of GPP and SIF simulated based on the optimized parameters were compared against those based on the default parameters. The simulated daily GPP based on the optimized parameters was more consistent with the flux tower GPP compared with that based on the default parameters. Similarly, the simulated SIF (757 nm, 771 nm, and 740 nm) based on the optimized parameters were also more consistent with the OCO-2 and GOME-2 SIF observations compared with that based on the default parameters. Figure 10 shows the scatter plots of predicted GPP against flux tower GPP from 2003 to 2017. The optimized model had improved performance in predicting GPP than the model with the default parameters at both hourly and daily timescales. We also found that the model with the optimized parameters had higher performance (lower RMSE and higher slope) in predicting SIF at GOME-2 SIF wavelength (740 nm) over the period 2007-2017 ( Figure 11).     Remote Sens. 2021, 13, x FOR PEER REVIEW Figure 11. Comparison of simulated SIF (in GOME-2 SIF wavelength) against GOME-2 SI vations from 2007 to 2017. The black points and blue points represent the SIF simulation fault parameters and that with the optimized parameters, respectively.

Evaluating the Capability of SIF in Estimating GPP
To evaluate the relationship between GPP and SIF at different timescales, w ined the relationships between our simulated SIF and flux tower GPP at differe lengths at both instantaneous and daily timescales for the Park Falls site (Figure model predicted SIF exhibited consistently strong linear correlation with the fl GPP for different wavelengths (757 nm, 771 nm, and 740 nm) at both timescales < 0.001). Meanwhile, for each wavelength, the RMSE was lower at the daily times at the hourly timescale. Figure 11. Comparison of simulated SIF (in GOME-2 SIF wavelength) against GOME-2 SIF observations from 2007 to 2017. The black points and blue points represent the SIF simulation with default parameters and that with the optimized parameters, respectively.

Evaluating the Capability of SIF in Estimating GPP
To evaluate the relationship between GPP and SIF at different timescales, we examined the relationships between our simulated SIF and flux tower GPP at different wavelengths at both instantaneous and daily timescales for the Park Falls site (Figure 12). The model predicted SIF exhibited consistently strong linear correlation with the flux tower GPP for different wavelengths (757 nm, 771 nm, and 740 nm) at both timescales (p value < 0.001). Meanwhile, for each wavelength, the RMSE was lower at the daily timescale than at the hourly timescale.
Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 22 Figure 11. Comparison of simulated SIF (in GOME-2 SIF wavelength) against GOME-2 SIF observations from 2007 to 2017. The black points and blue points represent the SIF simulation with default parameters and that with the optimized parameters, respectively.

Evaluating the Capability of SIF in Estimating GPP
To evaluate the relationship between GPP and SIF at different timescales, we examined the relationships between our simulated SIF and flux tower GPP at different wavelengths at both instantaneous and daily timescales for the Park Falls site (Figure 12). The model predicted SIF exhibited consistently strong linear correlation with the flux tower GPP for different wavelengths (757 nm, 771 nm, and 740 nm) at both timescales (p value < 0.001). Meanwhile, for each wavelength, the RMSE was lower at the daily timescale than at the hourly timescale. We also examined how well the simulated SIF by the optimized SCOPE model could be used to estimate GPP based on the relationship between these variables (Figure 13). We also examined how well the simulated SIF by the optimized SCOPE model could be used to estimate GPP based on the relationship between these variables ( Figure 13). Flux tower GPP had strong linear relationships with the satellite retrieved SIF and MODISderived EVI (p value < 0.001) (Figure 13a-d). The OCO-2 SIF at 757 nm was more strongly correlated with tower GPP than the OCO-2 SIF at 771 nm. The OCO-2 SIF and MODIS EVI data were more strongly correlated with tower GPP than GOME-2 SIF. We further used these relationships between satellite SIF and tower based GPP (Figure 13a-c) to predict the time series of GPP from our simulated SIF time series over the period 2003-2017 for different wavelengths (Figure 13e-g). For comparison purposes, we also calculated EVI from the surface reflectance simulated by the optimized SCOPE model for the period 2003-2017 and then used the relationship between tower GPP and MODIS EVI (Figure 13d) to predict GPP from EVI for the 15-year period (Figure 13h). The SIF-GPP models exhibited a better performance than the EVI-GPP model. Furthermore, the relationship between OCO-2 SIF and GPP had a better performance in predicting GPP than that between GOME-2 SIF and GPP. Our results showed that the SIF data simulated by the optimized SCOPE model can be used to predict GPP continuously. Flux tower GPP had strong linear relationships with the satellite retrieved SIF and MODIS-derived EVI (p value < 0.001) (Figure 13a-d). The OCO-2 SIF at 757 nm was more strongly correlated with tower GPP than the OCO-2 SIF at 771 nm. The OCO-2 SIF and MODIS EVI data were more strongly correlated with tower GPP than GOME-2 SIF. We further used these relationships between satellite SIF and tower based GPP (Figure 13ac) to predict the time series of GPP from our simulated SIF time series over the period 2003-2017 for different wavelengths (Figure 13e-g). For comparison purposes, we also calculated EVI from the surface reflectance simulated by the optimized SCOPE model for the period 2003-2017 and then used the relationship between tower GPP and MODIS EVI (Figure 13d) to predict GPP from EVI for the 15-year period (Figure 13h). The SIF-GPP models exhibited a better performance than the EVI-GPP model. Furthermore, the relationship between OCO-2 SIF and GPP had a better performance in predicting GPP than that between GOME-2 SIF and GPP. Our results showed that the SIF data simulated by the optimized SCOPE model can be used to predict GPP continuously.

Parameter Sensitivity of the SCOPE Model
We first identified the sensitive parameters of the SCOPE model (v1.7) based on the OCO-2 SIF observations and flux tower GPP using both the GSA (i.e., Saltelli's method) and OAT approaches. Sensitivity analysis often serves as a screening tool to reduce the parametric dimensionality and to identify the most important parameters on model performance for parameter optimization [31]. To accurately simulate and better understand the canopy fluorescence and photosynthesis dynamics, we identified and optimized the key physiological parameters in the SCOPE model (Vcmax, Rd, fqe). SIF is linked to the maximum carboxylation capacity, Vcmax, an important parameter that determines the capacity of photosynthesis [15,42], while their complex relationship is still not well understood. It is debatable to what extent SIF is sensitive to Vcmax across biomes [41,45]. A previous study [41] showed strong linear relationship between SIF and Vcmax using the SCOPE

Parameter Sensitivity of the SCOPE Model
We first identified the sensitive parameters of the SCOPE model (v1.7) based on the OCO-2 SIF observations and flux tower GPP using both the GSA (i.e., Saltelli's method) and OAT approaches. Sensitivity analysis often serves as a screening tool to reduce the parametric dimensionality and to identify the most important parameters on model performance for parameter optimization [31]. To accurately simulate and better understand the canopy fluorescence and photosynthesis dynamics, we identified and optimized the key physiological parameters in the SCOPE model (V cmax , Rd, fqe). SIF is linked to the maximum carboxylation capacity, V cmax , an important parameter that determines the capacity of photosynthesis [15,42], while their complex relationship is still not well understood. It is debatable to what extent SIF is sensitive to V cmax across biomes [41,45]. A previous study [41] showed strong linear relationship between SIF and V cmax using the SCOPE model and estimated seasonally-varying V cmax based on the GOME-2 SIF product. In contrast, other studies have shown that the relationship between V cmax and SIF within the SCOPE model was not stable and varied between different versions [19,20,39,45]. In this study, we used Saltelli's method and OAT approach to quantify the sensitivity of SIF and GPP to key parameters in SCOPE at two temperate forest sites. Both approaches showed that SIF had low sensitivity to V cmax , while GPP had high sensitivity to this parameter. Our results are inconsistent with the finding of a previous study [41] based on the old version of the SCOPE model that there is a strong linear relationship between SIF and V cmax . Our finding is instead consistent with the previous results of SIF simulation using SCOPE v1.7 [19,44,62], showing a relatively low sensitivity of SIF to photosynthesis in the global parameter space of the model.
SIF is a signal that is mainly related to the electron transport rate in the light reactions of photosynthesis and only partly related to the carboxylation rates [3,43], while carboxylation rates mainly depend on the leaf nitrogen content and fraction of nitrogen in RuBisCO as well as the internal CO 2 concentration [65]. It is noticed that the V cmax cannot be constrained by the optical data alone, and the GPP simulation was improved by combined use of optical and thermal infrared data [66]. Previous studies have demonstrated that V cmax only partly affected SIF [16,43]; our study also confirmed that SIF was less sensitive to V cmax .
In addition, our study also suggested the SIF simulation was more sensitive to fqe, a parameter that directly influences the SIF emissions and is often based on empirical values [46], while V cmax is the most important parameter for GPP simulation. The physiological parameters such as Cab and LAI are also sensitive parameters in SIF and GPP simulations, which is similar to the finding of a previous study [62]. Meanwhile, we also found that the SIF and GPP were constrained well by the concurrent use of OCO-2 SIF and flux tower GPP data. V cmax is a parameter that directly influences the photosynthetic capacity and indirectly influences the SIF simulation.

Parameter Estimation and Model Evaluation
Given the complexity of the SCOPE model, the parameter estimation strategy is of great importance for photosynthesis and SIF simulations. Previous studies used simplified approaches (e.g., using the linear relationship between SIF and V cmax or using a module of the SCOPE model) to optimize the parameters of the SCOPE model (e.g., [19,41]). As abovementioned, we found that the relationship between SIF and V cmax was not statistically significant in version 1.70 of the SCOPE model, while the simplified SCOPE module (i.e., only using a module of the SCOPE model) cannot reflect the linkages between SIF and photosynthesis processes. Therefore, we utilized the surrogate modeling-based approach to optimize the parameters of the model because this approach requires a limited number of model runs and is more computationally efficient. The surrogate approach has been proven to be efficient and effective in searching the approximate optimal parameters in previous studies (e.g., [37,40]). Our study showed that the surrogate approach is able to estimate parameters of a relatively complex model like SCOPE. However, it should be noted that the surrogate model is not the "real" model, and it is designed to use cheaper, approximate solutions to a simulation model. Its performance depends on multiple factors such as the surrogate method choice, optimization methods, and the initial sample size and complexity of the problems [37,40].
Previous studies have used inversion approaches to retrieve key parameters in controlling photosynthetic activity (e.g., [19,26,27,44]), and have demonstrated that optimizing the biochemical parameters such as the V cmax can effectively improve the estimation of photosynthetic capacity [19]. Previous studies have mainly utilized coarse-resolution SIF observations from GOME-2 [12,13,47], and few studies have used finer-resolution OCO-2 SIF to constrain parameters in the SCOPE model (e.g., [44]). We used both OCO-2 SIF and flux tower GPP data to constrain the SIF and GPP simulation by optimizing photosynthetic parameters and SIF related parameters. The optimized SCOPE model had satisfactory performance in simulating SIF and GPP. Similar to two previous studies [44,45], our study demonstrated that the OCO-2 SIF measurements constrained the SIF related parameters such as fqe well but not V cmax . Our study demonstrated that the OCO-2 SIF observations could improve the SIF simulation of the SCOPE model by optimizing uncertain parameters. The optimized SCOPE model was evaluated using flux tower GPP data and SIF observations from both OCO-2 and GOME-2. Our results showed that the optimized model captured the magnitude, seasonality, and interannual variability of both SIF and GPP fairly well. This implies that the model can be used to effectively predict SIF and GPP. The simulated SIF and GPP can be potentially used to examine the responses of fluorescence and photosynthesis to extreme events like drought. Unlike satellite data (e.g., GOME-2, OCO-2), the SCOPE model can be to simulate SIF and GPP for the past, present, and future. In addition, the predicted SIF by the SCOPE model can potentially be used to fill the data gaps for the OCO-2 SIF products to generate spatially and temporally continuous SIF products over the globe. The uncertainty/bias in the meteorological data used by the SCOPE model will certainly introduce uncertainty to the SIF estimates. The model can be run without the satellite input, which makes it possible to scale from instantaneous data into diurnal cycles, and calculate the SIF and carbon fluxes in cloudy and rainy days [15]. The produced SIF can also be further used to improve carbon cycle modeling at large scales based on the SIF-GPP relationships.
Previous studies showed significant relationships between satellite-derived SIF measurements (GOSAT, GOME-2, and OCO-2) and gridded GPP data (e.g., [5,12,13]). The relationship between SIF and GPP is sometimes complex and ecosystem-specific [4], and SIF contains information on canopy physiological processes [67]. The ground area of OCO-2 SIF soundings is much closer to the flux tower footprint, and several recent studies showed linear relationship between the OCO-2 SIF and flux tower GPP (e.g., [2,5,14,44]). Meanwhile, modeling studies have also demonstrated the linear SIF-GPP relationship using processed-based models like SCOPE. Similar to [44] and [20], we found strong linear correlation between SIF and GPP from both observations and models. Our results also showed that the SIF simulated by our optimized SCOPE model could be used to estimate GPP fairly well based on the SIF-GPP relationship established from satellite-derived SIF data and flux tower GPP. The long-term SIF simulations and the resulting GPP estimates can be used to examine the dynamics of terrestrial photosynthesis at various timescales.

Uncertainties and Future Perspectives
The default values of the parameters in terrestrial ecosystem models are typically based on field measurements, literature, assumptions, or calibrations [15,16], and some parameters are associated with substantial uncertainty. The uncertainty in these parameters can lead to significant uncertainty in modeled fluxes [27,29]. Evaluating the uncertainty of simulated fluorescence and photosynthesis is a persistent challenge. Several important factors can affect the SIF and GPP model performance. Aside from the uncertainty of model parameters, the uncertainty of flux tower GPP estimation (e.g., [27,68], SIF retrieval uncertainty (e.g., [12,13]), the mismatch of the spatial scales between satellite observation (e.g., the coarse resolution of GOME-2 SIF measurements), and flux tower footprint can also introduce uncertainty to the simulated SIF and GPP. For example, different NEE partitioning approaches will induce uncertainty in GPP estimates [68]. Our study found that the model parameters optimized with OCO-2 SIF were able to better estimate SIF than those optimized with GOME-2 SIF since the use of OCO-2 SIF observations can substantially reduce the scale mismatch between satellite grid cells and the flux tower footprint than the use of GOME-2 SIF data.
In this study, we only assessed the capacity of the OCO-2 SIF measurements and flux tower data in constraining the photosynthesis and SIF parameters within the SCOPE model with version 1.7. The recently released SCOPE 2.0 is likely to improve the modeling of SIF and photosynthesis. In addition, new SIF observations from the TROPOspheric Monitoring Instrument (TROPOMI), which provides global SIF measurements on a daily basis, might be valuable for constraining parameters because of its temporal continuity [69]. However, it should also be recognized that its grid cell (3.5~15 km × 7 km) is larger than the typical footprint of flux towers, leading to a sizeable scale mismatch. The next generation SIF mission, the FLEX ( [70]), will provide SIF data with a spectral range between 500 and 789 nm, temporal continuity, and small grid cells (300 m × 300 m), and these data will have greater potential for SIF and GPP modeling and parameter estimation than the currently available SIF observations. Additionally, the SCOPE model is based on a 1-D radiative transfer model (i.e., PROSAIL), and the model is thus limited to homogeneous vegetation canopies, and cannot reflect the inhomogeneity in vegetation properties of complex canopy (e.g., LAI in horizontal direction) well. The 3-D radiative transfer model (e.g., DART) can be expected to be used in the future for the SIF estimation of the vegetation with more complex canopies.

Conclusions
A SCOPE modeling study was conducted to examine how well OCO-2 SIF observations and flux tower GPP data could constrain uncertain parameters and thereby improve the capability of the model in simulating both SIF and GPP at two temperate forest sites. The sensitivity of GPP simulation to V cmax was greater than that of SIF simulation in SCOPE. The performances of SIF and GPP estimation were improved with the inverse estimation of key parameters in the model. Incorporating the OCO-2 SIF data constrained key uncertain parameters related to SIF simulation, but did not significantly improve the GPP simulation, suggesting that SIF observations have low capability in constraining photosynthetic capacity. The concurrent use of satellite SIF observations and flux tower GPP data improved the simulations of both SIF and GPP. The optimized model can be used to simulate SIF and GPP continuously for the past, present, and future. The simulated SIF can also potentially be used to fill the data gaps for satellite SIF products like OCO-2 to produce spatially and temporally continuous SIF products and to improve the diagnosis of carbon cycle at large scales.