An observational constraint on stomatal function in forests: evaluating coupled carbon and water vapor exchange with carbon isotopes in the Community Land Model (CLM4.5)

Abstract. Land surface models are useful tools to quantify contemporary and future climate impact on terrestrial carbon cycle processes, provided they can be appropriately constrained and tested with observations. Stable carbon isotopes of CO2 offer the potential to improve model representation of the coupled carbon and water cycles because they are strongly influenced by stomatal function. Recently, a representation of stable carbon isotope discrimination was incorporated into the Community Land Model component of the Community Earth System Model. Here, we tested the model's capability to simulate whole-forest isotope discrimination in a subalpine conifer forest at Niwot Ridge, Colorado, USA. We distinguished between isotopic behavior in response to a decrease of δ13C within atmospheric CO2 (Suess effect) vs. photosynthetic discrimination (Δcanopy), by creating a site-customized atmospheric CO2 and δ13C of CO2 time series. We implemented a seasonally varying Vcmax model calibration that best matched site observations of net CO2 carbon exchange, latent heat exchange, and biomass. The model accurately simulated observed δ13C of needle and stem tissue, but underestimated the δ13C of bulk soil carbon by 1–2 ‰. The model overestimated the multiyear (2006–2012) average Δcanopy relative to prior data-based estimates by 2–4 ‰. The amplitude of the average seasonal cycle of Δcanopy (i.e., higher in spring/fall as compared to summer) was correctly modeled but only when using a revised, fully coupled An − gs (net assimilation rate, stomatal conductance) version of the model in contrast to the partially coupled An − gs version used in the default model. The model attributed most of the seasonal variation in discrimination to An, whereas interannual variation in simulated Δcanopy during the summer months was driven by stomatal response to vapor pressure deficit (VPD). The model simulated a 10 % increase in both photosynthetic discrimination and water-use efficiency (WUE) since 1850 which is counter to established relationships between discrimination and WUE. The isotope observations used here to constrain CLM suggest (1) the model overestimated stomatal conductance and (2) the default CLM approach to representing nitrogen limitation (partially coupled model) was not capable of reproducing observed trends in discrimination. These findings demonstrate that isotope observations can provide important information related to stomatal function driven by environmental stress from VPD and nitrogen limitation. Future versions of CLM that incorporate carbon isotope discrimination are likely to benefit from explicit inclusion of mesophyll conductance.

Abstract. Land surface models are useful tools to quantify contemporary and future climate impact on terrestrial carbon cycle processes, provided they can be appropriately constrained and tested with observations. Stable carbon isotopes of CO 2 offer the potential to improve model representation of the coupled carbon and water cycles because they are strongly influenced by stomatal function. Recently, a representation of stable carbon isotope discrimination was incorporated into the Community Land Model component of the Community Earth System Model. Here, we tested the model's capability to simulate whole-forest isotope discrimination in a subalpine conifer forest at Niwot Ridge, Colorado, USA. We distinguished between isotopic behavior in response to a decrease of δ 13 C within atmospheric CO 2 (Suess effect) vs. photosynthetic discrimination ( canopy ), by creating a site-customized atmospheric CO 2 and δ 13 C of CO 2 time series. We implemented a seasonally varying V cmax model calibration that best matched site observations of net CO 2 carbon exchange, latent heat exchange, and biomass. The model accurately simulated observed δ 13 C of needle and stem tissue, but underestimated the δ 13 C of bulk soil carbon by 1-2 ‰. The model overestimated the multiyear (2006)(2007)(2008)(2009)(2010)(2011)(2012) average canopy relative to prior data-based estimates by 2-4 ‰. The amplitude of the average seasonal cycle of canopy (i.e., higher in spring/fall as compared to summer) was correctly modeled but only when using a revised, fully coupled A n −g s (net assimilation rate, stomatal conduc-tance) version of the model in contrast to the partially coupled A n − g s version used in the default model. The model attributed most of the seasonal variation in discrimination to A n , whereas interannual variation in simulated canopy during the summer months was driven by stomatal response to vapor pressure deficit (VPD). The model simulated a 10 % increase in both photosynthetic discrimination and water-use efficiency (WUE) since 1850 which is counter to established relationships between discrimination and WUE. The isotope observations used here to constrain CLM suggest (1) the model overestimated stomatal conductance and (2) the default CLM approach to representing nitrogen limitation (partially coupled model) was not capable of reproducing observed trends in discrimination. These findings demonstrate that isotope observations can provide important information related to stomatal function driven by environmental stress from VPD and nitrogen limitation. Future versions of CLM that incorporate carbon isotope discrimination are likely to benefit from explicit inclusion of mesophyll conductance.
B. Raczka et al.: An observational constraint on stomatal function in forests face (Le Quéré et al., 2015), but it is unclear how projected changes in temperature and precipitation will influence the future of this land carbon sink (Arora et al., 2013;Friedlingstein et al., 2006). A major source of uncertainty in climate model projections results from the disagreement in projected strength of the land carbon sink (Arora et al., 2013). Thus, it is critical to reduce this uncertainty to improve climate predictions, and to better inform mitigation strategies (Yohe et al., 2007).
An effective approach to reduce uncertainties in terrestrial carbon models is to constrain a broad range of processes using distinct and complementary observations. Traditionally, terrestrial carbon models have relied primarily upon observations of land-surface fluxes of carbon, water, and energy derived from eddy covariance flux towers to calibrate model parameters and evaluate model skill. Flux measurements best constrain processes that occur at diurnal and seasonal timescales (Braswell et al., 2005;Ricciuto et al., 2008). Traditional ecological metrics of carbon pools (e.g., leaf area index (LAI), biomass) are also commonly used to provide independent and complementary constraints upon ecosystem processes at longer timescales (Ricciuto et al., 2011;Richardson et al., 2010). However, neither flux nor carbon pool observations provide suitable constraints for the model formulation of plant stomatal function and the related link between the carbon and water cycles.
Stable carbon isotopes of CO 2 are influenced by stomatal activity in C3 plants (e.g., evergreen trees, deciduous trees), and thus provide a valuable but under-utilized constraint on terrestrial carbon models. Plants assimilate more of the lighter of the two major isotopes of atmospheric carbon ( 12 C vs. 13 C). This preference, termed photosynthetic discrimination ( canopy ), is primarily a function of two processes, CO 2 diffusion rate through the leaf boundary layer and into the stomata, and the carboxylation of CO 2 . The magnitude of canopy is controlled by CO 2 supply (depending on, for instance, atmospheric CO 2 concentration and stomatal conductance) and demand (depending on, for instance, photosynthetic rate; Flanagan et al., 2012). In general, environmental conditions favorable to plant productivity result in higher canopy during carbon assimilation compared to unfavorable conditions. Plants respond to unfavorable conditions by closing stomata and reducing the stomatal conductance which reduces canopy . Most relevant here, canopy responds to atmospheric moisture deficit (Andrews et al., 2012;Wingate et al., 2010), soil water content (McDowell et al., 2010), precipitation (Roden and Ehleringer, 2007), and nutrient availability (Cernusak et al., 2013). After carbon is assimilated, additional post-photosynthetic isotopic changes occur Brüggemann et al., 2011), but these impose a small influence on land-atmosphere isotopic exchange relative to photosynthetic discrimination.
The Niwot Ridge Ameriflux site, located in a subalpine conifer forest in the Rocky Mountains of Colorado, USA, has a long legacy of yielding valuable datasets to test carbon and water functionality of land surface models using stable isotopes. Niwot Ridge has a 17-year record of eddy covariance fluxes of carbon, water, and energy, as well as environmental data (Hu et al., 2010;Monson et al., 2002) and a 10-year record of δ 13 C of CO 2 in forest air (Schaeffer et al., 2008). From a carbon balance perspective, Niwot Ridge is representative of subalpine forests in western North America that, in general, act as a carbon sink to the atmosphere (Desai et al., 2011). Western forests make up a significant portion of the carbon sink in the United States (Schimel et al., 2002), yet this sink is projected to weaken with projected changes in temperature and precipitation (Boisvenue and Running, 2010).
The Community Land Model (CLM), the land subcomponent of the Community Earth System Model (CESM) has a comprehensive representation of biogeochemical cycling (Oleson et al., 2013) that can be applied across a range of temporal (hours to centuries) and spatial (site to global) scales. A mechanistic representation of photosynthetic discrimination based upon diffusion and enzymatic fractionation (Farquhar et al., 1989) was included in the latest release of CLM4.5 (Oleson et al., 2013), and is similar to the formulation implemented in other land surface models Scholze et al., 2003;Wingate et al., 2010;van der Velde et al., 2013). An early version of CLM simulated carbon (but not carbon isotope) dynamics at Niwot Ridge with reasonable skill (Thornton et al., 2002).
Here, we evaluate the performance of the 13 C / 12 C isotope discrimination submodel within CLM4.5 against a range of isotopic observations at Niwot Ridge, to examine what new insights an isotope-enabled model can bring upon ecosystem function. Specifically, we test whether CLM simulates the expected isotopic response to environmental drivers of CO 2 fertilization, soil moisture, and atmospheric vapor pressure deficit (VPD). A previous analysis at Niwot Ridge showed a seasonal correlation between VPD and photosynthetic discrimination (Bowling et al., 2014) suggesting that leaf stomata are responding to changes in VPD, and influencing discrimination. We use CLM to test whether VPD is the primary environmental driver of isotopic discrimination, as compared to soil moisture and net assimilation rate. Next, we determine whether site-specific boundary conditions (including, for instance, δ 13 C of atmospheric CO 2 ) combined with the representation of long-term (multidecadal to century) photosynthetic discrimination and simulated carbon pool turnover within the model, can accurately reproduce the measured δ 13 C in leaf tissue, roots and soil carbon. We then use CLM to determine if the increase in atmospheric CO 2 since 1850 has led to an increase in water-use efficiency (WUE), and whether net assimilation or stomatal conductance is the primary driver of such a change. Finally, we ask what distinct insights site-level isotope observations bring in terms of both model parameterization (i.e., stomatal conductance) and model structure as compared to the traditional observations (e.g., carbon fluxes, biomass). We focus the description of CLM4.5 (Sect. 2.1) upon photosynthesis, and its linkage to nitrogen, soil moisture, and stomatal conductance (Sect. 2.1.1). Next, we describe the model representation of carbon isotope discrimination by photosynthesis (Sect. 2.1.2). Because preliminary simulations demonstrated that model results were strongly influenced by nitrogen limitation, we used three separate nitrogen formulations (described in Sect. 2.1.2) to better diagnose model performance. Next, to provide context for subsequent descriptions of site-specific model adjustments we describe the field site, Niwot Ridge, including the site-level observations (Sect. 2.2) used to constrain and test the model. Patterns in plant growth and δ 13 C of biomass are strongly influenced by atmospheric CO 2 and δ 13 C of atmospheric CO 2 (δ atm ). Therefore, we designed a site-specific synthetic atmospheric CO 2 product (Sect. 2.3.1) and δ atm product (Sect. 2.3.2) for these simulations. The model setup and initialization procedure, intended to bring the system into steady state, is described in Sect. 2.3.3. This is followed by an explanation of the model calibration procedure that provided a realistic simulation of carbon and water fluxes (Sect. 2.4).

Community Land Model version 4.5
We used the Community Land Model version 4.5 (Oleson et al., 2013), which is the land component of the Community Earth System Model (CESM) version 1.2 (www.cesm. ucar.edu/models/cesm1.2/). Details regarding the Community Land Model can be found in Mao et al. (2016) and Oleson et al. (2013). Here, we emphasize the mechanistic formulation that controls photosynthetic discrimination ( canopy ) and factors that influence canopy including photosynthesis, stomatal conductance, water stress, and nitrogen limitation. A list of symbols is provided in Table 1.

Net Photosynthetic Assimilation
The leaf-level net carbon assimilation of photosynthesis, A n , is based on Farquhar et al. (1980) as where A c , A j , and A p are the enzyme (RuBisCO)-limited, light-limited, and product-limited rates of carboxylation, respectively, and Resp d is the leaf-level respiration. The enzyme limited rate is defined as where c i is the intercellular leaf partial pressure of CO 2 , o i = 0.209 P atm , P atm is atmospheric pressure, and K c , K o , and * are constants. The maximum rate of carboxylation at 25 • C, V cmax25 , is defined as where N a is the nitrogen concentration per leaf area, F LNR the fraction of leaf nitrogen within the RuBisCO enzyme, F NR the ratio of total RuBisCO molecular mass to nitrogen mass within RuBisCO, and a R25 is the specific activity of RuBisCO at 25 • C. The V cmax25 is adjusted for leaf temperature to provide V cmax in Eq.
(2), used in the final photosynthetic calculation. Both A j and A p are functions of V cmax as well (not shown). The variable β t represents the level of soil moisture availability, which influences both V cmax (Sellers et al., 1996), and stomatal conductance (Eq. 5). CLM calculates β t as a factor (0-1, high to low stress) by combining soil moisture, the rooting depth profile, and a plant-dependent response to soil water stress as where w i is a plant wilting factor for soil layer i and r i is the fraction of roots in layer i. The plant wilting factor is scaled according to soil moisture and water potential, depending on plant functional type (PFT). Soil moisture is predicted based upon prescribed precipitation and vertical soil moisture dynamics (Zeng and Decker, 2009). The root fraction in each soil layer depends upon a vertical exponential profile controlled by PFT-dependent root distribution parameters adopted from Zeng (2001). The carbon and water balance are linked through c i by the stomatal conductance to CO 2 , g s , following the Ball-Berry model (Ball et al., 1987) as defined by Collatz et al. (1991): where m is the stomatal slope, c s the partial pressure of CO 2 at the leaf surface, h s the relative humidity at the leaf surface, and b the minimum stomatal conductance when the leaf stomata are closed. The version of CLM used here has a two-layer (shaded, sunlit) representation of the vegetation (Oleson et al., 2013). Photosynthesis and stomatal conductance are calculated separately for the shaded and sunlit portion and the total canopy photosynthesis is the potential gross primary productivity (GPP), CF GPPpot : where LAI is the leaf area index and 12.011 −6 is a unit conversion factor. The total carbon available for new growth allocation (CF avail_alloc ) is defined as CF avail_alloc = CF GPPpot − CF GPP, mr − CF GPP, xs , where the maintenance respiration is derived either from recently assimilated photosynthetic carbon CF GPP, mr or, if photosynthesis is low or zero (e.g., night), the maintenance respiration is drawn from a carbon storage pool (CF GPP, xs ). In contrast, CF alloc , is the actual carbon allocated to growth calculated from the available nitrogen and fixed C : N ratios for new growth (e.g., stem, roots, leaves). The downregulation of photosynthesis from nitrogen limitation, f dreg , is given by The actual, nitrogen-limited GPP is defined as

Photosynthetic carbon isotope discrimination
The canopy-level fractionation factor α is defined as the ratio of 13 C / 12 C within atmospheric CO 2 (R a ) and the products of photosynthesis (R GPP ) as α = R a R GPP . The preference of C3 vegetation to assimilate the lighter CO 2 molecule during photosynthesis is simulated in CLM with two steps: diffusion of CO 2 across the leaf boundary layer and into the stomata, followed by enzymatic fixation to give the leaf-level fractionation factor: .
where c * i and c a are the intracellular and atmospheric CO 2 partial pressure, respectively. The numbers 4.4 and 22.6 represent the diffusional and enzymatic contributions to isotopic discrimination, respectively (Farquhar et al., 1989). The variable c * i (known in CLM as the "revised intracellular CO 2 partial pressure") is marked with an asterisk to indicate the inclusion of nitrogen downregulation defined as where g b is the leaf boundary layer conductance. Equation (11) is a general expression for c * i , where within the model c * i and discrimination are calculated for the sunlit and shaded layer of leaves separately and subject to the local environmental conditions unique to each layer (Oleson et al., 2013). The inclusion of the nitrogen downregulation factor f dreg reflects the two-stage process in which the potential photosynthesis (Eq. 6) and the actual photosynthesis (Eq. 9) are calculated within CLM and prevents a mismatch between the actual photosynthesis and the intracellular CO 2 . This mismatch is a result of the carbon-water (A n − g s ) coupling (Eq. 5) being imposed prior to the effect of nitrogen limitation (Eq. 9), and is an artifact of the model implementation. We also test a separate model formulation (described in detail in the next paragraph) specific to this analysis that imposes nitrogen limitation through the V cmax parameterization and removes the artifact of f dreg .
The sensitivity of preliminary model results to nitrogen limitation led us to test three distinct discrimination formulations ( Fig. 1; Table 2). The limited nitrogen formulation was based on the default version of CLM4.5 and included both nitrogen limitation and the nitrogen downregulation factor within the calculation of c * i as given in Eq. (11). The second, unlimited nitrogen formulation, which we created specifically for this analysis, also follows Eq. (11); however, the vegetation is allowed unlimited access to nitrogen (CF GPPpot = GPP, f dreg = 0) which ignores the nitrogen budget within CLM. We account for the increased productivity in the unlimited nitrogen model simulations by calibrating V cmax (Sect. 2.4). Finally, in the no-downregulation discrimination formulation (also created specifically for this analysis), we included nitrogen limitation, but removed the downregulation factor f dreg within the isotopic discrimination Eq. (11).
In the unlimited nitrogen formulation, we use a different modifier on V cmax25 ( Fig. 1; described in Sect. 2.4 and Figs. S1, S2 in the Supplement) in the calibrated runs to give similar carbon flux, water flux, and biomass as in the other two formulations, such that all three formulations have fluxes and biomass that are similar to what is observed at the site, and which presumably reflect nitrogen limitation. Thus, the distinction between these three formulations can be viewed entirely as when nitrogen limitation is imposed in relation to photosynthesis: (1) after photosynthesis via a downregulation between potential and actual GPP (Eq. 9) that feeds back on the c * i /c a used for isotopic discrimination but not on g s or A n in the limited nitrogen formulation; (2) before photosynthesis via V cmax , which limits photosynthetic capacity affecting both c * i /c a , g s and A n in the unlimited nitrogen formulation; and (3) after photosynthesis with no effect on either the c * i /c a for isotopic discrimination or g s or A n in the no-downregulation discrimination formulation. The downscaled portion of the carbon during nitrogen limitation (CF GPPpot − GPP) is removed from the system and does not appear as a respired flux (Fig. 1). In summary, the limited nitrogen (post-photosynthetic) formulation adjusts the photosynthetic rate by explicitly tracking N availability, whereas the unlimited nitrogen (pre-photosynthetic) formulation takes into account any N limitation through the V cmax parameterization. Because the limited nitrogen formulation reduces A n during the nitrogen downregulation step without explicitly solving for g s , the carbon-water cycle is partially coupled, whereas the unlimited nitrogen formulation is fully coupled.
Carbon isotope ratios are expressed by standard delta notation: where R x is the isotopic ratio of the sample of interest and R VPDB is the isotopic ratio of the Vienna Pee Dee Belemnite standard. The delta notation is dimensionless but expressed in parts per thousand (‰) where a positive (negative) value refers to a sample that is enriched (depleted) in 13 C / 12 C relative to the standard. Because this is the only carbon isotope Potential GPP Maintenance resp. flux (CF GPP,mr ) (eq. 7) Photosynthesis-Stomatal Conductance model: A n (eq. 1) V cmax (eq. S1), limited N V cmax (eq. S2), unlimited N g s (eq. 4)

Potential allocation
Actual allocation CF GPPpot (eq. 7) CF Avail_alloc (eq. 7) CF alloc (eq. 8) Limited N formulation reduces A n Figure 1. A simplified representation within CLM4.5 of assimilation and allocation of carbon for conifer species. The colored boxes and solid arrows represent carbon pools and carbon fluxes, respectively. The clear background boxes represent CLM submodels. N limitation is applied if the available N cannot meet the demand determined by the available carbon for allocation (CF avail_alloc ) and the C : N biomass ratio. The blue and red text and arrows represent the limited and unlimited nitrogen formulations, respectively. The no-downregulation discrimination formulation is exactly the same as the limited N formulation in this schematic.
ratio we are concerned with in this paper, the "13" superscript is omitted for brevity in subsequent definitions using the delta notation. The canopy-integrated photosynthetic discrimination, canopy , is defined as the difference between the δ 13 C of the atmospheric and assimilated carbon, The difference between δ 13 C of the total ecosystem respiration (ER) and GPP fluxes, called the isotope disequilibrium (Bowling et al., 2014), is defined as The ecosystem-level water-use efficiency (WUE) is defined as actual carbon assimilated (GPP) per unit water transpired (E T ) per unit land surface area: The intrinsic water-use efficiency (iWUE) from leaf-level physiological ecology is defined as where A n is the net carbon assimilated per unit leaf area and g s is the stomatal conductance. CLM calculates g s (Eq. 5) for shaded and sunlit portions of the canopy separately, therefore an overall conductance was calculated by weighting the conductance by sunlit and shaded leaf areas.

Niwot Ridge and site-level observations
Site-level observations and modeling were focused on the Niwot Ridge Ameriflux tower (US-NR1), a subalpine conifer forest located in the Rocky Mountains of Colorado, USA. The forest is approximately 110 years old and consists of lodgepole pine (Pinus contorta), Engelmann spruce (Picea engelmannii), and subalpine fir (Abies lasiocarpa). The site is located at an elevation of 3050 m above sea level, with mean annual temperature of 1.5 • C and precipitation of 800 mm, in which approximately 60 % is snow. More site details are available elsewhere (Hu et al., 2010;Monson et al., 2002). Flux and meteorological data were obtained from the Ameriflux archive (http://ameriflux.lbl.gov/). Net carbon exchange (NEE) observations were derived from flux tower measurements based on the eddy covariance method and were partitioned into component fluxes of GPP and ER according to two separate methods described by Reichstein et al. (2005) and Lasslop et al. (2010) using an online tool provided by the Max Planck Institute for Biogeochemistry (http://www.bgc-jena.mpg.de/~MDIwork/ eddyproc/). Seasonal patterns in δ GPP and δ ER were derived from measurements as described by Bowling et al. (2014). Observations of δ 13 C of biomass (Schaeffer et al., 2008) and carbon stocks (Bradford et al., 2008;Scott-Denton et al., 2003) were compared to model simulations. Schaeffer et al. (2008) reported soil, leaf, and root observations specific to each conifer species; however, the observed mean and standard error for all species were used for comparison because CLM treated all conifer species as a single PFT.
2.3 Atmospheric CO 2 , isotope forcing and initial vegetation state

Site-specific atmospheric CO 2 concentration time series
Global average atmospheric CO 2 concentrations increased roughly 40 % from 1850 to 2013 (from 280 to 395 ppm). The standard version of CLM4.5 includes an annually and globally averaged time series of this CO 2 increase; however, this does not capture the observed seasonal cycle of ∼ 10 ppm at Niwot Ridge (Trolier et al., 1996). Therefore, we created a site-specific atmospheric CO 2 time series (Fig. 2) to provide a seasonally realistic atmosphere at Niwot Ridge. From 1968 to 2013 the CO 2 time series was fit to flask observations (Dlugokencky et al., 2015) from Niwot Ridge. Prior to 1968, the CO 2 time series was created by combining the average multiyear seasonal cycle based on the Niwot Ridge flask data to the annual CO 2 product provided by CLM. More details are located in the Supplement.

Customized δ 13 C atmospheric CO 2 time series
As atmospheric CO 2 has increased, the δ 13 C of atmospheric CO 2 (δ atm ) has become more depleted (Francey et al., 1999), and this change continues at Niwot Ridge at −0.25 ‰ per decade (Bowling et al., 2014). The δ atm also varies seasonally, and depends on latitude (Trolier et al., 1996). However, CLM4.5 as released assigned a constant δ 13 C of −6 ‰. We therefore created a synthetic time series of δ atm from 1850 to 2013 (Fig. 2). From 1990 to 2013, the time series was fit to the flask observations (White et al., 2015) as described in Sect. 2.3.1. Prior to 1990, the interannual variation within the δ atm time series was fit to the ice core data from Law Dome (Francey et al., 1999; see also Rubino et al., 2013). This annual data product was then combined with the average seasonal cycle at Niwot Ridge as determined by the flask observations to create the synthetic product from 1850 to 1990. More details about the methods and the site-specific data set of atmospheric CO 2 and δ 13 C of CO 2 are located in the Supplement.

Model initialization
We performed an initialization to transition the model from near-bare ground conditions to present-day carbon stocks and LAI that allowed for proper evaluation of isotopic performance. This was implemented in four stages: (1) accelerated decomposition (1000 model years), (2) normal decomposition (1000 model years), (3) parameter calibration (1000 model years), and (4) transient simulation period . The first two stages were preset options within CLM with the first stage used to accelerate the equilibration of the soil carbon pools, which require a long period to reach steady state (Thornton and Rosenbloom, 2005). The parameter calibration stage was not a preset option but designed specifi-cally for our analysis. For this, we introduced a seasonally varying V cmax that scaled the simulated GPP and ecosystem respiration fluxes to present-day observations (Sect. 2.4). In the transient phase, we introduced time-varying atmospheric conditions from 1850 to 2013 including nitrogen deposition (CLM provided), atmospheric CO 2 , and δ atm (site-specific as described above). Environmental conditions of temperature, precipitation, relative humidity, radiation, and wind speed were taken from the Niwot Ridge flux tower observations from 1998 to 2013 and then cycled continuously for the entirety of the initialization process. We used a scripting framework (PTCLM) that automated much of the workflow required to implement several of these stages in a site-level simulation (Mao et al., 2016;Oleson et al., 2013).

Specific model details and model calibration
This version of CLM included a fully prognostic representation of carbon and nitrogen within its vegetation, litter, and soil biogeochemistry. We used the Century model representation for soil (three litter and three soil organic matter pools) with 15 vertically resolved soil layers (Parton et al., 1987). Nitrification and prognostic fire were turned off. Our initial simulations used prognostic fire, but we found that simulated fire was overactive leading to low simulated biomass compared to observations. Although Niwot Ridge has been subject to disturbance from fire and harvest in the past, ultimately our final simulations did not include either fire or harvest disturbance because the last disturbance occurred over 110 years ago (early 20th century logging; Monson et al., 2005).
Ecosystem parameter values (Table 3) used here were based upon the temperate evergreen needleleaf PFT within CLM. These values were based upon observations reported by White et al. (2000) intended for a wide range of temperate evergreen forests, and by Thornton et al. (2002) for Niwot Ridge. For this analysis, two site-specific parameter changes were made. First, the e-folding soil decomposition parameter was increased from 5 to 20 m. This parameter is a length scale for attenuation of decomposition rate for the resolved soil depth from 0 to 5 m where an increased value effectively increases decomposition at depth, thus reducing total soil carbon and more closely matching observations. Second, we performed an empirical photosynthesis scaling (Eq. 17, below) that reduced the simulated photosynthetic flux, as guided by eddy covariance observations (Figs. 3,  S1). Consequently, all downstream carbon pools and fluxes, including ecosystem respiration, aboveground biomass, and leaf area index, provided a better match to present-day observations. This approach also removed a systematic overestimation of winter photosynthesis. The model simulations without the photosynthetic scaling are referred to within the text and figures as the uncalibrated model, whereas model simulations that include the photosynthetic scaling are referred to as the calibrated model. We modified CLM for this where f df is the photosynthetic scaling factor, and all other parameters are identical to Eq. (3). These parameters were constant for the entirety of the simulations except for f df , an empirically derived time-dependent parameter ranging from 0 to 1. The value was set to zero to force photosynthesis to zero between 13 November and 23 March, consistent with flux tower observations where outside of this range GPP > 0 was never observed.
where the observed GPP was the daily average calculated from the partitioned flux tower observations (Reichstein et al., 2005) from 2006 to 2013, and the simulated GPP was the daily average of the unscaled value during the same time.
A polynomial was fit to Eq. (18) that represented f df for (1) both the limited nitrogen and no downregulation discrimination formulations and (2) the unlimited nitrogen formulation (Fig. S2). Note that CLM already includes a day length factor that also adjusts the magnitude of V cmax according to time of year; however, that default parameterization alone was not sufficient to match the observations. The lightlimited rate and product-limited rate of carboxylation (A j , A p ; Eq. 1) and maintenance leaf respiration are functions of V cmax (not shown) and are therefore subject to the same calibration. The CLM model (limited nitrogen simulation) was successful at simulating GPP, ER, and latent heat fluxes (Fig. 3), leaf area index (LAI), and aboveground biomass (Fig. 4), but only following site-specific calibration. Similar improvement was observed after calibration for the unlimited nitrogen run (not shown). The calibration also eliminated erroneous winter GPP. In general, terrestrial carbon models tend to overestimate photosynthesis during cold periods for temperate/boreal conifer forests (Kolari et al., 2007), including Niwot Ridge (Thornton et al., 2002). Although our calibration approach forced V cmax to zero during the winter, it did not solve the underlying mechanistic shortcoming. A more fundamental approach should address either cold inhibition (Zarter et al., 2006) of photosynthesis or soil water availability associated with snowmelt (Monson et al., 2005) to achieve the photosynthetic reduction. Nevertheless, within the confines of our study area, our calibration approach was sufficient to provide a skillful representation of photosynthesis and provided a sufficient testbed for evaluating carbon isotope behavior. We caution that because the f df parameter (Eq. 17) was calibrated specifically for Niwot Ridge, it would not be applicable outside this study area.

δ 13 C of carbon pools
The model performed better at simulating δ 13 C biomass of bulk needle tissue, roots, and soil carbon (Fig. 5) for the unlimited nitrogen and no downregulation discrimination cases as compared to the limited nitrogen case. When nitrogen limitation was included the model underestimated δ 13 C of sunlit needle tissue (1.8 ‰), bulk roots (1.0 ‰), and organic soil carbon (0.7 ‰). All simulations fell within the observed range of δ 13 C in needles that span from −28.7 ‰ (shaded) to −26.7 (sunlit). This vertical pattern in δ 13 C of leaves is common (Martinelli et al., 1998) and results from vertical differences in nitrogen allocation and photosynthetic capacity. The model results integrated the entire canopy and ideally should be closer to sun leaves (as in Fig. 5) given that the majority of photosynthesis occurs near the top of the canopy. Model simulations of δ 13 C of living roots were ∼ 1 ‰ more negative as compared to the structural roots. This range in δ 13 C results from decreasing δ atm with time (Suess effect, Fig. 2). The living roots had a relatively fast turnover time of carbon within the model, whereas the structural roots had a slower turnover time and reflected an older (more enriched δ atm ) atmosphere. The limited nitrogen simulation was a poor match to observations relative to the others (Fig. 5b).
There was an observed vertical gradient in δ 13 C of soil carbon (−24.9 to −26 ‰) with more enriched values at greater depth (Fig. 5c). This vertical gradient is commonly observed (Ehleringer et al., 2000). Simulated δ 13 C of soil carbon was most consistent with the organic horizon observations. There are a wide variety of post-photosynthetic fractionation processes in the soil system Brüggemann et al., 2011) that are not considered in the CLM4.5 model, so the match with observations is perhaps fortuitous.

Decadal changes in photosynthetic discrimination and driving factors
All modeled carbon pools showed steady depletion in δ 13 C since 1850 (coinciding with the start of the transient phase of simulations, Fig. 5). For the limited nitrogen run, there was a decrease in δ 13 C of 2.3 ‰ for needles, 2.3 ‰ for living for the limited nitrogen simulation. The observations are taken from the Ameriflux L2 processed eddy covariance flux tower data, partitioned into GPP and ER using the method of Reichstein et al. (2005). The uncalibrated simulation represents the CLM simulation without V cmax scaling and the calibrated simulation represents the CLM run using the V cmax scaling approach.
roots, and 0.1 ‰ for soil carbon. This occurred because of (1) decreased δ atm (Suess effect, Fig. 2) and (2) increased photosynthetic discrimination. We quantified the contribution of the Suess effect by performing a control run with constant δ atm , and kept other factors the same (Fig. 6). Approximately 70 % of the reduction in δ 13 C of needles occurred due to the Suess effect, and the remaining 30 % was caused by increased photosynthetic discrimination. This occurred as plants responded to CO 2 fertilization as illustrated in Fig. 7. The model indicated that plants responded to increased atmospheric CO 2 (∼ 40 % increase) by decreasing stomatal conductance (Eq. 5) by 20 % for the limited nitrogen run and 30 % for the unlimited nitrogen run (Fig. 7b) with associated change in c * i /c a (Fig. 7a). Other influences upon stomatal conductance were less significant, including A n (+10 % limited nitrogen, −10 % unlimited nitrogen; Fig. 7d), soil moisture availability (2-3 %; Fig. 7e), and negligible changes in relative humidity (multidecadal climate change effects are neglected due to methodological cycling of weather data). This finding that stomatal conductance responded to atmo- spheric CO 2 is consistent with both tree ring studies (Saurer et al., 2014) and site-level experiments (Ward et al., 2012).
The effect of CO 2 fertilization and associated response of stomatal conductance and net assimilation led to a multidecadal increase in c * i /c a for all model formulations (Fig. 7a). The c * i /c a increased from 0.71 to 0.76, 0.67 to 0.71, and 0.66 to 0.68 for the limited nitrogen, unlimited nitrogen, and no-downregulation discrimination formulations, respectively, from 1850 to 2013. All simulations therefore suggested an increase in photosynthetic discrimination. This increase in discrimination falls in between two hypotheses posed by Saurer et al. (2004) regarding stomatal response to increased CO 2 : (1) reduction in stomatal conductance causes c i to proportionally increase with c a keeping c i /c a constant and (2) minimal stomatal conductance response where c i increases at the same rate as c a (constant c a − c i ) causing c i /c a to increase. Our simulation generally agrees with the observed trend in c i /c a as estimated from tree ring isotope measurements from a network of European forests (Frank et al., 2015). When controlled for trends in climate, Frank et al. (2015) found that c i /c a was approximately constant during the last century. If the Niwot Ridge multidecadal warming trends in temperature and humidity (Mitton and Ferrenberg, 2012) were included in the CLM simulations (this analysis did not consider multidecadal climate change) the stom-  Table 2. Uncertainty bars for simulations represent 95 % confidence intervals of δ 13 C variation as a result of cycling the site-level meteorology observations. The observed values are from Schaeffer et al. (2008) with uncertainty bars representing standard error. Solid lines and dashed lines in middle panel represent living roots and structural roots, respectively. Figure 6. Simulation of δ 13 C of needle tissue using the limited nitrogen (default) CLM run. In the constant δ 13 C of CO 2 (δ atm ) simulation the model boundary condition was −6 ‰, whereas the transient δ atm simulation varied over time (Fig. 2). atal response may have been stronger, thereby holding c i /c a constant.
The simulated stomatal closure in response to CO 2 fertilization led to an increase in iWUE and WUE of approximately 20 and 10 %, respectively (Fig. 7f), from 1960 to 2000. This simulated increase in iWUE is consistent with the observation-based studies (Ainsworth and Long, 2005;Franks et al., 2013;Peñuelas et al., 2011) which indicate a 15-20 % increase in iWUE for forests during that time. The overall increase in WUE suggests that the vegetation at Niwot Ridge has some ability to maintain net ecosystem productivity when confronted with low soil moisture, low humidity conditions. Ultimately, whether Niwot Ridge maintains the current magnitude of carbon sink (Figs. 3, S1) will depend upon the severity of drought conditions, as improve-ments in WUE, in general, are only likely to negate weak to moderate levels of drought (Franks et al., 2013).
The limited nitrogen formulation simulated larger values of A n and g s , and smaller iWUE as compared to the unlimited nitrogen formulation (Fig. 7). This is because the unlimited nitrogen formulation was fully coupled (i.e., solved simultaneously) between A n and g s (Eq. 5). The limited nitrogen formulation, however, was only partially coupled because A n and g s were initially solved simultaneously through the potential A n (Eq. 1); however, under N limitation, A n becomes limited below its potential value (Eq. 9) through f dreg . Therefore, g s is calculated through the potential A n (Eq. 5) and not the nitrogen-limited A n .
The simultaneous increase in both simulated photosynthetic discrimination and iWUE conflicts with previous literature where increases in iWUE are typically linked with weakening discrimination (e.g., Saurer et al., 2004) using a linear model. In general, an increase in atmospheric CO 2 alone tends to increase iWUE because of reduced stomatal conductance; however, the impact upon discrimination is close to neutral because the increased supply of CO 2 external to the leaf is offset by reduced stomatal conductance (Saurer et al., 2004). The VPD likely plays an important role in determining the final trends for iWUE and discrimination, where an increasing VPD should further reduce stomatal conductance (VPD α 1 h s ; Eq. 5) thereby promoting the well-established relationship (increasing iWUE, decreasing discrimination). In contrast, a weak or decreasing trend in VPD should promote the opposite relationship (increasing iWUE, increasing discrimination).
The CLM model at present neglects mesophyll conductance (g m ). When Seibt et al. (2008) included g m in a model that linked iWUE to discrimination, they found there were certain conditions when iWUE and discrimination increased together. This is in part because mesophyll conductance, unlike stomatal conductance, does not respond as strongly to Figure 7. Diagnostic model variables that explain the discrimination trends (Fig. 5) for the three model formulations as described in Table 2 for (a) c * i /c a , (b) g s , (c) f dreg , (d) A n , (e) β t , and (f) the WUE and iWUE. Where the no-downregulation discrimination simulation is not shown, it was identical to the limited nitrogen simulation. Uncertainty bars represent 95 % confidence intervals of diagnostic variable variation as a result of cycling the site-level meteorology observations. The dashed lines represent WUE and the solid lines represent iWUE in (f). changes in VPD, yet has a significant impact upon c i /c a and discrimination (Flexas et al., 2006). Harvard Forest is an example of a site that was observed to show simultaneous increase in iWUE and discrimination over the last 2 decades, using data derived from tree rings (Belmecheri et al., 2014). In our model simulation, we do not consider multidecadal trends in climate or mesophyll conductance; therefore, increasing atmospheric CO 2 must be the primary driver for the modeled simultaneous increase in discrimination and iWUE at Niwot Ridge (Fig. 7). These trends in iWUE and discrimination have also been found in a fully coupled, isotope enabled, global CESM1.2 model run with climate simulated by CAM5 (Community Atmosphere Model) driven by CO 2 emissions (unpublished, K. Lindsay; Fig. S3). Specifically, a random sample of land model grid cells representing conifer species in British Columbia (lat: 52.3 • N, long: −122.5 • W) and Quebec (lat: 49.5 • N, long: −70.0 • W) all showed an increase in photosynthetic discrimination and a 10 % increase in WUE from 1850 to 2005. These randomly chosen grid cells are likely better analogs to the site-level simulations described here because they represent boreal conifer forests, whereas the grid cells that are in the Niwot Ridge area were heterogeneous in land cover (e.g., tundra, grassland, forest) and a poor representation of conifer forest.
The relationship between iWUE and discrimination in the global CESM1.2 model run (with model-simulated climate) suggest that the site-level trends are not isolated to the spe-cific conditions of Niwot Ridge, but are a function of the model formulation. There is a relationship between iWUE and c * i /c a (discrimination) as derived from Eq. (11) within the CLM model: The full derivation is provided in the Supplement. Note that according to Eq. (19) increasing iWUE can be consistent with weakening discrimination (decreasing c * i /c a (∼ α)) and therefore consistent with established understanding between trends in iWUE and discrimination. However, this can be moderated by increasing c a . During the course of our simulation (1850-2013), iWUE increased between 15 and 20 % (Fig. 7); however, c a increased by 40 %.

Magnitude of photosynthetic discrimination
The simulated photosynthetic discrimination (Fig. 8) was significantly larger than an estimate derived from observations and an isotopic mixing model (Bowling et al., 2014). For brevity, we refer to the estimates based on the Bowling et al. (2014) method as observed discrimination but highlight that they are derived from observations and not directly measured. On average, the simulated monthly growing season mean canopy discrimination was greater than observed values by 4.0, 2.3, and 1.8 ‰ for the limited nitrogen, unlimited nitrogen, and no-downregulation discrimina-tion formulations, respectively. The model-observation mismatch in discrimination, despite model-observation agreement to biomass, carbon, and latent heat flux tower observations (Fig. 3), highlights the independent and useful constraint isotopic observations provide for evaluating model performance. Specifically, the overestimation of discrimination may suggest the stomatal slope in the Ball-Berry model (m = 9 in Eq. 5) used for these simulations was too high. This is supported by Mao et al. (2016), who found a reduced stomatal slope (m = 5.6) was necessary for CLM4.0 to match observed δ 13 C in an isotope labeling study of loblolly pine forest in Tennessee. The stomatal slope was also important to match discrimination behavior in the ISOLSM model (Aranibar et al., 2006), a predecessor to CLM. A global analysis of stomatal slope inferred from leaf gas exchange measurements found that evergreen coniferous species, such as those at Niwot Ridge, had near the lowest values compared to other PFTs (Lin et al., 2015). In addition, they found that low stomatal slope values were characteristic of species with low stemwood construction costs per water transpired (high WUE), low soil moisture availability, and cold temperatures.
Alternatively, discrimination may be overestimated because CLM does not consider the resistance to CO 2 diffusion into the leaf chloroplast. The ability of CO 2 to diffuse across the chloroplast boundary layer, cell wall, and liquid interface is collectively known as the mesophyll conductance (g m ) (Flexas et al., 2008). Multiple studies suggest that g m is comparable in magnitude to g s , and responds similarly to environmental conditions (Flexas et al., 2008). CLM does not account for g m , and as a result assumes the intracellular CO 2 is the same as intercellular CO 2 , when it can be significantly lower (Di Marco et al., 1990;Sanchez-Rodriguez et al., 1999). The overestimation of c * i could have two important impacts upon our simulation. First, this may lead to unrealistically low values of V cmax in order to compensate for the overestimation of c * i . In fact, we reduced the default value of V cmax as much as 50 % in our simulation to match the eddy covariance flux tower observations (see Sect. 4.1). Second, the overestimation of c * i should cause an overestimation of discrimination (Eq. 10), which is also consistent with our simulations (Fig. 8). To determine whether the simulated discrimination bias is a model parameter calibration issue (g s ) or from excluding g m , we recommend a mechanistic representation of mesophyll conductance within CLM.
The mixing model approach estimate of canopy (17 ‰) (Bowling et al., 2014), combined with δ atm (−8.25 ‰) implies a δ 13 C of biomass between −26 and −25 ‰ (Fig. 8). This range is only slightly more enriched than the observed ranges of δ 13 C of needle and root biomass (−27 to −26 ‰). The fact that the different approaches to measure discrimination differ by only 1 ‰, whereas CLM simulates a canopy that is 1.8 to 4.0 ‰ greater than the mixing model discrimination, strongly suggests that the model has overestimated discrimination from 2006 to 2012. Therefore, what appeared to be a successful match between the simulated and observed δ 13 C biomass, may in fact have been fortuitous. A multidecadal time series of discrimination inferred from δ 13 C of tree rings (Saurer et al., 2014;Frank et al., 2015) would be useful to investigate this mismatch as a function of time, but these data are not presently available.
If the overestimation of modeled discrimination originates from a lack of response of stomatal conductance to environmental conditions, this could be a result of one or several of the following within the model: (1) the stomatal slope value is too high, (2) multidecadal trends in climate (e.g., VPD) have not been included in the simulation, (3) the model neglects g m , or (4) the Ball-Berry representation of g s is not sensitive enough to changes in environmental conditions (e.g., humidity, soil moisture). It has been shown that VPD may be an improved predictor of g s (Katul et al., 2000;Leuning, 1995) and discrimination (Ballantyne et al., 2010(Ballantyne et al., , 2011 as compared to relative humidity, currently used in CLM4.5. Future work should consider which of these scenarios is responsible for overestimation of discrimination.

Seasonal pattern of photosynthetic discrimination
The model formulations that did not explicitly consider the influence of nitrogen limitation upon discrimination (unlimited nitrogen, no downregulation discrimination) were most successful at reproducing the seasonality of discrimination (Figs. 8, S4). In general, the observed discrimination was stronger during the spring and fall and weaker during summer. This observed canopy seasonal range (excluding November) varied from 16.5 to 18 ‰ using Reichstein partitioning (Fig. 8), and was more pronounced using Lasslop partitioning (16.5 to 23 ‰) (Fig. S4). The nitrogen-limited simulated canopy had no seasonal trend, whereas the unlimited nitrogen and no-downregulation discrimination simulations ranged from 18.4 to 21.2 and 17.8 to 20.6 ‰, respectively.
The main driver of the seasonality of discrimination was the net assimilation (A n ) for the unlimited nitrogen formulation (Fig. 9). This was evident given the inversely proportional relationship between the simulated fractionation factor (α) and A n , consistent with Eq. (11). Stomatal conductance (g s ) also influenced the seasonal pattern. The most direct evidence for this was during the period between days 175 and 200 (Fig. 9), where A n descended from its highest value (favoring higher α), and g s abruptly ascended to its highest value (favoring higher α). The α responded to this increase in g s with an abrupt increase by approximately 0.003 (3 ‰). Similarly, the limited nitrogen simulation seasonal discrimination pattern was shaped by both A n and g s , although the magnitude for both was approximately 30 % higher during the summer months as compared to the unlimited nitrogen simulation. This was because the calibrated V cmax value for the limited nitrogen simulation was much higher than for the unlimited nitrogen simulation (Sect. 4.1). The difference in α between the two model formulations coincided with  the sharp increase in f dreg between days 125 and 275, providing strong evidence that the downregulation mechanism within the limited nitrogen formulation led to increased discrimination during the summer. Therefore, it follows that the nitrogen downregulation mechanism was the root cause of the small range in simulated seasonal cycle discrimination for the limited nitrogen formulation, which was inconsistent with the observations.  (Bowling et al., 2014). The horizontal lines represent δ 13 C of 17 ‰ and are included for reference.

Environmental factors influencing seasonality of discrimination
The simulated canopy was driven primarily by net assimilation (A n ), followed by vapor pressure deficit (VPD) (Fig. 10). The correlation between VPD and canopy was strongest for the unlimited nitrogen simulation, where the range in monthly average canopy spanned values from 18 to 22 ‰ (Fig. 10b, e, and h). This resembled the observed range in response based upon a fitted relationship from Bowling et al. (2014) that spanned from roughly 16 to 19 ‰ (Fig. 10a, b, and c), although with a consistent discrimination bias. The correlation between VPD and canopy , however, does not demonstrate causality. If that were the case, given that g s is a function of VPD (h s term in Eq. 5) and discrimination is a function of g s (Eqs. 10, 11), a similar relationship should have existed between g s and canopy . This, in fact, was not the case. Overall, the influence of g s (responding to VPD) (R value = −0.50) was secondary to A n (R value = −0.77) in driving changes in discrimination (Fig. 10). The model suggested that the range in seasonal discrimination (intraannual variation) was driven by the magnitude of A n based on the inverse relationship between A n and canopy , (Eq. 11) illustrated by the separation between months of low photosynthesis (October, May) vs. high photosynthesis (June, July, August). During times of relatively low photosynthesis, A n also drove the interannual variation in canopy . On the other hand, g s (VPD) was most influential in driving the interannual variation of discrimination during the summer months only, judging by the directly proportional relationship during the months of June, July, and August. Strictly speaking, g s is a function of h s (leaf relative humidity) and not atmospheric VPD in CLM. However, the two are closely related and the relationship between either variable (atmospheric VPD or simulated leaf humidity) to canopy was similar (Fig. S5). The limited nitrogen formulation did not produce as wide a range in discrimination as compared to the observations (Fig. 10a, d, and g). Part of this result was attributed to the lack of response between A n and canopy . In this case, the discrimination did not decrease with increasing A n because the signal was muted by the countering effect of f dreg . The limited nitrogen formulation was, however, able to reproduce the same discrimination response to g s as compared to the other model formulations. The tendency for the limited nitrogen model to simulate discrimination response to g s and not to A n may negatively impact its ability to simulate multidecadal trends in discrimination. This may not be a major detriment to sites such as Niwot Ridge which have maintained a consistent level of carbon uptake during the last decade, and is likely more susceptible to environmental impact upon stomatal conductance. However, sites that have shown a significant increase in assimilation rate (e.g., Harvard Forest; Keenan et al., 2013) are less likely to be well represented by this model formulation.
Given the dependence of forest productivity at Niwot Ridge on snowmelt (Hu et al., 2010), it was surprising that the model simulated minimal soil moisture stress (Fig. 9e) and therefore minimal discrimination response to soil moisture. However, this finding was consistent with Bowling et al. (2014), who did not find an isotopic response to soil mois-5198 B. Raczka et al.: An observational constraint on stomatal function in forests ture. In addition, lack of response to change in soil moisture may not be indicative of poor performance of the isotopic submodel performance, but rather an effect of the hydrology submodel (Duarte et al., 2016). However, a comparison of observed soil moisture at various depths at Niwot Ridge generally agrees with the CLM-simulated soil moisture (not shown), suggesting the lack of model response to soil moisture was not from biases in the hydrology model.

Discrimination formulations: implications for model development
The limited and unlimited model formulations tested in this study represented two approaches to account for nitrogen limitation within ecosystem models. The limited nitrogen formulation reduced photosynthesis, after the main photosynthesis calculation, so that the carbon allocated to growth was accommodated by available nitrogen. This allocation downscaling approach is common to a subset of models, for example, CLM (Thornton et al., 2007), DAYCENT (Parton et al., 2010) and ED2.1 (Medvigy et al., 2009). Another class of models limits photosynthesis based upon foliar nitrogen content and adjusts the photosynthetic capacity through nitrogen availability in the leaf through V cmax (e.g., CABLE, GDAY, LPJ-GUESS, OCN, SDVGM, TECO; see Zaehle et al. (2014). These foliar nitrogen models are similar to the unlimited nitrogen formulation of CLM because the scaling of photosynthesis was taken into account in the V cmax scaling methodology (see discussion in Sect. 2.1.2 and 2.4), prior to the photosynthesis calculation. In general, there were no categorical differences in behavior between these two classes of models during CO 2 manipulation experiments held at Duke Forest and ORNL (Zaehle et al., 2014). However, CLM4.0 was one of the few models in that study to consistently underestimate the NPP response to an increase of atmospheric CO 2 due to nitrogen limitation. This finding was attributed to a lower initial supply of nitrogen and too strong of a coupling between carbon and nitrogen that limited biomass production. Also within this experiment, it was found that models that had no or partial coupling (CLM4.0, DAYCENT) between A n and g s , generally predicted lower than observed WUE response to increases in CO 2 (De Kauwe et al., 2013). Similar to CLM4.0, the limited nitrogen formulation of CLM4.5 in this paper is partially coupled (see Sect The unlimited nitrogen formulation described in our study has similarities to a foliar nitrogen model, in that, the influence of nitrogen limitation is parameterized within V cmax . A true foliar nitrogen model, however, couples a dynamic ni-trogen cycle directly with the calculation of V cmax . This capability was recently developed within CLM (Ghimire et al., 2016) and is scheduled to be included in the next CLM release. Future work should test its functionality.
The performance of the unlimited nitrogen formulation was nearly identical to the no-downregulation discrimination formulation in terms of isotopic behavior despite the mechanistic differences. The no-downregulation discrimination formulation included nitrogen limitation within the bulk carbon behavior but ignored the impact of f dreg upon discrimination behavior. The relative high simulation skill with this formulation implied that the potential GPP linked to A n , was a more effective predictor of discrimination behavior than the downscaled GPP, which is linked to A n* (1 − f dreg ) (Eq. 11). There are several potential explanations for an unrealistically large value of f dreg . First, this could indicate that the V cmax parameter was too large, thereby requiring a large f dreg to compensate. As noted in Sect. 3.1, the default temperate evergreen V cmax25 was ∼ 62 µmol m −2 s −1 , much larger than what was found based on literature reviews (Monson et al., 2005;Tomaszewski and Sievering, 2007). We found to match the observed GPP we had to impose f dreg that had the same effect as reducing V cmax (Fig. S2) to values of 51 and 34 µmol m −2 s −1 for the limited nitrogen and unlimited nitrogen formulations, respectively. Alternatively, it could be that there are physiological processes that are acting to reduce nitrogen limitation (e.g., nitrogen storage pools or transient carbon storage as nonstructural carbohydrates), or that the current measurement techniques are underestimating GPP due to biases within the flux partitioning methods.

Disequilibrium, possible explanations of mismatch
Carbon cycle models (e.g., Fung et al., 1997) indicate that the steady decrease of δ atm (Suess effect, Fig. 2) should lead to a positive disequilibrium between land surface processes (δ 13 C difference between GPP and ER; Eq. 14). This is because the δ GPP reflects the most recent (δ 13 C-depleted) state of the atmosphere, whereas the δ ER reflects carbon (e.g., soil carbon) assimilated from an older (δ 13 C-enriched) atmosphere. This positive disequilibrium pattern promoted by the Suess effect was consistent with all CLM formulations for this study with an annual average disequilibrium of 0.8 ‰. In contrast, a negative disequilibrium (−0.6 ‰ ) was identified at Niwot Ridge based upon observations (Bowling et al., 2014) as well as in other forests Wehr and Saleska, 2015;Wingate et al., 2010). Bowling et al. (2014) hypothesized several reasons for this: (1) a strong seasonal stomatal response to atmospheric humidity, (2) decreased photosynthetic discrimination associated with CO 2 fertilization, (3) decreased photosynthetic discrimination associated with multidecadal warming and increased VPD, and (4) post-photosynthetic discrimination. We evaluated the first three hypotheses within the context of our CLM simulations.
The model results suggest a seasonal variation of discrimination that is a function of both VPD and A n . The simulated seasonal range in discrimination (Figs. 8,S4) varied by approximately 2 ‰, and this range in seasonal discrimination could contribute to a negative disequilibrium provided specific timing of assimilation, assimilate storage and respiration not currently considered in the model. For example, if a significant portion of photosynthetic assimilation was stored during the spring with relatively high discrimination and then respired during the summer, the net effect would deplete the δ ER and thereby promote negative disequilibrium during the summer months when discrimination is lower. Theoretically, this could be achieved by explicitly including carbohydrate storage pools within CLM. Isotopic tracer studies have shown assimilated carbon can exist for weeks to months within the vegetation and soil before it is finally respired (Epron et al., 2012;Hogberg et al., 2008). Although carbon storage pools are included in CLM, their allocation is almost always instantaneous for evergreen systems and could not provide the isotopic effect described above (Mao et al., 2016;Duarte et al., 2016).
The CO 2 fertilization effect tends to favor photosynthesis in plants and has been shown to simultaneously increase WUE and decrease stomatal conductance as inferred from δ 13 C in tree rings (Frank et al., 2015;Flanagan et al., 2012;Wingate et al., 2010). In general, a decrease in stomatal conductance and increase in WUE is associated with a decrease in C3 discrimination (Farquhar et al., 1982), which opposes the disequilibrium trend imposed by the Suess effect. The model simulation agrees with both these trends in WUE and stomatal conductance, yet simulates an increase in discrimination (Figs. 6, 7), which reinforces the Suess effect pattern upon disequilibrium. Although this appears to be a mismatch between forest processes and model performance, the model is operating within the limits of the discrimination parameterization (Eq. 17) in which the magnitude of photosynthetic discrimination is inversely proportional to the iWUE, but is also proportional to atmospheric CO 2 (see Sect. 3.2.1).
A multidecadal decrease in photosynthetic discrimination may also result from change in climate. Meteorological measurements at Niwot Ridge during the last several decades generally support conditions of higher VPD based upon a warming trend from an average annual temperature of 1.1 • C in the 1980s to 2.7 • C in the 2000s (Mitton and Ferrenberg, 2012) and no overall trend in precipitation. It is possible that a multidecadal trend in increasing VPD contributed to multidecadal weakening in photosynthetic discrimination given the observed (Bowling et al., 2014) and modeled (Fig. 10) correlation between canopy and VPD. The model meteorology only included the years 1998-2013 and did not include the rapid warming after the 1980s. It is unclear whether, if the full period of warming were to be included in the simulation, the simulated discrimination response to VPD would be enough to counter the Suess effect and lead to negative disequilibrium. Still, there is evidence that the model is over-estimating contemporary discrimination (Sect. 4.2) and the exclusion of the full multidecadal shift in VPD could be a significant reason why.
Finally, post-photosynthetic discrimination processes are likely to impact the magnitude and sign of the isotopic disequilibrium Brüggemann et al., 2011) at multiple temporal scales. None of these isotopic processes are currently modeled within CLM4.5, so at present the model cannot be used to examine them.

Conclusions
This study provides a rigorous test of the representation of C isotope discrimination within the mechanistic terrestrial carbon model CLM. CLM was able to accurately simulate δ 13 C in leaf and stem biomass and the seasonal cycle in canopy , but only when V cmax was calibrated to account for nitrogen limitation prior to photosynthesis (unlimited nitrogen formulation).
Although the unlimited nitrogen formulation (fully coupled carbon and water cycle) was able to match observed δ 13 C of biomass and seasonal patterns in discrimination, it still overestimated the contemporary magnitude of discrimination (2006)(2007)(2008)(2009)(2010)(2011)(2012). Future work should identify whether this overestimation was a result of parameterization (stomatal slope), exclusion of multidecadal shifts in VPD, limitations in the representation of stomatal conductance (Ball-Berry model), or absence of the representation of mesophyll conductance.
The model attributed most of the range in seasonal discrimination to variation in net assimilation rate (A n ) followed by variation in VPD, with little to no impact from soil moisture. The model suggested that A n drove the seasonal range in discrimination (across-month variation), whereas VPD drove the interannual variation during the summer months. This finding suggests that to simulate multidecadal trends in photosynthetic discrimination, response to assimilation rate and VPD must be well represented within the model.
The model simulated a positive disequilibrium that was driven by both the Suess effect and increased photosynthetic discrimination from CO 2 fertilization. It is possible that the negative disequilibrium that was inferred from observations (Bowling et al., 2014) was driven from the impacts of climate change and/or post-photosynthetic discrimination -not considered in this version of the model.
The model simulated a consistent increase in water-use efficiency as a response to CO 2 fertilization and decrease in stomatal conductance. The model simulated an increase in WUE despite an increase in discrimination; however, C3 plants typically express the opposite trends (increase in WUE, decrease in discrimination). Although CLM includes parameterization that promotes an increase in WUE with a decrease in discrimination, this trend was likely moderated by an increase in c a . Initial indications are that δ 13 C isotope data can bring additional constraint to model parameterization beyond what traditional flux tower measurements of carbon, water exchange, and biomass measurements. The isotope measurements suggested a stomatal conductance value generally lower than what was consistent with the flux tower measurements. Unexpectedly, the isotopes also provided guidance upon model formulation related to nitrogen limitation. The success of our empirical approach to account for nutrient limitation within the V cmax parameterization suggests that additional testing of foliar nitrogen models is worthwhile.

Information about the Supplement
All supplemental figures, derivations, and methodological details and synthetic atmospheric data (Sect. 2.3.1, 2.3.2) can be found in the Supplement.
The Supplement related to this article is available online at doi:10.5194/bg-13-5183-2016-supplement.