Long-term dynamics of peak photosynthesis timing and environmental controls in the Tibetan Plateau monitored by satellite solar-induced chlorophyll fluorescence

ABSTRACT

Vegetation phenological metrics have attracted increasing attentions during recent years. However, compared with the peak greenness timing (PGT), the peak photosynthesis timing (PPT) has been overlooked until now. As the world's third pole, the Tibetan Plateau is one of the most sensitive areas to global climatic change. In this study, long-term time-series solar-induced chlorophyll fluorescence (SIF) data were used to identify the PPT patterns, trends, and responses in the Tibetan Plateau over the past two decades. Our findings indicated that the SIF-derived PPT can better capture the characteristics from flux-measured gross primary productivity compared with the MODIS EVI-based PGT, exhibiting great potential to characterize large-scale vegetation PPT dynamics. The multi-year mean PPT mainly appeared in late July and early August with a gradually delayed trend from southeast to northwest. Meanwhile, early PPT occurred in many areas of the Tibetan Plateau (67.36%), with the overall trend advanced by about 0.89 day decade⁻¹. Interestingly, opposite trends of PPT were found with elevation below and above 2500 m asl. Preseason temperature, shortwave radiation, and precipitation jointly explained the strong spatial heterogeneity of PPT trend. All analyses help clarify the role of the ‘third pole’ in responding to environmental change and regional carbon cycles.

KEYWORDS:

• Vegetation phenology
• peak photosynthesis timing
• global change
• eddy covariance
• Tibetan Plateau

1. Introduction

Vegetation phenology refers to the annually recurring stages of vegetation development affected by climatic and other environmental factors, including the start, end, and length of the growing season. Shifts in vegetation phenology indicate the adaptation of vegetation to climate change, which can exert a substantial influence on terrestrial carbon, water cycles, and energy flow (Richardson et al. 2013). Therefore, a thorough understanding of the spatiotemporal dynamics of vegetation phenology and its response mechanism to global climate change is helpful to evaluate the terrestrial carbon sequestration potential (Berra and Gaulton 2021; Piao et al. 2019).

Traditionally, the observation of plant phenology relies on manual records in the field, providing long-term and accurate information on individual plant species (Zeng et al. 2020). Several observation networks have been established in parts of the Earth such as the USA National Phenology Network (Crimmins 2021), the Pan European Phenology Network (Templ et al. 2018), and the Chinese Phenological Observation Network (Wu et al. 2022). With the development of the Internet of Things in recent years, digital cameras as near-ground remote sensing sensors can track seasonal transitions of vegetation growth automatically and continuously with high frequency; thus, they have significant advantages in small-scale vegetation phenology observations. At present, a regional and global phenological camera network (the PhenoCam network) has been formed (Katal et al. 2022). However, the representativeness of such observations remains limited and spatially sparse, thereby constraining our understanding of spatial patterns in vegetation phenology.

Satellite-based remote sensing provides a feasible way to monitor vegetation dynamics at different spatial and temporal scales. A series of vegetation indices (VIs) that reflect the growth information via canopy greenness, such as the normalized difference vegetation index (NDVI), has been widely used to monitor phenological parameters. Several studies used MODIS-based VIs to map the peak photosynthesis timing (PPT) from regional (Wang and Wu 2019) to global scales (Gonsamo, Chen, and Ooi 2018). However, the discrepancy between these phenological metrics derived from greenness-based VIs and vegetation photosynthetic activity is large, such as the lag between the start of photosynthesis and leaf emergence, which reflects the shifts in canopy structure and ecosystem function, respectively (Jeong et al. 2017; Li et al. 2023).

The photosynthetic phenology indicates temporal dynamics of vegetation photosynthetic intensity, and it can be captured from time-series gross primary productivity (GPP) data. At the site level, GPP can be partitioned from flux-measured net ecosystem carbon exchange (NEE) based on the eddy covariance (EC) technique. At present, the global network of EC flux towers (FLUXNET) has more than 1000 sites across various ecosystems (Chu et al. 2017). Under the support of remote sensing data and ecosystem models, several satellite-based GPP products including MODIS (Running et al. 2004), BESS (Jiang and Ryu 2016), and VPM (Xiao et al. 2004) have been developed with long-term records, which provide a method to characterize vegetation photosynthetic activity. Although Yang et al. (2019) have used VPM-based GPP product to monitor PPT of temperate and alpine grasslands in China, great uncertainty remains in the accuracy of relevant satellite-based GPP products due to the algorithm, model structure, and input data, which restrict the phenological metrics of canopy photosynthesis (Wang et al. 2020). The validation results from FLUXNET showed that remote sensing-based GPP models can explain about 70% (53–78%) of site-level GPP changes (Lin et al. 2021).

Solar-induced chlorophyll fluorescence (SIF), as a proxy for vegetation photosynthesis, has been proven to be strongly correlated with GPP (Wang et al. 2023). Previous studies have shown that the SIF-derived photosynthetic phenology metrics are generally closer to phenology metrics extracted from flux-measured GPP than VI-based greenness phenology metrics (Liu et al. 2018; Zhang et al. 2020). In recent years, the emergence of global-scale SIF products with high temporal and spatial resolution based on discrete SIF soundings from satellites, such as OCO-2 and auxiliary data sets (Li and Xiao 2019; Zhang et al. 2018), has provided new insight into the extraction of large-scale photosynthetic phenology, which is significantly different from the conventional greenness phenology (Zhang et al. 2019). Previous studies about SIF-derived vegetation phenology generally focused on the start and end of growing season (Chang et al. 2019; Meng et al. 2021), whereas limited attention has been paid to PPT. As the change point of vegetation growth shifting from prosperity to decline, PPT represents the timing of maximum resource availability of vegetation and is closely linked with the transition stage in annual carbon uptake (Wang et al. 2020). Its variations including advance and delay are the response of vegetation to preseason hydrothermal conditions (preseason refers to the period prior to the timing of the phenological event) and will further affect ecosystem functions such as carbon, water, and energy balance between land surface and atmosphere (Zhang et al. 2023).

The Tibetan Plateau, known as the world's third pole, is one of the most sensitive areas to global climate change because of its unique terrain. The understanding of spatiotemporal dynamics of PPT, which reflects maximum vegetation photosynthetic activity and its environmental drivers, must be improved (Li et al. 2020). Here, we used a new global and fine-resolution SIF dataset (GOSIF, V2), along with the flux-measured GPP and MODIS EVI-based land phenology product. The goals of this study were as follows: (1) to examine the relationship between global reconstructed SIF and flux-measured GPP for different vegetation types; (2) to compare the consistency between site-level PPT from EC-based GPP (PPT_GPP) and SIF (PPT_SIF) or MODIS EVI (peak greenness timing, PGT_EVI), respectively; and (3) to evaluate the spatiotemporal dynamics of PPT in the Tibetan Plateau region and determine their response to climate factors during the past two decades. All analyses will improve our understanding and predictions of the effects of global climate change on ecosystem functions in alpine regions.

2. Methods and materials

2.1. Description of the study area

The Tibetan Plateau, located in southwestern China (Figure 1), covers a vast territory accounting for 26.8% of China's total land area, extending from the Pamir Plateau in the west, the Hengduan Mountains in the east, the southern margin of the Himalayas in the south, and the Kunlun Mountains to the north of the Qilian Mountains in the north (Yao et al. 2012). The terrain is quite complicated with an average altitude over 4000 m and a declining trend from west to east. The special geographical location at low-middle latitude and high altitude generates unique alpine climate and natural conditions. The climate is cold with multi-year mean temperature in most regions below 0°C (Chen et al. 2020). Influenced by the westerlies and the monsoons, annual precipitation ranges from over 1000 mm to below 100 mm, which decreases progressively from southeast to northwest (Wang, Pang, and Yang 2018). Solar radiation in the Tibetan Plateau is intense, which is far higher than that in other regions at similar latitudes in the mainland of China (Liu et al. 2012). In general, the vegetation types alter from forest in the southeast, to meadow and steppe in the central area, and desert in the northwest.

Figure 1. Spatial pattern of vegetation types in the Tibetan Plateau, which is reclassified from the 1:1,000,000 vegetation map of China.

2.2. Global reconstructed fluorescence data

The advent of global reconstructed SIF measurement provides an opportunity to monitor large-scale photosynthesis dynamics. Previous studies have shown a strong relationship between satellite-observed SIF from the Orbiting Carbon Observatory-2 (OCO-2) and terrestrial GPP (Li et al. 2018). However, the application of space–time discontinuous OCO-2 SIF retrievals at regional and even global scales remains a challenge. Fortunately, a global SIF product (GOSIF, V2) at a 0.05° spatial resolution and an 8-day interval is derived from individual OCO-2 SIF soundings, MODIS data, and meteorological reanalysis data (http://data.globalecology.unh.edu/). These SIF estimates from GOSIF performed well for GPP estimation at 91 EC-based flux sites globally (R² = 0.73, p < 0.001) (Li and Xiao 2019). Besides, the fine spatial resolution, globally continuous coverage, and long-term records make GOSIF useful to characterize terrestrial photosynthesis. This study used the time-series GOSIF dataset from 2001 to 2022 to extract PPT in the Tibetan Plateau region and analyze its spatiotemporal characteristics.

2.3. Eddy covariance-based GPP data

The EC technique directly measures NEE, which can then separate into GPP and ecosystem respiration (R_e). Therefore, terrestrial GPP can be used to generate photosynthetic phenology metrics and evaluate the performance of SIF-derived PPT. This study used GPP data from five EC flux sites in the Tibetan Plateau region covering three vegetation types (Table 1). GPP of each site was estimated from EC measurements by a standard procedure, including data quality control and gap filling (Li et al. 2023). The study used the night-time based flux-partitioning algorithm to separate GPP. Finally, we calculated the 8-day GPP to match the temporal interval of GOSIF and then detected the peak timing from GPP (PPT_GPP).

Table 1. Description of the flux sites used in the study.

Vegetation type	Site name	Location	MAT (°C)	MAP (mm)	ALT (m)	Dominant species	Period	References
Alpine Meadow	Damxung (DX)	91.066361E 30.497281N	1.3	476.8	4333	Stipa capillacea Carex montis-everestii Kobresia pygmaea	2004–2010	Shi et al. (2006)
	Naqu (NQ)	92.015333E 31.641833N	−1.9	380	4598	Kobresia pygmaea	2014 and 2017	Zhang et al. (2019)
Alpine Swamp	Haibei2 (HB2)	101.327333E 37.608500N	−1.7	570	3250	Kobresia tibetica Carex pamirensis Hippuris vulgaris Blysmus sinocompressus	2004–2008	Zhao et al. (2010)
	Shenzha (SZ)	88.707970E 30.946718N	0.5	388	4760	Carex Moocrofii Kobresia Tibetica	2017–2018	Qi et al. (2021)
Alpine Shrub	Haibei1 (HB1)	101.331167E 37.665167N	−1.6	560	3293	Helictotrichon tibeticum Elymus nutans	2003–2010	Zheng et al. (2013)

Notes: MAT and MAP represent mean annual temperature and mean annual total precipitation, respectively. ALT means altitude.

There are two alpine meadow flux sites (Damxung site [DX] and Naqu site [NQ]), two alpine swamp flux sites (Haibei swamp site [HB2] and Shenzha site [SZ]), and one alpine shrub flux site (Haibei shrub site [HB1]). The DX site and NQ site are situated in the hinterland of the Tibetan Plateau. The DX site is categorized as a region with plateau monsoon climate, mean annual temperature of 1.3°C, and annual mean precipitation of 476.8 mm. The NQ site has a sub-frigid monsoon climate and mean annual temperature of −1.9°C, which is the lowest among the five sites. Vegetation at the two sites is typical of alpine meadow. The HB2 site and SZ site are covered by alpine swamp. The HB2 site is located at the northeastern Tibetan Plateau, and the SZ site was built in the middle of Tibetan Plateau with an altitude of 4760 m, which is the highest among these sites. The HB1 site is close to HB2, which is covered by alpine shrub with characteristics of highland continental climate. Further details are available at the ChinaFLUX website (https://www.chinaflux.org) and Chinese Wetland Ecosystem Research Station Alliance website (https://wetland.ihb.cas.cn).

2.4. EVI-based phenology data

The EVI-based phenology data were obtained from the MODIS Land Cover Dynamics (MCD12Q2, Collection 6). Such MODIS land surface phenology product included the peak timing at 500 m spatial resolution and annual time step. The phenological events of vegetation within each year were estimated by modeling valid vegetation phases with a series of logistic functions based on the curvature change rate of smoothed EVI time series (Begue et al. 2014). Thus, this study extracted the peak layer of MCD12Q2 (https://lpdaacsvc.cr.usgs.gov/appeears/) as the MODIS-based peak greenness timing (PGT_EVI) according to the location of flux sites and relevant period of flux-measured GPP data and then aggregated it to 0.05° spatial resolution.

2.5. Auxiliary data

The climate data were obtained from China Meteorological Forcing Dataset (CMFD), which was developed from satellite-based products, reanalysis datasets, and ground-based observations (He et al. 2020). The CMFD was at high spatiotemporal resolution covering 1979 to 2018 with a spatial resolution of 0.1° and a temporal resolution of 3 h. Seven climatic factors are provided in the CMFD, including air temperature, downward shortwave radiation, and precipitation rate. This study interpolated these variables to the spatial resolution of 0.05° to match GOSIF data. We calculated monthly mean temperature, precipitation, and solar radiation to examine the response of PPT to climate dynamics.

To analyze the response of PPT to the variability of altitude in this area, the study also used the digital elevation model (DEM) data from the Shuttle Radar Topography Mission (SRTM; 90 m; https://srtm.csi.cgiar.org/).

2.6. Extraction method of PPT

We extracted PPT at two different scales, namely, site level and regional level, relying on TIMESAT software (v3.3) developed by Eklundh and Jönsson (2016). The adaptive Savitzky–Golay filter was used to lessen the noise in input data. The data were then smoothed by a least-squares fit to a quadratic polynomial, which could iteratively tighten the search window to capture rapid variations. The smoothing window was set to 5. According to a previous study (Li et al. 2023), the regional PPT was defined as a day of year halfway between the right 80% and left 80% of the input curve (Figure 2).

Figure 2. Extraction method of PPT based on TIIMESAT 3.3 from 8-day averaged GPP time series.

2.7. Statistical analysis

To assess the consistency between PPT_GPP and PPT_SIF or PGT_EVI, we calculated the root mean square error (RMSE) as follows: (1) $\begin{array}{c}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{(a_i-b_i)}^2}{n}}}\end{array}$(1) where $a_i$ or $b_i$ refers to the PPT_GPP and PPT_SIF or PGT_EVI at relevant year $i$, respectively; and $n$ refers to the number of data.

To analyze the long-term trend in PPT, we used the Mann–Kendall test and Theil–Sen slope at each pixel. If slope > 0, the PPT shifted toward autumn, and vice versa. The larger the absolute value of slope, the greater the intensity of change. Compared with the least squares regression, the Theil–Sen estimation is a non-parametric median-based slope estimator that is robust against outliers and can be computed as (2) $\begin{array}{c}\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}=\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\left({\displaystyle \frac{x_j-x_i}{j-i}}\right)\end{array}$(2) where $x_j$ and $x_i$ are the data values at times $j$ and $i$ ($1\le i<j\le 22$), respectively.

The Mann–Kendall test was implemented to evaluate the significance of the trend. As a non-parametric statistical method, the Mann–Kendall test has been widely used for trend analysis (Zhao et al. 2023). The calculation of the Mann–Kendall statistic (S) is as follows: (3) $\begin{array}{c}S=\sum _{i=1}^{n-1}\sum _{j=i+1}^n\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(x_j-x_i\right)\end{array}$(3) where $x_j$ and $x_i$ represent values at time $i$ and $j$ ($i<j$) in a time series of length $n$, respectively. The $\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(x_j-x_i\right)$ is computed as (4) $\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(x_j-x_i\right)=\{\begin{array}{c}1\mathrm{i}\mathrm{f}{x}_j>x_i\\ 0\mathrm{i}\mathrm{f}{x}_j=x_i\\ -1\mathrm{i}\mathrm{f}{x}_j<x_i\end{array}\end{array}$(4)

The variance of S was calculated using the following equation: (5) $\begin{array}{c}Var\left(S\right)={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)-{\sum }_{j=1}^n\left(t_j-1\right)\left(2t_j+5\right)}{18}}\end{array}$(5)

The standard test statistic (Z) was determined as follows: (6) $\begin{array}{c}Z=\{\begin{array}{cc}{\displaystyle \frac{S-1}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}\left(S\right)}}}& \mathrm{i}\mathrm{f}S>0\\ 0& \mathrm{i}\mathrm{f}S=0\\ {\displaystyle \frac{S+1}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}\left(S\right)}}}& \mathrm{i}\mathrm{f}S<0\end{array}\end{array}$(6)

We calculated the significance of 22 years of data and set the significance level at p < 0.1. The trend of PPT was considered significant with a 90% level of confidence when the absolute value of Z was greater than 1.64.

Wang and Wu (2019) revealed that the peak timing of China's terrestrial ecosystems was mostly determined by climatic factors within 0–1 month. Here, we first confirmed which month that vegetation photosynthetic activity in the Tibetan Plateau region peaked at and then set preseason as the month. Partial correlation analysis was performed to reveal the dominant role of climate factors during preseason (i.e. temperature, precipitation, and radiation) in the inter-annual PPT dynamics. Compared with correlation analysis, partial correlation analysis can measure the degree of association between two variables with other potential controlling variables removed. We further examined the most important environmental factor that determined PPT according to the significance level (p < 0.1).

For example, the calculation of the partial correlation coefficient of PPT to temperature when controlling precipitation and radiation is as follows: (7) $\begin{array}{c}p{r}_{12,3,4}=\sqrt{{\displaystyle \frac{R_{1\left(2,3,4\right)}^2-R_{1\left(3,4\right)}^2}{1-R_{1\left(3,4\right)}^2}}}\end{array}$(7) where variables 1–4 refer to PPT, temperature, precipitation, and radiation, respectively. $R_{1(2,3,4)}^2$ or $R_{1(3,4)}^2$ refers to the determination coefficient of the regression analysis of variable 1 (PPT) and variables 2–4 (temperature, precipitation, and radiation, respectively) or variables 3 and 4 (precipitation and radiation, respectively). The other partial correlation coefficients were calculated in a similar way.

3. Results

3.1. Relationships between SIF and site-level GPP

The scatter plots in Figure 3 showed a strong positive correlation between global reconstructed SIF and flux-measured GPP at these EC-based sites in the Tibetan Plateau region (p < 0.001). For each vegetation type, the correlation coefficient (r) ranged from 0.942 at the alpine meadow (Figure 3a) to 0.956 at the alpine swamp (Figure 3b), with the alpine shrub falling in between (Figure 3c). Furthermore, when mixing all vegetation types together, the correlation reached up to 0.877 (p < 0.001; Figure 3d). This result confirmed the strong potential of SIF in capturing photosynthetic activity for a variety of vegetation types.

Figure 3. Relationship between flux-measured GPP and global reconstructed SIF for different vegetation types in the Tibetan Plateau region at an 8-day interval: (a) alpine meadow, (b) alpine swamp, (c) alpine shrub, and (d) all vegetation types.

As illustrated in Figure 4, the seasonal dynamics of global reconstructed SIF generally agreed well with the flux-measured GPP at an 8-day interval. Both SIF and GPP exhibited similar climatology, with a single peak in the summertime. Specifically, the variability in SIF and GPP indicated that vegetation photosynthesis started around late March or April and ended in late October or early November. The PPT fell in approximately middle to late July, which was broadly consistent with the growth status of vegetation. Although a certain amount of discrepancies was observed between the global reconstructed SIF and the GPP observations at five flux sites, the time-series SIF showed highly consistent seasonality, especially in peak timing with GPP. Thus, the photosynthetic phenology in the Tibetan Plateau region could be feasibly monitored.

Figure 4. Climatology of global reconstructed SIF and the GPP observations at five flux sites in the Tibetan Plateau region. Each dot represents multi-year averaged SIF or GPP at an 8-day interval. Each curve represents fitted 8-day SIF or GPP. The vertical line represents the peak photosynthesis timing (PPT) from multi-year averaged SIF or GPP.

3.2. Site-level evaluation of PPT_SIF and PGT_EVI

To evaluate the accuracy of peak timing from global reconstructed SIF (PPT_SIF) and EVI (PGT_EVI) at the site level, we compared the consistency between PPT_SIF and PGT_EVI with PPT from the EC-based GPP (PPT_GPP) for all sites across different vegetation types (Figure 5). In comparison with PGT_EVI, PPT_SIF performed better with site-level PPT_GPP observations. The scatter plots between PPT_SIF and PPT_GPP were distributed close to the 1:1 line, with R² ranging from 0.52 to 0.90 and RMSE of 2.04 to 6.13 days, respectively. By contrast, PGT_EVI almost lagged behind PPT_GPP over all three vegetation types with relatively lower R² and higher RMSE values. Overall, PPT_SIF exhibited the best performance at the alpine shrub and swamp (RMSE = 2.04 and 2.64 days, respectively), followed by the alpine meadow ecosystem (RMSE = 6.13 days).

Figure 5. Comparisons of the peak photosynthesis timing from site-level GPP (PPT_GPP), global reconstructed SIF (PPT_SIF), and MODIS EVI (PGT_EVI) at the five flux sites in the Tibetan Plateau region: (a) and (b) alpine meadow; (c) and (d) alpine swamp; (e) and (f) alpine shrub. Solid line is the regression analysis between PPT_GPP and PGT_EVI or PPT_SIF. Dashed line represents the 1:1 line.

As shown in Figure 6, this study further compared multi-year mean PPTs including flux-measured PPT_GPP, global reconstructed PPT_SIF, and MODIS-based PGT_EVI at the five flux sites. In general, PPT_SIF exhibited good consistency with PPT_GPP at more sites than PGT_EVI. The average difference between multi-year mean PPT_GPP and PPT_SIF was 2.28 days (ranging from 0.81 to 5.08 days), which was far lower than that between PGT_EVI and PPT_GPP of 6.65 days (ranging from 1.37 to 13.56 days). Except for the SZ and NQ sites with limited site years of data, the biases between PPT_SIF and PPT_GPP were generally less than 2 days at the other sites, including the DX, HB2, and HB1 sites. However, such discrepancies were considerably larger between PGT_EVI and PPT_GPP. Therefore, the study used PPT_SIF to characterize the spatial and temporal patterns of PPT across the Tibetan Plateau region.

Figure 6. Site-level evaluation of multi-year mean peak photosynthesis timing from flux-measured GPP (PPT_GPP), global reconstructed SIF (PPT_SIF), and MODIS EVI (PGT_EVI) at the five flux sites. The dot and error bar mean multi-year mean and standard deviation, respectively.

3.3. Spatial pattern of PPT

As a result of complicated terrain factors in the Tibetan Plateau, the pattern of multi-year mean PPT exhibited strong spatial heterogeneity (Figure 7a). In general, we observed an increasing trend in PPT along the southeast–northwest direction. Early PPT was mainly located in the eastern borders of the Tibetan Plateau with low elevation, and it occurred in early July. However, PPT arrived in late July or even early August in the vast hinterland of the northwestern Tibetan Plateau area. Late PPT (>215) also occurred in the southwestern part of the Tibetan Plateau because of the effect of high elevation. An exception was a small area at the northwestern edge with low altitude where PPT arrived in mid-July. Overall, during the past two decades, the regional multi-year mean PPT mainly occurred on DOY 205–215, accounting for 56.91% of the Tibetan Plateau region. Thus, vegetation photosynthetic activity reached a peak at the end of July to the beginning of August.

Figure 7. Spatial pattern (a), latitudinal trend (b), and longitudinal trend (c) of multi-year mean peak photosynthesis timing (PPT) in the Tibetan Plateau region during 2001–2022. (d) and (e) Variations in elevation along latitude and longitude. (f) Distribution of PPT.

This study also investigated the latitudinal and longitudinal trends of PPT (Figure 7b,c) jointly with elevation (Figure 7d,e). We found an elevation-dependent spatial pattern of PPT. Overall, the variation in PPT along latitude and longitude was similar to elevation in the Tibetan Plateau area. PPT showed an obviously eastward-delayed trend in the longitude zone ranging from 82°E to 105°E, which was ascribed to the decreasing elevation. However, the variation in PPT in latitude was relatively complicated. PPT reached the earliest (during about DOY 201) at the lower and higher latitudes, which could be explained by the effect of relatively low elevation. By contrast, the lagging PPT was mainly concentrated between 28°N and 30°N with the highest elevation. It then fluctuated slightly from 30°N and 35.4°N around DOY 210. The study also found that PPT reached the latest on approximately DOY 214 at about 28.7°N.

3.4. Long-term trend in PPT

Figure 8 illustrates the long-term spatial trend in PPT over the past two decades in the Tibetan Plateau region, which also exhibited strong spatial heterogeneity. Overall, the regional PPT was advanced by 0.89 day decade⁻¹. The percentage of the advanced area was about 67.36% of the total area (Figure 8c), of which 11.37% passed the significance level (Figure 8d; p < 0.1). The advanced PPT was mainly distributed in the central hinterland of the region with high elevation. However, 32.64% of the total area experienced delayed PPT (Figure 8c), of which 2.21% passed the significance level (Figure 8d; p < 0.1), because of the sparse location in the borders of the area.

Figure 8. Spatial trend (a) and associated significant level (b) of the peak photosynthesis timing (PPT) in the Tibetan Plateau during 2001–2022. (c) Frequency histogram of the trend. (d) Proportion of different trends.

3.5. Effect of preseason climatic factors on PPT

Given the crucial role of preseason hydrothermal conditions to vegetation development, we evaluated the effects of preseason climatic factors including mean temperature, total precipitation, and total solar shortwave radiation on PPT by using spatial partial correlation analysis to account for the long-term trend in PPT during the past two decades. As illustrated in Figure 9, the dominant factors of PPT in the Tibetan Plateau area exhibited strong spatial heterogeneity. Among the three climate factors, preseason temperature was the most important factor dominating about 8.90% of all pixels, followed by solar shortwave radiation (8.15%) and precipitation (7.01%). PPT showed positive partial correlations with preseason temperature in 52.56% of all significant pixels (p < 0.1, Figure 9a), and it was mainly distributed in southwestern part of the Tibetan Plateau with high elevation. Pixels with significant negative partial correlations were much less, and they appeared in the central Tibetan Plateau area. Significant positive correlations between PPT and precipitation were mainly found in the central and northeastern Tibetan Plateau (about 61.78% of all significant pixels), indicating that additional precipitation would delay PPT (Figure 9b). The pixels with significant partial correlations between PPT and radiation were about 10.09% of the Tibetan Plateau region (Figure 9c), 53.61% of which exhibited negative effects (e.g. northeastern and central Tibetan Plateau). Thus, PPT coincided with the changes in radiation, and high levels of radiation led to the advancement in PPT. However, strong radiation may also delay PPT, which mainly occurred in the northeast low-elevation Tibetan Plateau area.

Figure 9. Spatial partial correlations between the peak photosynthesis timing (PPT) and climatic factors including (a) monthly mean temperature (Temp); (b) monthly total precipitation (Prec); (c) monthly total shortwave radiation (Srad); and (d) dominant factor in the Tibetan Plateau during 2001–2018. The significance level was set at p < 0.1.

4. Discussion

4.1. Differences among SIF, reconstructed SIF, and VIs in monitoring GPP

The study revealed a strong linear GOSIF–GPP relationship at the EC-based flux sites in the Tibetan Plateau area (Figure 3), which was confirmed by previous studies (Gao et al. 2021). The complicated relationships between SIF and other physiological processes, as well as the uncertainties in various SIF products, cannot be neglected. Although GOSIF is effective for assessing terrestrial photosynthesis (Tang et al. 2022), GOSIF is not a true SIF observation but reconstructed SIF product derived from the discrete OCO-2 SIF soundings and other auxiliary data. With the increasing demand of terrestrial vegetation photosynthesis monitoring, several reconstructed SIF data including gap-filled products, such as CSIF (Zhang et al. 2018), ${\overline{\mathrm{S}\mathrm{I}\mathrm{F}}}_{\mathrm{O}\mathrm{C}\mathrm{O}2\mathrm{\_}005}$ (Yu et al. 2019), and GOSIF (Li and Xiao 2019), and downscaled products like RSIF (Gentine and Alemohammad 2018), ${\overline{\mathrm{S}\mathrm{I}\mathrm{F}}}_{\mathrm{L}\mathrm{U}\mathrm{E}}$ (Duveiller et al. 2020), and ${\overline{\mathrm{S}\mathrm{I}\mathrm{F}}}_{005}$ (Wen et al. 2020) have been developed. The different input data and algorithm lead to uncertainties among such reconstructed SIF products, especially for the SIF signals derived from various satellite platforms. Bacour et al. (2019) compared the daily mean SIF products from OCO-2 and GOME-2, and they found disagreements in annual mean values and seasonal amplitude. Gao et al. (2021) found that the relative biases for GOSIF and CSIF are 7.14% and 11.06% in the estimated cumulative GPP across all EC towers, respectively. Moreover, uncertainty from the footprint mismatch between satellite-based reconstructed SIF product and EC-based GPP measurements still exists. However, a considerable proportion of spatial homogeneity of vegetation types in the Tibetan Plateau region could reduce the impact of such uncertainty to a certain extent. The strong linear relationship between SIF and GPP, as well as the consistent seasonal dynamics (Figures 3 and 4), also certified the feasibility to monitor regional photosynthetic phenology.

Compared with VIs, SIF is unique in capturing some physiological responses to stress that might not be reflected in VIs (Martini et al. 2022). Li et al. (2018) reported that OCO-2 SIF is more strongly correlated with tower GPP than conventional NDVI and EVI, EVI_BRDF, and the recently proposed NIRv. Zhang et al. (2021) found that EVI is not an accurate indicator of peak phenology in a warm-temperate mixed plantation. The notable separation between greenness-related VIs and photosynthesis-related SIF in a high-latitude ecosystem has been confirmed by Walther et al. (2018). Kong et al. (2020) also described asynchronization between VI-indicated greenness phenology and GPP-indicated photosynthetic phenology at 95 flux sites around the world. Yang et al. (2022) revealed the systemic inconsistency in phenological metrics between TROPOMI SIF and VIs in Tibetan Plateau grasslands. Recently, Dechant et al. (2022) found that canopy structure-related near-infrared reflectance of vegetation multiplied by incoming sunlight (NIR_VP) is a robust structural proxy for SIF and photosynthesis across scales.

Our study found that PPT_SIF was more consistent with PPT_GPP than PGT_EVI in alpine shrub, swamp, and meadow ecosystems (Figure 5), which could contribute to the mismatch between the timing of physiological parameter and timing of ecological functions, such as photosynthesis (D’Odorico et al. 2015; Jeong et al. 2017). Compared with the reflectance-based VIs, the global reconstructed SIF is more related to vegetation photosynthesis and not sensitive to cloud and atmospheric scattering, which accounted for the differences in accuracy. For large-area estimation, the remote sensing technique can only provide relatively accurate results of terrestrial ecosystems, which does not affect the long-term trend analysis. Similar to previous studies (Gao et al. 2021), the SIF exhibited strong potential in capturing the variability of GPP at the flux sites in the Tibetan Plateau region (Figure 4). Additionally, peak greenness emerged later than peak photosynthesis across years and vegetation types in circumpolar regions (Walther et al. 2018), which was also confirmed in our work. At the time when greenness is at a maximum, photosynthetic rates are decreasing. Therefore, our study provided an alternative framework to monitor global photosynthetic phenology by solely using the global reconstructed SIF dataset.

4.2. Responses of PPT to 3D environment and climate factors

Previous studies and Hopkins’ Bioclimatic Law reported the influence of 3D terrain on the spatial pattern of vegetation phenological events (Dai, Wang, and Ge 2014; Gao et al. 2020). In this study, we revealed the strong spatial heterogeneity and a northwestward delaying trend in PPT across the Tibetan Plateau region (Figure 7). Thus, our work investigated the responses of multi-year mean PPT during the past two decades to the variations in latitude, longitude, and elevation (Figure 10). Insignificant correlations were found between PPT and latitude (0.1°) with a low correlation coefficient (r = 0.121) and partial correlation coefficient (pr) lower than 0.01. However, we observed relatively stronger correlations between elevation (r = 0.568, p < 0.001) or longitude (r = −0.551, p < 0.001) and multi-year averaged PPT. When controlling the two other terrain factors, the correlation between multi-year averaged PPT and longitude weakened (pr = −0.373, p < 0.001), whereas elevation and multi-year averaged PPT still had a strong correlation (pr = 0.407, p < 0.001). Thus, the relatively strong correlation between longitude and multi-year averaged PPT was derived from the regional strong correlation between longitude and elevation. Despite complex effects exerted by longitude and elevation, our study demonstrated that these terrain factors mostly contributed to the spatial heterogeneity of PPT in the Tibetan Plateau region.

Figure 10. Distribution of multi-year averaged PPT across latitude, longitude, and elevation in the Tibetan Plateau.

Given the significant role of altitude, this study further discussed the response of multi-year averaged PPT to elevation with an interval of 100 m. As shown in Figure 11, an apparent inflection point was observed at about 2500 m asl. Specifically, when the elevation exceeded 2500 m in the Tibetan Plateau region, PPT exhibited a significant increase, which meant that peak timing was delayed by 4.99 days per 1000 m (R² = 0.87, p < 0.001). By contrast, a negative correlation between PPT and elevation was revealed when the altitude was below 2500 m, which meant that peak timing advanced by 4.49 days per 1000 m (R² = 0.80, p < 0.001). Ding et al. (2013) also reported the relationship between elevation from 2500 to 5500 m and associated phenological metrics including the start of growing season (−9 days per 1000 m), end of growing season (+1 day per 1000 m), and shortened length of growing period (−9 days per 1000 m) in the Tibetan Plateau. The opposite trend of PPT with elevation below 2500 m could be related to extensive human-induced activities (such as land management and land use difference) at the edges of the Tibetan Plateau (Ding et al. 2013).

Figure 11. Variations in averaged PPT along with the elevation gradients from the SRTM DEM dataset in the Tibetan Plateau. The averaged PPT values were calculated from multi-year PPT during 2001–2022 within each 100 m elevation belt.

In addition, climatic conditions are the direct factors affecting changes in vegetation phenology (Gong et al. 2023; Park et al. 2019; Shen et al. 2015). In this study, we systematically analyzed the influence of temperature, precipitation, and radiation on long-term PPT trends across various ecosystems. Our study revealed an overall advanced trend of PPT by approximately 0.89 day decade⁻¹ in the area (Figure 8). However, the long-term dynamics in vegetation PPT were eventually controlled by the alterations in the climate system (Wang et al. 2021). Zhang et al. (2018) also found that temperature is the dominant climatic variable driving the changes in vegetation phenology. Park et al. (2019) found a widespread warming-induced advance in PPT (−1.66 ± 0.30 day decade⁻¹, p < 0.001) across northern lands, and PPT in the boreal ecosystems showed a much closer alignment with peak radiation. Similar to the findings of Wang et al. (2015), our work also confirmed the influence exerted by precipitation on the peak timing. Indeed, dominant climate factors vary across large areas, and the impact of these climate factors on PPT can differ in various regions due to their varying sensitivity to environmental stress. For example, in the southwestern Tibetan Plateau, where spring is typically cold, a high summer temperature will increase the chance of drought and then lead to a delayed PPT (Wang and Wu 2019). Although the study explored the response of climatic factors to PPT, other environmental effects such as CO₂ fertilization and snow cover should be considered in the future.

5. Conclusion

PPT is an important but neglected phenological metric under the context of global climate change. Our study compared the global reconstructed SIF and flux-measured GPP at different EC-based flux sites in the Tibetan Plateau region. We found strong consistencies, which revealed the potential of SIF to characterize vegetation productivity. Moreover, PPT_SIF performed better than PGT_EVI in capturing the variability of PPT (PPT_GPP) at the site level. Thus, we used PPT_SIF to examine the spatiotemporal dynamics of PPT in the Tibetan Plateau region. Multi-year mean PPT mainly appeared in late July to early August, and a gradually delayed trend was revealed from southeast to northwest. This study also found that PPT over the Tibetan Plateau region experienced apparent advancement since 2000. Furthermore, the elevation exhibited a substantial effect on PPT from the perspective of 3D environment. Overall, the strong spatial heterogeneity of PPT was constrained by different climatic factors. Future work can focus on the regulatory mechanisms that underlie vegetation photosynthetic phenology. Our study contributes to the realization of interactions between global climate change and ecosystem responses in alpine regions.

Acknowledgements

We would like to express our gratitude to all scientists that have contributed to the multi-source datasets used in the study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All data can be accessed via the corresponding author upon request.

Additional information

Funding

This study was jointly supported by the Major Project of High Resolution Earth Observation System [grant number 30-Y60B01-9003-22/23] and the Special Project on National Science and Technology Basic Resources Investigation of China [grant number 2021FY100701].