Linking tree physiological constraints with predictions of carbon and water ﬂ uxes at an old-growth coniferous forest

. Old-growth coniferous forests of the Paci ﬁ c Northwest are among the most productive temperate ecosystems and have the capacity to store large amounts of carbon for multiple centuries. To date, there are considerable gaps in modeling ecosystem ﬂ uxes and their responses to physiological constraints in these old-growth forests. These model shortcomings limit our ability to understand and project how the old-growth forests of the Paci ﬁ c Northwest will respond to global climate change. This study applies the cohort-based Ecosystem Demography Model 2 (ED2) to the Wind River Experimental Forest (Washington, USA), a well-studied old-growth Douglas- ﬁ r – western hemlock ecosystem. ED2 is calibrated and validated using an extensive suite of forest inventory, eddy covariance, and biophysical observations. ED2 is able to reproduce observed forest composition and canopy structure, and carbon, water, and energy ﬂ uxes at the site. In the simulations, the effect of limited water supply on ecosystem carbon ﬂ uxes is mediated primarily by the forest ’ s gross primary productivity (GPP) response, rather than its heterotrophic respiration response. The simulation indicates that stomatal conductance is mainly determined by soil moisture during periods of low vapor pressure de ﬁ cit (VPD). However, when VPD is high, stomatal conductance is greatly reduced regardless of soil moisture status. During summer droughts, reduced soil moisture and increased VPD result in considerable stomatal closure and GPP reduction, which in turn decreases net carbon uptake. Cohort-based scheme integrates all canopy layers (species) that have distinct sensitivity to microclimate and respond distinctly to drought. This study is an initial ﬁ rst step to explore the potential importance of cohort-based model in simulating forest with complex structure, and to lay the foundation for applying cohort-based model at regional scales across the Paci ﬁ c Northwest.


INTRODUCTION
Old-growth coniferous forests, formerly widespread in the U.S. Pacific Northwest, have unique leaf physiological and conifer forest structural characteristics (White et al. 2000, Hessl et al. 2004. They are typically dominated by tall evergreen conifer species, such as Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco), which can live for 400-1000 yr (Franklin andDeBell 1988, Falk et al. 2005), and western hemlock (Tsuga heterophylla (Raf.) Sarg.). These forests potentially store more aboveground carbon than any other forest type in the world (Franklin and Waring 1980, Schimel et al. 2000, Smithwick et al. 2002, Suchanek et al. 2004, Hudiburg et al. 2009) and are extremely heterogeneous in canopy structure, tree size, and understory patchiness (Franklin et al. 1981). As one of the most productive temperate forest ecosystems (Franklin andWaring 1980, Field andKaduk 2004), old-growth coniferous forests have the capacity to both assimilate large amounts of carbon (Paw U et al. 2004) and store carbon for multiple centuries in live, woody debris, and soil carbon pools ). Estimates of late-successional and old-growth forests in the Northwest Forest Plan area (western Oregon and Washington, and northern California) put the estimate at 3.0 Mha, with less than a 3% decrease on federal lands by 2012 (Davis et al. 2015). Under a projected warming climate, more frequent drought-and heat-related physiological stresses can substantially affect carbon and water cycling of the old-growth forests and the feedback to atmosphere (Law and Waring 2015). Given the importance of these forests to regional carbon and water budgets and the provision of ecosystem services, along with their societal value, it is imperative to understand the mechanisms that control carbon and water dynamics and responses to environmental changes.
Given the structural complexities of the oldgrowth forest , species composition and size distribution may be important factors in shaping carbon and water dynamics between the ecosystem and atmosphere. The highly decoupled microclimate across different canopies drives distinct responses to physiological stresses among species/sizes. As climate warms in the region, potentially stronger water stress and higher mortality rates may alter forest composition (Littell et al. 2010, Spies et al. 2010, drive long-term change in forest productivity (Latta et al. 2010), and raise new challenges for future conservation planning (Carroll et al. 2010, Kim et al. 2017. Moreover, potential droughtrelated increase in future wildfire activities may dramatically change the forest structure in the region and release large amounts of carbon to the atmosphere Williams 2016, Riley andLoehman 2016). The ongoing transition to more shade-tolerant species-the replacement of Douglas-fir with western hemlock (Bible 2001) -may change the ecosystem-level water-use efficiency and carbon assimilation rate in the region. Therefore, models that can reproduce species composition, canopy structure, and competitions among species/sizes are best equipped to evaluate the forest responses to climate change.
To date, several terrestrial biome models have been applied to the evergreen coniferous forests in the Pacific Northwest (Waring et al. 2011, Shafer et al. 2015, Turner et al. 2015. Although these models have characterized important aspects of forest ecology and atmosphere-biosphere dynamics, their simplified representation of species composition and canopy structure may affect their simulation of natural forest stands of the Pacific Northwest in important ways. For example, in Biome-BGC, only one plant functional type (PFT) occupies each site, and no distinct age or size classes are represented (Turner et al. 2015(Turner et al. , 2016. 3-PG is a physiologically based stand growth model developed principally for plantation forestry applications but applied widely across the Pacific Northwest (Landsberg andWaring 1997, Waring et al. 2011). 3-PG simulates stem growth and leaf area dynamics of individual species, but it greatly simplifies or does not represent a range of forest processes and characteristics, including respiration, vertical canopy structure, and sub-monthly growth (Waring and McDowell 2002). The Lund-Potsdam-Jena dynamic global vegetation model (LPJ-DGVM;Shafer et al. 2015) captures a suite of biophysical processes but simplifies canopy structure: While multiple PFTs may occupy a site, each PFT is prescribed as only one single average big-leaf. Duarte et al. (2017) evaluated the performance of the Community Land Model version 4.5 (CLM4.5) at Wind River Experimental Forest, but the simulation represents the forest stand as a single needleleaf PFT. Pitfalls of oversimplification of vegetation composition and parameterization in CLM have been recognized by recent studies (Law and Waring 2015), and multiple PFTs with unique physiological parameters are required to represent forest types in the western United States (Buotte et al. 2019). At Wind River site, the old-growth forest has remarkably complex structure, with co-dominant but ecologically distinct conifer species. The potential importance of competition and succession dynamics among the cohorts occupying the multi-layered canopy at Wind River cannot be explored with the above models.
A second-generation dynamic global vegetation model (DGVM; Fisher et al. 2010Fisher et al. , 2015Fisher et al. , 2018 has the ability to simulate both sub-daily biophysical processes and longer-term demographic processes, thereby enabling exploration of the interaction of the shorter-and longer-term processes in response to climate drivers. Cohortbased models such as GAPPARD (Scherstjanoi et al. 2014) and TREEMIG  can better represent the heterogeneity of plants than big-leaf models and have lower computational cost than gap and individual-based models (Dietze and Latimer 2011, Christoffersen et al. 2016, Fischer et al. 2016. In this study, we applied a cohort-based model, the Ecosystem Demography Model 2 (ED2; Medvigy et al. 2009), at the Wind River Experimental Forest, an old-growth forest in Washington with a wealth of observation data and studies ). ED2 groups plants not only into multiple PFTs, but also into numerous size classes, which enables vertical interactions and resource competitions among different sizes of trees. We assumed that photosynthesis at this mature site was primarily limited by light and water availability (Brooks and Coulombe 2009, Tarvainen et al. 2016, Cornejo-Oviedo et al. 2017 and that variability in observed water and carbon fluxes during the 18yr study period (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) was primarily driven by climate. Therefore, we did not simulate nitrogen limitation during this exercise.
Our overall goal was to assess the ability of ED2 to simulate the sub-daily biophysical processes and the stand structure at Wind River, and to conduct an initial exploration of how stand structure interacts with the biophysical processes. In particular, we addressed the following questions: 1. How well can ED2 represent vegetation composition, stand structure, and ecosystem fluxes at Wind River? 2. How do carbon and water fluxes vary among the cohorts that comprise the simulated stand? 3. Are there significant differences in the way the cohorts respond to changes in the observed climate?
By pursuing these questions, we aim to establish the potential utility of using ED2 to explore long-term vegetation response under climate change scenarios (e.g., representative concentration pathways; Van Vuuren et al. 2011), and to lay the foundation for applying a secondgeneration DGVM across the Pacific Northwest.

Study site
Wind River Experimental Forest (45.82°N, 121.95°W) is located in the T.T. Munger Research Natural Area, in the southern Cascade Range of Washington State. The forest is situated at 371 m elevation within the relatively flat Wind River Valley on a shallow slope (~3.5%) between two hills, which protect the study area from severe windstorms ). The soils are generally coarse in texture, well drained, and nutrient poor. Winters are wet and mild with occasional snow, and summers are warm and dry (Falk et al. 2005(Falk et al. , 2008. During the period of 1998-2015, the total water-year (October-September) precipitation varied from~1269 mm to~2835 mm, with only~7% falling from June through September. Over the same period, the average annual air temperature varied from 8.1°C to 10.9°C. The winter monthly mean temperature in this period ranged from 0.8°C to 3.6°C, and the summer monthly mean temperature varied between 16.2°C and 20.1°C.
The forest at Wind River has a species composition characteristic of old-growth forests in the Pacific Northwest. It represents the upper boundary of several forest ecosystem gradients in age, biomass, stature, tree density, and structural complexity , Suchanek et al. 2004, Falk et al. 2005. The forest is composed of Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco), western hemlock (Tsuga heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata Donn ex D. Don), Pacific yew (Taxus brevifolia Nutt.), Pacific silver fir (Abies amabilis (Douglas ex Loudon) Douglas ex Forbes), noble fir (Abies procera Rehder), grand fir (Abies grandis (Douglas ex D. Don) Lindl.), and western white pine (Pinus monticola Douglas ex D. Don). The forest has a stand density of 427 trees/ha and basal area of 82.9 m 2 /ha. Douglas-fir dominates the forest in basal area (~43%) and wood volume (~50%). However, western hemlock, a late-successional species, comprises the majority of tree stems (~52%) while accounting for 32.4% of the total basal area . Deciduous trees comprise a small component of stand biomass (Thomas and Winner 2000).
The current stand is thought to have originated from a stand-replacing fire approximately 500 yr ago and has experienced minimal anthropogenic impact since ). Douglas-fir trees dominate the upper canopy (45-60 m), while western hemlocks dominate the middle canopy (15-45 m), along with occasional western redcedars. Western hemlocks have the greatest foliage area (~53%; Parker et al. 2004). Estimates of one-sided leaf area index (LAI) at the stand range from 6.3 to 9.2, with little seasonality (Parker et al. 2002, Roberts et al. 2004. The vertical distribution of LAI has a bottom-heavy pattern, with the maximum LAI occurring between 15 and 30 m aboveground Winner 2000, Parker et al. 2002). Total carbon storage at the site is estimated at 61.9 kg C/m 2 , with 39.8 kg C/m 2 stored in live biomass, 9.6 kg C/m 2 in woody detritus, 9.3 kg C/m 2 in mineral soil, and~3.2 kg C/m 2 in fine non-woody litter ).

ED2 model simulation
The ED2 model is a cohort-based terrestrial biosphere model that simulates a comprehensive suite of forest ecosystem processes, including fastresponse (half-hourly) biophysical and physiological processes, as well as long-term processes related to changes in forest structure, demography, and biogeochemistry (Medvigy et al. 2009). ED2 represents the landscape as a grid, and within each grid cell, a set of tree cohorts are simulated using partial differential equations to approximate average tree size and density. Species are grouped into plant functional types (PFTs) which are associated with a set of plant physiological traits that govern how a tree interacts with its environment. The default parameterization for the two conifer PFTs in ED2 does not adequately represent the species that currently occur in the Pacific Northwest. Therefore, we used field data to adjust six key PFT-specific parameters and to modify the dbh-toheight allometric functions for the northern pine and the late-successional conifer PFTs, to represent the two dominant tree species, Douglas-fir and western hemlock (Table 1), respectively. A complete list of tree species found at the study site and their associations with the two PFTs are given in Appendix S1: Table S1. The source code for the modified version of ED2 is provided in Data S1.
We ran ED2 for a single grid cell, representing the approximate measurement footprint of the flux tower at the Wind River Canopy Crane Research Facility. As with other ED2 studies Moorcroft 2012, Trugman et al. 2016), we initialized the simulation using forest inventory data . We initialized the simulation with three PFTs: a northern pine PFT representing Douglas-fir, with an initial average diameter at breast height (dbh) of 109.3 cm; a late-successional conifer PFT representing western hemlock, with an initial average dbh of 27.5 cm; and a latesuccessional broadleaf deciduous tree PFT representing the deciduous species at the site (e.g., vine maple; Appendix S1: Table S1), with an average dbh of 4.8 cm. The initial total soil carbon stock was set to 22.1 kg C/m 2 .
Simulations were driven with 1998-2015 halfhourly meteorological and radiation data from Wind River Crane Site (Wharton and Falk 2016), obtained from the AmeriFlux data portal (http:// ameriflux.lbl.gov/sites/siteinfo/US-Wrc, Wharton, S. 1998Wharton, S. -2016. We spun up the model to stabilize soil moisture values and density-dependent mortality. For the spin-up phase, we used a detrended version of the 1998-2015 meteorological and radiation data, repeating it three times, for a total of 54 simulated years. We assumed that the 18-yr time series was representative of the annual variability of the entire spin-up period.

Model evaluation and analysis
To evaluate the ED2 model skill for predicting stand structure at Wind River, we compared simulated plant densities and tree size distributions of two PFTs with the inventory data. To evaluate how well ED2 can simulate ecosystem dynamics, we compared model outputs against an extensive 18-yr dataset of concurrent observations at Wind River including (1) leaf area index (LAI) stratification within the canopy; (2) sensible and latent heat fluxes; (3) leaf temperature at different canopy heights; (4) soil temperature at 0 cm and 15 cm, and soil moisture at 30 cm and 50 cm; and (5) daily gross primary production (GPP), ecosystem respiration (R eco ), net ecosystem production (NEP), and averaged daily heterotrophic respiration (R H ) (Appendix S1: Table S2).
To quantify the variability of carbon and water fluxes among cohorts, we compared GPP among all the cohorts. For comparing stomatal conductance, leaf temperature, and leaf-level VPD, we grouped the cohorts into three thickness-even strata: lower (<20 m), middle (20-40 m), and upper (>40 m) canopies. To examine the ❖ www.esajournals.org variability of responses to drought, we compared stomatal conductance, leaf temperature, and VPD of across different heights of canopies between a typical summer (July 2014) and an exceptionally hot and dry summer (July 2015).
To characterize the response of the whole stand to drought, we selected four years with contrasting climates for further examination: a relatively wet and cool year (2012), an exceptionally dry year (2013), a hot year (2014), and an exceptionally hot and dry year (2015). We compared water and carbon fluxes of the four selected years to values from the baseline period of 2002-2011. We examined GPP, NEP, stomatal conductance, evapotranspiration (ET), and soil water content (SWC).

Model-observation comparison
ED2 was able to capture many aspects of PFT composition and canopy structure (Fig. 1). The modeled tree densities for both northern pine (NP) and late-successional conifer (LSC) PFTs agreed well with inventory data for Douglas-fir and western hemlock, respectively ( Fig. 1a, b).
The modeled slightly underestimated the density of large trees and failed to represent the mid-size ones. The small NP trees are mostly saplings, which have negligible carbon flux and water flux. In terms of tree height, ED2 simulation showed good agreement with measured tree heights of Douglas-fir at the site (Fig. 1c, d). For LSC trees, ED2 slightly underestimated the heights of western hemlock across all size classes. Total LAI (~6.3) simulated by ED2 was close to the value estimated by Parker et al. (2002) using the Beer-Lambert inversion method. However, the simulated vertical distribution of LAI is unlike the observed pattern from Parker et al. (2004), where most leaf biomass is concentrated at heights below 20 m (Fig. 2). This mismatch is primarily due to the fact that ED2 places all the leaves of a given cohort at the height of the cohort. For example, ED2 places all the leaves of a 20 m tall cohort at 20 m height.
ED2 reproduced the observed magnitude, variability, and seasonality of both sensible (H) and latent heat (LE) fluxes measured by an eddy covariance system above the canopy (Fig. 3). Linear regression analysis shows that the modeled daily H and LE have good (R 2 = 0.73) and ED2 captured both the magnitude and variability of observed GPP, R eco , and NEP on a seasonal timescale (Fig. 4). On a daily timescale, however, ED2 often underestimated the variability of GPP, R eco , and NEP, compared to the eddy covariance (EC) data. The discrepancies between measured and modeled GPP and R eco did not show a clear seasonal pattern. ED2 poorly models fluxes responding to temperature spikes in the summers of , 2007, and 2013). Both modeled and observed daily NEP values ranged between À5 and 5 g CÁm À2 Ád À1 . However, ED2 generally fails to capture peak NEP values in the summer. The slope of regression lines with the intercept forced through 0 was close to but slightly less than unity (GPP: 0.92, R eco : 0.90, NEP: 0.95), with relatively high R 2 (GPP: 0.64, R eco : 0.71, NEP: 050) and small RMSE (GPP: 1.31, R eco : 1.05, NEP: 1.21 g CÁm À2 Ád À1 ; Fig. 4). These lower-than-unity slopes represent a tendency by ED2 to overestimate high flux values (e.g., high NEP in mid-April) and underestimate low flux values (e.g., low GPP in winter).
The 18-yr averaged carbon fluxes from ED2 reproduced the seasonality of the observations, except in the case of simulated R eco , the timing of which appears to precede observations by 1-2 weeks in the spring (Fig. 5a). The averaged   . Positive values denote upward fluxes (i.e., ground to atmosphere). The intercepts of linear fits between model and eddy covariance data are forced to zero. In the right-hand panels, the thick black line represents the regression line, with the 95% confidence interval (thin black lines). In the linear fit for LE, points during 2009-2013 are ignored because of potential sensor calibration error. ❖ www.esajournals.org 7 April 2019 ❖ Volume 10(4) ❖ Article e02692 peaks in modeled and observed GPP were both about 7 g CÁm À2 Ád À1 in mid-June, with slightly lower peak R eco in mid-July. This asymmetric pattern of GPP and R eco resulted in net carbon uptake (i.e., positive NEP) starting from mid-February and ending around late June. The crossover of GPP and R eco occurred around late June; thereafter, the forest was a carbon source until early October when the rainy season started again (Fig. 5b). Therefore, the asymmetry of GPP and R eco determined the net C change throughout a year. This asymmetry is consistent with the eddy covariance-based analysis by Falk et al. (2008). In the spring, ED2 overestimates R eco , resulting in an underestimation of NEP in the spring. Compared with the measured NEP, ED2 slightly underestimated the daily net carbon sink in April by~0.8 g CÁm À2 Ád À1 and overestimated the net carbon source in late July through August by~0.4 g CÁm À2 Ád À1 . Annually, the modeled annual mean GPP (1417.6 AE 51.6 g CÁm À2 Áyr À1 ) and R eco (1382.3 AE 48.9 g CÁm À2 Áyr À1 ) summed to 35.3 AE 66.6 g C/m 2 NEP. These estimates are very close to the measured GPP, R eco , and NEP estimates of 1289.2 AE 215.4, 1255.6 AE 216.4, and 26.2 AE 88.1 g CÁm À2 Áyr À1 , respectively (Wharton and Falk 2016). In other words, both ED2 and EC data show that the Wind River old-growth forest functioned as a small annual carbon sink during the 18-yr measurement period. Additionally, ED2 captured the difference between the measured air temperature and top canopy temperature, yet poorly captured the temperature difference between top and bottom canopies (Appendix S1: Fig. S1). The model well reproduced soil temperature and soil moisture Fig. 4. Comparison of ED2-modeled and data-derived daily GPP, R eco , and NEP, from 1998 to 2015. Here, NEP = ÀNEE; therefore, positive NEP represents net ecosystem C uptake. The right panels represent linear fittings between modeled and measured fluxes, with intercept forced to zero. In the fitting panels, each season is represented in a unique color: winter (blue), spring (green), summer (orange), and autumn (brown). across soil depths (Appendix S1: Figs. S2 and S3). The modeled soil respiration also agreed well with the chamber measurements (Appendix S1: Fig. S4).

Heterogeneity in cohorts
Cohorts simulated by ED2 had distinctly different contributions to carbon fluxes (e.g., GPP; Fig. 6). On an annual basis, tall cohorts (e.g., dbh > 100 cm) accounted for 26% of the total stand GPP. Although the smaller cohorts had much lower GPP, their higher densities still led to the greatest proportion of total GPP. For example, cohorts with dbh < 50 cm contributed 56% of the total stand GPP, and cohorts with dbh between 50 and 100 cm contributed 18%.
The variations in stomatal conductance in responses to drought condition drove distinct changes in GPP across multiple cohorts (Fig. 8). For example, the extreme drought in July 2015 exerted a much stronger effect on the western hemlock occupying the middle-lower canopies (e.g., 10-40 m) than on the Douglas-fir at the upper canopies, reducing the entire stand GPP by about 36%. Although the extreme drought also reduced GPP of Douglas-fir in the upper canopies, the decrease was small relative to the base value.

Ecosystem-level responses to various hydrologic conditions
We compared results from 2012 to 2015, which were droughts years, to results from the baseline period of 2002-2011, to explore ecosystem-level response to drought. Integrating across all cohorts, the ecosystem-level photosynthetic capacity was increasingly sensitive to available water from 2012 to 2015 (Fig. 9). For instance, spring precipitation (March through May) in 2012 was 807 mm higher than the baseline 2002-2011 period, contributing to a 24 g C/m 2 per growing season increase in GPP, despite the average air temperature being~0.6°C lower than the baseline. In 2014, the warmer (+1.2°C) and wetter (+108 mm) spring also stimulated increased photosynthesis (+13 g CÁm À2 Áseason À1 ), thereby slightly compensating for the reduction in spring net carbon uptake (À15 g CÁm À2 Áseason À1 ) due to enhanced plant and soil respiration. In contrast, 2013 experienced an exceptionally dry (À144 mm) and slightly warmer (+0.7°C) growing season compared to the baseline period, and this led to continuous reduction in GPP from April through August (À24 g C/m 2 per growing season). The growing season of 2015 was even drier and warmer than that of 2013: Spring and summer temperatures were +2.2°C and +2.5°C warmer, respectively, than the baseline period. However, spring and summer precipitation amounts were below baseline averages, which resulted in a substantial reduction in growing season GPP (À133 g CÁm À2 Ágrowing season À1 ). The reduction in GPP in 2015, in tandem with a slight increase in ecosystem respiration (34 g CÁm À2 Ágrowing season À1 ), led to a net decrease in NEP (À167 g CÁm À2 Ágrowing season À1 ) for that year. Because ecosystem respiration was mostly sensitive to air temperature, the warmer growing season exacerbated the drought-induced reductions in NEP in 2013 and 2015 (Fig. 9).
The reduction of carbon uptake in dryer and warmer years was associated with anomalous water fluxes (Fig. 10), which in turn were driven by anomalous VPD and SWC. In 2013 and 2015, elevated VPD from April through August (+0.1 kPa and +0.4 kPa, respectively, relative to the baseline) reduced canopy stomatal conductance by 0.1 and 0.2 gÁm À2 Ás À1 , respectively. This reduction was partially responsible for the contemporaneous decrease in GPP (Fig. 9), and contributed to lower evapotranspiration (ET), as seen in year 2015 (À29 kgÁm À2 Ágrowing season À1 ). In 2013, reduced stomatal conductance did not lead to a net decrease in ET during the growing season (+7.9 kgÁm À2 Ágrowing season À1 ), likely due to higher interception loss of dew from foliar and trunk surfaces, given the modestly higher precipitation in April and May of that year. A similar phenomenon was simulated in 2014: Lower VPD reduced stomatal conductance in April and May, but it did not reduce modeled ET (+4.6 kgÁm À2 Ámonth À1 ) because of higher precipitation during those two months (+86 mm). In 2012, VPD was substantially lower than the baseline in June (À0.27 kPa) and July (À0.40 kPa), which drove higher stomatal conductance during those two months (+0.08 and +0.06 gÁm À2 Ás À1 , Fig. 6. The ED2-modeled GPP across different size of cohorts. The size of cohorts here is defined by dbh. ❖ www.esajournals.org respectively). Among the four years, only 2012 had greater SWC than the baseline for most of the growing season, despite stronger photosynthesis in the summer, because of soil water recharge from an anomalously wet spring. In the other three years, enhanced GPP and ET in early spring did not necessarily reduce SWC in late spring and summer because stomatal conductance was also lower than the baseline. SWC anomalies appear more sensitive to the timing of Fig. 7. Modeled stomatal conductance (gs), leaf temperature, and leaf vapor pressure deficit (VPD) across a series of tree height intervals in July 2014 (a-c) and July 2015 (d-f). A specific height does not represent the canopy corresponding to that height, but instead represents tree height. precipitation and the variability of air temperature than to ET withdrawals.
During the growing season, stomatal conductance was regulated by both SWC and VPD (Fig. 11a). For this analysis, we calculated the modeled leaf-specific VPD (VPD leaf ), since leaflevel measurements emphasize the strong dependence of stomatal conductance on leaf-to-air VPD (Schulze andHall 1982, Beer et al. 2009). Under low-VPD leaf conditions, stomatal conductance was regulated primarily by SWC, but as VPD leaf increased, the role of SWC diminished, with extreme VPD leaf levels largely shutting down stomatal conductance. For example, on days with midday VPD leaf below 2 kPa, SWC ranged from 0.13 to 0.34 m 3 /m 3 . In contrast, on days with midday VPD leaf higher than 8 kPa, SWC was limited to the range of 0.15-0.22 m 3 /m 3 . We illustrate the exceptional levels of drought stress in 2015 by comparing the distribution of midday VPD values for the growing season of 2015 with those taken from the entire 18-yr period (Fig. 11b). The growing season in 2015 was marked by notably more frequent extreme VPD leaf days (>6 kPa) than the entire 18-yr period, which reduced stomatal conductance, ET, and GPP.

DISCUSSION
The unique physiology and longevity of coniferous trees lead to their dominance in the Pacific Northwest (Bible 2001). Broadly generalized parameters for conifers may not work satisfactorily when modeling species in the northwestern United States (White et al. 2000, Hessl et al. 2004). Our use of site-specific parameterization for two coniferous PFTs in ED2 further supports this point. For example, considerably higher SLAs (specific leaf areas) for the two simulated conifers than the default values allow for more biomass allocation to metabolic rather than structural components of leaves. Increasing photosynthetic area per foliage mass of the lower canopy or understory conifer trees effectively increases light-use efficiency (Albani et al. 2006). Consequently, the density-dependent mortality associated with carbon starvation is significantly reduced in the model. Moreover, the other adjusted physiological parameters (Table 1)  together reduce the carbon cost and increase the productivity of coniferous PFTs, thereby increasing their competitive advantage against temperate deciduous PFTs.
With the improved values for key physiological and allometric parameters, ED2 was able to capture the general patterns of carbon, water, and energy fluxes at the site. The increase in summer stomatal conductance with tree height (Fig. 7) simulated by ED2 is at first puzzling given the opposite phenomenon reported by many (Bond et al. 2007, Woodruff et al. 2009). In the ED2 simulation, the upper-canopy foliage experiences more optimal leaf temperature and received more photosynthetically active radiation (PAR) in the summer than the lower canopy, resulting in higher photosynthetic rates and stomatal conductance. While ED2 may be overestimating the contrast between the upper and lower canopies, it remains possible that in highly shaded stands with sufficient soil moisture, the upper canopy maintains higher levels of productivity than the lower canopy during the summer. The mismatch in leaf skin temperature is partly due to ED2 locating all of a cohort's leaf biomass at a single height. The overestimation of radiative heating within canopy leads to overestimation of autotrophic respiration, which may indicate the model's over-sensitivity to temperature. As a consequence, the simulated R eco increases earlier in the beginning of growing season compared to the observation (Fig. 5).
In the following section, we discuss two primary findings of this study: (1) the importance of representing cohorts at the Wind River oldgrowth forest, rather than using big-leaf models that represent all the vegetation as a single leaf; and (2) the dynamic mechanisms of response to drought simulated by ED2. Finally, we discussed the limitation in current ED2 and potential effort to improve future simulations.

Importance of representing cohorts
This study highlights the utility of cohortbased models for forest stands with complex structure. Although ED2 does not reproduce spatial distribution of individual trees, representing stand composition and structure in terms of cohorts provides a useful tool to explore and evaluate the non-linear responses of different sizes of trees to climate anomalies. At the Wind River old-growth forest stand, we show that there is a clear divergence of microclimate between top and middle-low canopies. Further, we show that Douglas-fir and western hemlock trees respectively occupy upper and middlelower canopy layers and respond differently to the stratified microclimate. These results imply that the stand will not respond as a single bigleaf to climate change; instead, the results imply that the canopy layers occupied by different species will respond significantly differently to climate change. In particular, the top canopy, occupied by cohorts of mature Douglas-fir, are less sensitive to drought events, as demonstrated by observed and simulated responses to drought years (e.g., 2015). The large trees are able extract water from deep soil layers (>60 cm) during drought events, where the SWC remains relatively stable even in extremely dry summers (Warren et al. 2005). In contrast, the mid-canopy cohorts (western hemlock) with medium and small sizes of trees can only utilize soil water from shallower layers (<60 cm), where SWC is highly sensitive to atmospheric drying. Although the medium-size western hemlock cohorts have lower stomatal conductance and GPP rates than the large Douglas-fir cohorts, they contribute the most to ecosystem-level fluxes because of their dominance in tree density. This phenomenon is confirmed by our upscaled sap flow data at this site (Rastogi 2018).
A big-leaf ecosystem model can be and have been parameterized to exhibit the ecosystem dynamics at this site for average site conditions and climate. However, we hypothesize that those models will exhibit significant bias under extreme climate events or disturbance. For example, a big-leaf model may overestimate the decline of stomatal conductance and photosynthetic capacity of all trees at the site during drought events, even though mature Douglas-fir is likely to be resilient to a degree of drought. The difference between big-leaf models and cohort-based models will widen as we consider stand response to a rapidly warming climate in the future. Based on our current modeling results, we hypothesize that previously reported decline of large Douglas-fir trees (Shaw et al. 2005) may underestimate the resilience of Douglas-fir under a rapidly warming climate. Meanwhile, big-leaf models may underestimate the vulnerability of western hemlock that occupy the middle and lower canopies. The resilience of the mature Douglas-fir in the upper canopy may delay succession at the site. The timing of the potential decline of the top canopy (i.e., Douglasfir >45 m) under a rapidly warming climate, the ecosystem response under the dominance of the new dominant species (i.e., western hemlock), and feedback to the climate (e.g., evapotranspiration and albedo) can be more effectively captured by cohort-based model than big-leaf models.

Dynamic mechanisms of responses to drought
At the ecosystem scale, our simulation indicates that the control of soil moisture on carbon fluxes in the simulation is mainly from the GPP response rather than R H response, and the Wind River old-growth forest is simulated to be a partially water-limited ecosystem. The ecosystem respiration rate appears moderately resistant to water stress in the summer, and the seasonal reduction in summer soil moisture does not greatly affect the soil decomposition rates (Fig. 5). This is because deeper soil layers retain sufficient moisture, even though the surface soil dries out during summers. Furthermore, the quality of litter and coarse woody debris constrains decomposition rates in the system (Harmon et al. 2011). The abundant overstory leaves and leaf litter at this site help minimize soil water evaporation. Soil respiration rates are modest in the first place because of low organic matter in deeper soil layers, which partly limits the drought-induced change in soil respiration compared to the photosynthetic response.
Temperature control on productivity in the recent past may be a preview of future climate effects (Mariler et al. 2017). The exact effect of a warmer climate on forest productivity will likely depend on the type, timing, and magnitude of future climate anomalies (Sippel et al. 2016). This is illustrated by a comparison of simulation results from two droughts years: 2013 and 2015. Both 2013 and 2015 were exceptionally dry years, receiving low winter and spring precipitation.
However, the reduction in GPP was much larger in 2015 than in 2013, because the anomalously warm temperatures in spring 2015 depleted significantly more soil water than in 2013. The enhanced rates of photosynthesis and the depletion of soil water in the early spring further amplified the effect of extreme temperatures later in the summer. In particular, the high temperatures and VPD leaf in summer reduced stomatal conductance far more in 2015 than in 2013. While reductions in GPP are the primary driver of reductions in NEP, reductions in NEP may also reflect ED2's over-sensitivity to temperature for calculating R eco when temperatures warm in the spring. This partly explains the lower NEP in the spring of 2015 than in 2013.
A consistent VPD control on stomatal conductance has been observed both at the stand scale (Sulman et al. 2016) and at multiple biomes across the globe . Previous studies note that increases in VPD are always detrimental to transpiration and net carbon uptake, regardless of whether there is a co-occurring drought (Eamus et al. 2013). At Wind River, soil water deficit and high VPD do not always coincide. Duarte et al. (2017) applied a processbased model CLM4.5 at the Wind River forest and report that stomatal closure was primarily driven by VPD, rather than by soil moisture. This is likely because CLM4.5 simulates water limitation only in terms of available soil moisture and does not factor in water demand. Therefore, most of the variations in stomatal conductance displayed by CLM4.5 result from the varying effect of VPD on photosynthesis. In contrast, water limitation in ED2 is driven by both supply and demand. In ED2, when VPD is high, demand is low, so soil moisture limitation is low. Thus, drought may reduce stomatal conductance either by simulating stomatal closure at times of high VPD or by reducing available soil moisture. In our simulations, we observe that during periods of low VPD (<2 kPa), stomatal conductance is mainly determined by SWC. However, in exceptionally dry years (e.g., 2015) when VPD leaf is high (>6 kPa), reduction in stomatal conductance was controlled by VPD regardless of soil moisture status. This ability to resolve the different drivers of reductions in stomatal conductance further highlights the utility of ED2 for simulating climate change impacts in moist forests.

Modeling limitation and potential improvement
The multi-layered canopy structure in ED2 allows substantial light penetration into the lower canopy that facilitates light capture by the leaves. However, the use of top-flat scheme where the model locates all leaf biomass of a cohort at the cohort's height, along with the use of Beer-Lambert law for light distribution through the canopy, results in too much radiation penetrating to the lower canopy. As a consequence, lower-canopy leaves receive too much PAR and are heated excessively, which leads to overestimation of photosynthetic rates and underestimation of intercellular carbon dioxide concentrations (ci/ca ratios). A better scheme for vertical distribution of leaf biomass is therefore needed to improve the simulation of the lower canopy and understory physiology of forest ecosystems with complex canopy structures (Weng et al. 2015). Such model improvements are especially important for understanding the growth and demography of shade-tolerant tree species, such as the western hemlock occurring at Wind River.
Although mature stands are able to utilize deep soil water during drought periods, previous studies have found that foliage at the top of tall trees (e.g., height >50 m) often experiences very low water potentials and enhanced likelihood of xylem embolism because of the long transport distance between root layer and upper foliage layer Yoder 1997, Ryan et al. 2006). Not accounting for this long hydraulic pathway may cause overestimation in stomatal conductance in the upper canopy of Douglas-fir at Wind River, which has hydraulic path lengths that exceed 50 m (Wharton et al. 2009). Recent enhancements to ED2 include a mechanistic representation of water-limited photosynthesis, where the model tracks leaf and stem water potential and used it to solve for root zone water uptake, simulate transport of water vertically through the sapwood, and calculate ultimate transpiration rates. The hydrologically enhanced version has been shown to accurately capture ecosystem dynamics at diverse locations around the globe (Trugman et al. 2016. The complex dynamics of forest stand response to drought simulated by ED2-for example, the differences among the cohorts, the different ways that drought reduces net ecosystem productivity, and the importance of high spring temperatures-suggests earlier generations of dynamic vegetation models applied to the region may oversimplify important controls on plant response to climate change. For example, a rule-based biogeography algorithm (Neilson 1995) used in MC1 (Bachelet et al. 2001) and LPJ-DGVM (Sitch et al. 2003) relies on only minimum temperature for survival to rule out deciduous broadleaf trees in the Pacific Northwest (Fisher et al. 2015). However, these rules are built upon the bio-climate relationships diagnosed from historical quasi-steady-state distributions and may not adequately represent the actual mechanisms that may facilitate the establishment or the expansion of deciduous broadleaf species in the Pacific Northwest under a rapidly warming climate. More comprehensive and mechanistic representations of ecosystem dynamics and stand dynamics, such as provided by ED2, may be more instrumental in exploring consequences of climate change on forest succession in the Pacific Northwest.

CONCLUSIONS
This study provides a robust model evaluation for the ED2 at an approximately 500-yr-old Douglas-fir-western hemlock forest in the Pacific Northwest, USA. Using site-specific parameterization, the model captures carbon, water, and energy fluxes reported by observation studies and indicates that the old-growth forest was a small carbon sink during the 1998-2015 study period. ED2 simulates multiple cohorts that comprise the stand, and each cohort responds differently to the heterogeneous microclimate throughout the vertical canopy layers. Cohorts of mature Douglas-fir dominate the upper canopy, and by having access to soil water in deep layers, Douglas-fir exhibits greater resilience to drought. Cohorts of western hemlock occupy the understory and, taken together, dominate canopy conductance and productivity at the site. However, western hemlock exhibits greater vulnerability to drought due to its susceptibility to the drying of the shallow soil layers. By disentangling the climatic controls on carbon and water fluxes, our simulations show that precipitation-based drought causes only a modest reduction in forest productivity. At low and intermediate levels of VPD leaf , stomatal conductance and photosynthesis are controlled primarily by soil moisture. However, during periods of high VPD (>6 kPa), stomata close and conductance and photosynthesis are greatly reduced regardless of soil moisture status. Furthermore, the timing of warm spring temperatures may help deplete soil moisture, setting the stage for more extreme drought effects in the summer. This study emphasizes the utility and the importance of using the next generation of dynamic vegetation model in the Pacific Northwest, a model that has the ability to represent both fast ecosystem processes and the slower processes that govern stand composition and structure. Enhancements to ED2 and further explorations and applications to various sites in the region are expected to provide useful insights into the mechanisms and consequences of forest response to climate change.