Long-term, high-resolution GPP mapping in Qinghai using multi-source data and google earth engine

ABSTRACT

The terrestrial vegetation GPP of Qinghai Province is an important variable that characterizes the carbon cycling pattern. However, there is still a lack of a high-resolution GPP dataset for Qinghai Province. To address this issue, we processed all Landsat images of Qinghai from 1987 to 2021 using the GEE, and we combined multi-source auxiliary data to estimate GPP using the revised EC-LUE model. We compared our GPP dataset with flux observations to verify its accuracy. The results showed that our GPP dataset had a high correlation with the flux tower observations, with correlation coefficients of 0.984 at CF-AM site and 0.976 at CN-Ha2 site, respectively, and each site had an RMSE of $11.960\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot 16{\mathrm{d}}^{-1}$ and 12.986 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot 16{\mathrm{d}}^{-1}$, respectively. There are different deviations between our GPP dataset and the mainstream GPP datasets in various vegetation types, with the average correlation coefficient ranging from 0.431 to 0.943. By comparing with the flux observations and the related analysis, we demonstrated that our GPP dataset features better accuracy, higher spatial resolution, and more temporal coverage than mainstream GPP datasets. This study offers the first long-term high-resolution GPP dataset for Qinghai Province, and we believe that this dataset has important implications for ecological management and climate research.

KEYWORDS:

• GPP
• landsat
• revised EC-LUE model
• GEE
• NDVI

1. Introduction

Gross Primary Productivity (GPP) is the total amount of carbon dioxide ($\mathrm{C}{\mathrm{O}}_2$) fixed by green plants through photosynthesis (Beer et al. 2010). It drives many ecosystem processes and plays a crucial role in global carbon balance (Wild et al. 2022). Estimating its spatiotemporal variability is essential for assessing the carbon storage capacity of terrestrial ecosystems, predicting climate change, and studying carbon balance mechanisms (Chen et al. 2017; Ye et al. 2021; Zhang et al. 2021).

The Qinghai Plateau, also known as Qinghai Province, covers approximately 30% of the Qinghai-Tibetan Plateau and has a wide variety of natural resources, such as grasslands, forests, and wetlands. It is a sensitive area to global climate change, as its unique geographical environment influences and is influenced by the climate and environment of surrounding areas, such as the westerlies and Asian monsoons (Yang et al. 2014). Accurately quantifying the vegetation GPP of Qinghai Province and elucidating its spatiotemporal patterns are therefore critical to unveiling carbon cycling processes in this climatically sensitive region. This will also further our understanding of how vegetation productivity responds to climate change, facilitating ecological conservation and sustainable development across the plateau (Jingbin et al. 2020; Lin, Chang, and Feng 2021; Ma et al. 2018).

The eddy covariance (EC) technique is a direct method to monitor the carbon fluxes between the atmosphere and land surface, but it has some limitations. First, the measured values only reflect the carbon flux state within a small area around the EC towers, which ranges from a few hundred meters to several kilometers. EC measurements require a flat terrain and homogeneous vegetation within the footprint. Giannico et al. (2018) suggested that a high level of landscape heterogeneity would affect the accuracy of the flux measurements. Second, the number of flux tower sites in Qinghai Province is very limited. For example, only two sites in Qinghai Province are included in the FLUXNET2015 (Pastorello et al. 2020) dataset. The sparse distribution of flux tower sites leads to significant uncertainties when using the EC technique to monitor the spatial distribution and temporal trends of GPP. Therefore, researchers have applied remote sensing technology to estimate GPP by utilizing large-scale, spatially continuous land surface information (Lin, Li, and Liu 2018). Various GPP datasets have been developed to estimate global or regional GPP at different spatiotemporal scales using different data and models. For instance, Jiang et al. (2021) used a machine learning model and a reflectance fill and filter algorithm with MODIS data to produce a higher resolution (250 m, daily) GPP dataset for the contiguous United States, which improved the temporal and spatial resolution of existing MODIS GPP products. He et al. (2022) applied the Penman-Monteith-Leuning Version 2 water-carbon coupling model to estimate daily GPP and evapotranspiration data for China from 2000 to 2020, which provided valuable references for hydrological and ecological research. Zhang et al. (2017) employed an improved VPM model with MODIS and NCEP reanalysis data to estimate global GPP from 2000 to 2016 at an 8-day and 500 m spatiotemporal resolution. However, these datasets generally suffer from the trade-off between temporal and spatial resolution. Most of the existing global or regional scale GPP datasets have low spatial resolution (He et al. 2022; Jiang et al. 2021; Ma et al. 2018; Wild et al. 2022; Zhang et al. 2017). This means that the GPP estimation values in regions with strong heterogeneity at local scales, such as croplands, wetlands, pastures, and regions with complex vegetation composition, will have significant statistical uncertainty and error (Zheng et al. 2018). Moreover, the coarse spatial resolution restricts the application of these datasets in these regions, as they cannot provide sufficient spatial information to effectively differentiate between different land cover boundaries.

As the longest-running medium to high-resolution earth observation satellite (Wulder et al. 2022), Landsat has enormous potential for estimating GPP at regional and global scales. Many scholars have used Landsat data to estimate GPP. For example, Gitelson et al. (2012) found that the performance of MODIS GPP products in cropland areas was very poor due to the low spatial resolution and the lack of biome-specific parameters. Therefore, they used Landsat data to obtain some vegetation indices that represent chlorophyll content, and used these indices to estimate GPP. The results showed that the accuracy of Landsat-based GPP estimation was significantly improved, indicating that Landsat data can effectively handle the GPP estimation problem in areas with high heterogeneity. Knox et al. (2017) also used Landsat data to estimate wetland GPP, showing that Landsat can assess the carbon exchange status at a low cost. Wolanin et al. (2019) used Landsat 8 and Sentinel 2 data to estimate the GPP of C3 crops based on the radiative transfer model and machine learning algorithm. Even without using any local information corresponding to the flux sites, they were able to obtain good estimation results, demonstrating the potential of medium-to-high resolution earth observation satellites such as Landsat in mapping global crop productivity. These works have made great contributions to the field of productivity estimation at that time. However, they only applied Landsat data within a relatively small spatial range, mainly for model improvement or as auxiliary data. Applying Landsat data to obtain large-scale, spatially continuous GPP still faces many challenges, such as limited available observations due to cloud cover and the satellite's long repeat cycle, and large storage space and long processing times due to the high spatial resolution (30 m) of the data. With the emergence of Google Earth Engine (GEE), researchers have the possibility to address the aforementioned problems. He et al. (2018) used GEE in conjunction with time-series reflectance information from Landsat and MODIS to establish a pixel-level linear relationship between MODIS and Landsat NDVI, successfully estimating high-accuracy GPP data at a regional scale with an 8-day interval and 30 m spatial resolution. Robinson et al. (2018) also used Landsat data to produce a high-resolution GPP dataset covering the continental United States (CONUS). This is an innovative dataset that provides us with more spatial details of vegetation productivity on a large scale. However, this dataset is limited to the CONUS and does not cover Qinghai Province, which urgently needs vegetation productivity data. Therefore, we used GEE and a representative light use efficiency model revised EC-LUE (Yuan et al. 2007) to estimate the high spatiotemporal resolution GPP of Qinghai Province from 1987 to 2021, so that we can meet the demand for long-term and fine-resolution vegetation productivity monitoring in this region.

In this study, we leveraged Google Earth Engine to process the entire Landsat archive from 1987 to 2021 for Qinghai Province. The processed long-term Landsat-derived NDVI time series was used as input data, along with ERA5-Land reanalysis and auxiliary data. We employed the revised EC-LUE model to estimate vegetation GPP across Qinghai Province from 1987 to 2021. The main contributions of this work are: (1) We produced the first Landsat-based GPP dataset at 16-day temporal resolution and 30 m spatial resolution with a long time series (1987–2021) for Qinghai Province; and (2) We characterized the spatiotemporal variations in vegetation productivity across Qinghai Province based on the Landsat-derived GPP estimates.

2. Study area and data

2.1. Study area

Qinghai Province is located on the northeast Qinghai-Tibet Plateau, with an area of 724,000 ${\mathrm{k}\mathrm{m}}^2$. It spans from 89${}^{\circ }$E to 104${}^{\circ }$E and from 31${}^{\circ }$N to 40${}^{\circ }$N, and it has an average altitude of about 4000 m. Qinghai Province has a typical plateau continental climate, with long, cold winters and short, cool summers. The annual precipitation is low and unevenly distributed, mainly concentrated in summer and averaging between 156.6 and 669.1 mm. Temperatures range from 0.3${}^{\circ }$C to 7.2${}^{\circ }$C annually, with a large diurnal difference. Due to the influence of the Indian monsoon, Qinghai Province boasts world's most extensive and highest alpine meadows and grasslands (Li et al. 2016). The Province's unique alpine vegetation system, along with its vast glaciers and perennial snow cover, makes it an important ecological barrier and climate regulator in China (Figure 1).

Figure 1. Land cover map of Qinghai Province (2020) and the locations of four flux tower sites.

2.2. Data source

The datasets used for calculating GPP are presented in Table 1. For detailed information about each dataset, please refer to Sections 2.2.1 to 2.2.4.

Table 1. Input datasets used for generating GPP.

Dataset	Variable	Spatial resolution	Temporal resolution
Landsat Collection2 Level2 Dataset	Surface reflectance	30 m	16-day
ERA5-Land	2 m air temperature 2 m dew point temperature shortwave radiation	0.1${}^{\circ }$	Daily
SRTM digital elevation	elevation	90 m	–
CLCD	land cover type	30 m	Annual
ERSL ${\mathrm{C}\mathrm{O}}_2$	${\mathrm{C}\mathrm{O}}_2$ concentration	–	Monthly

2.2.1. Time-series landsat data

We used all Landsat-5, Landsat-7, and Landsat-8 Collection 2 Level 2 surface reflectance data in Qinghai province from 1987 to 2021 in GEE to estimate the FPAR parameter required for generating GPP. The temporal and spatial resolution of Landsat are 16-day and 30 m, respectively. All Landsat images have been orthorectified and atmospherically corrected by the USGS. The surface reflectance data of Landsat 5, 7, and 8 were generated by the LEDAPS and LaSRC algorithms (Masek et al. 2006; Vermote et al. 2016), respectively. According to the information provided by the Landsat Science Product Guide, we used a scale of 0.0000275 and an offset of $-0.2$ to convert the pixel values of surface reflectance products to actual surface reflectance. All pixels contaminated by cloud and cloud shadow were masked out using the QA band generated by the CFMASK algorithm (Foga et al. 2017). Since there are differences between different sensors, we used a conversion function to harmonize all Landsat data. Section 3.1 provides the details of this conversion function.

2.2.2. Meteorological and elevation data

The meteorological variables used in this study were obtained from the ERA5-Land daily reanalysis dataset in GEE, which provides high spatial resolution (0.1${}^{\circ }$) land information. The 0.1-degree spatial resolution daily 2m air temperature ($T_{2\mathrm{m}}$), 2 m dew point temperature ($T_d$), and shortwave radiation (Srad) were used to calculate the temperature and moisture stress factors on maximum light use efficiency of vegetation in the surrounding environment. The Shuttle Radar Topography Mission (SRTM) digital elevation Version 4 data with a spatial resolution of 90 m was resampled to 0.1-degree resolution using bilinear interpolation. The resampled elevation data, along with $T_{2m}$ and $T_d$, was used to calculate the vapor pressure deficit (VPD). The specific calculation process can be seen in 3.2 Meteorological Data Processing.

2.2.3. Land cover and $\mathrm{C}{\mathrm{O}}_2$ data

The land cover data used in this study is the 30 m spatial resolution land cover classification dataset (CLCD) developed by Yang and Huang (2021). This dataset is the first high-resolution annual land cover dataset in China developed using machine learning and cloud computing platforms. Moreover, it provides higher overall accuracy than commonly used classification datasets such as MCD12Q1 and GlobeLand30.

The carbon dioxide($\mathrm{C}{\mathrm{O}}_2$) concentration data sourced from the Earth System Research Laboratory (ERSL) (https://gml.noaa.gov/ccgg/trends/global.html) was used to estimate the $\mathrm{C}{\mathrm{O}}_2$ stress factor on maximum light use efficiency.

2.2.4. Validation and reference data

This study collected GPP flux measurements from four eddy covariance sites obtained from the FLUXNET2015 dataset (https://fluxnet.org/data/fluxnet2015-dataset/) and ChinaFlux (http://www.chinaflux.org). These four sites are CN-Ha2 and CN-HaM from the FLUXNET2015 dataset, and CF-SHR and CF-AM from ChinaFlux. Table S1 provides detailed information about these sites. These datasets have been widely used for vegetation ecophysiological mechanism studies and various model training and validation (Wang et al. 2021; Zhou et al. 2021). Creators of the FLUXNET2015 dataset used unified processing method and quality control to obtain the data, which ensured the consistency and comparability of the data from different sites (Pastorello et al. 2020). Contributors of ChinaFlux also used standardized quality control and processing methods to ensure high confidence of the flux observation data (Chang et al. 2023). In the FLUXNET2015 site data, the average of the day-time (GPP_DT_VUT_REF) and night-time (GPP_NT_VUT_REF) partition methods generated GPP were selected as the flux GPP. In ChinaFlux, we used the difference between daily ecosystem respiration (RE) and net ecosystem exchange (NEE) to calculate GPP (Chang et al. 2023), as shown in Equation (1(1) $GPP=RE-NEE$(1) ). After calculating the 16-day cumulative values of the flux data at each site, we used all the data from CN-HaM and CF-SHR sites as training data, and the data from CN-Ha2 and CF-AM sites as independent validation data, to train and validate the model. (1) $GPP=RE-NEE$(1) To evaluate the consistency and difference between the GPP estimated by this study and the existing GPP datasets, we compared our estimation to three widely used GPP datasets: the MOD17A2H (version 006) dataset, the PML-V2 GPP dataset, and the GLASS-GPP dataset. The MOD17A2H dataset is a cumulative 8-day 500 m GPP dataset based on the light use efficiency model proposed by Monteith (1972). It utilizes FPAR dataset (MOD15A2H), GMAO/NASA dataset, and biome look-up table to achieve near-real-time estimation of global GPP. This dataset is the most widely used dataset (Zhang et al. 2017). The PML-V2 GPP dataset is also an 8-day 500 m GPP dataset that used a water-carbon coupled model (referred to as PML-GPP). This dataset has been shown to have higher accuracy than other mainstream GPP products through its comparisons with 95 global flux sites. The GLASS (Global Land Surface Satellite) dataset (http://www.glass.umd.edu/GPP/AVHRR) is a global long-term land surface multi-parameter dataset. The dataset contains 8-day 0.05${}^{\circ }$ GPP products (referred to as GLASS-GPP). With the above three datasets were selected as reference datasets, we accumulated the 8-day GPP estimates from these datasets over 16-day intervals to match the temporal resolution of the GPP estimates from this study. It should be noted that although these three GPP datasets are widely used, there are still some potential errors in them. Existing studies have shown that the MODIS GPP product has errors mainly caused by meteorological data and vegetation parameter data (Heinsch et al. 2006; Zhang et al. 2019), and the GLASS-GPP product, which also uses a light use efficiency model, is likely to have the same errors. The values of different versions of MODIS GPP products show significant inconsistencies in some regions, possibly due to the use of different algorithms (Peter and Messina 2019). The PML-GPP product is also affected by the MCD12 land cover data it uses, which increases the uncertainty of its model parameters and thus affects its estimation accuracy (Zhang et al. 2019). Therefore, these GPP datasets can only be used as references for our estimated GPP data, and the flux observation GPP will be used to validate our data.

3. Methods

Generating accurate high-resolution GPP for the entire Qinghai Province faces two main challenges: (1) how to generate large-area high spatial-temporal resolution NDVI data under adverse factors such as cloud cover, which causes a large amount of missing data; and (2) how to accurately estimate the downward regulation scalars of temperature, moisture, and $\mathrm{C}{\mathrm{O}}_2$ on light use efficiency (LUE). To address these challenges, we developed an effective framework for estimating GPP accurately. Figure 2 shows the specific framework for calculating GPP. The first step is to generate high resolution NDVI data, which involves the following procedures: calculating NDVI from Landsat surface reflectance data, masking out clouds, shadows, and snow using the QA band in the Landsat surface reflectance products, synthesizing 16-day NDVI data, and filling and reconstructing NDVI data. The second step is to use the reconstructed NDVI data and meteorological data as inputs for the revised EC-LUE model to obtain GPP data. The third step is to compare and validate the estimated GPP with other GPP datasets and flux tower observations.

Figure 2. The flowchart of calculating GPP using the revised EC-LUE model.

3.1. NDVI filling and reconstruction

The vegetation index is a key input variable in the revised EC-LUE model, which reflects the vegetation's ability to absorb solar radiation through a linear transformation. We use Landsat surface reflectance data to calculate NDVI, which is defined as follows: (2) $NDVI=\frac{S{R}_{nir}-S{R}_{red}}{S{R}_{nir}+S{R}_{red}}$(2) where $S{R}_{nir}$ represents the surface reflectance in the near-infrared band of the corresponding sensor of Landsat satellites, and $S{R}_{red}$ represents the surface reflectance in the red band of the corresponding sensor of Landsat satellites.

To ensure the consistency of the time series NDVI, we apply a simple transformation function (Roy et al. 2016) to adjust the performance differences among Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI sensors. This function is originally designed in consideration of differences between Landsat-7 ETM+ and Landsat-8 OLI sensors. Since Landsat-5 TM and Landsat-7 ETM+ sensors have similar performance (Claverie et al. 2015), we extend the function to Landsat-5 TM sensor as well. The specific transformation function is as follows: (3) $NDV{I}_{OLI}=0.0235+0.9723\times NDV{I}_{TM,ETM+}$(3) where $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{O}\mathrm{L}\mathrm{I}}$ and $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{T}\mathrm{M},\mathrm{E}\mathrm{T}\mathrm{M}+}$ are the NDVI values calculated from Landsat-8 OLI and Landsat-5 TM or Landsat-7 ETM+ sensor data, respectively. One challenge we encountered was that the low temporal resolution of Landsat (about 16-day) and the atmospheric conditions limit the availability of Landsat-NDVI observations. To eliminate the effects of clouds and other factors on NDVI data, we first use the ‘QA_PIXEL’ band of Landsat surface reflectance products to evaluate the pixel quality and mask out clouds, shadows, and snow in the images. Then, we use the maximum value composite method to synthesize 16-day NDVI data for Qinghai Province starting from January 1st of each year. For the missing pixels caused by the masking process, we use the filling method proposed by Robinson et al. (2017) to select ‘clear’ pixels from the past five years based on the ‘QA_PIXEL’ band, and we use the median composite method to synthesize them. We choose a time scale of five years because the time span of the El Nio-Southern Oscillation (ENSO) is about 5 years, and ENSO has a huge impact on the biosphere. Selecting NDVI filling values beyond five years may introduce more uncertainty to the NDVI reconstruction (Running and Zhao 2021). To further smooth the filled NDVI data and reduce the noise, we use the harmonic regression model (Wilson, Knight, and McRoberts 2018) to decompose the NDVI time series into a limited number of harmonic terms to model NDVI. For further understanding of this method, please consult the supporting material.

3.2. Meteorological data processing

To match the temporal resolution of other data, we computed the 16-day average of daily $T_{2m}$ and $T_d$ in ERA5-Land, and the 16-day sum of Srad. We resampled all meteorological data to a spatial resolution of 30 m using bilinear interpolation. Since ERA5-Land does not include VPD data, we calculated VPD using the following formula (Fang et al. 2022): (4) $VPD=SVP-AVP$(4) (5) $SVP=6.112\times f_w\times {\exp }^{\frac{17.67T_{2\mathrm{m}}}{T_{2\mathrm{m}}+243.5}}$(5) (6) $AVP=6.112\times f_w\times {\exp }^{\frac{17.67T_d}{T_d+243.5}}$(6) where SVP and AVP are the saturated and actual vapor pressure in hPa, respectively. $T_d$ and $T_{2\mathrm{m}}$ are the dew point and air temperature in ${}^{\circ }$C, respectively. (7) $f_w=1+7\times {10}^{-4}+3.46\times {10}^{-6}{P}_{mst}$(7) (8) $P_{mst}=P_{msl}\times {\left(\frac{T_{2\mathrm{m}}+273.16}{T_{2\mathrm{m}}+273.16+0.0065Z}\right)}^{5.625}$(8) where $P_{msl}$ represents the mean sea level air pressure (1013.25 hPa), and Z represents the elevation (m).

3.3. Revised EC-LUE model

This paper uses the revised EC-LUE model to estimate GPP. However, our implementation of EC-LUE has one major difference compared with the original model, when considering the constraints of temperature and moisture on the maximum light use efficiency of vegetation, we did not follow the the original model's method of directly taking the minimum value of the temperature and moisture scalars. Instead, we assumed that temperature and moisture jointly affect the maximum light use efficiency of vegetation, hence we multiplied the temperature scalar with the moisture scalar when computing LUE of vegetation. This is because Qinghai Province is located on the Qinghai-Tibet Plateau, where the vegetation growing season is short, and the temperature is low in spring and autumn. If the original model's calculation framework is used, it would easily result in zero values in certain periods of estimated GPP, in particular spring and autumn. However, by observing the temporal variation of GPP at flux sites in Qinghai Province, we found that the GPP in spring and autumn account for about 10% of the annual GPP. Therefore, we made some modifications to the model to better adapt to the local conditions. (9) $GPP=FPAR\times PAR\times T_s\times W_s\times C_s\times {ϵ}_0$(9) where FPAR is the fraction of absorbed photosynthetically active radiation, which is derived from the reconstructed NDVI; PAR is the photosynthetically active radiation ($\mathrm{M}\mathrm{J}\cdot {\mathrm{m}}^{-2}$), which is obtained by multiplying Srad by 0.45; ${ϵ}_0$ represents the maximum light use efficiency of vegetation ($\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{M}{\mathrm{J}}^{-1}$). $T_s$, $W_s$, and $C_s$ are the temperature, moisture, and carbon dioxide downward regulation scalars for ${ϵ}_0$, which range from 0 to 1. (10) $FPAR=1.24\times NDVI-0.168$(10) (11) $T_s=\frac{(T-T_{min})(T-T_{max})}{(T-T_{min})(T-T_{max})-(T-T_{opt}{)}^2}$(11) (12) $W_s=\frac{VP{D}_0}{VPD+VP{D}_0}$(12) (13) $C_s=\frac{C_i-\phi }{C_i+2\phi }$(13) where NDVI is the normalized difference vegetation index obtained from Landsat data; T represents the average air temperature (${}^{\circ }$C); $T_{min}$, $T_{max}$, and $T_{opt}$ are the minimum, maximum, and optimum air temperature for vegetation photosynthesis, respectively. The higher average elevation and plateau continental climate on the Qinghai Plateau may lead to significant differences in the values of $T_{min}$, $T_{max}$, and $T_{opt}$ compared to other regions. Therefore, following the settings in previous studies (Huang et al. 2019; Wang et al. 2012), we set $T_{min}$, $T_{max}$, and $T_{opt}$ to 0${}^{\circ }$C, 35${}^{\circ }$C, and 13${}^{\circ }$C, respectively. $VP{D}_0$ represents the half-saturation coefficient of the vapor pressure deficit (kPa). $C_i$ is the intercellular carbon dioxide concentration (ppm); φ is the compensation point of carbon dioxide in the absence of dark respiration, which can be calculated using the following formula (Bi et al. 2022). (14) $\phi =4.22\times {\exp }^{\frac{37830(T_{2\mathrm{m}}-298.15)}{298.15R{T}_{2\mathrm{m}}}}$(14) (15) $C_i=C_a\times \chi$(15) (16) $\chi =\frac{\xi }{\xi +\sqrt{VPD}}$(16) (17) $\xi =\sqrt{\frac{356.51k}{1.6{\eta }^{\ast }}}$(17) (18) $K=K_c\times (1+\frac{P_0}{K_0})$(18) (19) $K_c=39.97\times {\exp }^{\frac{79.43(T_{2\mathrm{m}}-298.15)}{298.15R{T}_{2M}}}$(19) (20) $K_0=27840\times {\exp }^{\frac{36.38(T_{2\mathrm{m}}-298.15)}{298.15R{T}_{2M}}}$(20) where $C_a$ is the atmospheric $C{O}_2$ concentration (ppm); χ is the ratio of leaf-internal to ambient $C{O}_2$; K is the Michaelis-Menten coefficient of Rubisco; ${\eta }^{\ast }$ represents the water viscosity at 25${}^{\circ }$C (0.8903); $K_c$ and $K_0$ are the Michaelis-Menten constants of Rubisco for $C{O}_2$ and $O_2$, and $P_0$ represents the partial pressure of $O_2$ (21 kPa). R is the molar gas constant (8.314 $\mathrm{J}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}{\mathrm{K}}^{-1}$). To optimize the maximum light use efficiency and $VP{D}_0$ in the model, we used the Markov Chain-Monte Carlo method, which is widely applied to light use efficiency model parameter estimation. Based on previous experience (Bagnara et al. 2015), we set the prior distribution of the parameters as uniform distribution, and set the range of maximum light use efficiency and $VP{D}_0$ to 0–10 (since in most literature, these two parameters are set within 0–10). We chose the Metropolis-Hastings algorithm as the MCMC sampler (see Yuan et al. (2012) for a more detailed description). We used training data (CN-HaM and CF-SHR site data) corresponding to grassland and shrub vegetation cover types, so we only optimized the model parameters for these two vegetation types. We set the model parameters for other vegetation types based on the article by Yuan et al. (2019). The classification system used in the CLCD dataset is different from that used in the article by Yuan et al. (2019). In the CLCD dataset, vegetation is classified as cropland, forest, shrub, grassland, and wetland, while in the article by Yuan et al. (2019), forests are classified as deciduous broadleaf forest, evergreen needleleaf forest, evergreen broadleaf forest, and mixed forest. We set the parameter values for these four types of forests to the average value for forests. The specific parameter values are shown in Table 2.

3.4. Data validation and comparison

To verify the accuracy of the estimated GPP in this study, we calculated the root mean square error (RMSE), bias, and correlation coefficient (r) between the estimated GPP and the flux tower observed GPP. Similarly, we computed the RMSE, bias, and r between PML-GPP, GLASS-GPP, and MOD17A2H GPP and flux tower observed GPP for comparative analysis with our estimation. We also calculated the statistical metrics between the estimated GPP and the existing GPP products to evaluate their consistency. To match the different spatial resolutions of the GPP data, we resampled the high spatial resolution images using bilinear interpolation to the low spatial resolution before calculating these statistical metrics. When compared with PML-GPP and MOD17A2H, we resampled the estimated 30 m GPP to 500 m, and when compared with GLASS-GPP, we resampled the estimated 30 m GPP to 0.05${}^{\circ }$ (about 5566 m). The formulas for calculating these three statistical metrics are shown below: (21) $r=\frac{cov(\hat{\rho },\rho )}{\sqrt{VAR\left(\hat{\rho }\right)\times VAR\left(\rho \right)}}$(21) where $\hat{\rho }$ represents the estimated GPP; ρ represents the flux tower observed GPP; $cov(\hat{\rho },\rho )$ represents the covariance between the estimated GPP and the flux observation GPP; $VAR\left(\hat{\rho }\right)$ and $VAR\left(\rho \right)$ are the variance of the estimated GPP and the flux tower observed GPP, respectively. (22) $RMSE=\sqrt{\frac{1}{m}\sum _{i=1}^m(\widehat{{\rho }_i}-{\rho }_i{)}^2}$(22) (23) $Bias=\frac{\sum _{i=i}^m(\widehat{{\rho }_i}-{\rho }_i)}{m}$(23) where m represents the total number of observations; i represents the sequence number of the observation; $\widehat{{\rho }_i}$ and ${\rho }_i$ respectively represent the estimated GPP and the observed GPP of the i-th sample point.

Table 2. Revised EC-LUE model parameters for different vegetation types.

Vegetation Type	Cropland	Forest	Shrub	Grassland	Wetland
${ϵ}_0\left(\mathrm{g}\mathrm{C}{\mathrm{m}}^{-2}\mathrm{M}{\mathrm{J}}^{-1}\right)$	2.85	3.14	4.51	3.38	2.65
$VP{D}_0\left(KPa\right)$	0.89	0.48	0.46	0.23	1.16

4. Results

4.1. GPP comparison and validation

4.1.1. Validation against fluxnet GPP

We evaluated the accuracy of the estimated GPP in this study (referred to as $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$) by comparing it with flux tower observations and MOD17A2H, PML-GPP, and GLASS-GPP datasets. Figure 3 shows that most of the $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ data points at a 16-day time step are close to the 1:1 line with the flux tower observations, indicating a high level of agreement between the two. At the CF-AM site, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ has the second highest r (0.984) and slightly higher RMSE (11.960 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot 16{\mathrm{d}}^{-1}$) than PML-GPP, followed by MOD17A2H and GLASS-GPP, which have similar r values (around 0.98) but higher RMSE values. At the CN-Ha2 site, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ has a slightly lower r (0.976) than the other three datasets. The RMSE is 12.986 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot 16{\mathrm{d}}^{-1}$, higher than PML-GPP and GLASS-GPP, but lower than the MOD17A2H dataset. Overall, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ has the second highest correlation coefficient (0.978) and almost the lowest RMSE (12.312 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot 16{\mathrm{d}}^{-1}$). Figure 4 illustrates the time series data of remote sensing GPP datasets and flux tower observations. From the figure, we can see that $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ exhibits a consistent seasonal variation with the flux tower observations. Except for MOD17A2H, which underestimates GPP at both sites, there are relatively small differences among the different datasets. Overall, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ demonstrates similar accuracy to the other datasets. The high correlation coefficient and low RMSE indicate a small deviation between $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ and the flux tower observations.

Figure 3. FLUXNET GPP against $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$, MOD17A2H, PML-GPP and GLASS-GPP with 16-day interval, the first row corresponds to the data for the CF-AM site, the second row corresponds to the data for the CN-Ha2 site, and the third row corresponds to the data for both sites.

Figure 4. 16-day temporal variations of different GPP datasets and flux tower observations at the two sites in Qinghai Province: (a) CF-AM site and (b) CN-Ha2 site.

4.1.2. Comparison of different GPP datasets

The limited number of flux tower sites in Qinghai Province introduces some uncertainties in validation. Therefore, we also compared $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ with other datasets for different vegetation types to quantitatively assess the correlation and differences between $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ and the reference datasets (Figure 5). We found that $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ exhibits strong linear relationships with the other three datasets in most vegetation types. Specifically, in grassland types, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ has the highest average r (0.943) with PML-GPP, followed by MOD17A2H (0.891) and GLASS-GPP (0.868). The r values are concentrated, indicating high similarity between $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$-derived grassland productivity and the other datasets for the period of 2001-2018, with low RMSE and Bias suggesting minimal deviation. Overall, in the Shrub type, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ has lower r values with the other datasets, particularly with the MOD17A2H dataset, which has an average r of only 0.431. Both RMSE and Bias are higher for the MOD17A2H dataset compared to the PML-GPP and GLASS-GPP datasets. Among all reference datasets, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ shows the highest correlation and the least difference with PML-GPP. Specifically, the r values between PML-GPP and $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ range from 0.671 to 0.943 across different vegetation types, all higher than those of MOD17A2H and GLASS-GPP datasets. The RMSE and Bias range from 123.750 to 281.739 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-1}$ and −102.935 to 121.518 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-1}$, respectively, which are also the lowest compared to the other two reference datasets. Overall, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ exhibits high consistency with flux tower observations and reference datasets, indicating its high accuracy and suitability for applications.

Figure 5. Statistical metrics of the annual cumulative $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ compared to PML-GPP (PML), MOD17A2H (MOD), and GLASS-GPP (GLASS) from 2001 to 2018. These statistics were calculated based on 10,000 randomly generated sample points within each vegetation type.

We also compared $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ with other remote sensing GPP datasets for different regions and seasons in Qinghai Province. As shown in Figure 6, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ maintains consistent quarterly variations with the reference datasets in most regions, particularly in the northeastern, eastern, and southeastern parts of Qinghai Province. In the western part of the Sanjiangyuan region, where precipitation is low and vegetation is sparse, GPP values are also low, resulting in differences between $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ and the reference datasets in that area. Overall, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ shows high correlation with PML-GLASS, MOD17A2H, and GLASS-GPP, with average r values of 0.879, 0.901 and 0.839, respectively. Due to the differences in input data, model structures, and model parameterization methods, there may be variations in the magnitude of GPP estimated by different models (Zheng et al. 2018). Therefore, a better approach is to compare the interannual variations of different GPP anomalies (Figure 7). We found that all GPP datasets exhibit consistent interannual trends during the period from 1987 to 2021. From this figure, we can observe that each dataset shows a continuous increase in GPP in Qinghai.

Figure 6. Spatial pattern of pearson's r between $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ and (a) GLASS-GPP, (b) MOD17A2H, and (c) PML-GPP at 16-day intervals from 2001 to 2018. The corresponding frequency distributions of the pearson's r are shown in the (a), (b) and (c).

Figure 7. Time series of GPP anomalies for the $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ and PML-GPP, MOD17A2H, GLASS-GPP dataset

4.2. Spatial and temporal distribution patterns of GPP in qinghai province

We analyzed the seasonal variations of $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ in Qinghai Province. As shown in Figure 8, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ exhibits strong seasonal fluctuations. During the summer months (June to August), which are characterized by suitable temperature and rainfall for vegetation growth, GPP reaches its peak, accounting for the highest proportion (87.11%) of the annual value. In contrast, during the winter months (December to February), characterized by extremely low temperature and rainfall, vegetation photosynthetic capacity is severely weakened, resulting in near-zero GPP. GPP values in spring and autumn are similar, accounting for 6.14% and 6.75% of the annual value, respectively. From 1987 to 2021, the average annual GPP in Qinghai Province was 424.179 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-1}$. Over the course of 35 years, the annual total GPP shows a decreasing trend during the years of 1987 and 1992, followed by gradual increases in 1993–2005, 2006–2015, and a rapid increase after 2015. The peak GPP value was observed in 2020, reaching 537.64 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-1}$. Overall, GPP exhibited an increasing trend with a growth rate of 4.42 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-2}$ ($p<0.05$).

Figure 8. Distribution image of 16-day GPP in Qinghai Province. Due to the varying number of days in each year from 1987 to 2021, an index is used to represent the regional average cumulative GPP within each 16-day period of the year in Qinghai Province. For example, in the year 1987, a grid with an index of 1 represents the cumulative GPP regional average from January 1, 1987, to January 16, 1987, in Qinghai Province.

As Figure 9 illustrates, GPP in Qinghai Province exhibits a clear spatial heterogeneity. The eastern and southeastern regions of Qinghai Lake and the southern region of the Sanjiangyuan area have high GPP values. In the northeastern part of Qinghai Province, such as the Chaidamu Desert region, and the western part of the Sanjiangyuan area, GPP is extremely low due to factors such as scarce precipitation and severe soil salinization. Generally, all GPP datasets decreases gradually from south to north and increases from west to east. Regarding the temporal trend of GPP changes, GPP in Qinghai Province showed a general increase from 1987 to 2021, where 90.47% of the vegetation areas feature increasing GPP. The areas with the most pronounced GPP increase are near Qinghai Lake-Qilian Mountains, the agricultural and pastoral production areas in the western part of Qinghai Province, and the southern and southeastern regions of the Sanjiangyuan area, with a rate of up to 20 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-2}$. However, some areas in the central and southern parts of the Sanjiangyuan region exhibited a slight decreasing trend in GPP. The increasing trend of GPP from our observation is consistent with previous studies (Shen et al. 2015; Wei et al. 2022; Zhong et al. 2019), indicating that vegetation productivity in Qinghai Province, located on the Qinghai-Tibet Plateau, has been continuously increasing due to the warming and humidification of the plateau and the prolongation of the vegetation growing season.

Figure 9. Temporal trends and spatial distribution of GPP in Qinghai Province from 1987 to 2021. The White areas in (b) indicate no significant trend (p<0.05).

5. Discussion

5.1. Advancements and uncertainties

In this paper, we developed the first long-term, high-resolution 16-day GPP dataset, named $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$, so that we can aid Qinghai Province in addressing the need for accurate monitoring of vegetation productivity in the context of climate change. Compared to reference GPP datasets, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ featured the highest correlation and smaller deviations with flux tower observations. Moreover, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ provided vegetation productivity information at a higher spatial resolution, and exhibited a high level of spatiotemporal consistency with the reference datasets. The higher spatial resolution of $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ allows us to observe vegetation changes at finer scales, as Figure 10(e) illustrates, the values of $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ are significantly higher than those of other datasets, indicating that $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ can accurately capture the productivity changes of small patches of green vegetation within the vast Mugetan sandy land. In the data processing procedure, we employed transformation functions, data filling, and smoothing algorithms to adjust and reconstruct the time-series NDVI, accounting for the effects of atmospheric conditions, aerosols, and sensor variations on Landsat imagery. We also integrated the effects of hydro-thermal conditions and atmospheric $\mathrm{C}{\mathrm{O}}_2$ concentration on vegetation light-use efficiency by using the revised EC-LUE model, achieving accurate estimation of GPP in Qinghai.

Figure 10. Annual GPP of Mugetan sandy land in different spatial resolutions in 2017: (a) 0.05${}^{\circ }$; (b) 500 m; (c) 500 m; (d) 30 m. (e) represents the 16-day GPP time series data at the black dots in (a–d).

Despite the significant improvement in accuracy and the provision of long-term consistent high-resolution GPP observations, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ still has some limitations. Firstly, the input data used for GPP generation has inherent uncertainty. For example, Landsat NDVI can mitigate the effects of topography, atmosphere, and BRDF, but there are some known issues with Landsat products in the CFMASK algorithm used in this product struggles to identify thin clouds, and the detection of cloud shadows has not been fully validated (Foga et al. 2017), leading to some uncertainties in the time series NDVI. To minimize these errors, we employed data filling and smoothing algorithms. Secondly, the 30 m resolution landcover used in this paper is an annual data, but due to insufficient Landsat data, there was only one image before 1990, which nominally belongs to 1985. Therefore, we had to use the same landcover data for producing the GPP for 1987–1989, which may result in inaccurate classification information. This kind of error may also propagate to the GPP through the biome-specfic parameters in the revised EC-LUE model. Thirdly, the structure of revised EC-LUE model only considers the effects of hydro-thermal conditions and atmospheric $\mathrm{C}{\mathrm{O}}_2$ concentration on vegetation photosynthesis. Other factors affecting photosynthesis, such as soil nitrogen, potassium content and other soil properties, can be added to the model. However, due to the difficulty of obtaining soil properties and their mismatches of temporal coverage with the $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$, this limitation can only be improved by conducting research.

5.2. Future and benefits

In this paper, we successfully used tens of thousands of Landsat images over more than 30 years to map high-resolution vegetation productivity in Qinghai Province, which demonstrates the reliability of our proposed GPP processing workflow. Moreover, our work also provides useful reference for future high-resolution vegetation productivity estimation at continental and even global scale. Currently, there are many methods to monitor vegetation change information at high resolution, such as the widely used spatiotemporal fusion algorithms (Gao et al. 2006; Zhu et al. 2010) and the filtering methods with gap-filling and fusion (Cao et al. 2022; Chen et al. 2021). However, these methods have some drawbacks that are hard to overcome. First of all, these methods usually require fusing high spatial and low temporal resolution data with low spatial and high temporal resolution data. Generally speaking, low spatial and high temporal resolution data are chosen from MODIS or VIIRS data, since data with longer temporal coverage such as GIMMS NDVI have a huge difference in spatial resolution (about 8000 m) with Landsat (30 m), Sentinel (10 m) and other high spatial resolution data, which may lead to high uncertainty. The consequence of choosing MODIS and VIIRS data is that these methods cannot extend the observation period of high-resolution data to pre-2000. In the GPP processing workflow proposed in this paper, only Landsat satellite data is used, which is the longest global satellite observation record (Wulder et al. 2022). This allows us to extend the high-resolution GPP observation period to the 1980s, which is crucial for some climate change research and ecological management in climate-sensitive areas. Apart from the limitations discussed before, these methods are usually very time-consuming, such as spatiotemporal fusion algorithms. Although some scholars have successfully implemented this algorithm on the GEE platform (Nietupski et al. 2021), due to the low-end server cluster as the underlying foundation of GEE, and the limitation of the amount of data that can be passed into a single server in these clusters (Gorelick et al. 2017), GEE cannot effectively execute spatiotemporal fusion algorithms, and the computing efficiency is slow. The parallel processing mode of GEE itself makes the pixel-by-pixel calculation process used in this paper very efficient. This feature allows us to extend the processing workflow in this paper to any location in the world, so that we can help people quickly and conveniently obtain high-resolution GPP data.

6. Conclusion

We developed a long-term (1987–2021) and high spatiotemporal resolution (16-day, 30 m) GPP dataset for Qinghai Province using the Google Earth Engine (GEE) platform. By integrating multiple data sources, vegetation index filling and smoothing algorithms, and revised EC-LUE model, our $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ estimates agreed well with flux tower observations, with a high average correlation coefficient of 0.978 and the second lowest RMSE of 12.312 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot 16{\mathrm{d}}^{-1}$. This indicates a small deviation and consistent temporal variation between $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ and flux tower observations. Since flux tower sites are limited in Qinghai Province, we compared $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ with reference datasets under different vegetation types. Overall, $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ showed the highest correlation coefficients (0.868–0.943), lowest RMSE (123.750–276.978 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-1}$), and bias (−61.776–193.949 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-1}$) in grassland areas, which constitute the majority of Qinghai Province. These results demonstrate the high consistency of the $\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{L}}$ estimates with the reference datasets in grassland type. The validation results indicate that our estimated GPP for Qinghai Province is similar in accuracy to existing GPP datasets.

We analyzed the estimated GPP and found that vegetation productivity in Qinghai Province had a fluctuating growth over the past 35 years, with a growth rate of 4.42 $\mathrm{g}\mathrm{C}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{y}{\mathrm{r}}^{-2}$ (p<0.05). Approximately 90% of the vegetation areas exhibited an increasing trend in GPP. Qinghai Province's GPP had clear seasonality and spatial consistency, influenced by factors such as altitude and hydro-thermal conditions. Summer accounted for approximately 87% of the annual GPP, and most of the GPP was concentrated in the eastern and southeastern regions of Qinghai Province, where hydro-thermal conditions were favorable. In contrast, the GPP values were extremely low in the large desert areas of western and northwestern Qinghai Province.

Compared with existing datasets, this GPP dataset has some advantages or features. For example, our dataset is the first high resolution GPP dataset for Qinghai Province. Although some researchers have estimated high-resolution GPP using GEE (such as He et al. (2018); Robinson et al. (2018)), these datasets usually do not cover Qinghai Province, and our dataset has higher accuracy and reliability in Qinghai Province because we used independent site data located in Qinghai Province for training and validation. Secondly, our dataset has a longer time coverage, and the dataset is housed in GEE, which is convenient for researchers to analyze it. In summary, our GPP dataset can provide reliable observations of long-term GPP in Qinghai Province at a finer scale.

Disclosure statement

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

Data availability statement

The GPP dataset of this study are openly available in https://code.earthengine.google.com/3443ab7c12cf04448eda52376a44cc07

Additional information

Funding

This work was supported by the [National Natural Science Foundation of China] under Grant [U20A20981]; [Natural Science Foundation of Jiangsu Province] under Grant [BK20190785].