Articles | Volume 17, issue 9
Model description paper
08 May 2024
Model description paper |  | 08 May 2024

Terrestrial Ecosystem Model in R (TEMIR) version 1.0: simulating ecophysiological responses of vegetation to atmospheric chemical and meteorological changes

Amos P. K. Tai, David H. Y. Yung, and Timothy Lam

The newly developed offline land ecosystem model Terrestrial Ecosystem Model in R (TEMIR) version 1.0 is described here. This version of the model simulates plant ecophysiological (e.g., photosynthetic and stomatal) responses to varying meteorological conditions and concentrations of CO2 and ground-level ozone (O3) based on prescribed meteorological and atmospheric chemical inputs from various sources. Driven by the same meteorological data used in the GEOS-Chem chemical transport model, this allows asynchronously coupled experiments with GEOS-Chem simulations with unique coherency for investigating biosphere–atmosphere chemical interactions. TEMIR agrees well with FLUXNET site-level gross primary productivity (GPP) in terms of both the diurnal and monthly cycles (correlation coefficients R2>0.85 and R2>0.8, respectively) for most plant functional types (PFTs). Grass and shrub PFTs have larger biases due to generic model representations. The model performs best when driven by local site-level meteorology rather than reanalyzed gridded meteorology. Simulation using gridded meteorology agrees well for annual GPP in seasonality and spatial distribution with a global average of 134 Pg C yr−1. Application of Monin–Obukhov similarity theory to infer canopy conditions from gridded meteorology does not improve model performance, predicting an increase of +7 % in global GPP. Present-day O3 concentrations simulated by GEOS-Chem and an O3 damage scheme at high sensitivity show a 2 % reduction in global GPP with prominent reductions of up to 15 % in eastern China and the eastern USA. Regional correlations are generally unchanged when O3 is present and biases are reduced, especially for regions with high O3 damage. An increase in atmospheric CO2 concentration of 20 ppmv from the level in 2000 to the level in 2010 modestly decreases O3 damage due to reduced stomatal uptake, consistent with ecophysiological understanding. Our work showcases the utility of this version of TEMIR for evaluating biogeophysical responses of vegetation to changes in atmospheric composition and meteorological conditions.

1 Introduction

Terrestrial vegetation, as an integral part of the global biosphere, plays many vital roles in regulating the earth system. It facilities a substantial portion of the global land–atmosphere exchange of energy, momentum, and chemical species relevant for climate and atmospheric chemistry. It is a major sink for atmospheric carbon, sequestering an estimated 123 ± 8 Pg C of carbon dioxide (CO2) from the atmosphere annually through plant photosynthesis (Beer et al., 2010; Le Quéré et al., 2015), albeit with a relatively large observation-constrained range of 119–175 Pg C yr−1.

This vegetation-mediated process of CO2 sequestration, also known as gross primary productivity (GPP), is a key regulator of climate, and forests in particular are one of the largest providers of climate services (Bonan, 2008). Even before the industrial revolution, human perturbations of natural vegetation for agriculture, timber, and other uses had significant impacts on the natural carbon cycle. About a third of the total cumulative CO2 emission to date that is due to anthropogenic land cover change could have been emitted before the time of industrialization (Pongratz et al., 2009). Over the 20th century, widespread deforestation was estimated to result in a net warming of 0.13–0.15 °C due to biogeochemical warming (via carbon emission) partly offset by biogeophysical cooling (via higher albedo) (Pongratz et al., 2010). A reversal of historical land use trends, especially in the form of afforestation as well as careful management and preservation of existing forests, has the potential to help mitigate anthropogenic climate change, but the future carbon uptake capacity of forests can be substantially altered by an array of biogeochemical feedback mechanisms as forest ecosystems respond to changing climate and atmospheric composition (Arneth et al., 2010). Various global terrestrial ecosystem models have been employed, either stand-alone or coupled within an earth system model, to estimate future carbon budgets in response to global change. A multi-model comparison estimated that over the 21st century, the terrestrial biosphere can gain 0.2–1.5 Pg C for 1 part per million by volume (ppmv) increase in CO2 due to fertilization effects but lose 10–90 Pg C per degree increase in global surface temperature as forest ecosystems experience warming and more climatic stress (Arora et al., 2013).

An emerging research interest is the interactions between the terrestrial biosphere and atmospheric chemistry and the roles of short-lived atmospheric species in modulating terrestrial ecosystem functions. On the one hand, terrestrial ecosystems facilitate the removal of air pollutants from the atmosphere via the process of dry deposition, thus providing another important service for human benefits. The consequent health benefits are substantial: 17.4 × 106 t of air pollutants equivalent to USD 6.8 billion of public health cost were removed by forests in the contiguous USA in 2010 alone (Nowak et al., 2014), which is 6 % of the estimated total health cost of USD 109 billion (EUR 145 billion) due to air pollution in the USA in 2010 (Im et al., 2018). Globally it is estimated that dry deposition onto vegetated surfaces accounts for  20 % of the loss of tropospheric O3 (Wild, 2007), which is a major air pollutant detrimental to human health. On the other hand, the depositional uptake of O3 by leaves incurs substantial damage to vegetation, interfering with ecosystem functions and terrestrial biogeochemical cycling. In the process of dry deposition, O3 diffuses via leaf openings, known as stomata, into the leaf interior, where it impairs plant physiological functions and health. Stomatal uptake itself is responsible for 30 %–90 % of the deposition sink of O3 (Felzer et al., 2007; Ainsworth et al., 2012). O3 can significantly disrupt leaf photosynthesis rates, thereby hindering plant growth and reducing forest and crop productivity (Ainsworth et al., 2012). The O3-induced global yield losses for the key staple crops (wheat, rice, maize, and soybean) for 2000 were estimated to be worth USD 11–26 billion (Van Dingenen et al., 2009; Avnery et al., 2011). For natural vegetation and forests, observed GPP reductions average  10 % but could regionally be up to 30 % (e.g., Fares et al., 2013; Proietti et al., 2016; Moura et al., 2018). Modeling studies have estimated a 2 %–12 % decrease in GPP due to present-day O3, with large reductions of more than 20 % in the midlatitude regions of North America, Europe, and East Asia (Anav et al., 2011; Yue and Unger, 2014; Lombardozzi et al., 2015). O3 damage to plants in turn alters biosphere–atmosphere exchange, with ramifications for both climate and air quality. Models have estimated a 2 %–6 % decrease in global transpiration following O3 damage. The corresponding reductions in latent heat flux can regionally enhance temperature by up to 2–3 °C and alter rainfall (Li et al., 2016; Sadiq et al., 2017; Zhu et al., 2022).

Accurate predictions of both air quality and ecosystem functions, as well as their interactions, thus require proper representation of ecophysiological processes in the terrestrial ecosystems but are obscured by a complex array of nonlinear interactions between plant physiology, O3, CO2, and meteorological drivers. Elevated CO2 enhances photosynthesis and also induces stomatal closure (reducing stomatal conductance) over various timescales, likely reflecting the adaptation of plants to improve water use efficiency (Noormets et al., 2010; Franks et al., 2013). Sanderson et al. (2007) suggested that a doubling of CO2 could worsen O3 air quality by up to +8 ppbv (parts per billion by volume) due to reduced stomatal conductance and dry deposition. O3 damage to vegetation can potentially lead to a decline in the leaf area index (LAI) and stomatal uptake, which in turn would create strong positive feedback that would further enhance surface O3 by up to +6 ppbv (Sadiq et al., 2017; Zhou et al., 2018; Zhu et al., 2022). Furthermore, higher humidity generally promotes stomatal opening, while drought conditions often inhibit it (Dermody et al., 2008; Rhea et al., 2010; Monks et al., 2015). A modeling study by Emberson et al. (2013) suggested that the extended drought in association with the 2006 European heatwave might have shut down the dry depositional sink for O3 as plants closed their stomata to prevent excessive water loss, thereby leading to a greater number of O3-related premature human deaths. To complicate the matter further, O3 damage may cause stomata to respond more sluggishly to meteorological conditions. Under certain prolonged conditions (e.g., droughts) such sluggishness of stomatal response may cause them to be more open than without O3 damage (McLaughlin et al., 2007; Sun et al., 2012; Huntingford et al., 2018). These studies highlight the importance of considering the adaptive responses of plants to changing atmospheric composition and meteorological conditions in predicting future O3 air quality and ecosystem productivity, yet most atmospheric chemistry models to date rely on semi-empirical formulations for plant-mediated processes (e.g., dry deposition) without resolving ecophysiological processes that may evolve over time. Issues may also arise when coupling atmospheric chemistry and complex ecosystem models due to inconsistent driving inputs and model requirements (Clifton et al., 2020). As interpretation of model results depends largely on the underlying physiological processes, in-depth understanding of system behaviors is crucial yet lacking (Ganzeveld and Lelieveld, 1995; Hardacre et al., 2015).

A number of studies have taken advantage of the Earth System Modeling Framework (ESMF) to dynamically link dry deposition and O3 fluxes in atmospheric chemistry models to the photosynthetic and stomatal calculations in land surface models (e.g., Ganzeveld et al., 2010; Pacifico et al., 2012; Val Martin et al., 2014; Verbeke et al., 2015; Halladay and Good, 2017; Sadiq et al., 2017; Zhu et al., 2022; Bhattarai et al., 2023). These studies largely focused on long-term averages and trends rather than variability due to climate anomalies. Simulated climate is also often sensitive to land surface changes, and any simulated responses of meteorological variables to plant ecophysiological changes can further modify O3 through a cascade of feedbacks, potentially obscuring the importance and relative contribution from individual plant-mediated pathways. Fully coupled earth system models contain an intricate network of interdependencies among climate, atmospheric chemistry, and land surface and thus may not be ideal for calibrating specific model processes against observations. Stand-alone or coupled chemical transport models and ecosystem models driven by a consistent set of prescribed “offline” meteorology from observations and reanalysis datasets would be particularly useful to improve the understanding of O3–vegetation interactions in isolation and enhance model capability in predicting air quality under climate anomalies.

The Terrestrial Ecosystem Model in R version 1.0 (TEMIR v1.0), described in Sect. 2.2, is a stand-alone, multi-parameterization model designed to simulate important canopy and ecophysiological processes that are relevant for ecosystem exchange and atmospheric chemistry, including canopy radiative physics and aerodynamics, photosynthesis, stomatal behaviors, and dry deposition of different chemical species. It is designed to be entirely consistent with the GEOS-Chem global chemical transport model (CTM) in terms of model inputs and land surface representation. Driven by a consistent set of prescribed meteorological and surface flux inputs, asynchronously coupled GEOS-Chem–TEMIR experiments can be performed globally or regionally to simulate plant ecophysiological responses to changing atmospheric composition arising from, for example, O3 pollution and rising CO2, as well as to simulate a changing climate as simulated by climate models that have already been coupled to GEOS-Chem. It can also be used with user-defined meteorological and flux inputs (especially those directly from FLUXNET observations;, last access: 21 March 2024) to perform site-level simulations for various purposes, e.g., process investigation, predictions, model validation, and optimization with different parameterization schemes. (Versions of TEMIR with active biogeochemistry and crop biophysics are under development and not within the scope of this paper.) Validation and application of TEMIR to simulate O3 dry deposition and flux-based metrics of O3 damage to crops have been presented in several previous studies (i.e., Wong et al., 2019; Tai et al., 2021; Sun et al., 2022).

Developing an ecosystem model in the R programming language is beneficial to various ends. R is an increasingly popular tool for ecological research (R Core Team, 2022), especially in population and community ecology. Lai et al. (2019) surveyed more than 60 000 peer-reviewed ecology journal articles and found that the number of studies reported using R as their primary tool in data analysis increased from  10 % in 2008 to  60 % in 2017. However, ecosystem and earth system models are often written in low-level languages, such as Fortran, because the field of ecosystem and earth system modeling has close historical ties with geoscientific research due to the importance of representing the land cover and biogeochemical cycles in climate models, which are most often written in low-level languages that are less accessible to researchers outside of the field. Having a terrestrial ecosystem model in R may help enhance the accessibility to ecosystem modeling for ecological researchers who are more familiar with R, generate a common modeling framework across population, community, and ecosystem scales, and hopefully serve as a bridge between ecological and geoscientific fields to advance interdisciplinarity. Being an entirely free and open software, as well as a highly versatile and relatively user-friendly programming language, it may also help promote open science in environmental research and education, allowing the model to be more widely used as a policy-relevant assessment tool for practitioners such as those who need to assess the carbon uptake potential of tree planting or reforestation as means to achieve carbon neutrality.

2 Model description

2.1 GEOS-Chem model description

GEOS-Chem as a global CTM is widely used in research due to its versatility in tackling a multitude of atmospheric chemistry problems. We use GEOS-Chem v12.2.0 (, Bey et al., 2001) driven by assimilated meteorological observations from the Goddard Earth Observing System (GEOS) of the NASA Global Modeling Assimilation Office (GMAO). The driving meteorological data are available in 1-hourly and 3-hourly temporal resolutions with the finest horizontal resolution of 0.25° latitude by 0.3125° longitude and 72 hybrid vertical levels extending from the surface to 0.01 hPa, provided by the GEOS-Forward Processing product (GEOS-FP). Coarser resolutions of the reanalysis product Modern-Era Retrospective Analysis for Research and Applications v2 (MERRA-2) (Gelaro et al., 2017), which is a historical dataset spanning from 1980 to the present, can also be used. GEOS-Chem is equipped with detailed O3–NOx–VOC–aerosol chemical mechanisms (NOx: nitrogen oxides; VOC: volatile organic compound) that are used for simulating atmospheric chemistry and have been validated by many studies (e.g., Bey et al., 2001; Parrington et al., 2008; Zhang et al., 2010; Zhang and Wang, 2016; Hu et al., 2018). Anthropogenic missions of many species (e.g., CO, NOx, and non-methane VOCs) can be taken from global inventories (e.g., Community Emissions Data System – CEDS; Hoesly et al., 2018) and/or regional inventories through the Harmonized Emissions Component (HEMCO) v2.1 (Keller et al., 2014; Lin et al., 2021). Biogenic emissions are calculated online with the Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.1 (Guenther et al., 2012). An alternative photosynthesis-based isoprene emission (Pacifico et al., 2011) is also available (Lam et al., 2023). In this study, the meteorological input used for global simulations of both GEOS-Chem and TEMIR is the MERRA-2 product at a resolution of 2° × 2.5° latitude by longitude. The surface O3 concentrations simulated by GEOS-Chem are used as inputs for TEMIR to simulate the corresponding vegetation responses under the consistent set of MERRA-2 meteorology.

2.2 TEMIR description

TEMIR computes biogeophysical responses of terrestrial ecosystems to changes in the atmospheric (e.g., [CO2]) and terrestrial environment. Driven by the same consistent meteorological and land surface input data as GEOS-Chem, TEMIR is designed to be highly compatible with GEOS-Chem and can be coupled asynchronously with simulated atmospheric composition (e.g., [O3]) from GEOS-Chem, which is a vital aspect for the objectives of this study.

2.2.1 Plant type representation

Plant type categories considered in TEMIR follow the convention of the Community Land Model v4.5 (CLM4.5) (Oleson et al., 2013), embedded within the Community Earth System Model (CESM) v1.2.2. The plant type categories consist of 14 natural vegetation types (including generic C3 crops) (Lawrence and Chase, 2007) and 10 rainfed or irrigated crop types (Table S1 in the Supplement), giving a total of 24 different plant function types (PFTs), as well as one land type for unvegetated land or bare ground. Each model grid cell consists of a mosaic of natural or managed PFTs and/or bare ground, where only the natural PFTs share a single soil column, allowing them theoretically to compete for soil water. Each PFT or bare ground has a prescribed present-day fractional coverage in each grid cell, derived from Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data (Lawrence and Chase, 2007) according to climatic (temperature- and precipitation-based) rules (see Table 3 of Bonan et al., 2002), as well as managed crop distribution for non-generic crops (corn, temperate and winter cereals, and soybean) (Portmann et al., 2010). Each PFT has its own characteristic structural and physiological parameters (Table S2) as detailed in Oleson et al. (2013). The parameters used to represent vegetation structure include LAI, stem area index (SAI), and canopy height (h). This version of TEMIR lacks a full carbon cycle, and thus these structural parameters are prescribed as model input data. The monthly PFT-level LAI is derived from MODIS using the deaggregation methods described in Lawrence and Chase (2007), and the PFT-level SAI is derived from the LAI using the methods of Zeng et al. (2002). The PFT-level canopy heights are prescribed following Bonan et al. (2002). Users can specify any gridded total LAI input data, whereby the PFT-specific LAI that TEMIR requires is then scaled accordingly. This version of TEMIR does not dynamically simulate PFT coverage and structural parameters, and thus competition among different plant strategies or adaptation to environmental changes, such as climate change and air pollution, is not simulated. The effects of land use and land cover change (LULCC) or changing plant type distribution due to adaptation can only be included by user-modified prescribed PFT fractional coverage or LAI data obtained externally from other models or studies. These input data can, however, be updated every simulation year to represent continuous LULCC over interannual to multi-decadal timescales. Model development for a full carbon cycle for both natural vegetation and crops (Tai et al., 2021) is actively ongoing.

2.2.2 Canopy radiative transfer

Each PFT simulated per vegetated grid cell is represented as a single “big-leaf” canopy of sunlit and shaded leaves. We implement two alternative canopy radiative transfer schemes to calculate the sunlit and shaded LAI (LAIsun and LAIsha), absorbed photosynthetically active radiation (PAR) by sunlit and shaded leaves (ϕsun and ϕsha, in W m−2), canopy light extinction coefficient (Kb), surface albedo and other radiative variables as functions of direct beam and diffuse incident PAR reaching the canopy top (Idir and Idiff, in W m−2), cosine of the solar zenith angle (μ), and other vegetation parameters. The default scheme follows the two-stream approximation of Dickinson (1983) and Sellers (1985), which considers light attenuation by both leaves and stems. The details of the scheme are described in Sects. 3.1, 3.3, and 4.1 of Oleson et al. (2013). In brief, the absorbed PAR averaged over the sunlit and shaded canopy per unit plant area (leaf plus stem area) is


where fsun/sha,dir/diff is the fraction of direct or diffuse incident radiation absorbed by the sunlit or shaded leaves and stems as calculated by the two-stream approach. The sunlit and shaded plant area index (PAI = LAI + SAI) is


and Kb is calculated following the two-stream approximation. The sunlit and shaded LAIs ultimately used to calculate canopy photosynthesis are


An alternative simplified scheme that accounts for light attenuation by leaves only following the convention of the Zhang et al. (2002) dry deposition model as modified from Norman (1982), which is also implemented in TEMIR (see Sect. 2.2.6), is implemented as follows.


Here, exponent a=0.7 and exponent b=1 for LAI < 2.5 or downwelling shortwave radiation flux S<200 W m−2; otherwise, a=0.8 and b=0.8.

2.2.3 Canopy photosynthesis and conductance

Leaf photosynthesis of both C3 and C4 plants is represented by the well-established formulation that relates to Michaelis–Menten enzyme kinetics and photosynthetic biochemical pathways (Farquhar et al., 1980; von Caemmerer and Farquhar, 1981; Collatz et al., 1991, 1992), a formulation which considers three limiting regimes as given below.

  • i.

    The Rubisco-limited photosynthesis rate (Ac, in µmol CO2 m−2 s−1) captures the rate of carbon assimilation when substrate availability or enzyme activity is the limiting factor:

    (12) A c = V cmax c i - Γ c i + K c 1 + o i K o for C 3 plants V cmax for C 4 plants ,

    where ci (in Pa) is the intercellular CO2 partial pressure; Kc and Ko are the Michaelis–Menten constants for carboxylation and oxygenation (in Pa), respectively; oi (in Pa) is the intercellular oxygen partial pressure; Γ (in Pa) is the CO2 compensation point; and Vcmax (in µmol CO2 m−2 s−1) is the maximum rate of carboxylation.

  • ii.

    The RuBP-limited photosynthetic rate (Aj, in µmol CO2 m−2 s−1) defines the photosynthesis rate, as light intensity and thus RuBP regeneration are the limiting factors:

    (13) A j = J 4 c i Γ c i + 2 Γ for C 3 plants 0.23 ϕ for C 4 plants ,

    where J (in µmol m−2 s−1) is the electron transport rate and ϕ (in W m−2) is the absorbed PAR for either sunlit (ϕsun) or shaded (ϕsha) leaves as calculated by the canopy radiative transfer model (Sect. 2.2.2). For C3 plants, J is determined by ϕ as well and is determined as the smaller of the two roots of the following quadratic equation:

    (14) Θ PSII J 2 - I PSII + J max J + I PSII J max = 0 ,

    where Jmax (in µmol m−2 s−1) is the maximum potential rate of electron transport, ΘPSII=0.7 is the curvature parameter, and IPSII (in µmol m−2 s−1) is the light utilized in electron transport by photosystem II, determined by

    (15) I PSII = 2.3 Φ PSII ϕ ,

    where ΦPSII=0.85 is the quantum yield of photosystem II.

  • iii.

    The product-limited photosynthetic rate (Ap, in µmol CO2 m−2 s−1) represents the limitation from the regeneration rate of photosynthetic phosphate compounds:

    (16) A p = 3 T p for C 3 plants k p c i P atm for C 4 plants ,

    where Tp is the triose phosphate utilization rate (in µmol m−2 s−1), Patm (in Pa) is the ambient atmospheric pressure, and kp is the initial slope of the CO2 response curve for C4 plants.

The model considers co-limitation (Collatz et al., 1991, 1992), and the leaf-level gross photosynthesis rate (A, in µmol CO2 m−2 s−1) is given by the smaller root of the following equations:

(17) Θ c j A i 2 - A c + A j A i + A c A j = 0 , Θ ip A 2 - A i + A p A + A i A p = 0 .

where Θcj=0.98 for C3 plants and Θcj=0.80 for C4 plants, and Θip=0.95 for both. The net photosynthesis rate (An, in µmol CO2 m−2 s−1) is then

(18) A n = A - R d R d = 0.015 V cmax f R d T v f V cmax T v for C 3 plants 0.025 V cmax 1 + exp s 1 T v - s 2 1 + exp s 3 T v - s 4 1 + exp s 5 T v - s 6 for C 4 plants ,

where Rd (in µmol CO2 m−2 s−1) is the dark respiration rate; s1, s3, and s5 are 0.3, 0.2, and 1.3 K, respectively; s2, s4, and s6 are 313.15, 288.15, and 328.15 K−1, respectively; Tv is leaf temperature (in degrees K); and fRdTv and fVcmaxTv are functions to adjust for variations due to temperature (Bonan et al., 2011). All of the parameters (Vcmax, Jmax, Tp, Rd, Kc, Ko, Γ, and kp) are temperature-dependent and scale with their respective PFT-specific standard values at 25 °C by different formulations. Temperature acclimation of Vcmax and Jmax from the previous 10 d, as well as day-length dependence of Vcmax, is implemented as the default option. These are all detailed in Sects. 8.2 and 8.3 of Oleson et al. (2013).

The calculation of photosynthesis rates described above is coupled with that of stomatal conductance of water (gs, in m s−1) following the formulation of Ball et al. (1987) with m and b being the slope and intercept parameters derived from empirical data:

(19) g s = α m A n e s e sat c s P atm + b ,

where gs is controlled by the leaf surface CO2 partial pressure cs (in Pa), leaf surface water vapor pressure es (in Pa), and temperature-dependent saturation vapor pressure esat (in Pa); m=9 and b=10 000 µmol m−2 s−1 for C3 plants, and m=4 and b=40 000 µmol m−2 s−1 for C3 plants; and the factor α converts the unit of conductances from µmol H2O m−2 s−1, which is more common in ecophysiology literature, to m s−1, which is common in atmospheric science literature:

(20) α = 10 - 6 R uni θ atm P atm ,

where Runi=8.314468 J K−1 mol−1 is the universal gas constant and θatm (in degrees K) is the ambient atmospheric potential temperature. An alternative stomatal conductance scheme (Medlyn et al., 2011; Franks et al., 2017) is also implemented:

(21) g s = α 1.6 1 + m VPD A n c s P atm + b ,

where VPD = 0.001(esates) (in kPa) is the vapor pressure deficit, m has PFT-specific values consistent with CLM5.0 (Sect. 9.3 of Lawrence et al., 2020), and b=100µmol m−2 s−1. Photosynthesis and stomatal conductance are further related by the diffusive flux equations for CO2 and water vapor:


where ca (in Pa) and ea (in Pa) are the canopy air CO2 partial and water vapor pressure, ei (in Pa) is the saturation vapor pressure at the leaf temperature, E (in µmol H2O m−2 s−1) is the transpiration flux, and gb (in m s−1) is the leaf boundary layer conductance:

(24) g b = C v u d l .

The photosynthesis–stomatal conductance model considers limitation arising from soil water stress. A soil water stress factor (βt) scales the photosynthesis rate and stomatal conductance, being multiplied directly to A, Rd in Eq. (18) and b in Eq. (19) or Eq. (21) to account for soil water stress (Porporato et al., 2001; Verhoef and Egea, 2014). To compute βt, we consider a two-layer soil model consisting of a topsoil layer (0–5 cm) and a root zone beneath the top soil (5–100 cm), consistent with and constrained by the input soil moisture and model structure of MERRA-2. First, the soil matric potential in each layer i, ψi (in mm), that represents water availability in ecophysiological terms is evaluated as a function of soil wetness (si) and soil type:

(25) ψ i = ψ sat , i s i - B i ,

where ψsat,i and Bi refer to the saturated soil matric potential and soil water characteristic parameter, respectively, both depending on soil texture. A wilting factor, wi, is formulated as a function of ψi as well as ψc and ψo (Table S2), which refer to the matric potential at which stomatal closure and stomatal opening occur to the full extent, respectively:

(26) w i = 1 for ψ i > ψ o ψ c - ψ i ψ c - ψ o for ψ c ψ i ψ o 0 for ψ i < ψ c .

The function βt is then the average of the wilting factors weighted by the PFT-specific root fraction (ri) in each layer:

(27) β t = i = 1 2 w i r i .

A single-layer bulk soil formulation considering only the root zone (0–100 cm) is also implemented but is found to be inferior to the two-layer formulation in terms of reproducing observed GPP in semiarid locations (Lam and Tai, 2020).

The above equations calculate photosynthesis and conductance at the leaf level only, and appropriate scaling to account for vertical variation in leaf nitrogen content, light attenuation, and sunlit vs. shaded leaves is needed to obtain the canopy-level photosynthesis (i.e., GPP) and conductance. This is done by scaling Vcmax and other parameters as follows:


where Vcmax, sun and Vcmax, sha are the canopy-averaged leaf-level values for Vcmax for sunlit and shaded leaves, respectively, which are used to compute leaf-level photosynthesis and stomatal conductance for sunlit and shaded leaves separately (Asun, Asha, gssun, and gssha) from the equations above. Kn=0.30 is the canopy decay coefficient for nitrogen, calculated and calibrated to match an explicit multi-layer canopy (Bonan et al., 2012; Oleson et al., 2013). Vcmax, top is the PFT-specific value for Vcmax for the top of the canopy. Kb, LAIsun, and LAIsha are computed from the canopy radiative transfer model (Sect. 2.2.2). Other parameters, i.e., Jmax, Tp, kp, and Rd, scale similarly.

Canopy photosynthesis rate (i.e., GPP, in µmol CO2 m−2 s−1) and canopy conductance (gcan, in m s−1) per unit land area of a given PFT are then


The grid-cell-averaged values are obtained by weighting the PFT-level values by the PFT fractional coverage of the grid cell.

2.2.4 Canopy and surface layer aerodynamics

The canopy photosynthesis and conductance calculations above require micrometeorological variables (e.g., temperature T and specific humidity q) of the canopy air as inputs. The default approach for global and regional gridded simulations is to use reanalyzed meteorological variables at 2 m above the zero-plane displacement height (i.e., T2 m and q2 m) as the proxies for canopy air conditions. The default approach for site simulations is to directly use the measured micrometeorological variables regardless of the measurement height. We also implement an option to infer canopy air conditions from micrometeorological variables at any reference height (zref, in m) above the zero-plane displacement height (e.g., zref=10 m in the MERRA-2 reanalysis product) based on Monin–Obukhov similarity theory (MOST) (Monin and Obukhov, 1954), which relates the stability of the surface layer to the generation and suppression of turbulence through the Obukhov length:

(32) L Obuk = - ρ c p u 3 θ k g H ,

where ρ (in kg m−3) is the density of moist air (we use the value at zref), θ (in degrees K) is the potential temperature (we use T2 m as a proxy), H (in W m−2) is the sensible heat flux, cp (in J kg−1 K−1) is the heat capacity of air at constant pressure, k=0.4 is the von Kármán constant, and g=9.80616 m s−2 is the gravitational acceleration. The friction velocity u (in m s−1) is either provided as input or inferred iteratively from the wind speed at the reference height uref (in m s−1):

(33) u = k u ref ln z ref + z 0 m z 0 m - ψ m z ref + z 0 m L Obuk + ψ m z 0 m L Obuk - 1 ,

where z0m (in m) is the roughness length for momentum and the function ψm(x) for momentum flux follows the formulation of Zeng et al. (1998), consistent with the implementation in CLM4.5. The aerodynamic conductance (gah, in m s−1) for heat, water vapor, and other chemical species (e.g., ozone) between the reference height zref and the surface (treated as the zero-displacement height and where the canopy air is) is then

(34) g ah = k u ln z ref + z 0 m z 0 h - ψ h z ref + z 0 m L Obuk + ψ h z 0 h L Obuk - 1 ,

where z0h (m) is the roughness length for heat, water vapor, and other chemical species, and the function ψh(x) for heat and other material fluxes also follows the formulation of Zeng et al. (1998). Canopy air potential temperature (θa, in degrees K) and specific humidity (qa, in kg kg−1) can then be inferred as


where θref (in degrees K) and qref (in kg kg−1) are the potential temperature and specific heat capacity at zref, and E (in kg m−2 s−1) is the evapotranspiration flux.

For the computation of ozone damage and dry deposition fluxes (Sects. 2.2.5 and 2.2.6), gah is also needed and computed either using the default formulation above or an alternative formulation that is consistent with the default dry deposition scheme in GEOS-Chem (Wesely, 1989; Wang et al., 1998).

2.2.5 Ozone damage

Two ozone damage schemes are implemented in TEMIR, which considers the responses of vegetation in terms of photosynthesis and stomatal conductance. The first O3 damage scheme follows Sitch et al. (2007) and considers two levels of O3 sensitivity (high and low) for each of the five major plant groups, namely, “broadleaf”, “needleleaf”, “shrub”, “C3 grass”, and “C4 grass”, as defined by Karlsson et al. (2004) and Pleijel et al. (2004). These groups are mapped to the default TEMIR PFTs accordingly. The scheme represents O3 damage by an O3 impact factor (f) that is dependent on the instantaneous stomatal O3 flux into the leaf interior:

(37) f = 1 - a max O 3 g ah - 1 + g b - 1 + k O 3 g s - 1 - F crit , 0 ,

where [O3] (in nmol m−3) is the O3 concentration observed or of the lowest atmospheric model layer; the aerodynamic, leaf boundary layer, and stomatal conductances are calculated using the formulations in the previous sections; kO3=1.67 as defined by Sitch et al. (2007) is the ratio of the leaf resistance for O3 to that for water vapor; Fcrit represents a critical threshold accounting for O3 tolerance, below which instantaneous O3 exposure does not affect photosynthesis, and Fcrit=1.6 nmol m−2 s−1 for woody PFTs and Fcrit=5 nmol m−2 s−1 for grass PFTs; and the O3 sensitivity parameter a (in nmol−1 m2 s) is specific to the plant group and to the two levels of O3 sensitivity. Factor f is multiplied directly to the net photosynthesis rate An to represent O3 damage, which then indirectly affects gs via the coupling between An and gs. Since the calculation of f requires gs, the three variables, i.e., f, An, and gs, need to be solved together by numerical iterations. Noniterative methods give insignificant differences in performance.

The second scheme follows Lombardozzi et al. (2012, 2015) and considers three O3 sensitivity levels (high, average, and low) for each of the three major plant groups, namely, broadleaf, needleleaf, and grasses and crops, which are mapped to the default TEMIR PFTs accordingly. Unlike the Sitch et al. (2007) scheme, O3 damage is characterized by the cumulative uptake of O3 (CUO, in mmol m−2) instead of instantaneous O3 uptake, which is parameterized as the sum of the instantaneous stomatal O3 flux over the lifetime of the leaf:

(38) CUO = 10 - 6 O 3 g ah - 1 + g b - 1 + k O 3 g s - 1 - F crit Δ t ,

where Δt (in s) is the model time step, the critical threshold to account for O3 tolerance is set to Fcrit=0.8 nmol m−2 s−1, and CUO is only calculated when the LAI of the PFT concerned is larger than 0.5 to avoid unrealistically high CUO (Lombardozzi et al., 2012). Another important difference from Sitch et al. (2007) is that O3 damage alters photosynthesis and stomatal conductance separately using two different sets of O3 impact factors, fp and fc, respectively:


where the intercepts bp and bc, as well as the slopes ap and ac, are determined empirically for the three plant groups (Lombardozzi et al., 2015). Factors fp and fc are multiplied separately to An and gs, respectively, after the iterative calculation of An and gs.

2.2.6 Dry deposition

We implement two major dry deposition schemes: the Zhang et al. (2003) scheme used in several Canadian and American air quality models and the Wesely (1989) scheme widely used in many chemical transport models including WRF-Chem and GEOS-Chem. In each of the two schemes, the default stomatal conductance scheme is a semi-empirical formulation that is not coupled to plant ecophysiology. The default canopy radiative transfer and aerodynamic conductance also follow formulations that are different from the default TEMIR schemes described above. We implement options such that ecophysiology-based stomatal conductance (gs) computed from the photosynthesis model above (Sect. 2.2.3), as well as canopy radiative transfer (ϕsun, ϕsha, LAIsun, and LAIsha; Sect. 2.2.2) and aerodynamic conductance (gah; Sect. 2.2.4), can be used to replace the default options in the dry deposition schemes. A full evaluation of dry deposition velocities and fluxes computed by TEMIR using different combinations of schemes against O3 flux observations has been conducted by Sun et al. (2022). Tai et al. (2021) also evaluated how the dry depositional fluxes of O3 can affect global crop yields by integrating TEMIR with the Deposition of O3 for Stomatal Exchange (DO3SE) model (, last access: 31 March 2024).

3 Observational datasets for model evaluation

3.1 Site-level dataset

Site-level comparison utilizes the eddy covariance measurements from the flux tower sites of the FLUXNET network (, last access: 21 November 2018). The latest released dataset, FLUXNET2015, contains half-hourly or hourly measurements of carbon fluxes and various meteorological observations. Each site is classified to contain one PFT that follows the classes described in the International Geosphere Biosphere Program Data and Information System (IGBP-DIS) DISCover land cover dataset (Loveland and Belward, 1997).

Data after 2009 from FLUXNET are used for model validation, taking into account data quality using the FLUXNET quality flags for each meteorological variable measured. TEMIR follows the PFT classes of CLM4.5, which are defined differently to the IGBP-DIS scheme, so additional identification and matching of the PFT classes are performed based on the forest composition information provided as far as possible in the FLUXNET database, as shown in Table S1. Overall, five sites have mismatched PFTs where the TEMIR outputs do not contain the corresponding PFT classes specified by FLUXNET. Thus, a total of 49 sites were used for comparison as listed in Table S3 whereby most are in the Northern Hemisphere. Of these sites, 14 are evergreen needleleaf forests (ENF), 2 are evergreen broadleaf forests (EBF), and 6 are deciduous broadleaf forests (DBF), which together account for almost half of the total number of sites. The rest of the 27 sites are 1 open shrubland (OSH), 16 grasslands (GRA), and 10 croplands (CRO). The GPP products provided by FLUXNET are derived by partitioning from the net ecosystem exchange (NEE) that is directly measured by the flux towers, using two partitioning methods: one utilizing daytime data only as described by Lasslop et al. (2010) and the other utilizing the nighttime approach according to Reichstein et al. (2005). A major difference in the GPP products partitioned using these two methods is that the nighttime approach for some extreme events, such as droughts, gives negative GPP, which is not physical. The NEE of the dataset is also produced using two methods that consider constant vs. variable friction velocity thresholds to filter NEE accordingly. The FLUXNET GPP product used in our comparison is the mean GPP product partitioned using the daytime method with the variable friction velocity threshold denoted as GPP_DT_VUT_MEAN in the FLUXNET database.

3.2 Global-scale dataset

We use the solar-induced chlorophyll-fluorescence-inferred (SIF-inferred) GPP product derived by Li and Xiao (2019) for global validation. SIF-based GPP products rely on the correlation of SIF with GPP where there is a good representation (e.g., Shekhar et al., 2022). The global GPP dataset derived by Li and Xiao (2019) is based on the global SIF product from Orbiting Carbon Observatory-2 (OCO-2), namely GOSIF, and on the relationships between SIF and site-level observed GPP (Li et al., 2018). The resulting dataset has a spatial resolution of 0.05° and a monthly temporal resolution. For model validation, GOSIF GPP is regridded to a resolution of 2° × 2.5° and grid cells with zero GPP are excluded.

4 Model and simulation setup

4.1 Site-level simulations with TEMIR

For each of the 49 FLUXNET sites (Table S3) beginning from 2009, simulations shown in Table 1 are conducted for model–observation comparison, with ambient CO2 concentration kept constant at 390 ppmv. The first set is to use the default surface meteorological fields prescribed from MERRA-2 at 2° × 2.5° horizontal resolution consistent with the GEOS-Chem simulations described above. We also test turning on and off the option of using MOST to infer in-canopy micrometeorological variables from the prescribed meteorology at 2 m above displacement height, as described in Sect. 2.2.4. The second set is to use direct micrometeorological measurements available from the FLUXNET site towers as the driving meteorology and the default MERRA-2 meteorology used only to replace any missing or low-quality data. The results simulated are most relevant for evaluating model performance in reproducing the observed diurnal and seasonal cycle of GPP from each FLUXNET site.

Table 1TEMIR simulation settings for selected FLUXNET sites.

Download Print Version | Download XLSX

4.2 Global simulations with TEMIR

Global simulations from 2010 to 2015 are conducted under the same general setup as the site-level simulations, with ambient CO2 concentration fixed at 390 ppmv and driven by 2° × 2.5° MERRA-2 surface meteorology. A full list of MERRA-2 variables required for running gridded simulations of TEMIR is shown in Table S5. Functionalities of the model are tested for performance, as shown in Table 2, with the MOST option to infer in-canopy conditions and the Sitch O3 damage scheme with low and high sensitivity to assess the global O3 impact on vegetation. We use GEOS-Chem (Sect. 2.1) to simulate tropospheric O3, starting from 2009 to 2015 whereby the first year of simulation is considered as spin-up. The simulation uses the comprehensive chemistry scheme “tropchem”, which has tropospheric O3–NOx–VOC–aerosol chemistry accompanied by the default emission inventories (i.e., anthropogenic emissions from CEDS, Hoesly et al., 2018; and biogenic emissions from MEGAN, Guenther et al., 2012). The reanalyzed meteorological fields used are from MERRA-2 supplied by GMAO at 2° × 2.5° horizontal resolution, which are the identical meteorological inputs used in TEMIR. Simulated O3 concentrations at the lowest surface level are then fed into TEMIR as inputs accordingly. The results produced from these simulations (Table 2) are used to validate the spatial variability in seasonal and annual averages across the whole world against the GOSIF GPP dataset (Sect. 3.2). To retain representative results from all simulations (Table 2), grid cells with LAI < 0.5 are excluded. The transient LAI dataset for the simulation is from MODIS satellite data (Lawrence and Chase, 2007) and assimilated by Yuan et al. (2011). PFT distributional and structural data are regridded to the resolution of MERRA-2 data for TEMIR simulations. We also note that as the model mechanisms are essentially resolution-independent, the model can be straightforwardly modified to conduct simulations at higher resolutions as long as the corresponding input data are provided.

Table 2TEMIR simulation settings for global simulation for 2010–2015.

Download Print Version | Download XLSX

4.2.1 Global CO2–O3 factorial simulations with TEMIR

We perform factorial simulations (Table 3) to investigate the effects of CO2 fertilization, O3 damage, and their interactions on global primary productivity as an example to showcase the utility of the model. Global O3 surface concentrations of the year 2000 are simulated using GEOS-Chem (Sect. 2.1) with 1999 used as spin-up, and other settings are as described in Sect. 4.2. CO2 concentrations are changed in TEMIR as required for each simulation (Table 3) with reference to concentrations for 2000 and 2010 (Dlugokencky and Tans, 2022). The Sitch O3 damage scheme with high sensitivity is used for TEMIR simulations (Table 3) when global surface O3 concentrations are used as inputs; otherwise, no O3 damage scheme is used. Meteorological fields for 2000 in MERRA-2 are used for all GEOS-Chem and TEMIR simulations. The LAI is fixed at the year 2000 from MODIS for all GEOS-Chem and TEMIR simulations.

Table 3Factorial simulation settings. “N/A” for O3 year indicates that no O3 damage scheme is used.

Download Print Version | Download XLSX

5 Model evaluation

Statistics used for model validation are the adjusted coefficient of determination R2, the modified Nash–Sutcliffe model efficiency coefficient N, and the normalized mean bias B. R2 is a commonly used metric with a range of 0 to 1 (or 0 %–100 %) that represents the fraction of variability in observations that can be replicated by the model, whereby 1 indicates perfect correlation and 0 indicates no correlation. N addresses the sensitivity issues of R2 documented by Legates and McCabe (1999). With values from negative infinity to 1, it is a measure of the suitability of the model as a predictor instead of using the mean of the observations. When N=1, it indicates that the model perfectly replicates observations, and no preference is observed between the model and the mean of the observations as a predictor. Negative values in turn signify the incapability of the model in predicting system behaviors. B gives the relative difference of the magnitude of model results from the observations. The equations to compute these statistics are shown in Table 4.

Table 4Statistical metrics for TEMIR validation, where M and O respectively represent the simulated dataset and observational dataset, each containing n data points. M¯ and O¯ represent the means of the datasets in question.

Download Print Version | Download XLSX

5.1 Site-level validation

5.1.1 Validation on summer diurnal cycle

Figure 1 shows the statistical metrics (Table 4) used to perform a model–observation comparison for the diurnal GPP cycle calculated for each FLUXNET site taken from the second summer month of 2012, which corresponds to the month of July for sites in the Northern Hemisphere and January in the Southern Hemisphere (Table S3). Comparing the simulations with FLUXNET local meteorology and MERRA-2 meteorology (Fig. 1), statistical metrics generally do not differ substantially, but some significant differences can exist for some sites (e.g., CH-Cha and CH-Dav; see Fig. 3) where results from simulations with FLUXNET local meteorology show higher correlations. For simulations solely driven by MERRA-2 meteorology, inferring in-canopy meteorology using Monin–Obukhov similarity theory (Fig. 1c) gives insignificant differences for all sites.

Figure 1Statistical metrics (see Table 4) for model–observation comparison for the diurnal gross primary productivity from simulations (a) TEMIR_FLUX, (b) TEMIR_MO_off, and (c) TEMIR_ MO_on as described in Table 1.


Figure 2 shows the statistical metrics (Table 4) for model–observation comparison for the simulations using FLUXNET local meteorology as the driving meteorological input for different PFTs. The average correlations per PFT between observed and simulated GPP (Fig. 2) are high (R2>0.88), except for open shrubland (R2≈0.7). B shows a large variability due to various limitations of the model for each PFT. For forest sites, B generally has a smaller range and lower absolute mean values in comparison with the other PFTs, showcasing the better performance of TEMIR for forests. N shows similar behaviors, namely, the prediction for forest sites is satisfactory with mean values of N larger than 0.55, whereas N is mostly negative for grasslands, croplands, and open shrublands. All plots of diurnal cycles are shown in the Supplement, with relevant figures also included in the following discussion.

Figure 2Statistical metrics (see Table 4) for model–observation comparison for the diurnal gross primary productivity from the TEMIR_FLUX simulation (see Table 1) for each plant functional type listed in Table S1: (a) evergreen needleleaf forest (ENF), (b) evergreen broadleaf forest (EBF), (c) deciduous broadleaf forest (DBF), (d) open shrubland (OSH), (e) grassland (GRA), and (f) cropland (CRO).


We find that the correlations are above 90 % (Fig. 2e) for all grassland sites (e.g., AU-How in Fig. 3a, IT-MBo in Fig. 3b, and CZ-wet in Fig. 3c). Yet the spread of B is large, where we see absolute B values greater than +0.6 for 6 of the 16 sites, and the rest with absolute B less than 0.3. Overestimation of CH-Cha (Fig. 3d) is similar under FLUXNET meteorology, which is likely due to disturbances from intensive site management (i.e., cutting, slurry application, and grazing; Imer et al., 2013; Merbold et al., 2014); this is a shortcoming of simplistic model representation for crops. A possible explanation for the high B values is the fire-prone nature of these sites (i.e., AU-Stp in Fig. 3e) (Beringer et al., 2007, 2011; Hutley et al., 2011; Haverd et al., 2013) whereby the model is incapable of resolving such complexities as turnover and local disturbances. Another cause of overestimation is the simplistic and generic PFT classification for such biomes, which are usually sparsely populated yet with much diversity, as in the open shrubland site ES-LJu (Fig. 3f) (Serrano-Ortiz et al., 2009). Such generalization can also cause systematic inaccuracies in parameterization, where model parameters are better suited for European semiarid vegetation (e.g., CH-Fru, Imer et al., 2013; IT-MBo, Fig. 3b, Marcolla et al., 2011; CZ-wet, Fig. 3c;, Dušek et al., 2012) than similar sites in other regions (e.g., AU-Dry, Hutley et al., 2011; RU-Sam, Fig. 3g, Boike et al., 2013; US-SRG, Scott et al., 2015).

Simulated results for forest PFTs compare very well with observations, where N values are often greater than 0.5. TEMIR performs particularly well for evergreen needleleaf forests as seen in sites DE-Tha (Fig. 3h), FI-Hyy (Fig. 3i), and NL-Loo (Fig. 3j), which are mostly populated by mature Scots pine forests of over 70 years old. Sites CA-TP1 (Fig. 3k), CA-TP3 (Peichl et al., 2010; Arain et al., 2022), and DE-Lkb (Lindauer et al., 2014) are overestimated by the model as these forests are dominated by eastern white pine and Norway spruce that are less than 20 years old, so optimal productivity might not have been achieved. In comparison, the neighboring site CA-TP4 (Fig. 3l; Peichl et al., 2010) with over 70-year-old eastern white pine is better replicated by the model. The model has better performance with respect to site observations when using FLUXNET local meteorology (e.g., CH-Dav in Fig. 3m; Zielis et al., 2014), though the differences are insignificant for most sites. For deciduous broadleaf forest sites, although represented well overall, there is a systematic underestimation (e.g., FR-Fon in Fig. 3n), most likely due to inaccurate parameterization overcompensating for the uncertainties in satellite-derived LAI for broadleaf trees. The multi-year drought in the USA during the 2010s, which hinders plant productivity (Wolf et al., 2016; Xu et al., 2020), appears to improve model agreement by reducing the discrepancy (i.e., US-Oho and US-UMB) and even giving a positive model bias (i.e., US-MMS in Fig. 3o; Yi et al., 2017).

The correlation for croplands is high, but there is a spread in B giving varying N values. The range of model performance among cropland sites shows the limitation of the simplistic crop representation used in this version of TEMIR, whereby site-level settings, such as planting seasons and agricultural management (e.g., fertilizer usage, irrigation, and possible rotations between crop types and cultivars), are not considered. The generic crop representation fails to capture the maximum photosynthetic capacity of the planted crops. For example, site US-Ne1 has irrigated maize that has much higher GPP compared with the simulated generic crop as shown in Fig. 3p. Site DE-Kli (Fig. 3q) has a 5-year crop rotation with occasional fertilizer application (Prescher et al., 2010) and has higher productivity than that simulated by the model.

Figure 3Diurnal averaged gross primary productivity of selected sites representative of their respective vegetation types from simulations described in Table 1, with relevant site information annotated: (a) AU-How, (b) IT-MBo, (c) CZ-wet, (d) CH-Cha, (e) AU-Stp, (f) ES-LJu, (g) RU-Sam, (h) DE-Tha, (i) FI-Hyy, (j) NL-Loo, (k) CA-TP1, (l) CA-TP4, (m) CH-Dav, (n) FR-Fon, (o) US-MMS, (p) US-Ne1, and (q) DE-Kli. More details of these sites are given in Table S3.


5.1.2 Validation on seasonal cycle

Figure 4 shows the statistical metrics (Table 4) for monthly GPP averages of 2009–2013 to examine the model performance in seasonal GPP cycle. All plots of monthly cycles are shown in the Supplement, with plots of selected sites included in the following discussion. The model generally performs worse in capturing the seasonal cycle than the diurnal cycle. Between the different settings of meteorology used for simulations (Fig. 4), the differences in statistics are small. MERRA-2 meteorology shows good utility for most sites, with in-canopy meteorology inferred using Monin–Obukhov similarity theory improving correlation for some sites. FLUXNET local meteorology gives the smallest range of biases with performance similar to MERRA-2 meteorology simulations.

Figure 4Statistical metrics (see Table 4) of the monthly gross primary productivity from simulations (a) TEMIR_FLUX, (b) TEMIR_MO_off, and (c) TEMIR_MO_on as described in Table 1.


Comparing Figs. 2 and 5, model–observation comparison of monthly averages gives lower values of R2 for all sites in general. On the other hand, biases are distributed more evenly across the range with smaller extreme values compared with the biases from diurnal simulations. In terms of N, the model is less adept in reproducing seasonal variations (due to the reductions in correlation) regardless of the driving meteorology chosen. Figure 5c shows that deciduous broadleaf sites (e.g., US-MMS in Fig. 6a) give R2>0.85 with the maximum absolute B=+0.54 and minimum N=0.42 with a mean of 0.9. Monthly performance of forest sites shows a smaller range in B and lower absolute mean values of B in comparison with the other PFTs (Fig. 5). The prediction for forest sites is satisfactory with mean values of N larger than 0.65 (e.g., RU-Fyo in Fig. 6b; CA-TP4 in Fig. 6c), and monthly GPP inaccuracies for forest sites can be explained with similar reasoning as discussed in Sect. 5.1.1.

The correlation for grasslands is above 75 % for most sites (e.g., DE-Akm in Fig. 6d; IT-MBo in Fig. 6e), while sites AU-How (Fig. 6f) and AU-Stp have an R2 below 0.4. These sites are known to have fires occurring in the dry winter and spring from May to October, which corresponds to the low productivities (Beringer et al., 2007, 2011; Hutley et al., 2011; Haverd et al., 2013). Moreover, such disturbances on LAI with vegetation regrowth are complex and often overlooked by the model as shown in the simulation of the site AU-How with minimal productivity.

The simulation of monthly GPP of croplands most clearly shows the limitations of the generic model approach, as no site-specific crop phenology is available in this version of the model. The simulated seasonal cycle shows a typical annual peak usually in summer as dictated by meteorology that in general can yield good correlation (e.g., DE-Geb in Fig. 6g); yet, as many sites are intensively managed, the observed GPPs do not follow such a simplistic cycle, giving low correlations (e.g., IT-BCi in Fig. 6h). The parameters for a generic crop usually fail to represent the actual crop planted at the sites, and therefore large biases exist in the simulated GPP (e.g., BE-Lon in Fig. 6i; FR-Gri in Fig. 6j). GPP is more commonly underestimated at the US sites (e.g., US-Ne1 in Fig. 6k) where maize is usually planted and is more productive in comparison with other crops.

Figure 5Statistical metrics (see Table 4) of the monthly averaged gross primary productivity from the TEMIR_FLUX simulation (see Table 1) for each plant functional type listed in Table S1: (a) evergreen needleleaf forest (ENF), (b) evergreen broadleaf forest (EBF), (c) deciduous broadleaf forest (DBF), (d) open shrubland (OSH), (e) grassland (GRA), and (f) cropland (CRO).


Figure 6Monthly averaged gross primary productivity of selected sites representative of their respective vegetation types from simulations described in Table 1, with relevant site information annotated: (a) US-MMS, (b) RU-Fyo, (c) CA-TP4, (d) DE-Akm, (e) IT-MBo, (f) AU-How, (g) DE-Geb, (h) IT-BCi, (i) BE-Lon, (j) FR-Gri, and (k) US-Ne1.


5.2 Global validation

Figure 7 shows that the simulated annual averaged GPP for 2010–2015 is 134.7 Pg C yr−1 from the simulation using MERRA-2 meteorology (TEMIR_MO_off) (Fig. 7b), and the simulation with in-canopy meteorology inferred using MOST (TEMIR_MO_on) gives GPP over the same period as 144.7 Pg C yr−1 (Fig. 7c). Compared with the satellite-derived dataset (GOSIF GPP; Sect. 3.2) where the annual GPP of the same period is 128.4 Pg C yr−1 (Fig. 7a), TEMIR overestimates global GPP by  5 %–10 % depending on the input meteorology (Table 5). TEMIR performance is well within, and leans toward, the middle of the observation-constrained range in the literature of 119–175 Pg C yr−1. TEMIR closely agrees with models of similar design objectives, e.g., the Yale Interactive terrestrial Biosphere (YIBs) with GPP at 125±3 Pg C yr−1 (Yue and Unger, 2015) and JULES land surface model estimating GPP at 141 Pg C yr−1 (Slevin et al., 2017). TEMIR can largely reproduce the spatial distribution of GPP with respect to GOSIF GPP (Fig. 7), with grid cells with mixed savanna and forests showing larger discrepancies.

Table 5Gross primary productivity (GPP) product simulated and relevant simulation details of global TEMIR simulations described in Tables 2 and 3. “N/A” for O3 year indicates that no O3 damage scheme is used.

Download Print Version | Download XLSX

Figure 8 shows model–observation statistics (Table 4) for the model outputs of Fig. 7 in 12 regions. Global correlation of 6-year averaged GPP is around 83 % (Fig. 8b). In general, correlations are lower for regions closer to the Equator, where correlations for tropical regions are below 70 % but otherwise above 70 %, with correlation for Siberia close to 90 %. Correlation for the tropical Americas is  75 %, which is higher than other equatorial regions. Temperate North America shows a correlation of  60 %, which is lower than other regions at midlatitudes. We also see that GPP is underestimated in the grid cells with high crop density (Fig. 7b), which, as discussed in Sect. 5.1, is likely due to the generic crop representation of the TEMIR version giving poor model performance for this region. Simulated GPP driven by meteorology inferred with MOST gives small increases or decreases in regional correlation (e.g., correlation for North Africa and the Middle East drops from 48 % to 45 %, and correlation for temperate South America increases from 75 % to 78 %).

Figure 7Panel (a) shows the average global gross primary productivity for 2010–2015 from the GOSIF GPP product and differences in the simulated GPP from simulations. Panel (b) shows TEMIR_MO_off and (c) shows TEMIR_MO_on (see Table 1) for 2010–2015 compared with the GOSIF product.

Absolute biases are mostly within 25 % except for North Africa, the Middle East, and sub-Saharan Africa (except central Africa) when the driving meteorology is inferred from MOST. Simulated GPP driven by in-canopy meteorology inferred with MOST gives more positive biases for all regions, generally around +10 %. GPP for 6 of the 12 regions is overestimated by TEMIR and otherwise underestimated (Fig. 8b). Thus, in-canopy meteorology inferred with MOST results in regional bias changes from underestimation to overestimation for East and South Asia, Europe, temperate North America, and the tropical Americas. The bias for Europe is 6.2 % and +1.19 % when driving meteorology is inferred with MOST. The bias of +1.19 % is the smallest absolute bias of any region, which shows the possibility of in-canopy meteorology inferred with MOST improving GPP predictions for some but not all regions.

Figure 8Panel (a) shows regional division relevant for this study largely following Chen et al. (2017), and (b) shows the corresponding regional statistical metrics (see Table 4) of the averaged gross primary productivity for 2010–2015 from TEMIR simulations (see Table 2).

5.3 Effects of O3 and CO2 on global primary productivity

Figure 9 shows the simulated results where the Sitch O3 schemes (Sect. 2.2.5) of low sensitivity (Sl) and high sensitivity (Sh) are implemented for 2010–2015 using MERRA-2 meteorology. Figure 9a shows the mean daily 8 h averaged O3 concentration (MDA8), a common surface O3 metric, derived from the simulated hourly O3 concentration at the lowest model level affecting global vegetation under the Sitch O3 damage scheme. The global GPP values are 133.6 and 131.8 Pg C yr−1 for the Sitch O3 scheme at low and high sensitivity, respectively (Table 5), both of which are smaller than the 134.7 Pg C yr−1 from the simulation without O3 damage (Fig. 7b). These global GPP reductions are seemingly small (< 1 % to  2 %) and conceal larger regional changes. Figure 9c shows that the Sitch O3 damage scheme at high sensitivity leads to an up to 15 % reduction in GPP, whereas low sensitivity shows modest reductions of about half of those magnitudes. Particularly large O3-induced damage occurs in highly populated regions (e.g., the eastern USA, Europe, central Africa, northern India, and East Asia) associated with high anthropogenic emissions (NOx in particular). Many of these regions also contain arable lands, and thus O3 exposure can also affect food security (Feng et al., 2008; Avnery et al., 2011; Emberson et al., 2018; Ainsworth et al., 2020; Tai et al., 2021; Leung et al., 2022; Roberts et al., 2022).

Figure 8b shows the statistics of Table 4 per region (Fig. 8a) for simulations with O3 damage. The presence of O3 does not affect the model–observation correlations significantly for any region. When compared with the correlations of TEMIR_MO_off simulation results, correlations from O3-damaged GPP show small differences. O3 damage reduces the model overestimation with respect to GOSIF GPP. In particular, for eastern China and central Africa, implementing O3 damage reduces the positive model biases as seen in Fig. 7. Underestimation is worsened for the regions of temperate North America as well as East and South Asia where there is strong O3 damage (Fig. 9).

Figure 9Mean daily 8 h averaged (MDA8) O3 concentration of the lowest model layer averaged over 2010–2015 and percentage differences in average global GPP for 2010–2015 of the simulated results with the Sitch O3 damage scheme at (b) low sensitivity and (c) high sensitivity from the simulation TEMIR_MO_off (see Table 2).

Figure 10 shows the comparisons between simulations (Table 3) displaying the interplay of CO2 fertilization effects and O3 damage to GPP. CO2 fertilization (from 370 to 390 ppmv), shown in Fig. 10b, promotes regional productivity by up to 7 %. Global GPP enhancement is  2 % (Table 5), and thus simulations estimate that rising atmospheric CO2 concentration results in a global GPP increase of 0.126 % ppmv−1. As seen in Fig. 10c, O3-induced regional reductions are up to 15 % under the Sitch O3 damage scheme at high sensitivity, whereby results are similar to those in Fig. 9c. Figure 10d shows the differences in percentage of O3 damage to GPP of the simulation with O3 damage at a CO2 concentration of 390 ppmv from that at 370 ppmv (i.e., that of Fig. 9c). The positive values in Fig. 10d indicate that the O3-induced GPP reduction is smaller at a higher CO2 concentration, reflecting the additional benefits of CO2 fertilization from the reduced stomatal conductance, which improves water use efficiency and also decreases stomatal O3 uptake thus lessening O3-induced impacts.

Figure 10Plots showing results from simulations of Table 3: (a) gross primary productivity modeled for 2000 at a CO2 concentration of 370 ppmv (C0_O0), (b) percentage changes in GPP showing CO2 fertilization effects of year-2010 CO2 concentration at 390 ppm (100 % × (C1_O0 C0_O0)/C0_O0), (c) percentage changes in GPP due to O3 damage at high sensitivity of the Sitch O3 damage scheme for year-2000 modeled O3 concentration and CO2 concentration of 370 ppmv (100 % × (C0_O1 C0_O0)/C0_O0), and (d) differences in percentage of O3 damage at a CO2 concentration of 390 ppmv from that at 370 ppmv (100 % × (C1_O1 C1_O0)/C1_O0 100 % × (C0_O1 C0_O0)/C0_O0), whereby positive values indicate a reduction in percentage of O3 damage.

6 Discussion and conclusions

In this paper we provide a detailed model description of the newly developed Terrestrial Ecosystem Model in R (TEMIR) version 1.0, which simulates ecophysiological processes and functions (most importantly, photosynthesis and global primary productivity) of terrestrial ecosystems as represented by different PFTs, driven by prescribed meteorological conditions and atmospheric chemical composition. We specifically include the multiple parameterization schemes for stomatal O3 uptake and O3 damage to plants, as well as showcasing the utility of TEMIR in evaluating the responses of GPP to O3 damage, CO2 fertilization, and their interactions. The productivity simulated at site and global levels reproduces the observed diurnal and seasonal cycles well for evergreen needleleaf and deciduous broadleaf forests (especially those that are mature), with an annual average GPP of 134.7 Pg C yr−1 for 2010–2015 and a global reduction of up to 2 % when O3 damage is considered. This is validated against the productivity from the 49 FLUXNET sites and GOSIF GPP.

TEMIR-simulated global GPP lies well within the accepted range, but the associated large uncertainty is well acknowledged in the field (Bonan et al., 2011; Baldocchi et al., 2016; Zhang et al., 2017; Li and Xiao, 2019; Bi et al., 2022; Wild et al., 2022; Zhang and Ye, 2022), thus limiting the validity of global GPP model–observation comparison in this study (Sect. 5.2). Site-level validation may lend more credence by isolating certain PFTs for comparison, albeit being more limited in scope and scale unlike global comparisons. Our investigation suggests that possible PFT systematic biases exist generally for diurnal productivity, which reflect the limitations of having a set of prescribed parameters for generalized classes of plant functions (Harrison et al., 2021; Seiler et al., 2022; Liu et al., 2023; Wu et al., 2023a). For instance, there is a systematic underestimation for deciduous broadleaf forests, though it can be explained by the uncertainties in LAI datasets (Liu et al., 2018; Yang et al., 2023), and some regions show distinctive physiology and phenology of grasses and shrubs. Particularly for semiarid regions where the range of productivity is large, the model shows variable accuracy. In general, variability in prescribed LAI can be an important source of uncertainty in the model results. Single-site sensitivity simulations show that GPP generally linearly increases with LAI at low LAI, but as LAI becomes larger, GPP increases less than proportionately due to the canopy shading effect. Such nonlinearity of GPP responses to LAI changes is less important for small perturbations of LAI (e.g., < 20 %).

Simulating crops in ecosystem modeling remains particularly challenging (Deryng et al., 2016; Chopin et al., 2019; Muller and Martre, 2019; Boas et al., 2021), as it combines the nuances in phenology, physiology, coverage, and active human management with high spatiotemporal variations (Monfreda et al., 2008; Emberson et al., 2018; Ahmed et al., 2022; Gleason et al., 2022; Corcoran et al., 2023), which already exist for natural vegetation to a lesser degree. One particularly crucial aspect for improvement is to get crop LAI correct, which is typically more challenging to measure than trees with large canopies and often varies to greater extents with leaf orientations for different crops. More long-term ground-based and/or remote sensing measurements of crop LAI for different crop types across the world are particularly recommended, not only as input data but also for model validation in future development. Especially for site-level simulations, locally relevant crop physiological and structural parameters should also be measured and used. Ongoing development has already been attempting to enhance crop representation in a version of TEMIR with active crop biogeochemistry (Tai et al., 2021) to improve and reconcile model inaccuracies.

Incorporating site-level meteorology in simulations can improve performance for a few selected sites but otherwise is comparable to results from simulations with gridded assimilated meteorology as input. This highlights the fact that generalization and the coarse resolution of the MERRA-2 dataset used (due to computational limitations and necessary consistency with other input datasets) can drastically overlook regional and small-scale nuances. Furthermore, CO2 concentration was kept constant and spatially uniform in all simulations, which enables direct comparison with other modeling studies but ignores possible spatiotemporal variability in CO2 concentration (Cheng et al., 2022). Though such effects are usually minor on simulated GPP magnitudes (Lee et al., 2018; Tian et al., 2021), uncertainties should be minimized in any case; thus, it is recommended that users use the measured CO2 concentration, if available, as input, especially for site-level simulations. It is also recommended that users recalibrate relevant model parameters with site observations and available datasets (e.g., those of higher resolutions), such as LAI, Vcmax, PFT fractional coverage, and others, to yield the most accurate results. The non-dynamic representation of vegetation cover and parameterization is a shortcoming of TEMIR, and thus simulations overlook intricate and transient impacts of LULCC on land–atmosphere exchange (Ganzeveld et al., 2010; Pongratz et al., 2010; Prescher et al., 2010; Chen et al., 2018; Bastos et al., 2020; Hou et al., 2022). With the capacity of the current version of TEMIR, our simulations address these aspects by changing the input data of LAI and PFT fractions, derived from LULCC, for yearly and higher frequencies. LULCC can drive large regional changes, though recent LULCC mostly reduces GPP (due to urbanization, agricultural expansion, and deforestation), counteracted partly by CO2 fertilization effects (Wu et al., 2023b), and thus the validity of our results is likely unchanged. The assumption of sufficient nitrogen availability is a limitation (Sandor et al., 2018) as most non-tropical biomes experience varying nitrogen limitation (Davies-Barnard et al., 2022; Kou-Giesbrecht and Arora, 2023), thereby affecting photosynthetic capacities (Mason et al., 2022; Wang et al., 2022) and resource allocations in plants and soil (Zhang et al., 2020; Feng et al., 2023; Lu et al., 2023). Some models have N cycling (Yang et al., 2009; Gerber et al., 2010; Wiltshire et al., 2021; Hidy et al., 2022), but effects remains minor (Jain et al., 2009; O'Sullivan et al., 2019; Lin et al., 2023) in the recent decade and more relevant for assessing future global changes (Tharammal et al., 2018; Franz and Zaehle, 2021). Overall, TEMIR has great skill in capturing annual and seasonal GPP at the global scale as well as for some productive regions and certain PFTs, whereby the correlation is high in the range of 80 %–90 %, showcasing the utility of TEMIR at different scales. Caution should be taken with good knowledge of model preferences and the underlying theoretical assumptions for any given research question, especially when concerning multi-factor land–atmosphere interactions and vegetation responses to various environmental stresses (Kimmins et al., 2008; Zhao et al., 2022; Blanco and Lo, 2023; Rahman et al., 2023). Further development and validation of the model with detailed observations are crucial to provide more accurate vegetation parameterization for specific applications, e.g., to investigate vegetation responses to droughts and heatwave composition (e.g., Yan et al., 2022), especially at the regional and site levels.

The initial motivation and one of the most relevant applications of TEMIR is to address the impacts of O3 pollution and exposure on terrestrial ecosystem productivity, whereby an active Sitch O3 damage scheme improves model performance with respect to GPP. Concerning O3 damage to GPP, there is good agreement with previous studies in terms of both magnitudes and spatial variations (e.g., Lombardozzi et al., 2015; Sitch et al., 2007). For instance, the OCN model (Franz et al., 2017; Franz and Zaehle, 2021) simulated O3 as reducing GPP in Europe by  8 % and the JULES land surface model (Slevin et al., 2017) in the range of 10 %–20 % (Oliver et al., 2018). The Yale Interactive terrestrial Biosphere (YIBs) model (Yue and Unger, 2015) simulated O3 as reducing global GPP by 2 %–5 % with East Asia experiencing damage of 4 %–10 %. Yue and Unger (2014) also showed GPP reductions of 4 %–8 % in the eastern USA with high episodes giving a higher range of 11 %–17 %. YIBs has the capability of synchronous coupling (e.g., GEOS-Chem–YIBs; Lei et al., 2020), which reported similar ranges in GPP reductions: globally by 1.5 %–3.6 % and extremes of 11 %–14 % in the eastern USA and eastern China. This lends credence to the comparable performance of TEMIR v1.0, which has a more simplistic terrestrial ecosystem with prescribed ecosystem structure (noting that active biogeochemistry is in development).

O3 influences in the current version of TEMIR are limited to vegetation physiological and productivity responses. Intra- and interspecies differential sensitivity to O3 can cause competition (Agathokleous et al., 2020), affecting some species more than others in terms of biomass, flowering, and seed development, thus impacting community composition, PFT fractional coverage, and biodiversity (Calvete-Sogo et al., 2016; Fuhrer et al., 2016; Emberson, 2020). This can also be seen among functional groups. For example, perennial species retain more aboveground biomass than annual species, and angiosperms are more prone to O3 damage than gymnosperms, thus giving possible long-term biodiversity effects (Agathokleous et al., 2020). Such effects are further complicated by soil conditions (e.g., water and nitrogen content) and also spatial heterogeneity, whereby regional strategies might differ within functional groups, although this requires more studies to obtain observation-based parameterization.

Moreover, synchronous model coupling between a CTM or climate model and a fully prognostic biosphere model with active biogeochemistry is particularly suitable for examining O3–vegetation feedbacks (Danabasoglu et al., 2020; Franklin et al., 2020; Lam et al., 2023), especially for timescales long enough (e.g., multi-decadal) for ecosystem structure to co-evolve with the atmosphere. For instance, Sadiq et al. (2017) and Gong et al. (2021) showed that dynamic O3–vegetation interactions can lead to a long-term ecosystem decline and a positive feedback on O3 concentration in China and worldwide, respectively, worsening air quality. Zhu et al. (2022) and Jin et al. (2023) found similar positive O3–vegetation feedbacks in China with the coupled framework using WRF-Chem and Noah-MP. Yue et al. (2017) also investigated O3–aerosol–vegetation interactions in China. TEMIR can only be asynchronously coupled with GEOS-Chem and is not the best tool for investigating two-way O3–vegetation interactions, especially when such interactions relevantly happen within a model time step, but it is particularly suitable for estimating first-order effects of O3 pollution on vegetation in a computationally efficient manner. Zhou et al. (2018) indeed found that second-order effects of O3 pollution (i.e., additional effects of modified O3 concentrations after feedbacks are accounted for) on vegetation are negligible. Moreover, asynchronous coupling between TEMIR and GEOS-Chem, for example, and conducting factorial experiments with them, can help disentangle complex pathways and feedbacks that are often convoluted in fully coupled models.

We recognize that the O3 damage scheme in TEMIR does not account for sluggishness in stomatal responses (e.g., Clifton et al., 2020; Huntingford et al., 2018), which may modify further O3 uptake, although such effect is expected to be small at the resolution relevant for this study. O3 sensitivities also have crop-related inaccuracies due to the generic crop representation in this version of TEMIR. Such is a common practice in global-scale biosphere models, and Leung et al. (2020) suggested that if a study focuses on crop yields, species-specific calibration is required to reduce uncertainty and likely inaccuracies for the crops concerned. TEMIR v1.0 on a global scale is not suitable for any crop-focused investigations, but one may use the version of TEMIR implemented with additional crop functionalities, such as the calculation of phytotoxic O3 dose, taking advantage of the stomatal calculation in TEMIR and the subsequent estimation of O3–crop impacts (Tai et al., 2021). The utility of TEMIR in examining vegetation-mediated dry depositional sinks of O3 has also been demonstrated (Sun et al., 2022).

Mechanistic representations allow modeling for various meteorological conditions and are extremely useful to evaluate ecophysiological responses to a changing climate and intermittent climate extremes (e.g., Bonan, 2008, 2016; Cai and Prentice, 2020; Gang et al., 2022). Ciais et al. (2005) estimated a 30 % GPP reduction in Europe following the heatwave in 2003 and vegetation there became a net carbon source, attributable to the rainfall deficit and extreme summer heat. This was also found by Bamberger et al. (2017), whereby heat and drought impacts alter photosynthesis and vegetation state. More extreme events are projected for a future climate, which various models (e.g., O-CN and YIBs) have shown to decrease productivity (e.g., Franz and Zaehle, 2021; Yan et al., 2022). He et al. (2022), using models, showed that climate variability is the main factor controlling interannual GPP variability in grasslands in China. Such effect is the most prominent in summer, which is responsible for more than 40 % of decadal GPP variability in Chinese grasslands and the largest in comparison with effects from CO2 fertilization and nitrogen deposition. Similar to the case for O3–vegetation coupling discussed above, fully coupled climate–biosphere models can be particularly useful for examining two-way interactions and feedbacks and also long-term (multi-decadal to multi-centurial) co-evolution of climate and the biosphere. However, the embedded complex interactions may obscure the relative importance of different factors, making it a lot more difficult to attribute changes to specific factors. Offline modes, such as TEMIR, are therefore particularly useful for investigating and attributing biospheric variability and changes to prescribed changes in climatic variables.

In addition, we have demonstrated the utility of TEMIR in examining the direct and interactive effects of multiple atmospheric chemical species on global vegetation (i.e., CO2 and O3 concentrations). CO2 fertilization in TEMIR results in strong GPP enhancement as seen in many studies (e.g., Schimel et al., 2015; Cai and Prentice, 2020; Chen et al., 2022; Yang et al., 2022). Our simulations estimate that CO2 fertilization increases global GPP by 0.126 % ppmv−1, which is comparable to the value of 0.138 % ± 0.007 % ppmv−1 reported by Ueyama et al. (2020). It is noteworthy that some studies (e.g., Lee et al., 2018) suggest that overlooking spatiotemporal variability in atmospheric CO2 can lead to inaccuracies in seasonal and regional GPP estimation but only a minor influence on global GPP. Additional crop functionalities of TEMIR (Tai et al., 2021) can also address the CO2 fertilization effects on crops, although studies have also found that the resultant productivity increase gives larger yield quantity but does not necessarily translate to increased yield quality (Myers et al., 2014; Ebi et al., 2021; Xia et al., 2021). The competing effects of CO2 fertilization and O3 damage to vegetation have been well documented in field experiments, although magnitudes vary and are species-dependent (e.g., Oikawa and Ainsworth, 2016; Proietti et al., 2016; Karlsson et al., 2017; Moura et al., 2018; Zhang et al., 2018; Ainsworth et al., 2020; Xia et al., 2021). TEMIR shows that CO2 fertilization can reduce the percentage of O3 damage to vegetation ( 1 % globally), which is generally comparable to 1 %–2 % found by Oliver et al. (2018), whereas Sitch et al. (2007) simulated a higher range of 6 %–9 %. We note that while comparisons among models are useful, we must be mindful of the differences in model designs and setups (as mentioned in Sect. 1). Miner et al. (2017) cautioned that stomatal responses to CO2 can be highly species-dependent and variable under different soil conditions, adding more uncertainty to the parameterization of CO2–O3–vegetation interactions. The changing nitrogen deposition due to anthropogenic activities may likewise influence the interactions between vegetation, CO2, and O3 (Zhao et al., 2017; Liu et al., 2021), whereby nitrogen can limit CO2-promoted growth (Wang et al., 2020) or modify vegetation responses (e.g., for gs; Hu et al., 2021) with further implications on soil and nutrient cycling (Terrer et al., 2021). As atmospheric composition rapidly changes in the next century, these interactive mechanisms should be considered for modelers to more representatively and accurately model the future earth system (e.g., Bytnerowicz et al., 2007; Pu et al., 2017; Sicard et al., 2017; Franz and Zaehle, 2021; Leung et al., 2022).

All in all, the high adaptability of TEMIR, written in an freely open-source, widely used, and easy-to-learn programming language (R Core Team, 2022), is expected to facilitate fruitful contribution to research, at various spatiotemporal scales, on biosphere–atmosphere interaction. It also provides a readily available tool for policymakers, practitioners, and other stakeholders to assess the ecosystem services provided by vegetation in different regions or cities, as well as their sensitivities to future atmospheric changes, possibly enhancing the translational value of ecological and geoscientific research.

Code availability

The Terrestrial Ecosystem Model in R (version 1.0) source code is licensed and publicly available in the following repository: (Tai and Yung, 2023a); as well as on GitHub: (last access: 21 March 2024).

Data availability

The input dataset required to run Terrestrial Ecosystem Model in R (version 1.0) is licensed and publicly available in the following repository: (Tai and Yung, 2023b).


The supplement related to this article is available online at:

Author contributions

APKT conceived the study and developed TEMIR. DHYY developed additional model codes as well as performed the simulations and analysis. APKT and DHYY prepared the paper. TL developed multilayer soil representation and contributed to relevant sections of 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.


This work was supported by the National Natural Science Foundation of China (NSFC)/Research Grants Council (RGC) Joint Research Scheme (reference nos. N_CUHK440/20 and 42061160479) as well as the Collaborative Research Fund (reference no. C5062-21GF) of RGC, Hong Kong SAR, China.

Financial support

This research has been supported by the Research Grants Council, University Grants Committee (grant nos. N_CUHK440/20 and C5062-21GF).

Review statement

This paper was edited by Christoph Müller and reviewed by two anonymous referees.


Agathokleous, E., Feng, Z., Oksanen, E., Sicard, P., Wang, Q., Saitanis, C. J., Araminiene, V., Blande, J. D., Hayes, F., Calatayud, V., Domingos, M., Veresoglou, S. D., Peñuelas, J., Wardle, D. A., Marco, A. D., Li, Z., Harmens, H., Yuan, X., Vitale, M., and Paoletti, E.: Ozone affects plant, insect, and soil microbial communities: A threat to terrestrial ecosystems and biodiversity, Sci. Adv., 6, eabc1176,, 2020. 

Ahmed, Z., Gui, D., Qi, Z., Liu, Y., Liu, Y., and Azmat, M.: Agricultural system modeling: current achievements, innovations, and future roadmap, Arab. J. Geosci., 15, 363,, 2022. 

Ainsworth, E. A., Yendrek, C. R., Sitch, S., Collins, W. J., and Emberson, L. D.: The Effects of Tropospheric Ozone on Net Primary Productivity and Implications for Climate Change, Annu. Rev. Plant Biol., 63, 637–661,, 2012. 

Ainsworth, E. A., Lemonnier, P., and Wedow, J. M.: The influence of rising tropospheric carbon dioxide and ozone on plant productivity, Plant Biol. J., 22, 5–11,, 2020. 

Anav, A., Menut, L., Khvorostyanov, D., and Viovy, N.: Impact of tropospheric ozone on the Euro-Mediterranean vegetation, Glob. Change Biol., 17, 2342–2359,, 2011. 

Arain, M. A., Xu, B., Brodeur, J. J., Khomik, M., Peichl, M., Beamesderfer, E., Restrepo-Couple, N., and Thorne, R.: Heat and drought impact on carbon exchange in an age-sequence of temperate pine forests, Ecol. Process., 11, 7,, 2022. 

Arneth, A., Harrison, S. P., Zaehle, S., Tsigaridis, K., Menon, S., Bartlein, P. J., Feichter, J., Korhola, A., Kulmala, M., O'Donnell, D., Schurgers, G., Sorvari, S., and Vesala, T.: Terrestrial biogeochemical feedbacks in the climate system, Nat. Geosci., 3, 525–532,, 2010. 

Arora, V. K., Boer, G. J., Friedlingstein, P., Eby, M., Jones, C. D., Christian, J. R., Bonan, G., Bopp, L., Brovkin, V., Cadule, P., Hajima, T., Ilyina, T., Lindsay, K., Tjiputra, J. F., and Wu, T.: Carbon–Concentration and Carbon–Climate Feedbacks in CMIP5 Earth System Models, J. Climate, 26, 5289–5314,, 2013. 

Avnery, S., Mauzerall, D. L., Liu, J., and Horowitz, L. W.: Global crop yield reductions due to surface ozone exposure: 1. Year 2000 crop production losses and economic damage, Atmos. Environ., 45, 2284–2296,, 2011. 

Baldocchi, D., Ryu, Y., and Keenan, T.: Terrestrial Carbon Cycle Variability, F1000Research, 5, 2371,, 2016. 

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., and Biggins, J., Dordrecht, 221–224, 1987. 

Bamberger, I., Ruehr, N. K., Schmitt, M., Gast, A., Wohlfahrt, G., and Arneth, A.: Isoprene emission and photosynthesis during heatwaves and drought in black locust, Biogeosciences, 14, 3649–3667,, 2017. 

Bastos, A., O'Sullivan, M., Ciais, P., Makowski, D., Sitch, S., Friedlingstein, P., Chevallier, F., Rödenbeck, C., Pongratz, J., Luijkx, I. T., Patra, P. K., Peylin, P., Canadell, J. G., Lauerwald, R., Li, W., Smith, N. E., Peters, W., Goll, D. S., Jain, A. K., Kato, E., Lienert, S., Lombardozzi, D. L., Haverd, V., Nabel, J. E. M. S., Poulter, B., Tian, H., Walker, A. P., and Zaehle, S.: Sources of Uncertainty in Regional and Global Terrestrial CO2 Exchange Estimates, Global Biogeochem. Cy., 34, e2019GB006393,, 2020. 

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. 

Beringer, J., Hutley, L. B., Tapper, N. J., and Cernusak, L. A.: Savanna fires and their impact on net ecosystem productivity in North Australia, Glob. Change Biol., 13, 990–1004,, 2007. 

Beringer, J., Hutley, L. B., Hacker, J. M., Neininger, B., and Paw U, K. T.: Patterns and processes of carbon, water and energy cycles across northern Australian landscapes: From point to region, Agr. Forest Meteorol., 151, 1409–1416,, 2011. 

Bey, I., Jacob, D. J., Yantosca, R. M., Logan, J. A., Field, B. D., Fiore, A. M., Li, Q., Liu, H. Y., Mickley, L. J., and Schultz, M. G.: Global modeling of tropospheric chemistry with assimilated meteorology: Model description and evaluation, J. Geophys. Res.-Atmos., 106, 23073–23095,, 2001. 

Bhattarai, H., Tai, A. P. K., Martin, M. V., and Yung, D. H. Y.: Impacts of changes in climate, land use, and emissions on global ozone air quality by mid-21st century following selected Shared Socioeconomic Pathways, Sci. Total Environ., 906, 167759,, 2023. 

Bi, W., He, W., Zhou, Y., Ju, W., Liu, Y., Liu, Y., Zhang, X., Wei, X., and Cheng, N.: A global 0.05° dataset for gross primary production of sunlit and shaded vegetation canopies from 1992 to 2020, Sci. Data, 9, 213,, 2022. 

Blanco, J. A. and Lo, Y.-H.: Latest Trends in Modelling Forest Ecosystems: New Approaches or Just New Methods?, Current Forestry Reports, 9, 219–229,, 2023. 

Boas, T., Bogena, H., Grünwald, T., Heinesch, B., Ryu, D., Schmidt, M., Vereecken, H., Western, A., and Hendricks Franssen, H.-J.: Improving the representation of cropland sites in the Community Land Model (CLM) version 5.0, Geosci. Model Dev., 14, 573–601,, 2021. 

Boike, J., Kattenstroth, B., Abramova, K., Bornemann, N., Chetverova, A., Fedorova, I., Fröb, K., Grigoriev, M., Grüber, M., Kutzbach, L., Langer, M., Minke, M., Muster, S., Piel, K., Pfeiffer, E.-M., Stoof, G., Westermann, S., Wischnewski, K., Wille, C., and Hubberten, H.-W.: Baseline characteristics of climate, permafrost and land cover from a new permafrost observatory in the Lena River Delta, Siberia (1998–2011), Biogeosciences, 10, 2105–2128,, 2013. 

Bonan, G. B.: Forests and Climate Change: Forcings, Feedbacks, and the Climate Benefits of Forests, Science, 320, 1444–1449,, 2008. 

Bonan, G. B.: Forests, Climate, and Public Policy: A 500-Year Interdisciplinary Odyssey, Annu. Rev. Ecol. Evol. S., 47, 97–121,, 2016. 

Bonan, G. B., Levis, S., Kergoat, L., and Oleson, K. W.: Landscapes as patches of plant functional types: An integrating concept for climate and ecosystem models, Global Biogeochem. Cy., 16, 5-1–5-23,, 2002. 

Bonan, G. B., Lawrence, P. J., Oleson, K. W., Levis, S., Jung, M., Reichstein, M., Lawrence, D. M., and Swenson, S. C.: Improving canopy processes in the Community Land Model version 4 (CLM4) using global flux fields empirically inferred from FLUXNET data, J. Geophys. Res., 116, G02014,, 2011. 

Bonan, G. B., Oleson, K. W., Fisher, R. A., Lasslop, G., and Reichstein, M.: Reconciling leaf physiological traits and canopy flux data: Use of the TRY and FLUXNET databases in the Community Land Model version 4, J. Geophys. Res.-Biogeo., 117, G02026,, 2012. 

Bytnerowicz, A., Omasa, K., and Paoletti, E.: Integrated effects of air pollution and climate change on forests: A northern hemisphere perspective, Environ. Pollut., 147, 438–445,, 2007. 

Cai, W. and Prentice, I. C.: Recent trends in gross primary production and their drivers: analysis and modelling at flux-site and global scales, Environ. Res. Lett., 15, 124050,, 2020. 

Calvete-Sogo, H., Gonzalez-Fernandez, I., Sanz, J., Elvira, S., Alonso, R., Garcia-Gomez, H., Ibanez-Ruiz, M. A., and Bermejo-Bermejo, V.: Heterogeneous responses to ozone and nitrogen alter the species composition of Mediterranean annual pastures, Oecologia, 181, 1055–1067,, 2016. 

Chen, C., Riley, W. J., Prentice, I. C., and Keenan, T. F.: CO2 fertilization of terrestrial photosynthesis inferred from site to global scales, P. Natl. Acad. Sci. USA, 119, e2115627119,, 2022. 

Chen, L., Dirmeyer, P. A., Guo, Z., and Schultz, N. M.: Pairing FLUXNET sites to validate model representations of land-use/land-cover change, Hydrol. Earth Syst. Sci., 22, 111–125,, 2018. 

Chen, M., Rafique, R., Asrar, G. R., Bond-Lamberty, B., Ciais, P., and Zhao, F.; Regional contribution to variability and trends of global gross primary productivity. Environ. Res. Lett., 12, 105005,, 2017. 

Cheng, W., Dan, L., Deng, X., Feng, J., Wang, Y., Peng, J., Tian, J., Qi, W., Liu, Z., Zheng, X., Zhou, D., Jiang, S., Zhao, H., and Wang, X.: Global monthly gridded atmospheric carbon dioxide concentrations under the historical and future scenarios, Sci. Data, 9, 83,, 2022. 

Chopin, P., Bergkvist, G., and Hossard, L.: Modelling biodiversity change in agricultural landscape scenarios – A review and prospects for future research, Biol. Conserv., 235, 1–17,, 2019. 

Ciais, P., Reichstein, M., Viovy, N., Granier, A., Ogée, J., Allard, V., Aubinet, M., Buchmann, N., Bernhofer, C., Carrara, A., Chevallier, F., De Noblet, N., Friend, A. D., Friedlingstein, P., Grünwald, T., Heinesch, B., Keronen, P., Knohl, A., Krinner, G., Loustau, D., Manca, G., Matteucci, G., Miglietta, F., Ourcival, J. M., Papale, D., Pilegaard, K., Rambal, S., Seufert, G., Soussana, J. F., Sanz, M. J., Schulze, E. D., Vesala, T., and Valentini, R.: Europe-wide reduction in primary productivity caused by the heat and drought in 2003, Nature, 437, 529–533,, 2005. 

Clifton, O. E., Fiore, A. M., Massman, W. J., Baublitz, C. B., Coyle, M., Emberson, L., Fares, S., Farmer, D. K., Gentine, P., Gerosa, G., Guenther, A. B., Helmig, D., Lombardozzi, D. L., Munger, J. W., Patton, E. G., Pusede, S. E., Schwede, D. B., Silva, S. J., Sörgel, M., Steiner, A. L., and Tai, A. P. K.: Dry Deposition of Ozone Over Land: Processes, Measurement, and Modeling, Rev. Geophys., 58, e2019RG000670,, 2020. 

Collatz, G. J., Ball, J. T., Grivet, C., and Berry, J. A.: Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer, Agr. Forest Meteorol., 54, 107–136,, 1991. 

Collatz, G. J., Ribas-Carbo, M., and Berry, J.: Coupled Photosynthesis-Stomatal Conductance Model for Leaves of C4 Plants, Funct. Plant Biol., 19, 519–538,, 1992. 

Corcoran, E., Afshar, M., Curceac, S., Lashkari, A., Raza, M. M., Ahnert, S., Mead, A., and Morris, R.: Current data and modeling bottlenecks for predicting crop yields in the United Kingdom, Front. Sustain. Food Syst., 7, 1023169,, 2023. 

Danabasoglu, G., Lamarque, J. F., Bacmeister, J., Bailey, D. A., DuVivier, A. K., Edwards, J., Emmons, L. K., Fasullo, J., Garcia, R., Gettelman, A., Hannay, C., Holland, M. M., Large, W. G., Lauritzen, P. H., Lawrence, D. M., Lenaerts, J. T. M., Lindsay, K., Lipscomb, W. H., Mills, M. J., Neale, R., Oleson, K. W., Otto-Bliesner, B., Phillips, A. S., Sacks, W., Tilmes, S., van Kampenhout, L., Vertenstein, M., Bertini, A., Dennis, J., Deser, C., Fischer, C., Fox-Kemper, B., Kay, J. E., Kinnison, D., Kushner, P. J., Larson, V. E., Long, M. C., Mickelson, S., Moore, J. K., Nienhouse, E., Polvani, L., Rasch, P. J., and Strand, W. G.: The Community Earth System Model Version 2 (CESM2), J. Adv. Model. Earth Sy., 12, e2019MS001916,, 2020. 

Davies-Barnard, T., Zaehle, S., and Friedlingstein, P.: Assessment of the impacts of biological nitrogen fixation structural uncertainty in CMIP6 earth system models, Biogeosciences, 19, 3491–3503,, 2022. 

Dermody, O., Long, S. P., McConnaughay, K., and DeLucia, E. H.: How do elevated CO2 and O3 affect the interception and utilization of radiation by a soybean canopy?, Glob. Change Biol., 14, 556–564,, 2008. 

Deryng, D., Elliott, J., Folberth, C., Müller, C., Pugh, T. A. M., Boote, K. J., Conway, D., Ruane, A. C., Gerten, D., Jones, J. W., Khabarov, N., Olin, S., Schaphoff, S., Schmid, E., Yang, H., and Rosenzweig, C.: Regional disparities in the beneficial effects of rising CO2 concentrations on crop water productivity, Nat. Clim. Change, 6, 786–790,, 2016. 

Dickinson, R. E.: Land Surface Processes and Climate – Surface Albedos and Energy Balance, Adv. Geophys., 25, 305–353,, 1983. 

Dlugokencky, E. and Tans, P.: Trends in atmospheric carbon dioxide, National Oceanic and Atmospheric Administration, Earth System Research Laboratory (NOAA/ESRL), (last access: 21 March 2024), 2022. 

Dušek, J., Èížková, H., Stellner, S., Czerný, R., and Kvìt, J.: Fluctuating water table affects gross ecosystem production and gross radiation use efficiency in a sedge-grass marsh, Hydrobiologia, 692, 57–66,, 2012. 

Ebi, K. L., Anderson, C. L., Hess, J. J., Kim, S.-H., Loladze, I., Neumann, R. B., Singh, D., Ziska, L., and Wood, R.: Nutritional quality of crops in a high CO2 world: an agenda for research and technology development, Environ. Res. Lett., 16, 064045,, 2021. 

Emberson, L.: Effects of ozone on agriculture, forests and grasslands, Philos T. Roy. A-Math., 378, 20190327,, 2020. 

Emberson, L. D., Kitwiroon, N., Beevers, S., Büker, P., and Cinderby, S.: Scorched Earth: how will changes in the strength of the vegetation sink to ozone deposition affect human health and ecosystems?, Atmos. Chem. Phys., 13, 6741–6755,, 2013. 

Emberson, L. D., Pleijel, H., Ainsworth, E. A., van den Berg, M., Ren, W., Osborne, S., Mills, G., Pandey, D., Dentener, F., Büker, P., Ewert, F., Koeble, R., and Van Dingenen, R.: Ozone effects on crops and consideration in crop models, Eur. J. Agron., 100, 19–34,, 2018. 

Fares, S., Vargas, R., Detto, M., Goldstein, A. H., Karlik, J., Paoletti, E., and Vitale, M.: Tropospheric ozone reduces carbon assimilation in trees: estimates from analysis of continuous flux measurements, Glob. Change Biol., 19, 2427–2443,, 2013. 

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

Felzer, B. S., Cronin, T., Reilly, J. M., Melillo, J. M., and Wang, X.: Impacts of ozone on trees and crops, CR Geosci., 339, 784–798,, 2007. 

Feng, H., Guo, J., Peng, C., Kneeshaw, D., Roberge, G., Pan, C., Ma, X., Zhou, D., and Wang, W.: Nitrogen addition promotes terrestrial plants to allocate more biomass to aboveground organs: A global meta-analysis, Glob. Chang Biol., 29, 3970–3989,, 2023. 

Feng, Z., Kobayashi, K., and Ainsworth, E. A.: Impact of elevated ozone concentration on growth, physiology, and yield of wheat (Triticum aestivum L.): a meta-analysis, Glob. Change Biol., 14, 2696–2708,, 2008. 

Franklin, O., Harrison, S. P., Dewar, R., Farrior, C. E., Brannstrom, A., Dieckmann, U., Pietsch, S., Falster, D., Cramer, W., Loreau, M., Wang, H., Makela, A., Rebel, K. T., Meron, E., Schymanski, S. J., Rovenskaya, E., Stocker, B. D., Zaehle, S., Manzoni, S., van Oijen, M., Wright, I. J., Ciais, P., van Bodegom, P. M., Penuelas, J., Hofhansl, F., Terrer, C., Soudzilovskaia, N. A., Midgley, G., and Prentice, I. C.: Organizing principles for vegetation dynamics, Nat. Plants, 6, 444–453,, 2020. 

Franks, P. J., Adams, M. A., Amthor, J. S., Barbour, M. M., Berry, J. A., Ellsworth, D. S., Farquhar, G. D., Ghannoum, O., Lloyd, J., McDowell, N., Norby, R. J., Tissue, D. T., and von Caemmerer, S.: Sensitivity of plants to changing atmospheric CO2 concentration: from the geological past to the next century, New Phytol., 197, 1077–1094,, 2013. 

Franks, P. J., Berry, J. A., Lombardozzi, D. L., and Bonan, G. B.: Stomatal Function across Temporal and Spatial Scales: Deep-Time Trends, Land-Atmosphere Coupling and Global Models, Plant Physiol., 174, 583–602,, 2017. 

Franz, M. and Zaehle, S.: Competing effects of nitrogen deposition and ozone exposure on northern hemispheric terrestrial carbon uptake and storage, 1850–2099, Biogeosciences, 18, 3219–3241,, 2021. 

Franz, M., Simpson, D., Arneth, A., and Zaehle, S.: Development and evaluation of an ozone deposition scheme for coupling to a terrestrial biosphere model, Biogeosciences, 14, 45–71,, 2017. 

Fuhrer, J., Val Martin, M., Mills, G., Heald, C. L., Harmens, H., Hayes, F., Sharps, K., Bender, J., and Ashmore, M. R.: Current and future ozone risks to global terrestrial biodiversity and ecosystem processes, Ecol. Evol., 6, 8785–8799,, 2016. 

Gang, C., Wang, Z., You, Y., Liu, Y., Xu, R., Bian, Z., Pan, N., Gao, X., Chen, M., and Zhang, M.: Divergent responses of terrestrial carbon use efficiency to climate variation from 2000 to 2018, Global Planet. Change, 208, 103709,, 2022. 

Ganzeveld, L. and Lelieveld, J.: Dry deposition parameterization in a chemistry general circulation model and its influence on the distribution of reactive trace gases, J. Geophys. Res., 100, 20999,, 1995. 

Ganzeveld, L., Bouwman, L., Stehfest, E., van Vuuren, D. P., Eickhout, B., and Lelieveld, J.: Impact of future land use and land cover changes on atmospheric chemistry-climate interactions, J. Geophys. Res., 115, D23301,, 2010. 

Gelaro, R., McCarty, W., Suárez, 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. 

Gerber, S., Hedin, L. O., Oppenheimer, M., Pacala, S. W., and Shevliakova, E.: Nitrogen cycling and feedbacks in a global dynamic land model, Global Biogeochem. Cy., 24, GB1001,, 2010. 

Gleason, S. M., Barnard, D. M., Green, T. R., Mackay, S., Wang, D. R., Ainsworth, E. A., Altenhofen, J., Brodribb, T. J., Cochard, H., Comas, L. H., Cooper, M., Creek, D., DeJonge, K. C., Delzon, S., Fritschi, F. B., Hammer, G., Hunter, C., Lombardozzi, D., Messina, C. D., Ocheltree, T., Stevens, B. M., Stewart, J. J., Vadez, V., Wenz, J., Wright, I. J., Yemoto, K., and Zhang, H.: Physiological trait networks enhance understanding of crop growth and water use in contrasting environments, Plant Cell Environ., 45, 2554–2572,, 2022. 

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. 

Halladay, K. and Good, P.: Non-linear interactions between CO2 radiative and physiological effects on Amazonian evapotranspiration in an Earth system model, Clim. Dynam., 49, 2471–2490,, 2017. 

Hardacre, C., Wild, O., and Emberson, L.: An evaluation of ozone dry deposition in global scale chemistry climate models, Atmos. Chem. Phys., 15, 6419–6436,, 2015. 

Harrison, S. P., Cramer, W., Franklin, O., Prentice, I. C., Wang, H., Brannstrom, A., de Boer, H., Dieckmann, U., Joshi, J., Keenan, T. F., Lavergne, A., Manzoni, S., Mengoli, G., Morfopoulos, C., Penuelas, J., Pietsch, S., Rebel, K. T., Ryu, Y., Smith, N. G., Stocker, B. D., and Wright, I. J.: Eco-evolutionary optimality as a means to improve vegetation and land-surface models, New Phytol., 231, 2125–2141,, 2021. 

Haverd, V., Raupach, M. R., Briggs, P. R., Canadell, J. G., Isaac, P., Pickett-Heaps, C., Roxburgh, S. H., van Gorsel, E., Viscarra Rossel, R. A., and Wang, Z.: Multiple observation types reduce uncertainty in Australia's terrestrial carbon and water cycles, Biogeosciences, 10, 2011–2040,, 2013. 

He, P., Ma, X., and Sun, Z.: Interannual variability in summer climate change controls GPP long-term changes, Environ. Res., 212, 113409,, 2022. 

Hidy, D., Barcza, Z., Hollós, R., Dobor, L., Ács, T., Zacháry, D., Filep, T., Pásztor, L., Incze, D., Dencső, M., Tóth, E., Merganičová, K., Thornton, P., Running, S., and Fodor, N.: Soil-related developments of the Biome-BGCMuSo v6.2 terrestrial ecosystem model, Geosci. Model Dev., 15, 2157–2181,, 2022. 

Hoesly, R. M., Smith, S. J., Feng, L., Klimont, Z., Janssens-Maenhout, G., Pitkanen, T., Seibert, J. J., Vu, L., Andres, R. J., Bolt, R. M., Bond, T. C., Dawidowski, L., Kholod, N., Kurokawa, J.-I., Li, M., Liu, L., Lu, Z., Moura, M. C. P., O'Rourke, P. R., and Zhang, Q.: Historical (1750–2014) anthropogenic emissions of reactive gases and aerosols from the Community Emissions Data System (CEDS), Geosci. Model Dev., 11, 369–408,, 2018. 

Hou, H., Zhou, B. B., Pei, F., Hu, G., Su, Z., Zeng, Y., Zhang, H., Gao, Y., Luo, M., and Li, X.: Future Land Use/Land Cover Change Has Nontrivial and Potentially Dominant Impact on Global Gross Primary Productivity, Earth's Future, 10, e2021EF002628,, 2022. 

Hu, L., Keller, C. A., Long, M. S., Sherwen, T., Auer, B., Da Silva, A., Nielsen, J. E., Pawson, S., Thompson, M. A., Trayanov, A. L., Travis, K. R., Grange, S. K., Evans, M. J., and Jacob, D. J.: Global simulation of tropospheric chemistry at 12.5 km resolution: performance and evaluation of the GEOS-Chem chemical module (v10-1) within the NASA GEOS Earth system model (GEOS-5 ESM), Geosci. Model Dev., 11, 4603–4620,, 2018. 

Hu, Y., Schäfer, K. V. R., Zhu, L., Zhao, P., Zhao, X., Ni, G., Zhang, Y., Ye, H., Zhao, W., Shen, W., and Fu, S.: Impacts of Canopy and Understory Nitrogen Additions on Stomatal Conductance and Carbon Assimilation of Dominant Tree Species in a Temperate Broadleaved Deciduous Forest, Ecosystems, 24, 1468–1484,, 2021. 

Huntingford, C., Oliver, R. J., Mercado, L. M., and Sitch, S.: Technical note: A simple theoretical model framework to describe plant stomatal “sluggishness” in response to elevated ozone concentrations, Biogeosciences, 15, 5415–5422,, 2018. 

Hutley, L. B., Beringer, J., Isaac, P. R., Hacker, J. M., and Cernusak, L. A.: A sub-continental scale living laboratory: Spatial patterns of savanna vegetation over a rainfall gradient in northern Australia, Agr. Forest Meteorol., 151, 1417–1428,, 2011. 

Im, U., Brandt, J., Geels, C., Hansen, K. M., Christensen, J. H., Andersen, M. S., Solazzo, E., Kioutsioukis, I., Alyuz, U., Balzarini, A., Baro, R., Bellasio, R., Bianconi, R., Bieser, J., Colette, A., Curci, G., Farrow, A., Flemming, J., Fraser, A., Jimenez-Guerrero, P., Kitwiroon, N., Liang, C.-K., Nopmongcol, U., Pirovano, G., Pozzoli, L., Prank, M., Rose, R., Sokhi, R., Tuccella, P., Unal, A., Vivanco, M. G., West, J., Yarwood, G., Hogrefe, C., and Galmarini, S.: Assessment and economic valuation of air pollution impacts on human health over Europe and the United States as calculated by a multi-model ensemble in the framework of AQMEII3, Atmos. Chem. Phys., 18, 5967–5989,, 2018. 

Imer, D., Merbold, L., Eugster, W., and Buchmann, N.: Temporal and spatial variations of soil CO2, CH4 and N2O fluxes at three differently managed grasslands, Biogeosciences, 10, 5931–5945,, 2013. 

Jain, A., Yang, X., Kheshgi, H., McGuire, A. D., Post, W., and Kicklighter, D.: Nitrogen attenuation of terrestrial carbon cycle response to global environmental factors, Global Biogeochem. Cy., 23, GB4028,, 2009. 

Jin, Z., Yan, D., Zhang, Z., Li, M., Wang, T., Huang, X., Xie, M., Li, S., and Zhuang, B.: Effects of Elevated Ozone Exposure on Regional Meteorology and Air Quality in China Through Ozone-Vegetation Coupling, J. Geophys. Res.-Atmos., 128, e2022JD038119,, 2023. 

Karlsson, P. E., Uddling, J., Braun, S., Broadmeadow, M., Elvira, S., Gimeno, B. S., Le Thiec, D., Oksanen, E., Vandermeiren, K., Wilkinson, M., and Emberson, L.: New critical levels for ozone effects on young trees based on AOT40 and simulated cumulative leaf uptake of ozone, Atmos. Environ., 38, 2283–2294,, 2004. 

Karlsson, P. E., Klingberg, J., Engardt, M., Andersson, C., Langner, J., Karlsson, G. P., and Pleijel, H.: Past, present and future concentrations of ground-level ozone and potential impacts on ecosystems and human health in northern Europe, Sci. Total Environ., 576, 22–35,, 2017. 

Keller, C. A., Long, M. S., Yantosca, R. M., Da Silva, A. M., Pawson, S., and Jacob, D. J.: HEMCO v1.0: a versatile, ESMF-compliant component for calculating emissions in atmospheric models, Geosci. Model Dev., 7, 1409–1417,, 2014. 

Kimmins, J. P., Blanco, J. A., Seely, B., Welham, C., and Scoullar, K.: Complexity in modelling forest ecosystems: How much is enough?, Forest Ecol. Manage., 256, 1646–1658,, 2008. 

Kou-Giesbrecht, S. and Arora, V. K.: Compensatory Effects Between CO2, Nitrogen Deposition, and Nitrogen Fertilization in Terrestrial Biosphere Models Without Nitrogen Compromise Projections of the Future Terrestrial Carbon Sink, Geophys. Res. Lett., 50, e2022GL102618,, 2023. 

Lai, J., Lortie, C. J., Muenchen, R. A., Yang, J., and Ma, K.: Evaluating the popularity of R in ecology, Ecosphere, 10, e02567,, 2019. 

Lam, T. and Tai, A. P. K.: Simulating Gross Primary Productivity of vegetation under soil water stress using in-situ and reanalysis soil moisture data inputs, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-21397,, 2020. 

Lam, J. C. Y., Tai, A. P. K., Ducker, J. A., and Holmes, C. D.: Development of an ecophysiology module in the GEOS-Chem chemical transport model version 12.2.0 to represent biosphere–atmosphere fluxes relevant for ozone air quality, Geosci. Model Dev., 16, 2323–2342,, 2023. 

Lasslop, G., Reichstein, M., Papale, D., Richardson, A. D., Arneth, A., Barr, A., Stoy, P., and Wohlfahrt, G.: Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: critical issues and global evaluation, Glob. Change Biol., 16, 187–208,, 2010. 

Lawrence, D., Rosie, F., Charles, K., Keith, O., Sean, S., Mariana, V., Ben, A., Gordon, B., Bardan, G., Leo van, K., Daniel, K., Erik, K., Ryan, K., Peter, L., Fang, L., Hongyi, L., Danica, L., Yaqiong, L., Justin, P., William, R., William, S., Mingjie, S., Will, W., Chonggang, X., Ashehad, A., Andrew, B., Gautam, B., Patrick, B., Michael, B., Jonathan, B., Martyn, C., Tony, C., Kyla, D., Beth, D., Louisa, E., Josh, F., Mark, F., Pierre, G., Jan, L., Sam, L., Leung, L. R., William, L., Jon, P., Daniel, M. R., Ben, S., Jacquelyn, S., Andrew, S., Zachary, S., Jinyun, T., Ahmed, T., Quinn, T., Simone, T., Francis, V., and Xubin, Z.: Technical Description of version 5.0 of the Community Land Model (CLM), (last access: 21 March 2024), 2020. 

Lawrence, P. J. and Chase, T. N.: Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0), J. Geophys. Res., 112, G01023,, 2007. 

Le Quéré, C., Moriarty, R., Andrew, R. M., Canadell, J. G., Sitch, S., Korsbakken, J. I., Friedlingstein, P., Peters, G. P., Andres, R. J., Boden, T. A., Houghton, R. A., House, J. I., Keeling, R. F., Tans, P., Arneth, A., Bakker, D. C. E., Barbero, L., Bopp, L., Chang, J., Chevallier, F., Chini, L. P., Ciais, P., Fader, M., Feely, R. A., Gkritzalis, T., Harris, I., Hauck, J., Ilyina, T., Jain, A. K., Kato, E., Kitidis, V., Klein Goldewijk, K., Koven, C., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lenton, A., Lima, I. D., Metzl, N., Millero, F., Munro, D. R., Murata, A., Nabel, J. E. M. S., Nakaoka, S., Nojiri, Y., O'Brien, K., Olsen, A., Ono, T., Pérez, F. F., Pfeil, B., Pierrot, D., Poulter, B., Rehder, G., Rödenbeck, C., Saito, S., Schuster, U., Schwinger, J., Séférian, R., Steinhoff, T., Stocker, B. D., Sutton, A. J., Takahashi, T., Tilbrook, B., van der Laan-Luijkx, I. T., van der Werf, G. R., van Heuven, S., Vandemark, D., Viovy, N., Wiltshire, A., Zaehle, S., and Zeng, N.: Global Carbon Budget 2015, Earth Syst. Sci. Data, 7, 349–396,, 2015. 

Lee, E., Zeng, F.-W., Koster, R. D., Weir, B., Ott, L. E., and Poulter, B.: The impact of spatiotemporal variability in atmospheric CO2 concentration on global terrestrial carbon fluxes, Biogeosciences, 15, 5635–5652,, 2018. 

Legates, D. R. and McCabe, G. J.: Evaluating the use of “goodness-of-fit” Measures in hydrologic and hydroclimatic model validation, Water Resour. Res., 35, 233–241,, 1999. 

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. 

Leung, F., Williams, K., Sitch, S., Tai, A. P. K., Wiltshire, A., Gornall, J., Ainsworth, E. A., Arkebauer, T., and Scoby, D.: Calibrating soybean parameters in JULES 5.0 from the US-Ne2/3 FLUXNET sites and the SoyFACE-O3 experiment, Geosci. Model Dev., 13, 6201–6213,, 2020. 

Leung, F., Sitch, S., Tai, A. P. K., Wiltshire, A. J., Gornall, J. L., Folberth, G. A., and Unger, N.: CO2 fertilization of crops offsets yield losses due to future surface ozone damage and climate change, Environ. Res. Lett., 17, 074007,, 2022. 

Li, J., Mahalov, A., and Hyde, P.: Simulating the impacts of chronic ozone exposure on plant conductance and photosynthesis, and on the regional hydroclimate using WRF/Chem, Environ. Res. Lett., 11, 114017,, 2016. 

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., 11, 2563,, 2019. 

Li, X., Xiao, J., He, B., Altaf Arain, M., Beringer, J., Desai, A. R., Emmel, C., Hollinger, D. Y., Krasnova, A., Mammarella, I., Noe, S. M., Ortiz, P. S., Rey-Sanchez, A. C., Rocha, A. V., and Varlagin, A.: Solar-induced chlorophyll fluorescence is strongly correlated with terrestrial photosynthesis for a wide variety of biomes: First global analysis based on OCO-2 and flux tower observations, Glob. Change Biol., 24, 3990–4008,, 2018. 

Lin, H., Jacob, D. J., Lundgren, E. W., Sulprizio, M. P., Keller, C. A., Fritz, T. M., Eastham, S. D., Emmons, L. K., Campbell, P. C., Baker, B., Saylor, R. D., and Montuoro, R.: Harmonized Emissions Component (HEMCO) 3.0 as a versatile emissions component for atmospheric models: application in the GEOS-Chem, NASA GEOS, WRF-GC, CESM2, NOAA GEFS-Aerosol, and NOAA UFS models, Geosci. Model Dev., 14, 5487–5506,, 2021. 

Lin, S., Hu, Z., Wang, Y., Chen, X., He, B., Song, Z., Sun, S., Wu, C., Zheng, Y., Xia, X., Liu, L., Tang, J., Sun, Q., Joos, F., and Yuan, W.: Underestimated Interannual Variability of Terrestrial Vegetation Production by Terrestrial Ecosystem Models, Global Biogeochem. Cy., 37, e2023GB007696,, 2023. 

Lindauer, M., Schmid, H. P., Grote, R., Mauder, M., Steinbrecher, R., and Wolpert, B.: Net ecosystem exchange over a non-cleared wind-throw-disturbed upland spruce forest – Measurements and simulations, Agr. Forest Meteorol., 197, 219–234,, 2014. 

Liu, S., Yan, Z., Wang, Z., Serbin, S., Visser, M., Zeng, Y., Ryu, Y., Su, Y., Guo, Z., Song, G., Wu, Q., Zhang, H., Cheng, K. H., Dong, J., Hau, B. C. H., Zhao, P., Yang, X., Liu, L., Rogers, A., and Wu, J.: Mapping foliar photosynthetic capacity in sub-tropical and tropical forests with UAS-based imaging spectroscopy: Scaling from leaf to canopy, Remote Sens. Environ., 293, 113612,, 2023. 

Liu, X., Tai, A. P. K., and Fung, K. M.: Responses of surface ozone to future agricultural ammonia emissions and subsequent nitrogen deposition through terrestrial ecosystem changes, Atmos. Chem. Phys., 21, 17743–17758,, 2021. 

Liu, Y., Xiao, J., Ju, W., Zhu, G., Wu, X., Fan, W., Li, D., and Zhou, Y.: Satellite-derived LAI products exhibit large discrepancies and can lead to substantial uncertainty in simulated carbon and water fluxes, Remote Sens. Environ., 206, 174–188,, 2018. 

Lombardozzi, D., Levis, S., Bonan, G., and Sparks, J. P.: Predicting photosynthesis and transpiration responses to ozone: decoupling modeled photosynthesis and stomatal conductance, Biogeosciences, 9, 3113–3130,, 2012. 

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. 

Loveland, T. R. and Belward, A. S.: The IGBP-DIS global 1 km land cover data set, DISCover: First results, Int. J. Remote Sens., 18, 3289–3295,, 1997. 

Lu, X., Gilliam, F. S., Yue, X., Wang, B., and Kuang, Y.: Shifts in Above- Versus Below-Ground Carbon Gains to Terrestrial Ecosystems Carbon Sinks Under Excess Nitrogen Inputs, Global Biogeochem. Cy., 37, e2022GB007638,, 2023. 

Marcolla, B., Cescatti, A., Manca, G., Zorer, R., Cavagna, M., Fiora, A., Gianelle, D., Rodeghiero, M., Sottocornola, M., and Zampedri, R.: Climatic controls and ecosystem responses drive the inter-annual variability of the net ecosystem exchange of an alpine meadow, Agr. Forest Meteorol., 151, 1233–1243,, 2011. 

Mason, R. E., Craine, J. M., Lany, N. K., Jonard, M., Ollinger, S. V., Groffman, P. M., Fulweiler, R. W., Angerer, J., Read, Q. D., Reich, P. B., Templer, P. H., and Elmore, A. J.: Evidence, causes, and consequences of declining nitrogen availability in terrestrial ecosystems, Science, 376, eabh3767,, 2022. 

McLaughlin, S. B., Wullschleger, S. D., Sun, G., and Nosal, M.: Interactive effects of ozone and climate on water use, soil moisture content and streamflow in a southern Appalachian forest in the USA, New Phytol., 174, 125–136,, 2007. 

Medlyn, B. E., Duursma, R. A., Eamus, D., Ellsworth, D. S., Prentice, I. C., Barton, C. V. M., Crous, K. Y., De Angelis, P., Freeman, M., and Wingate, L.: Reconciling the optimal and empirical approaches to modelling stomatal conductance, Glob. Change Biol., 17, 2134–2144,, 2011. 

Merbold, L., Eugster, W., Stieger, J., Zahniser, M., Nelson, D., and Buchmann, N.: Greenhouse gas budget (CO2, CH4 and N2O) of intensively managed grassland following restoration, Glob. Change Biol., 20, 1913–1928,, 2014. 

Miner, G. L., Bauerle, W. L., and Baldocchi, D. D.: Estimating the sensitivity of stomatal conductance to photosynthesis: a review: The sensitivity of conductance to photosynthesis, Plant Cell Environ., 40, 1214–1238,, 2017. 

Monfreda, C., Ramankutty, N., and Foley, J. A.: Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000, Global Biogeochem. Cy., 22, GB1022,, 2008. 

Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Tr. Akad. Nauk SSSR Geofiz Inst., 24, 163–187, 1954. 

Monks, P. S., Archibald, A. T., Colette, A., Cooper, O., Coyle, M., Derwent, R., Fowler, D., Granier, C., Law, K. S., Mills, G. E., Stevenson, D. S., Tarasova, O., Thouret, V., von Schneidemesser, E., Sommariva, R., Wild, O., and Williams, M. L.: Tropospheric ozone and its precursors from the urban to the global scale from air quality to short-lived climate forcer, Atmos. Chem. Phys., 15, 8889–8973,, 2015. 

Moura, B. B., Alves, E. S., Marabesi, M. A., de Souza, S. R., Schaub, M., and Vollenweider, P.: Ozone affects leaf physiology and causes injury to foliage of native tree species from the tropical Atlantic Forest of southern Brazil, Sci. Total Environ., 610–611, 912–925,, 2018. 

Muller, B. and Martre, P.: Plant and crop simulation models: powerful tools to link physiology, genetics, and phenomics, J. Exp. Bot., 70, 2339–2344,, 2019. 

Myers, S. S., Zanobetti, A., Kloog, I., Huybers, P., Leakey, A. D. B., Bloom, A. J., Carlisle, E., Dietterich, L. H., Fitzgerald, G., Hasegawa, T., Holbrook, N. M., Nelson, R. L., Ottman, M. J., Raboy, V., Sakai, H., Sartor, K. A., Schwartz, J., Seneweera, S., Tausz, M., and Usui, Y.: Increasing CO2 threatens human nutrition, Nature, 510, 139–142,, 2014. 

Noormets, A., Kull, O., Sôber, A., Kubiske, M. E., and Karnosky, D. F.: Elevated CO2 response of photosynthesis depends on ozone concentration in aspen, Environ. Pollut., 158, 992–999,, 2010. 

Norman, J. M.: Simulation of microclimates, Biometeorology and Integrated Pest Management, 65–99,, 1982. 

Nowak, D., Hirabayashi, S., Bodine, A., and Greenfield, E.: Tree and forest effects on air quality and human health in the United States, Environ. Pollut., 193, 119–129,, 2014. 

Oikawa, S. and Ainsworth, E. A.: Changes in leaf area, nitrogen content and canopy photosynthesis in soybean exposed to an ozone concentration gradient, Environ. Pollut., 215, 347–355,, 2016. 

Oleson, K. W., Lawrence, D. M., Bonan, G. B., Drewniak, B., Huang, M., Levis, S., Li, F., Riley, W. J., Swenson, S. C., Thornton, P. E., Bozbiyik, A., Fisher, R., Heald, C. L., Kluzek, E., Lamarque, F., Lawrence, P. J., Leung, L. R., Muszala, S., Ricciuto, D. M., Sacks, W., Sun, Y., Tang, J., and Yang, Z.-L.: Technical Description of version 4.5 of the Community Land Model (CLM),, 2013. 

Oliver, R. J., Mercado, L. M., Sitch, S., Simpson, D., Medlyn, B. E., Lin, Y.-S., and Folberth, G. A.: Large but decreasing effect of ozone on the European carbon sink, Biogeosciences, 15, 4245–4269,, 2018. 

O'Sullivan, M., Spracklen, D. V., Batterman, S. A., Arnold, S. R., Gloor, M., and Buermann, W.: Have Synergies Between Nitrogen Deposition and Atmospheric CO2 Driven the Recent Enhancement of the Terrestrial Carbon Sink?, Global Biogeochem. Cy., 33, 163–180,, 2019. 

Pacifico, F., Harrison, S. P., Jones, C. D., Arneth, A., Sitch, S., Weedon, G. P., Barkley, M. P., Palmer, P. I., Serça, D., Potosnak, M., Fu, T.-M., Goldstein, A., Bai, J., and Schurgers, G.: Evaluation of a photosynthesis-based biogenic isoprene emission scheme in JULES and simulation of isoprene emissions under present-day climate conditions, Atmos. Chem. Phys., 11, 4371–4389,, 2011. 

Pacifico, F., Folberth, G. A., Jones, C. D., Harrison, S. P., and Collins, W. J.: Sensitivity of biogenic isoprene emissions to past, present, and future environmental conditions and implications for atmospheric chemistry, J. Geophys. Res.-Atmos., 117, D22302,, 2012. 

Parrington, M., Jones, D. B. A., Bowman, K. W., Horowitz, L. W., Thompson, A. M., Tarasick, D. W., and Witte, J. C.: Estimating the summertime tropospheric ozone distribution over North America through assimilation of observations from the Tropospheric Emission Spectrometer, J. Geophys. Res., 113, D18307,, 2008.Please provide article number or page range. 

Peichl, M., Brodeur, J. J., Khomik, M., and Arain, M. A.: Biometric and eddy-covariance based estimates of carbon fluxes in an age-sequence of temperate pine forests, Agr. Forest Meteorol., 150, 952–965,, 2010. 

Pleijel, H., Danielsson, H., Ojanperä, K., Temmerman, L. D., Högy, P., Badiani, M., and Karlsson, P. E.: Relationships between ozone exposure and yield loss in European wheat and potato – a comparison of concentration- and flux-based exposure indices, Atmos. Environ., 38, 2259–2269,, 2004. 

Pongratz, J., Reick, C. H., Raddatz, T., and Claussen, M.: Effects of anthropogenic land cover change on the carbon cycle of the last millennium, Global Biogeochem. Cy., 23, GB4001,, 2009. 

Pongratz, J., Reick, C. H., Raddatz, T., and Claussen, M.: Biogeophysical versus biogeochemical climate response to historical anthropogenic land cover change, Geophys. Res. Lett., 37, L08702,, 2010. 

Porporato, A., Laio, F., Ridol, L., and Rodriguez-Iturbe, I.: Plants in water-controlled ecosystems: active role in hydrologic processes and response to water stress III. Vegetation water stress, Adv. Water Resour., 24, 725–744, 2001. 

Portmann, F. T., Siebert, S., and Döll, P.: MIRCA2000-Global monthly irrigated and rainfed crop areas around the year 2000: A new high-resolution data set for agricultural and hydrological modeling, Global Biogeochem. Cy., 24, GB1011,, 2010. 

Prescher, A.-K., Grünwald, T., and Bernhofer, C.: Land use regulates carbon budgets in eastern Germany: From NEE to NBP, Agr. Forest Meteorol., 150, 1016–1025,, 2010. 

Proietti, C., Anav, A., De Marco, A., Sicard, P., and Vitale, M.: A multi-sites analysis on the ozone effects on Gross Primary Production of European forests, Sci. Total Environ., 556, 1–11,, 2016. 

Pu, X., Wang, T. J., Huang, X., Melas, D., Zanis, P., Papanastasiou, D. K., and Poupkou, A.: Enhanced surface ozone during the heat wave of 2013 in Yangtze River Delta region, China, Sci. Total Environ., 603–604, 807–816,, 2017. 

Rahman, M. H. U., Ahrends, H. E., Raza, A., and Gaiser, T.: Current approaches for modeling ecosystem services and biodiversity in agroforestry systems: Challenges and ways forward, Frontiers in Forests and Global Change, 5, 1032442,, 2023. 

R Core Team: R: A language and environment for statistical computing, Vienna, Austria, 2022. 

Reichstein, M., Falge, E., Baldocchi, D., Papale, D., Aubinet, M., Berbigier, P., Bernhofer, C., Buchmann, N., Gilmanov, T., Granier, A., Grunwald, T., Havrankova, K., Ilvesniemi, H., Janous, D., Knohl, A., Laurila, T., Lohila, A., Loustau, D., Matteucci, G., Meyers, T., Miglietta, F., Ourcival, J.-M., Pumpanen, J., Rambal, S., Rotenberg, E., Sanz, M., Tenhunen, J., Seufert, G., Vaccari, F., Vesala, T., Yakir, D., and Valentini, R.: On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm, Glob. Change Biol., 11, 1424–1439,, 2005. 

Rhea, L., King, J., Kubiske, M., Saliendra, N., and Teclaw, R.: Effects of elevated atmospheric CO2 and tropospheric O3 on tree branch growth and implications for hydrologic budgeting, Environ. Pollut., 158, 1079–1087,, 2010. 

Roberts, H. R., Dodd, I. C., Hayes, F., and Ashworth, K.: Chronic tropospheric ozone exposure reduces seed yield and quality in spring and winter oilseed rape, Agr. Forest Meteorol., 316, 108859,, 2022. 

Sadiq, M., Tai, A. P. K., Lombardozzi, D., and Val Martin, M.: Effects of ozone–vegetation coupling on surface ozone air quality via biogeochemical and meteorological feedbacks, Atmos. Chem. Phys., 17, 3055–3066,, 2017. 

Sanderson, M. G., Collins, W. J., Hemming, D. L., and Betts, R. A.: Stomatal conductance changes due to increasing carbon dioxide levels: Projected impact on surface ozone levels, Tellus B, 59, 404-411,, 2007. 

Sandor, R., Ehrhardt, F., Brilli, L., Carozzi, M., Recous, S., Smith, P., Snow, V., Soussana, J. F., Dorich, C. D., Fuchs, K., Fitton, N., Gongadze, K., Klumpp, K., Liebig, M., Martin, R., Merbold, L., Newton, P. C. D., Rees, R. M., Rolinski, S., and Bellocchi, G.: The use of biogeochemical models to evaluate mitigation of greenhouse gas emissions from managed grasslands, Sci. Total Environ., 642, 292–306,, 2018. 

Schimel, D., Stephens, B. B., and Fisher, J. B.: Effect of increasing CO2 on the terrestrial carbon cycle, P. Natl. Acad. Sci. USA, 112, 436–441,, 2015. 

Scott, R. L., Biederman, J. A., Hamerlynck, E. P., and Barron-Gafford, G. A.: The carbon balance pivot point of southwestern U.S. semiarid ecosystems: Insights from the 21st century drought, J. Geophys. Res.-Biogeo., 120, 2612–2624,, 2015. 

Seiler, C., Melton, J. R., Arora, V. K., Sitch, S., Friedlingstein, P., Anthoni, P., Goll, D., Jain, A. K., Joetzjer, E., Lienert, S., Lombardozzi, D., Luyssaert, S., Nabel, J. E. M. S., Tian, H., Vuichard, N., Walker, A. P., Yuan, W., and Zaehle, S.: Are Terrestrial Biosphere Models Fit for Simulating the Global Land Carbon Sink?, J. Adv. Model. Earth Sy., 14, e2021MS002946,, 2022. 

Sellers, P. J.: Canopy reflectance, photosynthesis and transpiration, Int. J. Remote Sens., 6, 1335–1372,, 1985. 

Serrano-Ortiz, P., Domingo, F., Cazorla, A., Were, A., Cuezva, S., Villagarcía, L., Alados-Arboledas, L., and Kowalski, A. S.: Interannual CO2 exchange of a sparse Mediterranean shrubland on a carbonaceous substrate, J. Geophys. Res., 114, G04015,, 2009. 

Shekhar, A., Buchmann, N., and Gharun, M.: How well do recently reconstructed solar-induced fluorescence datasets model gross primary productivity?, Remote Sens. Environ., 283, 113282,, 2022. 

Sicard, P., Anav, A., De Marco, A., and Paoletti, E.: Projected global ground-level ozone impacts on vegetation under different emission and climate scenarios, Atmos. Chem. Phys., 17, 12177–12196,, 2017. 

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. 

Slevin, D., Tett, S. F. B., Exbrayat, J.-F., Bloom, A. A., and Williams, M.: Global evaluation of gross primary productivity in the JULES land surface model v3.4.1, Geosci. Model Dev., 10, 2651–2670,, 2017. 

Sun, G., McLaughlin, S. B., Porter, J. H., Uddling, J., Mulholland, P. J., Adams, M. B., and Pederson, N.: Interactive influences of ozone and climate on streamflow of forested watersheds, Glob. Change Biol., 18, 3395–3409,, 2012. 

Sun, S., Tai, A. P. K., Yung, D. H. Y., Wong, A. Y. H., Ducker, J. A., and Holmes, C. D.: Influence of plant ecophysiology on ozone dry deposition: comparing between multiplicative and photosynthesis-based dry deposition schemes and their responses to rising CO2 level, Biogeosciences, 19, 1753–1776,, 2022. 

Tai, A. P. K., Sadiq, M., Pang, J. Y. S., Yung, D. H. Y., and Feng, Z.: Impacts of Surface Ozone Pollution on Global Crop Yields: Comparing Different Ozone Exposure Metrics and Incorporating Co-effects of CO2, Front. Sustain. Food Syst., 5, 534616,, 2021. 

Tai, A. P. K. and Yung, D. H. Y.: TEMIR v1.0 Public Release (v1.0), Zenodo [code],, 2023a. 

Tai, A. P. K. and Yung, D. H. Y.: Input Data for TEMIR v1.0 (1.0), Zenodo [data set],, 2023b. 

Terrer, C., Phillips, R. P., Hungate, B. A., Rosende, J., Pett-Ridge, J., Craig, M. E., van Groenigen, K. J., Keenan, T. F., Sulman, B. N., Stocker, B. D., Reich, P. B., Pellegrini, A. F. A., Pendall, E., Zhang, H., Evans, R. D., Carrillo, Y., Fisher, J. B., Van Sundert, K., Vicca, S., and Jackson, R. B.: A trade-off between plant and soil carbon storage under elevated CO(2), Nature, 591, 599–603,, 2021. 

Tharammal, T., Bala, G., Narayanappa, D., and Nemani, R.: Potential roles of CO2 fertilization, nitrogen deposition, climate change, and land use and land cover change on the global terrestrial carbon uptake in the twenty-first century, Clim. Dynam., 52, 4393–4406,, 2018. 

Tian, J., Zhang, Y., and Zhang, X.: Impacts of heterogeneous CO2 on water and carbon fluxes across the global land surface, International J. Digit. Earth, 14, 1175–1193,, 2021. 

Ueyama, M., Ichii, K., Kobayashi, H., Kumagai, T. O., Beringer, J., Merbold, L., Euskirchen, E. S., Hirano, T., Marchesini, L. B., Baldocchi, D., Saitoh, T. M., Mizoguchi, Y., Ono, K., Kim, J., Varlagin, A., Kang, M., Shimizu, T., Kosugi, Y., Bret-Harte, M. S., Machimura, T., Matsuura, Y., Ohta, T., Takagi, K., Takanashi, S., and Yasuda, Y.: Inferring CO2 fertilization effect based on global monitoring land-atmosphere exchange with a theoretical model, Environ. Res. Lett., 15, 084009,, 2020. 

Val Martin, M., Heald, C. L., and Arnold, S. R.: Coupling dry deposition to vegetation phenology in the Community Earth System Model: Implications for the simulation of surface O3, Geophys. Res. Lett., 41, 2988–2996,, 2014. 

Van Dingenen, R., Dentener, F. J., Raes, F., Krol, M. C., Emberson, L., and Cofala, J.: The global impact of ozone on agricultural crop yields under current and future air quality legislation, Atmos. Environ., 43, 604–618,, 2009. 

Verbeke, T., Lathière, J., Szopa, S., and de Noblet-Ducoudré, N.: Impact of future land-cover changes on HNO3 and O3 surface dry deposition, Atmos. Chem. Phys., 15, 13555–13568,, 2015. 

Verhoef, A. and Egea, G.: Modeling plant transpiration under limited soil water: Comparison of different plant and soil hydraulic parameterizations and preliminary implications for their use in land surface models, Agr. Forest Meteorol., 191, 22–32,, 2014. 

von Caemmerer, S. and Farquhar, G. D.: Some relationships between the biochemistry of photosynthesis and the gas exchange of leaves, Planta, 153, 376–387,, 1981. 

Wang, S., Zhang, Y., Ju, W., Chen, J. M., Ciais, P., Cescatti, A., Sardans, J., Janssens, I. A., Wu, M., Berry, J. A., Campbell, E., Fernández-Martínez, M., Alkama, R., Sitch, S., Friedlingstein, P., Smith, W. K., Yuan, W., He, W., Lombardozzi, D., Kautz, M., Zhu, D., Lienert, S., Kato, E., Poulter, B., Sanders, T. G. M., Krüger, I., Wang, R., Zeng, N., Tian, H., Vuichard, N., Jain, A. K., Wiltshire, A., Haverd, V., Goll, D. S., and Peñuelas, J.: Recent global decline of CO2 fertilization effects on vegetation photosynthesis, Science, 370, 1295–1300,, 2020. 

Wang, X., Chen, J. M., Ju, W., and Zhang, Y.: Seasonal Variations in Leaf Maximum Photosynthetic Capacity and Its Dependence on Climate Factors Across Global FLUXNET Sites, J. Geophys. Res.-Biogeo., 127, e2021JG006709,, 2022. 

Wang, Y., Jacob, D. J., and Logan, J. A.: Global simulation of tropospheric O3-NOx-hydrocarbon chemistry: 3. Origin of tropospheric ozone and effects of nonmethane hydrocarbons, J. Geophys. Res.-Atmos., 103, 10757–10767,, 1998. 

Wesely, M. L.: Parameterization of surface resistances to gaseous dry deposition in regional-scale numerical models, Atmos. Environ., 23, 1293–1304, 1989. 

Wild, B., Teubner, I., Moesinger, L., Zotta, R.-M., Forkel, M., van der Schalie, R., Sitch, S., and Dorigo, W.: VODCA2GPP – a new, global, long-term (1988–2020) gross primary production dataset from microwave remote sensing, Earth Syst. Sci. Data, 14, 1063–1085,, 2022. 

Wild, O.: Modelling the global tropospheric ozone budget: exploring the variability in current models, Atmos. Chem. Phys., 7, 2643–2660,, 2007. 

Wiltshire, A. J., Burke, E. J., Chadburn, S. E., Jones, C. D., Cox, P. M., Davies-Barnard, T., Friedlingstein, P., Harper, A. B., Liddicoat, S., Sitch, S., and Zaehle, S.: JULES-CN: a coupled terrestrial carbon–nitrogen scheme (JULES vn5.1), Geosci. Model Dev., 14, 2161–2186,, 2021. 

Wolf, S., Keenan, T. F., Fisher, J. B., Baldocchi, D. D., Desai, A. R., Richardson, A. D., Scott, R. L., Law, B. E., Litvak, M. E., Brunsell, N. A., Peters, W., and van der Laan-Luijkx, I. T.: Warm spring reduced carbon cycle impact of the 2012 US summer drought, P. Natl. Acad. Sci. USA, 113, 5880–5885,, 2016. 

Wong, A. Y. H., Geddes, J. A., Tai, A. P. K., and Silva, S. J.: Importance of dry deposition parameterization choice in global simulations of surface ozone, Atmos. Chem. Phys., 19, 14365–14385,, 2019. 

Wu, Q., Chen, S., Zhang, Y., Song, C., Ju, W., Wang, L., and Jiang, J.: Improved Estimation of the Gross Primary Production of Europe by Considering the Spatial and Temporal Changes in Photosynthetic Capacity from 2001 to 2016, Remote Sens., 15, 1172,, 2023a. 

Wu, Q., Wang, X., Chen, S., Wang, L., and Jiang, J.: Land Surface Greening and CO2 Fertilization More than Offset the Gross Carbon Sequestration Decline Caused by Land Cover Change and the Enhanced Vapour Pressure Deficit in Europe, Remote Sens., 15, 1372,, 2023b. 

Xia, L., Lam, S. K., Kiese, R., Chen, D., Luo, Y., van Groenigen, K. J., Ainsworth, E. A., Chen, J., Liu, S., Ma, L., Zhu, Y., and Butterbach-Bahl, K.: Elevated CO2 negates O3 impacts on terrestrial carbon and nitrogen cycles, One Earth, 4, 1752–1763,, 2021. 

Xu, B., Arain, M. A., Black, T. A., Law, B. E., Pastorello, G. Z., and Chu, H.: Seasonal variability of forest sensitivity to heat and drought stresses: A synthesis based on carbon fluxes from North American forest ecosystems, Glob. Change Biol., 26, 901–918,, 2020. 

Yan, M., Yue, X., Zhou, B., Sun, X., and Xin, N.: Projected changes of ecosystem productivity and their responses to extreme heat events in northern asia, Front. Earth Sci., 10, 970296,, 2022. 

Yang, R., Wang, J., Zeng, N., Sitch, S., Tang, W., McGrath, M. J., Cai, Q., Liu, D., Lombardozzi, D., Tian, H., Jain, A. K., and Han, P.: Divergent historical GPP trends among state-of-the-art multi-model simulations and satellite-based products, Earth Syst. Dynam., 13, 833–849,, 2022. 

Yang, X., Wittig, V., Jain, A. K., and Post, W.: Integration of nitrogen cycle dynamics into the Integrated Science Assessment Model for the study of terrestrial ecosystem responses to global change, Global Biogeochem. Cy., 23, GB4029,, 2009. 

Yang, X., Chen, X., Ren, J., Yuan, W., Liu, L., Liu, J., Chen, D., Xiao, Y., Song, Q., Du, Y., Wu, S., Fan, L., Dai, X., Wang, Y., and Su, Y.: A gridded dataset of a leaf-age-dependent leaf area index seasonality product over tropical and subtropical evergreen broadleaved forests, Earth Syst. Sci. Data, 15, 2601–2622,, 2023. 

Yi, K., Dragoni, D., Phillips, R. P., Roman, D. T., and Novick, K. A.: Dynamics of stem water uptake among isohydric and anisohydric species experiencing a severe drought, Tree Physiol, 37, 1379–1392,, 2017. 

Yuan, H., Dai, Y., Xiao, Z., Ji, D., and Shangguan, W.: Reprocessing the MODIS Leaf Area Index products for land surface and climate modelling, Remote Sens. Environ., 115, 1171–1187,, 2011. 

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., 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. 

Zeng, X., Zhao, M., and Dickinson, R. E.: Intercomparison of Bulk Aerodynamic Algorithms for the Computation of Sea Surface Fluxes Using TOGA COARE and TAO Data, J. Climate, 11, 2628–2644,<2628:IOBAAF>2.0.CO;2, 1998. 

Zeng, X., Shaikh, M., Dai, Y., Dickinson, R. E., and Myneni, R.: Coupling of the Common Land Model to the NCAR Community Climate Model, J. Climate, 15, 1832–1854,<1, 2002. 

Zhang, L., Moran, M. D., Makar, P. A., Brook, J. R., and Gong, S.: Modelling gaseous dry deposition in AURAMS: a unified regional air-quality modelling system, Atmos. Environ., 36, 537–560,, 2002. 

Zhang, L., Brook, J. R., and Vet, R.: A revised parameterization for gaseous dry deposition in air-quality models, Atmos. Chem. Phys., 3, 2067–2082,, 2003. 

Zhang, L., Jacob, D. J., Liu, X., Logan, J. A., Chance, K., Eldering, A., and Bojkov, B. R.: Intercomparison methods for satellite measurements of atmospheric composition: application to tropospheric ozone from TES and OMI, Atmos. Chem. Phys., 10, 4725–4739,, 2010. 

Zhang, L., Hoshika, Y., Carrari, E., Burkey, K. O., and Paoletti, E.: Protecting the photosynthetic performance of snap bean under free air ozone exposure, J. Environ. Sci., 66, 31–40,, 2018. 

Zhang, X., Ward, B. B., and Sigman, D. M.: Global Nitrogen Cycle: Critical Enzymes, Organisms, and Processes for Nitrogen Budgets and Dynamics, Chem. Rev., 120, 5308–5351,, 2020. 

Zhang, Y. and Wang, Y.: Climate-driven ground-level ozone extreme in the fall over the Southeast United States, P. Natl. Acad. Sci. USA, 113, 10025–10030,, 2016. 

Zhang, Y. and Ye, A.: Improving global gross primary productivity estimation by fusing multi-source data products, Heliyon, 8, e09153,, 2022. 

Zhang, Y., Xiao, X., Wu, X., Zhou, S., Zhang, G., Qin, Y., and Dong, J.: A global moderate resolution dataset of gross primary production of vegetation for 2000–2016, Sci. Data, 4, 170165,, 2017. 

Zhao, J., Liu, D., Zhu, Y., Peng, H., and Xie, H.: A review of forest carbon cycle models on spatiotemporal scales, J. Clean. Prod., 339,, 2022. 

Zhao, Y., Zhang, L., Tai, A. P. K., Chen, Y., and Pan, Y.: Responses of surface ozone air quality to anthropogenic nitrogen deposition in the Northern Hemisphere, Atmos. Chem. Phys., 17, 9781–9796,, 2017.  

Zhou, S. S., Tai, A. P. K., Sun, S., Sadiq, M., Heald, C. L., and Geddes, J. A.: Coupling between surface ozone and leaf area index in a chemical transport model: strength of feedback and implications for ozone air quality and vegetation health, Atmos. Chem. Phys., 18, 14133–14148,, 2018. 

Zhu, J., Tai, A. P. K., and Hung Lam Yim, S.: Effects of ozone–vegetation interactions on meteorology and air quality in China using a two-way coupled land–atmosphere model, Atmos. Chem. Phys., 22, 765–782,, 2022. 

Zielis, S., Etzold, S., Zweifel, R., Eugster, W., Haeni, M., and Buchmann, N.: NEP of a Swiss subalpine forest is significantly driven not only by current but also by previous year's weather, Biogeosciences, 11, 1627–1635,, 2014. 

Short summary
We have developed the Terrestrial Ecosystem Model in R (TEMIR), which simulates plant carbon and pollutant uptake and predicts their response to varying atmospheric conditions. This model is designed to couple with an atmospheric chemistry model so that questions related to plant–atmosphere interactions, such as the effects of climate change, rising CO2, and ozone pollution on forest carbon uptake, can be addressed. The model has been well validated with both ground and satellite observations.