Articles | Volume 17, issue 11
Model evaluation paper
12 Jun 2024
Model evaluation paper |  | 12 Jun 2024

Development and evaluation of the interactive Model for Air Pollution and Land Ecosystems (iMAPLE) version 1.0

Xu Yue, Hao Zhou, Chenguang Tian, Yimian Ma, Yihan Hu, Cheng Gong, Hui Zheng, and Hong Liao

Land ecosystems are important sources and sinks of atmospheric components. In turn, air pollutants affect the exchange rates of carbon and water fluxes between ecosystems and the atmosphere. However, these biogeochemical processes are usually not well presented in Earth system models, limiting the explorations of interactions between land ecosystems and air pollutants from regional to global scales. Here, we develop and validate the interactive Model for Air Pollution and Land Ecosystems (iMAPLE) by upgrading the Yale Interactive Terrestrial Biosphere Model with process-based water cycles, fire emissions, wetland methane (CH4) emissions, and trait-based ozone (O3) damage. Within iMAPLE, soil moisture and temperature are dynamically calculated based on the water and energy balance in soil layers. Fire emissions are dependent on dryness, lightning, population, and fuel load. Wetland CH4 is produced but consumed through oxidation, ebullition, diffusion, and plant-mediated transport. The trait-based scheme unifies O3 sensitivity of different plant functional types (PFTs) with the leaf mass per area. Validations show correlation coefficients (R) of 0.59–0.86 for gross primary productivity (GPP) and 0.57–0.84 for evapotranspiration (ET) across the six PFTs at 201 flux tower sites and yield an average R of 0.68 for CH4 emissions at 44 sites. Simulated soil moisture and temperature match reanalysis data with high R above 0.86 and low normalized mean biases (NMBs) within 7 %, leading to reasonable simulations of global GPP (R=0.92, NMB=1.3 %) and ET (R=0.93, NMB=-10.4 %) against satellite-based observations for 2001–2013. The model predicts an annual global area burned of 507.1 Mha, close to the observations of 475.4 Mha with a spatial R of 0.66 for 1997–2016. The wetland CH4 emissions are estimated to be 153.45 Tg [CH4] yr−1 during 2000–2014, close to the multi-model mean of 148 Tg [CH4] yr−1. The model also shows reasonable responses of GPP and ET to the changes in diffuse radiation and yields mean O3 damage of 2.9 % to global GPP. iMAPLE provides an advanced tool for studying the interactions between land ecosystems and air pollutants.

1 Introduction

As an important component on the Earth, land ecosystems regulate global carbon and water cycles. Every year, the terrestrial ecosystem assimilates ∼120 Pg (1 Pg=1015 g) of carbon from the atmosphere through vegetation photosynthesis (Beer et al., 2010). However, most of this carbon uptake returns to the atmosphere due to plant and soil respiration (Sitch et al., 2015), as well as other perturbations such as biomass burning and biogenic emissions (van der Werf et al., 2010; Carslaw et al., 2010), leading to a net carbon sink of only ∼2Pg C yr−1 during 1960–2021 (Friedlingstein et al., 2022). Meanwhile, land ecosystems affect atmospheric moisture and soil wetness through both physical (e.g., evaporation and runoff) and physiological (e.g., leaf transpiration and root hydrological uptake) processes. Observations show that transpiration accounts for 80 %–90 % of the terrestrial evapotranspiration (ET) (Jasechko et al., 2013) and makes significant contributions to land precipitation, especially over the tropical forests (Spracklen et al., 2012).

Different approaches have been applied to depict the spatiotemporal variations of ecosystem processes. The eddy covariance technique provides direct measurements of land carbon and water fluxes (Jung et al., 2011). However, the limited number and uneven distribution of ground sites result in large uncertainties in the upscaling of site-level fluxes to the global scale (Jung et al., 2020). Satellite retrieval provides a unique tool for the continuous representation of land fluxes in both space and time (Worden et al., 2021). However, most of the ecosystem variables (e.g., gross primary productivity, GPP) can only be derived using available signals from remote sensing through empirical relationships (Madani et al., 2017). As a comparison, process-based models build physical parameterizations based on field and/or laboratory experiments and validate against the available in situ and satellite-based observations (Niu et al., 2011; Castillo et al., 2012). These models can be further applied at different spatial (from site to global) and temporal (from days to centuries) scales to identify the main drivers of the changes in carbon and water fluxes (Sitch et al., 2015). For example, a total of 17 vegetation models were validated and combined to predict the land carbon fluxes in the past century (Friedlingstein et al., 2022); the ensemble mean of these models revealed a steadily increasing land carbon sink from 1960 with the dominant contribution by CO2 fertilization.

While many studies have quantified the ecosystem responses to the effects of CO2, climate, and human activities (Piao et al., 2009; Sitch et al., 2015), few have explored the interactions between air pollution and land ecosystems. Such biogeochemical processes are becoming increasingly important in the Anthropocene period with significant changes in atmospheric compositions. For example, observations found that nitrogen and phosphorus constrain the CO2 fertilization efficiency of global vegetation (Terrer et al., 2019), but such a limiting effect is ignored or underestimated in most of the current models (Wang et al., 2020). Tropospheric ozone (O3) damages plant photosynthesis and stomatal conductance, inhibiting carbon assimilation and the ET from the land surface (Sitch et al., 2007; Lombardozzi et al., 2015). Atmospheric aerosols can enhance photosynthesis through diffuse fertilization effects (Mercado et al., 2009) but meanwhile decrease photosynthesis by reducing precipitation (Yue et al., 2017). In turn, ecosystems act as both sources and sinks of atmospheric components. Biomass burning emits a large quantity of carbon dioxide, trace gases, and particulate matter, further influencing air quality (Chen et al., 2021), ecosystem functions (Yue and Unger, 2018), and global climate (Tian et al., 2022). Biogenic volatile organic compounds (BVOCs) are important precursors for both surface O3 and secondary organic aerosols (Wu et al., 2020), which can feed back to affect biogenic emissions (Yuan et al., 2016) and carbon assimilation (Rap et al., 2018). Wetland methane (CH4) emissions account for the dominant fraction of natural sources of CH4 and are projected to increase under global warming scenarios (Rosentreter et al., 2021; Zhang et al., 2017). On the other hand, stomatal uptake dominates the dry deposition of air pollutants over the vegetated land (Lin et al., 2020). Meanwhile, ET from forest results in the increase of water vapor in atmosphere (Spracklen et al., 2012), affecting the consequent rainfall and wet deposition of particles.

Currently, numerical models are in general developed separately for atmospheric chemistry and ecosystem processes. Chemical transport models are usually driven with prescribed emissions of biomass burning (Warneke et al., 2023) and wetland methane (Heimann et al., 2020), while ecosystem models often ignore the biogeochemical impacts of O3 and aerosols (Friedlingstein et al., 2022). In an earlier study, we developed and validated the Yale Interactive Terrestrial Biosphere (YIBs) model version 1.0 with a special focus on the interactions between atmospheric chemistry and land ecosystems (Yue and Unger, 2015). Thereafter, the YIBs model has been used offline to assess O3 vegetation damage (Yue et al., 2016), aerosol diffuse fertilization (Yue and Unger, 2017), and BVOC emissions (Cao et al., 2021a), as well as coupled to other models to investigate carbon–chemistry–climate interactions (Lei et al., 2020; Gong et al., 2021). The YIBs model joined the multi-model intercomparison project of TRENDY in the year 2020 and showed reasonable performance in the simulation of carbon fluxes (Friedlingstein et al., 2020). However, the YIBs model failed to predict typical hydrological variables such as ET and runoff due to missing carbon–water coupling modules. Furthermore, the model did not consider the nutrient limitation on plant photosynthesis and ignored some key exchange fluxes between the land and atmosphere.

In this study, we develop the interactive Model for Air Pollution and Land Ecosystems (iMAPLE) by coupling the process-based water cycle module from Noah-MP (Niu et al., 2011) to the carbon cycle in YIBs (Fig. 1). In addition, we update the original YIBs model with some major advances in biogeochemical processes including dynamic fire emissions, wetland CH4 emissions, nutrient limitations on photosynthesis, and trait-based O3 vegetation damage. Detailed descriptions of these updates are presented in the next section. iMAPLE is fully validated against available measurements in Sect. 3. The last section will summarize the model performance and rethink the prospective directions for future development.

Figure 1Illustration of biogeochemical processes in iMAPLE version 1.0. The carbon cycle is connected with the water cycle, wildfire emissions, biogenic volatile organic compounds (BVOCs) emissions, and wetland methane emissions and is affected by air pollutants including aerosols and ozone. The bold arrows indicate the directions of fluxes and air pollutants. The thin arrows indicate the influential pathways among different components. The dependences on key parameters are shown for some processes. Red font indicates new or updated processes in iMAPLE relative to the YIBs model. For detailed parameterizations please refer to Sect. 2.2.


2 Models and data

2.1 Main features of the YIBs model

The YIBs model is a process-based vegetation model predicting land carbon fluxes with dynamic changes in tree height, leaf area index, and carbon pools (Yue and Unger, 2015, hereafter YU2015). A total of nine plant functional types (PFTs) are considered including evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), tundra, shrubland, C3/C4 grassland, and C3/C4 cropland. At each grid box, a mixture of PFTs with each PFT fraction is used as model input, sharing the temperature or moisture information from the same soil column. Leaf photosynthesis is calculated using the well-established Michaelis–Menten enzyme kinetics scheme (Farquhar et al., 1980) and is coupled to stomatal conductance with modulations of air humidity and CO2 concentrations (Ball et al., 1987). The model applies a two-leaf approach to distinguish the irradiating states for sunlit and shading leaves and adopts an adaptive stratification for the radiative transfer processes within canopy layers (Spitters, 1986). The gross carbon assimilation is further regulated by the optimized plant phenology, which is mainly dependent on temperature and light for deciduous trees (Yue et al., 2015) but temperature and/or moisture for shrubland and grassland (YU2015). The assimilated carbon is allocated among the leaf, stem, and root to support autotrophic respiration and development, the latter of which is used to update plant height and leaf area (Cox, 2001). The input of litterfall triggers the carbon transition among 12 soil carbon pools and determines the magnitude of heterotrophic respiration with the joint effects of soil temperature, moisture, and texture (Schaefer et al., 2008). The net carbon uptake is then calculated by subtracting ecosystem respiration (plant and soil) and environmental perturbations (reforestation or deforestation) from the gross carbon assimilation (Yue et al., 2021). The YIBs model reasonably reproduces the observed spatiotemporal patterns of global carbon fluxes and makes contributions to the Global Carbon Project with long-term simulations of land carbon sink in the past century (Friedlingstein et al., 2020). The model specifically considers air pollution impacts on land ecosystems (Fig. 1), such as ozone vegetation damage (Yue and Unger, 2014) and the aerosol diffuse fertilization effect (Yue and Unger, 2017). YIBs implements two different schemes for BVOC emissions (Fig. 1), including the Model of Emissions of Gases and Aerosols from Nature (MEGAN, Guenther et al., 2012) and a photosynthesis-dependent (PS_BVOC) scheme (Unger et al., 2013).

2.2 New processes in iMAPLE

2.2.1 Process-based water cycles

The descriptions and units of all parameters used in this study are shown in Table S1 in the Supplement. We implement the hydrological module from Noah-MP into iMAPLE (Niu et al., 2011). The water budget closure is achieved by constructing water-balance equations among precipitation (P, Kgm-2s-1), evapotranspiration (ET, Kgm-2s-1), runoff, and terrestrial water storage change (ΔTWS) on each grid cell as follows:

(1) P = ET + runoff + Δ TWS .

Here, hourly P from MERRA-2 reanalyses is used as the input.

We then divide ET into three portions including plant transpiration (TRA), canopy evaporation (ECAN), and ground evaporation (EGRO):

(2) ET = TRA + ECAN + EGRO .

For vegetated grids, TRA is calculated as follows:

(3) TRA = ρ air × CP air × C tra × ( e sat - e ca ) PC ,

where ρair is air density, CPair is heat capacity of dry air, and PC is the psychrometric constant. esat is the saturated vapor pressure at the leaf temperature, eca is the vapor pressure of the canopy air, and Ctra is leaf transpiration conductance, which is calculated based on the Ball–Berry scheme of stomatal resistance (Yue and Unger, 2015). Meanwhile, ECAN is calculated as follows.


Here, Ccanopy,evap is the latent heat conductance from the wet leaf surface to canopy air. fwet is the wet fraction of canopy, which is a fraction of the maximum canopy precipitation interception capacity. EVAI is the effective vegetation area index and Rleaf,bdy is the bulk leaf boundary resistance. EGRO is calculated as follows:

(6) EGRO = C ground,evap ( e sat,ground RH - e ca ) .

Here, Cground,evap is the coefficient for latent heat at the ground, esat,ground is the saturated vapor pressure at the ground, and RH is the surface relative humidity.

Runoff includes surface (Rsrf) and subsurface (Rsub) components:

(7) runoff = R srf + R sub .

The surface runoff is calculated as follows:

(8) R srf = Q soil,srf - Q soil,in ,

where Qsoil,srf is the incident water in the soil surface and is the sum of the precipitation, snowmelt, and dewfall. Qsoil,in is the infiltration into the soil, which is derived from approximate solutions of Richards equations with considerations of the spatial variations in precipitation and infiltration capacity. Here, we assume exponential distributions of infiltration capacity in each grid cell following the approach by Schaake et al. (1996).


Here, Ic and Wd are the soil infiltration capacity of the model grid cell and the water deficit of the soil column, respectively. KΔt and Δt are the calibratable parameters and model time step. We assume free drainage processes in the soil column bottom, and thus the Rsub is calculated as follows:

(11) R sub = α slope × K 4 ,

where αslope=0.1 is the terrain slope index. K4 is the hydraulic conductivity in the bottom soil layer parameterized following the scheme in Clapp and Hornberger (1978) and is calculated using spatial soil profiles from Hengl et al. (2017).

Terrestrial water storage (TWS) is the sum of groundwater storage (Wgw), soil water content (Wsoil), and snow water equivalent (Wsnow):

(12) TWS = W gw + W snow + i = 1 N soil W soil .

Here, the soil module includes four layers (Nsoil=4) and Wsoil is calculated by the volumetric water content (Wi) as follows:

(13) W soil = ρ wat × W i × Δ Z i for  i = 1 , 2 , 3 , 4 ,

where water density (ρwat)=1000kg m−3, and ΔZi=0.1, 0.3, 0.6, and 1 m. Hourly Wi depends on variations of soil water diffusion (D) and hydraulic conductivity (K) as follows:

(14) W t = z D W z + K z .

Here, K and D are calculated following the parameterizations of Clapp–Hornberger curves (Clapp and Hornberger, 1978):


where φsat, Wsat, and Ksat are saturated soil capillary potential, volumetric water content, and hydraulic conductivity. Exponent b is an empirical constant depending on soil types. Soil moisture is calculated as the ratio of Ws to Wsat.

Soil temperature (Ts) is calculated through physical processes as follows:

(18) T s t = 1 C z K T T s z .

Here KT is soil specific heat capacity,

(19) K T = K e × K s - K dry , + K dry

where Ke, Ks, and Kdry are Kersten values as a function of soil wetness, saturated soil heat conductivity, and that under dry air conditions (Niu et al., 2011). C in Eq. (13) is the specific heat:

(20) C = W lip × C lip + W ice × C ice + ( 1 - W sat ) × C sat + ( W sat - W ) × C air .

Here, Wlip, Clip, and Wice, Cice indicate water content and heat capacity on soil water and ice. Csat and Cair are saturated and air heat capacity, which are empirical constants (Niu et al., 2011).

2.2.2 Dynamic fire emissions

We implement the active global fire parameterizations from Pechony and Shindell (2009) and Li et al. (2012) to iMAPLE. The fire emissions are determined by several key factors such as fuel flammability, natural ignitions, human activities, and fire spread. The fire count Nfire depends on flammability (Flam), fire ignition (including both natural ignition rate IN and anthropogenic ignition rate IA), and anthropogenic suppression (FNS):

(21) N fire = Flam × ( I N + I A ) × F NS .

Flam is a unitless metric representing conditions conducive to fire occurrence. It is parameterized as a function of vapor pressure deficit (VPD), precipitation (Prec), and leaf area index (LAI):

(22) Flam = VPD × e - 2 × Prec × LAI .

IN depends on the cloud-to-ground lightning, and IA can be expressed as

(23) I A = 0.03 × PD × k ( PD ) ,

where PD is population density. The empirical function of k(PD)=6.8×PD-0.6 stands for ignition potential by human activity. The fraction of non-suppressed fires FNS is derived as

(24) F NS = 0.05 + 0.95 × e - 0.05 × PD .

The burned area of a single fire (BAsingle) is typically taken to be elliptical in shape associated with length-to-breadth ratio (LB), head-to-back ratio (HB), and rate of fire spread (UP) as follows:

(25) BA single = π × UP 2 4 × LB × 1 + 1 HB 2 .

Then, LB and HB are related to changes in near-surface wind speed (U) as follows.


Meanwhile, UP is computed as the function of relative humidity (RH):

(28) UP = UP max × f RH × f θ × G ( U ) .

Here, UPmax is the maximum fire spread rate depending on PFTs, and fRH and fθ represent the dependence of fire spread on RH and on root-zone soil moisture, respectively. fθ is simply set to 0.5, and fRH is calculated as

(29) f RH = 1 , RH RH low RH up - RH RH up - RH low , RH low < RH < RH up 0 , RH RH up .

In this study, we set RHlow=30 % and RHup=70 % as the lower and upper thresholds of RH following the methods used in Li et al. (2012). If RH is higher than 70 %, natural fires will not occur or spread, and RH will no longer be a constraint for fire occurrence and spread if RH  30 %. G(U) is the limit of the fire spread:

(30) G ( U ) = LB 1 + 1 HB .

In general, the eccentricity of burned area is primarily influenced by near-surface wind speed, while the rate of fire spread is jointly regulated by near-surface wind speed and relative humidity. The shape of the fire is converted to a circular form when the near-surface wind speed reaches zero, and burning ceases to propagate once the relative humidity is above a specific threshold. The dependence of BAsingle on U and RH is shown in Fig. S1 in the Supplement.

Finally, the burned area (BA) is represented as

(31) BA = BA single × N fire .

The fire-emitted trace gases and aerosols (Emis) are calculated as

(32) Emis = BA × EF ,

where EF is the emission factors for different species (such as black carbon and organic carbon aerosols). It is important to note that the feedbacks from fire activities onto the terrestrial ecosystem have not been considered in the current version of iMAPLE due to the high complexity.

2.2.3 Wetland methane emissions

We implement process-based wetland CH4 emissions into iMAPLE. The anthropogenic sources of CH4 from Phase 6 of the Coupled Model Intercomparison Project (CMIP6,, last access: 5 June 2024) are also used as input for iMAPLE. For each soil layer, the flux of CH4 (FCH4) is calculated as the difference between production (PCH4) and consumption, which includes oxidation (OCH4), ebullition (ECH4), diffusion (DCH4), and plant-mediated transport through aerenchyma (ACH4) as follows:

(33) F CH 4 = P CH 4 - O CH 4 - E CH 4 - D CH 4 - A CH 4 .

The net methane emission to the atmosphere is the sum of ebullition, diffusion, and aerenchyma transport from the top soil layer.

The production of CH4 in soil depends on the quantity of carbon substrate and environmental conditions including soil temperature Ts, pH, and wetland inundation fraction fwetland as follows:

(34) P CH 4 = R h r f Ts f pH f wetland ,

where Rh is the heterotrophic respiration estimated at the grid cell (molCm-2s-1). r represents the release ratio of methane and carbon dioxide (Wania et al., 2010). We determine the dependence on Ts and soil pH in iMAPLE based on the parameterizations from the TRIPLEX-GHG model (Zhu et al., 2014). The impact factor of soil temperature fST can be calculated as follows (Zhang et al., 2002; Zhu et al., 2014).


Tmin, Tmax, and Topt represent the lowest, highest, and optimum temperature for the process of methane production and oxidation, respectively. In this study, Tmin=0°C, Tmax=45°C, and Topt=25°C (Zhu et al., 2014).

For the temperature dependence, the Q10 relationships are applied as follows:

(39) Q 10 = r b Q b T s - T base 10 .

Here rb is set to 3.0 and Qb is 1.33 with a base temperature (Tbase) of 25 °C (Zhu et al., 2014; Paudel et al., 2016). The inundation fraction of wetland at each cell describes the proportion of anaerobic conditions (Zhang et al., 2021). We ignore the impact of redox potential (Eh) because global observations are not available and the Eh-related processes are poorly characterized in current models (Wania et al., 2010).

The oxidation of CH4 is a series of aerobic activities related to temperature and CH4 concentrations:

(40) O CH 4 = [ CH 4 ] f Ts f CH 4 ,

where [CH4] is the methane amount in each soil layer (g C m−2 per layer). fCH4 is the CH4 concentration factor representing a Michaelis–Menten kinetic relationship:

(41) f CH 4 = [ CH 4 ] [ CH 4 ] + K CH 4 ,

where KCH4=5µmol L−1 is the half-saturation coefficient with respect to CH4 (Walter and Heimann, 2000). For temperature dependence of oxidation, the Q10 relationship with rb=2.0, Qb=1.9, and Tbase=12°C is adopted (Zhu et al., 2014; Paudel et al., 2016).

The diffusion of CH4 follows Fick's law with dependence on CH4 concentrations and the molecular diffusion coefficients of CH4 in the air (Da=0.2cm2 s−1) and water (Dw=0.00002cm2 s−1), respectively (Walter and Heimann, 2000). For each soil layer i, the diffusion coefficient Di can be calculated as follows:

(42) D i = D a × ( R sand × 0.45 + R silt × 0.2 + R clay × 0.14 ) × f tort × S poro × ( 1 - WFPS i ) + D w × WFPS i ,

where Rsand, Rsilt, and Rclay are the relative content of sand, silt, and clay in the soil, ftort=0.66 is the tortuosity coefficient, Sporo is soil porosity, and WFPS represents the pore space full of water (Zhuang et al., 2004).

The ebullition of CH4 occurs when the CH4 concentration is above the threshold of 0.5 mol CH4 m−3 (Walter et al., 2001). Since the process of ebullition occurs in a very short time, the bubbles will generate at once and all the flux will be released to the atmosphere if the concentration reaches the threshold. The plant-mediated transport of CH4 through aerenchyma is dependent on the concentration gradient of CH4 and plant-related factors (Zhu et al., 2014). The ACH4 is determined by the oxidation factor for roots and the aerenchyma factor for plants. The baseline value of the oxidation factor in roots is 0.5, with a regulatory range from 0.2 to 1.0 determined by the wetland plant types. The plant aerenchyma factor is calculated by the ratio of the plant root length density (typical value: 2.1 cm mg−1) and the root cross-sectional area (typical value: 0.0013 cm2), along with a plant root to atmosphere diffusion factor for methane which is modulated by plant type within a range of 0 to 1 (Zhang et al., 2002).

2.2.4 The downregulation on photosynthesis

We implement the downregulation parameterization from Arora et al. (2009) to indicate the nutrient limitations on leaf photosynthesis. A downregulating factor ε is calculated as a function of CO2 concentrations (C) as follows:

(43) ε ( C ) = 1 + γ gd ln ( C / C 0 ) 1 + γ g ln ( C / C 0 ) ,

where C0 is a reference CO2 concentration set to 288 ppm. The values of γgd=0.42 and γg=0.90 are derived from multiple measurements to constrain the CO2 fertilization. Then the downregulated photosynthesis is calculated by scaling the original value with the factor of ε.

2.2.5 Trait-based O3 vegetation damaging scheme

The YIBs model considers O3 vegetation damage using the flux-based scheme proposed by Sitch et al. (2007) (hereafter S2007), which determines the damaging ratio F of plant photosynthesis as follows:

(44) F = a PFT × max { f O 3 - t PFT , 0 } .

Here, fO3 denotes O3 stomatal flux (nmolm-2s-1) defined as

(45) f O 3 = [ O 3 ] r + k O 3 g p × ( 1 - F ) ,

where [O3] represents the O3 concentrations at the reference level (nmol m−3). r is the sum of boundary and aerodynamic resistance between the leaf surface and reference level (s m−1). gp is the potential stomatal conductance for H2O (m s−1). kO3=1.67 is a conversion factor of leaf resistance for O3 to that for water vapor. The level of O3 damage is then determined by the PFT-specific sensitivity aPFT and threshold tPFT, which are different among PFTs.

In iMAPLE, we implement the trait-based O3 vegetation damaging scheme to unify the inter-PFT sensitivities (Ma et al., 2023):

(46) a PFT = a LMA .

Here, a unified plant sensitivity a (nmol−1 g s) is scaled by leaf mass per area (LMA, g m−2) to derive the sensitivity of a specific PFT (aPFT). Accordingly, the damaging fraction F is modified as follows:

(47) F = a × max f O 3 LMA - t , 0 .

Here t (nmolg-1s-1) is a unified flux threshold for O3 vegetation damage. The fO3 in Eq. (45) is fed into Eq. (47) to build a quadratic equation for F. We solve the quadratic equation and select the F value within the range of [0,1]. The updated scheme considers the dilution effects of O3 dose through leaf cross-section by incorporating LMA. Plants with high LMA (e.g., ENF and EBF) usually have low sensitivities, and those with low LMA (e.g., DBF and crops) are more sensitive to O3 damage. The unified sensitivity a is set to 3.5 nmol−1 g s and threshold t is set to 0.019 nmolg-1s-1 by calibrating simulated F values with literature-based measurements (Ma et al., 2023).

2.3 Design of simulations

We perform four sensitivity experiments with iMAPLE. The baseline (BASE) simulation considers the two-way coupling between carbon and water cycles so that the prognostic soil meteorology drives canopy photosynthesis and evapotranspiration. A sensitivity run named BASE_NW is set up by turning off the water cycle in iMAPLE. In this simulation, the soil moisture and soil temperature are adopted from the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) reanalyses (Gelaro et al., 2017). The third and fourth runs turn on the O3 vegetation damage effect using either the LMA-based scheme (O3LMA) or the S2007 scheme (O3S2007). Surface hourly O3 concentrations are adopted from the chemical transport model simulations used in our previous study (Yue and Unger, 2018). For all simulations, iMAPLE is driven with hourly surface meteorology at a spatial resolution of 1° × 1° from the MERRA-2 reanalyses, including surface air temperature, air pressure, specific humidity, wind speed, precipitation, snowfall, and shortwave and longwave radiation. We run the model for the period of 1980–2021 using the initial conditions of the equilibrium soil carbon pool, tree height, and water fluxes from a spin-up run of 200 years driven with perpetual forcing for the year 1980.

iMAPLE is driven with observed CO2 concentrations from Mauna Loa (Keeling et al., 1976) and the land cover fraction of nine PFTs derived by combining satellite retrievals from both the Moderate Resolution Imaging Spectroradiometer (MODIS) (Hansen et al., 2003) and Advanced Very High Resolution Radiometer (AVHRR) (Defries et al., 2000). For fire emissions, we use Gridded Population of the World version 4 (, last access: 5 June 2024) to calculate human ignition and suppression. The lightning ignition is calculated using the flash rate from the Very High Resolution Gridded Lightning Climatology Data Collection Version 1 (, last access: 5 June 2024). For wetland CH4 emissions, we use the 2000–2020 global dataset of Wetland Area and Dynamics for Methane Modeling (WAD2M) derived from static datasets and remote sensing (Zhang et al., 2021), global soil pH from Hengl et al. (2017), and gridded soil texture from Scholes et al. (2011). For the LMA-based O3 damage scheme, we use gridded LMA from the trait-level dataset of TRY (Kattge et al., 2011) developed by extending field measurements with a random forest model (Moreno-Martínez et al., 2018).

2.4 Data for validations

We use observational datasets to validate the biogeochemical processes and related variables simulated by iMAPLE. For simulated carbon and water fluxes, site-level observations are collected from 201 sites from FLUXNET (Table S2 in the Supplement and Fig. 2). Among these sites, 95 have the ENF tree species as the major PFT and 106 are dominated by non-tree species, especially shrubland. Most (71 %) of the sites are located at the middle latitudes (30–60° N) of the Northern Hemisphere (NH), especially in the US and Europe. Compared to the previous evaluation in YU2015, we have many more sites in the tropics (22 in this study vs. 5 in YU2015), Asia (20 in this study vs. 1 in YU2015), and the Southern Hemisphere (28 in this study vs. 7 in YU2015) in this study. We also use the global gridded observations of GPP from the satellite retrievals including the solar-induced chlorophyll fluorescence (SIF) product GOSIF (Li and Xiao, 2019) and the Global Land Surface Satellite (GLASS) product (Yuan et al., 2010). The global observations of ET are adopted from the benchmark product of FLUXCOM (Jung et al., 2020) and the satellite-based GLASS product. For the dynamic fire module, we use monthly observed area burned from the Global Fire Emission Database version 4.1 with small fires (GFED4.1s) during 1997–2016 (van der Werf et al., 2010; Randerson et al., 2012). For methane emissions, we use site-level measurements of CH4 fluxes from FLUXNET-CH4 (Delwiche et al., 2021). We exclude monthly records with missing data for more than half of the days and calculate the long-term mean fluxes for the seasonal cycle. In total, we select 44 sites with at least 6 months of data available for the validations (Table S3 in the Supplement).

Figure 2Spatial distributions of 201 sites from global FLUXNET. The colors indicate various plant functional types (PFTs) including evergreen broadleaf forest (EBF, 13 sites), evergreen needleleaf forest (ENF, 57 sites), deciduous broadleaf forest (DBF, 25 sites), shrub (52 sites), grass (37 sites), and crop (17 sites). The black triangles indicate sites with at least 1 year of observations of diffuse radiation.

3 Model evaluations

3.1 Site-level evaluations

Simulated GPP shows correlation coefficients (R) of 0.59–0.86 for the six main PFTs with varied sample numbers (Fig. 3). The highest R is achieved for ENF, though the model underestimates the mean GPP magnitude by 20.62 % for this species. On average, simulated GPP is lower than observations for most PFTs. Compared to previous evaluation of the YIBs model (YU2015), iMAPLE with coupled water cycle improves the R of GPP simulations for ENF (from 0.65 to 0.86) and grassland (from 0.7 to 0.8) but worsens the predictions for other species such as EBF (from 0.65 to 0.59). The main cause of this degradation is the application of MERRA-2 reanalyses in the iMAPLE simulations instead of the site-level meteorology used in YU2015. The biases in the meteorological input may cause uncertainties in the simulation of GPP fluxes (Ma et al., 2021). In addition, the mismatch of vegetation cover and soil properties between the site location and the 1° × 1° grid in the simulation may further contribute to the modeling biases.

Figure 3Comparisons between observed and simulated monthly GPP from 201 FLUXNET sites. Each point indicates the average value of 1 month at a site. The red line represents linear regression between observations and simulations from the BASE experiment. The correlation coefficient (R), normalized mean bias, and numbers of samples (N) are shown in each panel. The comparisons are grouped into six PFTs including EBF, ENF, DBF, shrub, grass, and crop. The unit is gCm-2d-1.


Simulated ET matches observations with correlation coefficients of 0.57–0.84 at the FLUXNET sites (Fig. 4). Relatively better performance is achieved for ENF (R=0.83) and grassland (R=0.84), for which the model yields good predictions of GPP as well. In contrast, low correlations and high biases are predicted for shrubland and cropland. For the shrubland sites, different land types (e.g., closed shrublands, permanent wetlands, and woody savannas) share the same parameters in iMAPLE, resulting in the biases in depicting the site-specific carbon and water fluxes. For cropland, the prognostic phenology of grass species is applied in the model due to the missing plantation information for individual sites. Even with these deficits, iMAPLE in general captures the spatiotemporal variations of GPP and ET at most sites.

Figure 4The same as Fig. 3 but for ET. The unit is mm month−1.


We further compare the simulated wetland CH4 fluxes from the BASE experiment with observations at the FLUXNET-CH4 sites. Similar to the carbon flux sites, most of these CH4 flux sites are located in the NH (Fig. 5a). However, different from the carbon fluxes which usually range from 0 to 15 gCm-2d-1, the CH4 fluxes show a wide range across several orders of magnitude from 10−2 to 103g[CH4]m-2yr-1 (Fig. 5b). Such a large contrast requires a more realistic configuration of model parameters to distinguish the large gradient among sites. For example, US-Tw1 and US-Tw4 are two nearby sites within a distance of 1 km, where our simulations give a CH4 flux of 14.35 g[CH4]m-2yr-1 during 2011–2017. However, the average CH4 flux shows a difference of 3.7 times with 66.31 g[CH4]m-2yr-1 at US-Tw1 and 18.16 g[CH4]m-2yr-1 at US-Tw4 during 2011–2017. In the model, these two sites share the same land surface properties because they are located on the same grid. On average, simulated CH4 fluxes are correlated with observations at a moderate R of 0.68 and a normalized mean bias (NMB) of −28 %.

Figure 5(a) Spatial distribution of global FLUXNET-CH4 sites and (b) comparisons between observed and simulated monthly methane flux from the BASE experiment. Filled triangles indicate sites with at least 6 months of observations of wetland CH4 fluxes. Each point represents an average value of monthly methane emission at one site. The correlation coefficient (R), normalized mean bias, and numbers of samples (N) are shown in the right panel. The unit is g[CH4]m-2yr-1.

3.2 Grid-level evaluations

The coupling of the Noah-MP module enables the dynamic prediction of soil parameters by iMAPLE. We compare the simulated soil moisture and soil temperature from the BASE experiment with MERRA-2 reanalyses (Fig. 6). Both simulations (Fig. 6a) and observations (Fig. 6b) show low soil moisture over arid and semiarid regions with the minimum in North Africa. The model also captures the high soil moisture in tropical rainforest. However, the prediction underestimates soil moisture in boreal regions in the NH (Fig. 6c). On the global scale, simulated soil moisture matches observations with a high R of 0.86 and a low NMB of −6.9 %. These statistical metrics are further improved for the simulated soil temperature with R of 0.99 and NMB of 0.5 % against observations (Fig. 6f). The simulation reproduces the observed spatial pattern with a uniform warming bias.

Figure 6Comparisons of simulated (a) soil moisture (m3 m−3) and (d) soil temperature (K) from iMAPLE with (b, e) the MERRA-2 reanalyses. Both simulations from the BASE experiment and observations from MERRA-2 reanalyses are averaged for the period of 1980–2020. The spatial difference, correlation coefficient (R), normalized mean bias (NMB) between simulations and observations, and numbers of points (N) are shown in (c) and (f).

Driven with the prognostic soil moisture and temperature, iMAPLE predicts reasonable land carbon and water fluxes (Fig. 7). Simulated GPP (Fig. 7a) reproduces observed patterns (Fig. 7b) with high values in the tropical rainforest, moderate values in the boreal forests, and low values in the arid regions. On the global scale, our simulations yield a total GPP of 129.8 Pg C yr−1, similar to the observed amount of 125.4 Pg C yr−1. The predicted GPP is higher than observations over the tropical rainforest (Fig. 7c). However, such overestimation may instead be an indicator of biases in the ensemble observations, which are derived from the empirical models instead of direct measurements (Yuan et al., 2010; Running et al., 2004). Our site-level evaluations show that iMAPLE predicts reasonable GPP values at the EBF sites (Fig. 3). Despite this inconsistency, the model yields a high R of 0.92 and a small NMB of 1.3 % for GPP against observations on the global scale (Fig. 7c). Simulated ET (Fig. 7d) matches the observations (Fig. 7e) with high values in the tropical rainforest and secondary high values in the boreal forest. In general, the prediction is lower than observations except for the eastern US and eastern China (Fig. 7f). On average, iMAPLE shows R of 0.93 and NMB of −10.4 % in the simulation of ET compared to the ensemble of observations.

Figure 7Comparisons of simulated (a) gross primary productivity (GPP, gCm-2d-1) and (d) evapotranspiration (ET, mm month−1) with ensemble products from (b, e) observations. Simulated GPP and ET are performed by iMAPLE driven with meteorology from MERRA-2 reanalysis (BASE) during 2001–2013. Ensemble GPP products are from the average values of SIF-based GOSIF and satellite-based GLASS GPP products. Ensemble ET products include FLUXCOM and GLASS products during 2001–2013. The spatial difference, correlation coefficient (R), normalized mean bias (NMB) between simulations and observations, and numbers of points (N) are shown in (c) and (f). Only land grids with vegetation are shown in each panel, and their area-weighed values are shown in titles.

We further compare the simulated GPP with (BASE) or without (BASE_NW) a dynamic water cycle (Fig. 8). Relative to the simulations driven with MERRA-2 soil moisture and temperature, iMAPLE coupled with the Noah-MP water module predicts very similar GPP over the hotspot regions such as tropical rainforest and boreal forest (Fig. 8a). However, the coupled model predicts lower GPP for grassland in the tropics (e.g., South America and central Africa) but higher GPP in arid regions (e.g., South Africa and Australia). Since the baseline GPP is very low in arid regions, the relative changes are even larger than 100 % over those areas. These GPP differences are mainly driven by the changes in soil moisture, which increases over the arid regions with the dynamic water cycle (Fig. 6c). The reduction of soil moisture in the high latitudes of the NH shows limited impacts on the predicted GPP, likely because the boreal ecosystem is more dependent on temperature than moisture (Beer et al., 2010).

Figure 8Absolute (gCm-2d-1) and relative (%) differences of global GPP between simulations with (BASE) and without (BASE_NW) two-way carbon–water coupling processes. Simulation results are averaged for the period of 1980–2020.

3.3 Ecosystem perturbations to air pollution

Within the iMAPLE framework, the land ecosystem perturbs atmospheric components through emissions from biomass burning, wetland CH4, and BVOCs. We compare the simulated burned fraction and fire-emitted organic carbon (OC) emissions with observations from GFED4.1s (Fig. 9). The largest burned fraction is predicted over the Sahel region and the countries of Angola and Zambia, surrounding the low center of the Congo rainforest. Moderate burnings could be found in northern Australia and eastern South America. Most of these hotspots are located on grassland and shrubland in the tropics, where the high temperature and limited rainfall promote regional fire activities. The model reasonably captures the observed fire pattern with a spatial correlation of 0.66 and NMB of 6.05 % (Fig. 9c), though the model overestimates the area burned in southern Africa. The predicted fire area is used to derive biomass burning emissions of air pollutants (e.g., carbon monoxide, nitrogen oxides, black carbon, organic carbon, and sulfur dioxide) with the specific emission factors (Tian et al., 2023). Furthermore, we compare fire-emitted OC from the model with GFED4.1s. The spatial pattern of OC emissions is similar to that of burned area. The simulations yield a total of 16.8 Tg yr−1 for the global fire-emitted OC, slightly higher than the amount of 16.4 Tg yr−1 from GFED4.1s with some overestimations in tropical Africa (Fig. 9f).

Figure 9Comparisons of global burned fraction (%) and fire-emitted OC emissions (10−3kgkm-1yr-1) between (a, d) simulations and (b, e) observations. Simulations are performed using iMAPLE and observations are from GFED V4.1 fire emissions products. Both simulations from the BASE experiment and observations are averaged for the 1997–2016 period. The global total area burned is shown in (a) and (b), and total OC emissions are shown in (d) and (e). The spatial difference, correlation coefficient (R), and normalized mean biases between simulations and observations are shown in (c) and (f).

The wetland emissions of CH4 show hotspots over tropical rainforests (Fig. 10a), where the dense soil carbon provides abundant substrates for emissions and the warm climate promotes emission rates. The secondary hotspots are located in the boreal regions in the NH. This spatial pattern is very similar to the map of wetland CH4 emissions predicted by an ensemble of 13 biogeochemical models (Saunois et al., 2020). On the global scale, the total wetland emission is 153.45 Tg [CH4] yr−1 during 2000–2014, close to the average of 148 ± 25 Tg [CH4] yr−1 for 2000–2017 estimated by multiple models. As a comparison, anthropogenic sources of CH4 show a high amount in China and India due to the large emissions from fossil fuels and agriculture (Fig. 10b). On the global scale, the wetland emissions are equivalent to 45.3 % of the total anthropogenic emissions. As important factors driving CH4 emissions, heterotrophic respiration shows higher values over tropical regions and eastern China with a total amount of 73.2 Pg C yr−1 (Fig. 10c), and relatively high wetland coverage is found in boreal Asia and the Amazon (Fig. 10d).

Figure 10Global simulated CH4 emissions (g[CH4]m-2yr-1) from (a) wetland and (b) anthropogenic sources, (c) heterotrophic respiration (gCm-2d-1), and (d) fraction of wetland area. The simulations are from the BASE experiment. Anthropogenic sources are adopted from CMIP6 including the sectors of energy, agriculture, industrial, residential, shipping, solvents, and transportation. The global total emissions and heterotrophic respirations are shown in each panel. All variables are averaged for 2000–2014.

Isoprene emissions from the two schemes in iMAPLE show similar spatial distributions with hotspots over tropical rainforest (Fig. 11), where the warm climate and abundant light are favorable for biogenic emissions. Compared to the MEGAN scheme, the PS_BVOC scheme yields higher emissions in the tropical rainforest and boreal forest but lower emissions for shrubland and grassland in semiarid regions (Fig. 11c). Such differences are attributed to the varied processes as well as the emission factors. Our earlier study showed that the PS_BVOC scheme predicts stronger trends in isoprene emissions than MEGAN (Cao et al., 2021a) because the former considers both CO2 fertilization and inhibition effects, while the latter considers only the inhibition effects. On the global scale, isoprene emissions are 550 Tg yr−1 with PS_BVOC (Fig. 11a) and 611 Tg yr−1 with MEGAN (Fig. 11b). These amounts are higher than the ensemble mean of 448 Tg yr−1 from the CMIP6 models (Cao et al., 2021b) but in general within the range of 412–601 Tg yr−1 as summarized by Carslaw et al. (2010).

Figure 11Global isoprene emissions (mgCm-2d-1) from (a) MEGAN, (b) PS_BVOC schemes, and (c) their differences during 1980–2020. The simulations are from the BASE experiment. The global total emissions are shown in each panel.

3.4 Air pollution impacts on ecosystem fluxes

We assess the damaging effects of surface O3 on GPP with two schemes (O3LMA–BASE and O3S2007–BASE) (Fig. 12). Simulated GPP losses show similar patterns with high damage in the eastern US, western Europe, and eastern China, where the surface O3 level is high due to anthropogenic emissions. Limited GPP damage is predicted in the tropics though with abundant forest coverage due to the low level of O3 pollution. Compared to the S2007 scheme, predicted GPP loss is further alleviated in tropical rainforest with the LMA-based scheme because the latter scheme determines lower O3 sensitivity for evergreen trees due to their higher content of chemical resistance with larger LMA values (Ma et al., 2023). On the global scale, the average GPP loss is −2.9 % with the LMA scheme and −3.2 % with the S2007 scheme. Such damage to GPP is weaker than the estimate of −4.8 % in Ma et al. (2023) because of the differences in O3 concentrations, vegetation types, and photosynthetic parameters.

Figure 12Percentage changes in global GPP caused by ozone damage effects based on (a) LMA (O3LMA–BASE) and (b) S2007 (O3S2007–BASE) schemes. The ozone damage schemes include (a) trait leaf mass per area (LMA) based on the O3LMA experiment, (b) S2007 plant ozone sensitivity from the O3S2007 experiment, and (c) their differences.

Atmospheric aerosols cause perturbations to both direct and diffuse radiation, which have different efficiencies in enhancing plant photosynthesis. Here, we separate the diffuse (diffuse fraction >0.75) and direct (diffuse fraction <0.25) components using the observed diffuse fraction and solar radiation at six FLUXNET sites and aggregate the GPP and ET fluxes for different radiation periods at certain intervals (Fig. 13). At the six selected sites, observed GPP is higher and grows faster with more diffusive light than that under direct light conditions (Fig. 13a–f). Simulations in general reproduce such features with the comparable variability. In the earlier study, simulated diffuse fertilization efficiency for GPP (changes of GPP per unit diffuse radiation) was well validated against observations at more than 20 sites (Yue and Unger, 2018). Such amelioration of GPP suggests that moderate aerosol loading is beneficial for ecosystem carbon uptake (Yue and Unger, 2017). However, dense aerosol loading may instead weaken plant photosynthesis due to the large reduction in direct radiation.

Figure 13Observed and simulated responses of site-level (a–f) GPP and (g–l) ET to diffuse and direct radiation at the FLUXNET sites. Photosynthetically active radiation (PAR) reaching the surface is divided into diffuse (diffuse fraction >0.75) and direct (diffuse fraction <0.25) radiation at six FLUXNET sites with more than 10 years of observations. Observations (simulations) are grouped over PAR bins of 40 W m−2 with error bars (shading) indicating standard deviations of GPP and ET for each bin. Red (blue) and magenta (green) represent observed and simulated responses of GPP and ET to diffuse (direct) radiation. Units of GPP and ET are gCm-2d-1 and mm h−1, respectively.


We further evaluate the ET responses to diffuse and direct radiation from iMAPLE (Fig. 13g–l). Although ET is slightly higher in diffusive conditions, the growth rates are weaker than that of GPP. The main cause of such difference is related to the varied light dependence of ET components, which consist of canopy evaporation and transpiration. Transpiration is tightly coupled with photosynthesis and will increase by diffuse radiation at a similar rate. However, evaporation is more dependent on light quantity which will decrease with the extinction of aerosols. As a result, the weakened evaporation in part offsets the increased transpiration, leading to the smaller growth rate of ET than the responses of photosynthesis and the consequent enhancement in water use efficiency (Wang et al., 2023). iMAPLE reasonably captures the lower growth rates of ET than GPP in response to diffuse radiation at the selected sites.

4 Conclusions and discussion

We develop iMAPLE by coupling the Noah-MP water module with the YIBs vegetation model. Validations show that iMAPLE predicts a reasonable distribution of soil moisture and soil temperature. Driven with these prognostic soil conditions and meteorology from reanalyses, the model reasonably reproduces the observed spatiotemporal variations of both GPP and ET fluxes at 201 sites and on the global scale. We further update the biogeochemical processes in iMAPLE to extend the model's capability in quantifying interactions between air pollution and land ecosystems. The model reasonably predicts wetland CH4 emissions at 44 sites and yields a similar global map of CH4 emissions compared to an ensemble of 13 biogeochemical models. In addition, predicted biomass burning and biogenic emissions are consistent with either satellite retrievals or results from other models. We assess the impacts of surface O3 and aerosols on ecosystem fluxes. The LMA-based scheme links the O3 sensitivity with vegetation LMA and predicts a global map of GPP loss that is consistent with the traditional scheme using the PFT-specific sensitivity. The updated scheme effectively reduces modeling uncertainties by decreasing the number of parameters for O3 sensitivity and provides an option to apply the advanced LMA map from remote sensing. The model also reproduces the observed responses of GPP and ET to diffuse radiation with a lower growth rate for ET than GPP.

There are several limitations in the current version of iMAPLE. First, it does not include the dynamic nutrient cycle. Although we implement the downregulation from Arora et al. (2009) to constrain CO2 fertilization, this limitation is dependent only on the ambient CO2 concentrations and could not represent the heterogeneous distribution of nutrients. As a result, the model could not reveal the biogeochemical effects of nitrogen and phosphorus deposition on land ecosystems. Second, the feedback of fire activities to ecosystems is ignored. iMAPLE considers the impacts of fuel load on area burned at each modeling time step. However, these fire perturbations do not in turn change the vegetation distribution and composition. The vegetation model does not consider the competition among PFTs so that fire perturbations are not allowed to change vegetation coverage. As a result, the interactions between fire and ecosystems are underestimated in the current model framework, potentially leading to overestimations of wildfire activity due to remaining fuel loads. Third, iMAPLE does not consider the dynamic changes in wetland area for CH4 emissions. Although the Noah-MP module predicts runoff and underground water, the changes in hydrological cycles are not connected with wetland area in the model. Instead, a prescribed wetland dataset is applied to reduce the possible uncertainties but meanwhile limits the explorations of CH4 changes in the historical and future periods. Meanwhile, iMAPLE considers only dynamic soil water and temperature down to the 2 m level, which may affect the interactions between climate and the land terrestrial ecosystem, especially during drier conditions. These limitations will be the focus of model development in the next step.

iMAPLE inherits the good capability of the original YIBs model in the simulations of carbon cycle. Furthermore, iMAPLE upgrades the YIBs model with carbon–water coupling and more biogeochemical processes. With iMAPLE, we could assess the changes in carbon and water fluxes, as well as their coupling, in response to environmental perturbations (e.g., climate change, air pollution, land cover change). Meanwhile, by coupling iMAPLE with climate and/or chemical models, we could further quantify the changes in meteorology and atmospheric components in response to the biogeochemical and biogeophysical processes. For example, Lei et al. (2022) revealed the strong vegetation feedback to global surface O3 during drought periods using the YIBs model coupled to a chemical transport model. Xie et al. (2019) found a significant increase in atmospheric CO2 concentrations due to O3-induced vegetation damage using the YIBs model coupled with a regional climate–chemistry model. Gong et al. (2021) estimated a surface warming in polluted regions due to ozone–vegetation feedback using the YIBs model coupled with a global climate–chemistry model. These studies indicate that iMAPLE could be used either offline or online with other models to explore the interactions among climate, chemistry, and ecosystems.

Code availability

The code for iMAPLE version 1 is available at (Tian, 2023).

Data availability

All the validation data are available to download from the cited references or data links shown in Sect. 2.4. The simulation data of monthly output from the BASE experiment during 1980–2021 with iMAPLE are available at (Tian, 2023).


The supplement related to this article is available online at:

Author contributions

XY designed the research and wrote the paper. XY and HaZ optimized codes, performed simulations, and analyzed results. HaZ, CT, YM, YH, and CG implemented codes and collected data. HuZ helped with code implementations. All authors commented on and revised the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


The authors are grateful to the editor and two anonymous reviewers for their constructive comments that have improved this study.

Financial support

This work was jointly supported by the National Key Research and Development Program of China (grant no. 2019YFA0606802), the National Natural Science Foundation of China (grant no. 42275128), and the Natural Science Foundation of Jiangsu Province (grant no. BK20220031).

Review statement

This paper was edited by Fiona O'Connor and reviewed by two anonymous referees.


Arora, V. K., Boer, G. J., Christian, J. R., Curry, C. L., Denman, K. L., Zahariev, K., Flato, G. M., Scinocca, J. F., Merryfield, W. J., and Lee, W. G.: The Effect of Terrestrial Photosynthesis Down Regulation on the Twentieth-Century Carbon Budget Simulated with the CCCma Earth System Model, J. Climate, 22, 6066–6088,, 2009. 

Ball, J. T., Woodrow, I. E., and Berry, J. A.: A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions, in: Progress in Photosynthesis Research, edited by: Biggins, J., Nijhoff, Dordrecht, Netherlands, 221–224,, 1987. 

Beer, C., Reichstein, M., Tomelleri, E., Ciais, P., Jung, M., Carvalhais, N., Rodenbeck, C., Arain, M. A., Baldocchi, D., Bonan, G. B., Bondeau, A., Cescatti, A., Lasslop, G., Lindroth, A., Lomas, M., Luyssaert, S., Margolis, H., Oleson, K. W., Roupsard, O., Veenendaal, E., Viovy, N., Williams, C., Woodward, F. I., and Papale, D.: Terrestrial Gross Carbon Dioxide Uptake: Global Distribution and Covariation with Climate, Science, 329, 834–838,, 2010. 

Cao, Y., Yue, X., Lei, Y., Zhou, H., Liao, H., Song, Y., Bai, J., Yang, Y., Chen, L., Zhu, J., Ma, Y., and Tian, C.: Identifying the drivers of modeling uncertainties in isoprene emissions: schemes versus meteorological forcings, J. Geophys. Res., 126, e2020JD034242,, 2021a. 

Cao, Y., Yue, X., Liao, H., Yang, Y., Zhu, J., Chen, L., Tian, C., Lei, Y., Zhou, H., and Ma, Y.: Ensemble projection of global isoprene emissions by the end of 21st century using CMIP6 models, Atmos. Environ., 267, 118766,, 2021b. 

Carslaw, K. S., Boucher, O., Spracklen, D. V., Mann, G. W., Rae, J. G. L., Woodward, S., and Kulmala, M.: A review of natural aerosol interactions and feedbacks within the Earth system, Atmos. Chem. Phys., 10, 1701–1737,, 2010. 

Castillo, C. K. G., Levis, S., and Thornton, P.: Evaluation of the New CNDV Option of the Community Land Model: Effects of Dynamic Vegetation and Interactive Nitrogen on CLM4 Means and Variability, J. Climate, 25, 3702–3714,, 2012. 

Chen, G., Guo, Y., Yue, X., Tong, S., Gasparrini, A., Bell, M. L., Armstrong, B., Schwartz, J., Jaakkola, J. J. K., Zanobetti, A., Lavigne, E., Saldiva, P. H. N., Kan, H., Royé, D., Milojevic, A., Overcenco, A., Urban, A., Schneider, A., Entezari, A., Vicedo-Cabrera, A. M., Zeka, A., Tobias, A., Nunes, B., Alahmad, B., Forsberg, B., Pan, S.-C., Íñiguez, C., Ameling, C., De la Cruz Valencia, C., Åström, C., Houthuijs, D., Dung, D. V., Samoli, E., Mayvaneh, F., Sera, F., Carrasco-Escobar, G., Lei, Y., Orru, H., Kim, H., Holobaca, I.-H., Kyselý, J., Teixeira, J. P., Madureira, J., Katsouyanni, K., Hurtado-Díaz, M., Maasikmets, M., Ragettli, M. S., Hashizume, M., Stafoggia, M., Pascal, M., Scortichini, M., de Sousa Zanotti Stagliorio Coêlho, M., Ortega, N. V., Ryti, N. R. I., Scovronick, N., Matus, P., Goodman, P., Garland, R. M., Abrutzky, R., Garcia, S. O., Rao, S., Fratianni, S., Dang, T. N., Colistro, V., Huber, V., Lee, W., Seposo, X., Honda, Y., Guo, Y. L., Ye, T., Yu, W., Abramson, M. J., Samet, J. M., and Li, S.: Mortality risk attributable to wildfire-related PM2×5 pollution: a global time series study in 749 locations, The Lancet Planetary Health, 5, e579–e587,, 2021. 

Clapp, R. B. and Hornberger, G. M.: Empirical equations for some soil hydraulic properties, Water Resour. Res., 14, 601–604, 1978. 

Cox, P. M.: Description of the “TRIFFID” Dynamic Global Vegetation Model, Hadley Centre Technical Note 24, Berks, UK, (last access: 5 June 2024), 2001. 

Defries, R. S., Hansen, M. C., Townshend, J. R. G., Janetos, A. C., and Loveland, T. R.: A new global 1 km dataset of percentage tree cover derived from remote sensing, Global Change Biol., 6, 247–254,, 2000. 

Delwiche, K. B., Knox, S. H., Malhotra, A., Fluet-Chouinard, E., McNicol, G., Feron, S., Ouyang, Z., Papale, D., Trotta, C., Canfora, E., Cheah, Y.-W., Christianson, D., Alberto, Ma. C. R., Alekseychik, P., Aurela, M., Baldocchi, D., Bansal, S., Billesbach, D. P., Bohrer, G., Bracho, R., Buchmann, N., Campbell, D. I., Celis, G., Chen, J., Chen, W., Chu, H., Dalmagro, H. J., Dengel, S., Desai, A. R., Detto, M., Dolman, H., Eichelmann, E., Euskirchen, E., Famulari, D., Fuchs, K., Goeckede, M., Gogo, S., Gondwe, M. J., Goodrich, J. P., Gottschalk, P., Graham, S. L., Heimann, M., Helbig, M., Helfter, C., Hemes, K. S., Hirano, T., Hollinger, D., Hörtnagl, L., Iwata, H., Jacotot, A., Jurasinski, G., Kang, M., Kasak, K., King, J., Klatt, J., Koebsch, F., Krauss, K. W., Lai, D. Y. F., Lohila, A., Mammarella, I., Belelli Marchesini, L., Manca, G., Matthes, J. H., Maximov, T., Merbold, L., Mitra, B., Morin, T. H., Nemitz, E., Nilsson, M. B., Niu, S., Oechel, W. C., Oikawa, P. Y., Ono, K., Peichl, M., Peltola, O., Reba, M. L., Richardson, A. D., Riley, W., Runkle, B. R. K., Ryu, Y., Sachs, T., Sakabe, A., Sanchez, C. R., Schuur, E. A., Schäfer, K. V. R., Sonnentag, O., Sparks, J. P., Stuart-Haëntjens, E., Sturtevant, C., Sullivan, R. C., Szutu, D. J., Thom, J. E., Torn, M. S., Tuittila, E.-S., Turner, J., Ueyama, M., Valach, A. C., Vargas, R., Varlagin, A., Vazquez-Lule, A., Verfaillie, J. G., Vesala, T., Vourlitis, G. L., Ward, E. J., Wille, C., Wohlfahrt, G., Wong, G. X., Zhang, Z., Zona, D., Windham-Myers, L., Poulter, B., and Jackson, R. B.: FLUXNET-CH4: a global, multi-ecosystem dataset and analysis of methane seasonality from freshwater wetlands, Earth Syst. Sci. Data, 13, 3607–3689,, 2021. 

Farquhar, G. D., Caemmerer, S. V., and Berry, J. A.: A biochemical-model of photosynthetic CO2 assimilation in leaves of C-3 species, Planta, 149, 78–90,, 1980. 

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Hauck, J., Olsen, A., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Quéré, C., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S., Aragão, L. E. O. C., Arneth, A., Arora, V., Bates, N. R., Becker, M., Benoit-Cattin, A., Bittig, H. C., Bopp, L., Bultan, S., Chandra, N., Chevallier, F., Chini, L. P., Evans, W., Florentie, L., Forster, P. M., Gasser, T., Gehlen, M., Gilfillan, D., Gkritzalis, T., Gregor, L., Gruber, N., Harris, I., Hartung, K., Haverd, V., Houghton, R. A., Ilyina, T., Jain, A. K., Joetzjer, E., Kadono, K., Kato, E., Kitidis, V., Korsbakken, J. I., Landschützer, P., Lefèvre, N., Lenton, A., Lienert, S., Liu, Z., Lombardozzi, D., Marland, G., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S.-I., Niwa, Y., O'Brien, K., Ono, T., Palmer, P. I., Pierrot, D., Poulter, B., Resplandy, L., Robertson, E., Rödenbeck, C., Schwinger, J., Séférian, R., Skjelvan, I., Smith, A. J. P., Sutton, A. J., Tanhua, T., Tans, P. P., Tian, H., Tilbrook, B., van der Werf, G., Vuichard, N., Walker, A. P., Wanninkhof, R., Watson, A. J., Willis, D., Wiltshire, A. J., Yuan, W., Yue, X., and Zaehle, S.: Global Carbon Budget 2020, Earth Syst. Sci. Data, 12, 3269–3340,, 2020. 

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Gregor, L., Hauck, J., Le Quéré, C., Luijkx, I. T., Olsen, A., Peters, G. P., Peters, W., Pongratz, J., Schwingshackl, C., Sitch, S., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S. R., Alkama, R., Arneth, A., Arora, V. K., Bates, N. R., Becker, M., Bellouin, N., Bittig, H. C., Bopp, L., Chevallier, F., Chini, L. P., Cronin, M., Evans, W., Falk, S., Feely, R. A., Gasser, T., Gehlen, M., Gkritzalis, T., Gloege, L., Grassi, G., Gruber, N., Gürses, Ö., Harris, I., Hefner, M., Houghton, R. A., Hurtt, G. C., Iida, Y., Ilyina, T., Jain, A. K., Jersild, A., Kadono, K., Kato, E., Kennedy, D., Klein Goldewijk, K., Knauer, J., Korsbakken, J. I., Landschützer, P., Lefèvre, N., Lindsay, K., Liu, J., Liu, Z., Marland, G., Mayot, N., McGrath, M. J., Metzl, N., Monacci, N. M., Munro, D. R., Nakaoka, S.-I., Niwa, Y., O'Brien, K., Ono, T., Palmer, P. I., Pan, N., Pierrot, D., Pocock, K., Poulter, B., Resplandy, L., Robertson, E., Rödenbeck, C., Rodriguez, C., Rosan, T. M., Schwinger, J., Séférian, R., Shutler, J. D., Skjelvan, I., Steinhoff, T., Sun, Q., Sutton, A. J., Sweeney, C., Takao, S., Tanhua, T., Tans, P. P., Tian, X., Tian, H., Tilbrook, B., Tsujino, H., Tubiello, F., van der Werf, G. R., Walker, A. P., Wanninkhof, R., Whitehead, C., Willstrand Wranne, A., Wright, R., Yuan, W., Yue, C., Yue, X., Zaehle, S., Zeng, J., and Zheng, B.: Global Carbon Budget 2022, Earth Syst. Sci. Data, 14, 4811–4900,, 2022. 

Gelaro, R., McCarty, W., Suarez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G. K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.: The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), J. Climate, 30, 5419–5454,, 2017. 

Gong, C., Liao, H., Yue, X., Ma, Y., and Lei, Y.: Impacts of ozone-vegetation interactions on ozone pollution episodes in North China and the Yangtze River Delta, Geophys. Res. Lett., 48, e2021GL093814,, 2021. 

Guenther, A. B., Jiang, X., Heald, C. L., Sakulyanontvittaya, T., Duhl, T., Emmons, L. K., and Wang, X.: The Model of Emissions of Gases and Aerosols from Nature version 2.1 (MEGAN2.1): an extended and updated framework for modeling biogenic emissions, Geosci. Model Dev., 5, 1471–1492,, 2012. 

Hansen, M. C., DeFries, R. S., Townshend, J. R. G., Carroll, M., Dimiceli, C., and Sohlberg, R. A.: Global Percent Tree Cover at a Spatial Resolution of 500 Meters: First Results of the MODIS Vegetation Continuous Fields Algorithm, Earth Interact., 7, 1–15,<0001:GPTCAA>2.0.CO;2, 2003. 

Heimann, I., Griffiths, P. T., Warwick, N. J., Abraham, N. L., Archibald, A. T., and Pyle, J. A.: Methane Emissions in a Chemistry-Climate Model: Feedbacks and Climate Response, J. Adv. Model. Earth Sy., 12, e2019MS002019,, 2020. 

Hengl, T., de Jesus, J. M., Heuvelink, G. B. M., Gonzalez, M. R., Kilibarda, M., Blagotic, A., Shangguan, W., Wright, M. N., Geng, X. Y., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, Plos One, 12, e0169748,, 2017. 

Jasechko, S., Sharp, Z. D., Gibson, J. J., Birks, S. J., Yi, Y., and Fawcett, P. J.: Terrestrial water fluxes dominated by transpiration, Nature, 496, 347–350,, 2013. 

Jung, M., Reichstein, M., Margolis, H. A., Cescatti, A., Richardson, A. D., Arain, M. A., Arneth, A., Bernhofer, C., Bonal, D., Chen, J. Q., Gianelle, D., Gobron, N., Kiely, G., Kutsch, W., Lasslop, G., Law, B. E., Lindroth, A., Merbold, L., Montagnani, L., Moors, E. J., Papale, D., Sottocornola, M., Vaccari, F., and Williams, C.: Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and sensible heat derived from eddy covariance, satellite, and meteorological observations, J. Geophys. Res., 116, G00j07,, 2011. 

Jung, M., Schwalm, C., Migliavacca, M., Walther, S., Camps-Valls, G., Koirala, S., Anthoni, P., Besnard, S., Bodesheim, P., Carvalhais, N., Chevallier, F., Gans, F., Goll, D. S., Haverd, V., Köhler, P., Ichii, K., Jain, A. K., Liu, J., Lombardozzi, D., Nabel, J. E. M. S., Nelson, J. A., O'Sullivan, M., Pallandt, M., Papale, D., Peters, W., Pongratz, J., Rödenbeck, C., Sitch, S., Tramontana, G., Walker, A., Weber, U., and Reichstein, M.: Scaling carbon fluxes from eddy covariance sites to globe: synthesis and evaluation of the FLUXCOM approach, Biogeosciences, 17, 1343–1365,, 2020. 

Kattge, J., Diaz, S., Lavorel, S., Prentice, C., Leadley, P., Bonisch, G., Garnier, E., Westoby, M., Reich, P. B., Wright, I. J., Cornelissen, J. H. C., Violle, C., Harrison, S. P., van Bodegom, P. M., Reichstein, M., Enquist, B. J., Soudzilovskaia, N. A., Ackerly, D. D., Anand, M., Atkin, O., Bahn, M., Baker, T. R., Baldocchi, D., Bekker, R., Blanco, C. C., Blonder, B., Bond, W. J., Bradstock, R., Bunker, D. E., Casanoves, F., Cavender-Bares, J., Chambers, J. Q., Chapin, F. S., Chave, J., Coomes, D., Cornwell, W. K., Craine, J. M., Dobrin, B. H., Duarte, L., Durka, W., Elser, J., Esser, G., Estiarte, M., Fagan, W. F., Fang, J., Fernandez-Mendez, F., Fidelis, A., Finegan, B., Flores, O., Ford, H., Frank, D., Freschet, G. T., Fyllas, N. M., Gallagher, R. V., Green, W. A., Gutierrez, A. G., Hickler, T., Higgins, S. I., Hodgson, J. G., Jalili, A., Jansen, S., Joly, C. A., Kerkhoff, A. J., Kirkup, D., Kitajima, K., Kleyer, M., Klotz, S., Knops, J. M. H., Kramer, K., Kuhn, I., Kurokawa, H., Laughlin, D., Lee, T. D., Leishman, M., Lens, F., Lenz, T., Lewis, S. L., Lloyd, J., Llusia, J., Louault, F., Ma, S., Mahecha, M. D., Manning, P., Massad, T., Medlyn, B. E., Messier, J., Moles, A. T., Muller, S. C., Nadrowski, K., Naeem, S., Niinemets, U., Nollert, S., Nuske, A., Ogaya, R., Oleksyn, J., Onipchenko, V. G., Onoda, Y., Ordonez, J., Overbeck, G., Ozinga, W. A., Patino, S., Paula, S., Pausas, J. G., Penuelas, J., Phillips, O. L., Pillar, V., Poorter, H., Poorter, L., Poschlod, P., Prinzing, A., Proulx, R., Rammig, A., Reinsch, S., Reu, B., Sack, L., Salgado-Negre, B., Sardans, J., Shiodera, S., Shipley, B., Siefert, A., Sosinski, E., Soussana, J. F., Swaine, E., Swenson, N., Thompson, K., Thornton, P., Waldram, M., Weiher, E., White, M., White, S., Wright, S. J., Yguel, B., Zaehle, S., Zanne, A. E., and Wirth, C.: TRY – a global database of plant traits, Global Change Biol., 17, 2905–2935,, 2011. 

Keeling, C. D., Bacastow, R. B., Bainbridge, A. E., Ekdahl, C. A., Guenther, P. R., Waterman, L. S., and Chin, J. F. S.: Atmospheric carbon dioxide variations at Mauna Loa Observatory, Hawaii, Tellus A, 28, 538–551,, 1976. 

Lei, Y., Yue, X., Liao, H., Gong, C., and Zhang, L.: Implementation of Yale Interactive terrestrial Biosphere model v1.0 into GEOS-Chem v12.0.0: a tool for biosphere–chemistry interactions, Geosci. Model Dev., 13, 1137–1153,, 2020. 

Lei, Y., Yue, X., Liao, H., Zhang, L., Zhou, H., Tian, C., Gong, C., Ma, Y., Cao, Y., Seco, R., Karl, T., and Potosnak, M.: Global perspective of drought impacts on ozone pollution episodes, Environ. Sci. Technol., 56, 3932–3940, 2022. 

Li, F., Zeng, X. D., and Levis, S.: Corrigendum to “A process-based fire parameterization of intermediate complexity in a Dynamic Global Vegetation Model” published in Biogeosciences, 9, 2761–2780, 2012, Biogeosciences, 9, 4771–4772,, 2012. 

Li, X. and Xiao, J.: Mapping Photosynthesis Solely from Solar-Induced Chlorophyll Fluorescence: A Global, Fine-Resolution Dataset of Gross Primary Production Derived from OCO-2, Remote Sens.-Basel, 11, 2563,, 2019. 

Lin, M. Y., Horowitz, L. W., Xie, Y. Y., Paulot, F., Malyshev, S., Shevliakova, E., Finco, A., Gerosa, G., Kubistin, D., and Pilegaard, K.: Vegetation feedbacks during drought exacerbate ozone air pollution extremes in Europe, Nat. Clim. Change, 10, 444–451,, 2020. 

Lombardozzi, D., Levis, S., Bonan, G., Hess, P. G., and Sparks, J. P.: The Influence of Chronic Ozone Exposure on Global Carbon and Water Cycles, J. Climate, 28, 292–305,, 2015. 

Ma, Y., Yue, X., Zhou, H., Gong, C., Lei, Y., Tian, C., and Cao, Y.: Identifying the dominant climate-driven uncertainties in modeling gross primary productivity, Sci. Total Environ., 800, 149518,, 2021. 

Ma, Y., Yue, X., Sitch, S., Unger, N., Uddling, J., Mercado, L. M., Gong, C., Feng, Z., Yang, H., Zhou, H., Tian, C., Cao, Y., Lei, Y., Cheesman, A. W., Xu, Y., and Duran Rojas, M. C.: Implementation of trait-based ozone plant sensitivity in the Yale Interactive terrestrial Biosphere model v1.0 to assess global vegetation damage, Geosci. Model Dev., 16, 2261–2276,, 2023. 

Madani, N., Kimball, J. S., and Running, S. W.: Improving Global Gross Primary Productivity Estimates by Computing Optimum Light Use Efficiencies Using Flux Tower Data, J. Geophys. Res.-Biogeo., 122, 2939–2951,, 2017. 

Mercado, L. M., Bellouin, N., Sitch, S., Boucher, O., Huntingford, C., Wild, M., and Cox, P. M.: Impact of changes in diffuse radiation on the global land carbon sink, Nature, 458, 1014–1017,, 2009. 

Moreno-Martínez, Á., Camps-Valls, G., Kattge, J., Robinson, N., Reichstein, M., Bodegom, P. V., and Running, S. W.: Global maps of leaf traits at 3 km resolution, TRY File Archive,, 2018. 

Niu, G. Y., Yang, Z. L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., and Xia, Y. L.: The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements, J. Geophys. Res., 116, D12109,, 2011. 

Paudel, R., Mahowald, N. M., Hess, P. G. M., Meng, L., and Riley, W. J.: Attribution of changes in global wetland methane emissions from pre-industrial to present using CLM4.5-BGC, Environ. Res. Lett., 11, 034020,, 2016. 

Pechony, O. and Shindell, D. T.: Fire parameterization on a global scale, J. Geophys. Res.-Atmos., 114, D16115,, 2009. 

Piao, S. L., Ciais, P., Friedlingstein, P., de Noblet-Ducoudre, N., Cadule, P., Viovy, N., and Wang, T.: Spatiotemporal patterns of terrestrial carbon cycle during the 20th century, Global Biogeochem. Cy., 23, Gb4026,, 2009. 

Randerson, J. T., Chen, Y., van der Werf, G. R., Rogers, B. M., and Morton, D. C.: Global burned area and biomass burning emissions from small fires, J. Geophys. Res.-Biogeo., 117, G04012,, 2012. 

Rap, A., Scott, C. E., Reddington, C. L., Mercado, L., Ellis, R. J., Garraway, S., Evans, M. J., Beerling, D. J., MacKenzie, A. R., Hewitt, C. N., and Spracklen, D. V.: Enhanced global primary production by biogenic aerosol via diffuse radiation fertilization, Nat. Geosci., 11, 640–644,, 2018. 

Rosentreter, J. A., Borges, A. V., Deemer, B. R., Holgerson, M. A., Liu, S. D., Song, C. L., Melack, J., Raymond, P. A., Duarte, C. M., Allen, G. H., Olefeldt, D., Poulter, B., Battin, T. I., and Eyre, B. D.: Half of global methane emissions come from highly variable aquatic ecosystem sources, Nat. Geosci., 14, 225–230,, 2021. 

Running, S., Nemani, R., Heinsch, F., Zhao, M., Reeves, M., and Hashimoto, H.: A continuous satellite-derived measure of global terrestrial primary production, BioScience, 54, 547–560,[0547:ACSMOG]2.0.CO;2, 2004. 

Saunois, M., Stavert, A. R., Poulter, B., Bousquet, P., Canadell, J. G., Jackson, R. B., Raymond, P. A., Dlugokencky, E. J., Houweling, S., Patra, P. K., Ciais, P., Arora, V. K., Bastviken, D., Bergamaschi, P., Blake, D. R., Brailsford, G., Bruhwiler, L., Carlson, K. M., Carrol, M., Castaldi, S., Chandra, N., Crevoisier, C., Crill, P. M., Covey, K., Curry, C. L., Etiope, G., Frankenberg, C., Gedney, N., Hegglin, M. I., Höglund-Isaksson, L., Hugelius, G., Ishizawa, M., Ito, A., Janssens-Maenhout, G., Jensen, K. M., Joos, F., Kleinen, T., Krummel, P. B., Langenfelds, R. L., Laruelle, G. G., Liu, L., Machida, T., Maksyutov, S., McDonald, K. C., McNorton, J., Miller, P. A., Melton, J. R., Morino, I., Müller, J., Murguia-Flores, F., Naik, V., Niwa, Y., Noce, S., O'Doherty, S., Parker, R. J., Peng, C., Peng, S., Peters, G. P., Prigent, C., Prinn, R., Ramonet, M., Regnier, P., Riley, W. J., Rosentreter, J. A., Segers, A., Simpson, I. J., Shi, H., Smith, S. J., Steele, L. P., Thornton, B. F., Tian, H., Tohjima, Y., Tubiello, F. N., Tsuruta, A., Viovy, N., Voulgarakis, A., Weber, T. S., van Weele, M., van der Werf, G. R., Weiss, R. F., Worthy, D., Wunch, D., Yin, Y., Yoshida, Y., Zhang, W., Zhang, Z., Zhao, Y., Zheng, B., Zhu, Q., Zhu, Q., and Zhuang, Q.: The Global Methane Budget 2000–2017, Earth Syst. Sci. Data, 12, 1561–1623,, 2020. 

Schaake, J. C., Koren, V. I., Duan, Q.-Y., Mitchell, K., and Chen, F.: Simple water balance model for estimating runoff at different spatial and temporal scales, J. Geophys. Res., 101, 7461–7475,, 1996. 

Schaefer, K., Collatz, G. J., Tans, P., Denning, A. S., Baker, I., Berry, J., Prihodko, L., Suits, N., and Philpott, A.: Combined Simple Biosphere/Carnegie-Ames-Stanford Approach terrestrial carbon cycle model, J. Geophys. Res.-Biogeo., 113, G03034,, 2008. 

Scholes, R. J., Colstoun, E. B. D., Hall, F. G., Collatz, G. J., Meeson, B. W., Los, S. O., and Landis, D. R.: ISLSCP II Global Gridded Soil Characteristics, ORNL DAAC [data set], Oak Ridge, Tennessee, USA,, 2011. 

Sitch, S., Cox, P. M., Collins, W. J., and Huntingford, C.: Indirect radiative forcing of climate change through ozone effects on the land-carbon sink, Nature, 448, 791–794,, 2007. 

Sitch, S., Friedlingstein, P., Gruber, N., Jones, S. D., Murray-Tortarolo, G., Ahlström, A., Doney, S. C., Graven, H., Heinze, C., Huntingford, C., Levis, S., Levy, P. E., Lomas, M., Poulter, B., Viovy, N., Zaehle, S., Zeng, N., Arneth, A., Bonan, G., Bopp, L., Canadell, J. G., Chevallier, F., Ciais, P., Ellis, R., Gloor, M., Peylin, P., Piao, S. L., Le Quéré, C., Smith, B., Zhu, Z., and Myneni, R.: Recent trends and drivers of regional sources and sinks of carbon dioxide, Biogeosciences, 12, 653–679,, 2015. 

Spitters, C. J. T.: Separating the Diffuse and Direct Component of Global Radiation and Its Implications for Modeling Canopy Photosynthesis. 2. Calculation of Canopy Photosynthesis, Agr. Forest Meteorol., 38, 231–242,, 1986. 

Spracklen, D. V., Arnold, S. R., and Taylor, C. M.: Observations of increased tropical rainfall preceded by air passage over forests, Nature, 489, 282–285,, 2012. 

Terrer, C., Jackson, R. B., Prentice, I. C., Keenan, T. F., Kaiser, C., Vicca, S., Fisher, J. B., Reich, P. B., Stocker, B. D., Hungate, B. A., Penuelas, J., McCallum, I., Soudzilovskaia, N. A., Cernusak, L. A., Talhelm, A. F., Van Sundert, K., Piao, S. L., Newton, P. C. D., Hovenden, M. J., Blumenthal, D. M., Liu, Y. Y., Muller, C., Winter, K., Field, C. B., Viechtbauer, W., Van Lissa, C. J., Hoosbeek, M. R., Watanabe, M., Koike, T., Leshyk, V. O., Polley, H. W., and Franklin, O.: Nitrogen and phosphorus constrain the CO2 fertilization of global plant biomass, Nat. Clim. Change, 9, 684–689,, 2019. 

Tian, C.: interactive Model for Air Pollution and Land Ecosystems (iMAPLE) version 1.0, figshare [code and data set],, 2023. 

Tian, C., Yue, X., Zhu, J., Liao, H., Yang, Y., Lei, Y., Zhou, X., Zhou, H., Ma, Y., and Cao, Y.: Fire–climate interactions through the aerosol radiative effect in a global chemistry–climate–vegetation model, Atmos. Chem. Phys., 22, 12353–12366,, 2022. 

Tian, C., Yue, X., Zhu, J., Liao, H., Yang, Y., Chen, L., Zhou, X., Lei, Y., Zhou, H., and Cao, Y.: Projections of fire emissions and the consequent impacts on air quality under 1.5 °C and 2 °C global warming, Environ. Pollut., 323, 121311,, 2023. 

Unger, N., Harper, K., Zheng, Y., Kiang, N. Y., Aleinov, I., Arneth, A., Schurgers, G., Amelynck, C., Goldstein, A., Guenther, A., Heinesch, B., Hewitt, C. N., Karl, T., Laffineur, Q., Langford, B., A. McKinney, K., Misztal, P., Potosnak, M., Rinne, J., Pressley, S., Schoon, N., and Serça, D.: Photosynthesis-dependent isoprene emission from leaf to planet in a global carbon-chemistry-climate model, Atmos. Chem. Phys., 13, 10243–10269,, 2013. 

van der Werf, G. R., Randerson, J. T., Giglio, L., Collatz, G. J., Mu, M., Kasibhatla, P. S., Morton, D. C., DeFries, R. S., Jin, Y., and van Leeuwen, T. T.: Global fire emissions and the contribution of deforestation, savanna, forest, agricultural, and peat fires (1997–2009), Atmos. Chem. Phys., 10, 11707–11735,, 2010. 

Walter, B. P. and Heimann, M.: A process-based, climate-sensitive model to derive methane emissions from natural wetlands: Application to five wetland sites, sensitivity to model parameters, and climate, Global Biogeochem. Cy., 14, 745–765,, 2000. 

Walter, B. P., Heimann, M., and Matthews, E.: Modeling modern methane emissions from natural wetlands: 1. Model description and results, J. Geophys. Res., 106, 34189–34206,, 2001. 

Wang, B., Yue, X., Zhou, H., Lu, X., and Zhu, J.: Enhanced ecosystem water-use efficiency under the more diffuse radiation conditions, Global Biogeochem. Cy., 37, e2022GB007606,, 2023. 

Wang, S. H., Zhang, Y. G., Ju, W. M., Chen, J. M., Ciais, P., Cescatti, A., Sardans, J., Janssens, I. A., Wu, M. S., Berry, J. A., Campbell, E., Fernandez-Martinez, M., Alkama, R., Sftch, S., Friedlingstein, P., Smith, W. K., Yuan, W. P., He, W., Lombardozzi, D., Kautz, M., Zhu, D., Lienert, S., Kato, E., Poulter, B., Sanders, T. G. M., Kruger, I., Wang, R., Zeng, N., Tian, H. Q., Vuichard, N., Jain, A. K., Wiltshire, A., Haverd, V., Goll, D. S., and Penuelas, J.: Recent global decline of CO2 fertilization effects on vegetation photosynthesis, Science, 370, 1295–1300,, 2020. 

Wania, R., Ross, I., and Prentice, I. C.: Implementation and evaluation of a new methane model within a dynamic global vegetation model: LPJ-WHyMe v1.3.1, Geosci. Model Dev., 3, 565–584,, 2010. 

Warneke, C., Schwarz, J. P., Dibb, J., Kalashnikova, O., Frost, G., Al-Saad, J., Brown, S. S., Brewer, W. A., Soja, A., Seidel, F. C., Washenfelder, R. A., Wiggins, E. B., Moore, R. H., Anderson, B. E., Jordan, C., Yacovitch, T. I., Herndon, S. C., Liu, S., Kuwayama, T., Jaffe, D., Johnston, N., Selimovic, V., Yokelson, R., Giles, D. M., Holben, B. N., Goloub, P., Popovici, I., Trainer, M., Kumar, A., Pierce, R. B., Fahey, D., Roberts, J., Gargulinski, E. M., Peterson, D. A., Ye, X. X., Thapa, L. H., Saide, P. E., Fite, C. H., Holmes, C. D., Wang, S. Y., Coggon, M. M., Decker, Z. C. J., Stockwell, C. E., Xu, L., Gkatzelis, G., Aikin, K., Lefer, B., Kaspari, J., Griffin, D., Zeng, L. H., Weber, R., Hastings, M., Chai, J. J., Wolfe, G. M., Hanisco, T. F., Liao, J., Jost, P. C., Guo, H. Y., Jimenez, J. L., Crawford, J., and The FIREX-AQ Science Team: Fire Influence on Regional to Global Environments and Air Quality (FIREX-AQ), J. Geophys. Res., 128, e2022JD037758,, 2023. 

Worden, J., Saatchi, S., Keller, M., Bloom, A. A., Liu, J., Parazoo, N., Fisher, J. B., Bowman, K., Reager, J. T., Fahy, K., Schimel, D., Fu, R., Worden, S., Yin, Y., Gentine, P., Konings, A. G., Quetin, G. R., Williams, M., Worden, H., Shi, M. J., and Barkhordarian, A.: Satellite Observations of the Tropical Terrestrial Carbon Balance and Interactions With the Water Cycle During the 21st Century, Rev. Geophys., 59, e2020RG000711,, 2021. 

Wu, K., Yang, X. Y., Chen, D., Gu, S., Lu, Y. Q., Jiang, Q., Wang, K., Ou, Y. H., Qian, Y., Shao, P., and Lu, S. H.: Estimation of biogenic VOC emissions and their corresponding impact on ozone and secondary organic aerosol formation in China, Atmos. Res., 231, 104656,, 2020. 

Xie, X., Wang, T., Yue, X., Li, S., Zhuang, B., Wang, M., and Yang, X.: Numerical modeling of ozone damage to plants and its effects on atmospheric CO2 in China, Atmos. Environ., 217, 116970,, 2019. 

Yuan, W. P., Liu, S. G., Yu, G. R., Bonnefond, J. M., Chen, J. Q., Davis, K., Desai, A. R., Goldstein, A. H., Gianelle, D., Rossi, F., Suyker, A. E., and Verma, S. B.: Global estimates of evapotranspiration and gross primary production based on MODIS and global meteorology data, Remote Sens. Environ., 114, 1416–1431, 2010. 

Yuan, X. Y., Calatayud, V., Gao, F., Fares, S., Paoletti, E., Tian, Y., and Feng, Z. Z.: Interaction of drought and ozone exposure on isoprene emission from extensively cultivated poplar, Plant Cell Environ., 39, 2276–2287,, 2016. 

Yue, X. and Unger, N.: Ozone vegetation damage effects on gross primary productivity in the United States, Atmos. Chem. Phys., 14, 9137–9153,, 2014. 

Yue, X. and Unger, N.: The Yale Interactive terrestrial Biosphere model version 1.0: description, evaluation and implementation into NASA GISS ModelE2, Geosci. Model Dev., 8, 2399–2417,, 2015. 

Yue, X. and Unger, N.: Aerosol optical depth thresholds as a tool to assess diffuse radiation fertilization of the land carbon uptake in China, Atmos. Chem. Phys., 17, 1329–1342,, 2017. 

Yue, X. and Unger, N.: Fire air pollution reduces global terrestrial productivity, Nat. Commun., 9, 5413,, 2018. 

Yue, X., Unger, N., Keenan, T. F., Zhang, X., and Vogel, C. S.: Probing the past 30 year phenology trend of US deciduous forests, Biogeosciences, 12, 4693–4709,, 2015. 

Yue, X., Keenan, T. F., Munger, W., and Unger, N.: Limited effect of ozone reductions on the 20 year photosynthesis trend at Harvard forest, Global Change Biol., 22, 3750–3759,, 2016. 

Yue, X., Unger, N., Harper, K., Xia, X., Liao, H., Zhu, T., Xiao, J., Feng, Z., and Li, J.: Ozone and haze pollution weakens net primary productivity in China, Atmos. Chem. Phys., 17, 6073–6089,, 2017. 

Yue, X., Zhang, T., and Shao, C.: Afforestation increases ecosystem productivity and carbon storage in China during the 2000s, Agr. Forest Meteorol., 296, 108227,, 2021. 

Zhang, Y., Li, C., Trettin, C. C., Li, H., and Sun, G.: An integrated model of soil, hydrology, and vegetation for carbon dynamics in wetland ecosystems, Global Biogeochem. Cy., 16, 9-1–9-17,, 2002. 

Zhang, Z., Zimmermann, N. E., Stenke, A., Li, X., Hodson, E. L., Zhu, G. F., Huang, C. L., and Poulter, B.: Emerging role of wetland methane emissions in driving 21st century climate change, P. Natl. Acad. Sci. USA, 114, 9647–9652,, 2017. 

Zhang, Z., Fluet-Chouinard, E., Jensen, K., McDonald, K., Hugelius, G., Gumbricht, T., Carroll, M., Prigent, C., Bartsch, A., and Poulter, B.: Development of the global dataset of Wetland Area and Dynamics for Methane Modeling (WAD2M), Earth Syst. Sci. Data, 13, 2001–2023,, 2021.  

Zhu, Q., Liu, J., Peng, C., Chen, H., Fang, X., Jiang, H., Yang, G., Zhu, D., Wang, W., and Zhou, X.: Modelling methane emissions from natural wetlands by development and application of the TRIPLEX-GHG model, Geosci. Model Dev., 7, 981–999,, 2014. 

Zhuang, Q., Melillo, J. M., Kicklighter, D. W., Prinn, R. G., McGuire, A. D., Steudler, P. A., Felzer, B. S., and Hu, S.: Methane fluxes between terrestrial ecosystems and the atmosphere at northern high latitudes during the past century: A retrospective analysis with a process-based biogeochemistry model, Global Biogeochem. Cy., 18, GB3010,, 2004. 

Short summary
We develop the interactive Model for Air Pollution and Land Ecosystems (iMAPLE). The model considers the full coupling between carbon and water cycles, dynamic fire emissions, wetland methane emissions, biogenic volatile organic compound emissions, and trait-based ozone vegetation damage. Evaluations show that iMAPLE is a useful tool for the study of the interactions among climate, chemistry, and ecosystems.