Investigating old‐growth ponderosa pine physiology using tree‐rings, δ13C, δ18O, and a process‐based model

Abstract In dealing with predicted changes in environmental conditions outside those experienced today, forest managers and researchers rely on process‐based models to inform physiological processes and predict future forest growth responses. The carbon and oxygen isotope ratios of tree‐ring cellulose (δ13Ccell, δ18Ocell) reveal long‐term, integrated physiological responses to environmental conditions. We incorporated a submodel of δ18Ocell into the widely used Physiological Principles in Predicting Growth (3‐PG) model for the first time, to complement a recently added δ13Ccell submodel. We parameterized the model using previously reported stand characteristics and long‐term trajectories of tree‐ring growth, δ13Ccell, and δ18Ocell collected from the Metolius AmeriFlux site in central Oregon (upland trees). We then applied the parameterized model to a nearby set of riparian trees to investigate the physiological drivers of differences in observed basal area increment (BAI) and δ13Ccell trajectories between upland and riparian trees. The model showed that greater available soil water and maximum canopy conductance likely explain the greater observed BAI and lower δ13Ccell of riparian trees. Unexpectedly, both observed and simulated δ18Ocell trajectories did not differ between the upland and riparian trees, likely due to similar δ18O of source water isotope composition. The δ18Ocell submodel with a Peclet effect improved model estimates of δ18Ocell because its calculation utilizes 3‐PG growth and allocation processes. Because simulated stand‐level transpiration (E) is used in the δ18O submodel, aspects of leaf‐level anatomy such as the effective path length for transport of water from the xylem to the sites of evaporation could be estimated.


INTRODUCTION
Process-based tree growth models incorporate physiological principles that enable them to be widely applied to diverse species and sites, in contrast to empirical growth and yield models. This improves our understanding of how variable environmental conditions influence forest productivity and stand characteristics (Landsberg 2003). A widely used stand-level process model is Physiological Principles in Predicting Growth (3-PG), developed by Landsberg and Waring (1997) and since modified by numerous other investigators (Xenakis et al. 2008, Gonzalez-Benecke et al. 2014, Wei et al. 2014a, Almeida and Sands 2016, Forrester and Tang 2016, Meyer et al. 2018. The 3-PG model utilizes environmental conditions, stand characteristics, and species-specific physiological and allometric measurements to accurately predict growth and productivity in changing environmental conditions and on diverse forested stands (Coops et al. 1998, 2007, Waring and Gao 2016. 3-PG uses a simple light-use efficiency relationship to estimate carbon assimilation (i.e., gross primary productivity, GPP) and the original model version assumes a constant fraction of GPP (0.47) is allocated to net primary productivity (NPP), which is then partitioned into below-and aboveground biomass. The flexibility and simplicity of 3-PG make it advantageous and practical to use for diverse forest types and management applications while also allowing mechanistic processes to be easily examined.
Because carbon allocation is downstream of carbon assimilation, Wei et al. (2014a) used the carbon isotope composition (d 13 C) of tree rings as a new way to constrain 3-PG at the carbon assimilation step to more accurately represent allocation processes, thus improving 3-PG estimates of stand characteristics. This improvement occurred because the tree-ring d 13 C signal represents a balance between carbon and water fluxes during the assimilation process, enabling these fluxes to be integrated and mutually constrained within the model. The tree-ring d 13 C signal is related to both stomatal conductance (g s ) and photosynthesis (A), which together determine the ratio of intercellular to ambient [CO 2 ], and consequently 13 C discrimination (Farquhar et al. 1989). Thus, tree-ring d 13 C is commonly used to reconstruct gas exchange and physiological responses to environmental conditions (McCarroll andLoader 2004, Gessler et al. 2014).
In conjunction with the d 13 C of tree-ring cellulose (d 13 C cell ), the d 18 O of tree-ring cellulose (d 18 O cell ) has also been used to tease apart the relative contributions of A and g s to the d 13 C cell signal because the d 18 O cell signal primarily reflects water fluxes and processes affecting g s . Relying on a combination of both d 13 C cell and d 18 O cell is known as a "dual isotope approach" and is a useful proxy for leaf-level gas exchange (Scheidegger et al. 2000, Barnard et al. 2012). Together, tree-ring d 13 C cell and d 18 O cell can shed light on physiological responses to environmental conditions such as drought, temperature, relative humidity, fertilization, thinning, and pests (Williams et al. 2010, Brooks and Mitchell 2011, Marias et al. 2014, Saffell et al. 2014, Voelker et al. 2014a, b, Hartl-Meier et al. 2015. Therefore, the combination of a tree-ring d 18 O cell submodel with the existing d 13 C cell submodel (Wei et al. 2014a) can reveal more information about stand carbon dynamics and water use and improve our understanding of mechanisms underlying physiological responses to environmental variability.
However, our understanding of the mechanistic drivers of tree-ring d 18 O cell is incomplete (Gessler et al. 2014. Because 3-PG predicts growth and allocation (e.g., gross primary productivity [GPP], basal area, basal area increment [BAI], diameter at breast height [DBH], height) as modified by soil water balance and climate variations, the incorporation of a d 18 O cell submodel into 3-PG can reveal processes that underlie d 18 O cell predictions, as previous models of d 18 O cell have not explicitly included variation in growth or allocation in their predictions (Barbour andFarquhar 2000, Roden et al. 2000). The 3-PG model also provides monthly estimates of transpiration (E) required to determine the Peclet effect, which describes the mixing of isotopically unenriched water arriving at the evaporative sites via bulk flow driven by E with the isotopically enriched water diffusing back from the evaporative sites (Farquhar and Lloyd 1993). Although the Peclet effect has improved d 18 O cell predictions in multiple systems , Holloway-Phillips et al. 2016, the importance of the Peclet effect for estimating d 18 O has been debated (Ripullone et al. 2008, Xiao et al. 2012, B€ ogelein et al. 2017. This is partly because contributions of each Peclet effect component are challenging to disentangle, especially the effective path length (L), an elusive component of the Peclet effect that is not directly measurable (Song et al. 2013. 3-PG with the d 18 O cell submodel can advance our understanding of factors that drive d 18 O cell dynamics including estimates of L. The aim of this study is to test a further modification of 3-PG that includes both a d 13 C cell and d 18 O cell submodel to understand physiological drivers of stand characteristics using long-term trajectories of tree-ring growth, d 13 C cell , and d 18 O cell from old-growth Pinus ponderosa at the AmeriFlux Metolius site in the Oregon Cascades. This study builds upon previous application of 3-PG to P. ponderosa , Coops et al. 2005, Wei et al. 2014b) and is well suited to evaluate how the d 13 C cell and d 18 O cell submodels in 3-PG can improve our understanding of physiological drivers of tree-ring isotope ratios. For this work, we draw upon a combination of tree-ring growth and stable isotope observations, along with extensive meteorological and physiological measurements recorded at this site , 2001a, Warren et al. 2005. We parameterized the model using an "upland" set of oldgrowth P. ponderosa trees~1 km from the Metolius River located at the US-Me2 AmeriFlux site. We then tested the parameterized model on a nearby "riparian" set of similarly aged old-growth trees closer (~0.015 km) to the Metolius River to examine potential site and physiological factors driving observed differences in BAI and d 13 C cell between these two sites. The goals of this study were to (1) use long-term measured BAI, d 13 C cell , and d 18 O cell trajectories to evaluate the performance of the updated 3-PG model with the d 13 C cell and d 18 O cell submodels; (2) use this first test of 3-PG with a d 18 O cell submodel to improve our understanding of the mechanistic controls of d 18 O cell ; and (3) demonstrate how the model can be used to explore potential site and physiological differences between upland and riparian sets of trees.

Study site
This study was conducted in a P. ponderosa forest on the eastern side of the Cascade Mountains in central Oregon within the Metolius Research Natural Area 44.4957°N,121.6224°W) at an elevation of 915 m. The site consists of 27% old trees (~250 yr), 25% younger trees (~45 yr), and 48% mixed-age trees. Bitterbrush (Purshia tridentata) and bracken fern (Pteridium aquilinum) comprise the sparse understory. Precipitation is greatest between October and June with dry summer months. Winters are wet and cold, snow cover in winter is common, and freezing temperatures occur mostly at night and early morning. Soil is classified as sandy loam (73% sand, 21% silt, and 6% clay) and soil nutrients are low , Warren et al. 2005.

Climate data
AmeriFlux CDIAC climate and eddy flux (ecosystem gross primary productivity, GPP; ecosystem water vapor flux, LE) data for the intermediate ponderosa pine site (US-Me2) were available for 2002-2012. To extrapolate back in time, monthly minimum and maximum temperature (T min and T max , respectively) obtained from PRISM (http://www.prism.oregonstate.edu/) for 1895-2012 was compared with the site-level AmeriFlux data and that relationship (y = 1.27T min_PRISM + 2.46, R 2 = 0.92; y = 0.89T max_PRISM À 1.42, R 2 = 0.98) was used to correct PRISM data, which are based on 4-km grid cells (data available online). 6 Precipitation data were also obtained from PRISM. As PRISM climate data only extended to 1895, we focused on 1895-2002 in this study. Data on atmospheric [CO 2 ] and its d 13 C from Francey et al. (1999) were used in simulations of this period.
Appendix S1: Fig. S1 describes the required driving climate variables for 3-PG: annual average minimum air temperature (T min ), average air temperature (T av ), maximum air temperature (T max ), vapor pressure deficit (VPD), solar radiation, and precipitation from 1895 to 2002 (actual mean monthly climate data were used in the model, not the annual averages depicted in Appendix S1: Fig. S1). Mean monthly VPD was calculated as the difference between saturation vapor pressure at minimum and maximum temperature. Mean daily solar radiation was calculated from mean monthly T min and T max (Bristow and Campbell 1984, Thornton et al. 1997, Coops et al. 1998, 2000 as in Landsberg et al. (2003). The number of frost days per month (F) was calculated based on mean monthly T min : If T min > 6, then F was set to zero. This equation assumes photosynthesis does not occur on days with temperatures below À2°C (Waring 2000).

Tree-ring analyses
Tree cores were collected in early spring 2003 from two sets of trees: an upland set~1 km from the Metolius River at the AmeriFlux Metolius Intermediate Pine site (US-Me2), and a nearby riparian set within 0.015 km of the Metolius River located just north of Camp Sherman. The upland and riparian sites were <5 km apart. We sampled five upland and five riparian P. ponderosa trees of approximately the same stem diameter at 1.3 m height (upland, 87.8 AE 3.9 cm; riparian, 112.3 AE 8.2 cm) and age (mean age % 260 yr; Table 1). In spring 2003 prior to diameter growth, three 12 mm cores from each tree were collected for isotopic analysis, along with a separate 5 mm core that was collected as an archive. Cores were dried and sanded, and the 5 mm core was mounted. All cores were age dated, and ring widths were measured using a sliding stage incremental micrometer (Velmex, Bloomfield, New York, USA) with Measure J2X software (VoorTech Consulting, Holderness, New Hampshire, USA). Visual cross-dating was verified using the COFECHA program to identify false or missing rings (Holmes 1983) for all cores collected. Tree-ring widths were converted to basal area increment (BAI, cm 2 Átree À1 Áyr À1 ) by assuming a circular outline of stem cross-sections.

Sample preparation
The 12 mm cores were separated into annual increments spanning from 2002 to 1850 (152 yr). However, we focused on 1895-2002 in this study because PRISM climate data only went back to 1895. The annual increments from three cores per tree were combined for a single sample per tree per year. Each annual ring was ground with a ball mill to a fine powder. All samples were extracted for alpha-cellulose. Oils and resins were removed with toluene-ethanol and ethanol Soxhlet extractions (Leavitt and Danzer 1993). Holocellulose was isolated by delignification in an acetic acid-acidified sodium chlorite solution and converted to alpha-cellulose in sodium hydroxide (Sternberg 1989).
Approximately 0.8 mg of alpha-cellulose was loaded into tin capsules for C combustion and 0.4 mg into silver capsules for O pyrolysis for subsequent isotopic analysis by isotope ratio mass spectrometer (IRMS) at the Integrated Stable Isotope Research Facility at the Western Ecology Division of the U.S. EPA, Corvallis Oregon, USA. Samples analyzed for 13 C were flash combusted using an elemental analyzer (ECS 4010; Costech, Valencia, California, USA), and the resulting CO 2 analyzed by continuous-flow IRMS (Delta Plus XP, Finnigan MAT, now Fisher Scientific, Waltham, MA, USA). Each run was calibrated using three internal standards (NIST concentration standards of corn, bovine liver, and tomato) spanning the range of expected values, with an independent QC standard (cellulose) to calculate accuracy. Internal standards were routinely calibrated to international standards USGS42 (Tibetan hair), NIST 8542 sucrose, NIST 8573 and 8574 glutamic acid, and NIST 8514 graphite. Typical precision and accuracy was AE0.1& or better as determined by repeated measures of internal quality control standards and from sample replicates. Samples were analyzed for 18 O using a high temperature conversion elemental analyzer (TC/EA ThermoQuest Finnigan, now Fisher Scientific) interfaced to an IRMS (Thermo Electron Delta XL, now Fisher Scientific). Internal laboratory standards (NIST concentration standards of pine needles, sucrose, and corn) were used for calibration standards with an independent QC standard (cellulose) for accuracy estimates. IAEA-601 and IAEA-602 benzoic acid were used to routinely calibrate the internal standards. Typical error was AE0.2& or better as determined by repeated measures of internal quality control standards and from sample replicates. The C and O stable isotope ratios (R) of the heavier (i.e., 13 C, 18 O) to lighter isotope (i.e., 12 C, 16 O) were represented by delta (d) notation in parts per thousand (&) relative to the VPDB or VSMOW international standards (McCarroll and Loader 2004): 3-PG isotope model 3-PG is a simplified model that incorporates essential tree physiological, hydrologic, and growth processes to predict net primary productivity (NPP), biomass allocation, water use, soil water balance, stem mortality (self-thinning), litterfall, and root turnover on a monthly time step (Landsberg and Waring 1997). The input/driving data include mean monthly values of T min , T max , T av , precipitation, F, solar radiation, VPD, atmospheric [CO 2 ], d 13 C of the atmosphere, and d 18 O of source water. Stand initialization values include initial foliage biomass, root biomass, stem biomass, stand density (stocking), available soil water, and maximum available soil water and were based on previously reported values (Appendix S1: Table S1). The stand-level outputs include biomass pools for roots, stems, and foliage, GPP, NPP, transpiration (E), growth and stand characteristics, and d 13 C cell and d 18 O cell in this study. In all figures and tables, outputs are presented as averages for May-September when P. ponderosa physiological activity and radial growth at this site occurs (Fowells 1941).
3-PG is based on the light-use efficiency modeling paradigm that captures a positive linear relationship between plant growth and intercepted radiation. Specifically, 3-PG calculates GPP from absorbed photosynthetically active radiation (/ p.a , mol/m 2 ) and canopy quantum efficiency (a c , mol C/mol photon) and is constrained by factors that influence stomatal closure, including atmospheric VPD, soil moisture, temperature, frost, and site nutrient status: f ϴ , and f age are the temperature, frost, nutrition, VPD, soil water, and age modifiers, respectively, and a cx is the maximum canopy quantum efficiency.
The temperature modifier (f T ) incorporates the minimum, maximum, and optimum temperatures for growth. The frost modifier (f F ) is calculated using F. The nutrient modifier (f N ) is a function of site fertility rating (FR), ranging from 0 to 1, and is based on available soil nutrients. The VPD modifier (f D ) is a function of k g , a species-specific coefficient describing the strength of the response of canopy conductance (g c ) to VPD (Law et al. 2001a). The soil water modifier (f ϴ ) is calculated using the moisture ratio of current : available water and a soil water constant (c ϴ ) and power (r ϴ ) that reflect different soil types . For sandy loam at our study site, c ϴ is 0.4 and r ϴ is 7. The age modifier (f age ) accounts for reductions in hydraulic and stomatal conductance as stands age.
The 3-PG model the ratio of NPP to GPP is fixed at 0.47 . NPP is allocated to foliage, woody tissue, and root biomass pools based on speciesspecific partitioning algorithms, which also depend on site and growth conditions, litterfall, and root turnover . The model uses a simple relationship to determine root growth and turnover to estimate belowground biomass allocation. Allometric ratios are used to determine the allocation of biomass to stems and foliage. Stem growth, stand density, and stem mortality are calculated according to the self-thinning rule based on the negative relationship between tree density and stem mass (Landsberg and Waring 1997). Soil water balance is based on rainfall, irrigation, evapotranspiration, and runoff/drainage. Evapotranspiration (ET, JÁm À2 Ás À1 ) is determined from the Penman-Monteith equation and canopy conductance (Penman 1948, Monteith 1965, Monteith and Unsworth 2007: where e 20 is the ratio of the slope of the saturation-vapor pressure curve at T air = 20°C to the psychrometric constant (2.2); R n is the net radiation (W/m 2 ); q is density of air (1.2 kg/m 3 ); k is the latent heat of vaporization of water (2,460,000 J/kg); g bl is boundary layer conductance, and g c is canopy conductance (both in m/s); and e s À e a is the saturation vapor pressure deficit. E was converted to molÁm À2 Ás À1 to be used in the Peclet calculation in the d 18 O submodel. We updated the calculation of canopy conductance (g c ) by multiplying it by the frost modifier ( f F ) to prevent any E from occurring on days with frost. g c was calculated as: where TK 2 and TK 3 are temperature modifiers (0.244, 0.0368, respectively, (Wei et al. 2014a), g cmax is maximum canopy conductance, LAI is leaf area index, and LAI gcx is the LAI required for a stand to reach its g cmax (3.3 m 2 /m 2 ). The impact of this update is presented in Appendix S1: Fig. S2.

Allometric equation to estimate stem biomass
Diameter at breast height (DBH) and biomass measured in Pinus species (Gholz et al. 1979) were used to determine the stem constant (S c ) and stem power (S p ) used in 3-PG in (Wei et al. 2014b). Live branch mass, stem wood mass, and stem bark mass (Gholz et al. 1979) were summed to calculate total stem biomass. Total stem biomass was then plotted against stem DBH.
The relationship between DBH and total biomass (W) was described by an exponential function: where S c = 0.0273 and S p = 2.6405.
d 13 C cell theory and submodel The d 13 C of photosynthate (d 13 C plant ) is described in  as: where a is the kinetic fractionation effect associated with diffusion of CO 2 through stomata (4.4&), b is the net kinetic fractionation effect (27&) associated with discrimination against 13 C by the enzyme RUBISCO (ribulose bisphosphate carboxylase-oxygenase) during C fixation, and c i /c a is the weighted mean ratio of the intercellular CO 2 concentration (c i ) to that in the ambient air (c a ; Farquhar et al. , 1989. The c i can be estimated from c a , photosynthesis (A), and canopy conductance (g c ; Farquhar and Sharkey 1982): The value of 0.66 describes the ratio of diffusivities of CO 2 to water vapor in air (Wei et al. 2014a). Therefore, tree-ring d 13 C plant reflects factors that influence discrimination against 13 C during photosynthetic CO 2 fixation. These factors include the biochemical capacity to fix CO 2 (A), and the conductance (g) to CO 2 from the atmosphere to the sites of carboxylation. Although leaflevel g includes stomatal conductance (g s ) and mesophyll conductance, we assume the simplified equation from Farquhar et al. ( , 1989, assuming that canopy conductance (g c ) and c i are sufficient to model d 13 C plant (Cernusak et al. 2003).
The d 13 C submodel developed previously for 3-PG (Wei et al. 2014a) treats the canopy as a big leaf (Farquhar et al. 1989) and combines Eq. 8 and Eq. 9 so d 13 C plant is calculated as d 13 C plant % d 13 C air À a À ðb À aÞ 1 À A c a 0:66g : (10) To convert d 13 C plant of new photosynthate to d 13 C of tree-ring wood (d 13 C wood ), a constant offset (e sp ) of 1.99& was assumed (Wei et al. 2014a), similar to that observed in other Pinus species (Gessler et al. 2009, Wei et al. 2014b: This model was modified to include a constant offset (e wc ) of 1.5& observed in P. ponderosa (English et al. 2011) between the d 13 C wood and the d 13 C of tree-ring cellulose (d 13 C cell ):  (Craig andGordon 1965, Dongmann et al. 1974). Under steady state conditions where d 18 O es , d 18 O s , and d 18 O v represent the oxygen isotopic composition of leaf water at the site of evaporation, source water, and atmospheric water vapor, respectively. e a /e i is the ratio of ambient to saturation vapor pressure within the leaf, e* is the equilibrium fractionation between liquid water and vapor, and e k is the kinetic fractionation factor of vapor diffusion from the leaf to the atmosphere; d 18 O v was estimated using T av and was based on the assumption that vapor is in isotopic equilibrium with source water (Majoube 1971 Fig. S3) and selected the d 18 O s that resulted in the greatest Pearson correlation coefficient (R) between modeled and observed d 18 O cell for both upland and riparian sites. Based on the R comparisons, the d 18 O s calculated from the multiple linear regression model was selected for both the upland and riparian sites and was used for the presented results (Appendix S1: Table S2, Fig. S3).
Leaf water d 18 O (d 18 O lw ) heterogeneity can be explained further by the Peclet effect, which describes the ratio between the E-induced mass flow (advection) of unenriched source water to the evaporative sites and the back diffusion of isotopically enriched water from the sites of evaporation Lloyd 1993, Barbour 2007): P ¼ EL CD (15) where d 18 O lw is the steady state isotopic enrichment of mean leaf lamina water, P is the Peclet number describing the ratio of advection to diffusion, E is the leaf transpiration rate (molÁm À2 Ás À1 ), L is the scaled effective path length (m) for water movement from the veins to the site of evaporation, C is the molar density of water (55.56 9 10 3 mol/m 3 ), and D is the diffusivity of the heavy water isotopologue (H 2 18 O) in water (2.66 9 10 À9 m 2 /s). L is defined as the product of two components: l, the actual distance of the water pathway from xylem to the evaporative surface, and k, a scaling factor that accounts for the tortuosity of the path of water through a porous medium Lloyd 1993, Barbour et al. 2000).
Isotopic fractionation during the incorporation of the d 18 O lw signal into cellulose of plant tissue is described by the following equation (Roden et al. 2000): where f O is the proportion of oxygen atoms that exchange with source water during cellulose formation (0.42; Roden et al. 2000), and e o is a fractionation factor of +27& associated with the water/carbonyl interactions (Yakir et al. 1990 Wei et al. 2018). Data S1 includes the Python version of the 3-PG model with the d 13 C and d 18 O submodels. First, parameters were set to defaults used in previous applications of 3-PG at this Metolius site (Appendix S1: Table S1). Next, several parameters were calibrated following the approach of (Wei et al. 2014a) and . Maximum canopy conductance (g cmax ) and the coefficient describing the sensitivity of canopy conductance to VPD (k g ) were calibrated based on previously reported values of E and g for P. ponderosa at this site (Law et al. , 2001a, and the equation describing the relationship between g and VPD (Law et al. 2001a). Fertility rating (FR), foliage : stem partitioning ratio of tree diameter of 20 cm (pfs20), maximum root partitioning (p rx ), maximum tree stem mass likely in mature stands of 1,000 trees/ha (wSx1000), and maximum quantum efficiency (a cx ) were calibrated based on observed and previously reported values of LAI, basal area, BAI, stand density, d 13 C cell , and d 18 O cell with other parameters held constant (Wei et al. 2014a). L was calibrated to match modeled and observed d 18 O cell with a Peclet effect included.
The measured trajectories of d 13 C cell and d 18 O cell at the upland site were prioritized for calibration over previously reported stand characteristics, because we had 107 yr of measured tree-ring data; although the previously reported stand characteristics (Table 1) were measured at the same site, they were only for single years and not necessarily on the same trees for which we had the long-term measured tree-ring trajectories of d 13 C cell , d 18 O cell , and BAI. To do this, calibration of parameters as described above was determined first based on minimizing the root mean squared error (RMSE) between the measured and modeled values of d 13 C cell , d 18 O cell , and then BAI, and then previously reported stand characteristics.
The model trained on the upland site was then tested on the riparian site. Because the upland and riparian sites were <5 km apart, climate driving inputs and stand initialization values remained the same as those for the upland site, and the model was also run for the same time period as the upland site (i.e., from 1895 to 2002). Similar to the upland site, we compared the results of all three methods used to estimate monthly d 18 O s for the riparian site (Appendix S1: Table S2, Fig. S3) and selected the d 18 O s that resulted in the greatest R between modeled and observed d 18 O cell . The trained model did not predict the observed trajectories of BAI and d 13 C cell at the riparian site. Thus, we adjusted maximum available soil water (ASW), maximum canopy conductance (g cmax ), and wSx1000 to minimize the RMSE between the measured and simulated BAI, d 13 C cell , and d 18 O cell of riparian trees. By adjusting these parameters, we demonstrated how model parameterization can be used to identify site and physiological differences between upland and riparian trees. Other parameters that increased the RMSE between the measured and simulated BAI, d 13 C cell , and d 18 O cell of riparian trees were not adjusted.

Sensitivity analysis
We investigated the sensitivity of modeled BAI, d 13 C cell , d 18 O cell with Peclet, d 18 O es , E, g c , GPP, and LAI to AE20% and AE40% changes in the following parameters: a cx , FR, g cmax , k g , maximum ASW, pfs20, p rx , and wSx1000. These parameters were selected because they are known to influence d 13 C cell and biomass allocation (Wei et al. 2014a), but it is unknown how they might influence d 18 O cell . L was adjusted to evaluate the sensitivity of d 18 O cell with Peclet to changes in L. To conduct the sensitivity analysis, one parameter at a time was adjusted AE20% and AE40% of the original value while all other parameters were held constant. The % change in the aforementioned output response (BAI, d 13 C cell , d 18 O cell with Peclet, d 18 O es , E, g c , GPP, and LAI) was then quantified from the original output value. These % changes were averaged over the study period of 1895-2002. An output was considered "sensitive" if a change in parameter resulted in a ≥10% change in output.

Statistics
Modeled and observed values were compared using simple linear regression, Pearson correlation coefficient (R), the coefficient of determination (r 2 ), and root mean squared error (RMSE) in SigmaPlot 13.0 (Systat Software, San Jose, California, USA). A one-sample t test was used to test if mean difference for each year spanning 1895-2002 between observed and modeled BAI, d 13 C cell , and d 18 O cell , between upland and riparian sites, and between d 18 O cell with and without the Peclet effect was significantly different from zero.

RESULTS
Parameterizing 3-PG by prioritizing fit to measured BAI, d 13 C cell , and d 18 O cell , and secondarily previously reported stand characteristics, allowed the model to  predict stand characteristics at the upland site reasonably well (Figs. 1, 2). Predicted BAI was within range of observed BAI values (Fig. 3a) where the mean difference between observed and modeled BAI for each year spanning 1895-2002 was not significantly different from zero (P = 0.67; observed BAI = 30.3 AE 0.9 cm 2 /yr, modeled BAI = 30.8 AE 1.0 cm 2 /yr). The model also predicted d 13 C cell within range of observed d 13 C cell at the upland site (Fig. 3c)  Compared to upland trees, observed BAI of riparian trees was consistently and significantly greater (Fig. 3a,b) where the 1895-2002 mean observed BAI of riparian trees was 60.0 AE 1.9 cm 2 /yr, nearly twice that of upland trees, and the mean difference between upland and riparian observed BAI for each year spanning 1895-2002 was significantly different from zero (P < 0.001). Riparian trees also exhibited consistently and significantly lower (more negative) mean observed d 13 C cell than upland trees (Fig. 3c, d); the 1895-2002 mean observed d 13 C cell of riparian trees was À23.8 AE 0.04& which was 0.7& more negative than upland trees, and the mean difference between upland and riparian observed d 13 C cell for each year spanning 1895-2002 was significantly different from zero (P < 0.001). The observed d 13 C cell of the upland trees ranged from À24.8& to À22.3&, while observed d 13 C cell of the riparian trees ranged from À25.3& to À23.1&. In contrast to BAI and d 13 C cell , the mean difference between upland and riparian observed d  Table 2). The observed d 18 O of stem water and atmospheric water vapor were also similar between upland and riparian trees (Table 2). Because observed d 18 O data for river water, stem water, and water vapor were only available for one value per sampling date, we could not test for statistically significant differences.
The model trained on the upland site did not predict the observed BAI and d 13 C cell trajectories of riparian trees given their significantly greater BAI and lower d 13 C cell trajectories compared to upland trees (Fig. 3). However, we hypothesized that riparian trees had access to a greater water supply, which would increase canopy conductance and thus influence BAI and d 13 C cell . To test this hypothesis, we increased maximum available soil water (ASW) from 163 to 300 mm, maximum canopy conductance (g cmax ) from 0.014 to 0.0145 m/s, and decreased the maximum tree stem mass likely in mature stands of 1,000 trees/ha(wSx1000) from 110 to 100 kg/ tree. These parameter changes for the riparian trees reproduced the greater BAI and lower d 13 C cell values observed in riparian trees compared to upland trees without substantially influencing the modeled d 18 O cell of riparian trees (Fig. 3). The mean differences between riparian observed and modeled BAI, d 13 C cell , and d 18 O cell for each year spanning 1895-2002 were not significantly different from zero (P = 0.11, 0.20, 0.12, respectively; Fig. 3). By updating maximum ASW, g cmax , and wSx1000, we demonstrate how model parameterization can be used to identify site and physiological differences between upland and riparian trees.  Ecology,Vol. 100,No. 6 Although the mean differences between modeled and observed BAI for 1895-2002 at both the upland and riparian sites were not significantly different from zero (P > 0.05), modeled and observed BAI were not significantly correlated (simple linear regression, P > 0.05) and RMSE was 12.7 and 25.3, respectively (Fig. 5). In contrast, modeled and observed d 13 C cell at both the upland and riparian sites were significantly correlated (simple linear regression, P < 0.001, R = 0.64, 0.52, r 2 = 0.41, 0.27, respectively) and RMSE was 0.56 and 0.46, respectively (Fig. 5). Modeled and observed d 18 O cell with Peclet at both the upland and riparian sites were also significantly correlated (simple linear regression, P < 0.001, R = 0.47, 0.49, r 2 = 0.22, 0.24, respectively) and RMSE was 0.88 and 0.91, respectively (Fig. 5).
The sensitivity analysis showed that BAI, E, g c , GPP, and LAI were sensitive (i.e., each variable was changed by ≥10% of its original value) to changes in each of the tested parameters: a cx , FR, k g , pfs20, and p rx (Appendix S1: Table S3; Fig. 6). In contrast, d 13 C cell was only sensitive to changes in a cx , g cmax , and p rx . d 18 O es was not sensitive to any parameter, while d 18 O cell with Peclet was sensitive to a cx , FR, k g , pfs20, and p rx .
L was estimated using d 18 O cell with Peclet, where an L value of 0.010 m resulted in modeled d 18 O cell that best predicted the observed d 18 O cell (Fig. 7).

DISCUSSION
We used long-term tree-ring growth, d 13 C cell , and d 18 O cell to demonstrate for the first time the use of 3-PG with the newly added d 18 O submodel and to investigate physiological differences between old-growth upland and riparian P. ponderosa. This application of 3-PG in conjunction with measured BAI and the dual isotope approach across 107 yr is unprecedented. Model performance was better for both d 13 C cell and d 18 O cell than for BAI and stand growth characteristics (Figs. 1-3, 5). The d 18 O submodel with the Peclet effect improved predictions of d 18 O cell (Fig. 4; Appendix S1: Table S2) because it incorporates 3-PG's monthly E predictions, and also reflects growth and allocation processes, unlike previous d 18 O cell models based just on leaf water and cellulose isotope equations driven by precipitation and vapor isotopes and relative humidity. Our approach also provided a way to estimate L (Fig. 7), which is typically an unmeasurable component of the Peclet effect. By using the model to aid our understanding of the physiology driving the BAI, d 13 C cell , and d 18 O cell trajectories at the upland and riparian sites, we propose that upland and riparian trees were using the same source of water isotopically, but that greater water availability at the riparian site increased tree growth and lowered d 13 C cell values.

Model performance
Model parameterization based on previously reported parameters and stand characteristics for P. ponderosa resulted in overall reasonable values of stand characteristics with minor discrepancies (Fig. 1). Modeled height and LAI fell within the range of the mean of previously reported values. Modeled LAI reached~2.2 m 2 /m at a stand age of 30 yr (not shown), consistent with observed values for P. ponderosa at that same age (Law et al. 2001b), indicating the LAI was well predicted beyond our study period of 1895-2002 when the trees ranged from 155 to 262 yr in age. NPP was within range of previously reported NPP for P. ponderosa at the Metolius site . Modeled E and GPP were only slightly lower than the observed ecosystem water vapor flux (LE) and GPP AmeriFlux values (Fig. 2), but LE includes all forms of evaporation from the soil and understory, which is likely why observed LE was slightly greater than simulated E (i.e., transpiration only).
Minor discrepancies existed between modeled and previously reported stand growth characteristics. These discrepancies arose due to our prioritizing model calibration first on measured ring width and isotope values, and secondarily on previously reported stand-level values, which differed from our study trees. Previously reported stand densities ranged from 54 to 137 trees/ha (Table 1) and previously reported diameter at breast height (DBH) ranged from 55 to 63 cm (Law et al. 2001a, b, Youngblood et al. 2004, Warren et al. 2005) at this site. However, the trees cored in this study had DBH values ranging 88-112 cm (Table 1), as we preferentially selected old and large dominant trees. Thus, the model predicted fewer, larger diameter trees (Fig. 1), as 3-PG utilizes the established negative relationship between DBH and stand density (Meyer 1938). Because 3-PG was developed for even-aged stands, modeled basal area was greater than the basal area previously reported for younger or mixed-aged stands at the site (Fig. 1a). However, discrepancies between modeled and previously reported stand density and basal area are not surprising, as we prioritized tree-ring growth and isotope trajectories for model parameterization. This calibration method is consistent with Landsberg et al. (2003).
Although modeled BAI values were within range of observed values (Fig. 3a, b), the ability of 3-PG to predict variation in BAI over time was sub-optimal (i.e., larger RMSE; Fig. 5a, b). 3-PG was developed to simulate even-aged plantation stands, and thus stand growth is often challenging to model accurately in unmanaged forests (e.g., modeled DBH in Wei et al. 2014a). This is in  Ecology,Vol. 100,No. 6 part because microclimatic and site-specific conditions such as non-uniform tree spacing are not accounted for in a generalizable forest stand growth model like 3-PG ). This sub-optimal predictive ability of the model is also likely due to the relatively long period 107 yr) over which we applied the model compared to previous studies. Waring and Gao (2016) compared the tree-ring index of Picea crassifolia to normalized 3-PG-predicted diameter growth for 52 yr. Wei et al. (2014a) compared observed and predicted DBH of Abies grandis for 50 yr with sub-optimal success. The 107-yr timeframe distinguishes this study from previous tree-ring research using 3-PG; however, it is likely that 3-PG may more accurately predict observed growth across a shorter timeframe , Wei et al. 2014b). These latter studies applied 3-PG to P. ponderosa only for 25 and 2 yr, respectively. Although 3-PG's sub-optimal ability to predict BAI may potentially limit its capacity to model tree-ring isotopes, the model still simulated d 13 C cell and d 18 O cell well in this study (i.e.,small RMSE,Figs. 3,5). This is likely because the model predicted monthly E and GPP relatively well (Fig. 2). This suggests that 3-PG with the d 13 C and d 18 O submodels may be used to predict tree-ring isotopes over long timeframes if the model can reasonably predict at least one metric of productivity (e.g. GPP, BAI, basal area, DBH, height).  6. Sensitivity analysis results of the effect of a 40% increase in parameters: maximum quantum efficiency (a cx ), fertility rating (FR), maximum canopy conductance (g cmax ), sensitivity of canopy conductance to VPD (k g ), maximum available soil water (maxASW), foliage : stem partitioning ratio for tree diameter of 20 cm (pfs20), maximum root partitioning (p rx ), and maximum tree stem mass in mature stands of 1,000 trees/ha (wSx1000) on the percent change in output variables: basal area increment (BAI), d 13 C cell , d 18 O cell with Peclet, transpiration (E), canopy conductance (g c ), gross primary productivity (GPP), and leaf area index (LAI). Appendix S1: Table S3 shows all sensitivity analysis values including the effect of AE20% and AE40% shifts in parameters on percent change in output variables. and g cmax , before allocation processes, as indicated by the sensitivity analysis ( Fig. 6; Appendix S1:  Fig. 4; Appendix S1: Table S2). This also allowed L, the unmeasurable component of the Peclet effect, to be estimated (Fig. 7) and our estimate fell within the range of previously reported values for Pinus species (Song et al. 2013). This provides support for the use of a well-parameterized model based on measured and previously reported stand characteristics to estimate L and evaluate the impacts of variation in leaf-level physiology on the stand scale, although this needs to be tested further.

Comparing upland and riparian tree physiology
The model helped to explain the physiological mechanisms underlying the significant differences in BAI and d 13 C cell between upland and riparian trees. We hypothesized that close proximity of the riparian trees to the Metolius River would increase water availability and reduce drought stress compared to the upland trees, thereby altering canopy conductance, BAI, and d 13 C cell . By adjusting maximum ASW, g cmax , and wSx100 to better predict the greater BAI and lower d 13 C cell observed at the riparian site (Figs. 3, 5), we demonstrate how model parameterization can be used to investigate site and physiological differences between upland and riparian trees. This suggests that our hypothesis for the greater BAI and lower d 13 C cell was correct (Orwig and Abrams 1997, Adams and Kolb 2004). Because trees modulate stomatal conductance with water availability to maintain hydraulic function, the greater water availability also allowed riparian trees to maintain hydraulic function, gas exchange, and growth throughout more of the growing season (Panek and Goldstein 2001), consistent with the increase in g cmax for riparian trees. This resulted in greater 13 C discrimination due to lower relative stomatal constraints on A and thus greater CO 2 supply, imparting a lower d 13 C cell signal in tree rings compared to upland trees (McCarroll and Loader 2004). Because g cmax is related to hydraulic properties such as soil-to-leaf hydraulic conductance, greater g cmax is consistent with greater sapwood-specific native conductivity observed in P. ponderosa at a riparian site compared to an upland slope site (Stout and Sala 2003). The adjusted g cmax value of 0.0145 m/s is within range for this species and site , 2001a, Coops et al. 2005.
Given the greater measured growth and lower d 13 C cell of riparian trees and differing access to river water, we also expected d 18 O of the riparian trees to differ from that of the upland trees due to differences in source water and/or in leaf-level physiology. We had also expected that the d 18 O of source water (d 18 O s ) would be constant over time for riparian trees (Manga 2001), and their d 18 O cell variance would only be related to climate (e.g., relative humidity, VPD) and leaf-level physiology, while the d 18 O cell variance of the upland trees would be related to all three drivers. However, observed d 18 O cell did not substantially differ between the upland and riparian trees (Figs. 3, 5). To examine the mechanisms underlying the unexpectedly similar d 18 O cell trajectories between upland and riparian sets of trees, we discuss potential drivers of patterns in d 18 O cell : relative humidity, leaf-level physiology, and d 18 O s (Farquhar et al. 2007, Saugier et al. 2012) in terms of model performance and parameters.
First, differences in relative humidity and VPD can alter tree-ring d 18 O cell (Kahmen et al. 2011, Voelker et al. 2014a) as shown by a significant correlation between VPD and d 18 O cell (Pearson correlation coefficient R = 0.20, P = 0.04, data not shown). However, we expected that the upland and riparian trees experienced similar evaporative demand because they were <5 km from each other. The similarity of the d 18 O cell time series and high correlation between upland and riparian d 18 O cell trajectories (P < 0.001, R = 0.73) also is consistent with this assumption.
Second, leaf-level physiology may contribute to d 18 O cell patterns (Flanagan and Ehleringer 1991). The greater modeled g cmax , greater BAI, and lower d 13 C cell of the riparian trees suggested that riparian trees may have different leaf-level gas exchange compared to the upland trees. However, d 18 O cell was not considered sensitive to changes in g cmax in the model ( Fig. 6; Appendix S1: Table S3) and unexpectedly the d 18 O cell trajectories did not differ between sets of trees.
Finally, d 18 O s is considered to be a primary driver of d 18 O cell . The difference in proximity of the upland and riparian trees to the Metolius River Article e02656; page 14 DANIELLE E. M. ULRICH ET AL. Ecology,Vol. 100,No. 6 suggested that they may use different sources of water (Marshall and Monserud 2006) where riparian trees may rely primarily on river water and upland trees may rely on both river water and precipitation (Stout andSala 2003, Kerhoulas et al. 2013). This is why we expected that riparian trees would have a constant d 18 O s while upland trees would not. However, several lines of evidence suggest that upland and riparian trees were not using different sources of water: (1) observed d 18 O cell from the upland and riparian sites were both most correlated with modeled d 18 O cell with Peclet calculated using d 18 O s determined from the same method (multiple linear regression model, Appendix S1: Table S2); (2) the observed d 18 O signal of the tree-ring cellulose, stem water, and atmospheric water vapor were similar at both sites ( We found that a decrease in temperature, an increase in RH, an increase in L, an increase in E, and a decrease in d 18 O s can lower d 18 O cell (Appendix S1: Fig. S4). The d 18 O cell decline may have been caused by any combination of these factors unique to our study site. The presence of this d 18 O cell decline imprinted in the growth rings of trees from both sites provides more support that the upland and riparian sites and trees were likely using similar sources of water. However, the greater BAI, lower d 13 C cell , and increased maximum ASW and g cmax strongly suggested that the riparian trees had greater access to the same source of water (i.e., greater water availability) compared to upland trees.

CONCLUSIONS
We tested the 3-PG model with the d 13 C cell submodel and the newly added d 18 O cell submodel using long-term trajectories of measured growth, d 13 C cell , and d 18 O cell of old-growth P. ponderosa in central Oregon. The unprecedented use of a 107-yr period for which we had growth and isotope measurements revealed the model's strength in predicting d 13 C cell and d 18 O cell reasonably well across long timeframes but also highlighted the model's limitations in predicting certain stand growth characteristics. This first test of 3-PG with a d 18 O cell submodel improves our understanding of mechanistic drivers of d 18 O cell . Because d 18 O cell with the Peclet effect is calculated using stand-level E output predicted by 3-PG, the Peclet effect improved estimates of d 18 O cell and demonstrated that L and leaf-level physiology may be estimated using a well-parameterized model. The model helped to explain physiological drivers underlying the tree-ring growth, d 13 C cell , and d 18 O cell trajectories measured on the upland and riparian trees. The application of 3-PG with the d 13 C cell and d 18 O cell submodels to the upland and riparian sets of trees indicates the potential of such coupled models to be parameterized for diverse stands using site-and stand-specific information for examining the physiological mechanisms underlying forest responses to changes in climate.

ACKNOWLEDGMENTS
The authors would like to thank Dick Waring for help with 3-PG and for commenting on an earlier draft, Liang Wei, Carlos Gonzalez-Benecke, and Trent Seager for helpful discussions about 3-PG, Sean Hammond and Daniel Griffith for assistance with Python, and Linlin Gao and Steve Voelker for advice on tree-ring analyses. This work was supported by the NSF Graduate Research Fellowship Program and NSF grant IOS 11-46746. This manuscript has been subjected to Agency review and has been approved for publication. The views expressed in this paper are those of the author(s) and do not necessarily reflect the views or policies of the U.S. Environmental Protection Agency. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.