The Yale Interactive terrestrial Biosphere model version 1 . 0 : description , evaluation and implementation into NASA GISS

The land biosphere, atmospheric chemistry and climate are intricately interconnected, yet the modeling of carbon–climate and chemistry–climate interactions have evolved as entirely separate research communities. We describe the Yale Interactive terrestrial Biosphere (YIBs) model version 1.0, a land carbon cycle model that has been developed for coupling to the NASA Goddard Institute for Space Studies (GISS) ModelE2 global chemistry–climate model. The YIBs model adapts routines from the mature TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics) and CASA (Carnegie–Ames–Stanford Approach) models to simulate interactive carbon assimilation, allocation, and autotrophic and heterotrophic respiration. Dynamic daily leaf area index is simulated based on carbon allocation and temperatureand drought-dependent prognostic phenology. YIBs incorporates a semi-mechanistic ozone vegetation damage scheme. Here, we validate the present-day YIBs land carbon fluxes for three increasingly complex configurations: (i) offline local site level, (ii) offline global forced with WFDEI (WATCH Forcing Data methodology applied to ERA-Interim data) meteorology, and (iii) online coupled to the NASA ModelE2 (NASA ModelE2YIBs). Offline YIBs has hourly and online YIBs has halfhourly temporal resolution. The large observational database used for validation includes carbon fluxes from 145 flux tower sites and multiple satellite products. At the site level, YIBs simulates reasonable seasonality (correlation coefficient R > 0.8) of gross primary productivity (GPP) at 121 out of 145 sites with biases in magnitude ranging from −19 to 7 % depending on plant functional type. On the global scale, the offline model simulates an annual GPP of 125± 3 Pg C and net ecosystem exchange (NEE) of −2.5± 0.7 Pg C for 1982–2011, with seasonality and spatial distribution consistent with the satellite observations. We assess present-day global ozone vegetation damage using the offline YIBs configuration. Ozone damage reduces global GPP by 2–5 % annually with regional extremes of 4–10 % in east Asia. The online model simulates annual GPP of 123± 1 Pg C and NEE of −2.7± 0.7 Pg C. NASA ModelE2-YIBs is a useful new tool to investigate coupled interactions between the land carbon cycle, atmospheric chemistry, and climate change.


Introduction
The terrestrial biosphere interacts with the atmosphere through the exchanges of energy, carbon, reactive gases, water, and momentum fluxes.Forest ecosystems absorb an estimated 120 Pg C year −1 from the atmosphere (Beer et al., 2010) and mitigate about one-quarter of the anthropogenic carbon dioxide (CO 2 ) emissions (Friedlingstein et al., 2014).This carbon assimilation is sensitive to human-caused perturbations including climate change and land use change (Zhao and Running, 2010;Houghton et al., 2012) and is affected by atmospheric pollutants such as ozone and aerosols (Sitch et al., 2007;Mercado et al., 2009).Over the past 2-3 decades, a number of terrestrial biosphere models have been developed as tools to quantify the present-day global carbon budget in conjunction with available but sparse observations (e.g., Jung et al., 2009), to understand the relationships between terrestrial biospheric fluxes and environmental conditions (e.g., Zeng et al., 2005), to attribute drivers of trends in the carbon cycle during the anthropogenic era (e.g., Sitch et al., 2015), X. Yue and N. Unger: The Yale Interactive terrestrial Biosphere model version 1.0 and to project future changes in the land biosphere and the consequences for regional and global climate change (e.g., Friedlingstein et al., 2006).
Emerging research identifies climatically relevant interactions between the land biosphere and atmospheric chemistry (e.g, Huntingford et al., 2011).For instance, stomatal uptake is an important sink of tropospheric ozone (Val Martin et al., 2014) but damages photosynthesis, reduces plant growth and biomass accumulation, limits crop yields, and affects stomatal control over plant transpiration of water vapor between the leaf surface and atmosphere (Ainsworth et al., 2012;Hollaway et al., 2012).The indirect CO 2 radiative forcing due to the vegetation damage effects of anthropogenic ozone increases since the industrial revolution may be as large as +0.4 W m −2 (Sitch et al., 2007), which is 25 % of the magnitude of the direct CO 2 radiative forcing over the same period, and of similar magnitude to the direct ozone radiative forcing.Atmospheric oxidation of biogenic volatile organic compound (BVOC) emissions affects surface air quality and exerts additional regional and global chemical climate forcings (Scott et al., 2014;Unger, 2014a, b).Fine-mode atmospheric pollution particles affect the land biosphere by changing the physical climate state and through diffuse radiation fertilization (Mercado et al., 2009;Mahowald, 2011).Land plant phenology has experienced substantial changes in the last few decades (Keenan et al., 2014), possibly influencing both ozone deposition and BVOC emissions through the extension of growing seasons.These coupled interactions are often not adequately represented in current generation land biosphere models or global chemistry-climate models.Global land carbon cycle models often prescribe offline ozone and aerosol fields (e.g., Sitch et al., 2007;Mercado et al., 2009), and global chemistry-climate models often prescribe fixed offline vegetation fields (e.g., Lamarque et al., 2013;Shindell et al., 2013a).However, multiple mutual feedbacks occur between vegetation physiology and reactive atmospheric chemical composition that are completely neglected using these previous offline approaches.Model frameworks are needed that fully two-way couple the land carbon cycle and atmospheric chemistry and simulate the consequences for climate change.
Our objective is to present the description and presentday evaluation of the Yale Interactive terrestrial Biosphere (YIBs) model version 1.0 that has been developed for the investigation of carbon-chemistry-climate interactions.The YIBs model can be used in three configurations: (i) offline local site level, (ii) offline global forced with WFDEI (WATCH Forcing Data methodology applied to ERA-Interim data) meteorology, and (iii) online coupled to the latest frozen version of the NASA GISS (Goddard Institute for Space Studies) ModelE2 (Schmidt et al., 2014).The global climate model represents atmospheric gas-phase and aerosol chemistry, cloud, radiation, and land surface processes, and has been widely used for studies of atmospheric components, climate change, and their interactions (Schmidt et al., 2006;Koch et al., 2011;Unger, 2011;Shindell et al., 2013b;Miller et al., 2014).To our knowledge, this study represents the first description and validation of an interactive climate-sensitive closed land carbon cycle in NASA ModelE2.The impacts of the updated vegetation scheme on the chemistry and climate simulations in NASA ModelE2 will be addressed in other ongoing research.Section 2 describes the observational data sets used to evaluate YIBs land carbon cycle performance.Section 3 describes physical parameterizations of the vegetation model.Section 4 explains the model setup and simulations in three configurations.Section 5 presents the results of the model evaluation and Sect.6 summarizes the model performance.

YIBs design strategy
Many land carbon cycle models already exist (e.g., Sitch et al., 2015, and references therein;Schaefer et al., 2012, and references therein).We elected to build YIBs in a step-bystep process such that our research group has intimate familiarity with the underlying scientific processes, rather than adopting an existing model as a "black box".This unconventional interdisciplinary approach is important for discerning the complex mutual feedbacks between atmospheric chemistry and the land carbon sink under global change.The development of the YIBs land carbon cycle model has proceeded in three main steps.The first step was the implementation of vegetation biophysics, photosynthesis-dependent BVOC emissions and ozone vegetation damage that have been extensively documented, validated and applied in seven previous publications (Unger, 2013(Unger, , 2014a, b;, b;Unger et al., 2013;Unger and Yue, 2014;Yue and Unger, 2014;Zheng et al., 2015).The second step was the selection of the YIBs default phenology scheme based on rigorous intercomparison of 13 published phenological models (Yue et al., 2015a).This study represents the third step to simulate the closed climate-sensitive land carbon cycle: implementation of interactive carbon assimilation, allocation, autotrophic and heterotrophic respiration, and dynamic tree growth (changes in both height and LAI).For this third step, we purposefully select the mature, well-supported, well-established, readily available and accessible community algorithms: TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics; Cox, 2001;Clark et al., 2011) and the Carnegie-Ames-Stanford Approach (CASA) (Potter et al., 1993;Schaefer et al., 2008).TRIFFID has demonstrated previous usability in carbon-chemistry-climate interactions research., red)."mixed forests" are classified as ENF, "permanent wetlands", "savannas", and "woody savannas" as SHR.
The PFT of each site is described in Table S1 in the Supplement.
2 Observational data sets for validation

Site-level measurements
To validate the YIBs model, we use eddy covariance measurements from 145 flux tower sites (Fig. 1), which are collected by the North American Carbon Program (Schaefer et al., 2012) (K. Schaefer, personal communication, 2013) and downloaded from the FLUXNET (http://fluxnet.ornl.gov)network.Among these sites, 138 are located in the Northern Hemisphere, with 74 in Europe, 38 in the USA, and 24 in Canada (Table S1 in the Supplement).Sites on other continents are limited.Most of the sites have one dominant plant functional type (PFT), including 54 sites of evergreen needleleaf forests (ENF), 20 deciduous broadleaf forests (DBF), 9 evergreen broadleaf forests (EBF), 28 grasslands, 18 shrublands, and 16 croplands.We attribute sites with mixed forest to the ENF as these sites are usually at high latitudes.Each site's data set provides hourly or half-hourly measurements of carbon fluxes, including gross primary productivity (GPP) and net ecosystem exchange (NEE), and CO 2 concentrations and meteorological variables, such as surface air temperature, relative humidity, wind speed, and shortwave radiation.

Global measurements
We use global tree height, leaf area index (LAI), GPP, net primary productivity (NPP), and phenology data sets to validate the vegetation model.Canopy height is retrieved using 2005 remote sensing data from the Geoscience Laser Altimeter System (GLAS) aboard the ICESat satellite (Simard et al., 2011).LAI measurements for 1982-2011 are derived using the normalized difference vegetation index (NDVI) from Global Inventory Modeling and Mapping Studies (GIMMS; Zhu et al., 2013).Global GPP observations of 1982-2011 are estimated based on the upscaling of FLUXNET eddy covariance data with a biosphere model (Jung et al., 2009).This product was made to reproduce a model (LPJmL, Lund-Potsdam-Jena managed Land) using the fraction of absorbed PAR simulated in LPJmL.As a comparison, we also use GPP observations of 1982-2008 derived based on FLUXNET, satellite, and meteorological observations (Jung et al., 2011), which are about 10 % lower than those of Jung et al. (2009).
The NPP for 2000-2011 is derived using remote sensing data from Moderate Resolution Imaging Spectroradiometer (MODIS) (Zhao et al., 2005).We use the global retrieval of greenness onset derived from the Advanced Very High Resolution Radiometer (AVHRR) and the MODIS data from 1982 to 2011 (Zhang et al., 2014).All data sets are interpolated to the 1 • × 1 • offline model resolution for comparisons.

Vegetation biophysics
YIBs calculates carbon uptake for 9 PFTs: tundra, C3/C4 grass, shrubland, DBF, ENF, EBF, and C3/C4 cropland (Table 1).In the gridded large-scale model applications, each model PFT fraction in the vegetated part of each grid cell represents a single canopy.The vegetation biophysics simulates C3 and C4 photosynthesis with the well-established Michaelis-Menten enzyme-kinetics scheme (Farquhar et al., 1980;von Caemmerer and Farquhar, 1981) and the stomatal conductance model of Ball and Berry (Ball et al., 1987).The total leaf photosynthesis (A tot , µmol m −2 [leaf] s −1 ) is limited by one of three processes: (i) the capacity of the ribulose 1,5-bisphosphate (RuBP) carboxylase/oxygenase enzyme (rubisco) to catalyze carbon fixation (J c ), (ii) the capacity of the Calvin cycle and the thylakoid reactions to regenerate RuBP supported by electron transport (J e ), and (iii) the capacity of starch and sucrose synthesis to regenerate inorganic phosphate for photo-phosphorylation in C3 plants and phosphoenolpyruvate (PEP) in C4 plants (J s ).
The J c , J e , and J s are parameterized as functions of environmental variables (e.g., temperature, radiation, and CO 2 concentrations) and the maximum carboxylation capacity (V cmax , µmol m −2 s −1 ) (Collatz et al., 1991(Collatz et al., , 1992)): where c i and O i are the leaf internal partial pressure (Pa) of CO 2 and oxygen, * (Pa) is the CO 2 compensation point, and K c and K o (Pa) are Michaelis-Menten parameters for the carboxylation and oxygenation of rubisco.The parameters K c , K o , and * vary with temperature according to a Q 10 function.PAR (µmol m −2 s −1 ) is the incident photosynthetically active radiation, a leaf is leaf-specific light absorbance, and α is intrinsic quantum efficiency.P s is the ambient pressure and K s is a constant set to 4000 following Oleson et al. (2010).V cmax is a function of the optimal V cmax at 25 • C (V cmax25 ) based on a Q 10 function.Net carbon assimilation (A net ) of leaf is given by where R d is the rate of dark respiration set to 0.011 V cmax for C3 plants (Farquhar et al., 1980) and 0.025 V cmax for C4 plants (Clark et al., 2011).The stomatal conductance of water vapor (g s in mol [H 2 O] m −2 s −1 ) is dependent on net photosynthesis: where m and b are the slope and intercept derived from empirical fitting to the Ball and Berry stomatal conductance equations, RH is relative humidity, and c s is the CO 2 concentration at the leaf surface.In the model, the slope m is influenced by water stress, so that drought decreases photosynthesis by affecting stomatal conductance.Appropriate photosynthesis parameters for different PFTs are taken from Friend and Kiang (2005) and the Community Land Model (Oleson et al., 2010) with updates from Bonan et al. (2011) (Table 1).In future work, we will investigate the carbonchemistry-climate impacts of updated stomatal conductance models in YIBs (Berry et al., 2010;Pieruschka et al., 2010;Medlyn et al., 2011).The coupled equation system of photosynthesis, stomatal conductance and CO 2 diffusive flux transport equations form a cube in A net that is solved analytically (Baldocchi, 1994).A simplified but realistic representation of soil water stress β is included in the vegetation biophysics following the approach of Porporato et al. (2001).The algorithm reflects the relationship between soil water amount and the extent of stomatal closure ranging from no water stress to the soil moisture stress onset point (s * ) through to the wilting point (s wilt ).Stomatal conductance is reduced linearly between the PFTspecific values of s * and s wilt based on the climate model's soil water volumetric saturation in six soil layers (Unger et al., 2013).
The canopy radiative transfer scheme divides the canopy into an adaptive number of layers (typically 2-16) for light stratification.Each canopy layer distinguishes sunlit and shaded portions of leaves, so that the direct and diffuse photosynthetically active radiation (PAR) is used for carbon as-similation respectively (Spitters et al., 1986).The leaf photosynthesis is then integrated over all canopy layers to generate the GPP: GPP = LAI 0 A tot dL. (7)

Leaf phenology
Phenology determines the annual cycle of LAI.Plant phenology is generally controlled by temperature, water availability, and photoperiod (Richardson et al., 2013).For deciduous trees, the timing of budburst is sensitive to temperature (Vitasse et al., 2009) and the autumn senescence is related to both temperature and photoperiod (Delpierre et al., 2009).
For small trees and grasses, such as tundra, savanna, and shrubland, phenology is controlled by temperature and/or soil moisture, depending on the species type and locations of the vegetation (Delbart and Picard, 2007;Liu et al., 2013).In the YIBs model, leaf phenology is updated on a daily basis.
For the YIBs model, we build on the phenology scheme of Kim et al. (2015) and extend it based on long-term measurements of leaf phenology at five US sites (Yue et al., 2015a) and GPP at the 145 flux tower sties.A summary of the phenological parameters adopted is listed in Table 2.

Deciduous broadleaf forest (DBF)
We predict spring phenology of DBF using the cumulative thermal summation (White et al., 1997).The accumulative growing degree day (GDD) is calculated for the nth day from winter solstice if the 10-day average air temperature T 10 is higher than a base temperature T b : Here T b is set to 5 • C as in Murray et al. (1989).Similar to the approach outlined in Kim et al. (2015), the onset of greenness is triggered if the GDD exceeds a threshold value G b and a temperature-dependent phenological factor f T is calculated as follows: Following Murray et al. (1989), the threshold G b = a + b exp(r × NCD) is dependent on the number of chill days (NCD), which is calculated as the total days with < 5 • C from winter solstice.The autumn phenology is more uncertain than budburst because it is affected by both temperature and photoperiod (White et al., 1997;Delpierre et al., 2009).For the temperature-dependent phenology, we adopted the cumulative cold summation method (Dufrene et al., 2005;Richardson et al., 2006), which calculates the accumulative falling degree day (FDD) for the mth day from summer solstice as follows: where T s is 20 • C as in Dufrene et al. (2005).Similar to the budburst process, we determine the autumn phenological factor based on a fixed threshold F s : In addition, we assume photoperiod regulates leaf senescence as follows: where f P is the photoperiod-limited phenology.P is day length in minutes.P i and P x are the lower and upper limits of day length for the period of leaf fall.Finally, the autumn phenology of DBF is determined as the product of f T (Eq.11) and f P (Eq.12).Both the spring and autumn phenology schemes have been evaluated with extensive ground records over the USA in Y2015.

Shrubland
Shrub phenology is sensitive to temperature and/or water availability.We calculate correlation coefficients between Upper threshold of drought limit for herbs Dimensionless 0.9 Calibrated (Fig. S1) observed GPP and soil meteorology at 18 shrub sites (Fig. 2).
For 10 sites with annual mean soil temperature < 9 • C, the GPP-temperature correlations are close to 1 while the GPPmoisture correlations are all negative (Fig. 2a), suggesting that temperature is the dominant phenological driver for these plants.In contrast, for eight sites with an average soil temperature > 14 • C, GPP-moisture correlations are positive and usually higher than the GPP-temperature correlations, indicating that phenology is primarily regulated by water availability at climatologically warm areas.The wide temperature gap (9-14 • C) is due to the limit in the availability of shrub sites.Here, we select a tentative threshold of 12 • C to distinguish cold and drought species.We also try to identify phenological drivers based on soil moisture thresholds but find that both temperature-and drought-dependent phenology may occur at moderately dry conditions (Fig. 2b).
In the model, we apply the temperature-dependent phenology f T for shrubland if the site has annual mean soil temperature < 12 • C. We use the same f T as that for DBF (Eqs.9, 11), due to the lack of long-term phenology measurements at the shrub sites.However, if the soil temperature is > 12 • C, the plant growth is controlled by drought-limit phenology f D instead: where β 10 is the 10-day average water stress calculated based on soil moisture, soil ice fraction, and root fraction of each soil layer (Porporato et al., 2001).The value of β 10 changes from 0 to 1, with a lower value indicating drier soil.Two thresholds, β max and β min , represent the upper and lower thresholds that trigger the drought limit for woody species.The values of these thresholds are set to β max = 1 and β min = 0.4 so that the predicted phenology has the maximum correlations with the observed GPP seasonality (Fig. S1a in the Supplement).The shrub phenology applies for shrubland in tropical and subtropical areas, as well as tundra at the subarctic regions, though the phenology of the latter is usually dependent on temperature alone because the climatological soil temperature is < 12 • C.

Grassland
In the model, we consider temperature-dependent phenology for grassland based on soil temperature (ST) accumulation (White et al., 1997): where ST 10 is the 10-day average soil temperature and ST b = 0 • C. Similar to DBF, the onset of grass greenness is triggered if soil-temperature-based growing degree day (SGDD) is higher than a threshold value SG b (Kim et al., 2015): where SL g determines the grow length of grass.Both SG b and SL g are calibrated based on observed GPP seasonality at FLUXNET sites (Table 2).Grass phenology at warm sites is also sensitive to water stress (Fig. 2c).We apply the same drought-limit phenology f D as shrubland (Eq.13) for grassland but with calibrated thresholds β max = 0.9 and β min = 0.3 (Fig. S1b).Different from shrubland, whose phenology is dominated by drought when ST > 12 • C (Fig. 2a), grassland phenology is jointly affected by temperature and soil moisture (Fig. 2c).As a result, the final phenology for grassland at warm regions is the minimum of f T and f D .

Other PFTs
YIBs considers two evergreen PFTs, ENF at high latitudes and EBF in tropical areas.Observations do suggest that evergreen trees experience seasonal changes in LAI, following temperature variations and/or water availability (Doughty and Goulden, 2008;Schuster et al., 2014).However, due to the large uncertainty of evergreen phenology, we set a constant phenology factor of 1.0 for these species, following the approach adopted in other process-based vegetation models (Bonan et al., 2003;Sitch et al., 2003).We implement a parameterization for the impact of cold temperature (frost hardening) on the maximum carboxylation capacity (V cmax ) so as to reduce cold injury for ENF during winter (Hanninen and Kramer, 2007).EBF may experience reduced photosynthesis during the dry season through the effects of water stress on stomatal conductance (Jones et al., 2014).
Crop phenology depends on planting and harvesting dates.In YIBs, we apply a global data set of crop planting and harvesting dates (Sacks et al., 2010;Unger et al., 2013).Crop budburst occurs at the planting date and the crop continues to grow for a period of 30 days until reaching full maturity (f = 1).The crop leaves begin to fall 15 days prior to the harvest date, after which phenology is set to 0. A similar treatment has been adopted in the CLM (Community Land Model; Bonan et al., 2003).Thus, crop productivity but not crop phenology is sensitive to the imposed meteorological forcings.

Carbon allocation
We adopt the autotrophic respiration and carbon allocation scheme applied in the dynamic global vegetation model (DGVM) TRIFFID (Cox, 2001;Clark et al., 2011).On a daily basis, the plant LAI is updated as follows: where f is the phenological factor, and LAI b is the biomassbalanced (or available maximum) LAI related to tree height.LAI b is dependent on the vegetation carbon content C veg , which is the sum of carbon from leaf (C l ), root (C r ), and stem (C w ): where each carbon component is a function of LAI b : Here σ l is the specific leaf carbon density.a wl and b wl are PFT-specified allometric parameters (Table 1).The vegetation carbon content C veg is updated every 10 days based on the carbon balance of assimilation, respiration, and litter fall.
The NPP is the net carbon uptake: Here GPP is the total photosynthesis rate integrated over LAI.Autotrophic respiration (R a ) is split into maintenance (R am ) and growth respiration (R ag ) (Clark et al., 2011): The maintenance respiration is calculated based on nitrogen content in leaf (N l ), root (N r ), and stem (N w ) as follows: where R d is the dark respiration of leaf, which is dependent on leaf temperature and is integrated over whole canopy LAI.The factor of 0.012 is the unit conversion (from mol CO 2 m −2 s −1 to kg C m −2 s −1 ) and β is water stress representing soil water availability.The nitrogen contents are given by Here n 0 is leaf nitrogen concentration, n rl and n wl are ratios of nitrogen concentrations of root and stem to leaves, η is a factor scaling live stem mass to LAI and tree height H .We adopt the same values of n 0 , n rl , n wl and η as that of the TRIFFID model (Table 1) except that n rl is set to 0.5 following observations of deciduous trees by Sugiura and Tateno (2011).The growth respiration is dependent on the residual between GPP and R am based on a ratio r g set to 0.2 for all PFTs (Knorr, 2000): The λ in Eq. ( 19) is a partitioning coefficient determining the fraction of NPP used for spreading: where LAI min and LAI max are minimum and maximum LAI values for a specific PFT (Table 1).In the current model version, we turn off the fractional changes by omitting λNPP in the carbon allocation but feeding it as input for the soil respiration.The litter fall rate l in Eq. ( 19) consists of contributions from leaf, root, and stem as follows: Here γ l , γ r , and γ w are turnover rate (year −1 ) for leaf, root, and stem carbon respectively.The leaf turnover rate is calculated based on the phenology change every day.The root and stem turnover rates are PFT-specific constants (Table 1), derived based on the meta-analysis by Gill and Jackson (2000) for root and Stephenson and van Mantgem (2005) for stem.

Soil respiration
The soil respiration scheme is developed based on the CASA model (Potter et al., 1993;Schaefer et al., 2008), which considers carbon flows among 12 biogeochemical pools.Three live pools, including leaf C l , root C r , and wood C w , contain biomass carbon assimilated from photosynthesis.Litterfall from live pools decomposes and transits in nine dead pools, which consist of one coarse woody debris (CWD) pool, three surface pools, and five soil pools.The CWD pool is composed of dead trees and woody roots.Both surface and soil have identical pools, namely structural, metabolic, and microbial pools, which are distinguished by the content and functions.The structural pool contains lignin, the metabolic pool contains labile substrates, and the microbial pool represents microbial populations.The remaining two soil pools, the slow and passive pools, consist of organic material that decays slowly.The full list of carbon flows among different pools has been illustrated by Schaefer et al. (2008) (cf.their Fig. 1).When carbon transfers from pool j to pool i, the carbon loss of pool j is where C j is the carbon in pool j , k j is the total carbon loss rate of pool j , and f j 2i is the fraction of carbon lost from pool j transferred to pool i.The coefficient k j is dependent on soil temperature, moisture, and texture.Meanwhile, the carbon gain of pool i is, where e j 2i is the ratio of carbon received by pool i to the total carbon transferred from pool j .The rest of the transferred carbon is lost due to heterotrophic respiration: As a result, the carbon in the ith pool is calculated as The total heterotrophic respiration (R h ) is the summation of R j 2i for all pair pools where carbon transitions occur.The total soil carbon is the summation of carbon for all dead pools: The net ecosystem productivity (NEP) is calculated as where NEE is the net ecosystem exchange, representing net carbon flow from land to atmosphere.YIBs does not yet account for NEE perturbations due to dynamic disturbance.

Ozone vegetation damage effects
We apply the semi-mechanistic parameterization proposed by Sitch et al. (2007) to account for ozone damage to photosynthesis through stomatal uptake.The scheme simulates associated changes in both photosynthetic rate and stomatal conductance.When photosynthesis is inhibited by ozone, stomatal conductance decreases accordingly to resist more ozone molecules.We employed an offline regional version of YIBs to show that present-day ozone damage decreases GPP by 4-8 % on average in the eastern USA and leads to larger decreases of 11-17 % in east coast hotspots (Yue and Unger, 2014).In the current model version, the photosynthesis and stomatal conductance responses to ozone damage are coupled.In future work, we will update the ozone vegetation damage function in YIBs to account for decoupled photosynthesis and stomatal conductance responses based on recent extensive metadata analyses (Wittig et al., 2007;Lombardozzi et al., 2013).

BVOC emissions
YIBs incorporates two independent leaf-level isoprene emission schemes embedded within the exact same host model framework (Zheng et al., 2015).The photosynthesis-based isoprene scheme simulates emission as a function of the electron transport-limited photosynthesis rate (J e , Eq. 3), canopy temperature, intercellular CO 2 (c i ) and * (Arneth et al., 2007;Unger et al., 2013).The MEGAN (Model of Emissions of Gases and Aerosols from Nature) scheme applies the commonly used leaf-level functions of light and canopy temperature (Guenther et al., 1993(Guenther et al., , 1995(Guenther et al., , 2012)).Both isoprene schemes account for atmospheric CO 2 sensitivity (Arneth et al., 2007).Long-term increases (decreases) in atmospheric CO 2 decrease (increase) isoprene emissions (Unger et al., 2013).The CO 2 sensitivity is higher under lower atmospheric CO 2 levels than at present day.Leaf-level monoterpene emissions are simulated using a simplified temperaturedependent algorithm (Lathière et al., 2006).The leaf-level isoprene and monoterpene emissions are integrated over the multiple canopy layers in the exact same way as GPP to obtain the total canopy-level emissions.

Implementation of YIBs into NASA ModelE2 (NASA ModelE2-YIBs)
NASA ModelE2 has a spatial resolution of 2 • × 2.5 • latitude by longitude with 40 vertical levels extending to 0.1 hPa.In the online configuration, the global climate model provides the meteorological drivers to YIBs and the land-surface hydrology submodel provides the soil characteristics (Rosenzweig and Abramopoulos, 1997; Schmidt et al., 2014).Recent relevant updates to NASA ModelE2 include a dynamic fire activity parameterization from Pechony and Shindell ( 2009) and climate-sensitive soil NO x emissions based on Yienger and Levy (1995) (Unger and Yue, 2014).Without the YIBs implementation, the default NASA ModelE2 computes dry deposition using fixed LAI and vegetation cover fields from Olson et al. (2001), which are different from the climate model's vegetation scheme (Shindell et al., 2013b).With YIBs embedded in NASA ModelE2, the YIBs model provides the vegetation cover and LAI for the dry deposition scheme.The online-simulated atmospheric ozone and aerosol concentrations influence terrestrial carbon assimilation and stomatal conductance at the 30 min integration time step.In turn, the online vegetation properties, and water, energy and BVOC fluxes affect air quality, meteorology and the atmospheric chemical composition.The model simulates the interactive deposition of inorganic and organic nitrogen to the terrestrial biosphere.However, the YIBs biosphere currently applies fixed nitrogen levels and does not yet account for the dynamic interactions between the carbon and nitrogen cycles, and the consequences for carbon assimilation, which are highly uncertain (e.g., Thornton et al., 2007;Koven et al., 2013;Thomas et al., 2013;Zaehle et al., 2014;Houlton et al., 2015).

Site-level simulations (YIBs-site)
We perform site-level simulations with the offline YIBs model at 145 eddy covariance flux tower sites for the corresponding PFTs (Fig. 1).Hourly in situ measurements of meteorology (Sect.2.1) are used as input for the model.We gap-filled missing measurements with the Global Modeling and Assimilation Office (GMAO) Modern Era-Retrospective Analysis (MERRA) reanalysis (Rienecker et al., 2011), as described in Yue and Unger (2014).All grasslands and most croplands are considered as C3 plants, except for some sites where corn is grown.Meteorological measurements are available for a wide range of time periods across the different sites ranging from the minimum of 1 year at some sites (e.g., BE-Jal) to the maximum of 16 years at Harvard Forest (US-HA1).The soil carbon pool initial conditions at each site are provided by the 140-year spinup procedure using YIBsoffline (Supplement).An additional 30-year spinup is conducted for each site-level simulation using the initial height H 0 for corresponding PFT (Table 1) and the fixed meteorology and CO 2 conditions at the first year of observations.Then, the simulation is continued with year-to-year forcings at the specific site for the rest of the measurement period.For all grass and shrub sites, two simulations are performed.One applies additional drought controls on phenology as described in Sects.3.2.2 and 3.2.3,while the other uses only temperature-dependent phenology.By comparing results of these two simulations, we assess the role of drought phenology for plants in arid and semi-arid regions.

Global offline simulation (YIBs-offline)
The global offline YIBs applies the CLM land cover data set (Oleson et al., 2010).Land cover is derived based on retrievals from both MODIS (Hansen et al., 2003) and AVHRR (Defries et al., 2000).Fractions of 16 PFTs are aggregated into 9 model PFTs (Table 1).The soil carbon pool and tree height initial conditions are provided by the 140-year spinup procedure using YIBs-offline (Supplement).The global offline YIBs model is driven with WFDEI meteorology (Weedon et al., 2014) at 1 • × 1 • horizontal resolution for the period of 1980-2011.Observed atmospheric CO 2 concentrations are adopted from the fifth assessment report (AR5) of the Intergovernmental Panel on Climate Change (IPCC; Meinshausen et al., 2011).We evaluate the simulated longterm 1980-2011 average tree height/LAI and carbon fluxes with available observations and recent multimodel intercomparisons.Attribution of the decadal trends in terrestrial carbon fluxes are explored in a separate follow-on companion study (Yue et al., 2015b).

Global online simulation in NASA ModelE2-YIBs
The global land cover data is identical to that used in YIBsoffline (Sect.4.2) based on the CLM cover.Because our major research goal is to study short-term (seasonal, annual, decadal) interactions between vegetation physiology and atmospheric chemistry, we elect to prescribe the PFT distribution in different climatic states.We perform an online atmosphere-only simulation representative of the present day (∼ 2000s) climatology by prescribing fixed monthlyaverage sea surface temperature (SST) and sea ice temperature for the 1996-2005 decade from the Hadley Center as the boundary conditions (Rayner et al., 2006) soil carbon pools and tree heights are provided by the 140year spinup process described in the Supplement using YIBsoffline but for year 2000 (not 1980) fixed WFDEI meteorology and atmospheric CO 2 conditions.The NASA ModelE2-YIBs global carbon-chemistry-climate model is run for an additional 30 model years.The first 20 years are discarded as the online spinup and the last 10-year results are averaged for the analyses including comparisons with observations and the YIBs-offline.

Ozone vegetation damage simulation (YIBs-ozone)
We perform two simulations to quantify ozone vegetation damage with the offline YIBs model based on the high and low ozone sensitivity parameterizations (Sitch et al., 2007).Similar to the setup in Yue and Unger ( 2014), we use offline hourly surface ozone concentrations simulated with the NASA ModelE2 based on the climatology and precursor emissions of the year 2000 (Sect.4.3).In this way, atmospheric ozone photosynthesis damage affects plant growth, including changes in tree height and LAI.We compare the simulated ozone damage effects with the previous results in Yue and Unger ( 2014) that used prescribed LAI.For this updated assessment, we do not isolate possible feedbacks from the resultant land carbon cycle changes to the surface ozone concentrations themselves, for instance through concomitant changes to BVOC emissions and water fluxes.The importance of these feedbacks will be quantified in future research using the online NASA ModelE2-YIBs framework.

Site-level evaluation
The simulated monthly-average GPP is compared with measurements at 145 sites for different PFTs (Fig. 3).GPP simulation biases range from −19 to 7 % depending on PFT.The highest correlation of 0.86 is achieved for DBF, mainly contributed by the reasonable phenology simulated at these sites (Fig. S2).The correlation is also high for ENF sites even though phenology is set to a constant value of 1.0.A relatively low correlation of 0.65 is modeled for EBF sites (Fig. S2).However, the site-specific evaluation shows that the simulations reasonably capture the observed magnitude and seasonality, including the minimum GPP in summer due to drought at some sites (e.g., FR-Pue and IT-Lec).Predictions at crop sites achieve a medium correlation of 0.77, because the prescribed crop phenology based on the planting and harvesting dates data set matches reality for most sites with some exceptions (e.g., CH-Oe2).Measured GPP at shrub and grass sites show varied seasonality.For most sites, the maximum carbon fluxes are measured in the hemispheric summer season.However, for sites with arid or Mediterranean climate, the summer GPP is usually the lowest during the year (e.g., ES-LMa and US-Var in Fig. S2) while the peak flux is observed during the wet season when the climate is cooler and moister.Implementing the drought-dependent phenology helps improve the GPP seasonality and decrease the root-mean-square error (RMSE) at most warm climate shrub and grass sites (Fig. S3).
A synthesis of the site-level evaluation is presented in Fig. 4.Among the 145 sites, 121 have correlations higher than 0.8 for the GPP simulation (Fig. 4a).Predictions are better for PFTs with larger seasonal variations.For example, high correlations of > 0.8 are achieved at 95 % for the ENF and DBF sites, but only 70 % for grass and 45 % for EBF sites.Low relative biases (−33 to 50 %) are achieved at 94 sites (Fig. 4b).For most PFTs, a similar fraction (65 %) of the sites have low biases falling into that range, except for cropland, where only seven sites (45 %) have the low biases.The RMSE is lower than 3 g [C] day −1 for 107 out of 145 sites (Fig. 4c).The highest RMSE is predicted for crop sites, where the model misses the large interannual variations due to crop rotation at some sites (e.g., BE-Lon, DE-Geb, and US-Ne2).The YIBs model performs simulations at the PFT level while measurements show large uncertainties in the carbon fluxes among biomes/species within the same PFT (Luyssaert et al., 2007).The simulated intraspecific variations (in the form of standard deviation) are smaller than the measured/derived values for most PFTs (Table S2), likely because of the application of fixed photosynthetic parameters for each PFT (Table 1).
Compared with GPP, the NEE simulations have smaller correlations with measurements because of the limited sea- sonality in the observations at most sites (Fig. S4).In total, 74 sites (51 %) have correlation coefficients higher than 0.6 (Fig. 4d) and 75 sites (52 %) have absolute biases within ± 0.5 g [C] day −1 (Fig. 4e).For most ENF sites, the maximum net carbon uptake (the minimum NEE) is observed in spring or early summer, when GPP begins to increase while soil respiration is still at a low rate due to the cool and wet conditions (e.g., CA-Ojp and ES-ES1).Compared with other PFTs, the DBF trees usually have larger seasonality with the NEE peak in the early summer.Such seasonality helps promote correlations between model and measurements, resulting in high R (> 0.8) for 17 out of 20 sites (Fig. 4d).For shrub and grass sites, the observed seasonality of NEE is not regular, though most show maximum carbon uptake in spring or early summer.Implementation of drought-dependent phenology helps improve the simulated NEE seasonality at some sites of these PFTs (e.g., ES-LMa and IT-Pia); however, such improvement is limited for others (Fig. S4).Simulated crop NEE reaches its maximum magnitude in summer at most sites, consistent with observations, leading to a high R (> 0.8) in 10 out 16 sites (Fig. 4d).The RMSE of simulated NEE is larger for crops relative to other PFTs because the model does not treat crop rotation (Fig. 4f).

Evaluation of YIBs-offline
YIBs-offline forced with WFDEI meteorology simulates reasonable spatial distributions for tree height, LAI, and GPP, all of which show maximums in the tropical rainforest biome and medium values in the Northern Hemisphere high latitudes (Fig. 5).Compared with the satellite observations, the simulated height is underestimated by 30 % on the annual and global mean basis (Fig. 5b).Regionally, the prediction is larger by only 4 % for tropical rainforest and temperate DBF, but by 27 % for boreal ENF, for which the model assumes a constant phenology of 1.0 all the year round.However, for the vast areas covered with grass and shrub PFTs, the simulated height is lower by 41 % with maximum underestimation in eastern Siberia, where the model land is covered by short tundra.The modeled LAI is remarkably close to observations on the annual and global mean basis (Fig. 5c, d).
However, there are substantial regional biases in model LAI.
Model LAI prediction is higher by 0.8 m 2 m −2 (70 %) for boreal ENF and by 0.1 m 2 m −2 (5 %) for tropical rainforest.In contrast, the simulation underestimates LAI of tropical C4 grass by 0.4 m 2 m −2 (30 %) and shrubland by 0.2 m 2 m −2 (30 %).The GPP simulation is lower than the FLUXNETderived value by 5 % on the global scale, which is contributed by the minor underestimation for all PFTs except for tropical rainforest, where model predicts 9 % higher GPP than observations (Fig. 5f).
The model simulates reasonable seasonality for LAI and land carbon fluxes (Fig. 6).Tree height shows limited sea- sonal variations, especially for DBF, ENF, and EBF trees.LAI, GPP, and NPP also exhibit small seasonality over tropical areas, such as the Amazon, central Africa, and Indonesia.However, for temperate areas, such as North America, Europe and east Asia, these variables show large seasonal variations with a minimum in winter and maximum in summer.The LAI is overestimated by 20 % in the Amazon during the December-January-February season but underestimated by 25 % in Indonesia during summer (Fig. 6b).For GPP and NPP, the positive bias in Indonesia is even larger at 45 % during summer (Fig. 6c, d).
On the global scale, YIBs-offline simulates a GPP of 124.6 ± 3.3 Pg C a −1 and NEE of −2.5 ± 0.7 Pg C a −1 for 1982-2011.These values are consistent with estimates upscaled from the FLUXNET observations (Jung et al., 2009(Jung et al., , 2011;;Friedlingstein et al., 2010) and simulations from 10 other carbon cycle models (Piao et al., 2013) (Fig. 7).The net biome productivity (NBP) is in opposite sign to NEE.Tropical areas (23 • S-23 • N) account for 63 % of the global GPP, including 27 % from the Amazon rainforest, 21 % from central Africa, and 5 % from Indonesia forest (Table 3).A lower contribution of 57 % from the tropics is predicted for both NPP and heterotrophic respiration.However, for NEE, only 40 % of the land carbon sink is contributed by tropical forests and grasslands, while 56 % is from temperate forests and grasslands in North America, Europe, and east Asia.
We compare the simulated budburst dates with observations from satellite retrieval (Fig. 8).The model captures the basic spatial pattern of spring phenology with earlier to later budburst dates from lower to higher latitudes.On average, the observed budburst date in the Northern Hemisphere is DOY 133 (13 May) and the simulated is DOY 132 (12 May).Such close estimate results from the regional delay of 10 days (DOY 119 vs. 129) in Europe and advance of 4 days (DOY 140 vs. 136) in east Asia.In Y2015, extensive (∼ 75 000 records) ground-based measurements have been used to validate the simulated spring and autumn phenology in the USA and both the spatial distribution and interannual variation of the simulation are reasonable.

Evaluation of NASA ModelE2-YIBs
NASA ModelE2-YIBs simulations of global land carbon fluxes show similar spatial distribution and magnitude as the YIBs-offline model (Figs.S6-S8).However, due to differences in the meteorological forcings (Figs.S9-S12), regional discrepancies between the two configurations occur.The predicted LAI with NASA ModelE2-YIBs is lower by 20 % in the Amazon region than with YIBs-offline (Fig. S6), following the similar magnitude of differences in regional GPP and NPP (Figs.S7, S8).We performed driver attribution sensitivity simulations in which the YIBs-offline con-   figuration is driven with the same meteorological forcings simulated by NASA ModelE2 except for one selected field from the WFDEI reanalysis.We found that the anomalously warmer climate over the Amazon in the global climate model (Fig. S9) causes the lower GPP in that region in NASA ModelE2-YIBs.The temperature optimum for C3 photosynthesis is around 30 • C, above which the maximum rate of electron transport (Eq. 3) decreases dramatically (Farquhar et al., 1980).As a result, the higher NASA ModelE2-YIBs surface temperature in the tropical rainforest results in the lower photosynthesis rates there.With the exception of the Amazon, the NASA ModelE2-YIBs June-July-August GPP and NPP show low biases in central Africa and high latitudes in North America and Asia but high biases in Europe, western USA, and eastern China (Figs.S7, S8).The sensitivity tests attribute these discrepancies to differences in canopy humidity (Fig. S11) and soil wetness (Fig. S12).Low soil wetness decreases water stress β, reduces the slope m of the Ball-Berry equation (Eq.6), and consequently limits photosynthesis by declining stomatal conductance in combination with low humidity.On the global scale, the ModelE2-YIBs simulates annual GPP of 122.9 Pg C, NPP of 62 Pg C, and NEE of −2.7 Pg C, all of which are close to the YIBs-offline simulation (Table 3) and consistent with results from observations and model intercomparison (Fig. 7).

Assessment of global ozone vegetation damage
Ozone dampens GPP and consequently affects tree growth and LAI.In North America, the annual average reductions range from 2 to 6 %, depending on the plant sensitivity to ozone damage (Table 3).Locally, average damages reach as high as 5-11 % in the eastern USA with maximums up to 11-23 % (Fig. 9a, b).These values are higher than the estimate of 4-8 % (maximum 11-17 %) by Yue and Unger (2014) because the latter used prescribed LAI in the simulation and did not consider the LAI reductions due to ozone damage (Fig. 9c, d).The YIBs model predicts a similar magnitude of damages in Europe compared to North America but almost doubled its effects in east Asia (Table 3) due to the high ozone concentrations there, especially in boreal summer (Fig. S5).Ozone-induced GPP reductions are limited in tropical areas (Fig. 5e) because the surface ozone levels there are very low, for example, especially over the Amazon forest (Fig. S5).The damage to LAI generally follows the pattern of GPP reductions but the response signal is weaker than that of GPP (Fig. 9c, d).

Conclusions and discussion
We describe and evaluate the process-based YIBs interactive terrestrial biosphere model.YIBs is embedded into the NASA ModelE2 global chemistry-climate model and is an important urgently needed development to improve the biological realism of interactions between vegetation, atmospheric chemistry and climate.We implement both temperature-and drought-dependent phenology for DBF, shrub, and grass species.The model simulates interactive ozone vegetation damage.The YIBs model is fully validated with land carbon flux measurements from 145 ground stations and global observations of canopy height, LAI, GPP, NPP, and phenology from multiple satellite retrievals.
There are several limitations in the current model setup.The vegetation parameters, V cmax25 , m, and b (Table 1) are fixed at the PFT level, which may induce uncertainties in the simulation of carbon fluxes due to intraspecific variations (Kattge et al., 2011).The model does not yet include a dynamic treatment of nitrogen and phosphorous availability because current schemes suffer from large uncertainties (Thornton et al., 2007;Zaehle et al., 2014;Houlton et al., 2015).Phenology is set to a constant value of 1 for ENF and EBF, which is not consistent with observations (O' Keefe, 2000;Jones et al., 2014).The ozone damage scheme of Sitch et al. (2007) considers coupled responses of photosynthesis and stomatal conductance while observations suggest a decoupling (Lombardozzi et al., 2013).
Despite these limitations, the YIBs model reasonably simulates global land carbon fluxes compared with both sitelevel flux measurements and global satellite observations.YIBs is primed for ongoing development, for example, incorporating community dynamics including mortality, establishment, seed transport and dynamic fire disturbance (Moorcroft et al., 2001).NASA ModelE2-YIBs is available to be integrated with interactive ocean and atmospheric carbon components to offer a full global carbon-climate model, for example, for use in interpreting and diagnosing new satellite data sets of atmospheric CO 2 concentrations.In the current form, NASA ModelE2-YIBs provides a useful new tool to investigate the impacts of air pollution on the carbon budget, water cycle, and surface energy balance, and, in turn, the impacts of changing vegetation physiology on the atmospheric chemical composition.Carbon-chemistry-climate interactions, a relatively new interdisciplinary research frontier, are expected to influence the evolution of Earth's climate system on multiple spatiotemporal scales.

Figure 2 .
Figure 2. Correlations between monthly GPP and soil variables at (a, b) shrub and(c, d) grass sites.For each site, we calculate correlation coefficients of GPP-soil temperature (red points) and GPPsoil moisture (blue squares).These correlation coefficients are then plotted against the annual mean (a, c) soil temperature ( • C) or (b, d) soil moisture (fraction) at each site.

Figure 3 .
Figure 3.Comparison between observed and simulated monthly GPP (in g C m −2 day −1 ) from FLUXNET and NACP networks grouped by PFTs.Each point represents the average value of 1 month at one site.The red lines indicate linear regression between observations and simulations.The regression fit, correlation coefficient, and relative bias are shown on each panel.The PFTs include evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), shrubland (SHR), grassland (GRA), and cropland (CRO).The detailed comparison for each site is shown in Fig. S2.

Figure 4 .
Figure 4. Bar charts of (a, d) correlation coefficients (R), (b, e) biases, and (c, f) RMSE for monthly (a-c) GPP and (d-f) NEE between simulations and observations at 145 sites.Each bar represents the number of sites where the R, bias, or RMSE of simulations fall between the specific ranges as defined by the x axis intervals.The minimum and maximum of each statistical metric are indicated as the two ends of the x axis in the plots.The values of the x axis are not even.The absolute biases instead of relative biases are shown for NEE because the long-term average NEE (the denominator) is usually close to zero at most sites.The PFT definitions are as in Fig. 2. Detailed comparisons at each site are shown in Figs.S2 and S4.

Figure 5 .
Figure 5. Simulated (a) tree height, (c) LAI, and (e) GPP and their differences relative to observations (b, d, f).GPP data set is from Jung et al. (2009).Simulations are performed with WFDEI reanalysis.Statistics are the annual average for the period 1982-2011.The boxes in panel (a) represent six regions used for seasonal comparison in Fig. 6.

Figure 6 .
Figure 6.Comparison of annual (a) tree height and seasonal (b) LAI, (c) GPP, and (d) (NPP between simulations and observations for the six regions shown in Fig. 5a.GPP data set is from Jung et al. (2009).Values at different regions are marked using different symbols, with distinct colors indicating seasonal means for winter (blue, December-February), spring (green, March-May), summer (red, June-August), and autumn (magenta, September-November).

Figure 7 .
Figure 7.Comparison of simulated global GPP and net biome productivity (NBP) from (red) YIBs-offline and (blue) ModelE2-YIBs models with 10 other carbon cycle models for 1982-2008.Each black symbol represents an independent model as summarized in Piao et al. (2013).Error bars indicate the standard deviations for interannual variability.The gray shading represents the global residual land sink (RLS) calculated in Friedlingstein et al. (2010).The green line at the top represents the range of GPP for 1982-2008 estimated by Jung et al. (2011) and the magenta line represents GPP for 1982-2011 from Jung et al. (2009).

Figure 8 .
Figure 8.Comparison of simulated budburst dates in the Northern Hemisphere with remote sensing.Simulated phenology in each grid square is the composite result from DBF, tundra, shrubland, and grassland based on PFT fraction and LAI in that grid box.Both simulations and observations are averaged for the period 1982-2011.Results for the Southern Hemisphere are not shown due to the limited coverage of deciduous forests and cold grass species.

Figure 9 .
Figure 9. Percentage of ozone vegetation damage to (top) GPP and (bottom) LAI with (a, c) low and (b, d) high sensitivity.Both damages of GPP and LAI are averaged for 1982-2011.Offline surface ozone concentrations (Fig. S5) are simulated by the GISS ModelE2 with climatology of the year 2000.

Table 1 .
Photosynthetic and allometric parameters for the vegetation model.

Table 2 .
Phenological parameters for the vegetation model.

Table 3 .
Summary of carbon fluxes and ozone vegetation damage in different domains and for the tropics (23 • S-23 • N).