Improved Estimation of the Gross Primary Production of Europe by Considering the Spatial and Temporal Changes in Photosynthetic Capacity from 2001 to 2016

Abstract

The value of leaf photosynthetic capacity (V_cmax) varies with time and space, but state-of-the-art terrestrial biosphere models rarely include such V_cmax variability, hindering the accuracy of carbon cycle estimations on a large scale. In particular, while the European terrestrial ecosystem is particularly sensitive to climate change, current estimates of gross primary production (GPP) in Europe are subject to significant uncertainties (2.5 to 8.7 Pg C yr⁻¹). This study applied a process-based Farquhar GPP model (FGM) to improve GPP estimation by introducing a spatially and temporally explicit V_cmax derived from the satellite-based leaf chlorophyll content (LCC) on two scales: across multiple eddy covariance tower sites and on the regional scale. Across the 19 EuroFLUX sites selected for independent model validation based on 9 plant functional types (PFTs), relative to the biome-specific V_cmax, the inclusion of the LCC-derived V_cmax improved the model estimates of GPP, with the coefficient of determination (R²) increased by 23% and the root mean square error (RMSE) decreased by 25%. V_cmax values are typically parameterized with PFT-specific V_cmax calibrated from flux tower observations or empirical V_cmax based on the TRY database (which includes 723 data points derived from V_cmax field measurements). On the regional scale, compared with GPP, using the LCC-derived V_cmax, the conventional method of fixing V_cmax using the calibrated V_cmax or TRY-based V_cmax overestimated the annual GPP of Europe by 0.5 to 2.9 Pg C yr⁻¹ or 5 to 31% and overestimated the interannually increasing GPP trend by 0.007 to 0.01 Pg C yr⁻² or 14 to 20%, respectively. The spatial pattern and interannual change trend of the European GPP estimated by the improved FGM showed general consistency with the existing studies, while our estimates indicated that the European terrestrial ecosystem (including part of Russia) had higher carbon assimilation potential (9.4 Pg C yr⁻¹). Our study highlighted the urgent need to develop spatially and temporally consistent V_cmax products with a high accuracy so as to reduce uncertainties in global carbon modeling and improve our understanding of how terrestrial ecosystems respond to climate change.

1. Introduction

The terrestrial ecosystem offsets approximately 30.5% of the carbon dioxide (CO₂) released due to anthropogenic activity [1] and plays a prominent role in regulating global carbon cycling [2]. As a quantitative indicator of the total amount of carbon assimilated via photosynthesis, terrestrial gross primary production (GPP) serves as the initial driver of the global carbon cycle [3]. On the continental scale, the terrestrial ecosystems in Europe have been proven to be ecologically fragile and particularly sensitive to climate change [4,5,6]. According to long-term-recorded remotely sensed data, the terrestrial ecosystem has demonstrated widespread greening since the 1980s, especially in the Northern Hemisphere [7]. Modeling studies have suggested that the total GPP values for Europe vary in a wide range from 2.5 to 8.7 Pg C yr⁻¹ [8,9,10,11,12,13,14]. While there are some discrepancies between study regions (e.g., those including or excluding part of Russia), estimates of European GPP are subject to significant uncertainties, hindering our understanding of the role of the European terrestrial ecosystems in mitigating climate change.

Attempts to model GPP on a global scale using remote sensing data fall into three general categories. The first is empirical approaches based on statistical models or machine learning. These studies empirically link GPP with the spectral vegetation index [15,16,17], leaf area index (LAI) [18], and sun-induced chlorophyll fluorescence (SIF) [19,20]. Some studies use machine learning methods to estimate local or global GPP [21,22]. The second category is the widely used light use efficiency (LUE) model [23]. These models assume that GPP is a product of the fraction of absorbed photosynthetic active radiation (APAR) and LUE reduced by modifying factors. Examples are BIOMASS [24], CASA [25], C-Fix [8], 3-PG [26], VPM [27], EC-LUE [28], the P-model [29], and CCW [30] models. The third category is process-based terrestrial biosphere models (TBMs), such as CENTURY [31], TEM [24], Biome-BGC [32], BESS [33], BEPS [34], and FGM [35]. Biologically, GPP is a product of leaf-scale photosynthesis. Thus, GPP is related to both internal and environmental factors, including rapid leaf-level biochemical reactions, stomatal conductance [36], canopy structure [37], climatic factors, soil moisture [38], and slower environmental acclimation processes [39]. However, empirical models, machine learning methods, and LUE models have a limited ability to simulate the response of GPP to complicated environmental and internal biological factors due to their inadequate representation of the mechanisms that regulate the physiological process of photosynthesis.

TBMs have proven to be particularly useful for estimating GPP due to their inclusion of the biochemical processes of photosynthesis. TBMs commonly include the mechanistic leaf photosynthesis model developed by Farquhar et al. (1980). To estimate the spatiotemporal patterns of GPP on a large scale, TBMs require a range of forcing data, such as meteorological data, land cover, the leaf area index, the clumping index, and leaf trait information. Most importantly, the leaf photosynthesis rate simulated by the Farquhar model is particularly sensitive to the parameterization of leaf photosynthetic capacity at 25 °C (V_cmax) ($\mathsf{\mu }\mathrm{mol}$ CO₂ m⁻² s⁻¹) [2,40,41], reflecting the active amount and kinetic activity of the Rubisco enzyme in leaves. Inadequate constraints on V_cmax lead to substantial uncertainties in GPP estimation [42]. Traditionally, V_cmax is estimated by measuring the net photosynthesis rate (A_n) relative to internal CO₂ pressure (C_i) (i.e., A_n–C_i curve) at different CO₂ concentrations with saturating irradiance or using a modified ‘one-point method’ [43,44]. Measuring one A_n–Ci curve can take up to one hour, generating only one V_cmax value based on a small number of leaf samples. Thus, field measurements of V_cmax on the leaf scale [45] are laborious, time-consuming, and, most importantly, mismatched with the footprints (~100 m–450 m) of eddy covariance (EC) flux towers [46].

Due to the lack of spatiotemporal information on V_cmax, the state-of-the-art TBMs generally assume a constant V_cmax for a specific plant functional type (PFT). More specifically, PFT-specific V_cmax values are typically parameterized using two types of data: (1) optimal V_cmax data calibrated from a GPP derived from eddy covariance flux tower measurements [35,47,48]; and (2) empirical V_cmax compiled from field measurements reported in the literature [49]. However, many studies have proven the existence of variations in V_cmax with space and time, even for plants with the same PFT [50,51,52,53,54]. Consequently, the conventional parameterization methods that fix V_cmax as a PFT-specific constant can lead to substantial bias in GPP estimations. In addition, previous studies have mostly been conducted on the global scale, and TBMs with parameters calibrated based on all flux tower observations worldwide may have limited accuracy on the regional scale, as in Europe. Thus, European GPP could possibly be improved by including spatially and temporally explicit V_cmax values combined with a model calibration method based on EC observations based in Europe alone.

Recent advances in remote sensing have made it possible to derive dynamic V_cmax information on a global scale [55]. Global V_cmax products are generally estimated based on the strong correlation between V_cmax and two major biochemistry properties (i.e., the leaf chlorophyll content (LCC) and leaf nitrogen content) [56], which can be estimated from remotely sensed hyperspectral land surface reflectance and solar-induced chlorophyll fluorescence (SIF) data. For example, two V_cmax products are estimated from GOME-2 SIF data [57] and GOME-2/OCO-2 SIF data [51]. However, the spatial resolution of SIF-derived V_cmax products is generally too coarse (i.e., 36 km–1°) for regional studies. In contrast, the V_cmax products derived from LCC have a higher spatial resolution (i.e., 500 m to 1$°$) [58,59]. While some LCC-based V_cmax products have a relatively low update frequency and a short accumulation time (i.e., less than a decade with a one-month interval) [58], a few LCC-based V_cmax products provide approximately twenty years of comprehensive V_cmax estimation at a 500 m spatial resolution and 8-day temporal resolution [59]. Thus, the newly developed remote sensing V_cmax products provide an opportunity to improve the estimation of GPP in Europe by including spatial and temporal variations in V_cmax.

While a high value has been placed on the modeling of V_cmax dynamics on a global scale, the concomitant increase in our understanding of the V_cmax change effect on GPP has only partially been realized. This brings us to the crux of our study: the quantitative analysis of uncertainties in GPP using the conventional constant V_cmax parameterization in TBMs. In this study, we hypothesize that by considering changes in V_cmax, we can improve the estimation of spatial and temporal variations in GPP. Specifically, we address the following three scientific questions: (1) How much carbon has been assimilated by the terrestrial ecosystem in Europe? (2) Can GPP estimation be improved by including the spatiotemporal dynamics of V_cmax compared with the conventional method of fixing V_cmax as a PFT-specific constant? (3) How much uncertainty lies in European GPP estimates that do not consider changes in V_cmax? By answering these questions, our study offers an improved estimation of the carbon assimilated by the terrestrial ecosystem in Europe and can help us to better understand the role of the Northern Hemisphere in mitigating climate change.

2. Materials and Methods

2.1. Study Regions and Flux Towers

The study region covers the mainland of Europe and part of Russia, excluding England and parts of Siberia. Flux towers provide direct measurements of ecosystem carbon fluxes. In this study, 40 sites located in Europe (Figure 1) were selected from FLUXNET 2015 (https://fluxnet.org/data/fluxnet2015-dataset/) [60] based on the availability of V_cmax data. In addition, we screened out EuroFLUX sites with inconsistent profiles of LAI and GPP derived from eddy covariance data. Given the limited pool of shrubland sites in EuroFLUX, the SH site from AmeriFLUX was also included, since only one CSH site (i.e., US-KS2) is available in AmeriFLUX. The land cover types in the study region were extracted from the MODIS International Geosphere–Biosphere Programme (IGBP) classification product (Figure 1), which has a spatial resolution of 500 m. The land cover types were grouped into a total of ten PFTs, including croplands (CRO), closed shrublands (CSH), deciduous broadleaf forest (DBF), deciduous needleleaf forest (DNF), evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), grasslands (GRA), mixed forest (MF), open shrublands (OSH), and wetland (WET).

2.2. Methods

2.2.1. A Process-Based Farquhar GPP Model (FGM)

We recently developed a large-scale Farquhar GPP model (FGM) based on eddy covariance data and remote sensing data [35]. The FGM model was initially developed from a stand-level GPP model based on Song et al. (2009). Derived from the Song et al. (2009) model, the FGM estimates GPP by integrating the Farquhar leaf-level biochemical photosynthesis model [40] with a two-leaf radiation interception simulation method. In the Song et al. (2009) model, GPP is solved based on three complex equations: Fick’s law, the Farquhar photosynthesis model, and a model for stomatal conductance. This approach is computationally expensive. To reduce the computational need for large-scale GPP estimation at a high spatiotemporal resolution, we introduced the optimal stomatal conductance theory to compute GPP more efficiently with the FGM.

The FGM simulates carbon assimilation using the Farquhar, von Caemmerer, and Berry (i.e., FvCB) [40] enzyme kinetic model, which couples electron transport and the Calvin-Benson cycle. The function of the leaf photosynthetic rate takes the minimum of the V_cmax-limited photosynthesis rate (i.e., ${\mathrm{A}}_{\mathrm{v}}$) and light-limited photosynthesis rate (i.e., ${\mathrm{A}}_{\mathrm{j}}$). Some fundamental equations for GPP estimation with the FGM are described here:

$${\mathrm{A}}_{\mathrm{n}}=\min \left\{\begin{array}{c}{\mathrm{A}}_{\mathrm{v}}\\ {\mathrm{A}}_{\mathrm{j}}\end{array}\right\}$$

(1)

$${\mathrm{A}}_{\mathrm{v}}=\frac{{\mathrm{V}}_{\mathrm{cmax}}{(\mathrm{C}}_{\mathrm{i}}{-\mathsf{\Gamma }}^{*})}{{\mathrm{C}}_{\mathrm{i}}{+\mathrm{K}}_{\mathrm{C}}{(1+\mathrm{O}/\mathrm{K}}_{\mathrm{O}})}{-\mathrm{R}}_{\mathrm{d}}$$

(2)

$${\mathrm{A}}_{\mathrm{j}}=\frac{{\mathrm{J}(\mathrm{C}}_{\mathrm{i}}{-\mathsf{\Gamma }}^{*})}{4{.5\mathrm{C}}_{\mathrm{i}}+10{.5\mathsf{\Gamma }}^{*}}{-\mathrm{R}}_{\mathrm{d}}$$

(3)

where ${\mathrm{A}}_{\mathrm{n}}$ is the net photosynthesis rate; ${\mathrm{A}}_{\mathrm{v}}$ is the minimum of the ${\mathrm{V}}_{\mathrm{cmax}}$-limited photosynthesis rate; ${\mathrm{A}}_{\mathrm{j}}$ is the light-limited photosynthesis rate; R_d is the dark respiration rate; ${\mathsf{\Gamma }}^{*}$ is the CO₂ compensation point; ${\mathrm{K}}_{\mathrm{C}}$ and ${\mathrm{K}}_{\mathrm{O}}$ are the Michaelis–Menten constants of Rubisco for CO₂ and O₂, respectively; $\mathrm{O}$ is the intercellular oxygen partial pressure in the leaves; and J is the rate of electron transport.

According to the optimal stomatal conductance theory, plants adjust their stomata to minimize the combined unit costs of transpiration and carbon assimilation [61]. Assuming the residual conductance parameter ${\mathrm{g}}_0$ equals zero, the ratio of the intercellular CO₂ concentration (${\mathrm{C}}_{\mathrm{i}}$) to ambient CO₂ concentration (${\mathrm{C}}_{\mathrm{a}}$) is regulated by the atmospheric vapor pressure deficit based on optimal stomatal theory [62,63]:

$$\frac{{\mathrm{C}}_{\mathrm{i}}}{{\mathrm{C}}_{\mathrm{a}}} \approx \frac{{\mathrm{g}}_1}{{\mathrm{g}}_1+\sqrt{\mathrm{D}}}$$

(4)

where $\mathrm{D}$ is the vapor pressure deficit in kPa, and ${\mathrm{g}}_1$ is an empirical parameter in kPa^0.5. According to its theoretical interpretation, the parameter ${\mathrm{g}}_1$ increases with the marginal water cost of carbon $\mathsf{\lambda }$ and the CO₂ compensation point ${\mathsf{\Gamma }}^{*}$. Thus, species with a high ${\mathrm{g}}_1$ will have a low instantaneous water use efficiency, i.e., a lower ratio of photosynthesis to the transpiration rate.

The FGM estimates the mean carboxylation capacity of a unit sunlit leaf (${\mathrm{V}}_{\mathrm{cmax}25\_\mathrm{sunlit}}$) and shaded leaf (${\mathrm{V}}_{\mathrm{cmax}25\_\mathrm{shaded}}$) area with the following models [53,64,65]:

$${\mathrm{V}}_{\mathrm{cmax}25\_\mathrm{sunlit}}=\frac{{\mathsf{\Omega }\mathrm{LV}}_{\mathrm{cmax}25} \left(1.0-\exp \left(-{\mathrm{k}}_{\mathrm{n}}-{\mathrm{K}}_{\mathrm{b}}\left({\mathsf{\theta }}_{\mathrm{z}}\right)\mathsf{\Omega }\mathrm{L}\right)\right)}{{(\mathrm{k}}_{\mathrm{n}}{+\mathrm{K}}_{\mathrm{b}}\left({\mathsf{\theta }}_{\mathrm{z}}\right)\mathsf{\Omega }\mathrm{L})}{/\mathrm{L}}_{\mathrm{sunlit}}$$

(5)

$${\mathrm{V}}_{\mathrm{cmax}25\_\mathrm{shaded}}=\frac{{\mathsf{\Omega }\mathrm{LV}}_{\mathrm{cmax}25}}{{\mathrm{L}}_{\mathrm{shaded}}}\left[\frac{1.0-\exp \left(-{\mathrm{k}}_{\mathrm{n}}\right)}{{\mathrm{k}}_{\mathrm{n}}}-\frac{1.0-\exp \left(-{\mathrm{k}}_{\mathrm{n}}-\left({\mathsf{\theta }}_{\mathrm{z}}\right)\mathsf{\Omega }\mathrm{L}\right)}{{(\mathrm{k}}_{\mathrm{n}}{+\mathrm{K}}_{\mathrm{b}}\left({\mathsf{\theta }}_{\mathrm{z}}\right)\mathsf{\Omega }\mathrm{L})}\right]$$

(6)

where $\mathsf{\Omega }$ is the clumping index; L is the total LAI; ${\mathrm{L}}_{\mathrm{sunlit}}$ and ${\mathrm{L}}_{\mathrm{shaded}}$ are the LAIs for sunlit leaves and shaded leaves, respectively; ${\mathrm{k}}_{\mathrm{n}}$ is the coefficient of leaf nitrogen allocation; ${\mathsf{\theta }}_{\mathrm{z}}$ is the solar zenith angle; ${\mathrm{K}}_{\mathrm{b}}\left({\mathsf{\theta }}_{\mathrm{z}}\right)$ is the light extinction coefficient; and ${\mathrm{V}}_{\mathrm{cmax}25}$ is the maximum carboxylation rate standardized to 25 °C for sunlit leaves. ${\mathrm{L}}_{\mathrm{sunlit}}$ and ${\mathrm{L}}_{\mathrm{shaded}}$ are estimated using Beer’s law as follows:

$${\mathrm{L}}_{\mathrm{sunlit}}=\frac{-{\mathrm{e}}^{-{\mathrm{K}}_{\mathrm{b}}{(\mathsf{\theta }}_{\mathrm{z}})\mathsf{\Omega }\mathrm{L}}}{{\mathrm{K}}_{\mathrm{b}}{(\mathsf{\theta }}_{\mathrm{z}})}$$

(7)

$${\mathrm{L}}_{\mathrm{shaded}}=\mathrm{L}-{\mathrm{L}}_{\mathrm{sunlit}}$$

(8)

In this study, V_cmax represents the maximum carboxylation rates, standardized to $25 °\mathrm{C}$ hereafter (i.e., V_cmax25). In the FGM, the temperature effect on V_cmax is modeled as follows:

$${\mathrm{V}}_{\mathrm{cmax}}=\frac{{\mathrm{V}}_{\mathrm{cmax}25} \exp \left(\mathrm{a}1\right(\mathrm{T} 25\left)\right)}{(1+\exp (\mathrm{a}2(\mathrm{T}-41\left)\right)}$$

(9)

where T is air temperature and a1 (0.051) and a2 (0.205) are empirical parameters based on measurements [65,66].

The canopy total GPP is the sum of the GPP for both the sunlit and the shaded leaves:

$${\mathrm{GPP}=\mathrm{A}}_{{\mathrm{n}}_{\mathrm{sunlit}}}{\mathrm{L}}_{\mathrm{sunlit}}{+\mathrm{A}}_{{\mathrm{n}}_{\mathrm{shade}}}{\mathrm{L}}_{\mathrm{shade}}$$

(10)

where ${\mathrm{A}}_{{\mathrm{n}}_{\mathrm{sunlit}}}$ and ${\mathrm{A}}_{{\mathrm{n}}_{\mathrm{shade}}}$ are the net photosynthesis rate for a unit sunlit leaf area index and a unit shaded leaf area index, respectively.

Additional details of the theoretical framework, default parameter values for the FGM, and model calibration method can be found in our previous study [35].

2.2.2. Model Calibration and Validation Methods

In the original version of the FGM, both ${\mathrm{g}}_1$ and ${\mathrm{V}}_{\mathrm{cmax}}$ were biome-specific and calibrated together using the FLUXNET2015 dataset [60]. In this study, because a new parameterization method for ${\mathrm{V}}_{\mathrm{cmax}}$ was adopted, we only needed to calibrate the parameter ${\mathrm{g}}_1$ for various PFTs in the FGM. Here, we randomly split the sites for each PFT at a 1:1 ratio for independent model calibration and validation (

Table 1. Sites for model calibration and validation and the calibrated ${\mathrm{g}}_1$ for nine plant functional types (PFTs).

PFTs	${\mathbf{g}}_1$	$\mathbf{Sites} \mathbf{for} \mathbf{Calibration} (\mathbf{Latitude}°$$, \mathbf{Longitude}°)$	$\mathbf{Sites} \mathbf{for} \mathbf{Validation} (\mathbf{Latitude}°$$, \mathbf{Longitude}°)$
CSH 1	1.14	US-KS2 (28.61, −80.67) [60,67]	RU-Vrk (67.05, 62.94)
CRO	10	IT-BCi (40.52, 14.96) [68] IT-Cas (45.07, 8.72) DE-Seh (50.87, 6.45) [69] DE-Geb (51.10, 10.91) [70]	FR-Gri (48.84, 1.95) [71] BE-Lon (50.55, 4.75) [72] DE-Kli (50.89, 13.52) [73]
DBF	1.66	IT-Col (41.85, 13.59) [74] FR-Fon (48.48, 2.78) [75] DE-Lnf (51.33, 10.37) [70]	IT-Ro1 (42.41, 11.93) [76] IT-Ro2 (42.39, 11.92) [77] DE-Hai (51.08, 10.45) [78]
ENF	0.62	IT-Ren (46.59, 11.43) [79] CZ-BK1 (49.50, 18.54) [80] NL-Loo (52.17, 5.74) [81] RU-Fyo (56.46, 32.92) [82]	IT-Lav (45.96, 11.28) [83] CH-Dav (46.82, 9.86) [84] DE-Obe (50.79, 13.72) DE-Tha (50.96, 13.57) [85]
EBF	0.62	FR-Pue (43.74, 3.60) [86]	IT-Cpz (41.71, 12.38) [87]
GRA	1.14	IT-Tor (45.84, 7.58) [88] CH-Cha (47.21, 8.41) [89] CH-Oe1 (47.29, 7.73) [90] CZ-BK2 (49.49, 18.54)	IT-MBo (46.01, 11.05) [91] CH-Fru (47.12, 8.54) [92] DE-Gri (50.95, 13.51) [73]
MF	0.62	CH-Lae (47.48, 8.36) [93] BE-Bra (51.03, 6.0) [94]	BE-Vie (50.30, 6.0) [95]
OSH	10	ES-LgS in 2008 (36.93, −2.75) [96]	ES-Lgs in 2007 (36.93, −2.75) [96]
WET	0.62	CZ-wet (49.02, 14.77) [97] DE-Spw (51.89, 14.03) DE-Zrk (53.88, 12.89) [98]	DE-SfN (47.81, 11.33) [99] DE-Akm (53.87, 13.68)

¹ Abbreviations: closed shrublands (CSH), croplands (CRO), deciduous broadleaf forest (DBF), evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), grasslands (GRA), mixed forest (MF), open shrublands (OSH), and wetland (WET). Site locations are shown in Figure 1. Please refer to the official website (https://fluxnet.org/sites/site-list-and-pages/, accessed on 19 February 2023) for more detailed descriptions.

). We calibrated the parameter ${\mathrm{g}}_1$ using data collected at calibration sites for nine PFTs, and the remaining sites were reserved for independent validation. For most PFTs, the sites split for calibration covered a wide spatial distribution and ensured representative EC data for the calibration. Moreover, the independent validation reduced the number of uncertainties caused by the inclusion of prior information using reduplicate sites for the calibration and validation. The independent model calibration and validation ensured a reliable and objective assessment of the model’s performance.

Given that only one OSH site (i.e., ES-Lgs) is available in Europe, the flux data for ES-Lgs were separated by year for the model calibration and validation, respectively. When calibrating the values of ${\mathrm{g}}_1$ for different PFTs, V_cmax was set to the seasonal dynamics derived from the LCC. Table 1 lists the results of the calibrated ${\mathrm{g}}_1$ values for different PFTs. There are no flux towers for deciduous needleleaf forests (DNF) in Europe. Thus, we adopted the default values of ${\mathrm{g}}_1$ for DNF [35].

2.2.3. Simulation Experiments

Four simulation scenarios were designed depending on the types of V_cmax used in the FGM. First, we defined a reference setup that included changes in all the factors (i.e., simulation “All”) to drive the FGM in a straightforward manner. In the “All” simulation, the FGM was operated with the spatially and temporally explicit V_cmax derived from the LCC in a 500 m grid cell and 8-day temporal interval from 2001 to 2016 (Table 2). By considering the spatial and temporal changes in V_cmax, the European GPP for the contemporary climate was estimated based on the “All” simulations.

While the other input data were kept fixed for all the runs, three other simulation scenarios (Table 2) were designed with different types of photosynthetic capacity (PC) parameterization methods, as follows: (a) In simulation “PC1”, the FGM was operated with 500 m, 8-day V_cmax data derived from the LCC in 2001 by considering spatial and seasonal changes in V_cmax, neglecting the interannual variations in V_cmax from 2001 to 2016. (b) In simulation “PC2”, the FGM was parameterized with PFT-specific V_cmax constants retrieved from flux tower observations. This typical V_cmax parameterization method is widely adopted for TBMs. (c) In simulation “PC3”, the FGM was parameterized with PFT-specific V_cmax constants provided by the TRY database (which includes 723 data points of V_cmax field measurements) [49]. Kattge et al. (2009) compiled data on qualitative traits, climate, and soil to subdivide terrestrial vegetation into PFTs and set V_cmax to different empirical values for different PFTs.

We then compared the simulations “PC1”, “PC2”, and “PC3” with the reference to quantify the effects of the changed photosynthesis capacity on the magnitude and spatial pattern of, as well as the temporal variation in, GPP. Thus, the simulation difference between “All” and “PC1” (i.e., “All”−“PC1”) represents the impacts of interannual changes in V_cmax on GPP. The simulation difference between “All” and “PC2” (i.e., “All”−“PC2”) or “PC3” (i.e., “All”–“PC3”) represents the uncertainties regarding GPP using two typical parameterization methods by fixing V_cmax as a PFT-specific constant.

Table 2. Scenario designs used to quantify the effects of changes in photosynthesis capacity (PC) on GPP based on the FGM. The symbol ‘△’ indicates that the input variable changes over time, while the symbol ‘▲’ indicates that the seasonality of V_cmax on a large scale is included. The symbol ‘$◆$’ indicates that the input variable is fixed at a biome-specific constant.

Simulation	All (LCC-Derived Vcmax)	PC1 (LCC-Derived Vcmax)	PC2 (PFT-Specific Vcmax)	PC3 (PFT-Specific Vcmax)
LULC	△	△	△	△
LAI	△	△	△	△
Ω	△	△	△	△
Vcmax	△ (8-day, 2001 to 2016)	▲ (8-day, 2001)	$◆$ (Optimal constants retrieved from eddy covariance data)	$◆$ (Empirical constants based on TRY 1 database)
DSR	△	△	△	△
Ta	△	△	△	△
RH	△	△	△	△
CO2	△	△	△	△

¹ See Kattge et al. (2009).

Table 2. Scenario designs used to quantify the effects of changes in photosynthesis capacity (PC) on GPP based on the FGM. The symbol ‘△’ indicates that the input variable changes over time, while the symbol ‘▲’ indicates that the seasonality of V_cmax on a large scale is included. The symbol ‘$◆$’ indicates that the input variable is fixed at a biome-specific constant.

Simulation	All (LCC-Derived Vcmax)	PC1 (LCC-Derived Vcmax)	PC2 (PFT-Specific Vcmax)	PC3 (PFT-Specific Vcmax)
LULC	△	△	△	△
LAI	△	△	△	△
Ω	△	△	△	△
Vcmax	△ (8-day, 2001 to 2016)	▲ (8-day, 2001)	$◆$ (Optimal constants retrieved from eddy covariance data)	$◆$ (Empirical constants based on TRY 1 database)
DSR	△	△	△	△
Ta	△	△	△	△
RH	△	△	△	△
CO2	△	△	△	△

¹ See Kattge et al. (2009).

2.2.4. Data Analysis Methods

(1)

Model performance evaluation during the calibration process

We conducted Monte Carlo simulations to calibrate parameter g1 and V_cmax for each biome. The optimal values of g1 for different biomes were set as driving parameters for the FGM, while the optimal values of V_cmax for different biomes were only used for the quantification of the V_cmax change effect on GPP. The Monte Carlo simulations identified optimal parameters when the simulated data reached the strongest agreement with the observed data using a weighted R² (wR²) as a performance indicator [30,100]:

$${\mathrm{wR}}^2=\{\begin{array}{c}\left|\mathrm{b}\right|{\mathrm{R}}^2, \mathrm{b}<1\\ {\left|\mathrm{b}\right|}^{-1}{\mathrm{R}}^2, \mathrm{b}\ge 1\end{array}$$

(11)

where b and R² are the slope and coefficient of determination for the regression of the modeled GPP and GPP estimates derived from eddy covariance data (EC-GPP) when the intercept is forced to zero, respectively [101]. A coefficient of determination (R²) equal to one and b equal to one (wR² = 1) indicate a perfect model performance. The range of R² is 0 to 1, which describes the proportion of the observed dispersion that is explained by the prediction.

(2)

Quantification of the accuracy of the FGM GPP

To quantify the accuracy of the FGM GPP directly, first, we compared the simulated GPP with the EC-GPP collected during the daytime at the 19 validation sites. The R² and root mean square error (RMSE) were estimated using the regression lines between the modeled GPP and EC-GPP to evaluate the FGM model performance.

Second, we quantified the accuracy of the spatial pattern of European GPP predicted by the FGM. The multi-year mean of the European vegetation photosynthesis rate was calculated based on different GPP products and the remotely sensed sun-induced chlorophyll fluorescence (SIF) data for the same study area. The results of the multi-year means of the annual European GPP and SIF data with coarse spatial resolutions were resampled to the target spatial resolution of 500 m by nearest neighbor interpolation. We calculated the spatial correlation matrix between the different GPP products and the SIF data. The correlation matrix provides the correlation coefficients between each combination of two inputs using Person’s correlation (r) metric as an indicator. It is calculated as:

$${\mathrm{r}}_{\mathrm{X},\mathrm{Y}}=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{X}}_{\mathrm{i}}- \overline{\mathrm{X}}\right)\left({\mathrm{Y}}_{\mathrm{i}}- \overline{\mathrm{Y}}\right)}{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{X}}_{\mathrm{i}}- \overline{\mathrm{X}}\right)}^2}\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Y}}_{\mathrm{i}}- \overline{\mathrm{Y}}\right)}^2}}$$

(12)

where n is the number of vegetated pixels in the study area, i is the grid cell index, X_i is the estimates of GPP based on the FGM, and Y_i is the estimates of GPP based on other GPP products or the SIF data.

Third, we quantified the accuracy of the interannual changes in the FGM GPP. We calculated the annual GPP of Europe by accumulating 8-day GPP predictions with a yearly temporal resolution from 2001 to 2016 based on the FGM and other GPP products. During the study period, the annual GPP of Europe predicted by the FGM was evaluated against other global GPP products using the correlation coefficient r and interannual trend b as two quantitative indicators. Here, r describes the temporal correlation coefficient between two GPP products based on the regression of interannual GPP dynamics from 2001 to 2016, in turn based on the FGM GPP and other GPP products, and b is the slope (Pg C yr⁻²) for the regression of interannual GPP dynamics (Pg C yr⁻¹) from 2001 to 2016.

(3)

Quantification of the V_cmax change effects on GPP

On a large scale, we quantify the V_cmax change effect on GPP by measuring the magnitude in percent as the mean absolute difference between the pixel-based means of the simulations, on the one hand, and “PC2” and “PC3”, on the other, relative to the mean of the reference simulation with “All”, as:

$${\mathrm{Effect}}_{\mathrm{Magnitude}}=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}|\overline{\mathrm{AR}}-\overline{\mathrm{REF}}|}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\overline{{\mathrm{REF}}_{\mathrm{i}}}} \times 100$$

(13)

where REF is the reference modeling setup using the LCC-derived V_cmax, and AR is an alternative realization where the V_cmax of the reference setup has changed. The single overbar denotes the grid-cell-based temporal mean.

2.3. Data

2.3.1. Flux Data

FLUXNET2015 provides gap-filled EC-GPP and corresponding meteorological data on a daily time scale. We excluded sites without V_cmax data. A total of forty sites in EuroFLUX and one site in AmeriFLUX were selected from FLUXNET 2015 [60]. Time series LAI data were extracted from the GLASS LAI product for the pixels in which the flux towers were located. Because of the potential for uncertainties in both GPP and LAI, we excluded flux data with inconsistent temporal profiles of EC-GPP and LAI. This screening can reduce the amount of noise from the spatial mismatch between the remotely sensed data and field observations. In addition, nighttime flux data were removed. A total of 188 site years were selected for model calibration and validation purposes. The site-level daily EC-GPP, shortwave radiation, temperature, VPD, and LAI were used to drive the FGM.

2.3.2. Forcing Datasets for the FGM

Model inputs related to vegetation and environmental forcing data are listed in

Table 3. Vegetation and environmental inputs for the FGM from 2000 to 2016.

Parameter	Source	Time	Temporal Resolution	Spatial Resolution	Reference
Land use and land cover (LULC)	MODIS C6	2001 to 2016	yearly	500 m	[107]
Leaf area index (LAI)	GLASS V5	2001 to 2016	8-day	500 m	[108,109]
Clumping index (Ω)	MODIS BRDF-derived	2006	8-day	500 m	[102]
Photosynthetic capacity (Vcmax)	Chlorophyll content	2001 to 2016	8-day	500 m	[59,103,104,105,106]
Downward shortwave radiation (DSR)	GLASS V5	2001 to 2016	daily	5 km	[110,111]
Air temperature (Ta)	CRUNCEP	2001 to 2016	6 h	0.5°	[112]
Vapor pressure deficit (VPD)	CRUNCEP	2001 to 2016	6 h	0.5°	[112]
Ambient CO2 concentration	MLO	2001 to 2016	daily	site	http://www.esrl.noaa.gov

. PFTs on a large scale were determined using MODIS land use and land cover data according to the MODIS IGBP classification protocol. The FGM uses meteorological data (downward shortwave radiation (DSR), mean air temperature, and vapor pressure deficit (VPD)) from the Climatic Research Unit-NCEP (CRUNCEP), LAI data from the Global Land Surface Satellite (GLASS) product, and ambient CO₂ concentration from the Mauna Loa Observatory (MLO) as inputs. We further included the spatially resolved model inputs for the clumping index (CI) and V_cmax. The CI data were from the global CI map derived from the MODIS bidirectional reflectance distribution function (BRDF) product [102]. Temporally and spatially continuous global V_cmax maps were estimated based on the remotely sensed chlorophyll content [103,104] using Rubisco–chlorophyll relationships between vegetation types via meta-analyses [105,106]. All data were processed to a target spatial resolution of 500 × 500 m and temporal resolution of 8 days.

2.3.3. Global GPP Products for Intercomparison

To examine the spatial pattern and interannual dynamics of the FGM GPP, a total of nine popular global GPP products (

Table 4. Information on nine GPP products for intercomparison.

GPP	Spatial Resolution	Temporal Resolution	Method	Time Period	Reference
GOSIF	$0.05°$	Annual	Empirical model	2001–2016	[113,120]
BEPS	0.073°	Daily	TBM	2001–2016	[34,37,121,122]
GLASS (v6)	500 m	Annual	LUE model	2001–2016	[28,117]
MODIS (c6)	500 m	Annual	LUE model	2001–2016	[115]
VPM	500 m	Annual	LUE model	2001–2016	[116]
CCW	$0.05°$	Annual	LUE model	2001–2016	[30]
FLUXCOM	$0.5°$	Annual	Machine Learning	2001–2016	[118,119]

) were estimated using different methods, including one empirical model, four light use efficiency (LUE) models, three machining learning methods, and one process-based biophysical model. The GOSIF GPP product was derived from the empirical relationship between GPP and SIF [113,114]. The LUE-based GPP products included the CCW [30], MOD17 [115], VPM [116], and GLASS [28,117] products. Three machine-learning-based GPP products, including an artificial neural network (ANN), the multivariate adaptive regression splines method (MARS), and the random forest method (RF), were derived from FLUXCOM GPP [118,119]. In addition, a global process-based GPP estimation derived from the BEPS was also included [34]. We evaluated the FGM GPP with the collected GPP products and GOSIF data on a yearly scale. In addition, we also evaluated the FGM GPP with the GOSIF data, derived from the SIF soundings of the Orbiting Carbon Observatory-2 (OCO-2), MODIS data, and meteorological reanalysis data [35,113].

3. Results

3.1. Model Evaluation

3.1.1. Including Dynamic V_cmax Information Improved GPP Estimation at EuroFLUX Sites

We first compared the performance of the FGM between the 19 validation sites (Table 1) using the LCC-derived dynamic V_cmax (“All”) and TRY-based constant V_cmax (“PC3”). While “PC3” used only the LAI to describe changes in vegetation status, “All” considered the variability in both the LAI and the V_cmax in the estimation of GPP. Compared with the TRY-based V_cmax, the FGM improved the estimation of daily GPP, with the R² increased from 0.52 to 0.64 and RMSE decreased by 25%. When forcing the intercept to zero, the R² was much higher (0.95 to 0.96). Meanwhile, the scatters between the estimated GPP and EC-GPP from “All” were closer to the 1:1 line, with smaller biases than “PC3” after including the LCC-derived V_cmax (Figure 2).

The FGM model generally showed a reasonable performance in predicting daily GPP compared with in situ GPP measurements across different biomes in Europe (Figure 3 and Figure 4). Overall, the FGM captured an average of 87.5% of the variation in the EC-GPP from the validation dataset. We randomly split the forty EuroFLUX sites into 1:1 by biome to conduct independent model calibration and validation. Our independent validation indicated that the R² between the modeled GPP and EC-GPP ranged from 0.79 to 0.93, with RMSE values ranging from 0.9 to 2.73 g C m⁻² d⁻¹ (Figure 3). The independent validation analysis indicated the strongest model performances for MF and DBF (Figure 3b,f) and less satisfactory agreement for CRO (Figure 3a). The agricultural sites had a relatively high RMSE of 2.72 g C m⁻² d⁻¹ and a relatively low R² of 0.79, mainly based on deviations in amplitude and growing season periods (Figure S1). For example, in the case of DE_Seh, the FGM predicted lower GPP values than those of EC-GPP (Figure S1). The modeled GPP values were close to the values of EC-GPP in 2008 but showed lower productivity in the vegetation phase compared to the values of EC-GPP in 2007 for the same site (e.g., DE_Seh), probably due to uncertainties in the remote sensing products, such as the LAI (Figure S2). Despite some discrepancies, the FGM generally simulated the variations in the GPP on the daily time scale effectively.

3.1.2. FGM GPP Estimations Matched with GOSIF and Other GPP Products

In addition to the flux tower measurements, we further introduced nine global GPP products and remotely sensed GOSIF data to examine the large-scale pattern and temporal dynamics of European GPP estimated by the FGM.

The FGM effectively simulated the general pattern of GPP along the temperature gradient across Europe (Figure 5a). The mean GPP increased from the boreal to temperate regions and decreased from the temperate regions to the Mediterranean regions. The spatial pattern of the annual mean GPP modeled by the FGM was well correlated with that of the other GPP products (r = 0.61–0.8) (Figure 5b–f) and GOSIF data (r = 0.77) (Figure 5g), although there were regional discrepancies in magnitude.

The multiyear mean of the annual GPP of Europe estimated by the FGM (9.4 Pg C yr⁻¹) was reasonable compared with the other GPP products (5.9 to 9.2 Pg C yr⁻¹) (Figure 6a). Moreover, the annual total GPP showed a significant increasing trend from 2001 to 2016 (+0.051 Pg C yr⁻², R² = 0.76, p < 0.01), which was in accordance with the other GPP products (Figure 6a). The annual total GPP across Europe increased from 9.09 Pg C yr⁻¹ in 2001 to 9.94 Pg C yr⁻¹ in 2016. In addition, the interannual dynamics of the FGM GPP correlated well with those of the GOSIF GPP, GLASS GPP, VPM GPP, BEPS GPP, and CCW GPP (Figure 6b). These evaluation results indicated that the FGM GPP was reasonable and could be used to further quantify the V_cmax change effect on GPP.

3.2. Impacts of V_cmax Change on GPP across Europe

3.2.1. Dynamic V_cmax Information Is Important for the Accurate Estimation of GPP Seasonality

We further evaluated the seasonal variation in V_cmax using three types of V_cmax data: (a) V_cmax derived from the LCC (i.e., LCC-derived V_cmax); (b) V_cmax retrieved by the model calibration method (i.e., Calibrated V_cmax); and (c) V_cmax based on the TRY database (i.e., TRY-based V_cmax) (Figure 7). In comparison with the calibrated V_cmax, the LCC-derived V_cmax showed major differences, without a consistent bias in any one direction. In spring and autumn, the LCC-derived V_cmax reduced the overestimation of V_cmax for most PFTs and reduced the underestimation for EBF; in summer, the LCC-derived V_cmax reduced the overestimation of V_cmax for WET, GRA, DBF, and DNF and the underestimation for CRO, SH, MF, EBF, and ENF. In comparison with the V_cmax from the TRY database, we found that the TRY-based V_cmax was consistently larger than the LCC-derived V_cmax throughout the seasons for all biomes, especially for CRO and GRA. The seasonal pattern of the LCC-derived V_cmax was very similar to those of the LAI and GPP.

Consequently, these differences in V_cmax (Figure 8a,c) led to synchronous changes in the GPP estimated by the FGM (Figure 8b,d). In comparison with the GPP from “PC2” (Figure 8b), in spring and autumn, the inclusion of the LCC-derived V_cmax reduced the overestimation of GPP for all the PFTs, especially for MF, DBF, DNF, and ENF. In summer, the LCC-derived V_cmax reduced the overestimation of GPP for GRA and the underestimation of GPP for CRO, SH, MF, and EBF. In comparison with the GPP from “PC3”, the inclusion of the LCC-derived V_cmax reduced the overestimation of GPP using the TRY-based V_cmax (Figure 8d).

3.2.2. Including Dynamic V_cmax Information Improved the Estimation of the GPP Spatial Pattern

The conventional method of fixing V_cmax using the model calibration method and TRY database overestimated the European GPP by 0.5 Pg C yr⁻¹ (Figure 9a) and 2.9 Pg C yr⁻¹ (Figure 9b), respectively. Compared with “PC2”, using spatiotemporally explicit V_cmax information, the terrestrial ecosystem productivity mainly increased for regions in marine climate zones between 43°N–60°N and 0°E–30°E but decreased for most of the regions in other climate zones (Figure 9c). In contrast, the FGM driven by the TRY-based V_cmax overestimated the GPP for almost all the regions, especially for cropland (Figure 9d). Including dynamic V_cmax information improved the FGM’s performance in simulating the variability in European GPP based on the spatial correlations between the estimated annual GPP and other GPP products or GOSIF data (Figure 9e).

3.2.3. Interannual Changes in V_cmax Only Have a Minor Effect on GPP in a Limited Period of 16 Years

We further evaluated the V_cmax change effect on GPP on annual time scales (Figure 10). Compared with “PC2”, the annual GPP of cropland, grassland, and forests demonstrated reductions of 8%, 21%, and 9%, respectively, with the inclusion of the LCC-derived V_cmax. However, we also noticed a slight 8% increase in the GPP for shrubland. Using the TRY-based V_cmax, the FGM overestimated the annual GPP of cropland, grassland, shrubland, and forests by 49%, 34%, 26%, and 15%, respectively.

European vegetation appeared markedly more productive from 2001 to 2016. However, if we fixed V_cmax in the FGM using the calibrated V_cmax or TRY-based V_cmax, both the magnitude of and the rate of increase in GPP were overestimated (Figure 10c). Compared with “All”, “PC2” and “PC3” overestimated the annual increasing GPP trend by 14% and 20%, respectively. The magnitude and interannual dynamics of GPP simulated by the FGM for “PC1” (9.42 Pg C yr⁻¹) were close to that for “All” (9.40 Pg C yr⁻¹). In contrast, the annual GPPs for “PC2” (9.90 Pg C yr⁻¹) and “PC3” (12.30 Pg C yr⁻¹) were significantly higher than that for “All”. Thus, the inclusion of spatial and seasonal variations in V_cmax improved the GPP estimation, while interannual changes in V_cmax contributed little to the GPP in the limited study period of sixteen years.

4. Discussion

4.1. Effects of V_cmax Change on GPP Estimation

Many studies have demonstrated that V_cmax changes across both space and time. Leaf chlorophyll abundance is closely linked to photosynthesis capacity [56,59,123]. The V_cmax derived from the LCC (i.e., LCC-derived V_cmax) showed strong seasonality and significant spatial variation across Europe (Figure S3a) but only slight interannual variation over the limited study period of 16 years (Figure S3b). A comparison with the EC-GPP and other GPP or SIF products supported our hypothesis: the consideration of spatiotemporal changes in V_cmax provided more reliable GPP estimations for Europe. Compared with the GPP estimations obtained by fixing the V_cmax using the TRY database, the inclusion of temporally and spatially explicit V_cmax using the satellite-derived LCC product reduced the bias in the estimated daily GPP (Figure 2) and increased the spatial consistency between the FGM GPP and other GPP products or GOSIF data (Figure 9).

The positive impact of the LCC on GPP simulations on the site level in Europe that found in our study is comparable with the results of previous studies on single sites (with R² enhanced by 10–12% and RMSE decreased by 24–32%) [124,125] or across multiple sites with different PFTs (with R² enhanced by 9–22% and RMSE decreased by 15–32%) [123]. In this study, we found a 23% increase in R² and a 25% decrease in RMSE (Figure 2) between the modeled GPP and EC-GPP for 19 EuroFLUX sites across 9 PFTs using the LCC-derived V_cmax to estimate GPP. On the regional scale, we found a 17% decrease in the annual GPP across the 19 EuroFLUX sites and a 24% decrease in the annual GPP for Europe. The maximum leaf photosynthesis capacity is known to change with the seasons under the influences of multiple factors, such as leaf development [53], changes in climatic variables [126], and drought conditions [127]. In this study, we found a higher V_cmax during summer than in spring and autumn (Figure 7). However, the V_cmax field measurements were generally collected close to the peak growing seasons. Thus, assuming a constant V_cmax based on the TRY database led to the overestimation of V_cmax in spring and autumn, which further resulted in an overestimation of GPP. The overestimation of the TRY-based GPP was in agreement with previous findings observed on a global scale. However, Luo et al. (2019) found only a 7% decrease in the annual GPP across 124 sites and a 7% decrease in the global GPP [123], which is much lower than the regional V_cmax change effect on GPP across Europe (17–24%).

While the V_cmax across Europe showed a small (but not significant) increasing trend from 2001 to 2016 (Figure S5), we found that including the interannual changes in V_cmax had only a minor impact on the interannual GPP change trend for the limited study period of 16 years (Figure 10). However, we could not neglect the interannual changes in V_cmax, since plants may continue to acclimate their leaf chemistry and photosynthesis capacity in response to climate change, especially in response to continued global warming and elevated CO₂ concentrations [128]. According to optimality theory, rising CO₂ and warming can reduce the global canopy demand for Rubisco and result in reductions in V_cmax in the long term [129]. In contrast, we found that the LCC-derived V_cmax showed an increasing trend (that was not statistically significant) across Europe (Figure S3), which is also revealed by the V_cmax estimated from SIF data [51,57]. Plants in high arctic regions are sensitive to changes in temperature [130]. During the study period, the interannual mean air temperature in Europe showed a significant increasing trend (+0.029 °C yr⁻¹, R² = 0.55, p = 0.01) (Figure S4b). However, we also observed enhanced VPD and water stress caused by global warming in Europe (Figure S4c). Plants may adapt to combined changes in different environmental factors, such as radiation brightening, warming temperatures, and enhanced VPD (Figure S4), by increasing their V_cmax to match the light-limited rate of photosynthesis and optimize carbon fixation [49,131,132].

4.2. Comparison with other GPP Products

The increasing interannual trend in the GPP predicted by the FGM (+0.55% yr⁻¹) was in the range of that estimated using other GPP products (+0.47% yr⁻¹ to +0.92% yr⁻¹) (Figure 6). From 2001 to 2016, terrestrial ecosystem productivity showed a significant increasing trend (p < 0.01) in Europe according to five previous GPP products (i.e., BEPS, MODIS, GLASS, GOSIF, and VPM) (Figure 6). The annual total GPP across Europe from 2001 to 2016 predicted by the VPM [116] showed an increasing trend of +0.92% yr⁻¹, which is almost double the predictions of the FGM (Figure 6). Other studies reported that the increasing GPP trend detected by the VPM may be an overestimate [2], since the VPM is not strictly calibrated using field observations at FLUXNET sites [116]. In contrast, the CCW GPP products failed to detect the increasing GPP trend across Europe while successfully capturing the increasing trend of GPP on the global scale [2]. In the case of LUE models, the model parameters, especially those related to the fraction of photosynthetically active radiation (FPAR) and LUE, may have uncertainties and lead to errors in model estimations. Previous studies that estimated GPP dynamics were based mainly on LUE models and process-based models and, in most cases, did not include the spatial and temporal dynamics of V_cmax [51]. These models are unlikely to produce reliable simulations of photosynthesis–climate interactions at a fine temporal resolution or on a large scale.

Prior studies highlighted continuous increases in global terrestrial production during the last two to three decades based on remote sensing data [2,133]. In particular, enhanced GPP mainly occurs in the boreal and temperate regions, where widespread greening and climate warming occur [130]. In this study, we also found that European vegetation showed a significant ‘greening’ trend (i.e., increases in the LAI) from 2001 to 2016 (+0.62% yr⁻¹). In addition, we further found that the GPP increasing trend throughout the European terrestrial ecosystem estimated by the FGM (+0.55% yr⁻¹) was proportional to the greening rate detected by the LAI and other GPP products (0.47–0.67% yr⁻¹) (Figure 11). With the help of spatially and temporally continuous V_cmax maps, we can expect that process-based models will help us to better understand the driving forces of enhanced carbon assimilation in Europe. However, determining how land surface greening, climate change, and other factors contribute to the increase in GPP observed across Europe is beyond the scope of this study and warrants further investigation.

4.3. Uncertainties in V_cmax Data and Implications for Photosynthesis Simulations

To examine the accuracy of the LCC-based V_cmax products used in this study, we first built an observational dataset of V_cmax by compiling field measurements collected at nine sites covering four PFTs (i.e., DBF, EBF, ENF, and GRA) across Europe (Figure 12) [127,134,135,136,137,138,139,140,141]. Then, we compared the mean V_cmax seasonality derived from the LCC, PFT-specific V_cmax, and field measurements of V_cmax for these sites (Figure 12a–i). The distribution of the monthly mean value of V_cmax measurements across the different sites was comparable to that of the corresponding V_cmax derived from the LCC (Figure 12k). The monthly averaged V_cmax values for site NOIT0-03, derived from the LCC, were well correlated with the field data collected (r = 0.64) during the growing season (i.e., from April to October) (Figure 12l). When we calibrated the V_cmax according to the PFTs, the FGM overestimated the V_cmax during spring, autumn, and winter at most sites, the same result as that which we obtained for the whole study area. It is worth noting that the V_cmax used in this study represents the maximum carboxylation rate at a standardized temperature of 25 °C (i.e., V_cmax25), which is a proxy for the Rubisco content rather than a realized rate at ambient temperatures. Thus, here, we observed larger V_cmax values (Figure S3a) than the V_cmax rates at the average growing season temperature reported in other studies [142].

Recently, several global-scale V_cmax products of different spatial and temporal resolutions have been distributed [51,57,58] with contrasting patterns. Large uncertainties still exist regarding the current V_cmax products, especially in terms of the temporal variations during the growing season [53]. Because a limited quantity of field-measured V_cmax data are available for validation, remote sensing V_cmax products have not yet been fully tested. Although efforts have been made to predict V_cmax on the global scale using remote sensing data [55,57,58,131,143], the mechanisms driving the spatiotemporal variability in plant photosynthetic production (e.g., environmental acclimation, leaf age effect) are still ongoing [53,54,139].

We found that the seasonal pattern of the LCC-derived V_cmax was very similar to that of the LAI. Thus, an alternative strategy for the LCC-based V_cmax is to use a PFT-specific V_cmax and scale it according to the LAI seasonality. This LAI-based V_cmax scaling approach was used for some models (e.g., BESS) [33,144] and could result in a robust performance, at least regarding the seasonality aspect. The V_cmax change effect on the GPP investigated here highlights the need for detailed studies using multi-source V_cmax datasets on different scales (e.g., site-level field measurements and large-scale remote sensing retrievals).

Although we improved European GPP estimations by including the spatiotemporal dynamics of V_cmax, there are still some uncertainties regarding the GPP estimated by the FGM. Uncertainties regarding the input parameter datasets are one possible source. Another possible limitation is that the algorithm of the FGM, as presented here, does not use precipitation or soil moisture data directly and implements only VPD, which is partly related to soil moisture. We assumed that VPD would be able to replace soil moisture for the assessment of the influence of drought on GPP. This is an example of how a future model could be improved.

5. Conclusions

In this study, by including the spatial and temporal variations in the maximum photosynthetic capacity rate (i.e., V_cmax) derived from the leaf chlorophyll metric using a remote-sensing-driven process-based model (i.e., FGM), we improved the estimation of the European GPP dynamics from 2001 to 2016 at 8-day time intervals and a 500 m spatial resolution. Compared with the traditional method of fixing the V_cmax as a PFT-specific constant using the empirical parameterization method, we obtained an improved model performance by modeling GPP considering the spatial and temporal variations in V_cmax. The FGM predictions revealed a greening and more productive Europe, consistent with the existing global-scale GPP products and recent literature reports of enhanced carbon sinks in boreal and temperate regions. Our reanalysis suggests that a process-based GPP model using Farquhar’s photosynthesis model requires the careful parameterization of V_cmax to accurately represent the photosynthetic capacity of terrestrial ecosystems. This study contributes to a better understanding of the role of European vegetation in the global carbon cycle.