Evaluation of annual trends in carbon cycle variables simulated by CMIP6 Earth system models in China

Systematic evaluation of the carbon cycle physical and biological variables simulated in Earth System Model (ESM) participating in the Coupled Model Intercomparison Project Phase 6 (CMIP 6) is fundamental to the understanding of terrestrial ecosystems, as well as to future projections. Leaf Area Index (LAI), Gross Primary Productivity (GPP), Net Primary Productivity (NPP), Net Ecosystem Productivity (NEP) and Land Surface Temperature (LST) as key indicators of carbon cycle performance in ESM outputs, play a critical role in evaluating ecosystem functions. Assessing these metrics can provide valuable insights into the biases in model-simulated ecosystems and offer guidance for model improvement. In this study, we assessed the interannual trends performance of LAI, GPP, NPP, NEP and LST simulated by 12 CMIP6 ESMs during the historical period by using satellite LAI, NPP, NEP, LST and CSIF data as observations. The findings indicate that: (1) There are significant uncertainties in the overall trends and interannual variability in LAI, NPP, and LST captured by the CMIP6 ESM. Meanwhile, simulated GPP and NEP trends were lower than observations with discrepancies reaching 0.03 yr⁻¹ for GPP and 2.46 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for NEP. (2) Spatially, CMIP6 ESMs exhibited widespread underestimation of trends in LAI, GPP, NPP, and NEP across China. The MME underestimated these variables in 46.29 % (LAI), 43.47 % (GPP), 49.81 % (NPP), and 61.34 % (NEP) of the study area. Meanwhile, the simulated LST trend is underestimated in northern China, while its overestimations in western and southern China. (3) ESMs inadequate responsiveness to anthropogenic and environmental forcing and incomplete mechanistic representation of plant respiration pathways struggled accurate simulation of trends in LAI, GPP, NPP, NEP and LST.

1 Introduction

Since the Industrial Revolution, human activities have significantly altered the structure of the atmosphere and terrestrial biosphere, leading to profound changes in the structure and functioning of terrestrial ecosystems (Allen et al., 2018). Vegetation, as a key component of the biosphere, regulates the exchange of carbon, water, momentum, and energy between the land and the atmosphere (Wu et al., 2020). Vegetations maximize water-use efficiency at the leaf scale by dynamically regulating stomatal conductance to effectively respond to rising atmospheric CO₂ concentrations and global warming (Fu et al., 2022). However, due to its high sensitivity to environmental conditions (Jung et al., 2017; Nemani et al., 2003), vegetation has undergone notable changes. The response of vegetation to global change has attracted significant attention in recent years (Chen et al., 2019; Hovenden et al., 2019; Zhu et al., 2016).

To capture the responses and feedbacks of terrestrial ecosystems to global environmental change, several Earth System Models (ESMs) in Coupled Model Intercomparison Project Phase 6 (CMIP6) simulate terrestrial vegetation structure and photosynthetic capacity (Song et al., 2021). Four key outputs of ESMs, which describe canopy structure and photosynthetic capacity – Leaf Area Index (LAI), Gross Primary Productivity (GPP), Net Primary Productivity (NPP), Net Ecosystem Productivity (NEP) – are critical components of the terrestrial carbon cycle (Le Quéré et al., 2009; Wang et al., 2015), and the main physical factors – Land Surface Temperature (LST) – controlling the land carbon balance (Anav et al., 2013; Piao et al., 2009). The relationship between LAI, GPP, and environmental factors is complex, as environmental drivers significantly influence LAI dynamics and collectively determine vegetation photosynthetic performance (Zhao et al., 2022, 2020). Vegetation productivity can be directly estimated using canopy parameters, resulting in a high coupling between GPP and LAI (Walther et al., 2019; Tian et al., 2024).

The LAI, defined as the total one-sided leaf area per unit ground surface area, is a key parameter of vegetation canopy structure that directly influences light interception, transpiration, and the spatial heterogeneity of GPP. GPP refers to the total amount of carbon dioxide fixed into organic compounds by vegetation through photosynthesis, serving as a core indicator of an ecosystem's carbon sequestration capacity. NPP represents the net carbon accumulation after subtracting autotrophic respiration from total photosynthetic fixation, reflecting the primary production potential and health of ecosystems as influenced by GPP and plant physiological regulation. NEP denotes the net carbon exchange between the ecosystem and the atmosphere by further subtracting heterotrophic respiration from NPP, making it a crucial measure for assessing regional carbon source/sink status under the influence of environmental factors such as atmospheric CO₂ and climate (Fang et al., 2001). LST is the thermodynamic temperature at the land-atmosphere interface, playing a key role in surface energy and water exchange while jointly affecting ecological processes through interactions with solar radiation, soil properties, vegetation, and atmospheric conditions (Li et al., 2023). In-depth research on LST facilitates a deeper understanding of surface-atmosphere exchange processes at global and regional scales and provides high-quality quantitative indicators of surface conditions for scientific applications. Consequently, LST has been designated as an indispensable observation indicator for the International Geosphere and Biosphere Program (IGBP) and the Global Climate Observing System (GCOS) (Townshend et al., 1994; Hollmann et al., 2013).

Both remote sensing and ecosystem models can be effective in investigating large-scale vegetation dynamics (Zhao et al., 2020). Ecosystem models offer an advantage in attributing vegetation growth changes by addressing limitations inherent in traditional statistical models, such as the inability to distinguish correlation from causation and to capture nonlinear relationships (Piao et al., 2006; Zhu et al., 2016). Earth System Models (ESMs), which simulate the physical, chemical, and biological processes and their interactions within the atmosphere, land, ocean, biosphere, and cryosphere (Ziehn et al., 2017; Wu et al., 2019; Bao et al., 2020; Zhang et al., 2020), are critical tools for modeling terrestrial ecosystems (Li et al., 2019a; Lawrence, 2020).

Understanding vegetation responses to past global environmental changes is crucial for improving predictions of future spatiotemporal changes in global vegetation (Zhao et al., 2020). However, numerous studies have evaluated and projected variables simulated by ESMs, revealing a general overestimation trend for model simulations (Anav et al., 2013; Mahowald et al., 2016; Kim et al., 2018; Zeng et al., 2016; Park and Jeong, 2021). Anav et al. (2013) reported that half of the models overestimated LAI trends from 1986 to 2005. Ziehn et al. (2017) found that the ACCESS-ESM1 model overestimated the seasonal mean and peak amplitude of global LAI during historical periods. Zhu et al. (2019) observed that ELMv1-ECA significantly overestimated LAI across most terrestrial surfaces, particularly in the tropics, while underestimating GPP in boreal forest systems and tropical ecosystems. Song et al. (2021) evaluated multiple ESMs and concluded that global LAI overestimation primarily stemmed from overestimated LAI in non-forest areas, with peak LAI in some regions occurring 1–2 months later than observed.

China's carbon sink accounts for approximately 10 %–31 % of the global terrestrial ecosystem carbon sink (Piao et al., 2022). More accurate simulations of terrestrial ecosystems using CMIP ESM models are highly valuable for both scientific research and policy development. Currently, CMIP has progressed to its sixth phase (CMIP6). Compared to CMIP5, CMIP6 models feature significantly improved resolution, incorporate more comprehensive biogeochemical processes (Eyring et al., 2016), and consider more plausible future scenarios (O'Neill et al., 2016), an advancement that will enhance the potential for simulating terrestrial ecosystems. However, research evaluating the outputs of CMIP6 models remains limited compared to CMIP5 evaluations. In existing evaluation studies, researchers have made significant progress in assessing the performance of ESMs. For example, Zhao et al. (2020) focused on the temporal variability of LAI in CMIP6 models, while Song et al. (2021) evaluated global LAI and tree height performance. However, limitations in spatial resolution and geographical scope have left gaps in the systematic evaluation of ESM simulations for specific regions, such as China. Similarly, Sun et al. (2023) conducted evaluations of LAI, GPP and NPP across Asia but lacked a comprehensive methodology and finer-scale analyses for historical periods. These gaps highlight the need for more detailed and region-specific assessments of CMIP6 ESMs outputs.

Based on the current lack of multi-variable evaluations in ESM assessment studies and the absence of systematic regional evaluations in China, this study utilizes remote sensing data as observational references to evaluate the interannual trends performance of LAI, GPP, NPP, NEP and LST simulated by 12 CMIP6 ESMs (ACCESS-ESM1-5, BCC-CSM2-MR, CanESM5, CESM2-WACCM, EC-Earth3-Veg, INM-CM4-8, INM-CM5-0, IPSL-CM6A-LR, MPI-ESM1-2-HR, MPI-ESM1-2-LR, NorESM2-LM, and NorESM2-MM) in China. A 3×3 sliding window approach is applied spatially to provide a more systematic assessment of trends performance. Finally, remote sensing data of tree cover were used to evaluate the simulated trends of tree distribution patterns.

2 Data and method

2.1 The study area

China is situated at the intersection of East and Central Asia, with a geographical range approximately between 18 and 53° N latitude and 73 and 135° E longitude. China exhibits diverse topographical features, ranging from expansive plains and hills in the east to high plateaus and mountains in the west, with an overall terrain gradient sloping from west to east. The monsoon region in the east has abundant precipitation and a large temperature difference from north to south, which is suitable for a wide variety of plants to grow, while the arid region in the west has scarce precipitation and sparse vegetation. Against this climatic background, China has developed a diversity of vegetation types ranging from humid broad-leaved evergreen forests to arid desert grasslands, with significant regional differences between seasonal and inter-annual variations in vegetation (Fig. 1).

Figure 1 Distribution of vegetation types in China.

2.2 CMIP6 Earth system models

The Coupled Model Intercomparison Project (CMIP), organized under the auspices of the World Climate Research Programme (WCRP) Working Group on Coupled Modelling (WGCM), which enables the simulation of difficult-to-observe early global coupled climate models by conducting experiments using atmospheric models coupled to the dynamical ocean, a simple land surface, and thermodynamic sea-ice (Meehl et al., 1997). CMIP6 adopts a novel, more coordinated organizational structure (Eyring et al., 2016).

This study utilizes the outputs of LAI and GPP variables from 12 CMIP6 models (ACCESS-ESM1-5, BCC-CSM2-MR, CanESM5, CESM2-WACCM, EC-Earth3-Veg, INM-CM4-8, INM-CM5-0, IPSL-CM6A-LR, MPI-ESM1-2-HR, MPI-ESM1-2-LR, NorESM2-LM, and NorESM2-MM) under five scenarios: historical, SSP1-26, SSP2-45, SSP3-70, and SSP5-85 as the simulated values of the model. To ensure consistency, the outputs from the first realization (r1i1p1) of each model were selected. Monthly data outputs from all models were uniformly resampled to a spatial resolution of 0.5°×0.5°. Additionally, based on whether the models incorporate a dynamic vegetation module, the 12 models were classified into groups with and without dynamic vegetation integration (Table 1).

Table 1 Spatial resolution of selected CMIP6 models and variables information.

To provide a more robust and integrated assessment than any single model can offer, we primarily utilize the Multi-Model Ensemble (MME) mean, calculated as the arithmetic mean of all available models, as a central benchmark for evaluating the simulated spatial patterns and interannual variability of vegetation and climatic variables over China. The MME mean-based approach, an established statistical integration technique, synthesizes outputs from diverse models through averaging (Zeng et al., 2016). This approach combines the true and noise signals, where the inter-model errors are mitigated through model averaging, thereby amplifying the underlying true signals.

2.3 Benchmark data

The 2003–2014 LAI data from the reprocessed MODIS product in China (Yuan et al., 2011), were used to validate the LAI outputs of 12 ESMs under the historical scenario. The reprocessed MODIS product refines the original MODIS LAI data (MOD15A2H and MYD15A2H) using modified temporal-spatial filtering (mTSF) technique to fill data gaps and correct low-quality data. Subsequently, the TIMESAT SG filtering technique was applied for post-processing to produce the final product. The reprocessed remote sensing LAI data were then resampled to a 0.5°×0.5° monthly grid for consistency with the model outputs. The analysis period spans 2003 to 2014. This timeframe was selected to utilize the stable post-launch era of the MODIS Terra-Aqua dual-satellite constellation for vegetation monitoring, commencing after potential initial sensor calibration transients, and to align with the end year of the CMIP6 historical experiment simulations used for comparison.

Solar-Induced Chlorophyll Fluorescence (SIF), which integrates the complex physiological functions of plants and can directly reflect the dynamic changes in plants' actual photosynthetic process, exhibits a strong linear relationship with GPP, thereby serving as a direct observational indicator for GPP (Mohammed et al., 2019; Frankenberg et al., 2014; Walther et al., 2016). Compared with GOSIF-GPP and FLUXCOM-GPP, CSIF exhibits higher temporal resolution and more robust spatial data gap-filling capability. Therefore, the CSIF dataset was used to validate the GPP outputs of the models under the historical scenario. CSIF data are generated by using surface reflectance data from MODIS C6 (MCD43C1) as input (Zhang et al., 2018), trained with daily SIF observations derived from the Orbiting Carbon Observatory-2 using a machine learning algorithm. CSIF effectively captures the seasonal dynamics of satellite-observed SIF and shows a strong correlation with ecosystem GPP, making it a suitable proxy for GPP in vegetation phenology analysis. The CSIF data were aggregated from 4 d to monthly scale data, and the spatial resolution was sampled into a 0.5°×0.5° grid using mean aggregation.

The NPP dataset was derived from the MOD17A3HGF Version 6 product of MODIS, featuring a spatial resolution of 500 m×500 m and a temporal resolution of 8 d. Annual NPP values were calculated as the cumulative sum of 8 d Net Photosynthesis (PSN) products for each year. The reprocessed remote sensing NPP data were then resampled to a 0.5°×0.5° monthly grid for consistency with the model outputs.

The NEP dataset was obtained from the reprocessed annual NEP data based on MOD17A3HGF V6 products, provided by the Loess Plateau Subcenter of the National Earth System Science Data Center. This dataset retains a spatial resolution of 500 m×500 m, the study resampled the remote sensing NEP data to 0.5°×0.5° monthly grid.

The LST data were sourced from the Global Daily 0.05° Spatiotemporally Continuous LST dataset hosted by the National Tibetan Plateau Data Center. The dataset reconstructs clear-sky LST by applying an Empirical Orthogonal Function (EOF) interpolation method to MODIS Terra/Aqua LST products. It further integrates ERA5-Land climate reanalysis data using a Cumulative Distribution Function (CDF) matching method, generating high-quality global spatiotemporally continuous LST data under both clear-sky and all-weather conditions. The final product has a spatial resolution of 0.05°×0.05°, while the remote sensing LST dataset were resampled to 0.5°×0.5° monthly grid in this study.

Tree cover fraction (TCF) were extracted from the MODIS Terra Vegetation Continuous Fields (VCF) product (https://www.earthdata.nasa.gov/data/catalog/lpcloud-mod44b-061, last access: 14 November 2025), providing global tree cover fraction at a 500 m×500 m spatial resolution. To consistency with the model output, the remote sensing TCF data were resampled to a 0.5°×0.5° monthly grid.

To verify the accuracy of CMIP6 ESMs in simulating climate data (temperature, precipitation, and solar radiation), we employed historical temperature and precipitation data sourced from the reanalyzed CRUNCEP dataset. Specifically, we utilized monthly data from the atmospheric stress component of CRUNCEP (https://gdex.ucar.edu/datasets/d314003/, last access: 14 November 2025) and standardized its spatial resolution to 0.5°×0.5°. Given that CRUNCEP does not include solar radiation data, this study also uses monthly-scale net solar radiation data from ERA5 (https://cds.climate.copernicus.eu/datasets/reanalysis-era5-land-monthly-means?tab=download, last access: 14 November 2025), with its spatial resolution resampled to 0.5°×0.5°.

2.4 Auxiliary data

The IGBP classification scheme from the MODIS MCD12Q1 product was used to categorize the vegetation in China into four land cover types: forest (Tree Cover), grassland (Grassland), cropland (Crop), and non-vegetated areas (Non-Vegetable). Specifically, evergreen coniferous forest, evergreen broadleaf forest, deciduous coniferous forest, deciduous broadleaf forest, mixed forest, woody savanna and savanna are combined to be classified as Tree Cover; closed shrubland, open shrubland, grassland and permanent wetland are combined to be classified as Grassland; farmland and agricultural land (with natural vegetation) are classified as Cropland; and urban and built-up land, permanent snow and ice, unutilized land, and water are combined to be classified as Non-Vegetable areas. The specific classification details are summarized in Table S1 in the Supplement.

2.5 Analysis

The Theil-Sen Median method was employed to estimate the interannual trends of observed and simulated LAI, GPP, NPP, NEP and LST from 2003 to 2014, while the Mann-Kendall (MK) test was applied to assess the statistical significance of these trends. Additionally, the performance of each model in simulating the interannual LAI, GPP, NPP, NEP and LST time series was quantified using the standard deviation (SD) of interannual time series, correlation coefficient (r), and root mean square error (RMSE).

Observed LAI and CSIF trends were used to validate the spatial trends of simulated LAI and GPP, providing a comprehensive spatial assessment of model performance. At the pixel scale, a 3×3 sliding window approach was applied to extract paired pixel value sequences from model simulations and satellite observations. These paired sequences were statistically compared using a paired t-test (also known as a dependent or correlated t-test). This method compares the means and standard deviations of two related groups to determine whether significant differences exist between them. The hypotheses for the paired t-test are formulated as follows: H₀: m=0 (null hypothesis, where m is the mean of differences between paired samples); H₁: m≠0 (alternative hypothesis).

$$\begin{array}{}\text{(1)}& {\displaystyle }t={\displaystyle \frac{m}{S/\sqrt{n}}}\text{(2)}& {\displaystyle }S=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(d_i-m{)}^{\mathrm{2}}}{n-\mathrm{1}}}}\text{(3)}& {\displaystyle }m={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{d}_i}{n}}\end{array}$$

Here, $i=\mathrm{1},\mathrm{2},\mathrm{\dots },n$. Where m represents the mean of paired sample differences, S denotes the standard deviation of these differences, and n is the number of paired samples. The critical t-value is determined by degrees of freedom (df) and significance level (α). If the calculated t-value is less than the critical value at α-level (P), this indicates P>α, leading to acceptance of H₀ and rejection of H₁ at the α significance level. Conversely, if the calculated t-value exceeds the critical value, H₀ is rejected in favor of H₁.

To enable a comparison between the trends of CSIF and GPP, it is necessary to normalize the annual-scale CSIF data and the GPP outputs from the ESMs. Normalization ensures that both datasets are brought to a comparable scale, facilitating a more accurate and meaningful trend analysis.

$$\begin{array}{}\text{(4)}& X={\displaystyle \frac{X_i-\text{mean}}{\text{std}}}\end{array}$$

where X is the annual normalized CSIF or GPP data, X_i is the annual unnormalized data, mean is the average of the multi-year annual data, and std is the standard deviation of the multi-year annual data.

3 Results

3.1 Trend of observed variables in China

From 2003 to 2014, most regions across China exhibited increasing trends in MODIS LAI (82.60 % of the area; Fig. S1 in the Supplement) and CSIF (86.50 %; Fig. S2), with significant increases covering 44.10 % for LAI and 48.34 % for CSIF concentrated in southern forests, northeastern China, and the Loess Plateau. Notably, CSIF demonstrated stronger rising trends than LAI in the North China Plain. MODIS NPP trends (Fig. S3) showed minimal significant changes, as 84.57 % of the area exhibited non-significant variability, while significant increases and decreases covered only 15.23 % and 0.20 % respectively. These patterns closely mirrored MODIS NEP trends (Fig. S4), where 85.07 % of the area showed non-significant changes with merely 0.28 % significant decreases. For MODIS LST during 2003–2019 (Fig. S5), 47.34 % and 52.66 % of the study area exhibited increasing and decreasing LST trends, respectively, yet 97.89 % of the area showed non-significant changes, with significant warming limited to 0.78 %.

Spatially coherent trends emerged across northwestern China (Xinjiang grasslands/croplands, central Inner Mongolia, Tibetan Plateau), where declining LAI and CSIF aligned with significant LST increases, while vegetation productivity declines were also observed in parts of southern Tibet. Conversely, northeastern China consistently exhibited significant increases in LAI, CSIF, NPP, and NEP, coinciding with declining LST trends. The Loess Plateau (northern Shaanxi/Ningxia) similarly demonstrated concurrent increases in LAI, CSIF, NPP, and NEP. Southern forested regions showed LAI and CSIF gains, though central Yunnan experienced notable CSIF reductions alongside NPP and NEP declines. Contrasting patterns characterized eastern China: the North China Plain had significant CSIF increases but prominent CSIF reductions later, alongside NEP declines and LST warming, while southeastern provinces (Guangdong, Fujian) and the Yangtze River Delta featured NPP and NEP reductions alongside localized CSIF decreases.

3.2 Evaluating the overall trends and interannual variability of CMIP6 variables

As shown in the bar chart comparing observed LAI trends with MME simulated LAI (Fig. 2a), the MODIS LAI exhibited a trend of 0.11 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, which falls within the MME uncertainty range (0.0060 ± 0.0058 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). This indicates that, overall, the models are capable of capturing the observed LAI trend. However, it is worth noting that the trend values simulated by individual models vary significantly, ranging from −0.0021 yr⁻¹ (INM-CM5-0) to 0.019 yr⁻¹ (NorESM2-MM), contributing to substantial uncertainties in trends of simulated LAI.

Figure 2 Overall annual average (a) LAI and (c) GPP trends in China during 2003–2014. The asterisk (*) indicates the significant trend (p<0.05). The Taylor diagrams compare the remotely-sensed and model-simulated annual mean (b) LAI and (d) GPP for the historical period (2003–2014). The standard deviation shows the interannual variability of the observed and simulated LAI. The dashed green lines show centered root mean square difference (RMSD) between model simulations and satellite observations.

A similar analysis was conducted for the normalized interannual trends of observed CSIF and simulated GPP in China during 2003–2014 (Fig. 2c). The observed CSIF exhibited a significant trend of 0.13 yr⁻¹, which substantially exceeded all model simulations, including the MME. Simulated GPP trends ranged from −0.020 yr⁻¹ (INM-CM5-0) to 0.10 yr⁻¹ (MME). These results demonstrate a systematic underestimation of GPP trends in China by DGVMs.

The Taylor diagram (Fig. 2b and d) is employed to evaluate and compare the interannual variations of LAI and GPP among individual models, the MME, and observations. According to the reference observational data, the standard deviation (SD) of the interannual time series for LAI and CSIF over the Chinese region during 2003–2014 is 0.041 and 0.49, respectively. The SD of ESM-simulated LAI ranges from 0.011 (MPI-ESM1-2-LR) to 0.089 (NorESM2-MM), while that of GPP ranges from 0.139 (ACCESS-ESM1-5) to 0.515 (INM-CM4-8). The MME-simulated SD for LAI and GPP are 0.021 and 0.425, respectively. The root mean square deviation (RMSD) of the models for LAI ranges from 0.024 (NorESM2-LM) to 0.063 (NorESM2-MM), and for GPP, from 0.334 (EC-Earth3-Veg) to 0.697 (INM-CM5-0). The MME-simulated RMSD for LAI and GPP are 0.022 and 0.249, respectively. Generally, the MME-simulated SD and RMSD for both LAI and GPP are smaller than those of most individual models. The correlation coefficients between the ESM-simulated interannual LAI time series and the observed LAI range from −0.409 (INM-CM5-0) to 0.855 (NorESM2-LM). The correlation coefficient for the MME-simulated LAI is 0.922, exceeding that of any individual model. Similarly, the correlation coefficients between the ESM-simulated GPP and observed CSIF time series range from −0.207 (INM-CM5-0) to 0.809 (NorESM2-MM), with the MME-simulated GPP showing a correlation coefficient of 0.847, which is higher than that of all individual models.

Figure 3a illustrates the observed NPP exhibited a trend of 2.70 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, which lies within the uncertainty range of the MME trend (1.19 ± 1.53 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), indicating that models broadly captured the direction of NPP changes. However, three models – ACCESS-ESM1-5 (−0.49 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), INM-CM5-0 (−1.42 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), and MPI-ESM1-2-LR (−0.27 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) – simulated declining NPP trends. In contrast, MODIS NPP showed a statistically significant increasing trend during 2003–2014, whereas only INM-CM4-8, NorESM2-LM, and the MME produced significant positive trends among the models. The observed NEP exhibited a statistically significant increasing trend (2.66 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), whereas simulated NEP trends were systematically lower than observations (Fig. 3c). The MME simulated NEP trend of 0.20 ± 0.64 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ – substantially below the observed value – with exceptionally large model simulations uncertainties. Specifically, ACCESS-ESM1-5 (−0.05 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), BCC-CSM2-MR (−0.21 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), CanESM5 (−0.76 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), and MPI-ESM2-LR (−0.14 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) simulated declining NEP trends, highlighting the models' underestimation of carbon sink intensification.

Figure 3 Overall annual average (a) NPP and (c) NEP trends in China during 2003–2014. The asterisk (*) indicates the significant trend (p<0.05). The Taylor diagrams compare the remotely-sensed and model-simulated annual mean (b) NPP and (d) NEP for the historical period (2003–2014). The standard deviation shows the interannual variability of the observed and simulated LAI. The dashed green lines show centered root mean square difference (RMSD) between model simulations and satellite observations.

The Taylor diagram (Fig. 3b and d) evaluates interannual variations of NPP and NEP across individual models, the MME, and observations. The SD of ESM-simulated NPP ranges from 7.70 (ACCESS-ESM1-5) to 24.95 (BCC-CSM2-MR), while RMSD values span from 8.82 (INM-CM4-8) to 25.66 (INM-CM5-0). For NEP, SD ranges from 9.19 (CanESM5) to 20.24 (BCC-CSM2-MR) with RMSD values between 13.55 (IPSL-CM6A-LR) and 21.41 (BCC-CSM2-MR). The MME demonstrates notably lower variability, with SD values of 5.65 for NPP and 6.19 for NEP, and RMSD values of 9.71 for NPP and 11.28 for NEP. Generally, the MME-simulated SD and RMSD for both variables are lower than those of all individual models.

The observed LST exhibited a declining trend of −0.0045 $\mathrm{°}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in China during 2003–2014 (Fig. 4a), which lies within the uncertainty range of the MME simulated a positive LST trend (0.0088 ± 0.035 $\mathrm{°}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Notably, only five models – ACCESS-ESM1-5, CESM2-WACCM, INM-CM5-0, MPI-ESM1-2-HR, and NorESM2-MM – captured negative LST trends, whereas all others simulated positive LST trends. The NorESM2-LM model produced a statistically significant positive LST trend.

Figure 4 (a) Overall annual average LST trends in China during 2003–2014. The asterisk (*) indicates the significant trend (p<0.05). (b) The Taylor diagrams compare the remotely-sensed and model-simulated annual mean LST for the historical period (2003–2014). The standard deviation shows the interannual variability of the observed and simulated LAI. The dashed green lines show centered root mean square difference (RMSD) between model simulations and satellite observations.

The Taylor diagram evaluates interannual variations of LST across individual models, the MME, and observations (Fig. 4). The SD of ESM-simulated LST ranges from 0.20 (INM-CM4-8) to 0.48 (IPSL-CM6A-LR), while RMSD values span from 0.30 (MPI-ESM1-2-HR) to 0.59 (IPSL-CM6A-LR). The MME demonstrates notably lower variability and error, with an SD of 0.10 and RMSD of 0.36.

3.3 Assessment of spatial patterns and trends in CMIP6 variables

From a spatial perspective, a comparison between simulated and observed LAI trends reveals that most models exhibit significant overestimation and underestimation in China (Fig. 5). The MME overestimated, underestimated, and simulated well with observed LAI trends in 21.71 %, 46.29 %, and 32.00 % in China, respectively. Specifically, the MME predominantly underestimated LAI trends in regions with statistically significant observed increases – particularly southern China, the Loess Plateau, and northeastern China – while overestimating trends in areas with declining or non-significant LAI changes (e.g., northwestern Xinjiang, the Tibetan Plateau, and parts of the North China Plain). Among individual models, EC-Earth3-Veg exhibited relatively better performance, but even this model aligns with 40.36 % of its simulated LAI trends matching observations (Fig. 5f). However, most models showed dominant underestimation patterns: CanESM5 (62.97 %), INM-CM4-8 (60.09 %), and INM-CM5-0 (63.50 %) had the largest proportions of underestimated LAI trends.

Figure 5 Evaluation of the LAI trend performance of the CMIP6 ESMs from 2003 to 2014 in China. A paired t-test with a sample size of 9 was conducted using a 3×3 sliding window to determine whether the model simulated a trend that was not significantly different from the observed data (Good), significantly smaller than the observed data (Underestimate), or significantly larger than the observed data (Overestimate). U, G, and O represent “Underestimate”, “Good”, and “Overestimate”, respectively.

The spatial distribution of normalized GPP trend discrepancies between simulated and observed CSIF exhibited similar patterns to that of the model LAI performance (Fig. 6). For the MME, overestimated, underestimated, and simulated well GPP trends accounted for 20.26 %, 43.47 %, and 36.26 % in China, respectively (Fig. 6a). Regions with overestimated GPP trends were predominantly distributed in Yunnan Province, the Sichuan-Tibet Plateau, and northwestern border areas. Among individual models, INM-CM4-8 exhibited the largest proportion of overestimated GPP trends (18.06 %), alongside 55.62 % underestimated areas (Fig. 6g). In contrast, MPI-ESM1-2-LR showed extreme underestimation biases, with 75.47 % of its simulated GPP trends being lower than observations – the highest underestimation proportion among all models – while only 7.68 % of areas were overestimated (Fig. 6k).

Figure 6 Evaluation of the normalized GPP trend performance of the CMIP6 ESMs from 2003 to 2014 in China. A paired t-test with a sample size of 9 was conducted using a 3×3 sliding window to determine whether the model simulated a trend that was not significantly different from the observed data (Good), significantly smaller than the observed data (Underestimate), or significantly larger than the observed data (Overestimate). U, G, and O represent “Underestimate”, “Good”, and “Overestimate”, respectively.

The spatial patterns of discrepancies between simulated and MODIS NPP trends during 2003–2014 are shown in Fig. 7. Simulated NPP trends were generally overestimated in southwestern and southern China but underestimated across most other regions (Fig. 7). For the MME, overestimated, underestimated, and simulated well NPP trends occupied 28.62 %, 49.81 %, and 21.57 % of the study area, respectively. Overestimated regions were predominantly concentrated in the Central Plains, Sichuan-Tibet Plateau, and Yunnan Province (Fig. 7a). Notably, the NorESM2-LM model showed the highest proportion of overestimated NPP trends (47.46 %) among all models, with underestimated and simulated well areas covering 33.84 % and 18.70 %, respectively. Its overestimations spanned southwestern China, southern China, the North China Plain, the middle Yangtze River Basin, and parts of northeastern China (Fig. 7l). Two other models – MPI-ESM1-2-HR and NorESM2-MM – also exhibited substantial overestimation (42.88 % and 46.07 %, respectively), while their underestimation proportions reached 29.86 % and 28.30 %. In contrast, CanESM5 and INM-CM5-0 displayed most underestimation biases, with 62.27 % and 69.85 % of their simulated NPP trends falling below observations.

Figure 7 Evaluation of the NPP trend performance of the CMIP6 ESMs from 2003 to 2014 in China. A paired t-test with a sample size of 9 was conducted using a 3×3 sliding window to determine whether the model simulated a trend that was not significantly different from the observed data (Good), significantly smaller than the observed data (Underestimate), or significantly larger than the observed data (Overestimate). U, G, and O represent “Underestimate”, “Good”, and “Overestimate”, respectively.

Figure 8 displays spatial discrepancies between simulated and MODIS NEP trends across China during 2003–2014. The MME underestimated, overestimated, and simulated well NEP trends in 61.34 %, 16.60 %, and 22.06 % in China, respectively (Fig. 8a), with spatial patterns similar to that of the simulated NPP trend spatial performance, indicating systemic uncertainties in coupled carbon flux simulations. Among individual models, NorESM2-LM exhibited the highest proportion of overestimated NEP trends (29.72 %), alongside 49.38 % underestimation and 20.90 % simulated well (Fig. 8i), while NorESM2-MM showed the lowest underestimation rate (35.38 %) with 27.12 % overestimation and 37.51 % simulated well (Fig. 8j). In contrast, all other models underestimated NEP trends in over 50 % of the region, with CanESM5 displaying the most severe underestimation bias (65.61 % underestimated, 12.33 % overestimated) (Fig. 8e).

Figure 8 Evaluation of the NEP trend performance of the CMIP6 ESMs from 2003 to 2014 in China. A paired t-test with a sample size of 9 was conducted using a 3×3 sliding window to determine whether the model simulated a trend that was not significantly different from the observed data (Good), significantly smaller than the observed data (Underestimate), or significantly larger than the observed data (Overestimate). U, G, and O represent “Underestimate”, “Good”, and “Overestimate”, respectively.

Figure 9 illustrates spatial differences between simulated and observed LST trends in China during 2003–2014. Most models overestimated LST trends, with the MME overestimating, underestimating, and simulating well observations in 54.57 %, 28.02 %, and 17.42 % of the study area, respectively (Fig. 9a). Spatially, the MME overestimated LST trends in southern China, northeastern China, Xinjiang, the Tibetan Plateau, and the Shaanxi-Gansu region compared to observations. While some models' proportion of underestimated LST trend – including ACCESS-ESM1-5 (45.36 %), CESM2-WACCM (46.19 %), INM-CM5-0 (63.99 %), and NorESM2-MM (59.59 %) – showed substantial underestimation biases, the majority of models still exhibited larger overestimated areas. Notably, NorESM2-LM displayed the most extreme overestimation, with 84.38 % of its simulated LST trends exceeding observations, except in localized underestimations over the North China Plain, central Inner Mongolia, and parts of Yunnan (Fig. 9l).

Figure 9 Evaluation of the LST trend performance of the CMIP6 ESMs from 2003 to 2014 in China. A paired t-test with a sample size of 9 was conducted using a 3×3 sliding window to determine whether the model simulated a trend that was not significantly different from the observed data (Good), significantly smaller than the observed data (Underestimate), or significantly larger than the observed data (Overestimate). U, G, and O represent “Underestimate”, “Good”, and “Overestimate”, respectively.

4 Discussion

4.1 ESM captures CMIP6 variables trends in China

Previous studies have found that, whether for CMIP5 or CMIP6, models have not effectively captured the long-term trends of vegetation (Anav et al., 2013; Song et al., 2021). Overestimations of LAI, GPP, and NPP are more common in the mid-latitude regions of Asia (Sun et al., 2023). This study also reveals that the models have misrepresented the trends of LAI and GPP in the China region, with significant underestimation of LAI and GPP growth in the forested areas of central and southern China, and overestimation of LAI and GPP trends in grassland areas of eastern Inner Mongolia and the Tibetan Plateau, where the growth or decline of LAI and GPP is less significant. Ecosystem models suggest that the global greening of vegetation is mainly driven by the increase in atmospheric CO₂ concentrations, while climate change and land-use changes are strong regional drivers of vegetation dynamics (Chen et al., 2019; Zhu et al., 2016). The misrepresentation of LAI and GPP trends in the models can be attributed to unreasonable parameter settings and missing mechanisms (Song et al., 2021). The overestimation of the growing season length (GSL) (Verger et al., 2016) in the models is a key factor contributing to the widespread overestimation of LAI in northern and temperate regions. Additionally, in the models, the reduction in water limitation mechanisms results in the overestimation of LAI in non-forest areas (Song et al., 2021). The overestimation of maximum carboxylation rates in the models leads to excessive biomass allocation to leaves, further overestimating LAI in some regions, compounded by missing mechanisms (Li et al., 2019a). The lack of adequate mechanisms for soil moisture and nutrient limitations in the models results in an overestimation of the positive drivers for increased LAI caused by elevated atmospheric CO₂ concentrations, temperature rise, and woody plant invasions (Huang et al., 2017; Reich et al., 2014; Wang et al., 2014).

Despite substantial progress in understanding vegetation autotrophic respiration through field experiments and observational studies, significant challenges remain in simulating respiratory processes and their responses to environmental changes due to inherent physiological complexities (Smith and Dukes, 2013; Yuan et al., 2024). Empirical studies suggest that the ratio of NPP to GPP often converges toward a constant value across regions, and vegetation adapts photorespiration to changes in photosynthesis through self-regulation (Gifford, 1995; Dewar et al., 1998; Collalti and Prentice, 2019). However, respiratory processes are inherently harder to quantify than photosynthesis, leading models to rely heavily on empirical functions for simulating NPP and NEP. Such parameterizations often fail to capture real-world vegetation respiratory acclimation. Spatial discrepancies in ESM simulated NPP trends (Fig. 7) mirror those of GPP trends (Fig. 6), with systematic underestimations in northeastern forested areas, the Loess Plateau, Inner Mongolian grasslands, and the middle-lower Yangtze River Basin, alongside overestimations over Tibetan Plateau grasslands. These biases indicate that ESMs primarily misrepresent trends in photosynthetic capacity, which leads to misestimation in simulated NPP. Satellite observations reveal pronounced NPP declines in Guangdong and Fujian provinces during 2003–2014 (Fig. 2c), yet models paradoxically overestimate NPP trends in these regions. Critically, the NPP and GPP ratio is not static but varies across ecosystems and responds to environmental drivers (Zhang et al., 2009; Amthor, 2000; DeLucia et al., 2007). For instance, temperature, precipitation, and rising CO₂ levels can enhance plant carbon-use efficiency by reducing photorespiration, thereby increasing NPP (Slot and Kitajima, 2015; Schuur et al., 2001; Aspinwall et al., 2017), leads to changes of the ratio of NPP and GPP (Zhang et al., 2009). Schuur et al. (2001) demonstrated that excessive rainfall in warm-humid regions (e.g., southern China) may leach nutrients, reduce light availability, and impair soil oxygen diffusion – all suppressing NPP. Current ESMs inadequately represent these nuanced respiratory responses to environmental, particularly in regions like Guangdong and Fujian.

ESMs face critical limitations in representing soil carbon-atmosphere interactions and respiratory mechanisms, largely due to insufficient understanding of large-scale vegetation autoregulatory respiration and soil heterotrophic respiration processes (Zhao et al., 2019; Tang et al., 2020). These models often lack mechanistic descriptions of soil respiration (Tang et al., 2020) and struggle to simulate soil carbon-climate feedbacks (Lehmann and Kleber, 2015). The spatial discrepancies between simulated and observed NEP trends in China (Fig. 8) closely resemble those of NPP trends (Fig. 7), as NEP – calculated as the difference between NPP and soil heterotrophic respiration – inherits biases from both components. Models fail to capture regional heterogeneity in soil respiration, and misestimation in simulated NPP (e.g., underestimation in afforestation zones or overestimation in grasslands) propagate directly into NEP biases. For instance, systematic underestimation of NPP in northeastern forests and the Loess Plateau (Fig. 8), combined with oversimplified soil respiration parameterizations, leads to misestimation in simulated NEP. This highlights a fundamental issue: current models cannot disentangle uncertainties in photosynthesis and respiration, nor represent their spatially divergent responses to environmental drivers. These limitations underscore the urgent need for mechanistic advancements in simulating soil-vegetation-atmosphere continuum dynamics.

In southern China and northeastern forested regions, simulated LST trends were significantly higher than observed values (Fig. 9), coinciding with systematic underestimation of LAI trends in these areas (Fig. 5). This contradiction highlights a critical model bias: increased LAI in forested ecosystems reduces LST through enhanced shortwave radiation interception, modified albedo, and altered aerodynamic and surface resistances (Li et al., 2015, 2019b). Vegetation growth typically cools surface temperatures by improving canopy shading and evapotranspiration efficiency – mechanisms poorly represented in current parameterizations. The models' failure to capture LAI increases thus leads to simulated LST warming trends, underscoring the necessity for improved representation of biophysical feedbacks in land surface.

The large inter-model spread and the model-observation mismatch in carbon cycle trends can be attributed to two primary sources: biases in the simulated climate and differences in model structure/parameterization. Our evaluation of the models' climatic outputs (Fig. S9 and Table S2) reveals that models did not accurately reproduce the observed climate over China for 2003–2014, generally overestimating mean temperature while underestimating mean precipitation and solar radiation. The high RMSD values, particularly for precipitation and radiation, indicate substantial errors in the simulated climatic drivers that propagate into the carbon cycle simulations. Furthermore, parameterization and model structure are fundamental for ecosystem models to generate realistic projections, playing a critical role in their accuracy (Luo et al., 2016). As synthesized from previous studies (Table S3), the selected ESMs exhibit considerable diversity in their key land surface components and related parameters (Spafford and Macdougall, 2021; Arora et al., 2020; Pan et al., 2025). This structural and parametric heterogeneity is a major factor contributing to the divergent performances in simulating the trends of carbon-cycle variables among the models. However, it is important to note that while this study focused on evaluating model performance in capturing interannual trends and variability, the assessment of seasonal cycle dynamics (e.g., phenological timing, amplitude of seasonal fluctuations) across China's diverse regions represents a critical avenue for future research to further refine our understanding of model capabilities in simulating terrestrial carbon cycle processes.

4.2 Potential anthropogenic impact for ESM inaccuracy

A series of ecological engineering initiatives, such as afforestation, reforestation, and grassland restoration in southern China and the Loess Plateau, have greatly promoted vegetation growth. However, the models may fail to capture the effects of these human activities. ESMs likely fail to adequately represent these anthropogenic impacts.

To better investigate the promoting effects of human activities, such as afforestation and the Grain-for-Green Program, on vegetation growth, we analysed the spatial trends of MODIS TCF changes in China during 2003–2014 (Fig. S6). The results indicate that 64.19 % and 32.49 % of the study area exhibited increasing and decreasing TCF trends, respectively. Among these, statistically significant increases and decreases accounted for 14.57 % and 1.28 % of the total area, while 84.15 % showed no significant trend. Most tree-abundant regions in eastern and central China displayed TCF gains. In contrast, declining TCF trends were observed in western China, where natural conditions are unsuitable for tree growth and land cover is predominantly grassland or non-vegetated areas. Significant TCF increases were concentrated in the Shaanxi-Gansu-Ningxia region, the Greater Khingan Mountains region in northeastern China, and Guizhou, Guangdong, and Guangxi provinces.

The observed TCF exhibited a significant increasing trend of 0.14 % yr⁻¹ (Fig. S7a), whereas simulated trends ranged from −0.0041 % yr⁻¹ (INM-CM4-8) to 0.096 % yr⁻¹ (MPI-ESM1-2-HR), with the MME yielding a trend of 0.039 % ± 0.032 % yr⁻¹. Critically, all models systematically underestimated TCF changes, and none captured the significant increase observed in satellite data.

Figure S8 illustrates spatial differences between simulated and observed TCF trends across China. While most models showed reasonable TCF trend simulations in the sparsely vegetated northwestern regions, they overestimated trends on the Tibetan Plateau and underestimated trends in eastern China. The MME underestimated, overestimated, and simulated well TCF trends in 43.72 %, 21.47 %, and 34.81 % in China, respectively. CanESM5 and MPI-ESM1-2-HR demonstrated relatively better simulated TCF trend with observations, well simulating TCF trends in 41.81 % and 42.79 % of the area, respectively, while overestimating 22.90 % and 21.76 % and underestimating 35.28 % and 35.45 % (Fig. S8d and i). However, both models still overestimated TCF changes in eastern coastal regions and parts of the Central Plains. Notably, MPI-ESM1-2-LR exhibited the highest proportion of overestimated TCF trends (27.92 %) among all models, with 36.24 % underestimation. In contrast, INM-CM5-0 showed the most underestimation, underestimating TCF trends in 53.29 % of the area while overestimating 12.46 %. The underestimation in regions with observed significant TCF increases (e.g., afforestation zones in southern China and the Loess Plateau) suggests that CMIP6 models likely fail to adequately represent anthropogenic contributions to vegetation growth, such as afforestation policies and ecological engineering impacts.

4.3 Uncertain of the observed datasets in China

Accurate observational data are essential for determining and improving the precision of models (Luo et al., 2016). In this study, two remote sensing datasets that better reflect actual vegetation growth, namely the reprocessed MODIS dataset and the CSIF dataset, were selected for evaluation. However, satellite remote sensing data in China are subject to considerable uncertainty due to various factors. High cloud cover during the rainy season and snow cover in high-latitude areas during winter can introduce inaccuracies. Additionally, satellite sensors are prone to degradation over time, leading to reduced sensitivity. Although the reprocessed MODIS and CSIF datasets utilize spatiotemporal filtering and machine learning techniques to enhance data quality (Zhang et al., 2018; Yuan et al., 2011), significant uncertainties remain. Similarly, limitations in satellite remote sensing-based carbon accounting (Araza et al., 2023), contribute to substantial uncertainties in MODIS NPP and NEP products (Sun et al., 2021; Huang et al., 2018; Ma et al., 2016). Concurrently, studies attribute uncertainties in MODIS LST primarily to spatial inconsistencies and surface emissivity uncertainties, the latter resulting from inadequate global representativeness in land cover classification (Ma et al., 2021; Wan et al., 2002; Duan et al., 2019).

5 Conclusion

This study evaluates the trends in performance of CMIP6 ESM-simulated LAI, GPP, NPP, NEP and LST in China from 2003 to 2014 using MODIS (LAI, NPP, NEP and LST) and CSIF, as observational references. Overall, ESMs fail to capture the trends of CMIP6 ESM-simulated variables in China. In particular, the simulations overestimate and underestimate that are spatially distributed over a large part of the area. The 12 CMIP6 ESMs revealed substantial errors in overall trends in LAI, NPP, and LST, along with significant underestimations in overall trends of GPP and NEP. Discrepancies between observed and simulated trends reached 0.03 yr⁻¹ for GPP and 2.46 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for NEP. Spatially, CMIP6 ESMs underestimated trends in LAI, GPP, NPP, and NEP across China. The MME underestimated these variables in 46.29 % (LAI), 43.47 % (GPP), 49.81 % (NPP), and 61.34 % (NEP) of the study area. Concurrently, models predominantly overestimated LST trends in southern China, while predominantly underestimating trends in northern China (MME overestimation: 54.57 % and underestimation: 28.02 %).

The CMIP6 ESMs exhibit significant misestimation of trends in carbon cycle physical and chemical variables, primarily attributed to inadequate response of environmental factor dynamics and insufficient consideration of anthropogenic influences. This limitation is compounded by the model's failure to adequately represent ecosystem and soil respiration mechanisms, which affects the simulation accuracy of NPP and NEP. Furthermore, systematic misestimation of LAI trend prevents accurate reproduction of the inhibitory effect of LAI on LST. The model simulations exhibit pronounced regional biases, including a systematic overestimation of ecological variables on the cold-arid Tibetan Plateau, a general underestimation over the extensively vegetated Loess Plateau likely linked to misrepresented ecological restoration processes, and widespread simulation errors in the rapidly urbanizing Pearl River Basin, potentially due to unaccounted-for anthropogenic pressures.