Evaluating the Community Land Model (CLM4.5) at a coniferous forest site in northwestern United States using flux and carbon-isotope measurements

Abstract

Droughts in the western United States are expected to intensify with climate change. Thus, an adequate representation of ecosystem response to water stress in land models is critical for predicting carbon dynamics. The goal of this study was to evaluate the performance of the Community Land Model (CLM) version 4.5 against observations at an old-growth coniferous forest site in the Pacific Northwest region of the United States (Wind River AmeriFlux site), characterized by a Mediterranean climate that subjects trees to water stress each summer. CLM was driven by site-observed meteorology and calibrated primarily using parameter values observed at the site or at similar stands in the region. Key model adjustments included parameters controlling specific leaf area and stomatal conductance. Default values of these parameters led to significant underestimation of gross primary production, overestimation of evapotranspiration, and consequently overestimation of photosynthetic ${}^{\mathrm{13}}$C discrimination, reflected in reduced ${}^{\mathrm{13}}$C : ${}^{\mathrm{12}}$C ratios of carbon fluxes and pools. Adjustments in soil hydraulic parameters within CLM were also critical, preventing significant underestimation of soil water content and unrealistic soil moisture stress during summer. After calibration, CLM was able to simulate energy and carbon fluxes, leaf area index, biomass stocks, and carbon isotope ratios of carbon fluxes and pools in reasonable agreement with site observations. Overall, the calibrated CLM was able to simulate the observed response of canopy conductance to atmospheric vapor pressure deficit (VPD) and soil water content, reasonably capturing the impact of water stress on ecosystem functioning. Both simulations and observations indicate that stomatal response from water stress at Wind River was primarily driven by VPD and not soil moisture. The calibration of the Ball–Berry stomatal conductance slope ($m_{\mathrm{bb}}$) at Wind River aligned with findings from recent CLM experiments at sites characterized by the same plant functional type (needleleaf evergreen temperate forest), despite significant differences in stand composition and age and climatology, suggesting that CLM could benefit from a revised $m_{\mathrm{bb}}$ value of 6, rather than the default value of 9, for this plant functional type. Conversely, Wind River required a unique calibration of the hydrology submodel to simulate soil moisture, suggesting that the default hydrology has a more limited applicability. This study demonstrates that carbon isotope data can be used to constrain stomatal conductance and intrinsic water use efficiency in CLM, as an alternative to eddy covariance flux measurements. It also demonstrates that carbon isotopes can expose structural weaknesses in the model and provide a key constraint that may guide future model development.

Introduction

The frequency, duration, and severity of droughts are expected to increase in the 21st century with climate change (Burke et al., 2006; Sheffield and Wood, 2008; Dai, 2013; Prein et al., 2016). In the western United States in particular, the combination of warmer temperature, larger vapor pressure deficit (VPD), reduced snowfall and snow pack, earlier snow melt, and extended growing season length is expected to lead to an intensification of water stress during the summer (Boisvenue and Running, 2010; Spies et al., 2010; Swain and Hayhoe, 2015; Fyfe et al., 2017). In this drying scenario, an accurate representation of ecosystem response to water stress in land models is critical for projecting carbon dynamics (and climate) into the future.

The land carbon and water cycles are coupled by the plant stomata through CO${}_{\mathrm{2}}$ uptake (photosynthesis) and water vapor loss (transpiration). While stomatal conductance responds to atmospheric vapor pressure deficit, soil moisture, and various other environmental factors, its modeling still represents a major challenge for the scientific community (Damour et al., 2010). Many stomatal conductance models have been proposed, including different approaches to account for water stress, but each model is subject to its own limitations (Damour et al., 2010; Miner et al., 2017; Sperry et al., 2017). Traditionally, stomatal conductance models have been calibrated through leaf- to canopy-level observations of water exchange.

Stable carbon isotopes provide an alternative observation to constrain stomatal conductance and offer an opportunity for model evaluation and improvement. During photosynthesis, plants discriminate against the heavier stable isotope of carbon (${}^{\mathrm{13}}$C) in favor of the lighter, more abundant ${}^{\mathrm{12}}$C stable isotope. This discrimination in C${}_{\mathrm{3}}$ plants, expressed as $\mathrm{\Delta }=\left[\right(R_{\mathrm{air}}/R_{\mathrm{plant}})-\mathrm{1}]\times \mathrm{1000}$ (‰), where $R_{\mathrm{air}}$ and $R_{\mathrm{plant}}$ are the ${}^{\mathrm{13}}$C : ${}^{\mathrm{12}}$C isotope ratios of atmospheric CO${}_{\mathrm{2}}$ and plant assimilated carbon, respectively, can be estimated according to the model proposed by Farquhar and Richards (1984) as $$\mathrm{\Delta }=a+\left(b,-,a\right)c_{\mathrm{i}}/c_{\mathrm{a}},$$ where $c_{\mathrm{i}}/c_{\mathrm{a}}$ is the ratio of intracellular CO${}_{\mathrm{2}}$ concentration to atmospheric CO${}_{\mathrm{2}}$ concentration, $a$ is the ${}^{\mathrm{13}}$C discrimination associated with the process of CO${}_{\mathrm{2}}$ diffusion through the stomata, and $b$ is the ${}^{\mathrm{13}}$C discrimination associated with the process of assimilation of CO${}_{\mathrm{2}}$ via rubisco ($a\approx \mathrm{4.4}$ ‰ and $b\approx \mathrm{27}$ ‰; Farquhar et al., 1989). The $c_{\mathrm{i}}/c_{\mathrm{a}}$ ratio correlates negatively with leaf intrinsic water use efficiency (iWUE), defined as the ratio of net leaf assimilation to stomatal conductance (Farquhar et al., 1989). Under water stress, C${}_{\mathrm{3}}$ plants tend to reduce stomatal conductance and increase water use efficiency, leading to reductions in $c_{\mathrm{i}}/c_{\mathrm{a}}$ and ${}^{\mathrm{13}}$C discrimination, affecting the carbon isotope ratio (${\mathit{\delta }}^{\mathrm{13}}$C) of photosynthesis and consequently of carbon pools and respiration. Experimental studies have shown, for instance, correlations between the ${\mathit{\delta }}^{\mathrm{13}}$C of ecosystem respiration and soil water content (SWC), atmospheric vapor pressure deficit, and precipitation. Bowling et al. (2008) and Brüggemann et al. (2011) present extensive reviews of experimental results on the link between environmental factors and the isotopic signature of carbon pools and fluxes, demonstrating that isotopic measurements provide insights into the response of stomatal conductance and iWUE to water stress. Furthermore, stable carbon isotopes have been used to partition photosynthetic and respiration fluxes from flux tower data (e.g., Wehr and Saleska, 2015) and to identify the strength of land and ocean sinks (e.g., Alden et al., 2010; van der Velde et al., 2013).

Photosynthetic ${}^{\mathrm{13}}$C discrimination is represented in biospheric models including the stable-isotope-enabled land surface model (ISOLSM) (Riley et al., 2002), the Simple Biosphere (SiB) model (Suits et al., 2005), the Lund–Potsdam–Jena (LPJ) dynamic global vegetation model (Scholze et al., 2003, 2008), the Land Surface Processes and Exchanges model of the University of Bern (LPX-Bern) (Spahni et al., 2013; Stocker et al., 2013), the hybrid SiB–CASA (combining biophysics from SiB and biochemistry from the Carnegie–Ames–Stanford Approach model) (van der Velde et al., 2013, 2014), and the Community Land Model (CLM) (Oleson et al., 2013).

Modeling studies have shown that stable carbon isotopes provide a constraint upon stomatal conductance (Aranibar et al., 2006; Raczka et al., 2016; Mao et al., 2016). Aranibar et al. (2006) evaluated the performance of ISOLSM at the Metolius Old Pine AmeriFlux site and were able to calibrate the slope of the stomatal conductance equation ($m_{\mathrm{bb}}$ in the Ball–Berry stomatal conductance model; see Eq. 2) with the aid of foliar ${\mathit{\delta }}^{\mathrm{13}}$C data measured at the site. Raczka et al. (2016) evaluated photosynthetic ${}^{\mathrm{13}}$C discrimination in CLM version 4.5 (CLM4.5) against ${\mathit{\delta }}^{\mathrm{13}}$C observations of photosynthesis and biomass at the Niwot Ridge AmeriFlux site and found the model to perform poorly with its default nitrogen limitation approach, resulting in overestimation of stomatal conductance and ${}^{\mathrm{13}}$C discrimination. By using an alternative approach in which a nitrogen downscaling factor is directly applied to $V_{\mathrm{cmax}\mathrm{25}}$ (maximum rate of carboxylation at 25 ${}^{\circ }$C), they found significant improvement in the simulations, but with results still suggesting that a smaller $m_{\mathrm{bb}}$ value (they used the default C${}_{\mathrm{3}}$ value, $m_{\mathrm{bb}}=\mathrm{9})$ would better simulate the site observations. Mao et al. (2016) evaluated CLM4.0 at a loblolly pine site in Tennessee, USA, and were able to adequately simulate the observed biomass ${\mathit{\delta }}^{\mathrm{13}}$C values with an optimized $m_{\mathrm{bb}}$ value of 5.6. Keller et al. (2017) used a global tree-ring ${\mathit{\delta }}^{\mathrm{13}}$C data set to evaluate the 20th century trend in photosynthetic ${}^{\mathrm{13}}$C discrimination and iWUE as modeled by CLM4.5 and LPX-Bern. LPX-Bern was found to perform well, while CLM simulated a significantly stronger increase (decrease) in iWUE (${}^{\mathrm{13}}$C discrimination) than that indicated by the tree-ring data. The default CLM parameterization and configuration were used in their study. Keller et al. (2017) suggested that the model–data mismatch was associated with the stomatal conductance parameterization ($m_{\mathrm{bb}}$ in particular) and the shortcomings of the nitrogen limitation scheme.

The present study focuses on CLM – the land component of the Community Earth System Model (CESM), a fully coupled global climate model widely used by the scientific community (http://www.cesm.ucar.edu/publications/) – and further evaluates the performance of its latest release (CLM4.5 – hereafter referred simply as CLM; Oleson et al., 2013; http://www.cesm.ucar.edu/models/cesm1.2/) against observations at a coniferous forest site in the Pacific Northwest region of the United States, with particular attention to the simulation of stomatal conductance and its response to water stress. The study site, Wind River, is part of the AmeriFlux eddy covariance network. Wind River is an old-growth forest ($\sim$ 500 years) characterized by a Mediterranean climate, due to which trees are naturally subject to water stress each summer. The combination of long-term measurements of energy and carbon fluxes, meteorology, biological variables, and stable carbon isotope ratios makes the site a good choice for evaluating carbon cycle and carbon isotope components of CLM. In addition to energy flux observations that allow for the estimation of canopy conductance, this study leverages the recent inclusion of photosynthetic ${}^{\mathrm{13}}$C discrimination within CLM and also uses ${\mathit{\delta }}^{\mathrm{13}}$C observations to better diagnose the simulation of stomatal conductance at the site. We test whether a reduced stomatal conductance at similar needleleaf evergreen temperate forest sites (Mao et al., 2016; Raczka et al., 2016) is appropriate for Wind River. This study also provides further investigation on the nitrogen limitation issue identified by Raczka et al. (2016) and the adequacy of the default parameters used in CLM, especially those regulating stomatal conductance. We test whether the calibration scheme (optimized parameters, nitrogen limitation) proposed by Raczka et al. (2016) for Niwot Ridge is appropriate for Wind River. By comparing the results at Wind River against those at different sites characterized by the same plant functional type (needleleaf evergreen temperate tree) but with different stand composition and age and climatology (Mao et al., 2016; Raczka et al., 2016), this study also seeks to identify general improvements in model parameterization.

Material and methods

This section provides a description of CLM, focusing on key formulations of relevance to the present study (Sect. 2.1), followed by a description of the study site (Sect. 2.2), the eddy covariance and meteorological data sets used to drive and assess the model (Sect. 2.3), the carbon isotope data sets used to assess the photosynthetic ${}^{\mathrm{13}}$C discrimination in CLM (Sect. 2.4), and also a description of the CLM configuration, simulations performed, and calibration of model parameters (Sect. 2.5 and 2.6). Section 2.7 describes the methodology used in the calculation of canopy conductance values from eddy covariance observations, which are compared against simulated values as a way to assess the model skill in simulating leaf stomatal conductance.

Model description

This section focuses on describing CLM's approach to the simulation of stomatal conductance and photosynthetic ${}^{\mathrm{13}}$C discrimination, key aspects of this study. For a full description of the model, the reader is referred to Oleson et al. (2013).

In CLM, leaf stomatal conductance ($g_{\mathrm{s}})$ is calculated based on the Ball–Berry model as described by Collatz et al. (1991) and implemented by Sellers et al. (1996) in the SiB2 model: $$g_{\mathrm{s}}=m_{\mathrm{bb}}{\displaystyle \frac{A_{\mathrm{n}}\left({\mathit{\beta }}_t\right)}{c_{\mathrm{s}}/P_{\mathrm{atm}}}}{h}_{\mathrm{s}}+b_{\mathrm{bb}}{\mathit{\beta }}_t,$$ where $A_{\mathrm{n}}\left({\mathit{\beta }}_t\right)$ is the potential net leaf photosynthesis (without nitrogen limitation) as a function of a soil moisture stress factor (${\mathit{\beta }}_t)$, $c_{\mathrm{s}}$ is the CO${}_{\mathrm{2}}$ partial pressure at the leaf surface, $P_{\mathrm{atm}}$ is the atmospheric pressure, $h_{\mathrm{s}}$ is the relative humidity at the leaf surface (defined as the ratio of vapor pressure at the leaf surface to saturation vapor pressure inside the leaf at vegetation temperature $T_{\mathrm{v}})$, $m_{\mathrm{bb}}$ is a slope coefficient, and $b_{\mathrm{bb}}$ corresponds to the minimum stomatal conductance in the original Ball–Berry model. The soil moisture stress factor ${\mathit{\beta }}_t$ is defined as $${\mathit{\beta }}_t={\sum }_i{w}_i{r}_i,$$ where $r_i$ is the root fraction at soil layer $i$ and $w_i$ is a corresponding plant wilting factor. The former is defined as (Oleson et al., 2013) $$r_i=\mathrm{0.5}\left(e^{-r_{\mathrm{a}}{z}_{h,i-\mathrm{1}}},+,e^{-r_{\mathrm{b}}{z}_{h,i-\mathrm{1}}}\right)-\mathrm{0.5}\mathit{\alpha }\left(e^{-r_{\mathrm{a}}{z}_{h,i}},+,e^{-r_{\mathrm{b}}{z}_{h,i}}\right),$$ where $z_{h,i}$ (m) is the depth from the soil surface to the interface between layers $i$ and $i+\mathrm{1}$ ($z_{h,\mathrm{0}}=\mathrm{0}$ corresponds to the soil surface), $r_{\mathrm{a}}$ and $r_{\mathrm{b}}$ are root distribution parameters (m${}^{-\mathrm{1}})$, $\mathit{\alpha }=\mathrm{1}$ for $\mathrm{1}\le i<N_{\mathrm{levsoi}}$, and $\mathit{\alpha }=\mathrm{0}$ for $i=N_{\mathrm{levsoi}}$ ($N_{\mathrm{levsoi}}$ is the number of soil layers). The plant wilting factor for soil layer $i$ is defined as (Oleson et al., 2013) $$w_i=\left\{\begin{array}{ll}{\displaystyle \frac{{\mathrm{\Psi }}_{\mathrm{c}}-{\mathrm{\Psi }}_i}{{\mathrm{\Psi }}_{\mathrm{c}}-{\mathrm{\Psi }}_{\mathrm{o}}}}\left[{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{sat},i}-{\mathit{\theta }}_{\mathrm{ice},i}}{{\mathit{\theta }}_{\mathrm{sat},i}}}\right]\le \mathrm{1}& \mathrm{for}{T}_i>T_{\mathrm{f}}-\mathrm{2}\\ & \mathrm{and}{\mathit{\theta }}_{\mathrm{liq},i}>\mathrm{0}\\ \mathrm{0}& \mathrm{for}{T}_i\le T_{\mathrm{f}}-\mathrm{2}\\ & \mathrm{or}{\mathit{\theta }}_{\mathrm{liq},i}=\mathrm{0}\end{array}\right.,$$ where ${\mathrm{\Psi }}_i$is the soil water matric potential, ${\mathrm{\Psi }}_{\mathrm{c}}$ and ${\mathrm{\Psi }}_{\mathrm{o}}$ are the soil water potential when stomata are fully closed or fully open, respectively (${\mathrm{\Psi }}_{\mathrm{c}}=-\mathrm{255}\mathrm{000}$ mm and ${\mathrm{\Psi }}_{\mathrm{o}}=-\mathrm{66}\mathrm{000}$ mm for the needleleaf evergreen temperate tree plant functional type, hereafter referred simply as NETT PFT), ${\mathit{\theta }}_{\mathrm{sat},i}$ is the saturated volumetric water content, ${\mathit{\theta }}_{\mathrm{ice},i}$ is the volumetric ice content, ${\mathit{\theta }}_{\mathrm{liq},i}$ is the volumetric liquid water content, $T_i$ is the soil layer temperature, and $T_{\mathrm{f}}=\mathrm{273.15}$ K is the freezing temperature of water. The sum in Eq. (3) is defined over the entire soil column, resulting in ${\mathit{\beta }}_t$ values from 0 (maximum soil moisture stress) to 1 (no soil moisture stress). In CLM's implementation of the Ball–Berry model (Eq. 2), ${\mathit{\beta }}_t$ is used to downscale $b_{\mathrm{bb}}$, directly impacting $g_{\mathrm{s}}$. ${\mathit{\beta }}_t$ also indirectly impacts $g_{\mathrm{s}}$ through the $A_{\mathrm{n}}$ term, as ${\mathit{\beta }}_t$ is used to downscale the maximum rate of carboxylation (${\mathit{\beta }}_t{V}_{\mathrm{cmax}})$ and also leaf respiration (${\mathit{\beta }}_t{R}_d)$ (Oleson et al., 2013).

Stomatal conductance ($g_{\mathrm{s}})$ and $A_{\mathrm{n}}$ are solved separately for sunlit and shaded leaves. Canopy conductance ($G_{\mathrm{c}})$ is given by $$G_{\mathrm{c}}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{b}}+r_{\mathrm{s}}^{\mathrm{sun}}}}{\mathrm{LAI}}^{\mathrm{sun}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{b}}+r_{\mathrm{s}}^{\mathrm{sha}}}}{\mathrm{LAI}}^{\mathrm{sha}}$$ and potential gross primary production (${\mathrm{GPP}}_{\mathrm{pot}}$, without nitrogen limitation) by $${\mathrm{GPP}}_{\mathrm{pot}}=(A_{\mathrm{n}}^{\mathrm{sun}}+R_d^{\mathrm{sun}}){\mathrm{LAI}}^{\mathrm{sun}}+(A_{\mathrm{n}}^{\mathrm{sha}}+R_d^{\mathrm{sha}}){\mathrm{LAI}}^{\mathrm{sha}},$$ where $r_{\mathrm{b}}$ is the leaf boundary layer resistance, $r_{\mathrm{s}}=\mathrm{1}/g_{\mathrm{s}}$ is the leaf stomatal resistance, $\mathrm{LAI}$ is the leaf area index, and $R_d$ is the leaf-level respiration (sun and sha superscripts denote sunlit and shaded leaves, respectively). Photosynthetic parameters such as $V_{\mathrm{cmax}\mathrm{25}}$ are solved separately for sunlit and shaded leaves and their canopy scaling scheme is detailed in Oleson et al. (2013, Sect. 8.3).

Based on nitrogen availability and nitrogen requirements for allocation of new carbon tissue, CLM calculates actual gross primary production ($\mathrm{GPP}$) as $$\mathrm{GPP}={\mathrm{GPP}}_{\mathrm{pot}}(\mathrm{1}-d).$$ The nitrogen down-regulation factor ($d)$ is defined as $$d={\displaystyle \frac{{\mathrm{CF}}_{\mathrm{avail}\mathrm{\_}\mathrm{alloc}}-{\mathrm{CF}}_{\mathrm{alloc}}}{{\mathrm{GPP}}_{\mathrm{pot}}}},$$ where ${\mathrm{CF}}_{\mathrm{avail}\mathrm{\_}\mathrm{alloc}}$ is the carbon flux from photosynthesis, which is available to new-growth allocation and ${\mathrm{CF}}_{\mathrm{alloc}}$ is the actual carbon allocation to new growth (limited by nitrogen availability). This implementation of nitrogen down-regulation makes CLM a partially coupled model with respect to net leaf photosynthesis and stomatal conductance. Note that GPP is calculated via down-regulation (Eq. 8) after the solution for $A_{\mathrm{n}}$ and $g_{\mathrm{s}}$ is obtained. Modeled $g_{\mathrm{s}}$ remains consistent with $A_{\mathrm{n}}$ (potential, not actual net leaf photosynthesis).

The original implementation of ${}^{\mathrm{13}}$C in CLM was developed in consultation with Neil Suits (Suits et al., 2005) and is described in Oleson et al. (2013). Photosynthetic ${}^{\mathrm{13}}$C discrimination in CLM for C${}_{\mathrm{3}}$ plants follows the model proposed by Farquhar and Richards (1984) (cf. Eq. 1): $$\mathrm{\Delta }=\mathrm{4.4}+\mathrm{22.6}{c}_{\mathrm{i}}/c_{\mathrm{a}}.$$ CLM calculates the intracellular-to-atmospheric CO${}_{\mathrm{2}}$ concentration ratio, $c_{\mathrm{i}}/c_{\mathrm{a}}$, in Eq. (10) as $${\displaystyle \frac{c_{\mathrm{i}}}{c_{\mathrm{a}}}}=\mathrm{1}-{\displaystyle \frac{A_{\mathrm{n}}(\mathrm{1}-d)}{c_{\mathrm{a}}}}\left[{\displaystyle \frac{\mathrm{1.4}}{g_{\mathrm{b}}}},+,{\displaystyle \frac{\mathrm{1.6}}{g_{\mathrm{s}}}}\right],$$ where $g_{\mathrm{b}}$ is the leaf boundary layer conductance. CLM does not account for mesophyll conductance (intracellular CO${}_{\mathrm{2}}$ is assumed to be the same as intercellular CO${}_{\mathrm{2}})$. Assuming $g_{\mathrm{b}}\gg g_{\mathrm{s}}$ (typically true for coniferous needles), Eq. (11) can be approximated by $$\begin{array}{cc}{\displaystyle \frac{c_{\mathrm{i}}}{c_{\mathrm{a}}}}& {\displaystyle }\cong \mathrm{1}-{\displaystyle \frac{\mathrm{1.6}(\mathrm{1}-d)}{c_{\mathrm{a}}}}\left[{\displaystyle \frac{A_{\mathrm{n}}}{g_{\mathrm{s}}}}\right]\\ & {\displaystyle }& {\displaystyle }\cong \mathrm{1}-{\displaystyle \frac{\mathrm{1.6}(\mathrm{1}-d)}{c_{\mathrm{a}}}}\mathrm{iWUE},\end{array}$$ where $\mathrm{iWUE}=A_{\mathrm{n}}/g_{\mathrm{s}}$ is the intrinsic water use efficiency. Note that $c_{\mathrm{i}}/c_{\mathrm{a}}$ and consequently $\mathrm{\Delta }$ correlate negatively with iWUE. All other terms being constant in Eq. (12), an increase in iWUE is expected to result in a reduction of the photosynthetic ${}^{\mathrm{13}}$C discrimination, i.e., an increase in the assimilation of the heavier ${}^{\mathrm{13}}$C stable isotope relative to the lighter, more abundant ${}^{\mathrm{12}}$C stable isotope. Note also that $A_{\mathrm{n}}$ is multiplied by $(\mathrm{1}-d)$ in Eqs. (11) and (12), making $c_{\mathrm{i}}$ consistent with the actual, nitrogen-limited $\mathrm{GPP}$. However, it is important to highlight that $g_{\mathrm{s}}$ is consistent with $A_{\mathrm{n}}$ (potential net assimilation), not $A_{\mathrm{n}}(\mathrm{1}-d)$ (actual net assimilation). The implications of this mismatch to the simulation of $\mathrm{\Delta }$ are discussed in Raczka et al. (2016) and later in the present paper.

Location and view of the Wind River AmeriFlux site, US-Wrc (satellite image from Google Earth).

[Figure omitted. See PDF]

The carbon isotope ratio of the $\mathrm{GPP}$ flux (${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}})$ is calculated in CLM based on the prescribed ${\mathit{\delta }}^{\mathrm{13}}$C of atmospheric CO${}_{\mathrm{2}}$, the carbon assimilation and photosynthetic ${}^{\mathrm{13}}$C discrimination by sunlit and shaded leaves, and their respective LAIs. The ${\mathit{\delta }}^{\mathrm{13}}$C of newly allocated carbon is the same as ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$. The ${\mathit{\delta }}^{\mathrm{13}}$C of the leaf carbon pool, for instance, depends on the allocation flux and its ${\mathit{\delta }}^{\mathrm{13}}$C (${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}})$ and the turnover time of the pool. CLM does not include any representation of post-photosynthetic ${}^{\mathrm{13}}$C discrimination.

Site description

The site for this study (Wind River) is part of the AmeriFlux eddy covariance network (Baldocchi et al., 2001), with a long record of meteorological, biological, surface flux (energy and carbon), and carbon isotope measurements for model assessment (1998–present). The site is located in the Pacific Northwest region of the United States, in the state of Washington (45.8205, $-$121.9519; 371 m elevation; see Fig. 1). Wind River is characterized by an old-growth conifer forest dominated by Douglas fir (Pseudotsuga menziesii) and western hemlock (Tsuga heterophylla) trees, with a mean canopy height of 56 m. Douglas fir trees are about 40–65 m high, corresponding to about 50 % of the wood volume of the stand and 33 % of the leaf area, while western hemlock trees are more numerous and smaller, corresponding to about 53 % of the leaf area of the stand (Unsworth et al., 2004; Parker et al., 2004). No significant disturbances have occurred at the site in the past $\sim$ 450–500 years. The local climate is strongly seasonal, marked by dry summers and wet winters. The climate summary reported by Shaw et al. (2004) indicates a mean annual precipitation of 2223 mm, with only $\approx$ 5 % falling during June, July, and August. During winter, much of the precipitation falls as snow, and the average snowpack depth exceeds 100 mm. The mean annual, January, and July air temperatures are $\mathrm{8.7}\pm \mathrm{6.5}$, $\mathrm{0.1}\pm \mathrm{2.3}$, and $\mathrm{17.7}\pm \mathrm{1.7}$ ${}^{\circ }$C, respectively.

Eddy covariance and meteorological data

Air temperature, relative humidity, wind speed, incident shortwave radiation, incident long-wave radiation, atmospheric pressure, and precipitation observed at the Wind River site from 1998 to 2006 were used to drive CLM. The time series were gap-filled using data from nearby towers and climate stations or interpolated in case of missing data. The gap-filled data product used to drive CLM in this study was created at Oak Ridge National Laboratory following the methodology described in Barr et al. (2013).

The L4 data set based on the eddy covariance observations (version V002, daily averages) was downloaded from the AmeriFlux repository at ftp://cdiac.ornl.gov/pub/ameriflux/data/ (AmeriFlux data are currently hosted at http://ameriflux.lbl.gov, see Wharton, 1998–2016). This data set contains friction-velocity-filtered, gap-filled, and partitioned fluxes and was used to assess the simulated surface fluxes of sensible heat ($H$), latent heat (LE), and carbon, including GPP and ecosystem respiration (ER). The ER product was estimated according to the short-term temperature response of measured nighttime net ecosystem exchange (NEE) (Reichstein et al., 2005), and GPP was estimated as the difference between ER and NEE: i.e., ER – NEE. The gap-filled NEE values (and derived GPP and ER) using the marginal distribution sampling method (Reichstein et al., 2005) were used in this study.

Eddy covariance and meteorological data from the AmeriFlux L2 data product (version V007, 30 min averages) were used to calculate canopy conductance ($G_{\mathrm{c}}$; see Sect. 2.7) and atmospheric VPD. In the analysis, 30 min surface flux data were rejected during periods when the wind direction was in the [45${}^{\circ }$ : 135${}^{\circ }$] sector (same criterion used by Wharton et al., 2012), as the northeast-to-southeast wind sector is characterized by heterogeneous (age-fragmented) land cover. The data were averaged hourly prior to $G_{\mathrm{c}}$ and VPD calculation. L2 SWC data were also used in the analysis. Missing SWC data from the L2 data set in the year 2002 were replaced by respective L1 data (version Apr2013).

The AmeriFlux L2 data product (version V007, 30 min averages) was also used to assess the energy balance closure at the site. The energy balance ratio, $\mathrm{EBR}=(H+\mathrm{LE})/(R_{\mathrm{n}}-G)$, where $R_{\mathrm{n}}$ is net radiation and $G$ is soil heat flux, was calculated for dry season months (June to September) using 10:00–14:00 PST data and rejecting periods with rain or unfavorable wind direction ([45${}^{\circ }$ : 135${}^{\circ }$] sector). With the available data, $\mathrm{EBR}$ could be calculated for the years of 1998–2001, 2004, and 2006.

Carbon isotope data

Estimated ${\mathit{\delta }}^{\mathrm{13}}$C values of ER (Lai et al., 2005) and observed ${\mathit{\delta }}^{\mathrm{13}}$C values of leaf tissue and soil organic matter (SOM) (Fessenden and Ehleringer, 2003) at Wind River were used to assess the photosynthetic ${}^{\mathrm{13}}$C discrimination in CLM. Lai et al. (2005) used an automated air sampling system, with inlets at 0.5 m above ground level and at 0.5 canopy height, collecting 15 flasks weekly during the growing season. Most of the flasks (13 out of 15) were dedicated to nighttime sampling (over a single night). The Keeling plot method was used to infer the weekly ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ using the CO${}_{\mathrm{2}}$ and ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ observations (for simplicity, the resulting ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ values are referred to as observations in the text). The monthly averages (June–November) from 2001 to 2003 reported by Lai et al. (2005) were used as references in the present study. Fessenden and Ehleringer (2003) conducted measurements of ${\mathit{\delta }}^{\mathrm{13}}$C of bulk organic tissue from current-year needles of Tsuga heterophylla trees and seedlings at the top (55 m), middle (25 m), and bottom (2 m) of the canopy. They also conducted vertical profile measurements of ${\mathit{\delta }}^{\mathrm{13}}$C of bulk soil organic carbon down to 20 cm depth. The measurements were performed on a 1-month to 2-month time interval. The values reported by Fessenden and Ehleringer (2003) for the growing season in 1999 and 2000 were used as references.

In the present study, both observed and modeled carbon isotope ratios were expressed as ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=(\frac{R_x}{R_{\mathrm{std}}}-\mathrm{1})\times \mathrm{1000}$ (‰), where $R_x$ is the ${}^{\mathrm{13}}$C : ${}^{\mathrm{12}}$C isotope ratio of the carbon pool or flux of interest and $R_{\mathrm{std}}$ is the ${}^{\mathrm{13}}$C : ${}^{\mathrm{12}}$C isotope ratio of a standard reference material (Vienna Pee Dee Belemnite standard).

CLM configuration and simulations

CLM was run at site level using the Point CLM (PTCLM) scripting framework (see Kluzek, 2013), as in recent studies (e.g., Mao et al., 2016; Raczka et al., 2016). Land cover was defined as the needleleaf evergreen temperate tree plant functional type. The model was configured to use CLM4.5 physics and CLM4.5 carbon–nitrogen (CN) biogeochemistry. The vertical soil carbon profile option was turned on, and the CENTURY carbon model was selected for the decomposition parameters. The nitrification and denitrification sub-model was switched off, as preliminary simulations indicated an excess of nitrogen availability and forest productivity when the respective module was active. Given that the Wind River site is characterized by an old-growth mature forest, no land-cover disturbance was considered in the simulations.

The model was spun-up in a two-stage process, using a preindustrial component set with a constant, preindustrial atmospheric CO${}_{\mathrm{2}}$ concentration and ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ of 285 ppmv and $-$6.5 ‰, respectively. The model was run in accelerated decomposition mode for 600 years (first stage) and then in normal decomposition mode for 1000 years (second and final stage), using the local observed meteorological data (Sect. 2.3) from 1998 to 2006 to drive the model (continuously cycled). Following the spin-up process, a transient run (1850–2006) was performed with prescribed nitrogen deposition, atmospheric CO${}_{\mathrm{2}}$ concentration, and atmospheric ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$.

The transient atmospheric CO${}_{\mathrm{2}}$ concentrations used in this study were based on the CMIP5 recommendations for annual global mean values (Meinshausen et al., 2011). The transient atmospheric ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ values used here were based on ice-core and flask measurements reported by Francey et al. (1999) (annual values in their spline fitting from 1850 to 1981) and flask measurements in Mauna Loa (annual averages from 1981 to 2006) by the Scripps CO${}_{\mathrm{2}}$ program (Keeling et al., 2005), following a similar methodology as in Raczka et al. (2016) (note that, unlike in Raczka et al., 2016, here a seasonal cycle was not superimposed onto the time series). As in the spin-up process, the local observed meteorological data from 1998 to 2006 were cycled during the transient run. The driver data and model years were aligned in a way to guarantee a perfect match between them during the final 9 years of the simulation (1998–2006).

Summary of changes in CLM parameters during the calibration process. The parameters listed, excluding $M$, $Q_{\mathrm{10}}$, $m_{\mathrm{bb}}$, $b_{\mathrm{bb}}$, and soil hydraulic parameters, correspond to the needleleaf evergreen temperate tree plant functional type (NETT PFT).

Parameter	Description	CLM name	Default CLM	Calibrated CLM
			value	value
$A_{\text{s : l}}$	New stem C: new leaf C ratio (g C g C${}^{-\mathrm{1}})$	stem_leaf	Dynamic	2
$A_{\text{fr : l}}$	New fine root C: new leaf C ratio (g C g C${}^{-\mathrm{1}})$	froot_leaf	1	0.385
CN${}_{\mathrm{l}}$	Leaf C : N ratio (g C g N${}^{-\mathrm{1}})$	leafcn	35	38.2
CN${}_{\mathrm{llit}}$	Leaf litter C : N ratio (g C g N${}^{-\mathrm{1}})$	lflitcn	70	76.4
CN${}_{\mathrm{fr}}$	Fine root C : N ratio (g C g N${}^{-\mathrm{1}})$	frootcn	42	64.7
${\mathit{\tau }}_{\mathrm{l}}$	Leaf longevity (yr)	leaf_long	3	5
$r_{\mathrm{b}}$	Root distribution parameter (m${}^{-\mathrm{1}})$	rootb_par	2	1
${\mathrm{SLA}}_{\mathrm{0}}$	Specific leaf area at canopy top (m${}^{\mathrm{2}}$ leaf g C${}^{-\mathrm{1}})$	slatop	0.010	0.006
$m$	$\mathrm{SLA}\left(x\right)$ slope (m${}^{\mathrm{2}}$ ground g C${}^{-\mathrm{1}})$	dsladlai	0.00125	0.000985
$M$	Plant mortality rate (% yr${}^{-\mathrm{1}})$	am	2	1.5
$Q_{\mathrm{10}}$	Temperature sensitivity coefficient of maintenance	q10	1.5	2.5
	respiration and decomposition (–)
Soil hydraulic	Version used	origflag (namelist	0 (CLM4.5)	1 (CLM4.0)
parameters		variable)
$m_{\mathrm{bb}}$	Ball–Berry equation slope (–)	mbbopt	9	6
$b_{\mathrm{bb}}$	Ball–Berry equation intercept ($\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ leaf s${}^{-\mathrm{1}})$	bbbopt	10 000	5000

CLM calibration

Initial simulations using the default parameters from CLM resulted in a poor representation of the carbon dynamics at the Wind River site (Figs. A5 and A6). GPP and forest biomass were significantly underestimated. The seasonality of ER was poorly represented and the simulated late-summer GPP was impacted by an underestimation of SWC, resulting from a poor representation of soil hydrology at the site. Furthermore, the modeled evapotranspiration values were significantly overestimated for the given values of GPP, i.e., the simulated water use efficiency was much lower than that observed. As a result of CLM's poor performance in the simulation of GPP and evapotranspiration, the modeled photosynthetic ${}^{\mathrm{13}}$C discrimination was found to be overestimated, reflected in reduced ${}^{\mathrm{13}}$C : ${}^{\mathrm{12}}$C ratios of carbon fluxes and pools.

In order to improve the representation of the carbon dynamics at the site, key model parameters were calibrated as detailed in Appendix A and summarized in Table 1. The adjusted parameters were primarily based on biological measurements at Wind River or at similar stands in the Pacific Northwest. Parameters controlling specific leaf area and stomatal conductance were found to be critical to the simulation of GPP and evapotranspiration and were manually adjusted in a way to minimize the differences between model output and site observations (eddy covariance fluxes). The default soil hydraulic parameters used in CLM4.5 were found to be inadequate at Wind River, leading to severe underestimation of SWC and unrealistic soil moisture stress in the model during late summer. These parameters were reverted back to their default values in CLM4.0, with significant improvement in the representation of soil hydrology at the site. In an additional measure to reduce the unrealistic late-summer soil moisture stress in the model, root distribution was adjusted based on CLM's default parameter values for the broadleaf evergreen temperate tree plant functional type, shifting roots towards deeper soil layers (justified based on physical understanding of the site – see Appendix A8).

Bayesian parameter calibration is a common approach used in modeling studies to account for both the prior parameter distributions and more recent observations. In this case, a Bayesian calibration approach would be complicated by the current lack of prior parameter distributions within CLM in order to create a model ensemble and the computational expense of running a calibration. Commonly used techniques such as Markov chain Monte Carlo (MCMC) are prohibitively expensive with long CLM simulations, and more advanced techniques for calibration (e.g., using surrogate modeling approaches) are still under development. The simpler approach used here proved to be an effective method to improve model performance at the Wind River AmeriFlux site. The reader is referred to Appendix A for a more complete description of the parameters that were adjusted and the calibration approach used. All model results presented and discussed in Sects. 3 and 4, unless noted otherwise, are based on the optimized model.

Canopy conductance

Observed canopy conductance ($G_{\mathrm{c}}$, m s${}^{-\mathrm{1}})$ was calculated by combining hourly tower data (see Sect. 2.3) with the Penmann–Monteith equation (Monteith, 1964) as in Wharton et al. (2012): $$G_{\mathrm{c}}={\left[{\displaystyle \frac{\mathit{\rho }{c}_p\mathrm{VPD}}{\mathit{\gamma }\mathrm{LE}}},+,{\displaystyle \frac{\left(\frac{{\mathrm{\Delta }}_{\mathrm{sat}}}{\mathit{\gamma }}\right)\left(\frac{H}{\mathrm{LE}}\right)-\mathrm{1}}{G_{\mathrm{a}}}}\right]}^{-\mathrm{1}},$$ where $\mathit{\rho }$ and $c_p$ are the density and specific heat of air, respectively (kg m${}^{-\mathrm{3}}$, J kg${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}})$, VPD is the atmospheric vapor pressure deficit (kPa), LE is the latent heat flux (W m${}^{-\mathrm{2}})$, ${\mathrm{\Delta }}_{\mathrm{sat}}$ is the slope of the saturation vapor pressure curve as a function of air temperature (kPa K${}^{-\mathrm{1}})$, $\mathit{\gamma }$ is the psychrometric constant (kPa K${}^{-\mathrm{1}})$, $H$ is the sensible heat flux (W m${}^{-\mathrm{2}})$, and $G_{\mathrm{a}}=u_{\ast }^{\mathrm{2}}/U$ is the aerodynamic conductance for momentum transfer (m s${}^{-\mathrm{1}})$, where $u_{\ast }$ is the friction velocity and $U$ is the wind speed. Atmospheric pressure and air temperature data and the ideal gas law were later used to convert the $G_{\mathrm{c}}$ values to mmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. The calculation of $G_{\mathrm{c}}$ was restricted to daytime hours (10:00–16:00 PST) and to the months of June to September (dry season). Rain events and periods with LE $<$ 5 W m${}^{-\mathrm{2}}$ or relative humidity $>$ 80 % were disregarded. $G_{\mathrm{c}}$ values outside the interval of 0 to 1000 mmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ were also disregarded.

Modeled leaf area index (a), carbon stocks (b–f), and ${\mathit{\delta }}^{\mathrm{13}}$C of leaf (g) and soil organic matter (h) during the transient run (lines) compared against site observations (points and error bars). Modeled values in panels (a–f) correspond to annual averages. Modeled ${\mathit{\delta }}^{\mathrm{13}}$C values in panels (g) and (h) were calculated from annual averages of the respective ${}^{\mathrm{13}}$C and ${}^{\mathrm{12}}$C pools. Observations in panels (a–e) (average $\pm$ SD) are from the AmeriFlux database (based on Thomas and Winner, 2000 and Harmon et al., 2004). Observations in panels (g) and (h) correspond to the average $\pm$ SD of the measurements reported by Fessenden and Ehleringer (2003) in their Figs. 2b and 3 (leaf ${\mathit{\delta }}^{\mathrm{13}}$C at canopy top (55 m), middle (25 m), and bottom (2 m) and SOM ${\mathit{\delta }}^{\mathrm{13}}$C at 20 cm depth).

[Figure omitted. See PDF]

For comparison against observations, modeled canopy conductance values were calculated using the same methodology described above, but using hourly CLM output ($H$, LE, $u_{\ast })$ instead. An alternative would be to calculate canopy conductance directly by upscaling CLM's leaf stomatal conductance and leaf boundary layer conductance using LAI (Eq. 6). Canopy conductance values derived from both approaches were found to be strongly correlated. The Penmann–Monteith method was ultimately selected for the calculation of $G_{\mathrm{c}}$ in order to allow a more direct comparison between modeled and observed values. This comparison was done as a way to assess the performance of CLM in the simulation of leaf stomatal conductance.

Results

Carbon pools and isotopic signatures

Figure 2 shows modeled LAI, carbon stocks (leaf, fine root, coarse root, tree wood, and SOM carbon), and ${\mathit{\delta }}^{\mathrm{13}}$C of leaf and SOM pools throughout the transient run (1850–2006). Before the transient run, the model was spun-up and successfully equilibrated under the defined preindustrial scenario, with LAI, carbon stocks, and leaf and SOM carbon isotope ratios reaching steady state (results not shown). The cyclic behavior exhibited in Fig. 2 is related to the driving meteorological data set, which was cycled throughout the simulation period (Sect. 2.5).

From 1850 to 2006, modeled LAI and carbon stocks (Fig. 2a–f) increased due to CO${}_{\mathrm{2}}$ fertilization and increasing nitrogen deposition. Average values of LAI, leaf carbon, and tree wood carbon were in agreement with the reference values reported in the AmeriFlux database for the Wind River site. Modeled fine root and coarse root carbon were underestimated, but within 2 standard deviations from the reference values.

The ${\mathit{\delta }}^{\mathrm{13}}$C of leaves and SOM was initialized in the model with a value of $-$6 ‰ (default value in CLM, close to the preindustrial atmospheric ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ value of $-$6.5 ‰ used in this study). During the model spin-up, in which constant preindustrial atmospheric ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ and CO${}_{\mathrm{2}}$ concentration values were prescribed, the ${\mathit{\delta }}^{\mathrm{13}}$C values stabilized at $\approx -\mathrm{26}$ ‰. During the transient simulation, the ${\mathit{\delta }}^{\mathrm{13}}$C of both leaves (Fig. 2g) and SOM (Fig. 2h) decreased (the pools became isotopically “lighter”), mostly due to the decreasing atmospheric ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ values associated with the Suess effect (Keeling, 1979) but also due to the increasing atmospheric CO${}_{\mathrm{2}}$ concentration values. The ${\mathit{\delta }}^{\mathrm{13}}$C of leaves declined faster over the years than the ${\mathit{\delta }}^{\mathrm{13}}$C of SOM, given the fact that leaves have a significantly shorter turnover time than SOM and therefore present a faster response to the changes in atmospheric ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ and CO${}_{\mathrm{2}}$ concentration. Modeled ${\mathit{\delta }}^{\mathrm{13}}$C of leaves compared well against the site observations for top- and mid-canopy leaves ($-$0.8 and $+$0.8 ‰ difference, respectively), and modeled ${\mathit{\delta }}^{\mathrm{13}}$C of SOM (top 1 m of soil) compared well against site observations for SOM at 20 cm below ground ($-$0.4 ‰ difference).

It is important to clarify that CLM has leaf properties that vary continuously with canopy depth and that two leaf categories (sunlit and shaded leaves) are estimated dynamically at every time step, as a function of canopy structure and solar elevation angle (Thornton and Zimmerman, 2007). The modeled leaf ${\mathit{\delta }}^{\mathrm{13}}$C output corresponds to the isotopic signature of the entire leaf carbon pool, which is calculated from both sunlit and shaded portions of the leaf canopy (see Sect. 2.1). The observed leaf ${\mathit{\delta }}^{\mathrm{13}}$C values in Fig. 2g correspond to measurements at canopy top (55 m), middle (25 m), and bottom (2 m). As pointed out by Fessenden and Ehleringer (2003), the decrease in the observed leaf ${\mathit{\delta }}^{\mathrm{13}}$C values (i.e., increase in photosynthetic ${}^{\mathrm{13}}$C discrimination) with canopy depth can be explained by light reduction within the canopy. In principle, the observed mid-canopy values are expected to better represent the isotopic composition of leaves for the whole canopy, in comparison with the observed values at the two canopy extremes, especially given the larger amount of leaf biomass in the mid-canopy. However, considering how light is reduced within the canopy, the top-canopy ${\mathit{\delta }}^{\mathrm{13}}$C value should still be representative of a significant fraction of the canopy as well; thus, the whole canopy ${\mathit{\delta }}^{\mathrm{13}}$C is expected to lay somewhere in between the top- and mid-canopy values. As shown in Fig. 2g, the modeled ${\mathit{\delta }}^{\mathrm{13}}$C of the leaf carbon pool was the average between the observed values at canopy top and middle.

The overall agreement between the observed and modeled carbon isotope ratios indicates that CLM had skill in simulating the balance between assimilation and stomatal conductance and the associated photosynthetic ${}^{\mathrm{13}}$C discrimination. The adjustment of the parameters controlling stomatal conductance in the model ($m_{\mathrm{bb}}$ and $b_{\mathrm{bb}}$ – see Sect. 2.6, Table 1 and Appendix A9) to improve the simulation of evapotranspiration had a significant impact on the simulation of photosynthetic ${}^{\mathrm{13}}$C discrimination. When using the default parameter values (resulting in significantly higher stomatal conductance values), the modeled values of ${\mathit{\delta }}^{\mathrm{13}}$C in leaves and SOM were generally 2–3 ‰ lower (Fig. A1), departing from site observations.

Energy and carbon fluxes

Modeled energy and carbon fluxes are compared against daily-averaged observations in Fig. 3 for the period between 1998 and 2006. “Observed” GPP and ER values were obtained from applying a partitioning model to NEE measurements (Sect. 2.3), but are referred to as observations in the text.

Modeled LE values were close to observations, with a mean bias error (MBE) of $\approx$ $-$3 W m${}^{-\mathrm{2}}$ and a RMSE of $\approx$ 20 W m${}^{-\mathrm{2}}$. The adjustment of the stomatal conductance parameters $m_{\mathrm{bb}}$ and $b_{\mathrm{bb}}$ (Table 1) was fundamental in modifying the LE simulations. When using the default parameter values, the modeled evapotranspiration was significantly overestimated, with summer values exceeding observations by almost 100 % (Fig. A2a).

In 1998–2003 the model overestimated $H$ (MBE $\approx$ 32 W m${}^{-\mathrm{2}}$, RMSE $\approx$ 40 W m${}^{-\mathrm{2}})$, while in 2004–2006 the modeled values were closer to observations (MBE $\approx$ 10 W m${}^{-\mathrm{2}}$, RMSE $\approx$ 36 W m${}^{-\mathrm{2}})$. The modeled $H$ values did not present significant interannual variability in 1998–2006; however, the observations showed significantly smaller fluxes in 1998–2003 than in 2004–2006. Such changes in the magnitude of $H$ were reported as a potential data issue in the Wind River site documentation available in the AmeriFlux repository.

Modeled sensible heat flux (a), latent heat flux (b), gross primary production (c), and ecosystem respiration (d) vs. site observations. Orange and red (cyan and blue) lines correspond to modeled (observed) values. For a clearer visualization, the daily averages (thin lines) were smoothed with a Bézier algorithm (thick lines).

[Figure omitted. See PDF]

The overall mean $\mathrm{EBR}$ calculated from site observations was 0.88 (see Sect. 2.3 for calculation approach). The energy balance closure for years 2004 and 2006 was high (mean EBR $=$ 1.01 and 1.09, respectively). The model bias of $H$ and LE was relatively small in those years (Fig. 3a, b). In the years 1998, 2000, and 2001, mean $\mathrm{EBR}$ was significantly lower (0.63, 0.69, and 0.76, respectively). Modeled $H$ presented a relatively large positive bias in those years (Fig. 3a). As discussed above, the observed $H$ values in 1998–2003 were significantly smaller than in 2004–2006, while the LE observations showed approximately the same pattern over the years. The low $\mathrm{EBR}$ for the years 1998, 2000, and 2001 supports the remark included in the AmeriFlux documentation regarding a potential data issue with $H$ and suggests that the observed values were biased low in 1998–2003. Mean $\mathrm{EBR}$ in 1999 was relatively high (0.92), in which the reduced $H$ values (Fig. 3a) were compensated for by larger LE values (Fig. 3b). In that year, modeled $H$ (LE) had a positive (negative) bias in respect to the observations.

Modeled GPP resembled observed values, with small differences (MBE $\approx$ 0.23 g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$, RMSE $\approx$ 1.60 g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}})$. Modeled ER exhibited closer correspondence with measurements during the spring and summer months in general (MBE $\approx$ 0.82 g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$, RMSE $\approx$ 1.85 g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}})$, with summer peaks especially close to measured values. In the colder months, modeled ER was significantly overestimated (MBE $\approx$ 1.46 g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$, RMSE $\approx$ 1.77 g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}})$.

Despite the significant improvement in the seasonal behavior of ER after the $Q_{\mathrm{10}}$ adjustments discussed in Appendix A6, the results indicate that further adjustments also including the base rate of maintenance respiration and the base decomposition rates for each litter and SOM pool within CLM would be necessary to better simulate the observed ER at Wind River. The results suggest that lower base rates and higher $Q_{\mathrm{10}}$ values would improve the simulations at the site.

Isotopic signatures of GPP and ER

Diurnal cycle

Modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ exhibited a well-defined diurnal cycle (Fig. 4), with minimum values in the early morning and late afternoon and a peak value typically in midafternoon, reflecting diurnal changes in the simulated iWUE (see Eqs. 10 and 12). Modeled ${\mathit{\delta }}^{\mathrm{13}}$C values of the heterotrophic component of ecosystem respiration (HR) were approximately constant, with a $\approx \mathrm{0.2}$ ‰ change over the entire period of study (1998–2006). Conversely, modeled ${\mathit{\delta }}^{\mathrm{13}}$C values of the autotrophic component (AR) were found to be virtually equal to modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ values during daytime. At nighttime, modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{AR}}$ was found to change abruptly towards values closer to modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{HR}}$. Because AR was the major component of the total ecosystem respiration (ER $=$ AR $+$ HR; see Fig. 4a), modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ exhibited a similar behavior compared to modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{AR}}$ (Fig. 4b).

Mean diurnal cycle of modeled carbon fluxes (a) and their respective carbon isotope ratios (b) for the summer months (June–September) of the years 1998–2006. Fluxes include gross primary production (black circles), ecosystem respiration (purple squares), autotrophic respiration (blue crosses), and heterotrophic respiration (red diamonds). Bars correspond to $\pm$1 SD.

[Figure omitted. See PDF]

In CLM, newly assimilated carbon is first allocated to meet the total maintenance respiration demand of live plant tissues (top priority). When this demand exceeds the supply of carbon via photosynthesis (e.g., during nocturnal periods, wintertime, stress periods), carbon is drawn from a storage pool (excess maintenance respiration pool; CS${}_{\mathrm{xs}})$, which is allowed to run a deficit state. The reason CLM allows this deficit state is to avoid the requirement of knowing the size of the total storage pool available to plants and thus the possibility of vegetation dying in a given location if the storage pool is depleted. When negative, CS${}_{\mathrm{xs}}$ is gradually replenished with newly assimilated carbon at a potential rate of $-{\mathrm{CS}}_{\mathrm{xs}}/{\mathit{\tau }}_{\mathrm{xs}}$, where ${\mathit{\tau }}_{\mathrm{xs}}$ is a time constant (set to 30 days in CLM). The carbon allocation flux to replenish CS${}_{\mathrm{xs}}$ receives second priority in the model, while the carbon allocation fluxes to support plant growth have third priority. Given this allocation structure, ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{AR}}$ will follow ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ during daytime (assuming GPP is enough to meet the maintenance respiration demand) and the ${\mathit{\delta }}^{\mathrm{13}}$C of the “excess maintenance respiration flux” (XSMR) during nighttime. CLM does not calculate the isotopic signature of XSMR from CS${}_{\mathrm{xs}}$, but from bulk vegetation tissues (total vegetation carbon, TOTVEGC). This is done because CS${}_{\mathrm{xs}}$ is not a physical quantity but a construct of CLM. Note that XSMR borrows carbon from the CS${}_{\mathrm{xs}}$ pool, which is allowed to run a deficit state. This debt is paid in the future with the replenishment of the CS${}_{\mathrm{xs}}$ pool with newly assimilated carbon. This construct makes the ${\mathit{\delta }}^{\mathrm{13}}$C of CS${}_{\mathrm{xs}}$ nonphysical; therefore, the approximation that ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{XSMR}}={\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{TOTVEGC}}$ is more physically realistic. This approximation makes the nocturnal ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{AR}}$ follow ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{TOTVEGC}}$, explaining the low sensitivity of the nocturnal ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{AR}}$ to recent ${}^{\mathrm{13}}$C discrimination in the results shown in Fig. 4b.

Autotrophic respiration at Wind River is likely fueled by a mixture of stored and recently fixed carbon, as indicated by ${}^{\mathrm{14}}$C measurements from root respiration (RR) at the site (Taylor et al., 2015). This process cannot be appropriately modeled by CLM with the current carbon allocation scheme, impacting the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$. An explicit representation of carbohydrate storage pools within CLM to support the maintenance respiration demand would improve the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$. The need for a better representation of carbohydrate storage pools within CLM was also highlighted by the ${}^{\mathrm{13}}$CO${}_{\mathrm{2}}$-labeling study conducted by Mao et al. (2016).

It is important to highlight that, unlike models such as SiB (Sellers et al., 1996; Vidale and Stöckli, 2005), CLM does not have a prognostic canopy airspace where ${\mathit{\delta }}^{\mathrm{13}}$CO${}_{\mathrm{2}}$ is impacted by photosynthetic and respiratory fluxes; thus, the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ is not affected by the limitations in the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ described above.

Seasonal cycle

Modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ exhibited a well-defined seasonal pattern, peaking during the summer as a result of a decrease in the photosynthetic ${}^{\mathrm{13}}$C discrimination associated with higher iWUE values (Fig. 5; see also Eqs. 10 and 12). The summer peak in iWUE was linked to changes in stomatal conductance in response to increased VPD and reduced SWC during the dry summer season.

On a monthly scale, roughly indicated by the smoothed curve in Fig. 5, the modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ values presented a similar seasonal pattern in comparison with the ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ observations by Lai et al. (2005) at the site. Differences between ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ and ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ are obviously expected, as ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ depends on the contribution of recently assimilated carbon to AR, the AR : ER ratio, and also post-photosynthetic fractionation (Bowling et al., 2008; Brüggemann et al., 2011). The seasonal pattern in the observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ (Fig. 5) could be partially attributed to an eventual spring-to-summer decrease in the AR : ER ratio (assuming ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{HR}}$ > ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{AR}})$. ${}^{\mathrm{14}}$C measurements from belowground respiration components at Wind River reported by Taylor et al. (2015) do indicate a spring-to-summer decrease in the contribution of RR towards total soil respiration (SR $=$ RR $+$ HR). The similarity of the seasonal patterns of observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ and modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ suggests that stomatal response to water stress could also be driving the seasonal pattern in the observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ at the site. The broader implication is that ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$, which can be more easily measured than ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$, can be reasonably used as a surrogate to indicate forest response to water stress at Wind River.

Modeled ${\mathit{\delta }}^{\mathrm{13}}$C of gross primary production (lines) and observed ${\mathit{\delta }}^{\mathrm{13}}$C of ecosystem respiration (circles). Thin orange line corresponds to daily averages using 10:00–16:00 PST data only. For a clearer visualization, this curve was smoothed with a Bézier algorithm (thick red line). Blue circles correspond to site observations (monthly averages) reported by Lai et al. (2005).

[Figure omitted. See PDF]

Due to the limitations in the carbon allocation scheme used in CLM (Sect. 3.3.1), the simulated ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ values were found to be inconsistent with the site observations, with nocturnal values approximately constant throughout the entire period of study (1998–2006), exhibiting little sensitivity to recent photosynthetic ${}^{\mathrm{13}}$C discrimination. Diurnal values, however, were found to be strongly correlated with ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$, given the fact that in CLM current photosynthate directly fuels AR (results not shown).

The adjustment of the stomatal conductance parameters $m_{\mathrm{bb}}$ and $b_{\mathrm{bb}}$ to improve the simulation of evapotranspiration (Sect. 2.6, Table 1 and Appendix A9) led to a significant change in the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$. When the default parameter values were used, modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ was generally 2–3 ‰ lower due to higher photosynthetic ${}^{\mathrm{13}}$C discrimination (Fig. A2b), also presenting a considerable reduction in the amplitude of the seasonal cycle. The difference between modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ and observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ was significantly larger. As discussed in Sect. 3.1, site observations of leaf and SOM ${\mathit{\delta }}^{\mathrm{13}}$C support the notion that the default stomatal conductance parameters are inadequate at Wind River, resulting in excessive photosynthetic ${}^{\mathrm{13}}$C discrimination.

Ecosystem response to water stress

Overall, CLM was able to reasonably capture the observed interannual variability in GPP at the study site (Fig. 3c). The behavior of observed GPP in 2002 and 2006 stands out, showing an early-season peak followed by a quick reduction, suggesting strong water stress in those years, especially in 2002. Among the years studied here, 2002 and 2006 had the lowest summer precipitation. Spring precipitation was also low in 2006 but normal in 2002. Observed canopy conductance during the spring and summer of 2006 was smaller than in 2002, but a stronger attenuation of GPP was observed in 2002, suggesting that water stress was not the main reason for the attenuated GPP values in 2002. CLM was able to simulate the observed GPP behavior in 2006 but not in 2002. The reason for the model–data mismatch in the spring–summer of 2002 is currently unclear. Despite the fact that meteorological forcing data from 1998 to 2006 were continuously cycled throughout the transient run (1850–2006), meaning that the impact of any slow secular change in the forcing data was not captured in the simulation, the simulated GPP still compared reasonably well against observations.

Throughout the simulation period (1998–2006), CLM predicted a few periods when the ecosystem was under the influence of soil moisture stress (Fig. 6). As indicated by the ${\mathit{\beta }}_t$ parameter (Eq. 3), which varies from 0 (maximum soil moisture stress) to 1 (no soil moisture stress) (see Sect. 2.1), those periods included the summers of 1998, 2006, 2003, and 2002, in decreasing order of stress intensity. The departures from ${\mathit{\beta }}_t=\mathrm{1}$ typically occurred when modeled SWC (top five soil layers, 0–27 cm) decreased below $\approx \mathrm{20}$ %. Note that, at Wind River, SWC at permanent wilting point and at field capacity is 14 and 30 %, respectively (Wharton et al., 2009).

Hourly soil water content and CLM's soil moisture stress parameter, ${\mathit{\beta }}_t$ (black line). Observed SWC (blue line) corresponds to the integrated value for the top 30 cm of soil, while modeled SWC (red line) corresponds to the integrated value for the top five soil layers in CLM (0–27 cm). At Wind River, SWC at permanent wilting point and at field capacity is 14 and 30 %, respectively (Wharton et al., 2009).

[Figure omitted. See PDF]

With the adjustment of soil hydraulic parameters (Appendix A7), CLM was able to adequately simulate SWC throughout most of the years within the study period (Fig. 6), especially during the summer months, with an overall summer MBE of 3.24 %. However, the simulated SWC significantly departed from observations in 1999–2002. CLM, which was driven by observed precipitation at the site, indicated higher SWC than observations in 1999–2002, particularly during the summer months, with a summer MBE of 8.05 %. For the remaining years, summer MBE was $-$0.27 %. The SWC observations starting on the second year of the site records (1999) up to the data gap in 2002 presented a different pattern in comparison with the remaining years, showing an apparent negative offset of near 10 %. It is likely that the apparent shift in the time series of observed SWC was instrument-related. In 1999–2002, soil moisture monitoring at the site consisted of two, two-pronged time-domain reflectometer (TDR) probes instead of six, three-pronged TDR probes, likely resulting in less-accurate data collection.

Hourly observed canopy conductance vs. observed VPD for the summer months (June–September) of the years 1999–2006, restricted to 10:00–16:00 PST (additional restrictions were imposed to the calculation of $G_{\mathrm{c}}$; see Sect. 2.7). The years 1998 and 2005 were not included due to missing data. Data points were segregated according to observed SWC in panel (a) and according to modeled SWC in panel (b) (see Fig. 6). Lines correspond to the linear regression between $\log G_{\mathrm{c}}$ and VPD using all data points (solid lines) and using only points within the lowest SWC bin (red circles, dashed lines).

[Figure omitted. See PDF]

Observed canopy conductance was found to be strongly dependent on VPD, following a decreasing exponential relationship (Fig. 7). In order to investigate the additional dependence on soil moisture stress, the data points were divided into four bins according to the observed values of SWC (Fig. 7a). The linear regression fit between $\log G_{\mathrm{c}}$ and VPD for the points corresponding to the lowest SWC bin (SWC < 17.5 %, $\approx \mathrm{22}$ % of all data points) was virtually the same as the linear regression considering all data points. If the forest were under soil moisture stress at those low SWC levels, the former regression curve with data points from the lowest SWC bin would be expected to be found below the latter. Instead, the SWC < 17.5 % regression curve was very similar – even slightly above the regression curve using all data points.

The lack of sensitivity of observed $G_{\mathrm{c}}$ to observed SWC (Fig. 7a) was likely associated with a negative bias in SWC in 1999–2002. Observed $G_{\mathrm{c}}$ was found to respond to modeled SWC (driven by observed precipitation) (Fig. 7b). As discussed above, the observed SWC values in 1999–2002 were suspected to have a negative bias, i.e., drier than reality. $G_{\mathrm{c}}$ values in Fig. 7a corresponding to the summer of 1999–2002 were tagged as belonging to the lowest SWC bin, but in reality, they could be associated with wetter, non-moisture stress conditions. Assuming that CLM's summer simulated SWC (driven by observed precipitation) was not as biased as the observed SWC might be, we instead used the modeled SWC values to probe the $G_{\mathrm{c}}$ vs. VPD relationship under different SWC regimes in Fig. 7b. Interestingly, with this approach, a distinct pattern emerged for the data points within the lowest SWC bin. The regression curve considering all data points was $\log G_{\mathrm{c}}=-\mathrm{0.59}\mathrm{VPD}+\mathrm{6.06}$ ($r=-\mathrm{0.60})$ and when considering only the data points from the lowest bin (modeled SWC < 21.25 %, $\approx \mathrm{24}$ % of all points), the regression curve was $\log G_{\mathrm{c}}=-\mathrm{0.50}\mathrm{VPD}+\mathrm{5.71}$ ($r=-\mathrm{0.56})$. The latter regression curve corresponded to reasonably lower $G_{\mathrm{c}}$ values, especially at low VPD levels, which is compatible with a moisture stress scenario. The result supports the suspicion of a negative bias in the observed SWC data in 1999–2002.

Similar to observations, modeled canopy conductance was also found to be strongly dependent on VPD (Fig. 8). This is expected given the Ball–Berry stomatal conductance model used in CLM (Eq. 2). The Ball–Berry model has a direct dependence on leaf relative humidity (leaf RH), not leaf VPD, but these variables are strongly correlated. The correlation between modeled $G_{\mathrm{c}}$ and RH was found to be slightly higher than between modeled $G_{\mathrm{c}}$ and VPD, while observed $G_{\mathrm{c}}$ correlated slightly better with VPD than RH (results not shown). The results indicate that a direct dependence on leaf VPD, rather than leaf RH, in CLM's stomatal conductance model would lead to a more accurate representation of stomatal functioning at Wind River, but overall, for the period analyzed in the present study, such improvement is expected to be small. The general dependence of modeled canopy conductance on VPD was very similar in comparison with observations, as indicated by the linear regression curve between $\log G_{\mathrm{c}}$ and VPD in Fig. 8 using all data points ($\log G_{\mathrm{c}}=-\mathrm{0.59}\mathrm{VPD}+\mathrm{6.04}$; compare with $\log G_{\mathrm{c}}=-\mathrm{0.59}\mathrm{VPD}+\mathrm{6.06}$ in Fig. 7b). The correlation between observed $\log G_{\mathrm{c}}$ and VPD, however, was lower than for the model results ($r=-\mathrm{0.60}$ and $r=-\mathrm{0.91}$, respectively).

Hourly modeled canopy conductance vs. observed VPD for the summer months (June–September) of the years 1999–2006, restricted to 10:00–16:00 PST (additional restrictions were imposed to the calculation of $G_{\mathrm{c}}$; see Sect. 2.7). Note that observed air temperature and relative humidity were used to drive CLM. The years 1998 and 2005 were not included for consistency with Fig. 7. Data points were segregated according to the soil moisture stress parameter ${\mathit{\beta }}_t$ in panel (a) and according to modeled SWC in panel (b) (see Fig. 6). Lines correspond to the linear regression between $\log G_{\mathrm{c}}$ and VPD using all data points (solid lines) and using only points within the lowest SWC bin (red circles, dashed line).

[Figure omitted. See PDF]

The impact of soil moisture stress on $G_{\mathrm{c}}$ was reasonably captured in CLM (Fig. 8b; cf. Fig. 7b). The impact of soil moisture stress on modeled $G_{\mathrm{c}}$ is clearly visible in Fig. 8a, in which the data points were binned according to ${\mathit{\beta }}_t$. With increasing soil moisture stress (decreasing ${\mathit{\beta }}_t$ values), the modeled $G_{\mathrm{c}}$ values still maintained a strong dependence on VPD, but were shifted downward, particularly at low VPD levels. In order to allow a more direct comparison against Fig. 7b, the data points were binned according to modeled SWC in Fig. 8b. The points in the lowest SWC bin (SWC < 21.25 %, $\approx \mathrm{22}$ % of all points) roughly corresponded to the periods under soil moisture stress (${\mathit{\beta }}_t<\mathrm{1})$. The regression curve for the SWC < 21.25 % group lay reasonably below the regression curve considering all data points ($\log G_{\mathrm{c}}=-\mathrm{0.53}\mathrm{VPD}+\mathrm{5.80}$, $r=-\mathrm{0.90}$, and $\log G_{\mathrm{c}}=-\mathrm{0.59}\mathrm{VPD}+\mathrm{6.04}$, $r=-\mathrm{0.91}$, respectively). The regression curves associated with SWC < 21.25 % were similar for the observed and modeled results (Figs. 7b and 8b), indicating that CLM could reasonably simulate soil moisture stress at Wind River, although with a small underestimation (i.e., a small overestimation of $G_{\mathrm{c}}$; note the $G_{\mathrm{c}}$ intercepts at 301 and 331 mmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ in Figs. 7b and 8b, respectively). It is important to point out, however, that modeled SWC was used to segregate the observations in Fig. 7b due to the potential bias in the SWC observations discussed above.

Modeled ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ and $G_{\mathrm{c}}$ values were highly correlated ($r=-\mathrm{0.88}$, $p<\mathrm{0.001}$; Fig. 9b). Modeled $G_{\mathrm{c}}$ generally decreased into the summer season, leading to an increase in water use efficiency and a decrease in photosynthetic ${}^{\mathrm{13}}$C discrimination, resulting in higher ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ values. Observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ was found to have a low negative correlation with observed $G_{\mathrm{c}}$, but it was not statistically significant ($r=-\mathrm{0.27}$, $p=\mathrm{0.396}$; Fig. 9a). The low correlation was likely a result of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ reflecting constraints of prior environmental drivers in comparison with the more rapid response of $G_{\mathrm{c}}$ to more recent environmental drivers. Another possible explanation is that the monthly ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ values in Fig. 9a were obtained by averaging up to four discrete weekly observations (see Sect. 2.4), in contrast with the calculation of monthly $G_{\mathrm{c}}$, which used daytime values for each day of the month. It is important to mention that the observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ values show a clear seasonal pattern (Fig. 5), with values peaking during summer likely in response to changes in $g_{\mathrm{s}}$ and iWUE associated with increasing water stress (see discussion in Sect. 3.3.2), but the present results indicate a lag in this response.

Observed ${\mathit{\delta }}^{\mathrm{13}}$C of ecosystem respiration vs. observed canopy conductance (a) and modeled ${\mathit{\delta }}^{\mathrm{13}}$C of gross primary production vs. modeled canopy conductance (b) for the summer months of 2001–2003. Except for the observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$, data points correspond to monthly averages of daytime (10:00–16:00 PST) data (additional restrictions were imposed to the calculation of $G_{\mathrm{c}}$; see Sect. 2.7). Observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ corresponds to the monthly averages reported by Lai et al. (2005). Numbers at the center of each point indicate the month.

[Figure omitted. See PDF]

Discussion

Ecosystem response to water stress

We found that the major cause of water stress leading to stomatal response at Wind River during summer was the elevated VPD, and not the reduced soil moisture (Sect. 3.4). Observed canopy conductance values at the site strongly decreased at moderate VPD levels, regardless of soil moisture conditions (Fig. 7b). The high sensitivity of stomatal response to changes in VPD was also shown and discussed in Wharton et al. (2009). As pointed out in their study, “even under moderate VPD levels, foliage at the tops of tall evergreen conifer trees often reach near critical values for cavitation due to a long path distance between the water table and the hydraulic capacity of the xylem, and as a result shut their stomata frequently (Ryan and Yoder 1997)”. They also point out that soil moisture depletion is usually not limiting at the site because the mature trees are capable of tapping water from deeper soil layers. This is generally consistent with our findings (Sect. 3.4); however, we also found that stomatal conductance responded to soil moisture stress during periods of more severe SWC depletion and low VPD (Fig. 7b).

Overall, CLM was able to simulate the observed response of canopy conductance to VPD and SWC, reasonably capturing the impact of water stress on ecosystem functioning (Fig. 8). Similarly to observations, VPD exerted a strong limitation on modeled $G_{\mathrm{c}}$, while SWC was usually not limiting. Note that ${\mathit{\beta }}_t$ was equal to 1 (no soil moisture stress) throughout most of the period of study (Fig. 6), in alignment with the explanation by Wharton et al. (2009) that the mature trees at the site are capable of accessing water from deeper soil layers. Note also that the default NETT PFT root distribution in CLM was shifted towards deeper soil layers (Appendix A8), aiming to improve the simulation of ${\mathit{\beta }}_t$. Despite the good overall model–data agreement ($G_{\mathrm{c}}$ dependency on VPD and SWC) after calibration, the results indicate a small underestimation of soil moisture stress in CLM, as discussed in Sect. 3.4. Calibration of the parameters controlling the plant wilting factor (Eq. 5) and additional calibration of the root distribution parameters could improve the results but are out of scope here.

An obvious but important point that must be highlighted is that in order to adequately simulate soil moisture stress, CLM must first adequately simulate SWC. Even when the model is driven by observed precipitation data (the case of the present study), this task is not trivial. As discussed in Appendix A7, CLM's hydrology sub-model performed poorly at Wind River when the default soil hydraulic parameters were used, leading to a strong dry bias in SWC. The original parameters used in the previous version of CLM (version 4.0) were found to perform much better at the site, likely due to a reduction of subsurface runoff and consequent increase in water retention in the soil column. As the default parameter values are intended for global simulations, it is natural to expect site-to-site variation in model performance (see Sect. 4.2). Raczka et al. (2016), for instance, did not find issues with the default soil hydraulic parameters in their CLM4.5 simulation at the Niwot Ridge AmeriFlux site. This difference in impact between the sites may have resulted from unique soil properties or differences in precipitation and evaporative demand between the sites during the summer.

As pointed out in Sect. 3.4, the results of the present study indicate that a direct dependence on leaf VPD in CLM's stomatal conductance model, rather than leaf RH, would lead to a more accurate representation of stomatal functioning at Wind River, but overall, for the period analyzed in the present study, such improvement is expected to be small. It is important to emphasize that this expectation refers to the results presented here only. In the case of model predictions under future climate scenarios, in which atmospheric VPD is predicted to change while RH stays the same (as discussed in Sato et al., 2015), a direct dependence on leaf VPD in the stomatal conductance model becomes critical. The next CLM release (version 5) is expected to replace the Ball–Berry model with the Medlyn model (Medlyn et al., 2011), which directly depends on leaf VPD. This modification is expected to be more relevant for climate change simulations. Note that the present analysis is based on a hindcast simulation using a stable climate.

Calibration of CLM

Substantial calibration of model parameters was necessary to simulate the observed energy and carbon dynamics at Wind River, an old-growth ($\sim$ 500 years old) coniferous forest site dominated by Douglas fir and western hemlock trees and characterized by a Mediterranean climate. This is not surprising given that the default parameters used in CLM are intended for global simulations. Thus, model performance at particular sites is expected to vary greatly, requiring site-specific calibration in order to adequately simulate the observations. This is also demonstrated in the studies by Raczka et al. (2016) and Mao et al. (2016). Raczka et al. (2016) investigated the performance of CLM at the Niwot Ridge AmeriFlux site, a $\sim$ 110-year-old subalpine coniferous forest site in Colorado, USA, consisting of lodgepole pine, Engelmann spruce, and subalpine fir, while Mao et al. (2016) evaluated CLM in a 10-year-old loblolly pine stand in Tennessee, USA. In both cases significant site-specific specification and calibration of model parameters were also necessary. Note that these sites fall into the same PFT category as Wind River (NETT PFT). Despite the significant differences between the three sites, the results presented here and in Raczka et al. (2016) and Mao et al. (2016) converge in respect to the calibration of the Ball–Berry stomatal conductance slope, $m_{\mathrm{bb}}$. It is promising that despite the range in stand age and climate conditions amongst these sites, there appears to be a consensus that reduced stomatal conductance is required across all sites. This bodes well when upscaling to regional simulations.

A reduction of $m_{\mathrm{bb}}$ from 9 (default) to 6 was necessary to simulate the observed GPP, LE, and ${\mathit{\delta }}^{\mathrm{13}}$C values (leaf, SOM) at Wind River. This aligns with the results by Mao et al. (2016), as they were able to simulate the observations at their Tennessee site, including biomass ${\mathit{\delta }}^{\mathrm{13}}$C values, with an optimized $m_{\mathrm{bb}}$ of 5.6. However, as discussed in Appendix A9, the present results show that the significant reduction of $m_{\mathrm{bb}}$ from 9 to 6 may represent a tradeoff with model representation of nitrogen limitation. When using CLM's default nitrogen limitation scheme and $m_{\mathrm{bb}}$ value, Raczka et al. (2016) found significant overestimation of ${}^{\mathrm{13}}$C discrimination at Niwot Ridge due to excessive stomatal conductance, similar to the present study. When using an alternative nitrogen limitation scheme based on $V_{\mathrm{cmax}\mathrm{25}}$ down-regulation, maintaining the coupling between net leaf assimilation and $g_{\mathrm{s}}$, Raczka et al. (2016) found significant improvement in the simulations. This alternative scheme was also tested here while keeping the default $m_{\mathrm{bb}}$ value, and the results were similar compared to the model run with the default nitrogen limitation scheme and $m_{\mathrm{bb}}=\mathrm{6}$ (Appendix A9).

The results in the present study indicate that it is possible to account for the partial coupling between net leaf assimilation and stomatal conductance in CLM through the adjustment of $m_{\mathrm{bb}}$ to achieve reasonable carbon and energy exchange behavior, including ${}^{\mathrm{13}}$C discrimination. This is also supported by the results in Mao et al. (2016). A more detailed evaluation of model skill in simulating ${}^{\mathrm{13}}$C discrimination with this approach, in comparison with the $V_{\mathrm{cmax}\mathrm{25}}$ down-regulation approach (fully coupled CLM), would depend on high-frequency observations of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ as in Raczka et al. (2016). These data were not available at Wind River. Note that ${}^{\mathrm{13}}$C discrimination at Wind River was inferred from ${\mathit{\delta }}^{\mathrm{13}}$C measurements of leaves and SOM.

The results in Raczka et al. (2016), Mao et al. (2016), and the present study indicate that $m_{\mathrm{bb}}=\mathrm{9}$ is excessive when the default nitrogen limitation implementation is used in the simulations, with the latter two studies indicating that $m_{\mathrm{bb}}\approx \mathrm{6}$ is a more appropriate value to simulate the site observations. This agreement at three very distinct sites is encouraging and suggests that CLM could possibly benefit from a revised $m_{\mathrm{bb}}$ value of 6 for the NETT PFT, keeping in mind that such adjustment to improve model skill would also account for structural error. At the same time, the results presented here and in Raczka et al. (2016) indicate that the default $m_{\mathrm{bb}}=\mathrm{9}$ is reasonable for simulations when the $V_{\mathrm{cmax}\mathrm{25}}$ down-regulation scheme is implemented in the model, although Raczka et al. (2016) still found a small overestimation of ${}^{\mathrm{13}}$C discrimination at Niwot Ridge, suggesting that a smaller $m_{\mathrm{bb}}$ value would better simulate the site dynamics. It is important to point out, however, that the experimental literature indicates generally lower $m_{\mathrm{bb}}$ values for coniferous forests (see for example the surveys by Williams et al., 2004, Table 6.3, and Miner et al., 2017, Fig. 1). The SiB model (Sellers et al., 1996), for instance, uses $m_{\mathrm{bb}}=\mathrm{6}$ for conifers and $m_{\mathrm{bb}}=\mathrm{9}$ for other C${}_{\mathrm{3}}$ plants, while CLM uses $m_{\mathrm{bb}}=\mathrm{9}$ for all C${}_{\mathrm{3}}$ plants. Further investigation of the applicability of the revised $m_{\mathrm{bb}}$ value (or the current default value while using the $V_{\mathrm{cmax}\mathrm{25}}$ down-regulation scheme as in Raczka et al., 2016) at other NETT PFT sites is recommended for future studies.

Recommendations for structural improvement within CLM

The results of the present study demonstrate that ${\mathit{\delta }}^{\mathrm{13}}$C observations can be used to constrain stomatal conductance and iWUE in CLM as an alternative to eddy covariance flux measurements, leveraging the recent implementation of photosynthetic ${}^{\mathrm{13}}$C discrimination within the model. The adjustments made on the parameters controlling stomatal conductance within the model, originally aiming to improve the simulation of evapotranspiration, were critical to simulating the observed photosynthetic ${}^{\mathrm{13}}$C discrimination at Wind River, inferred from ${\mathit{\delta }}^{\mathrm{13}}$C measurements of leaves and SOM. As discussed in Sect. 4.2, these adjustments to improve model skill interacted strongly with the nitrogen limitation scheme. A possible interpretation of results from this and other recent studies is that growth limitation due to restricted nitrogen availability does not operate instantaneously upon photosynthesis (e.g., through nitrogen downscaling in the default version of CLM4.5) but is accounted for further “downstream” during the allocation of carbon.

For example, Metcalfe et al. (2017) proposed a revised model structure in which GPP is not instantaneously down-regulated during photosynthesis, but the excess photosynthate, which cannot be allocated to structural pools due to insufficient nitrogen supply, is allocated to a new nonstructural carbohydrate storage pool within the model. Carbon from this pool is able to return to the atmosphere via the inclusion of a single additional respiration term within the model. This new model structure provides a solution for the issue regarding the partial coupling between net leaf assimilation and stomatal conductance. Alternatively, a foliar nitrogen model could be used to account for nitrogen limitation directly within the estimation of photosynthetic capacity (Ghimire et al., 2016), removing the requirement for nitrogen downscaling. A similar approach is planned to be included in the next release of CLM (version 5.0).

The use of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ observations as a strong constraint upon CLM is hindered by the lack of an explicit representation of carbohydrate storage pools within the model to support autotrophic respiration (Fig. 4). The results from the ${}^{\mathrm{13}}$C-labeling study by Mao et al. (2016) also illustrate the issue and highlight the need for structural improvements in CLM's carbon allocation scheme. One implication of this issue is that it prevents a more direct use of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ observations – which are easier to obtain and more frequently available than ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}}$ observations – for evaluation of ${}^{\mathrm{13}}$C discrimination in CLM. It may also limit the applicability of CLM for global atmospheric ${}^{\mathrm{13}}$C budget studies focusing on land–ocean flux partitioning (e.g., van der Velde et al., 2013), as errors in the simulation of the land isotopic disequilibrium (${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}-{\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{GPP}})$ can propagate to the estimation of the land–ocean partitioning and the estimation of variability in each sink (van der Velde et al., 2014). Van der Velde et al. (2014) were able to reasonably simulate mean observed ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ values for a selection of sites from the Biosphere–Atmosphere Stable Isotope Network (BASIN; Pataki et al., 2003) using a modified version of the SiB–CASA model including representation of ${}^{\mathrm{13}}$C isotopes and modified carbon storage pools. The original SiB–CASA model (Schaefer et al., 2008) has a single storage pool representing sugars and starch, with only the sugar portion being readily available for plant growth and maintenance. The effective pool turnover rate in this configuration is $\sim$ 70 days. In the modified model, sugar and starch allocation are simulated separately with two distinct pools, with prescribed turnover rates of 7 days (sugar to starch) and 63 days (starch to sugar). Van der Velde et al. (2014) found significant improvement in the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ with the new carbon allocation approach. We recommend that CLM adopt a similar carbon allocation scheme, moving away from the deficit-based accounting scheme (Sect. 3.3.1) towards an explicit representation of carbohydrate storage pools such as in the SiB–CASA model (van der Velde et al., 2014).

Another shortcoming in CLM is the fact that mesophyll conductance ($g_{\mathrm{m}})$ is not simulated, i.e., intracellular and intercellular CO${}_{\mathrm{2}}$ values are assumed to be equal. As demonstrated here and in Raczka et al. (2016) and Mao et al. (2016), CLM is able to reasonably simulate ${}^{\mathrm{13}}$C discrimination by either adjusting the stomatal conductance slope parameter or using an alternative nitrogen limitation scheme ($V_{\mathrm{cmax}\mathrm{25}}$ down-regulation), but the impact of not including $g_{\mathrm{m}}$ in the simulations must be investigated. Mesophyll conductance was recently incorporated in CLM by Sun et al. (2014); however, it still has to be linked to the carbon isotope submodel. This is another area in which ${}^{\mathrm{13}}$C observations can be used for model evaluation and development.

Conclusions

After substantial calibration of model parameters, CLM was able to simulate energy and carbon fluxes, leaf area index, and carbon stocks at an old-growth coniferous forest (Wind River AmeriFlux site) in general agreement with site observations. Overall, the calibrated CLM was able to simulate the observed response of canopy conductance to atmospheric vapor pressure deficit and soil water content, reasonably capturing the impact of water stress on ecosystem functioning. Key model adjustments to simulating observed flux and carbon stock patterns included (1) parameters controlling the variation in specific leaf area through the forest canopy (SLA${}_{\mathrm{0}}$, $m$), with significant impact on GPP, (2) parameters controlling stomatal conductance ($m_{\mathrm{bb}}$, $b_{\mathrm{bb}})$, with significant impact on the simulated latent heat flux and water use efficiency, and (3) soil hydraulic parameters, with impact on soil water content and on the soil moisture stress parameter, ${\mathit{\beta }}_t$.

The calibrated CLM was able to simulate carbon isotope ratios of leaves and soil organic matter at Wind River, in general agreement with site observations. The adjustments made on the parameters controlling stomatal conductance within the model, originally aiming to improve the simulation of evapotranspiration, were critical to simulating the observed photosynthetic ${}^{\mathrm{13}}$C discrimination at the site, inferred from ${\mathit{\delta }}^{\mathrm{13}}$C measurements of leaves and soil organic matter. This demonstrates that stable carbon isotopes can serve as an alternative to eddy covariance flux measurements for constraining stomatal conductance. The simulation of nocturnal ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ was found to be inconsistent with site observations, with results showing little sensitivity to recent photosynthetic ${}^{\mathrm{13}}$C discrimination. The inclusion of explicit carbohydrate storage pools within CLM (and removal of the current deficit-based carbon accounting system) to support the maintenance respiration demand from live plant tissues would improve the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$.

We found that an optimized stomatal slope value ($m_{\mathrm{bb}}=\mathrm{6})$ was necessary at Wind River, consistent with previous CLM experiments from distinct needleleaf evergreen temperate forest sites. This suggests that this parameterization could apply to broader-scale simulations of this PFT. We also found a tradeoff between adjustment of stomatal slope and changes to the nitrogen limitation scheme. The best long-term solution may be to replace this nitrogen scheme with alternative approaches.

The hydrology submodel within CLM and its parameterization deserve special attention because the simulation of soil water content has a direct impact on ${\mathit{\beta }}_t$, and thus on stomatal conductance. Wind River required a unique calibration to achieve reasonable soil moisture, which was not consistent across other sites. This suggests that simulation of soil moisture in regional studies should be used with caution.

The recent inclusion of the photosynthetic ${}^{\mathrm{13}}$C discrimination functionality in CLM opens a new opportunity for model testing and development. The results presented here demonstrate that carbon isotopes can expose structural weaknesses in the model, such as the deficit-based accounting system in CLM's carbon allocation scheme and the partial coupling between net leaf photosynthesis and stomatal conductance caused by the nitrogen limitation scheme. ${\mathit{\delta }}^{\mathrm{13}}$C observations provide a key constraint that may guide future CLM development.

CLM4.5 (Oleson et al., 2013) can be obtained at http://www.cesm.ucar.edu/models/cesm1.2/. AmeriFlux data are currently available at http://ameriflux.lbl.gov (see Wharton, 1998–2016). Model output data presented in this paper are available upon request to the corresponding author.

CLM calibration

Most of the adjustments were performed on parameters particular to the needleleaf evergreen temperate tree plant functional type in CLM. For brevity, this plant functional type is referred to as NETT PFT in the following sections.

Carbon allocation ratios

By default, CLM uses a dynamic new-stem-carbon-to-new-leaf-carbon allocation ratio ($A_{\text{s : l}}$, g C g C${}^{-\mathrm{1}})$ for the NETT PFT, which rises with increasing net primary production. A survey by White et al. (2000) indicates an average $A_{\text{s : l}}$ of $\mathrm{2.2}\pm \mathrm{0.89}$ g C g C${}^{-\mathrm{1}}$ for needleleaf evergreen forests. Measurements reported by Hudiburg et al. (2013) for a region close to the Wind River site and characterized by forests of similar species composition vary between approximately 1 and 3.5 g C g C${}^{-\mathrm{1}}$ (their Fig. A1 – Mesic sites). A fixed value of $A_{\text{s : l}}=\mathrm{2}$ g C g C${}^{-\mathrm{1}}$ (value also used by Thornton et al., 2002 in their Biome-BGC simulations for the Wind River site) was found to improve the simulated forest biomass and was adopted in this study for the NETT PFT.

The new-fine-root-carbon-to-new-leaf-carbon allocation ratio parameter ($A_{\text{fr : l}}$, g C g C${}^{-\mathrm{1}})$ for the NETT PFT was also changed based on observations at the Wind River site reported in the AmeriFlux database, indicating $A_{\text{fr : l}}=\mathrm{0.385}$ g C g C${}^{-\mathrm{1}}$ rather than the default value of 1 g C g C${}^{-\mathrm{1}}$. The change meant a significantly greater carbon investment to leaves, helping to increase the modeled GPP towards the site observations.

Carbon : nitrogen ratios

Leaf litter C : N ratio (CN${}_{\mathrm{llit}}$, g C g N${}^{-\mathrm{1}})$ for the NETT PFT was adjusted based on measurements at the Wind River site (Klopatek, 2007) to 76.4 g C g N${}^{-\mathrm{1}}$ (mean observed value). Based on the mean observed CN${}_{\mathrm{llit}}$ and assuming a nitrogen retranslocation efficiency of 50 % (survey by Parkinson, 1983, indicates efficiencies around 50 % for conifer trees and 36–69 % for Douglas fir in particular), the leaf C : N ratio (CN${}_{\mathrm{l}}$, g C g N${}^{-\mathrm{1}})$ for NETT PFT was adjusted to 38.2 g C g N${}^{-\mathrm{1}}$. The updated parameters differ little from the default values (CN${}_{\mathrm{llit}}=\mathrm{70}$ g C g N${}^{-\mathrm{1}}$, CN${}_{\mathrm{l}}=\mathrm{35}$ g C g N${}^{-\mathrm{1}})$.

Fine-root C : N ratio (CN${}_{\mathrm{fr}}$, g C g N${}^{-\mathrm{1}})$ for the NETT PFT was also adjusted based on measurements at the Wind River site (Klopatek, 2007). The value was adjusted from 42 g C g N${}^{-\mathrm{1}}$ (default) to 64.7 g C g N${}^{-\mathrm{1}}$ (mean observed value), meaning a significantly smaller nitrogen investment in fine roots resulting in more nitrogen for investment in leaves. This change helped to increase the modeled GPP towards the site observations.

Modeled ${\mathit{\delta }}^{\mathrm{13}}$C of leaf (a) and soil organic matter (b), calculated from annual averages of the respective ${}^{\mathrm{13}}$C and ${}^{\mathrm{12}}$C pools during the transient run (lines). Results from two model configurations are presented: CLM (calibrated model, solid lines) and CLM${}^{\ast }$ (calibrated model using the default stomatal conductance parameters ($m_{\mathrm{bb}}$ and $b_{\mathrm{bb}}$; see Table 1), dashed lines). Site observations (average $\pm$ SD, blue points and error bars) are also shown (see caption of Fig. 2 for details).

[Figure omitted. See PDF]

Leaf longevity

Measurements reported by Hudiburg et al. (2013) for a region near the Wind River site and characterized by forests of similar species composition indicate leaf longevity (${\mathit{\tau }}_{\mathrm{l}})$ of 5 years. This value was adopted for the NETT PFT, replacing the default value of 3 years. This change contributed particularly to an increase in the modeled leaf area index.

Specific leaf area

In CLM, specific leaf area (SLA, m${}^{\mathrm{2}}$ leaf g C${}^{-\mathrm{1}})$ is assumed to be linear with canopy depth $x$ (expressed as overlying leaf area index, m${}^{\mathrm{2}}$ leaf m${}^{-\mathrm{2}}$ ground) (Thornton and Zimmermann, 2007): $$\mathrm{SLA}\left(x\right)=\mathrm{SL}{A}_{\mathrm{0}}+mx,$$ where SLA${}_{\mathrm{0}}$ is the specific leaf area at the top of canopy and $m$ is a linear coefficient (m${}^{\mathrm{2}}$ ground g C${}^{-\mathrm{1}})$. Integrating this equation over the canopy, a relationship can be established in which leaf area index (LAI, m${}^{\mathrm{2}}$ leaf m${}^{-\mathrm{2}}$ ground) is calculated as a function of leaf carbon (C${}_{\mathrm{l}}$, g C m${}^{-\mathrm{2}}$ ground), knowing the parameters SLA${}_{\mathrm{0}}$ and $m$ (Thornton and Zimmermann, 2007): $$\mathrm{LAI}={\displaystyle \frac{\mathrm{SL}{A}_{\mathrm{0}}\left(e^{m{\mathrm{C}}_{\mathrm{l}}},-,\mathrm{1}\right)}{m}}.$$ The default NETT PFT values in CLM for SLA${}_{\mathrm{0}}$ and $m$ are 0.01 m${}^{\mathrm{2}}$ leaf g C${}^{-\mathrm{1}}$ and 0.00125 m${}^{\mathrm{2}}$ ground g C${}^{-\mathrm{1}}$, respectively. These values were found to be too large for the Wind River site. Using them in Eq. (A2) with a C${}_{\mathrm{l}}$ of 941 g C m${}^{-\mathrm{2}}$ ground (mean observation at the Wind River site reported in the AmeriFlux database) results in an LAI of $\approx \mathrm{18}$ m${}^{\mathrm{2}}$ leaf m${}^{-\mathrm{2}}$ ground, instead of $\approx \mathrm{9}$ m${}^{\mathrm{2}}$ leaf m${}^{-\mathrm{2}}$ ground according to the observations at the Wind River site (AmeriFlux database).

In CLM, the maximum rate of carboxylation at 25 ${}^{\circ }$C ($V_{\mathrm{cmax}\mathrm{25}})$ is proportional to the area-based leaf nitrogen concentration defined as ${\mathrm{N}}_{\mathrm{a}}=\mathrm{1}/\left({\mathrm{CN}}_{\mathrm{l}}{\mathrm{SLA}}_{\mathrm{0}}\right)$, i.e., $V_{\mathrm{cmax}\mathrm{25}}\propto \mathrm{1}/{\mathrm{SLA}}_{\mathrm{0}}$. Using the default NETT PFT values for SLA${}_{\mathrm{0}}$ and $m$ led to the development of large and thin leaves with reduced ${\mathrm{N}}_{\mathrm{a}}$ and $V_{\mathrm{cmax}\mathrm{25}}$, resulting in excessive LAI and significant down-regulation of GPP. Smaller SLA${}_{\mathrm{0}}$ values were attempted (manual trial and error), with $m$ values constrained by Eq. (A2), the SLA${}_{\mathrm{0}}$ value, and the site observations of LAI and C${}_{\mathrm{l}}$ mentioned above, aiming to minimize model errors in the simulation of GPP and LAI. $\mathrm{SL}{A}_{\mathrm{0}}=\mathrm{0.006}$ m${}^{\mathrm{2}}$ leaf g C${}^{-\mathrm{1}}$ and $m=\mathrm{0.000985}$ m${}^{\mathrm{2}}$ ground g C${}^{-\mathrm{1}}$ were found to significantly improve the simulations and were adopted instead of the default values. Measurements reported by Woodruff et al. (2004) indicate that the ratio of leaf dry mass to leaf area reaches 263 g m${}^{-\mathrm{2}}$ leaf near the canopy top at Wind River (their Fig. 6). Assuming that the mass of carbon is 50 % of the dry mass, the observed value corresponds to 131.5 g C m${}^{-\mathrm{2}}$ leaf, i.e., an SLA${}_{\mathrm{0}}$ value of 0.0076 m${}^{\mathrm{2}}$ leaf g C${}^{-\mathrm{1}}$, indicating that the optimized SLA${}_{\mathrm{0}}$ value moved in the right direction from the default NETT PFT value (0.0100 down to 0.0060 m${}^{\mathrm{2}}$ leaf g C${}^{-\mathrm{1}})$, but ended up slightly lower than the observed value.

Tree mortality

Results reported by van Mantgem et al. (2009) indicate an increasing trend in plant mortality rates ($M$, yr${}^{-\mathrm{1}})$ for Pacific Northwest forests, with $M$ growing from $\approx \mathrm{1}$ % yr${}^{-\mathrm{1}}$ in 2000 towards $\approx \mathrm{1.5}$ % yr${}^{-\mathrm{1}}$ in 2010. In CLM, a default rate of $M=\mathrm{2}$ % yr${}^{-\mathrm{1}}$ is used for all vegetation types, which was found to be excessive at Wind River, leading to a reduced modeled forest biomass. $M=\mathrm{1.5}$ % yr${}^{-\mathrm{1}}$ was found to yield results closer to site observations and was therefore adopted in this study.

Temperature sensitivity coefficient ($Q_{\mathrm{10}})$

The effect of temperature on maintenance respiration (component of autotrophic respiration) in CLM is calculated via a $Q_{\mathrm{10}}$ formulation, in which the base rate of maintenance respiration is multiplied by $Q_{\mathrm{10}}^{(T_{\mathrm{a}}-T_{\mathrm{ref}})/\mathrm{10}}$, where $Q_{\mathrm{10}}$ is a temperature sensitivity coefficient, $T_{\mathrm{a}}$ is air temperature, and $T_{\mathrm{ref}}$ is a reference temperature. For the maintenance respiration cost for live fine roots, soil temperature at the respective soil layer ($T_{\mathrm{s},i})$ is used instead of $T_{\mathrm{a}}$. Similarly, the effect of temperature on decomposition (and therefore on heterotrophic respiration) is also calculated via a $Q_{\mathrm{10}}$ formulation, in which the base rates of decomposition are multiplied by $Q_{\mathrm{10}}^{(T_{\mathrm{s},i}-T_{\mathrm{ref}})/\mathrm{10}}$. In CLM, a default $Q_{\mathrm{10}}$ of 1.5 is used for both maintenance respiration and decomposition. However, nighttime CO${}_{\mathrm{2}}$ flux measurements above the canopy at Wind River, which would include the sum of autotrophic and heterotrophic respiration, indicate a $Q_{\mathrm{10}}$ of 2.49 (Misson et al., 2007). By adjusting CLM's $Q_{\mathrm{10}}$ to 2.5 for both maintenance respiration and decomposition, the seasonal behavior of ecosystem respiration better corresponded with observed values. This was especially the case for heterotrophic respiration, reducing the model overestimation during winter and the model underestimation during summer.

Soil hydraulic properties

Initial runs indicated poor performance of CLM in the simulation of soil water content at the Wind River site (strong dry bias), which resulted in an unrealistic down-regulation of GPP due to soil moisture stress late in the dry summer season. When using the original soil hydraulic properties from CLM4.0 the results were greatly improved, with a wetter soil and a reduction of the unrealistic soil moisture stress. The observed improvement was likely related to a smaller subsurface runoff in CLM4.0 and consequently greater water retention in the soil. In CLM, subsurface runoff is proportional to a term representing the maximum drainage when the water table depth is at the surface ($q_{\mathrm{drai}}^{max})$. In CLM4.0, $q_{\mathrm{drai}}^{max}=\mathrm{0.0055}$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, while in CLM4.5 $q_{\mathrm{drai}}^{max}=\mathrm{10}\sin \mathit{\beta }$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, where $\mathit{\beta }$ is the mean grid cell topographic slope. Even for a small 1${}^{\circ }$ slope, $q_{\mathrm{drai}}^{max}$ is significantly larger than in CLM4.0 (0.1745 kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$. The soil hydraulic properties from CLM4.0 were therefore used in this study.

Root distribution

In CLM, root distribution over soil depth is calculated as in Eq. (4). Root fraction ($r_i)$ in combination with a plant wilting factor ($w_i$, Eq. 5) for each soil layer $i$ is used to calculate an integrated soil moisture stress parameter in CLM, ${\mathit{\beta }}_t$ (Eq. 3), which down-regulates stomatal conductance in the model (Eq. 2).

Shaw et al. (2004) provide a good description of rooting depth at Wind River: “Plant roots are concentrated above 50 cm in soil profiles; however, roots as deep as 2.05 m have been observed in younger forests growing on nearly identical soils (T. Hinckley, personal communication). Many coarse roots of Douglas fir extend to depths greater than 1.0 m. Tip-up mounds of windthrown western hemlock trees typically have a classic flat root plate indicative of shallow rooting” (Douglas fir and western hemlock are the dominant species at the site). With the default NETT PFT root distribution parameters in Eq. (4) ($r_{\mathrm{a}}=\mathrm{7}$ m${}^{-\mathrm{1}}$ and $r_{\mathrm{b}}=\mathrm{2}$ m${}^{-\mathrm{1}})$, the total root fraction in the top 46 and 130 cm of soil is 78 and 96 %, respectively (note the small fraction of roots at depths below 1.3 m, 4 %). The site description above (Shaw et al., 2004) suggests that the default parameters are inadequate at Wind River, resulting in a too-shallow rooting profile.

In this study the NETT PFT $r_{\mathrm{b}}$ parameter was changed to 1 m${}^{-\mathrm{1}}$ (default CLM value for broadleaf evergreen temperate tree PFT), shifting roots towards deeper soil layers, in order to make water stored at deeper soil layers available to the trees and, along with the changes in the soil hydraulic properties discussed in Appendix A7, reduce the excessive late-summer soil moisture stress and down-regulation of GPP in the model. With the adjusted $r_{\mathrm{b}}$ parameter, the total root fraction in the top 46 and 130 cm of soil is 67 and 86 %, respectively (14 % below 1.3 m), which seems more reasonable based on Shaw et al. (2004) and the fact that Douglas fir trees at the site are about 500 years old and 40–65 m tall. The adjustment of soil moisture stress in CLM via root distribution was therefore physically justified.

The plant wilting factor, $w_i$, offers an additional path for adjustment of the simulated soil moisture stress, but it was not investigated in this study.

Modeled latent heat flux (a) and ${\mathit{\delta }}^{\mathrm{13}}$C of gross primary production (b, lines) for 1998–2006. The curves presented correspond to Bézier-smoothed daily averages as in Figs. 3 and 5. Results from two model runs are presented: CLM (calibrated model, solid red lines) and CLM${}^{\ast }$ (calibrated model using the default stomatal conductance parameters ($m_{\mathrm{bb}}$ and $b_{\mathrm{bb}}$; see Table 1), dashed red lines). The blue line and circles correspond to site observations. The circles in panel (b) are the monthly averages of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ reported by Lai et al. (2005).

[Figure omitted. See PDF]

Stomatal conductance

In CLM, leaf stomatal conductance is calculated based on the Ball–Berry model as described by Collatz et al. (1991) and implemented by Sellers et al. (1996) in the SiB2 model (see Eq. 2). The default values set for the parameters $m_{\mathrm{bb}}$ and $b_{\mathrm{bb}}$ in CLM for C${}_{\mathrm{3}}$ plants (9 and 10 mmol m${}^{-\mathrm{2}}$ leaf s${}^{-\mathrm{1}}$, respectively) were found to be inadequate at Wind River, leading to a significant overestimation of latent heat fluxes due to excessive plant transpiration (after the adjustments discussed in the aforementioned sections, which resulted in higher forest productivity). These default parameter values were established based on the values used in the SiB2 model (Sellers et al., 1996). In SiB2, however, a distinction was made for coniferous forests ($m_{\mathrm{bb}}=\mathrm{6})$ but was not carried over to CLM. Observations reported in the literature support this lower $m_{\mathrm{bb}}$ value for conifers (see for example the survey by Williams et al., 2004, Table 6.3, and Miner et al., 2017, Fig. 1). Conversely, $b_{\mathrm{bb}}$ values reported in the literature are highly variable (1–400 mmol m${}^{-\mathrm{2}}$ leaf s${}^{-\mathrm{1}}$ in the survey by Barnard and Bauerle, 2013, for a broad range of plant species). In CLM4.0, the default $b_{\mathrm{bb}}$ for C${}_{\mathrm{3}}$ plants is significantly smaller than in CLM4.5 (2 vs. 10 mmol m${}^{-\mathrm{2}}$ leaf s${}^{-\mathrm{1}})$ (Oleson et al., 2010). Values of $m_{\mathrm{bb}}=\mathrm{6}$ and $b_{\mathrm{bb}}=\mathrm{5}$ mmol m${}^{-\mathrm{2}}$ leaf s${}^{-\mathrm{1}}$ were found to greatly improve the modeled latent heat fluxes at the Wind River site and were therefore adopted in this study. The updated values also resulted in a great improvement in the simulation of ${\mathit{\delta }}^{\mathrm{13}}$C of leaves, SOM, and GPP. Figures A1 and A2 illustrate the impact of the stomatal conductance parameters on model performance, particularly in regards to latent heat fluxes and photosynthetic ${}^{\mathrm{13}}$C discrimination.

Modeled fraction of potential GPP (GPP/GPP${}_{\mathrm{pot}})$. Data points correspond to daily means averaged over 1850–2006 (calibrated CLM simulation). The fitted curve is $f_n(x$ $=$ day of year) $=$ $-\mathrm{1.39697}\times {\mathrm{10}}^{-\mathrm{14}}{x}^{\mathrm{6}}+\mathrm{1.71948}\times {\mathrm{10}}^{-\mathrm{11}}{x}^{\mathrm{5}}-\mathrm{8.26883}\times {\mathrm{10}}^{-\mathrm{9}}{x}^{\mathrm{4}}+\mathrm{1.90682}\times {\mathrm{10}}^{-\mathrm{6}}{x}^{\mathrm{3}}-\mathrm{1.97639}\times {\mathrm{10}}^{-\mathrm{4}}{x}^{\mathrm{2}}+\mathrm{0.0055728}x+\mathrm{0.966272}$ for 31 $\le x\le$ 332 and $f_n\left(x\right)=\mathrm{1}$ elsewhere. Note that GPP/GPP${}_{\mathrm{pot}}=\mathrm{1}-d$, where $d$ is the nitrogen down-regulation factor as defined in Eq. (9) within text.

[Figure omitted. See PDF]

Modeled gross primary production (a), latent heat flux (b), and ${\mathit{\delta }}^{\mathrm{13}}$C of gross primary production (c) for 1998–2006. The curves presented correspond to Bézier-smoothed daily averages as in Figs. 3 and 5. Results from two model runs are presented: CLM (calibrated model, red lines) and CLM${}^{\mathit{\#}}$ (calibrated model using $m_{\mathrm{bb}}=\mathrm{9}$ and the alternative nitrogen limitation scheme discussed in Appendix A9, black lines). In CLM${}^{\mathit{\#}}$, $V_{\mathrm{cmax}\mathrm{25}}$ was multiplied by a seasonally varying nitrogen down-regulation factor, calculated based on the mean (1850–2006) seasonal cycle of $\mathrm{GPP}/{\mathrm{GPP}}_{\mathrm{pot}}=\mathrm{1}-d$ in the CLM run ($f_n\left(x\right)$ in Fig. A3) subtracted by 0.35 (manual adjustment applied to avoid excessive productivity during the transient simulation). Blue lines and circles correspond to site observations. The circles in panel (c) are the monthly averages of ${\mathit{\delta }}^{\mathrm{13}}$C${}_{\mathrm{ER}}$ reported by Lai et al. (2005).

[Figure omitted. See PDF]

Comparison of CLM performance at Wind River when using default, “out-of-the-box” parameters (black lines) and calibrated parameters (red lines). Modeled values correspond to annual averages. Observations (average $\pm$ SD, blue points and error bars) are from the AmeriFlux database (based on Thomas and Winner, 2000, and Harmon et al., 2004).

[Figure omitted. See PDF]

Comparison of CLM performance at Wind River when using default, “out-of-the-box” parameters (black lines) and calibrated parameters (red lines). Observations (blue lines) are from the AmeriFlux database. For a clearer visualization, the data presented correspond to Bézier-smoothed daily averages as in Fig. 3.

[Figure omitted. See PDF]

It is important to highlight that the default nitrogen limitation scheme was used in the simulations. As discussed in Sect. 2.1, this scheme makes CLM a partially coupled model in respect to net leaf photosynthesis and stomatal conductance: while the potential GPP is down-regulated in response to nitrogen availability, stomatal conductance remains consistent with potential net leaf photosynthesis ($A_{\mathrm{n}})$. With this structure, CLM is expected to overestimate plant transpiration and photosynthetic ${}^{\mathrm{13}}$C discrimination. The calibration of the Ball–Berry stomatal conductance parameters discussed above, especially the significant reduction of $m_{\mathrm{bb}}$ from 9 to 6, must also have compensated for this structural issue within the model. Note that nitrogen down-regulation is significant at Wind River, peaking at $\sim$ 0.25 (GPP/GPP${}_{\mathrm{pot}}$ $=$ 0.75) in May (Fig. A3).

When using the default nitrogen limitation scheme in CLM, the modeled ${}^{\mathrm{13}}$C discrimination values reported by Raczka et al. (2016) for the Niwot Ridge AmeriFlux site (also a coniferous forest site) were significantly overestimated, i.e., ${\mathit{\delta }}^{\mathrm{13}}$C values of GPP and biomass were significantly smaller than observations. To improve the simulation, Raczka et al. (2016) removed the post-photosynthetic nitrogen down-regulation of $A_{\mathrm{n}}$ and ${\mathrm{GPP}}_{\mathrm{pot}}$ ($d=\mathrm{0}$; see Eq. 9) and included a foliar nitrogen-limiting factor in the calculation of $V_{\mathrm{cmax}\mathrm{25}}$, making the model fully coupled in respect to net leaf photosynthesis and stomatal conductance. With this configuration, their simulation of ${}^{\mathrm{13}}$C discrimination improved significantly, but the values still presented a small overestimation in respect to the site observations. According to Raczka et al. (2016), overestimation of $g_{\mathrm{s}}$ due to an inadequate $m_{\mathrm{bb}}$ value (too high) could be a reason for the mismatch (they used the default value of 9 in their simulation).

The alternative nitrogen limitation scheme (via $V_{\mathrm{cmax}\mathrm{25}}$ down-regulation, as in Raczka et al., 2016) was also investigated here. The simulation of LE, GPP, and ${}^{\mathrm{13}}$C discrimination when using this configuration and the default $m_{\mathrm{bb}}$ value of 9 was found to be similar to the results when using the default nitrogen limitation scheme and $m_{\mathrm{bb}}=\mathrm{6}$ (Fig. A4). The results in Fig. A4 indicate that the calibration of $m_{\mathrm{bb}}$ from 9 to 6 represents a tradeoff with the approach to nutrient limitation, compensating for elevated, nitrogen-unlimited (potential) net leaf photosynthesis used in the calculation of $g_{\mathrm{s}}$.

CLM performance: default vs. calibrated parameters

In order to illustrate the effect of altering the model parameters discussed in this Appendix (see summary of changes in Table 1), Figs. A5 and A6 compare the performance of CLM for key model outputs when using “out-of-the-box” parameters and calibrated parameters. Note the significant improvement in the simulation of LAI, biomass, and CO${}_{\mathrm{2}}$ $/$ H${}_{\mathrm{2}}$O fluxes.

The authors declare that they have no conflict of interest.

Acknowledgements

This research was supported by the US Department of Energy's Office of Science, Terrestrial Ecosystem Science Program, under award number DE-SC0010624. BMR and DRB were supported by the US Department of Energy's Office of Science, Terrestrial Ecosystem Science Program, under award number DE-SC0010625. PET was supported by the US Department of Energy's Office of Science, Biological and Environmental Research, Accelerated Climate Modeling for Energy project. We acknowledge the Wind River Field Station AmeriFlux site (US-Wrc, PIs: Kenneth Bible, Sonia Wharton) for its data records. Funding for AmeriFlux data resources was provided by the US Department of Energy's Office of Science. Data and logistical support were also provided by the US Forest Service Pacific Northwest Research Station. We gratefully acknowledge the comments and suggestions of the editor and three anonymous reviewers, which have greatly improved our paper. Edited by: Andreas Ibrom Reviewed by: three anonymous referees