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

Abstract

Subtropical forests can sequester a larger amount of atmospheric carbon dioxide (CO₂) 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 fluorescence (SIF) opens a new window for quantification. 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⁻² d⁻¹) than that in the evergreen broadleaf forest site (5.71 ± 0.31 g C m⁻² d⁻¹). 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 confirmed the good performance of the SIF-based GPP estimates. These findings 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.

1. 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 near-surface 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, contributed a lot to this warming [4,5]. Fortunately, as an important part of terrestrial ecosystems, forests make up 30% of the Earth’s land surface, which plays a crucial role in regulating atmospheric CO₂ concentrations under the influence of both photosynthesis and respiration [6,7,8]. Essentially, the carbon sequestration potential of terrestrial forests is jointly affected by climate change and human-induced perturbation through land-use changes including deforestation and afforestation [9,10,11]. It is therefore imperative to accurately evaluate the dynamics of terrestrial forest carbon cycles in order to predict the feedback on Earth’s climate system.

With the development of the eddy covariance-based flux measurements (EC) in recent decades, the ongoing land–atmosphere carbon and water exchanges can be directly observed [12]. Currently, over 900 EC-based flux towers on the Earth are operated in the long run and jointly form the global network FLUXNET [13]. These sites are located in various ecosystems such as tropical and subtropical forests, temperate conifer and broadleaf (evergreen and deciduous) forests, croplands, wetlands, grasslands, and tundra. Forest flux sites occupy over a half of theFLUXNET2015 dataset. These EC-based carbon flux measurements (NEE) can be further separated into gross primary production (GPP) and ecosystem respiration (R_e) through the daytime-based [14] approach and the nighttime-based [15] partitioning method. Nevertheless, these tower-based fluxes were generally observed at site-level with small areas (<1 km²), and long-term financial support is necessary to construct and maintain these instruments [16,17]. In order to quantify the spatial and temporal patterns of terrestrial GPP across large areas, these spatially sparse flux measurements need to be upscaled to spatially continuous estimates.

Various methods, including the satellite-based vegetation indices or leaf area index and process-based and light-use efficiency GPP models, were developed and compared to evaluate the spatio-temporal dynamics of vegetation productivity [18,19]. Despite that, large uncertainties still exist in these GPP estimates, owing to the model structure and different input data sets such as the fractional photosynthetically active radiation (FPAR) product, the potential error contributions from land cover classification and maximum light use efficiency (ε_max) [20,21]. In recent years, satellite-based sun-induced vegetation chlorophyll fluorescence (SIF) observations, like the space-borne SIF retrievals of GOSAT [22], GOME-2 [23], and OCO-2 [24], have been recognized as a new proxy for vegetation GPP over various terrestrial ecosystems. Because of its inherent connection with photosynthesis, as well as the simple and robust relationship between SIF signals and GPP observations [25,26,27], these SIF products are crucial to examining the dynamics of ecosystems and responses to climatic change. Even so, the tremendous applicable values were severely constrained by such spatially sparse SIF signals from site to regional or global levels. A high-resolution, spatially contiguous global SIF product from space-borne SIF retrievals can fulfill the need for repeatedly and precisely quantifying vegetation productivity over different time and space scales [28,29], which can bridge the gaps between site-level EC-based flux towers and satellite-driven SIF data sets.

Subtropical forest ecosystems in southern China play a critical role in mitigating climate change from regional to even global scales due to the strong carbon fixation capability [30,31], and the main forest types are comprised of evergreen needleleaf forest and evergreen broadleaf forest depending on the contrasting tree species. Harris et al. [32] found that the tropical and subtropical forests sequestered more atmospheric CO₂ emissions relative to the temperate and boreal forests (−8.6 vs. −4.4 and −2.5 Gt CO₂ year⁻¹, respectively). However, considerable uncertainties still exist in these satellite-derived GPP products. The study used two subtropical forest flux sites, jointly with the satellite-derived GPP estimates based on ecosystem process-based and light-use efficiency models as well as a new globally reconstructed SIF datasets, (1) to analyze the seasonal dynamics in site-level GPP observations, and the associated environmental controls across the two subtropical forest flux sites including evergreen needleleaf forest and evergreen broadleaf 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.

2. Materials and Methods

2.1. 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 CO₂ and H₂O fluxes from autumn 2002 via the EC-based 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 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⁻² m⁻² [37].

2.2. 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₂/H₂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₂O and CO₂, 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.

2.3. 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.

2.4. 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).

2.5. 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₃ and C₄ 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).

2.6. 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²) and the root mean square error (RMSE) were used to assess the model’s performance. The equations of R² and RMSE are as follows:

$$R^2=1-\raisebox{1ex}{${\displaystyle \sum {(\stackrel{⌢}{y}-y)}^2}$}\!\left/ \!\raisebox{-1ex}{${\displaystyle \sum y^2}$}\right.$$

(1)

$$RMSE=\sqrt{\raisebox{1ex}{${\displaystyle \sum {(\stackrel{⌢}{y}-y)}^2}$}\!\left/ \!\raisebox{-1ex}{$n$}\right.}$$

(2)

where y and $\stackrel{⌢}{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² 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 satellite-based GPP products.

3. Results

3.1. 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 Figure 2, Figure 3 and Figure 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⁻² d⁻¹) than that at the DHS site (5.71 ± 0.31 g C m⁻² d⁻¹). 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. Correlations between site-level GPP observations and the environmental factors including solar radiation (R_g), air temperature (T_a), vapour pressure deficit (VPD), and precipitation (P) at the two subtropical forest flux sites, respectively.

	Site	Rg (W m−2)	Ta (°C)	VPD (h Pa)	P (mm)
GPP (g C m−2 d−1)	QYZ (Evergreen needleleaf forest)	0.855 **	0.908 **	0.793 **	0.051
	DHS (Evergreen broadleaf forest)	0.790 **	0.665 **	0.677 **	−0.069

Note: ** means that the correlation is significant at the 0.01 level (p < 0.01).

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.

3.2. 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 satellite-derived 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.

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. 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	R2	RMSE (g C m−2 d−1)
GPP (g C m−2 d−1)	QYZ (Evergreen needleleaf forest)	2003	0.78	0.98
		2004	0.94	0.90
		2005	0.90	0.71
		2006	0.86	0.88
		2007	0.82	0.89
		2008	0.87	0.85
		2009	0.83	0.94
		2010	0.75	1.11
	DHS (Evergreen broadleaf forest)	2003	0.40	0.80
		2004	0.61	0.89
		2005	0.42	1.03
		2006	0.30	0.98
		2007	0.39	0.99
		2008	0.45	0.91
		2009	0.31	1.01
		2010	0.28	0.92

). At the evergreen needleleaf forest, the regressions of site-level GPP observations and GPP estimates in the cross-validation implied R² changed between 0.75 and 0.94, and the associated RMSE values varied from 0.71 to 1.11 g C m⁻² d⁻¹ for the GPP estimates by withholding years. Meanwhile, R² changed between 0.28 and 0.61, and the associated RMSE values varied from 0.80 to 1.03 g C m⁻² d⁻¹ 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² = 0.82 and RMSE = 0.93 g C m⁻² d⁻¹ at the QYZ site (Figure 6a), and R² = 0.35 and RMSE = 0.95 g C m⁻² d⁻¹ 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.

3.3. Model Performance against Satellite-Based GPP Products

Furthermore, this study evaluated the model performance (GPP_SIF) against satellite-based GPP products including MODIS GPP (GPP_MOD) and BESS GPP (GPP_BESS) across the two subtropical forests (Figure 6, Figure 7 and Figure 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 Figure 6 and Figure 7 at the QYZ site, GPP_MOD consistently underestimated the tower-measured GPP with large discrepancies (Figure 6a; RMSE = 2.92 g C m⁻² d⁻¹), whereas GPP_BESS agreed well the variability in the measured GPP with RMSE of 1.10 g C m⁻² d⁻¹ (Figure 6c). As illustrated in Figure 6 and Figure 8 at the DHS site, GPP_BESS severely overestimated the tower-measured GPP in the evergreen broadleaf forest with RMSE of 2.69 g C m⁻² d⁻¹ (Figure 6d), whereas GPP_MOD provided slightly better estimates (Figure 6b; RMSE = 1.44 g C m⁻² d⁻¹). Nevertheless, this study also found that GPP_SIF seemed unable to capture the dramatic decline in the coniferous site in summer.

3.4. 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⁻² y⁻¹, which was over twice of the estimate by GPP_MOD (812 ± 39 g C m⁻² y⁻¹), and close to the estimates by GPP_SIF (1800 ± 70 g C m⁻² y⁻¹) and GPP_BESS (1795 ± 83 g C m⁻² y⁻¹). At the evergreen broadleaf forest site, multi-year mean tower-measured GPP was 1411 ± 89 g C m⁻² y⁻¹, which was quite close to the estimate by GPP_SIF (1415 ± 26 g C m⁻² y⁻¹), 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.

4. 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₂ 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 long-term 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² and RMRE being 0.76 and 1.66 g C m⁻² d⁻¹, respectively, for the evergreen needleleaf forest, and R² and RMRE being 0.35 and 2.45 g C m⁻² d⁻¹, 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 sun-induced 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 (Figure 5 and Figure 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⁻² y⁻¹ and 1411 ± 89 g C m⁻² y⁻¹, respectively, which were apparently larger than that of the boreal black spruce (812 g C m⁻² y⁻¹) and Scots pine forest (about 1000 g C m⁻² y⁻¹) [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.

5. 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 MODIS and BESS GPP products, the SIF-based GPP estimates showed better performance at the seasonal and annual time scales. The results can help improve our understanding of the carbon cycles in the subtropical forests and pave a novel way toward monitoring GPP accurately for evaluating climate change mitigation performance.