Improving the Gross Primary Production Estimate by Merging and Downscaling Based on Deep Learning

A reliable estimate of the gross primary productivity (GPP) is crucial for understanding the global carbon balance and accurately assessing the ability of terrestrial ecosystems to support the sustainable development of human society. However, there are inconsistencies in variations and trends in current GPP products. To improve the estimation accuracy of GPP, a deep learning method has been adopted to merge 23 CMIP6 data to generate a monthly GPP merged product with high precision and a spatial resolution of 0.25◦, covering a time range of 1850–2100 under four climate scenarios. Multi-model ensemble mean and the merged GPP (CMIP6DL GPP) have been compared, taking GLASS GPP as the benchmark. Compared with the multi-model ensemble mean, the coefficient of determination between CMIP6DL GPP and GLASS GPP was increased from 0.66 to 0.86, with the RMSD being reduced from 1.77 gCm−2d−1 to 0.77 gCm−2d−1, which significantly reduced the random error. Merged GPP can better capture long-term trends, especially in regions with dense vegetation along the southeast coast. Under the climate change scenarios, the regional average annual GPP shows an upward trend over China, and the variation trend intensifies with the increase in radiation forcing levels. The results contribute to a scientific understanding of the potential impact of climate change on GPP in China.


Introduction
Global climate change, especially the impact of climate warming on the human living environment, has attracted more and more attention [1]. With the continuous development of the global industry, the use of fossil fuels is increasing. This process releases large amounts of carbon dioxide and changes the proportions of greenhouse gases, causing global climate problems such as the intensification of the greenhouse effect. The gross primary productivity (GPP), formed by terrestrial vegetation through photosynthesis to fix CO2 and energy in the biosphere, is the initial stage of entering the carbon cycle process and the basis of the ecosystem carbon cycle [2]. Terrestrial ecosystems can influence global climate change through the carbon cycle [3]. Therefore, the research on the carbon cycle of global terrestrial ecosystems has received widespread attention and has become a hot spot in global change research [4].
An essential indicator of terrestrial ecosystem carbon cycle research is GPP. The gross primary productivity, that is, the total photosynthetic uptake or carbon assimilation by plants, is a key component of terrestrial carbon balance [5]. It is the main driving factor of the terrestrial carbon sink, which is responsible for absorbing about 30% of anthropogenic carbon dioxide emissions [6]. In general, it can be said that GPP plays a vital role in the global carbon cycle [7][8][9][10]. Characterizing the spatiotemporal changes and trends of GPP is essential for a deeper understanding of the carbon cycle between terrestrial ecosystems and excellent results [53]. It establishes the relationship between reference data and multi-model datasets in the historical period to generate a merged GPP product with the distribution characteristics of reference data. Furthermore, GPP-merged data in the projection period under multiple climate scenarios can be obtained based on 23 CMIP6 model projections if the relationship built during the historical period is applied to the projection period. The evaluation in the historical period shows that the quality of the merged data is significantly improved compared with the multi-model ensemble mean, indicating that the merged data are reliable. Therefore, the variation trends in GPP under multiple climate scenarios in the future were further analyzed based on the merged data.
This study aims to produce an improved high-resolution merged GPP dataset by adopting a new method to combine the advantages of 23 CMIP6 GPP datasets, without relying on any prior knowledge, to analyze the long-term variation trends in GPP under multiple climate scenarios in the future based on the improved merged data. Based on the historical CMIP6 GPP data, the deep learning method called Ynet [52] neural network is used to train the deep learning model with merging and downscaling functions. The 23 groups of monthly GPP datasets are integrated into this high-resolution long-series dataset called CMIP6 DL GPP. The period and the spatial resolution are 1850-2100 and 0.25 degrees. The future projection period includes four climate change scenarios, such as SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8. 5. The performance of the merged product has been evaluated based on eddy covariance tower GPP data, and the effectiveness of the merging method has been discussed. Finally, the variation trends of GPP over China under future climate change scenarios have been analyzed based on CMIP6 DL GPP.

Study Area
This study takes China as the research region. The climate characteristics of China are complex and various, with a monsoon climate in the east, a temperate continental climate in the northwest, and an alpine climate in the Qinghai-Tibet Plateau. In general, the climate in China is characterized by a monsoon climate with a high temperature, a rainy summer and cold winter with little rain. The high-temperature period is consistent with the rainy period in the study region. There is apparent spatial heterogeneity in GPP due to the differences in topography, vegetation types, and hydrothermal conditions in the different areas [54]. The growth of vegetation is closely related to the climate of the region where it is located. Therefore, according to the Coben climate division and the classic climate division of China, China is divided into four climatic parts for regional analysis: the arid, semi-arid, semi-humid, humid, and Qinghai-Tibet Plateau regions.

GLASS GPP Data
Global Land Surface Satellite (GLASS) gross primary production (GPP) is a global surface remote sensing product with a long time series and high precision based on multisource remote sensing data and ground observation [45]. The Bayesian multi-algorithm integration method integrates eight widely used light use efficiency (LUE) models, including CASA, CFix, CFlux, EC-LUE, MODIS, VPM, VPRM, and two-leaf. The spatial and temporal resolution of GLASS GPP are 5 km and eight days, from 1982 to 2018. The unit of the data is gCm −2 d −1 . The development and validation of the algorithm are based on the data from global fluxnet, which contains nine types of terrestrial ecosystems, such as evergreen broad-leaved forest, evergreen coniferous forest, deciduous broad-leaved forest, mixed forest, temperate grassland, tropical savanna, shrub, farmland, and tundra. GLASS GPP can be downloaded from http://www.geodata.cn (accessed on 23 April 2023).

CMIP6 GPP Data
The monthly GPP data from 23 available CMIP6 models are used in the study, with the unit of kgCm −2 s −1 . The period is over 1850-2100, including the historical and future projection period. More specifically, the projection period of 2015-2100 contains the four most interesting scenarios: SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5. The model and GLASS GPP data are first uniformly interpolated to 2 • due to the inconsistent spatial resolution. Meanwhile, it is necessary to unify the unit of model simulation data and GLASS GPP to gCm −2 d −1 . The first member, r1i1p1f1, is selected if possible (r1i1p1f2 is used if unavailable). CMIP6 GPP can be downloaded from https://esgf-node.llnl.gov/projects/cmip6/ (accessed on 9 June 2023). The information on the 23 models used is listed in the Table 1 below.

Deep Learning Network
In this study, a GPP-merged product conforming to the distribution characteristics of the GLASS GPP dataset is developed by merging the simulated data of CMIP6 models. A neural network named Ynet [51] is adopted in the data merging part, which combines data merging and image super-resolution technology. The deep learning model establishes relationships between the reference and the CMIP6 GPP data in the historical period and is then applied in future projections. Compared with the multi-model ensemble mean, the merged data have higher accuracy and spatial resolution over 1850-2100 in China. The multi-model ensemble mean is a simple average of the simulation results of the 23 CMIP6 models. It uses the simplest strategy of assuming that all 23 model data have the same uncertainty; thus, the multi-model ensemble mean is a simple average of each model. The merged data in this study are generated based on the deep learning method, which compensates for this deficiency. The novelty of the deep learning method used in this study is that it combines data merging and image super-resolution technology, which combines the advantages of 23 CMIP6 GPP datasets without relying on any prior knowledge.
The model architecture is shown in Figure 1. This structure consists of three parts. The first part is a symmetric encoder-decoder structure with a skip connection. Similar to the residual encoder-decoder network, this structure solves the problem of gradient disappearance when the error calculated by the loss function in the model training process is propagated back. This structure captures the abstract features of low-resolution images with noise and outputs a "cleaner" image, which can be regarded as a feature extractor. In addition to the 30 convolutional layers, this part also contains 15 deconvolution layers. It may introduce "checkerboard artefacts", resulting in lower image output quality [55]. Therefore, to reduce the problem, a convolutional layer is added after each deconvolution layer to compensate for upsampling. Specifically, this part's input to the model is lowresolution CMIP6 GPP. The second part of the model is upsampling, which is mainly used for downscaling. It consists of one upper sampling layer using bilinear interpolation and two convolution layers with the same feature depth as the input channel. As in the first part, the last two convolution layers were added to eliminate the checkerboard effect.
The third part is data merging, which consists of three input datasets, including upsampling output, auxiliary data, and unsampled CMIP6 GPP. The three datasets calculated by the two convolution layers are joined to obtain high-resolution data. Among them, auxiliary data are used to help improve the results, which remained constant throughout different months over the training and testing periods. The loss function can be calculated according to the following formula: Among them, θ is the network parameter that needs to be optimized. f (X i , θ) represents the learned function. In addition, n is the total number of training samples. X and Y represent the input data and the target at position i.
The model architecture is shown in Figure 1. This structure consists of three parts. The first part is a symmetric encoder-decoder structure with a skip connection. Similar to the residual encoder-decoder network, this structure solves the problem of gradient disappearance when the error calculated by the loss function in the model training process is propagated back. This structure captures the abstract features of low-resolution images with noise and outputs a "cleaner" image, which can be regarded as a feature extractor. In addition to the 30 convolutional layers, this part also contains 15 deconvolution layers. It may introduce "checkerboard artefacts", resulting in lower image output quality [55]. Therefore, to reduce the problem, a convolutional layer is added after each deconvolution layer to compensate for upsampling. Specifically, this part's input to the model is lowresolution CMIP6 GPP. The second part of the model is upsampling, which is mainly used for downscaling. It consists of one upper sampling layer using bilinear interpolation and two convolution layers with the same feature depth as the input channel. As in the first part, the last two convolution layers were added to eliminate the checkerboard effect.
The third part is data merging, which consists of three input datasets, including upsampling output, auxiliary data, and unsampled CMIP6 GPP. The three datasets calculated by the two convolution layers are joined to obtain high-resolution data. Among them, auxiliary data are used to help improve the results, which remained constant throughout different months over the training and testing periods. The loss function can be calculated according to the following formula: Among them, θ is the network parameter that needs to be optimized. f(Xi, θ) represents the learned function. In addition, n is the total number of training samples. X and Y represent the input data and the target at position i.

Evaluation Indicators
The evaluation indicators, including mean absolute error (MAE), root mean square deviation (RMSD), unbiased root mean square deviation (ubRMSD), and coefficient of determination (r 2 ), have been used to assess the performance of the merged data objectively and find out if it captures the distribution characteristics of the reference data. The calculation formulae are as follows:

Evaluation Indicators
The evaluation indicators, including mean absolute error (MAE), root mean square deviation (RMSD), unbiased root mean square deviation (ubRMSD), and coefficient of determination (r 2 ), have been used to assess the performance of the merged data objectively and find out if it captures the distribution characteristics of the reference data. The calculation formulae are as follows: where n represents the total number of samples. M i and R i represent the merged and reference data at i, respectively. M and R represent the average of the merged and reference data.

Deep Learning Merging Model Training
Statistical analysis ( Figure S1) shows that GPP data belong to heavy-tail distribution. The characteristic of the skewed distribution data is that there are more samples with little influence and fewer samples with substantial impact. The concrete manifestation in GPP data is that the number of pixels near the value of 0 is far more than that of other values. In addition, the larger the value, the smaller the number of pixels. Such datasets will make the deep learning network perform well in the low-value range, while the efficiency is relatively low in the high-value range, resulting in overall reduced accuracy. Therefore, the log1p function converts GPP data to make it more obedient to Gaussian distribution during data preprocessing.
The period from 1982 to 2014, when CMIP6 GPP data and GLASS GPP data coincide, was selected to establish the relationship between input data and reference data using the deep learning network. The GPP data simulated by 23 CMIP6 models and the reference data in the same month were matched as a group, totaling 396 samples. All samples were divided into three datasets according to time: training (1982-2010), validation (2011-2012), and test dataset (2013-2014). The number of iterations and the initial learning rate were 150 times and 1 × 10 −4 during the training stage, respectively. Figure 2 shows loss function curves for training and validation datasets. As shown in the figure, the loss reduced sharply in the first ten iterations of the training process. Subsequently, it shows a fluctuant reduction. The model is stable after 77 iterations, with a decreased loss. Loss no longer changes after 112 iterations, indicating that the model has been trained.
where n represents the total number of samples. Mi and Ri represent the merged and reference data at i, respectively. and represent the average of the merged and reference data.

Deep Learning Merging Model Training
Statistical analysis ( Figure S1) shows that GPP data belong to heavy-tail distribution. The characteristic of the skewed distribution data is that there are more samples with li le influence and fewer samples with substantial impact. The concrete manifestation in GPP data is that the number of pixels near the value of 0 is far more than that of other values. In addition, the larger the value, the smaller the number of pixels. Such datasets will make the deep learning network perform well in the low-value range, while the efficiency is relatively low in the high-value range, resulting in overall reduced accuracy. Therefore, the log1p function converts GPP data to make it more obedient to Gaussian distribution during data preprocessing.
The period from 1982 to 2014, when CMIP6 GPP data and GLASS GPP data coincide, was selected to establish the relationship between input data and reference data using the deep learning network. The GPP data simulated by 23 CMIP6 models and the reference data in the same month were matched as a group, totaling 396 samples. All samples were divided into three datasets according to time: training (1982-2010), validation (2011-2012), and test dataset (2013-2014). The number of iterations and the initial learning rate were 150 times and 1 × 10 −4 during the training stage, respectively. Figure 2 shows loss function curves for training and validation datasets. As shown in the figure, the loss reduced sharply in the first ten iterations of the training process. Subsequently, it shows a fluctuant reduction. The model is stable after 77 iterations, with a decreased loss. Loss no longer changes after 112 iterations, indicating that the model has been trained.

Quality Evaluation of Merging Results
The trained merging model has been tested. The results of comparing GLASS GPP, CMIP6 DL GPP, and multi-model ensemble mean with fluxnet data in China are shown in Figure 3. However, the results only represent the accuracy of the ten sites, as there are only ten unevenly distributed flux towers in China. As can be seen from the figure, all three datasets have a wide distribution range in the four indicators. The distribution range in MAE, RMSD and r 2 of the multi-model ensemble mean is wider than that of the other two datasets, indicating that it has a higher uncertainty. Although the distribution range is relatively small on ubRMSD, many outliers exist in the multi-model ensemble mean. The CMIP6DL GPP, and multi-model ensemble mean with fluxnet data in China are shown in Figure 3. However, the results only represent the accuracy of the ten sites, as there are only ten unevenly distributed flux towers in China. As can be seen from the figure, all three datasets have a wide distribution range in the four indicators. The distribution range in MAE, RMSD and r 2 of the multi-model ensemble mean is wider than that of the other two datasets, indicating that it has a higher uncertainty. Although the distribution range is relatively small on ubRMSD, many outliers exist in the multi-model ensemble mean. The uncertainty of CMIP6DL GPP is lower compared with the multi-model ensemble mean, while the median performs the best in MAE, RMSD, and ubRMSD. The merged data have improved overall (Figure 3a-c), indicating that data are closer to the actual value than the multi-model ensemble mean.  Table 2 shows the model's performance on the test dataset taking GLASS GPP as the benchmark, with the results showing that the merged GPP has improved in the aspect of the four indices and has a higher accuracy in China as a whole compared with the multimodel ensemble mean. Table 2. Mean absolute error (MAE), root mean square deviation (RMSD), unbiased root mean square difference (ubRMSD), and coefficient of determination (r 2 ) between merged GPP, the multimodel ensemble mean, and GLASS GPP. The normalized Taylor diagram obtained by calculating correlation coefficients, ubRMSD, and standard deviation ratios is used to evaluate further the CMIP6DL GPP, CMIP6 GPP, and muti-model ensemble mean. GLASS GPP has been used as the reference due to its relatively low uncertainties. Figure 4 shows the distribution of correlation coefficient, ubRMSD, and standard deviation ratio between 25 datasets (including CMIP6DL GPP, the multi-model ensemble mean, and 23 CMIP6 datasets) and the reference data in 2013-2014. As shown in Figure 4, the red dot represents GLASS GPP. Since it is used as reference data, the correlation coefficient and standard deviation ratio are equal to 1, and  Table 2 shows the model's performance on the test dataset taking GLASS GPP as the benchmark, with the results showing that the merged GPP has improved in the aspect of the four indices and has a higher accuracy in China as a whole compared with the multi-model ensemble mean. Table 2. Mean absolute error (MAE), root mean square deviation (RMSD), unbiased root mean square difference (ubRMSD), and coefficient of determination (r 2 ) between merged GPP, the multi-model ensemble mean, and GLASS GPP. The normalized Taylor diagram obtained by calculating correlation coefficients, ubRMSD, and standard deviation ratios is used to evaluate further the CMIP6 DL GPP, CMIP6 GPP, and muti-model ensemble mean. GLASS GPP has been used as the reference due to its relatively low uncertainties. Figure 4 shows the distribution of correlation coefficient, ubRMSD, and standard deviation ratio between 25 datasets (including CMIP6 DL GPP, the multi-model ensemble mean, and 23 CMIP6 datasets) and the reference data in 2013-2014. As shown in Figure 4, the red dot represents GLASS GPP. Since it is used as reference data, the correlation coefficient and standard deviation ratio are equal to 1, and ubRMSD is 0. The black dot represents CMIP6 DL GPP data, which is the closest to the red dot in the figure, indicating that the overall performance of CMIP6 DL GPP data is better than that of all datasets simulated by other models. The correlation coefficient of all CMIP6 GPP is less than 0.82, and ubRMSD is greater than 0.5. There are significant differences among the 23 model data, of which CAS-ESM2-0 performs the best, even better than the multi-model ensemble mean. The multi-model ensemble mean is still superior to most single CMIP6 data, though it is affected by the high error of model data. In general, CMIP6 DL GPP data perform the best in all aspects. dot in the figure, indicating that the overall performance of CMIP6DL GPP data is be er than that of all datasets simulated by other models. The correlation coefficient of all CMIP6 GPP is less than 0.82, and ubRMSD is greater than 0.5. There are significant differences among the 23 model data, of which CAS-ESM2-0 performs the best, even be er than the multi-model ensemble mean. The multi-model ensemble mean is still superior to most single CMIP6 data, though it is affected by the high error of model data. In general, CMIP6DL GPP data perform the best in all aspects. To evaluate the spatial accuracy of the merged data, the spatial absolute error distributions of CMIP6DL GPP, the multi-model ensemble mean, and GLASS GPP are plo ed, respectively. Figure 5a shows the MAE between CMIP6DL GPP and GLASS GPP in the test dataset, with the results showing that the regions with MAE less than 0.5 account for 73.91%. The errors in Northwest, Northeast, and North China are generally minor, less than 0.21 gCm −2 d −1 . However, the effects of CMIP6DL GPP in Central China and the eastern part of Southwest China are not ideal, with the error ranging from 0.5 to 1.5 gCm −2 d −1 , with a total of 7.4% of the regional errors being even more significant than 1 gCm −2 d −1 . Relatively large errors mainly occur in Yunnan Province, with two-thirds of the regional errors being more significant than 0.5 gCm −2 d −1 . Figure 5b shows MAE between the multimodel ensemble mean and GLASS GPP. Except for Northwest China, the error in most regions is beyond 1 gCm −2 d −1 . Specifically, the regions with an error greater than 1 gCm −2 d −1 account for 45.17%, even more significant than 2 gCm −2 d −1 in Central, East, and South China. In summary, CMIP6DL GPP generated by merging the model Ynet is highly reliable. To evaluate the spatial accuracy of the merged data, the spatial absolute error distributions of CMIP6 DL GPP, the multi-model ensemble mean, and GLASS GPP are plotted, respectively. Figure 5a shows the MAE between CMIP6 DL GPP and GLASS GPP in the test dataset, with the results showing that the regions with MAE less than 0.5 account for 73.91%. The errors in Northwest, Northeast, and North China are generally minor, less than 0.21 gCm −2 d −1 . However, the effects of CMIP6 DL GPP in Central China and the eastern part of Southwest China are not ideal, with the error ranging from 0.5 to 1.5 gCm −2 d −1 , with a total of 7.4% of the regional errors being even more significant than 1 gCm −2 d −1 . Relatively large errors mainly occur in Yunnan Province, with two-thirds of the regional errors being more significant than 0.5 gCm −2 d −1 . Figure 5b shows MAE between the multi-model ensemble mean and GLASS GPP. Except for Northwest China, the error in most regions is beyond 1 gCm −2 d −1 . Specifically, the regions with an error greater than 1 gCm −2 d −1 account for 45.17%, even more significant than 2 gCm −2 d −1 in Central, East, and South China. In summary, CMIP6 DL GPP generated by merging the model Ynet is highly reliable.

Assessment of the Variations in the Merged GPP
GLASS GPP is a dataset generated by the light use efficiency model integrating CO2, radiation, and vapor pressure deficit (VPD), providing a reliable long-term GPP estimation. Figure 6 shows the long-term changes in GLASS GPP, CMIP6DL GPP, and multimodel ensemble mean from 1982 to 2014. As shown in Figure 6a, the variations in GLASS GPP are not apparent and almost remain unchanged in the northwest arid regions. There

Assessment of the Variations in the Merged GPP
GLASS GPP is a dataset generated by the light use efficiency model integrating CO 2 , radiation, and vapor pressure deficit (VPD), providing a reliable long-term GPP estimation. Figure 6 shows the long-term changes in GLASS GPP, CMIP6 DL GPP, and multi-model ensemble mean from 1982 to 2014. As shown in Figure 6a, the variations in GLASS GPP are not apparent and almost remain unchanged in the northwest arid regions. There is an increasing trend in Central China and a significant decrease in the eastern part of Inner Mongolia, the southern part of East China, and South China. The long-term changes in the multi-model ensemble mean are shown in Figure 6c. There is a sharp increasing trend throughout China, except for the western regions. The characteristics of the decrease in the southern coastal areas and Inner Mongolia are not captured, which is inconsistent with the spatial characteristics of the variations in GPP over China. Figure 6b describes the spatial distribution of the long-term variations of the merged GPP data, similar to GLASS GPP. However, the variations are gentler than the reference data since the deep learning merging model is built using a convolution neural network with a "smoothing" effect. That is, it is impossible for the deep learning merging model to learn the extreme value. Compared with the multi-model ensemble mean, the data quality has been significantly improved in eastern Inner Mongolia and South China, dramatically increasing the merged dataset's reliability in the projection period.

Assessment of the Variations in the Merged GPP
GLASS GPP is a dataset generated by the light use efficiency model integrating CO2, radiation, and vapor pressure deficit (VPD), providing a reliable long-term GPP estimation. Figure 6 shows the long-term changes in GLASS GPP, CMIP6DL GPP, and multimodel ensemble mean from 1982 to 2014. As shown in Figure 6a, the variations in GLASS GPP are not apparent and almost remain unchanged in the northwest arid regions. There is an increasing trend in Central China and a significant decrease in the eastern part of Inner Mongolia, the southern part of East China, and South China. The long-term changes in the multi-model ensemble mean are shown in Figure 6c. There is a sharp increasing trend throughout China, except for the western regions. The characteristics of the decrease in the southern coastal areas and Inner Mongolia are not captured, which is inconsistent with the spatial characteristics of the variations in GPP over China. Figure 6b describes the spatial distribution of the long-term variations of the merged GPP data, similar to GLASS GPP. However, the variations are gentler than the reference data since the deep learning merging model is built using a convolution neural network with a "smoothing" effect. That is, it is impossible for the deep learning merging model to learn the extreme value. Compared with the multi-model ensemble mean, the data quality has been significantly improved in eastern Inner Mongolia and South China, dramatically increasing the merged dataset's reliability in the projection period.  To evaluate the differences in the merged results in the time dimension, the grids are randomly selected in four climatic regions over China (arid, semi-arid and semi-humid, humid region, and Qinghai-Tibet Plateau) to test the time series. As can be seen from Figure 7, there are similar seasonal fluctuations in the three datasets. The fluctuation range of the test dataset is highly consistent in the CMIP6 DL GPP and the reference data, while it is more significant in the multi-model ensemble mean. GPP has consistently been overestimated throughout the period in the multi-model ensemble mean. As the GPP increases, the error of the multi-model ensemble mean increases. Specifically, the high value of GPP is concentrated in June and July, when the error of the multi-model ensemble mean is the largest. The impacts of temperature, water, solar radiation, and vegetation growth are the main causes of seasonal variation. The temperature, precipitation, vegetation leaf area, and solar radiation of each region reach the highest values during in the year in summer, and then gradually decrease. The comprehensive impact of these factors ultimately leads to an inverted U-shaped change in the average monthly GPP of each vegetation coverage area.
is concentrated in June and July, when the error of the multi-model ensemble mean is the largest. The impacts of temperature, water, solar radiation, and vegetation growth are the main causes of seasonal variation. The temperature, precipitation, vegetation leaf area, and solar radiation of each region reach the highest values during in the year in summer, and then gradually decrease. The comprehensive impact of these factors ultimately leads to an inverted U-shaped change in the average monthly GPP of each vegetation coverage area.

Temporal and Spatial Distribution of Merged GPP in Historical Period
The GPP data over the historical period (1850-2014) and the four climate scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5) are merged using the trained deep learning merging model based on the data simulated by 23 CMIP6 models. The spatial distribution of GPP in different seasons is shown in Figure 8 to study the seasonal variation characteristics during the historical period (1850-2014). There is a significant spatial distribution gradient of GPP in various seasons in China, mainly due to land use and cover conditions and temperature and water environmental factors. In southern China, the climate is warm and humid, with the best hydrothermal conditions, sufficient sunlight and high light utilization efficiency, high vegetation coverage, and the highest GPP in the growing season. GPP is high in northeast forest regions in summer. However, photosynthesis is inhibited by low temperature in other seasons, resulting in relatively low GPP. In the deserts of the northwest, most areas of the Qinghai-Tibet Plateau, and the grasslands in central and western Inner Mongolia, the plant growing season is relatively short due to water stress and low-temperature constraints, which leads to low productivity and even the lowest GPP. Overall, GPP decreases from southeast to northwest and coastal to inland, consistent with current research results [56][57][58][59]. The spatial distribution of GPP is highly compatible with precipitation in summer, indicating that precipitation is an important factor affecting the distribution of GPP in China. In addition, the GPP near the Tianshan Mountains in Xinjiang is relatively high in summer. The prevailing westerly winds from the Atlantic Ocean and the airflow from the Arctic Ocean encounter the uplift of mountain slopes and

Temporal and Spatial Distribution of Merged GPP in Historical Period
The GPP data over the historical period (1850-2014) and the four climate scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5) are merged using the trained deep learning merging model based on the data simulated by 23 CMIP6 models. The spatial distribution of GPP in different seasons is shown in Figure 8 to study the seasonal variation characteristics during the historical period (1850-2014). There is a significant spatial distribution gradient of GPP in various seasons in China, mainly due to land use and cover conditions and temperature and water environmental factors. In southern China, the climate is warm and humid, with the best hydrothermal conditions, sufficient sunlight and high light utilization efficiency, high vegetation coverage, and the highest GPP in the growing season. GPP is high in northeast forest regions in summer. However, photosynthesis is inhibited by low temperature in other seasons, resulting in relatively low GPP. In the deserts of the northwest, most areas of the Qinghai-Tibet Plateau, and the grasslands in central and western Inner Mongolia, the plant growing season is relatively short due to water stress and low-temperature constraints, which leads to low productivity and even the lowest GPP. Overall, GPP decreases from southeast to northwest and coastal to inland, consistent with current research results [56][57][58][59]. The spatial distribution of GPP is highly compatible with precipitation in summer, indicating that precipitation is an important factor affecting the distribution of GPP in China. In addition, the GPP near the Tianshan Mountains in Xinjiang is relatively high in summer. The prevailing westerly winds from the Atlantic Ocean and the airflow from the Arctic Ocean encounter the uplift of mountain slopes and generate intense precipitation, which promotes a high vegetation coverage in the west and northwest of the Tianshan Mountains [60]. At the same time, the high mountain snowmelt in the Tianshan Mountains also plays a promoting role.
The spatial distribution of GPP during the four seasons is significantly different, mainly due to the comprehensive influences of the monsoon climate and vegetation type. Throughout the whole year, GPP is the highest in summer (June to August), followed by spring (March to May), autumn (September to November), and the lowest in winter, which entirely fits within the ecological definition of the GPP itself. In spring, the high value of GPP is concentrated in the area south of 28 • N. As the summer monsoon moves northward, the temperature in most parts increases, with the average GPP in the historical period reaching 13.6 gCm −2 d −1 . The spatial distribution of GPP is highly affected by hydrothermal conditions. Thus, the high value is mainly distributed in Yunnan and Guizhou. In autumn, the spatial distribution is similar to that in spring due to the decrease in temperature in northern China. GPP is mostly lower than 1 gCm −2 d −1 in winter, and even 0 in some regions. Due to the minimum precipitation and temperature in winter, most of the vegetation is in the non-growing season, which makes the vegetation coverage the lowest, resulting in the lowest GPP throughout the year. Although crops are sown in the north, they grow slowly due to the low temperature. On the contrary, the GPP in southern Yunnan remains high even in winter (Figure 8d). mainly due to the comprehensive influences of the monsoon climate and vegetation type. Throughout the whole year, GPP is the highest in summer (June to August), followed by spring (March to May), autumn (September to November), and the lowest in winter, which entirely fits within the ecological definition of the GPP itself. In spring, the high value of GPP is concentrated in the area south of 28° N. As the summer monsoon moves northward, the temperature in most parts increases, with the average GPP in the historical period reaching 13.6 gCm −2 d −1 . The spatial distribution of GPP is highly affected by hydrothermal conditions. Thus, the high value is mainly distributed in Yunnan and Guizhou. In autumn, the spatial distribution is similar to that in spring due to the decrease in temperature in northern China. GPP is mostly lower than 1 gCm −2 d −1 in winter, and even 0 in some regions. Due to the minimum precipitation and temperature in winter, most of the vegetation is in the non-growing season, which makes the vegetation coverage the lowest, resulting in the lowest GPP throughout the year. Although crops are sown in the north, they grow slowly due to the low temperature. On the contrary, the GPP in southern Yunnan remains high even in winter (Figure 8d). The global climate has undergone significant changes in the past few decades, including frequent extreme climate disasters and sustained warming. Ecosystems and terrestrial carbon sinks will be significantly affected by climate change. On the other hand, GPP is greatly influenced by human activity [61]. Although large-scale afforestation benefits vegetation restoration and increases GPP, population growth and rapid economic development have exacerbated urbanization and disrupted the balance of terrestrial ecosystems. These human interventions significantly impact the formation of GPP dynamics [62,63]. Figure 9 shows the spatial distribution of variation trends in GPP in all seasons over China during the historical period (1850-2014). There are increasing trends in most regions; however, changes in these trends exhibit significant heterogeneity in space and time. The trend is the strongest in summer, followed by spring and autumn, and the most gentle in winter. In summer, the variation trend of GPP in more than one-third of the The global climate has undergone significant changes in the past few decades, including frequent extreme climate disasters and sustained warming. Ecosystems and terrestrial carbon sinks will be significantly affected by climate change. On the other hand, GPP is greatly influenced by human activity [61]. Although large-scale afforestation benefits vegetation restoration and increases GPP, population growth and rapid economic development have exacerbated urbanization and disrupted the balance of terrestrial ecosystems. These human interventions significantly impact the formation of GPP dynamics [62,63]. Figure 9 shows the spatial distribution of variation trends in GPP in all seasons over China during the historical period (1850-2014). There are increasing trends in most regions; however, changes in these trends exhibit significant heterogeneity in space and time. The trend is the strongest in summer, followed by spring and autumn, and the most gentle in winter. In summer, the variation trend of GPP in more than one-third of the regions is relatively large, concentrated in Northeast China, the south of North China, and the southwest regions. The annual mean value of GPP in these regions is significant, with the variation trends also being substantial. It is mainly related to the superior local water and heat conditions and the suitable climate environment, which results in the overall growth of vegetation. The GPP in the south of Central China and East China shows a weak decreasing trend. It is mainly due to the increase in temperature in the southern region since the 1990s, which increases the frequency of heat waves and drought events [64,65]. There is a negative impact on the terrestrial carbon cycle, resulting in a decline in GPP. Compared with summer, the regions with a substantial increase in autumn tend to reduce; meanwhile, the trend is gentler. In spring, the increasing trends concentrate in Anhui and Guizhou. In winter, there is a trend of widespread decline in Inner Mongolia, Jilin, and Liaoning, accounting for 31.57% of the whole region, while there is an increasing trend in some regions in the south. since the 1990s, which increases the frequency of heat waves and drought events [64,65]. There is a negative impact on the terrestrial carbon cycle, resulting in a decline in GPP. Compared with summer, the regions with a substantial increase in autumn tend to reduce; meanwhile, the trend is gentler. In spring, the increasing trends concentrate in Anhui and Guizhou. In winter, there is a trend of widespread decline in Inner Mongolia, Jilin, and Liaoning, accounting for 31.57% of the whole region, while there is an increasing trend in some regions in the south.

Long-Term Changes in Merged GPP in Projection Period under Multiple Climate Scenarios
To explore the difference between the changes in GPP in the near future (2021-2040) and the far future (2071-2100), the spatial distribution of the multi-year mean of GPP in each season under different climate scenarios in historical baseline, near, and far future is described in Figures 10 and S2, respectively. Among them, the distributions under different scenarios in the future describe the change in the multi-year mean compared with that during baseline. The positive value represents the increase in GPP, while the negative one represents the decrease. Figure 10a1-a4 shows that the spatial distribution of the mean annual GPP in each season of the baseline GPP increased significantly in the four seasons compared with the seasonal average over the period 1850-2014 ( Figure 8). In spring, the average annual GPP greater than 2 gCm −2 d −1 occupied 22.93% of the region, concentrated in the south of 30° N. In summer, the area with an annual average GPP of more than 2 gCm −2 d −1 doubled compared with spring, reaching 48.97%, with the maximum yearly peak at 14.19 gCm −2 d −1 . There are similar peaks and spatial distribution in autumn and spring, while it remains almost the same in winter.
There are significant seasonal differences in the variations in GPP based on the historical baseline under multiple scenarios in the near future. In spring, GPP increases in

Long-Term Changes in Merged GPP in Projection Period under Multiple Climate Scenarios
To explore the difference between the changes in GPP in the near future (2021-2040) and the far future (2071-2100), the spatial distribution of the multi-year mean of GPP in each season under different climate scenarios in historical baseline, near, and far future is described in Figure 10 and Figure S2, respectively. Among them, the distributions under different scenarios in the future describe the change in the multi-year mean compared with that during baseline. The positive value represents the increase in GPP, while the negative one represents the decrease. Figure 10a1-a4 shows that the spatial distribution of the mean annual GPP in each season of the baseline GPP increased significantly in the four seasons compared with the seasonal average over the period 1850-2014 ( Figure 8). In spring, the average annual GPP greater than 2 gCm −2 d −1 occupied 22.93% of the region, concentrated in the south of 30 • N. In summer, the area with an annual average GPP of more than 2 gCm −2 d −1 doubled compared with spring, reaching 48.97%, with the maximum yearly peak at 14.19 gCm −2 d −1 . There are similar peaks and spatial distribution in autumn and spring, while it remains almost the same in winter.
There are significant seasonal differences in the variations in GPP based on the historical baseline under multiple scenarios in the near future. In spring, GPP increases in most eastern regions, and rises more with the enhancement of radiation forcing levels. However, GPP decreased significantly in the east of the plain of Henan Province under the high emission scenario. In summer, GPP increased substantially in nearly half of the regions of the country, including the Northeast, North, Southwest, and south of Northwest China. GPP accelerated with the increase in radiation forcing levels, with more than 20% of the regions having values greater than 0.4 gCm −2 d −1 under SSP5-8.5. In autumn, it increased significantly in Yunnan, while there was a slight weakening trend in the north of Henan. The spatial characteristics distribution maintains consistency under multiple climate scenarios. Winter is the only season when GPP decreases in a large area. Under the four scenarios, the GPP in 21% of regions is lower than the baseline, mainly distributed in the west, Inner Mongolia, and Heilongjiang. The reduction in GPP is the same under different scenarios. Except for the decrease in GPP in Fujian, it generally increases in the south.
increased significantly in Yunnan, while there was a slight weakening trend in the north of Henan. The spatial characteristics distribution maintains consistency under multiple climate scenarios. Winter is the only season when GPP decreases in a large area. Under the four scenarios, the GPP in 21% of regions is lower than the baseline, mainly distributed in the west, Inner Mongolia, and Heilongjiang. The reduction in GPP is the same under different scenarios. Except for the decrease in GPP in Fujian, it generally increases in the south. The change in GPP becomes increasingly intense with the enhancement of radiative forcing levels in the far future, as seen from Figure S2. The variations in GPP are studied with the Heihe-Tengchong line as the boundary. In spring, it can be found that GPP increases more as radiative forcing levels are enhanced in the eastern part of the dividing line. The area with an increase of more than 0.8 gCm −2 d −1 rises from 0.62% to 23.99%. The change in the west of the boundary is gentle, with an average rise of 0.09 gCm −2 d −1 . The area of GPP increasing by more than 0.8 gCm −2 d −1 expands from 0.62% to 23.99% of the regions. The change in the western part of the dividing line is relatively gentle, with an average increase of 0.09 gCm −2 d −1 . Summer is the peak photosynthesis season for plants in northern ecosystems compared to other seasons, leading to a more extensive range and broader areas of increase. GPP increases in about 94% of the regions under SSP1-2.6 relative to the baseline. With the enhancement of radiation forcing levels, the rise of GPP in Northeast, North, and the east of Southwest China becomes more and more intense, while that in Henan is turning from decrease to increase. The magnitude of the increase is further accelerated under SSP5-8.5. The reduction in GPP in Henan in autumn only occurs in the low-emission scenario, which differs from spring. In winter, GPP reduces in parts of the west and northeast regions, roughly similar to the scenario in the near future (Figure The change in GPP becomes increasingly intense with the enhancement of radiative forcing levels in the far future, as seen from Figure S2. The variations in GPP are studied with the Heihe-Tengchong line as the boundary. In spring, it can be found that GPP increases more as radiative forcing levels are enhanced in the eastern part of the dividing line. The area with an increase of more than 0.8 gCm −2 d −1 rises from 0.62% to 23.99%. The change in the west of the boundary is gentle, with an average rise of 0.09 gCm −2 d −1 . The area of GPP increasing by more than 0.8 gCm −2 d −1 expands from 0.62% to 23.99% of the regions. The change in the western part of the dividing line is relatively gentle, with an average increase of 0.09 gCm −2 d −1 . Summer is the peak photosynthesis season for plants in northern ecosystems compared to other seasons, leading to a more extensive range and broader areas of increase. GPP increases in about 94% of the regions under SSP1-2.6 relative to the baseline. With the enhancement of radiation forcing levels, the rise of GPP in Northeast, North, and the east of Southwest China becomes more and more intense, while that in Henan is turning from decrease to increase. The magnitude of the increase is further accelerated under SSP5-8.5. The reduction in GPP in Henan in autumn only occurs in the low-emission scenario, which differs from spring. In winter, GPP reduces in parts of the west and northeast regions, roughly similar to the scenario in the near future (Figure 10b1-e4). There is a considerable area reduction in GPP in the eastern part of the Qinghai-Tibet Plateau under the SSP1-2.6 scenario. With the increase in radiation forcing levels, the area of GPP reduction gradually decreases. GPP increased south of 30 • N, with the maximum increase reaching 2.48 gCm −2 d −1 in Yunnan under the SSP5-8.5 scenario. Compared with the near future, GPP growth will be more dramatic in the far future.
The multi-year variations in GPP in each pixel can better reflect the characteristics of spatial changes in addition to the regional average. Figure 11 describes the spatial distribution of GPP variation trends in China under four scenarios from 2015 to 2100. GPP generally shows an increasing trend in the projection period, with the trend becoming more and more intense as the radiation forcing levels increase. There is no significant change in GPP in the northwest arid regions under different scenarios, which maintains a low growth. A slight decreasing trend can be observed visually at the junction of Shandong, Henan, Jiangsu, and Anhui provinces under the SSP1-2.6 scenario, which is related to the environmental effects of aerosols [66]. The spatial pattern of GPP variations under the medium and high emission scenario is similar, with GPP increasing in each pixel. Likewise, the variation trends become more intense as radiative forcing levels increase. The most significant increase in GPP occurs in Yunnan and Guizhou under SSP2-4.5, SSP3-7.0, and SSP5-8.5. The regions where GPP is more sensitive to radiation forcing levels are concentrated in the southwest, the east of the northwest, and the northeast. GPP increases rapidly as radiative forcing levels increase. It is found that the increase in the high-emission scenario is even four times that of the low-emission scenario.
spatial changes in addition to the regional average. Figure 11 describes the spatial distribution of GPP variation trends in China under four scenarios from 2015 to 2100. GPP generally shows an increasing trend in the projection period, with the trend becoming more and more intense as the radiation forcing levels increase. There is no significant change in GPP in the northwest arid regions under different scenarios, which maintains a low growth. A slight decreasing trend can be observed visually at the junction of Shandong, Henan, Jiangsu, and Anhui provinces under the SSP1-2.6 scenario, which is related to the environmental effects of aerosols [66]. The spatial pa ern of GPP variations under the medium and high emission scenario is similar, with GPP increasing in each pixel. Likewise, the variation trends become more intense as radiative forcing levels increase. The most significant increase in GPP occurs in Yunnan and Guizhou under SSP2-4.5, SSP3-7.0, and SSP5-8.5. The regions where GPP is more sensitive to radiation forcing levels are concentrated in the southwest, the east of the northwest, and the northeast. GPP increases rapidly as radiative forcing levels increase. It is found that the increase in the high-emission scenario is even four times that of the low-emission scenario. Considering the characteristics of the future variations in GPP, it is divided into three periods for analysis, including 2021-2040 (near future), 2041-2070 (medium future), and 2071-2100 (far future) due to the long projection period of the merged data. As can be seen from Figure S3, the increasing trend of GPP turns into a decreasing trend as time goes on under the SSP1-2.6 scenario. In the near future, GPP shows a significant increasing trend except for the junction of Shandong, Jiangsu, Henan, and Anhui provinces. In the medium future, the changes in average regional GPP are relatively stable, though there are scattered pixels with increasing and decreasing trends of GPP across the country. In the far future, the regions with reduced trends of GPP nationwide account for about 93% of the country. Under the low emission scenario, the air temperature rises first and remains stable in the far future, while the carbon dioxide emissions begin to decrease in 2020 [67]. This shows that GPP is greatly affected by air temperature. The amount of carbon dioxide Considering the characteristics of the future variations in GPP, it is divided into three periods for analysis, including 2021-2040 (near future), 2041-2070 (medium future), and 2071-2100 (far future) due to the long projection period of the merged data. As can be seen from Figure S3, the increasing trend of GPP turns into a decreasing trend as time goes on under the SSP1-2.6 scenario. In the near future, GPP shows a significant increasing trend except for the junction of Shandong, Jiangsu, Henan, and Anhui provinces. In the medium future, the changes in average regional GPP are relatively stable, though there are scattered pixels with increasing and decreasing trends of GPP across the country. In the far future, the regions with reduced trends of GPP nationwide account for about 93% of the country. Under the low emission scenario, the air temperature rises first and remains stable in the far future, while the carbon dioxide emissions begin to decrease in 2020 [67]. This shows that GPP is greatly affected by air temperature. The amount of carbon dioxide directly affects vegetation's carbon sequestration when the temperature is stable, leading to GPP decline. GPP increases in the near future under the SSP2-4.5 scenario; however, the variation trends slow down over time. Under SSP3-7.0 and SSP5-8.5 scenarios, GPP increased over time in all regions, with the trends strengthening gradually.

Discussion
It can be seen from Figure 10 and Figure S2 that the variations in GPP in the north and south of China are opposite in winter while increasing in all regions in spring, summer, and autumn. However, this cannot reflect the impact of the reduction in winter on the change in the average GPP across the country. Therefore, Figure 12 shows the inter-annual changes in GPP in China under the historical baseline and different scenarios in the future. The thick grey line represents the line between the baseline and the future. GPP shows an increasing trend from 1.76 gCm −2 d −1 to 1.84 gCm −2 d −1 from 1995 to 2014. It may result from a combination of factors such as appropriate climate change, the fertilization effect of CO2, nitrogen deposition, and human activities such as afforestation and agricultural irrigation [68][69][70][71]. From 2015 to 2035, GPP increases steadily yearly with similar increasing trends under the four scenarios. From 2036, the difference in the variations in GPP in the four scenarios is noticeable. Among them, the variation speed of GPP is the fastest under SSP5-8.5, with the average reaching 2.49 gCm −2 d −1 at the end of the 21st century. The second fastest occurs under SSP3-7.0, which is slightly slower than that of the highemission scenario. The annual average is roughly the same before 2055 as under SSP2-4.5. However, the gap between the two scenarios gradually increases from 2067, when the growth rate under SSP2-4.5 slows down. From 2081 to 2100, it is almost stable and begins declining slowly in 2097. Among the four scenarios, GPP decreases significantly only in the SSP1-2.6 scenario. Although there is a decreasing trend at the end of the century under the SSP2-4.5 scenario (Figure 12), GPP shows an overall increasing trend in each pixel over the period 2015-2100 (Figure 11b). More specifically, it keeps increasing, as do the other three scenarios in the near future, followed by a slight fluctuation around 1.96 gCm −2 d −1 for a long time in the medium term. It is worth noting that GPP declines significantly in the far future, especially from 2075 onwards. Meteorological factors are the dominant factors in the interannual variation in global GPP. In the context of climate warming, the frequent occurrence of extreme events such as high temperatures and droughts is expected to impact vegetation growth, leading to fluctuations in ecosystem GPP. According to the definition of severe climate, the Intergovernmental Panel on Climate Change [72] pointed out that climate change has led to changes in the frequency, intensity, spatial scope, duration, and time of extreme climate, which may have an unprecedented impact on the terrestrial carbon cycle. In addition, due to global warming, climate change is expected to further increase the frequency, persistence, and intensity of extreme weather in the mid to late 21st century [73][74][75], making the impact of future climate change on terrestrial ecosystems even more uncertain [76,77].
autumn. However, this cannot reflect the impact of the reduction in winter on the change in the average GPP across the country. Therefore, Figure 12 shows the inter-annual changes in GPP in China under the historical baseline and different scenarios in the future. The thick grey line represents the line between the baseline and the future. GPP shows an increasing trend from 1.76 gCm −2 d −1 to 1.84 gCm −2 d −1 from 1995 to 2014. It may result from a combination of factors such as appropriate climate change, the fertilization effect of CO2, nitrogen deposition, and human activities such as afforestation and agricultural irrigation [68][69][70][71]. From 2015 to 2035, GPP increases steadily yearly with similar increasing trends under the four scenarios. From 2036, the difference in the variations in GPP in the four scenarios is noticeable. Among them, the variation speed of GPP is the fastest under SSP5-8.5, with the average reaching 2.49 gCm −2 d −1 at the end of the 21st century. The second fastest occurs under SSP3-7.0, which is slightly slower than that of the high-emission scenario. The annual average is roughly the same before 2055 as under SSP2-4.5. However, the gap between the two scenarios gradually increases from 2067, when the growth rate under SSP2-4.5 slows down. From 2081 to 2100, it is almost stable and begins declining slowly in 2097. Among the four scenarios, GPP decreases significantly only in the SSP1-2.6 scenario. Although there is a decreasing trend at the end of the century under the SSP2-4.5 scenario (Figure 12), GPP shows an overall increasing trend in each pixel over the period 2015-2100 (Figure 11b). More specifically, it keeps increasing, as do the other three scenarios in the near future, followed by a slight fluctuation around 1.96 gCm −2 d −1 for a long time in the medium term. It is worth noting that GPP declines significantly in the far future, especially from 2075 onwards. Meteorological factors are the dominant factors in the interannual variation in global GPP. In the context of climate warming, the frequent occurrence of extreme events such as high temperatures and droughts is expected to impact vegetation growth, leading to fluctuations in ecosystem GPP. According to the definition of severe climate, the Intergovernmental Panel on Climate Change [72] pointed out that climate change has led to changes in the frequency, intensity, spatial scope, duration, and time of extreme climate, which may have an unprecedented impact on the terrestrial carbon cycle. In addition, due to global warming, climate change is expected to further increase the frequency, persistence, and intensity of extreme weather in the mid to late 21st century [73][74][75], making the impact of future climate change on terrestrial ecosystems even more uncertain [76,77].  Large-scale studies using merged GPP can reduce the uncertainty of terrestrial ecosystem carbon cycle research. However, the merging process will also introduce certain uncertainties, mainly affected by input errors, scale effects, and merging algorithms. Input errors are caused by the CMIP6 datasets, GLASS GPP, and flux tower observation GPP. Although model simulation provides an essential research tool for the carbon cycle of large-scale terrestrial ecosystems, there are still significant differences in the simulation results of the models at regional and global scales due to the differences in the structure, parameter settings, input data, and spatial resolution of each CMIP6 model. Schaefer et al. [5] compared GPP simulation differences in North America across 26 models to explore the responses of different model structures and environmental factors to simulations based on data from 39 flux stations. Zhang et al. [78] studied the impact of different parameter settings and meteorological data in the CEVSA2 model on the GPP simulation results in the Changbai Mountain region of China from 2003 to 2005, indicating that the GPP difference caused by parameter setting in the model is 5% to 8%. The reference data GLASS GPP is generated by integrating multiple algorithms using the Bayesian integration method, reducing the uncertainty of a single algorithm. However, there are differences in spatial scales between ground and satellite data, and the algorithm fails to distinguish potential differences in photosynthesis between C3 and C4 crops. In addition, the problem of mixed pixels can also affect the accuracy of GLASS GPP estimation [79]. The GPP observed by the flux tower is considered the benchmark, while the high heterogeneity of the ground can also make the representativeness of the flux station data poor. Therefore, significant uncertainties in the GPP benchmark lead to a lack of consensus on the global GPP distribution [15,80]. Chen et al. [81] showed that the flux-tower-estimated GPP represents 50% to 60% of the actual situation in the case of surface heterogeneity. There is still a problem of insufficient spatial representation compared to the range of flux contribution areas assumed by eddy covariance.
The mismatch between the spatial resolution of the model and the floor area of the flux tower causes the uncertainty imposed by scale effects. Even with the same model and source of input data, the simulation results may still be different due to spatial heterogeneity caused by different spatial resolutions [1]. At the regional scale, improving spatial resolution can reduce the uncertainty of GPP simulation [82]. However, in practical applications, models often need to simulate large-scale and long-time series data. When considering factors such as simulation time, data storage capacity, and analysis difficulty, they often reduce the spatial resolution to improve simulation efficiency. In GPP estimation at different spatial scales, there will be different judgments on the coverage type, leading to incorrect maximum light energy utilization, further resulting in errors. In mid-and low-resolution GPP remote sensing products, most of the pixels under heterogeneous ground surfaces are mixed. For example, there is a significant difference between MODIS surface coverage with 1 km and global 30 m resolution. At different scales, heterogeneity characteristics have an important impact on surface vegetation parameter characteristics [83]. Taking the mixed pixel of forest and grassland as an example, there is a significant difference in the a priori maximum light energy utilization rate between the two types. It may bring about significant errors using a single type of maximum weak energy utilization rate to replace pixel values. Wang et al. [1] discussed the impact of two spatial resolutions of 1 km and 8 km input data on GPP simulation, with the results showing that the difference was mainly due to the difference in LAI caused by mixed pixels within the range of 8km. The difference in simulation results with different spatial resolutions at forest stations was more significant than that at grassland stations. Therefore, it is feasible to appropriately reduce spatial resolution to improve the model's simulation efficiency. However, it is necessary to minimize the error of low spatial resolution in the GPP simulation of forest ecosystems and growing seasons.
In the long sequence data merging, considering time information can effectively enhance model learning capabilities and improve the quality of merged data, which helps to reduce the data uncertainty in the projection period. Since the deep learning merging model focuses more on spatial features and no individual modules can capture temporal characteristics, the temporal information is not fully utilized. In addition, there are other factors contributing to uncertainties. For example, complex vegetation structures, significant dry and wet seasons, and cloud cover can also bring uncertainties [84]. In highly productive regions, methods based on optical remote sensing often underestimate GPP due to the saturation of reflectance measurements in high-density canopy [85]. In addition, optical remote sensing is strongly affected by cloud cover, resulting in data gaps and high uncertainty in areas with frequent cloud cover and high GPP, such as tropical forests [16].

Conclusions
In this study, a GPP-merged dataset with high resolution in China has been developed using the deep learning method. It combines the advantages of 23 CMIP6 GPP datasets without relying on any prior knowledge. The lower resolution of the multi-model data has been improved to 0.25 • of the merged GPP. From the perspective of evaluation indicators, the performance of merged data is superior to the others in all aspects. In China, the error has been kept at a low level. The merged dataset has significantly improved in temporal and spatial dimensions compared with every single model and the multi-model mean. The merged data are closer to the actual value than the multi-model mean according to the fluxnet data. The merged data greatly reduce the uncertainties and improve the reliability of GPP estimation in the projection period. A further analysis of the spatiotemporal change pattern of GPP in China has resulted in the following findings.
During the baseline, GPP is the highest in summer, followed by spring and autumn, and the lowest in winter due to the combined effects of factors such as the monsoon climate and vegetation types. GPP exhibits significant spatial heterogeneity due to differences in topography, vegetation types, and hydrothermal conditions in different regions. Generally, GPP increases gradually from the northwest inland to the southeast coast. GPP in the southern part of Yunnan remains high in winter. According to the long-term variation trends, GPP is increasing in most regions of China. It is worth noting that the characteristics of the decrease in the southern coastal areas and Inner Mongolia are not captured by the multi-model ensemble mean; however, they are captured by the merged GPP. There are seasonal differences in the variations in GPP relative to the baseline under different scenarios in the future projection period, while the variations are almost the same across different scenarios in the near future. In summer, GPP increased rapidly with the increase in radiation forcing levels. In winter, GPP decreases in a large region mainly distributed in the west, Inner Mongolia, and Heilongjiang. The long-term variation trends of GPP in China under different scenarios are significantly different. The growth rate of GPP has reached the maximum under the SSP5-8.5 scenario since 2036. The increasing trend becomes slightly slower under the SSP3-7.0 scenario than under the high-emission scenario. Furthermore, GPP decreased significantly under the SSP1-2.6 scenario. The annual GPP in each pixel in China shows an increasing trend, which is more intense with the increase in radiation forcing levels. GPP changes from increasing to decreasing over time, since GPP is affected by temperature and carbon dioxide under the SSP1-2.6 scenario. Under the SSP2-4.5 scenario, GPP will increase in the near future; however, the trend will slow down as time goes by. Under the SSP3-7.0 and SSP5-8.5 scenarios, GPP increases with time; the increasing trend rises gradually.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.