Remotely Monitoring Vegetation Productivity in Two Contrasting Subtropical Forest Ecosystems Using Solar-Induced Chlorophyll Fluorescence

: Subtropical forests can sequester a larger amount of atmospheric carbon dioxide (CO 2 ) relative to other terrestrial ecosystems through photosynthetic activity and act as an important role in mitigating global climate warming. Compared with the model-based gross primary production (GPP) products, satellite-derived solar-induced ﬂuorescence (SIF) opens a new window for quantiﬁcation. Here, we used the remotely sensed SIF retrievals, two satellite-driven GPP products including MODIS (GPP MOD ) and BESS (GPP BESS ), and tower-based GPP measurements at two contrasting subtropical forests to provide a systematic analysis. Our results revealed that GPP and the associated environmental factors exhibited distinct seasonal patterns. However, the peak GPP values had large differences, with stronger GPP in the evergreen needleleaf forest site (8.76 ± 0.71 g C m − 2 d − 1 ) than that in the evergreen broadleaf forest site (5.71 ± 0.31 g C m − 2 d − 1 ). The satellite-derived SIF retrievals showed great potential in quantifying the variability in GPP, especially for the evergreen needleleaf forest with r reaching up to 0.909 ( p < 0.01). GPP MOD and GPP BESS showed distinctly different performances for the two subtropical forests, whereas the GPP estimates by exclusive use of satellite-based SIF data promised well to the tower-based GPP observations. Multi-year evaluation again conﬁrmed the good performance of the SIF-based GPP estimates. These ﬁndings will provide an alternative framework for quantifying the magnitude of forest GPP and advance our understanding of the carbon sequestration capacity of subtropical forest ecosystems. sum, the present study exhibited the strong linear GPP–SIF relationships across the two subtropical forest ecosystems and great potential to track the variability in tower-based GPP by solely using the satellite-derived SIF product. Compared with two widely-used


Introduction
Earth's climate has changed dramatically in the past century and will keep changing during the next few centuries. The IPPC AR6 report revealed that the global annual nearsurface temperature has been rising steadily from the latter part of the 19th century and places the year 2020 as one of the three warmest years on record [1]. Such changes have exerted substantial effects on terrestrial ecosystem processes and functions [2,3], whereas human-induced activities, especially greenhouse gases emitted from fossil fuel combustion, forest; (2) to develop an alternative method to predict the GPP dynamics by solely using the satellite-derived SIF retrievals via leave-one-out validation; and (3) to evaluate the model performance against satellite-based GPP products. All analyses will be helpful to enhance our understanding of the carbon cycles in the subtropical conifer and broadleaved forests.

Description of the Flux Sites
As illustrated in Figure 1, the study used two forest sites with long-term continuous flux measurements to represent the typical evergreen needleleaf forest and evergreen broadleaf forest in the subtropical humid region of Southeast China, which also belonged to Chinese Terrestrial Ecosystem Flux Research Network (ChinaFLUX). Both sites started to continuously measure the long-term CO2 and H2O fluxes from autumn 2002 via the ECbased method. The Qianyanzhou (QYZ) Ecological Research station (115°03′29.2″ E, 26°44′29.1″ N) is located in Jiangxi Province on the typical red-soil hilly area in the mid-subtropical monsoon climate region with planted coniferous forests. The terrain of the underlying surface is gentle, with the slope gradient ranging from 2.8° to 13.5°. Multi-year mean temperature, precipitation, and relative humility at this site are 17.8 °C, 1471.2 mm, and 84%, respectively [33,34]. The soil type is comprised of red soil, fluvo-aquic soil, paddy soil, and meadow soil, and the parent materials generally include rotten rock, mudstone, glutinite, and river alluvium. Coniferous forests planted inabout 1985 are the dominant vegetation type. The major tree species consist of masson pine, splash pine, and cedarwood, with some Schima superba andorange. The canopy is nearly close with evergreen vegetation occupying approximately 76% of total land area.
The Dinghushan (DHS) flux site (112°32′3.8″ E, 23°10′24″ N) is situated in the evergreen broadleaf forest, in the kernel area of the Dinghushan Biosphere Reserve, Guangdong Province, within the southern subtropical region. This site has abundant resources of radiation, heat, and rainfall. Multi-year mean temperature, precipitation, and relative The Qianyanzhou (QYZ) Ecological Research station (115 • 03 29.2" E, 26 • 44 29.1" N) is located in Jiangxi Province on the typical red-soil hilly area in the mid-subtropical monsoon climate region with planted coniferous forests. The terrain of the underlying surface is gentle, with the slope gradient ranging from 2.8 • to 13.5 • . Multi-year mean temperature, precipitation, and relative humility at this site are 17.8 • C, 1471.2 mm, and 84%, respectively [33,34]. The soil type is comprised of red soil, fluvo-aquic soil, paddy soil, and meadow soil, and the parent materials generally include rotten rock, mudstone, glutinite, and river alluvium. Coniferous forests planted inabout 1985 are the dominant vegetation type. The major tree species consist of masson pine, splash pine, and cedarwood, with some Schima superba andorange. The canopy is nearly close with evergreen vegetation occupying approximately 76% of total land area.
The Dinghushan (DHS) flux site (112 • 32 3.8" E, 23 • 10 24" N) is situated in the evergreen broadleaf forest, in the kernel area of the Dinghushan Biosphere Reserve, Guangdong Province, within the southern subtropical region. This site has abundant resources of radiation, heat, and rainfall. Multi-year mean temperature, precipitation, and relative humidity are 20.9 • C, 1956 mm, and 82%, respectively [35,36]. The terrain is basically flat, especially in the prevailing northeast wind direction. The soil type is mostly lateritic, and a humus layer is covered. Monsoon evergreen broadleaf forest is the typical vegetation, with the dominant species of Castanopis chinensis, Schim asuperba, Cryptocarya chnensis. The stand is about 400 years old, with complicated forest structure including six layers, namely, three sub-layers of arbor, sapling shrub layer, herbaceous seedling layer, and interlayer plant layer. The mean canopy height is about 17 m, and the leaf area index (LAI) spans from 3.8 to 4.2 m −2 m −2 [37].

Data Processing of the EC-Based Fluxes
As an indispensable part of the global-scale long-term network of micrometeorological flux observations, ChinaFLUX addressed to reveal the variability in land-atmosphere carbon and water fluxes, as well as the underlying mechanisms controlling it [38,39]. In this study, the EC-based flux measurements between 2003 and 2010 at the two subtropical forest flux sites (QYZ and DHS) were used for analysis. At these flux sites, the eddy covariance system together with an automatic meteorological station (AWS) were mounted and continuously observe the site-scale net carbon exchanges (NEE) and the related meteorological data. These synchronously observed climate data consisted of near-surface air temperature (T a ), solar radiation (R g ), soil temperature and water content (T s and SWC), relative humidity, vapor pressure deficit (VPD) and the amount of precipitation (P). The open-path EC system contains two sensors including a fast response infrared CO 2 /H 2 O gas analyser (Li-7500, LI-COR Biosciences, NE, USA) and a 3-D sonic anemometer (CSAT3, Campbell Scientific, UT, USA). The raw data are sampled at the frequency of 10 Hz by the data logger (CR5000, Campbell Scientific, UT, USA). Then, the procedures including spike detection and despiking, two-dimensional coordinate rotation, time delay removal of H 2 O and CO 2 , virtual temperature correction, density effects (WPL correction) and frequency response corrections were completed to produce a half-hour flux dataset [40,41]. More details about the processing procedures can be seen in Tang et al. [42].
Furthermore, owing to the power failure, instrument malfunctions, and severe weather conditions, a few flux observations were lost. These gaps were interpolated with a standardized algorithm developed by ChinaFlux [40,41]. Meanwhile, these site-level NEE measurements were separated into GPP and R e using a short-term temperature dependent method [15,43]. The R-based software (REddyProc, MPI, Jena, Germany) developed by the Max Planck Institute (MPI) for Biogeochemistry was applied to complete the gap-filling and flux partitioning procedures [44]. Finally, the half-hourly GPP data were processed into the daily time scale and averaged over each 8-day period to match the composites of satellite-based SIF and GPP products.

SIF Data Set
The space-borne SIF observation offers the possibility to accurately quantify terrestrial ecosystem productivity. Recent research by Sun et al. [26] and Li et al. 2018 [45] proved that the satellite-derived SIF retrievals from the NASA's Orbiting Carbon Observatory-2 (OCO-2) have become an alternative proxy for vegetation GPP across multiple spatial and temporal scales. Nevertheless, the tremendous applicable values were severely constrained by the spatially-sparse nature of OCO-2 data, particularly for regional and global scales. Thus, this study used a newly released global contiguous SIF product with the spatial resolution of 0.05 degree and temporal resolution of 8 days (GOSIF V2, http://globalecology.unh.edu, accessed on 20 October 2021). These SIF retrievals were derived from OCO-2, MERRA-2, and MODIS data sets and agreed well with the site-level GPP observations from more than 90 flux sites spanning various vegetation types. The GOSIF product exhibits amounts of advantages such as better spatial resolution, longer observation periods, and spatially continuous coverage in comparison to the coarse-resolution SIF product directly from discrete SIF signals. More descriptions of the features on the GOSIF data can be found in Li & Xiao [28]. Here, we used the GOSIF time series between 2003 and 2010 to extract the SIF retrievals at the two subtropical forest flux sites and explore its potential to capture the variability in site-level GPP.

MODIS GPP Data
The study used the satellite-driven MOD17A2 GPP data (V6) based on the concept of a light use efficiency model for comparison, which has been widely used to repeatedly monitor global vegetation conditions, as well as the associated carbon sequestration potentials [46], and as inputs to terrestrial land surface models to infer the energy and ecohydrological processes [17]. The satellite remote-sensing-based vegetation-related dynamic input data of the model used the MOD15A2H data set, including the composite of photosynthetically active radiation (FPAR) and leaf area index (LAI) for 8-day periods. The missing data owing to cloudiness in this product are gap-filled according to the MODIS quality assessment (QA) data layers. Meanwhile, the model used the meteorological data from 6-h National Center for Environmental Prediction/Department of Energy (NCEP/DOE) reanalysis II products to drive the algorithm. The relevant Biome Property Look-up Tables (BPLUT) were also updated across the plant functional types. Currently, these GPP products have several time scales including 8-day, monthly, and annually, and this study used the 8-day mean values, which can be accessed through the website https://lpdaac.usgs.gov/products/mod17a2hv006/ (accessed on 20 October 2021).

BESS GPP Product
The Breathing Earth System Simulator (BESS) algorithm used seven MODIS-driven atmosphere and land data, four different satellite products, four reanalysis data, and three ancillary input data sets to produce global-scale GPP products. As an ecological process-based model, the modules including atmosphere and canopy radiative transfer algorithm, canopy photosynthesis algorithm, transpiration algorithm, and energy balance algorithm were coupled together [47]. The coupled carbon-water module composed of Farquhar's algorithm for biochemical photosynthetic pathways such as C 3 and C 4 plants, a two-leaf canopy radiative transfer algorithm, and the quadratic Penman-Monteith model is integrated to quantify terrestrial GPP for sunlit and shade canopy by iterative steps. Jiang & Ryu [48] improved the model algorithm and evaluated BESS GPP against FLUXNET 2015 dataset with MPI-BGC product across multiple spatio-temporal scales. The performance has also been evaluated by previous research [21,49]. Now, the global BESS GPP dataset between 2000 and 2015 with 1 km spatial resolution for 8-day periods can be accessed at the Ecological Sensing Artificial Intelligence (AI) Laboratory at the Soul National University (the website: http://environment.snu.ac.kr/, accessed on 22 October 2021).

Statistical Analysis
In total, 16 site-years of flux measurements at the two contrasting subtropical forest sites were used to explore the main environmental factors determining the dynamics in GPP of the evergreen needleleaf forest and evergreen broadleaf forest. The Pearson correlation coefficients (r) jointly with relevant p-value between these in situ GPP measurements, and the vegetation and climate-related factors including SIF, R g , T a , VPD, and P were calculated. When the r value is close to 1, a strong positive linear relationship is indicated when it passes the significance test (p < 0.01).
Then, this study aimed to develop an alternative method to retrieve GPP by exclusive use of satellite-based SIF data and evaluated its performance against satellite-based MODIS and BESS GPP products. Here, we used the leave-one-out cross-validation approach to test the reliability and robustness of the developed model. For each forest type, there were a total of 8 years of GPP measurements. One year of the GPP observations were withheld successively as a testing dataset, and then the data from the remaining seven years were used as training samples for model development. Then, we compared the model-estimated GPP including the SIF-based GPP (GPP SIF ), MODIS GPP (GPP MOD ), and BESS GPP (GPP BESS ) with the tower-measured GPP, respectively.
Two widely used indicators including the coefficient of determination (R 2 ) and the root mean square error (RMSE) were used to assess the model's performance. The equations of R 2 and RMSE are as follows: where y and y represented the observed and predicted values of GPP; n is the total number of observations (eight-day periods). Generally, the optimal GPP estimates had the largest R 2 and least RMSE. In this study, all statistical analyses are implemented by SPSS 19.0 (IBM, Chicago, IL, USA). When the tower-measured and model-estimated GPP regression was close to the 1:1 line, it indicated a high precision of the model estimation. However, the scatter plots between site-level GPP and modeled GPP can only reflect the overall consistency. Thus, in order to examine the performance of the modeled GPP over time, this study also compared the seasonal consistency between tower-based GPP observations (GPP EC ) at the two subtropical forest flux sites and three different satellitebased GPP products.

Dynamics in GPP as Well as Environmental Factors
Seasonal variations in site-level GPP and climate-related factors including R g , T a , VPD, and P across the two contrasting subtropical conifer and broadleaved forest ecosystems are illustrated in Figures 2-4. For each forest type, a total of 8 years of flux observations was used for analysis. It implied that the GPP of both the QYZ and DHS sites showed apparent single-peak patterns. As the temperature rose in springtime (Figure 3), GPP enhanced increasingly with the plant growth status, peaked in summertime, and decreased again in autumn. However, the peak timing appeared earlier in the evergreen needleleaf forest at the end of June or early July, by comparison with the evergreen broadleaf forest during about middle to late July. In addition, the peak values of GPP had large differences, with larger photosynthetic carbon fixation potential at the QYZ site (8.76 ± 0.71 g C m −2 d −1 ) than that at the DHS site (5.71 ± 0.31 g C m −2 d −1 ). The distribution of rainfall is mainly concentrated in summer with large fluctuations throughout the measurement period ( Figure 3). As a proxy for plant water stress, the VPD of the QYZ site also showed distinct seasonal dynamics, whereas the VPD of the DHS site fluctuated severely around the year.      Table 1 compared the correlations between these site-level GPP observations and the environmental factors including R g , T a , VPD, and P at the two subtropical forest flux sites, respectively. It revealed that, except P, almost all variables strongly affected the dynamics in forest GPP. In the evergreen needleleaf forest, GPP was proven to exert a stronger influence by T a (r = 0.908), followed by R g and VPD of 0.855 and 0.793, respectively. Contrastingly, R g was the dominant environmental factor of GPP in the evergreen broadleaf forest (r = 0.79). Then, both VPD and T a also strongly and positively affected the variability in GPP with similar correlations of 0.677 and 0.655, respectively.

Model Development by Exclusive Use of SIF Data
The study also examined the relationships between these 8-day tower-based GPPs and the associated SIF products at the two subtropical forest sites. As shown in Figure 5, the GPP at both the QYZ and DHS sites exhibited strong correlations to these satellitederived SIF retrievals, revealing the great potential to monitor the dynamics of GPP by the remote-sensing-based SIF technique. The correlation between GPP and SIF reached up to 0.909 (p < 0.01) in the evergreen needleleaf forest. Although the r value was slightly lower (r = 0.589) in the evergreen broadleaf forest, it also exhibited a strong linear relationship with SIF.

Model Development by Exclusive Use of SIF Data
The study also examined the relationships between these 8-day tower-based GPPs and the associated SIF products at the two subtropical forest sites. As shown in Figure 5, the GPP at both the QYZ and DHS sites exhibited strong correlations to these satellitederived SIF retrievals, revealing the great potential to monitor the dynamics of GPP by the remote-sensing-based SIF technique. The correlation between GPP and SIF reached On this basis, the work proposed the models for quantifying the seasonal dynamics of GPP at the two subtropical forest ecosystems by exclusive use of satellite-based SIF data. Then, this study evaluated the model-estimated GPP performance using the in situ measured GPP through leave-one-out cross-validation method ( Table 2). At the evergreen needleleaf forest, the regressions of site-level GPP observations and GPP estimates in the cross-validation implied R 2 changed between 0.75 and 0.94, and the associated RMSE values varied from 0.71 to 1.11 g C m −2 d −1 for the GPP estimates by withholding years. Meanwhile, R 2 changed between 0. 28  On this basis, the work proposed the models for quantifying the seasonal dynamics of GPP at the two subtropical forest ecosystems by exclusive use of satellite-based SIF data. Then, this study evaluated the model-estimated GPP performance using the in situ measured GPP through leave-one-out cross-validation method ( Table 2). At the evergreen needleleaf forest, the regressions of site-level GPP observations and GPP estimates in the cross-validation implied R 2 changed between 0.75 and 0.94, and the associated RMSE values varied from 0.71 to 1.11 g C m −2 d −1 for the GPP estimates by withholding years. Meanwhile, R 2 changed between 0.28 and 0.61, and the associated RMSE values varied from 0.80 to 1.03 g C m −2 d −1 for the GPP estimates by withholding years at the evergreen broadleaf forest. For all site-years measurements, the regression of in situ measured GPP and the estimated GPP in the cross-validation implied that R 2 = 0.82 and RMSE = 0.93 g C m −2 d −1 at the QYZ site (Figure 6a), and R 2 = 0.35 and RMSE = 0.95 g C m −2 d −1 at the DHS site ( Figure 6b). As illustrated in Figure 6, the scatterplots between GPP observations and GPP estimates were close to the 1:1 regression line, demonstrated by the great performance of such SIF-based method. Table 2. Accuracy assessment of the model-estimated GPP relying on time-series SIF data at the two subtropical forest flux sites through leave-one-out cross-validation method.

Site
Validation Year   Table 2. Accuracy assessment of the model-estimated GPP relying on time-series SIF data at the two subtropical forest flux sites through leave-one-out cross-validation method.

Model Performance against Satellite-Based GPP Products
Furthermore, this study evaluated the model performance (GPPSIF) against satellitebased GPP products including MODIS GPP (GPPMOD) and BESS GPP (GPPBESS) across the two subtropical forests (Figures 6-8). Overall, among these GPP estimates, GPPSIF exhibited the optimal performance in promising the variability in site-level GPP observations. However, GPPMOD and GPPBESS showed distinctly different performances for the evergreen needleleaf forest and evergreen broadleaf forest, respectively. As illustrated in Figures 6  and 7 at the QYZ site, GPPMOD consistently underestimated the tower-measured GPP with large discrepancies (Figure 6a; RMSE = 2.92 g C m −2 d −1 ), whereas GPPBESS agreed well the variability in the measured GPP with RMSE of 1.10 g C m −2 d −1 (Figure 6c). As illustrated in Figures 6 and 8 at the DHS site, GPPBESS severely overestimated the tower-measured GPP in the evergreen broadleaf forest with RMSE of 2.69 g C m −2 d −1 (Figure 6d), whereas

Model Performance against Satellite-Based GPP Products
Furthermore, this study evaluated the model performance (GPP SIF ) against satellitebased GPP products including MODIS GPP (GPP MOD ) and BESS GPP (GPP BESS ) across the two subtropical forests (Figures 6-8). Overall, among these GPP estimates, GPP SIF exhibited the optimal performance in promising the variability in site-level GPP observations. However, GPP MOD and GPP BESS showed distinctly different performances for the evergreen needleleaf forest and evergreen broadleaf forest, respectively. As illustrated in Figures 6 and 7 at the QYZ site, GPP MOD consistently underestimated the tower-measured GPP with large discrepancies (Figure 6a; RMSE = 2.92 g C m −2 d −1 ), whereas GPP BESS agreed well the variability in the measured GPP with RMSE of 1.10 g C m −2 d −1 (Figure 6c). As illustrated in Figures 6 and 8 at the DHS site, GPP BESS severely overestimated the towermeasured GPP in the evergreen broadleaf forest with RMSE of 2.69 g C m −2 d −1 (Figure 6d), whereas GPP MOD provided slightly better estimates (Figure 6b; RMSE = 1.44 g C m −2 d −1 ). Nevertheless, this study also found that GPP SIF seemed unable to capture the dramatic decline in the coniferous site in summer.
GPPMOD provided slightly better estimates (Figure 6b; RMSE = 1.44 g C m −2 d −1 ). Nevertheless, this study also found that GPPSIF seemed unable to capture the dramatic decline in the coniferous site in summer.   GPPMOD provided slightly better estimates (Figure 6b; RMSE = 1.44 g C m −2 d −1 ). Nevertheless, this study also found that GPPSIF seemed unable to capture the dramatic decline in the coniferous site in summer.

Yearly Evaluation of GPP
As illustrated in Figure 9, the study work examined the multi-year mean tower-based GPP observations (GPP EC ) and three different satellite-based GPP products including the GPP SIF , GPP MOD and GPP BESS at the two subtropical forest flux sites. At the evergreen needleleaf forest site, multi-year mean tower-measured GPP was 1799 ± 95 g C m −2 y −1 , which was over twice of the estimate by GPP MOD (812 ± 39 g C m −2 y −1 ), and close to the estimates by GPP SIF (1800 ± 70 g C m −2 y −1 ) and GPP BESS (1795 ± 83 g C m −2 y −1 ). At the evergreen broadleaf forest site, multi-year mean tower-measured GPP was 1411 ± 89 g C m −2 y −1 , which was quite close to the estimate by GPP SIF (1415 ± 26 g C m −2 y −1 ), whereas GPP MOD and GPP BESS exhibited a larger overestimation by about 22.8% and 55.0%, respectively. All analysis again confirmed the good performance by the SIF-based GPP estimates.

Yearly Evaluation of GPP
As illustrated in Figure 9, the study work examined the multi-year mean tower-based GPP observations (GPPEC) and three different satellite-based GPP products including the GPPSIF, GPPMOD and GPPBESS at the two subtropical forest flux sites. At the evergreen needleleaf forest site, multi-year mean tower-measured GPP was 1799 ± 95 g C m −2 y −1 , which was over twice of the estimate by GPPMOD (812 ± 39 g C m −2 y −1 ), and close to the estimates by GPPSIF (1800 ± 70 g C m −2 y −1 ) and GPPBESS (1795 ± 83 g C m −2 y −1 ). At the evergreen broadleaf forest site, multi-year mean tower-measured GPP was 1411 ± 89 g C m −2 y −1 , which was quite close to the estimate by GPPSIF (1415 ± 26 g C m −2 y −1 ), whereas GPPMOD and GPPBESS exhibited a larger overestimation by about 22.8% and 55.0%, respectively. All analysis again confirmed the good performance by the SIF-based GPP estimates.

Discussion
Climate-induced issues must be mitigated by diverse stakeholders such as the policymakers, investors, companies, scientists, and the public [32]. In response to the changing climate, currently, more and more governments are implementing their aims to combat global climate warming [50]. In the fall of 2020, the Chinese government also made a specific target to reach its carbon peak by 2030 and then realizing carbon neutrality by 2060. Carbon neutrality means that one region or country gradually achieves net-zero CO2 emissions by adopting diverse management measures. Terrestrial vegetation GPP represents the largest carbon flux amongst global carbon cycles and supplies important ecosystem services for human society [51,52]. Both the Paris Agreement and the IPCC reports admitted that the targets of climate change mitigation could not be realized without the important role of terrestrial forests, whereas the rate of contribution by forests in reducing atmospheric greenhouse gas concentrations remains a challenge [53,54].
The traditional vegetation monitoring from the satellite-based Earth observation (EO) technique is usually constrained by the vegetation structure and the sensitivity to aboveground biomass density, especially for the complicated forest ecosystems [55][56][57],

Discussion
Climate-induced issues must be mitigated by diverse stakeholders such as the policymakers, investors, companies, scientists, and the public [32]. In response to the changing climate, currently, more and more governments are implementing their aims to combat global climate warming [50]. In the fall of 2020, the Chinese government also made a specific target to reach its carbon peak by 2030 and then realizing carbon neutrality by 2060. Carbon neutrality means that one region or country gradually achieves net-zero CO 2 emissions by adopting diverse management measures. Terrestrial vegetation GPP represents the largest carbon flux amongst global carbon cycles and supplies important ecosystem services for human society [51,52]. Both the Paris Agreement and the IPCC reports admitted that the targets of climate change mitigation could not be realized without the important role of terrestrial forests, whereas the rate of contribution by forests in reducing atmospheric greenhouse gas concentrations remains a challenge [53,54].
The traditional vegetation monitoring from the satellite-based Earth observation (EO) technique is usually constrained by the vegetation structure and the sensitivity to aboveground biomass density, especially for the complicated forest ecosystems [55][56][57], leading to the spatial pattern and long-term changes of forest GPP being poorly understood. Vegetation GPP can be quantified using satellite remote sensing and ecological models, but Tang et al. [42] revealed that the widely used MODIS product was unable to capture the longterm dynamics of global karst ecosystem GPP. In this study, the contrasting performance of the MODIS GPP and BESS GPP products were revealed when monitoring GPP of the evergreen needleleaf forest and evergreen broadleaf forest ( Figure 6). MODIS GPP exhibited better performance in the broadleaved forest (Figure 6b), whereas a large underestimation was found for the coniferous forest (Figure 6a). Wang et al. [20] systematically evaluated the latest MODIS GPP products across multiple biomes using global eddy covariance flux data and revealed that the ε max of MOD17A2H have the largest underestimation reaching up to 31% for the evergreen needleleaf forest. In contrast, BESS GPP showed better performance in the coniferous forest (Figure 6c), which was consistent with the results by Lin et al. [58], with the R 2 and RMRE being 0.76 and 1.66 g C m −2 d −1 , respectively, for the evergreen needleleaf forest, and R 2 and RMRE being 0.35 and 2.45 g C m −2 d −1 , respectively, for the evergreen broadleaf forest.
Fortunately, scientists found that when vegetation chlorophyll is exposed to solar radiation, SIF will be emitted as a by-product of vegetation photosynthesis activity, which is different from the reflectance-based biophysical parameters. The satellite-based suninduced fluorescence technique has been developed in recent years and recognized as a rapidly advancing front on the terrestrial vegetation science [26,59]. SIF showed more advantages in monitoring terrestrial GPP by comparison with other remotely sensed indicators such as vegetation indices. Because SIF is related to both the amount of absorbed photosynthetic active radiation (APAR) as well as the efficiency of APAR used to carry out photosynthesis [60], satellite SIF retrievals thus have great potential to monitor photosynthetic activity, which has also been confirmed by our study (Figures 5 and 6). Nevertheless, the spatial scale mismatch between the actual flux tower footprint and the pixel of GOSIF product used in the comparisons is an issue. The coarse spatiotemporal sampling and low signal-to-noise ratio of the raw SIF data hamper their application over small or fragmented ecosystems. Fortunately, with the development of satellite SIF observations, TROPOMI within the TROPOSIF project funded by the European Space Agency can provide daily global measurements at a much denser spatial and temporal sampling than earlier satellite instruments. Recently, a new global SIF dataset between May 2018 and April 2021 was produced [61], which could reduce the uncertainty of scale mismatch to a certain extent. Although this study developed the accurate SIF-based models to capture the variability in tower-based GPP, more forest flux sites across different climate zones and stand age structures are necessary to enhance the empirical model's robustness. In addition, a recent study [59] revealed that extreme weather events such as heatwave broke down the linearity between SIF and GPP of evergreen broadleaved trees. As illustrated in Figure 7, this work also found that the SIF lacked sensitivity to the influence of summer drought in 2003 on forest GPP in the QYZ site. However, the process-based BESS GPP captured the dramatic decline.
Multi-year mean GPP observations in the subtropical conifer and broadleaved forests in the study reached up to 1799 ± 95 g C m −2 y −1 and 1411 ± 89 g C m −2 y −1 , respectively, which were apparently larger than that of the boreal black spruce (812 g C m −2 y −1 ) and Scots pine forest (about 1000 g C m −2 y −1 ) [62,63]. Recent work has also implied that the implementation of a series of afforestation and reforestation programmes could significantly increase vegetation cover condition and enhance the vegetation carbon sequestration capacity [42,64]. Alternatively, deforestation could severely reduce the forest coverage, as well as the carbon fixation [65]. However, the explicit effect of various human-induced activities on regional or global carbon budget remained unclear. Our study provided a new framework to quantify the magnitude and dynamic of forest carbon fixation.

Conclusions
In sum, the present study exhibited the strong linear GPP-SIF relationships across the two subtropical forest ecosystems and great potential to track the variability in tower-based GPP by solely using the satellite-derived SIF product. Compared with two widely-used