Ecosystem climate sensitivities drive the divergence in aerosol-induced carbon uptake across CMIP6 models

Anthropogenic aerosols significantly affect the terrestrial carbon cycle. Many models have been developed to simulate the effects of aerosols on regional ecosystem productivity. However, the differences among models in simulating the impacts of aerosols on gross primary production (GPP) remain unclear. To investigate the response of GPP to aerosol loadings among different models, we analyzed historical and hist-piAer simulations from five Earth System Models (ESMs) in Coupled Model Intercomparison Project Phase 6 (CMIP6). The results showed that all models captured the decrease in GPP from 1850 to 2014 (mean: −0.0586 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and the magnitudes of aerosol-induced GPP changes varied greatly (−0.0234 to −0.1151 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). To analyze the roles of aerosol representations and model sensitivities to climatic factors across ESMs, we developed a biophysical attribution framework. Our results showed that inter-model discrepancies in simulating the effects of aerosols on GPP were primarily driven by the differences in ecosystem climate sensitivities across ESMs, especially the response of photosynthesis to radiation and temperature. These findings provide critical insights into understanding the impacts of anthropogenic aerosols on the terrestrial ecosystem carbon cycle.

1 Introduction

Terrestrial gross primary production (GPP) is the largest carbon flux in the global carbon cycle (Anav et al., 2015; Lai et al., 2024). Understanding the response of GPP to various environmental factors is critical for accurately simulating the photosynthesis of terrestrial ecosystems (Piao et al., 2008; Huang et al., 2019). Atmospheric aerosol loadings have significantly increased since the Industrial Revolution due to the increased combustion of fossil fuels (Liu et al., 2022). The increased aerosol loadings significantly affect the amount of solar radiation reaching the Earth's surface (Tan et al., 2023), cloud properties (Manshausen et al., 2022), and regional climate (Najafi et al., 2015; Leung and Van Den Heever, 2023). Aerosols also play an important role in the photosynthesis of terrestrial ecosystems by altering the vegetation growing environment, such as radiation and temperature (Zhang et al., 2023b, 2019).

Atmospheric aerosols can affect GPP through four pathways. First, increasing aerosols can reduce the incoming radiation by absorbing and scattering sunlight (Wu et al., 2025). Second, aerosol loadings also increase the fraction of diffuse radiation (DF) reaching the Earth's surface. The increased DF can enhance the canopy light-use efficiency (LUE) (Gu et al., 2003, 2002). Third, aerosols can also influence the total radiation and DF reaching the surface by affecting the cloud properties (Khatri et al., 2021). Furthermore, aerosols also influence the terrestrial ecosystem photosynthesis through altering the air temperature and precipitation (Wang et al., 2018; Zhang et al., 2021a). To quantify the effect of aerosols on GPP, ground-based measurements and model simulations have been widely used.

Ground-based measurements provided some insights into the effect of aerosols on GPP at site scale. These studies showed that the increased aerosol loadings enhanced the canopy LUE by reducing the light saturation in the upper layers and enhancing photosynthesis in the lower canopy layers (Gu et al., 2002, 2003). Niyogi et al. (2004) showed that the effect of diffuse radiation induced by clouds and aerosols on canopy LUE varied with the vegetation types due to the canopy structure. Ground-based measurements also showed that the enhanced photosynthetic rates of sunlit and shaded leaves under high aerosol loadings conditions were driven by different environmental factors. The enhanced photosynthesis for sunlit and shaded leaves is induced by lower vapor pressure deficit (VPD) and higher diffuse radiation, respectively (Wang et al., 2018).

To investigate the changes of regional GPP induced by aerosols, model simulations were conducted. For example, Mercado et al. (2009) and Rap et al. (2018) showed that anthropogenic aerosols enhanced land carbon uptake due to the diffuse fertilization effect (DFE). Yue and Unger (2017) showed that the aerosol-induced change in net primary productivity (NPP) over China was from −3 % to 6 % depending on the local aerosol optical depth (AOD). However, these studies did not account for indirect aerosol radiative effects and aerosol climatic effects. To comprehensively understand the impact of aerosols, many other modelling studies were conducted. For example, Zhang et al. (2021a) reported that aerosols enhanced vegetation carbon dioxide sink since 1850 due to the DFE and cooling effects induced by aerosols. Zhang et al. (2023a) found that aerosols caused 0.43 % reduction in net biome production from 1980 to 2014 using the Community Earth System Model (CESM, version 2.1.3) and the dominant variable is the changes of temperature. Zhou et al. (2024) simulated the impact of the Clean Air Action plan on ecosystem carbon assimilation and found that aerosol reductions led to NPP increase of 20.1 ± 10.9 Tg C yr⁻¹ and the aerosol climatic effects on NPP are twice those of the aerosol radiative effects. These studies indicated that there were still large uncertainties in simulating the effects of aerosols on GPP (Liu et al., 2021; Zhang et al., 2021b).

The uncertainties in aerosol-induced GPP changes could be induced by aerosol direct and indirect effects and model sensitivities to climatic factors (defined as the ecosystem climate sensitivity). Bellouin et al. (2020) reported that there were large uncertainties in simulating aerosol radiative forcing. Additionally, many studies demonstrated that the parameterization of vegetation photosynthesis within Earth System Models (ESMs) also has large uncertainties (Hu et al., 2022; Gier et al., 2024). Liu et al. (2021) showed that current LUE models have large bias in estimating the DFE. However, it remains unclear whether the impact of aerosols on GPP simulated by different ESMs is consistent and the dominant factors driving divergence among different ESMs remain unclear. In this study, we used simulations with and without anthropogenic aerosol emissions from the Coupled Model Intercomparison Project Phase 6 (CMIP6). Our objectives of this study were: (1) to quantify the consistency among CMIP6 models in estimating the impacts of aerosols on terrestrial GPP; (2) to explore the contributors for inter-model differences. This multi-model assessment will enhance our understanding of the interactions between anthropogenic aerosols, climate, and terrestrial ecosystems.

2 Data and method

2.1 CMIP6 simulations

During the 1850–2014 period, the main anthropogenic aerosols simulated by CMIP6 models are Sulfate, Organic Carbon (OC), and Black Carbon (BC) (Zhang et al., 2022). The variations of anthropogenic aerosols are driven by anthropogenic and biomass burning emissions. To investigate the effect of anthropogenic aerosols on terrestrial GPP, we used the paired simulations from the Aerosol and Chemistry Model Intercomparison Project (AerChemMIP), a CMIP6-endorsed activity (Collins et al., 2017). We selected five Earth System Models (ESMs), including BCC-ESM1, IPSL-CM6A-LR, NorESM2-LM, MPI-ESM-1-2-HAM, and UKESM1-0-LL. These models have a diverse range of land surface components. Four of the five models considered the differential effects of direct and diffuse radiation on canopy photosynthesis (Table 1). For each model, we compared the historical experiment against the hist-piAer experiment from 1850 to 2014. The historical experiment is driven by all time-evolving natural and anthropogenic forcings, while the hist-piAer experiment is run in parallel with the historical experiment but fixes the anthropogenic aerosol emissions at pre-industrial levels. This experimental design can be used to calculate the variations of GPP induced by aerosols.

Table 1 CMIP6 Earth system models (ESMs) used in this study. For each model, the land component model and whether the model accounts for the diffuse fertilization effect (DFE) on canopy photosynthesis or not are listed.

The monthly GPP, surface downwelling shortwave radiation (rsds), near-surface air temperature (tas), top-of-atmosphere incident shortwave radiation (rsdt), precipitation (pr), total cloud cover percentage (clt), aerosol optical depth at 550 nm (od550aer) from historical and hist-piAer experiments were used in this study. We also used the leaf area index (LAI) from historical and hist-piAer experiments to calculate the fraction of absorbed photosynthetically active radiation (PAR) by using Beer-Lambert law. The model simulations can be downloaded from Earth System Grid Federation (ESGF). Only NorESM2-LM and UKESM1-0-LL historical experiments provide diffuse radiation datasets. To illustrate the impact of DF on vegetation photosynthesis, we calculated the clearness index (CI, rsds$/$rsdt). Previous studies demonstrated that CI was strongly correlated with DF (Zhang et al., 2023c). We also show the scatter plots of DF against CI from these two models (Fig. S1 in the Supplement). The results also indicate that there is a very high correlation between these two variables (R²=0.727). All data were regridded to a resolution of 1.25°×2.5° (latitude by longitude). The impact of aerosols on GPP was isolated by comparing historical and hist-piAer scenarios (historical-hist-piAer).

2.2 Observational data for model validation

In this study, monthly eddy covariance flux measurements from FLUXNET were used to assess the performance of GPP from CMIP6 ESMs. FLUXNET is a global network of eddy covariance towers, which can provide measured data on energy, water, and carbon dioxide exchanges between the biosphere and atmosphere (Pastorello et al., 2020). In this study, we used datasets from FLUXNET2015, which has over 1500 site-years of measurements from 212 locations (Lasslop et al., 2010). We utilized data records containing more than 80 % of measured values and good quality gap-filled data (NEE_VUT_REF_QC≥0.8) to ensure the quality of GPP. Figure S2 shows monthly GPP from five CMIP6 models (a–e) against FLUXNET site observations. The results reveal a systematic underestimation of high GPP (slopes 0.406–0.632) and low coefficient of determination (R²=0.305–0.438) at the site scale. Additionally, we also used the FLUXCOM-X products to evaluate the simulated GPP from CMIP6 ESMs. FLUXCOM-X is the global terrestrial GPP and evapotranspiration (ET) products derived from a newly data-driven scaling framework (X-BASE) (Nelson et al., 2024). Nelson et al. (2024) demonstrated that the X-BASE dataset was significantly improved compared to previous versions of FLUXCOM. The FLUXCOM-X products were also regridded into 1.25° in latitude and 2.5° in longitude. Figure S3 shows the performance of monthly GPP from five CMIP6 ESMs (a–e) against the FLUXCOM-X GPP. The coefficients of determination (R²) for these models range from 0.517 to 0.678, with root mean square errors (RMSEs) between 1.642 and 2.563 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The spatial distribution of observed and simulated GPP from 2001 to 2014 were also shown in Fig. S4.

2.3 Attribution Framework of Inter-Model Spread

The inter-model spread is attributed to discrepancies in simulated aerosol radiative and climatic effects and the sensitivities of model to climatic factors. To quantify the sources of uncertainty in aerosol-induced GPP changes, we developed an attribution framework based on the method of Yu and Huang (2023). The framework is based on the biophysical principle that GPP is the product of photosynthetically active radiation (PAR), fraction of absorbed PAR (fPAR) and LUE. GPP can be calculated as follows:

$$\begin{array}{}\text{(1)}& \text{GPP}=\text{PAR}\cdot \text{fPAR}\cdot \text{LUE}(\text{tas},\text{pr},\text{CI})\end{array}$$

where LUE is dependent on environmental conditions including tas, pr, and CI. To mathematically represent the aerosol-induced anomaly, a first-order Taylor expansion is applied to Eq. (1):

$$\begin{array}{}\text{(2)}& \mathit{\delta }\text{GPP}\approx {\displaystyle \frac{\partial \text{GPP}}{\partial \text{PAR}}}\cdot \mathit{\delta }\text{PAR}+{\displaystyle \frac{\partial \text{GPP}}{\partial \text{fPAR}}}\cdot \mathit{\delta }\text{fPAR}+{\displaystyle \frac{\partial \text{GPP}}{\partial \text{LUE}}}\cdot \mathit{\delta }\text{LUE}.\end{array}$$

δGPP is the mean difference in GPP induced by anthropogenic aerosols during 1850–2014 and can be calculated by subtracting the GPP of hist-piAer from that of historical for each ESM. Equation (2) can be rewritten as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathit{\delta }\text{GPP}& \approx (\text{fPAR}\cdot \text{LUE})\cdot \mathit{\delta }\text{PAR}+(\text{PAR}\cdot \text{LUE})\cdot \mathit{\delta }\text{fPAR}\\ & +(\text{PAR}\cdot \text{fPAR})\cdot \mathit{\delta }\text{LUE}.\end{array}\end{array}$$

Bloomfield et al. (2022) showed that a generalized linear mixed-effects model could well represent the response of LUE to environmental factors. In addition, aerosol-induced changes in climatic variables are small. Therefore, the change in LUE (δLUE) can be approximated linearly:

$$\begin{array}{}\text{(4)}& \mathit{\delta }\text{LUE}\approx {\displaystyle \frac{\partial \text{LUE}}{\partial \text{tas}}}\mathit{\delta }\text{tas}+{\displaystyle \frac{\partial \text{LUE}}{\partial \text{pr}}}\mathit{\delta }\text{pr}+{\displaystyle \frac{\partial \text{LUE}}{\partial \text{CI}}}\mathit{\delta }\text{CI}.\end{array}$$

Substituting Eq. (4) into Eq. (3) can get the full decomposition:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathit{\delta }\text{GPP}& \approx \text{fPAR}\cdot \text{LUE}\cdot \mathit{\delta }\text{PAR}+\text{PAR}\cdot \text{LUE}\cdot \mathit{\delta }\text{fPAR}\\ & +{\displaystyle \frac{\partial \text{LUE}}{\partial \text{tas}}}\cdot \text{fPAR}\cdot \text{PAR}\cdot \mathit{\delta }\text{tas}\\ & +{\displaystyle \frac{\partial \text{LUE}}{\partial \text{pr}}}\cdot \text{fPAR}\cdot \text{PAR}\cdot \mathit{\delta }\text{pr}\\ & +{\displaystyle \frac{\partial \text{LUE}}{\partial \text{CI}}}\cdot \text{fPAR}\cdot \text{PAR}\cdot \mathit{\delta }\text{CI}.\end{array}\end{array}$$

A multivariate regression model was constructed for specific plant functional type (PFT) and ESM to capture the impacts of climatic drivers and systematic model biases. The regression equation for a specific model (m) is defined as:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathit{\delta }{\text{GPP}}_m& \approx {\mathit{\beta }}_{\mathrm{0},m}+{\mathit{\beta }}_{\mathrm{1},m}\mathit{\delta }{\text{PAR}}_m+{\mathit{\beta }}_{f,m}{\text{PAR}}_{\text{clim}}\cdot \mathit{\delta }{\text{fPAR}}_m\\ & +{\mathit{\beta }}_{\mathrm{2},m}{\text{PAR}}_{\text{clim}}\cdot \mathit{\delta }{\text{tas}}_m+{\mathit{\beta }}_{\mathrm{3},m}{\text{PAR}}_{\text{clim}}\cdot \mathit{\delta }{\text{pr}}_m\\ & +{\mathit{\beta }}_{\mathrm{4},m}{\text{PAR}}_{\text{clim}}\cdot \mathit{\delta }{\text{CI}}_m.\end{array}\end{array}$$

Here, δGPP_m represents the aerosol-induced anomaly of GPP (Historical-Hist-piAer) from model m. PAR_clim is the climatological baseline PAR (0.45⋅rads). β₁ is the product of fPAR and LUE, while β_f represents the baseline LUE. β₂₋₄ represents the product of fPAR and the partial derivatives of LUE to climatic factors, while the intercept β₀ represents the systematic bias of the model. To address multicollinearity, standardized ridge regression was used for the specific PFT.

To quantify the inter-model divergence, we calculated the deviation of model m from the multi-model ensemble mean (mmm) by using the Eq. (7).

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }\left(\mathit{\delta }\text{GPP}\right)=\mathit{\delta }{\text{GPP}}_m-\mathit{\delta }{\text{GPP}}_{\text{mmm}}\end{array}$$

By substituting the regression equations into Eq. (7) and rearranging terms, we derived the final equation:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{\Delta }\left(\mathit{\delta }\text{GPP}\right)& =\underset{\text{State contribution}}{\underbrace{{\sum }_i{\mathit{\beta }}_{i,\text{mmm}}(X_{i,m}-X_{i,\text{mmm}})}}\\ & +\underset{\text{Sensitivity contribution}}{\underbrace{\left\{\begin{array}{c}{\sum }_i{X}_{i,m}({\mathit{\beta }}_{i,m}-{\mathit{\beta }}_{i,\text{mmm}})\\ +({\mathit{\beta }}_{f,m}{X}_{f,m}-{\mathit{\beta }}_{f,\text{mmm}}{X}_{f,m})\\ +({\mathit{\beta }}_{\mathrm{0},m}-{\mathit{\beta }}_{\mathrm{0},\text{mmm}})\end{array}\right\}}}\end{array}\end{array}$$

where X_i represents the independent variables (including interaction terms). This equation decomposes the model spread into two components:

1. State contribution: The divergence arising from differences in the simulated aerosol radiative and climatic effects (e.g., differences in simulated pr: $X_{i,m}-X_{i,\text{mmm}}$), weighted by the mean pr sensitivity (β_i,mmm).

2. Sensitivity contribution: The divergence arising from ecosystem climate sensitivities across ESMs. This term represents the contributions from dynamic photosynthesis sensitivity differences (${\mathit{\beta }}_{i,m}-{\mathit{\beta }}_{i,\text{mmm}}$), structural feedback divergence driven by LAI simulations (${\mathit{\beta }}_{f,m}{X}_{f,m}-{\mathit{\beta }}_{f,\text{mmm}}{X}_{f,m}$) and the systematic bias differences (${\mathit{\beta }}_{\mathrm{0},m}-{\mathit{\beta }}_{\mathrm{0},\text{mmm}}$).

3 Results

3.1 The changes of global GPP induced by aerosols

Figure 1 shows the spatial patterns of aerosol-induced changes in GPP from five CMIP6 ESMs and their multi-model mean from 1850 to 2014. In the multi-model ensemble mean (Fig. 1a), aerosol loadings lead to a reduction in GPP (0.0586 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), with 70.31 % of global land areas experiencing decreased GPP. Notably, the Northern and Southern Hemispheres exhibit contrasting responses: most areas in the Northern Hemisphere show decreased GPP, whereas some regions in the Southern Hemisphere show modest increases. The results suggest that the changes of GPP induced by aerosols in Northern and Southern Hemispheres are asymmetric. Large positive GPP anomalies can be observed around 30° S, while a pronounced decline in GPP is shown around 70° N.

Figure 1 The spatial pattern of changes in ecosystem GPP ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) induced by aerosols from CMIP6 models (a multi-model mean, b BCC-ESM1, c IPSL-CM6A-LR, d MPI-ESM-1-2-HAM, e NorESM2-LM, f UKESM1-0-LL).

Among the individual models, BCC-ESM1 shows the general spatial pattern of the multi-model ensemble mean but exhibits stronger regional variability. More areas show positive GPP anomalies in BCC-ESM1 than in the multi-model ensemble mean over the Northern Hemisphere. Weak aerosol effects are simulated by IPSL-CM6A-LR with small changes in GPP. In contrast, MPI-ESM1-2-HAM model reveals large GPP reductions across central and southern China, western Europe, and the eastern United States. All these regions have relatively high aerosol emissions. NorESM2-LM shows increased GPP in some areas of South America and decreased GPP over Europe. UKESM1-0-LL shows increased GPP in the eastern United States, central Africa, and parts of South America and decreased GPP across northern Europe and Asia.

In these models, decreased GPP can be found in northern mid-to-high latitudes and eastern China. The result shows a robust signal of aerosol-induced suppression of photosynthesis in these regions. Reduced GPP can be observed in 58.18 %, 68.64 %, 77.91 %, 53.86 %, and 72.94 % of the global regions for BCC-ESM1, IPSL-CM6A-LR, MPI-ESM-1-2-HAM, NorESM2-LM, and UKESM1-0-LL, respectively. The mean aerosol-induced GPP change is −0.0398, −0.0234, −0.0959, −0.0298, and −0.1151 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ from BCC-ESM1, IPSL-CM6A-LR, MPI-ESM-1-2-HAM, NorESM2-LM, and UKESM1-0-LL, respectively. Although all models show the reduction in GPP, the magnitude and spatial distribution of GPP changes vary greatly among models. These discrepancies suggest the uncertainty in quantifying aerosol impacts on the terrestrial carbon cycle in current ESMs.

In all four seasons, the CMIP6 models consistently show negative GPP anomalies induced by aerosols. However, there are large differences in the magnitude and spatial distribution (Figs. S5–S8). In March–May, the multi-model ensemble mean reveals widespread GPP reductions over the mid- and high-latitudes of the Northern Hemisphere, particularly across East Asia and Europe. The differences among the models are significant, with BCC-ESM1, MPI-ESM-1-2-HAM, and UKESM1-0-LL simulating stronger reductions, while IPSL-CM6A-LR and NorESM2-LM show weaker responses of GPP to aerosols. During the period of June–August, the variations of GPP are greater and there are more regions showing positive anomalies. The differences among models become more pronounced, especially in the low- and mid-latitudes. For example, simulations from BCC-ESM1 and UKESM1-0-LL show that the impacts of aerosols on GPP are positive in most regions of the United States and Europe, whereas MPI-ESM-1-2-HAM reveals that the GPP anomalies are negative in half of these regions. Meanwhile, the IPSL-CM6A-LR model simulation indicates that the changes of GPP are negative over the United States. In September–November, the negative anomalies are also shown but weaker than those during the period of June–August. Over Australia, the aerosol effects on GPP are positive from BCC-ESM1 and NorESM2-LM, but negative from IPSL-CM6A-LR and MPI-ESM-1-2-HAM. In December–February, aerosols consistently exhibit a small negative effect on GPP in the Northern Hemisphere, whereas model simulations show large discrepancies in the Southern Hemisphere. In all, these results demonstrate that aerosols generally suppress global GPP, but with significant differences in the amplitude and spatial distribution among models.

Figure 2a shows the time series of the ten-year average changes of GPP induced by aerosols and AOD variations from 1850 to 2014. The results show that global GPP and AOD have experienced large variations. The aerosol-induced GPP changes show a decreasing trend from 1850 to the mid-20th century, with a marked shift around the 1950s. This reduction was induced by the increasing aerosol emissions. From 1850 to 1890, some models show a positive impact of aerosols on GPP. However, the increment in ΔAOD during this period is negligible across all models. This indicates that anthropogenic aerosols have little impacts on the change of GPP from 1850 to 1890. The absence of a consistent directional GPP response suggests that these variations might be related to the internal climate variability noise. After 1980, an increase in GPP can be observed. This aligns with the decreasing aerosol loadings. There are notable differences between the models in the magnitude and timing of GPP and AOD changes. MPI-ESM-1-2-HAM and UKESM1-0-LL exhibit a larger variation of GPP compared to the other models. Figure 2b shows the scatter plots of the annual mean of GPP changes induced by aerosols against the AOD variations. Aerosol-induced changes in GPP are significantly related to the AOD (p<0.001). However, the sensitivities of GPP to aerosol loadings are different among models. These discrepancies highlight the uncertainty in simulating atmospheric aerosol loadings and the impact of aerosols on global productivity.

Figure 2 (a) Time series of aerosol-induced GPP changes (solid lines) and AOD variations (dashed lines) from 1850 to 2014 with a ten-years moving window; (b) The scatter plots of annual mean of GPP changes induced by aerosols against the AOD variations.

3.2 Changes of aerosols and meteorological factors

Analysis of AOD differences at 550 nm (od550aer) between historical and hist-piAer experiments reveals significant discrepancies among CMIP6 ESMs (Fig. 3). The multi-model ensemble mean differences in AOD (Fig. 3a) show significant increase in AOD across northern mid-latitudes, especially in major industrial regions including North America, Europe, and East Asia. In contrast, IPSL-CM6A-LR (Fig. 3c) shows decreased AOD in some regions. MPI-ESM-1-2-HAM (Fig. 3d) and UKESM1-0-LL (Fig. 3f) show high aerosol loadings in industrialized regions of North America and Eurasia. NorESM2-LM (Fig. 3e) shows a relatively modest aerosol increase. The spatial distribution of AOD reveals substantial inter-model discrepancies in simulating the global aerosol loadings.

Figure 3 The spatial pattern of mean differences of aerosol optical depth (AOD) at 550 nm (od550aer) between historical and hist-piAer experiments over the period 1850–2014 (a multi-model mean, b BCC-ESM1, c IPSL-CM6A-LR, d MPI-ESM-1-2-HAM, e NorESM2-LM, f UKESM1-0-LL).

Analysis of mean differences in rsds between historical and hist-piAer scenarios further demonstrates pronounced inter-model variability (Fig. 4). The multi-model ensemble mean differences in rsds (Fig. 4a) show reduction in most of the regions, especially in northern mid-latitudes. Models such as BCC-ESM1 (Fig. 4b) and UKESM1-0-LL (Fig. 4f) show similar spatial distribution of variations in rsds. Among all models, MPI-ESM-1-2-HAM (Fig. 4d) has the largest area where rsds increases. NorESM2-LM (Fig. 4e) show widespread decreases. IPSL-CM6A-LR (Fig. 4c) presents a more complex spatial distribution, highlighting reductions in shortwave radiation over parts of Eurasia and Africa but also notable regional increases. This is consistent with the spatial distribution of AOD. The result shows that there are large differences among ESMs in simulating radiation.

Figure 4 The spatial distribution of mean differences of shortwave radiation (rsds, W m⁻²) between historical and hist-piAer experiments from 1850 to 2014 (a multi-model mean, b BCC-ESM1, c IPSL-CM6A-LR, d MPI-ESM-1-2-HAM, e NorESM2-LM, f UKESM1-0-LL).

The spatial distribution of mean differences in the CI between historical and hist-piAer scenarios is consistent with that of rsds. BCC-ESM1 (Fig. 5b) and UKESM1-0-LL (Fig. 5f) show the decreases in CI in most regions. The CI decreases significantly in northern mid-latitude regions with the high aerosol loadings. In the NorESM2-LM (Fig. 5e), CI shows a decreasing trend in almost all regions. IPSL-CM6A-LR (Fig. 5c) reveals heterogeneous changes, with localized areas showing increased CI amidst predominant decreases. The result also shows large discrepancies in simulating the variations of CI induced by the aerosols.

Figure 5 The spatial pattern of mean differences of clearness index (CI) between historical and hist-piAer experiments from 1850 to 2014 (a multi-model mean, b BCC-ESM1, c IPSL-CM6A-LR, d MPI-ESM-1-2-HAM, e NorESM2-LM, f UKESM1-0-LL).

The spatial distribution of mean differences in tas between historical and hist-piAer scenarios reveals consistent cooling trends in response to increased aerosol loadings across these models in CMIP6 (Fig. 6). The multi-model ensemble mean differences of tas (Fig. 6a) shows that aerosols can induce the decrease in tas. BCC-ESM1 (Fig. 6b), MPI-ESM-1-2-HAM (Fig. 6d), and UKESM1-0-LL (Fig. 6f) prominently show widespread cooling in the Northern Hemisphere. The tas in high-latitude areas of the Northern Hemisphere decreases more than that in mid- and low-latitude regions. The temperature in the Southern Hemisphere is less affected than that in the Northern Hemisphere. NorESM2-LM (Fig. 6e) and IPSL-CM6A-LR (Fig. 6c) also exhibit decreased tas, but with weaker magnitude.

Figure 6 The spatial pattern of mean differences of near-surface air temperature (tas, °C) between historical and hist-piAer experiments from 1850 to 2014 (a multi-model mean, b BCC-ESM1, c IPSL-CM6A-LR, d MPI-ESM-1-2-HAM, e NorESM2-LM, f UKESM1-0-LL).

Analysis of mean differences in pr between historical and hist-piAer scenarios reveals that there are also large differences in simulating the impact of aerosols on pr across CMIP6 models (Fig. 7). The multi-model ensemble mean differences of pr show that aerosols induce a reduction in pr in most regions (Fig. 7a). BCC-ESM1 (Fig. 7b) shows modest decrease in pr across mid-latitudes and tropical regions. In contrast, IPSL-CM6A-LR (Fig. 7c), MPI-ESM-1-2-HAM (Fig. 7d), NorESM2-LM (Fig. 7e), and UKESM1-0-LL (Fig. 7f) show complex spatial distribution of pr changes, with both pronounced regional increases and decreases. The response of pr to aerosols suggests that there are some uncertainties in aerosol-cloud interactions. This may also induce the uncertainties in simulating regional hydrological cycles.

Figure 7 The spatial pattern of mean differences of precipitation (pr, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) between historical and hist-piAer experiments from 1850 to 2014 (a multi-model mean, b BCC-ESM1, c IPSL-CM6A-LR, d MPI-ESM-1-2-HAM, e NorESM2-LM, f UKESM1-0-LL).

3.3 Attribution of Inter-Model Spread in aerosol-induced GPP changes

We applied the attribution framework to quantify the drivers of inter-model spread in aerosol-induced GPP anomalies. The framework decomposes the total spread into contributions from aerosol radiative and climatic effects (“state”) and ecosystem climate sensitivities (“sensitivity”). The framework captures the variability of aerosol-induced GPP changes for all ESMs, with the R² ranging from 0.514 to 0.788 (Fig. 8c, red line). The framework can explain more than 50 % of GPP changes. We also showed the performance of the attribution framework across ESMs per PFT (Table S1 in the Supplement). The coefficient of determination (R²) exceeds 0.6 across different ESMs and PFTs. These results suggest that the framework can be used for analyzing the contribution of “state” and “sensitivity”.

Figure 8 Decomposition of drivers governing the inter-model spread in aerosol-induced GPP anomalies. Contributions ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) arising from differences in (a) aerosol radiative and climatic effects (State; including ΔPAR and interaction terms) and (b) ecosystem climate sensitivities (Sensitivity; including dynamic Δβ and residual). (c) Aggregated relative contributions (bars) and coefficient of determination (R², line). Components: State (blue), Sensitivity (orange), Residual (green).

The decomposition of total spread reveals the dominant driver (Fig. 8c, bars). The GPP anomalies in NorESM2-LM are dominated by state contributions (80.5 %). This indicates that the aerosol radiative and climatic effects simulated by this model have a large difference with the ensemble mean. In contrast, the main driver of BCC-ESM1, IPSL-CM6A-LR, MPI-ESM-1-2-HAM, and UKESM1-0-LL is the sensitivity contributions, accounting for 143.6 %, 165.2 %, −71.8 %, and −138.2 %, respectively. This implies a divergence in model ecological parameterization. Four of five models show that the sensitivity contribution is higher than the state contribution in the inter-model spread. The inter-model discrepancies in aerosol-induced GPP changes are driven by the parameterization of canopy photosynthesis in the ESMs.

Figure 8a shows the contribution of discrepancies in simulated aerosol radiative and climatic effects across models. IPSL-CM6A-LR shows that radiation anomalies have a large positive contribution (0.0282 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). This indicates that this model simulates a weaker aerosol dimming effect than other models (Fig. 4c). This is the primary driver for positive GPP anomalies of IPSL-CM6A-LR. For UKESM1-0-LL, radiation shows a strong negative contribution (−0.0174 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). This means that the aerosol-induced dimming effect simulated by UKESM1-0-LL is stronger than those of other models (Fig. 4f). For NorESM2-LM, the interaction between PAR and CI plays a prominent role (0.0198 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). The interaction between PAR and CI also plays an important role in IPSL-CM6A-LR, but with a negative effect (−0.0237 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$).

Figure 8b shows the contribution of ecosystem climate sensitivities to the inter-model spread in GPP anomalies. The parameterization of model canopy photosynthesis will induce the divergence in ecosystem climate sensitivities to temperature, precipitation, and radiation. Large difference in the response of photosynthesis to radiation can be observed. This indicates that the assumption of canopy radiative transfer in the ESMs introduces some differences in the response of photosynthesis to radiation. β_PAR and β_CI represent the sensitivity of vegetation photosynthesis to light quantity and light quality (light composition and spectral distribution), respectively. The sensitivities of photosynthesis to CI (β_CI) from BCC-ESM1 (0.0250 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and UKESM1-0-LL (0.0662 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) show positive contributions to GPP anomalies. This suggests that these models simulate a stronger diffuse fertilization effect than the multi-model ensemble mean. For UKESM1-0-LL, the diffuse fertilization effect partially reduces by the covaried decreased PAR (β_PAR: −0.0667 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). IPSL-CM6A-LR shows the opposite pattern with a positive contribution from PAR sensitivity (0.0151 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and a negative contribution from CI sensitivity (−0.0150 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). This indicates that the assumption of canopy radiative transfer in this model is insensitive to the variations of light quality. The sensitivity of photosynthesis to temperature reveals the different thermal adaptation strategies used in the ESMs. MPI-ESM-1-2-HAM shows the highest positive contribution of temperature (β_tas: 0.0365 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). This suggests that this model has a larger GPP gain under aerosol-induced cooling compared to the multi-model ensemble mean. Conversely, BCC-ESM1 shows a negative contribution of temperature sensitivity (−0.0168 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). This indicates that the photosynthesis in the model gains less than other ESMs. The sensitivity of photosynthesis to precipitation (β_pr) links to the soil hydrology schemes and the response of stomatal conductance to water stress. The contribution of precipitation sensitivity to GPP anomalies is generally smaller than that of other environmental factors across most models. This suggests that the response of ESMs to precipitation is relatively consistent. BCC-ESM1 and IPSL-CM6A-LR show a moderate positive contribution of precipitation sensitivity (0.0134 and 0.0098 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), while other three models show a small negative contribution. We also calculate the intercept (residual) term, which can represent the systematic deviation in baseline productivity from the ensemble mean. For MPI-ESM-1-2-HAM and NorESM2-LM, the baseline bias shows large negative contributions (−0.0297 and −0.0111 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively). This demonstrates that the offsets in base-state parameterizations (e.g., lower Vcmax compared to the MMM) is also a major source of uncertainty in simulating the effect of aerosols on plant productivity.

To analyze the spatial distribution of dominant driver governing the inter-model spread in aerosol-induced GPP anomalies, we showed the state and sensitivity contributions onto the global spatial grid (Fig. S9). The results from spatial distribution further confirmed that the sensitivity contribution is the primary driver of GPP anomalies induced by aerosols. The contribution of ecosystem climate sensitivities ranges from 57.5 % in NorESM2-LM to 78.7 % in UKESM1-0-LL. Although the primary driver for each model is the sensitivity contribution, the spatial maps reveal significant inter-model divergence in the directionality of the sensitivities. MPI-ESM-1-2-HAM and UKESM1-0-LL show widespread negative sensitivity contributions in most regions, while BCC-ESM1, MPI-ESM-1-2-HAM, and UKESM1-0-LL exhibit widespread negative state contributions in most regions.

4 Discussion

4.1 Divergent aerosol impacts on terrestrial GPP in CMIP6 models

This study revealed significant inter-model divergence among CMIP6 ESMs in simulating the effect of aerosols on terrestrial GPP. In this study, we showed that the anthropogenic aerosols decreased the terrestrial GPP in all ESMs from 1850 to 2014 and the decreased GPP increased with the AOD. Zhang et al. (2023a) and Zhou et al. (2024) also showed that the anthropogenic aerosols caused a reduction of terrestrial carbon sink. This is generally consistent with the results in this study. However, Zhang et al. (2021a) reported that anthropogenic aerosols increased terrestrial carbon sink by 22.6 Pg C using ORCHIDEE_DF land components and IPSL-CM6A-LR climate and aerosol forcing data. This bias might be induced through two pathways. First, AOD from IPSL-CM6A-LR is lower than that from other models. Second, a new development of ORCHIDEE trunk (ORCHIDEE_DF) with two-stream canopy light transmission model was used, which can better capture the diffuse radiation fertilization effect than the original IPSL-CM6A-LR model (Zhang et al., 2020).

4.2 Ecosystem climate sensitivities of ESMs in CMIP6

Our attribution framework reveals that the substantial inter-model spread in aerosol-induced GPP anomalies is not merely caused by the divergent aerosol radiative and climatic effects, but is also governed by the terrestrial ecosystem climate sensitivities. This finding challenges the traditional view of improving the simulation of aerosol radiative and climatic effects alone and suggests that ecosystem climate sensitivities of ESMs are also a dominant source of uncertainties in simulating aerosol-induced GPP changes. Inter-model differences in ecosystem climate sensitivities primarily came from the radiation and temperature sensitivities of photosynthesis in the ESMs.

A critical source of divergence lies in the representation of canopy radiative transfer, specifically the response to light quality. Compared to the MMM, UKESM1-0-LL shows strong positive anomalies in the sensitivity to diffuse radiation. The sensitivity coefficients for CI from UKESM1-0-LL are lowest in 8 out 11 PFTs (Table S6). This behavior is controlled by the land component (JULES). The model uses a multi-layer canopy scheme with explicit light interception calculations for sunlit and shaded leaves at each depth (Sellar et al., 2019; Clark et al., 2011). This canopy radiative transfer model allows the diffuse radiation to reach the deeper canopy and can capture the DFE. The IPSL-CM6A-LR shows the lowest DFE among all five models, because the model uses a standard “big leaf” approach and does not consider the DFE (Cheruy et al., 2020; Zhang et al., 2020). Table S6 shows that sensitivity coefficients for CI from IPSL-CM6A-LR approach zero across almost all PFTs. NorESM2-LM (CLM5) integrates a revised two-stream approximation with the Medlyn stomatal conductance model (Lawrence et al., 2019). This combination also can capture the DFE (β_CI<0 in major PFTs). However, this model shows lower DFE than other models and this might be induced by the nutrient limitation.

Most land components of ESMs use the traditional Farquhar-Berry-Collatz framework for simulating photosynthesis (Arora, 2003; Clark et al., 2011; Reick et al., 2021; Boucher et al., 2020). However, our results showed that there were large differences in the temperature sensitivities (β_tas) and structural feedbacks (fPAR) across different PFTs (Tables S3 and S4). This is because there are large differences in how they consider the influence of temperature on GPP. The impact of temperature on photosynthesis can be divided into three parts, including chemical limits, adaptation to heat, and phenology. For the level of immediate chemical reactions, the models have some differences in estimating the response of V_cmax and J_max to temperature. MPI-ESM-1-2-HAM (JSBACH 3.2) adopts a strict physical chemistry approach. The method uses the Arrhenius equation and a specific inhibition function to calculate the impact of temperature on key rates (Reick et al., 2021). When the temperature is higher than about 55 °C, the photosynthetic rates will be zero. In contrast, UKESM1-0-LL (JULES) adopts a formula based on Q₁₀ factors (Clark et al., 2011). These divergent responses are captured by the β_Tas values (Table S4). At the level of adaptation to heat, NorESM2 (CLM5) introduces a mechanism for thermal acclimation based on the LUNA module (Lawrence et al., 2019). BCC-ESM1 and IPSL-CM6A-LR use a traditional method, which assumes that the response of plants to temperature is fixed (Boucher et al., 2020; Li et al., 2019). CLM5 can adjust the nitrogen-use efficiency based on past environmental conditions. Therefore, the photosynthetic capacity (V_cmax25) in the model also changes dynamically. This will improve the accuracy of GPP in the area with big seasonal changes and the photosynthesis in the model might be more sensitive to warming than the model used traditional method. The plant phenology in BCC-ESM1 (BCC_AVIM2.0) is controlled by temperature (Li et al., 2019). The model adopts a method based on accumulated heat to determine when leaves grow. This is different from JSBACH, which uses chill days to break dormancy (Reick et al., 2021). Therefore, BCC-ESM1 uses heat to support growth, while JSBACH focuses on the end of cold days. This difference will lead to some divergences in predicting when the growing season starts, especially in cold regions. Our results are fully consistent with the theory. Table S3 showed that the structural feedback of BCC-ESM1 is generally the highest among all ESMs across almost all PFTs. This suggests that heat-accumulation phenology scheme of BCC-ESM1 makes the vegetation phenology more sensitive to temperature changes induced by aerosols. In summary, these models use different methods to estimate the impact of temperature on photosynthesis and this will lead to the different response of GPP to temperature.

4.3 Limitations and implications for future projections

Validation against FLUXNET observations reveals a systemic underestimation of GPP across all five models (Fig. S2). However, this underestimation is largely attributed to the scale mismatch between ground-based eddy covariance measurements and model grid pixels. Comparisons of CMIP6 ESMs against the FLUXCOM-X did not show this underestimation (Fig. S3). This indicates that ESMs can capture the global GPP magnitudes. Our attribution framework incorporated a systematic bias term to account for the bias of GPP. In addition, our framework was developed to investigate the sources of inter-model spread, rather than to evaluate the accuracy of model simulations. This also mitigates the direct impact of GPP bias on our analysis.

Although the attribution framework in this study can capture more than 50 % of inter-model GPP spreads, there are still large remaining unexplained variations. First, the framework omits higher-order non-linear climate interactions. For example, heatwaves accompanied by concurrent precipitation deficits exert an exponentially influence on VPD and vegetation stomatal conductance (Wang et al., 2025). Second, we utilized precipitation as the primary hydrological driver. However, precipitation might not adequately reflect the actual water stress on vegetation photosynthesis in some regions (Song et al., 2022). Third, the unexplained inter-model spread is highly related to unaccounted biogeochemical constraints, particularly the nutrition limitations in the model. Furthermore, it is impossible to fully disentangle internal climate variability by using single ensemble member per model. This limitation mainly occurrs during the early industrial period, when the AOD changes are small and is comparable in magnitude to the internal climate variability induced noise in the GPP difference between historical and hist-piAer simulations.

In this study, our framework identified structural and physiological mechanisms driving the inter-model spread in simulating the response of photosynthesis to aerosols. This also highlights the pathways for future model evaluation and reducing the inter-model spread. First, canopy radiative transfer module should be evaluated to make sure that the model can capture the impacts of diffuse and direct radiation accurately. Some parameters that affect canopy radiation transfer module should be evaluated and incorporated into the model. For example, Li et al. (2023) reported that the light environment within canopy was affected by the clumping index. However, many ESMs do not incorporate the clumping index and this will induce some uncertainties in simulating the canopy light environment (Fang, 2021). Second, the responses of photosynthesis to soil moisture and air temperature require rigorous validation (Gabele et al., 2026). Although some models have incorporated the acclimation of photosynthesis via various approaches, most of the approaches were not sufficient to capture the non-linear impact of air temperature and precipitation on photosynthesis. Eco-evolutionary optimality (EEO) theories offer an opportunity to re-evaluate and constrain the underlying physiological responses (Ren et al., 2025). Using these observation-supported theories to evaluate the ESMs is the key to reducing the uncertainties in simulating the response of photosynthesis to temperature and precipitation.

5 Summary

Anthropogenic aerosol loadings have increased significantly since 1850. The increased aerosols can significantly affect the terrestrial carbon cycle through reducing the total shortwave radiation, increasing the diffuse radiation fraction, and altering the temperature and precipitation. Many models were developed for simulating the effects of aerosols on regional terrestrial carbon cycle. However, divergence among models in simulating the effects of aerosols on gross primary production (GPP) still need further investigation. In this study, we investigated the differences in simulating the aerosol-induced GPP changes among the models and the driving factors using five Earth System Models (ESMs) from CMIP6, including BCC-ESM1, IPSL-CM6A-LR, MPI-ESM-1-2-HAM, NorESM2-LM, and UKESM1-0-LL. Our results indicated that all five models simulated a reduction in global GPP. However, there are large uncertainties in the magnitude and spatial distribution of these changes. Our results showed that inter-model spread was mainly caused by terrestrial ecosystem climate sensitivities, rather than atmospheric aerosol radiative and climatic effects in ESMs. Specifically, the divergence was mainly induced by the different assumptions of canopy radiative transfer and thermal acclimation. Our findings indicated that refining atmospheric aerosol optical properties alone was insufficient to reduce inter-model spread in simulating aerosol-induced GPP changes. Future efforts should be used to improve the response of photosynthesis to climatic factors.