# Leveraging 35 years of Pinus taeda research in the southeastern US to constrain forest carbon cycle predictions: regional data assimilation using ecosystem experiments

R. Quinn Thomas, Evan B. Brooks, Annika L. Jersild, Eric J. Ward, Randolph H. Wynne, Timothy J. Albaugh, Heather Dinon-Aldridge, Harold E. Burkhart, Jean-Christophe Domec, Thomas R. Fox, et al.

► To cite this version:

R. Quinn Thomas, Evan B. Brooks, Annika L. Jersild, Eric J. Ward, Randolph H. Wynne, et al. Leveraging 35 years of Pinus taeda research in the southeastern US to constrain forest carbon cycle predictions: regional data assimilation using ecosystem experiments. Biogeosciences, 2017, 14 (14), pp.3525-3547. 10.5194/bg-14-3525-2017. hal-01600941

HAL Id: hal-01600941

https://hal.science/hal-01600941v1

Submitted on 27 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Leveraging 35 years of *Pinus taeda* research in the southeastern US to constrain forest carbon cycle predictions: regional data assimilation using ecosystem experiments

R. Quinn Thomas<fnref n="1" />, Evan B. Brooks<fnref n="1" />, Annika L. Jersild<fnref n="1" />, Eric J. Ward<fnref n="2" />, Randolph H. Wynne<fnref n="1" />, Timothy J. Albaugh<fnref n="1" />, Heather Dinon-Aldridge<fnref n="3" />, Harold E. Burkhart<fnref n="1" />, Jean-Christophe Domec<fnref n="4" />,<fnref n="5" />, Thomas R. Fox<fnref n="1" />, Carlos A. Gonzalez-Benecke<fnref n="6" />, Timothy A. Martin<fnref n="7" />, Asko Noormets<fnref n="8" />,ª, David A. Sampson<fnref n="9" />, and Robert O. Teskey<fnref n="10" />

*Correspondence to:* R. Quinn Thomas (rqthomas@vt.edu)

Received: 14 February 2017 – Discussion started: 16 February 2017

Revised: 22 May 2017 – Accepted: 19 June 2017 – Published: 26 July 2017

**Abstract.** Predicting how forest carbon cycling will change in response to climate change and management depends on the collective knowledge from measurements across environmental gradients, ecosystem manipulations of global change factors, and mathematical models. Formally integrating these sources of knowledge through data assimilation, or model-data fusion, allows the use of past observations to constrain model parameters and estimate prediction uncertainty. Data assimilation (DA) focused on the regional scale has the opportunity to integrate data from both environmental gradients and experimental studies to constrain model parameters. Here, we introduce a hierarchical Bayesian DA approach (Data Assimilation to Predict Productivity for Ecosystems and Regions, DAPPER) that uses observations of carbon stocks, carbon fluxes, water fluxes, and vegetation dynamics from loblolly pine plantation ecosystems across the southeastern US to constrain parameters in a modified version of the Physiological Principles Predict-

ing Growth (3-PG) forest growth model. The observations included major experiments that manipulated atmospheric carbon dioxide (CO₂) concentration, water, and nutrients, along with nonexperimental surveys that spanned environmental gradients across an 8.6 × 10⁵ km<fnref n="2" /> region. We optimized regionally representative posterior distributions for model parameters, which dependably predicted data from plots withheld from the data assimilation. While the mean bias in predictions of nutrient fertilization experiments, irrigation experiments, and CO₂ enrichment experiments was low, future work needs to focus modifications to model structures that decrease the bias in predictions of drought experiments. Predictions of how growth responded to elevated CO₂ strongly depended on whether ecosystem experiments were assimilated and whether the assimilated field plots in the CO₂ study were allowed to have different mortality parameters than the other field plots in the region. We present predictions of stem biomass productivity under elevated CO₂, decreased precip-

<footnote n="1">¹Department of Forest Resources and Environmental Conservation, Virginia Tech, Blacksburg, VA, USA</footnote>

<footnote n="2">²Climate Change Science Institute and Environmental Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA</footnote>

<footnote n="3">³State Climate Office of North Carolina, North Carolina State University, Raleigh, NC, USA</footnote>

<footnote n="4">⁴Bordeaux Sciences Agro, UMR 1391 INRA-ISPA, Gradignan CEDEX, France</footnote>

<footnote n="5">⁵Nicholas School of the Environment, Duke University, Durham, NC, USA</footnote>

<footnote n="6">⁶Department of Forest Engineering, Resources and Management, Oregon State University, Corvallis, OR, USA</footnote>

<footnote n="7">⁷School of Forest Resources and Conservation, University of Florida, Gainesville, FL, USA</footnote>

<footnote n="8">⁸Department of Forestry and Environmental Resources, North Carolina State University, Raleigh, NC, USA</footnote>

<footnote n="9">⁹Decision Center for a Desert City, Arizona State University, Tempe, AZ, USA</footnote>

<footnote n="10">¹⁰Warnell School of Forestry and Natural Resources, University of Georgia, Athens, Athens, GA, USA</footnote>

<footnote>ªcurrent address: Department of Ecosystem Science and Management, Texas A&M University, College Station, TX, USA</footnote>
itation, and increased nutrient availability that include esti-
mates of uncertainty for the southeastern US. Overall, we
(1) demonstrated how three decades of research in southeast-
ern US planted pine forests can be used to develop DA tech-
niques that use multiple locations, multiple data streams, and
multiple ecosystem experiment types to optimize parameters
and (2) developed a tool for the development of future pre-
dictions of forest productivity for natural resource managers
that leverage a rich dataset of integrated ecosystem observa-
tions across a region.

# 1 Introduction

Forest ecosystems absorb and store a large fraction of an-
thropogenic carbon dioxide (CO2) emissions (Le Quéré et
al., 2015; Pan et al., 2011) and supply wood products to a
growing human population (Shvidenko et al., 2005). There-
fore, predicting future carbon sequestration and timber sup-
ply is critical for adapting forest management practices to
future environmental conditions and for using forests to as-
sist with the reduction in atmospheric CO2 concentrations.
The key sources of information for developing these predic-
tions are results from global change ecosystem manipulation
experiments, observations of forest dynamics across environ-
mental gradients, and process-based ecosystem models. The
challenge is integrating these three sources into a common
framework for creating probabilistic predictions that provide
information on both the expected future state of the forest
and the probability distribution of those future states.

Data assimilation (DA), or data-model fusion, is an in-
creasingly used framework for integrating ecosystem obser-
vations into ecosystem models (Luo et al., 2011; Niu et al.,
2014; Williams et al., 2005). DA integrates observations with
ecosystem models through statistical, often Bayesian, meth-
ods that can generate probability distributions for ecosys-
tem model parameters and initial states. DA allows for the
explicit accounting of observational uncertainty (Keenan et
al., 2011), the incorporation of multiple types of observa-
tions with different timescales of collection (MacBean et
al., 2016; Richardson et al., 2010), and the representation of
prior knowledge through informed parameter prior distribu-
tions or specific relationships among parameters (Bloom and
Williams, 2015).

Using DA to parameterize ecosystem models with obser-
vations from multiple locations that leverage ecosystem ma-
nipulation experiments and environmental gradients will al-
low for predictions to be consistent with the rich history of
global change research in forest ecosystems. Ecosystem ma-
nipulation experiments provide a controlled environment in
which data collected can be used to describe how forests ac-
climate and operate under altered environmental conditions
(Medlyn et al., 2015) and can potentially allow for the op-
timization of model parameters associated with the altered

environmental factor in the experiment. Furthermore, the as-
similation of data from ecosystem manipulation experiments
may increase parameter identifiability (reducing equifinality;
Luo et al., 2009), where two parameters have compensating
controls on the same processes, by isolating the response to
a manipulated driver. Observations that span environmental
gradients include measures of forest ecosystem stocks and
fluxes across a range of climatic conditions, nutrient avail-
abilities, and soil water dynamics. These studies leverage
time and space to quantify the sensitivity of forest dynam-
ics to environmental variation. However, covariation of en-
vironmental variation can pose challenges separating the re-
sponses to individual environmental factors. Overall, assimi-
lating observations from a region that includes environmental
gradients and manipulation experiments is a useful extension
of prior DA research focused on DA at a single site with mul-
tiple types of observations (Keenan et al., 2012; Richardson
et al., 2010; Weng and Luo, 2011).

Southeastern US planted pine forests are ideal ecosystems for exploring the application of DA to carbon cycle and forest production predictions. These ecosystems are dominated by loblolly pine (Pinus taeda L.), thus allowing for a single parameter set to be applicable to a large region containing many soil types and climatic gradients. Loblolly pine represents more than one half of the standing pine volume in the southern United States (11.7 million ha) and is by far the single most commercially important forest tree species for the region, with more than 1 billion seedlings planted annually (Fox et al., 2007; McKeand et al., 2003). There is also a rich history of experimental research located across the region focused on global change factors that have included nutrient addition (Albaugh et al., 2016; Carlson et al., 2014; Raymond et al., 2016), water exclusion (Bartkowiak et al., 2015; Tang et al., 2004; Ward et al., 2015; Will et al., 2015), and water addition experiments (Albaugh et al., 2004; Allen et al., 2005; Samuelson et al., 2008). The region also includes a multiyear ecosystem CO2 enrichment study (McCarthy et al., 2010). Furthermore, many of these experiments are multi-factor with water exclusion by nutrient addition (Will et al., 2015), water addition by nutrient addition (Albaugh et al., 2004; Allen et al., 2005; Samuelson et al., 2008), and CO2 by nutrient addition treatments (McCarthy et al., 2010; Oren et al., 2001). Beyond experimental treatments, southeastern US loblolly pine ecosystems include at least two eddy-covariance sites with high-frequency measurements of C and water fluxes along with biometric observations over many years (Noormets et al., 2010; Novick et al., 2015) and sites with multiyear sap flow data (Ewers et al., 2001; Gonzalez-Benecke and Martin, 2010; Phillips and Oren, 2001). Finally, there are studies that include plots that span the regional environmental gradients and extend back to the 1980s (Burkhart et al., 1985). Overall, the multi-decadal availability of observations of C stocks (or biomass), leaf area index (LAI), C fluxes, water fluxes, and vegetation dynamics in plots with experimental manipulation and plots
across environmental gradients, is well suited to potentially constrain model parameters and predictions of how carbon cycling responds to environmental change.

Using loblolly pine plantations across the southeastern US as a focal application, our objectives were to (1) develop and evaluate a new DA approach that integrates diverse data from multiple locations and experimental treatments with an ecosystem model to estimate the probability distribution of model parameters, (2) examine how the predictive capacity and optimized parameters differ between an assimilation approach that only uses environmental gradients and an assimilation approach that uses both environmental gradients and ecosystem manipulations, and (3) demonstrate the capacity of the DA approach to predict, with uncertainty, regional forest dynamics by simulating how forest productivity responds to drought, nutrient fertilization, and elevated atmospheric CO₂ across the southeastern US.

## 2 Methods

### 2.1 Observations

We used 13 different data streams from 294 plots at 187 unique locations spread across the native range of loblolly pine trees to constrain model parameters (Table 1; Fig. 1). The data streams covered the period between 1981 and 2015. The Forest Modeling Research Cooperative (FMRC) thinning Study provides the largest number of plots that span the region (Burkhart et al., 1985). In this study, we only used the control plots that were not thinned. The Forest Productivity Cooperative (FPC) Region-wide 18 (RW18) study included control and nutrient fertilization addition plots that span the region (134.4 kg ha⁻¹ N + 13.44 kg ha⁻¹ P biannually) (Albaugh et al., 2015). The Pine Integrated Network: Education, Mitigation, and Adaptation Project (PINEMAP) study included four locations dispersed across the region that included a replicated factorial experiment with control, nutrient fertilization (224 kg ha⁻¹ N + 27 kg ha⁻¹ P + micronutrients once at project initiation), throughfall reduction (30% reduction), and fertilization by throughfall treatments (Will et al., 2015). The Southeast Tree Research and Education Site (SETRES) study was located at a single location and included replicated control, irrigation (~650 mm of added water per year), nutrient fertilization (~100 kg N ha⁻¹ + 17 kg P ha⁻¹ with micronutrients applied annually with absolute amount depending on foliar nutrient ratios), and fertilization by irrigation treatments (Albaugh et al., 2004). The Waycross study was a single site with a non-replicated fertilization treatment. The annual application of nutrient fertilization was focused on satisfying the nutrient demand by the trees and resulted in one of the most productive stands in the region (Bryars et al., 2013). These five studies included data streams of stand stem biomass (defined as the sum of stem wood, stem bark, and branches) and live stem density. Waycross and SETRES

included LAI measurements from litterfall traps (Waycross) or estimates from LI-COR LAI-2000 (SETRES). SETRES also included fine root and coarse root measurements. In the PINEMAP, SETRES, and RW18 studies we only used foliage biomass estimates from the control plots. We excluded the foliage biomass estimates from the treatment plots because they were derived from allometric models that may not have captured changes in allometry due to the experimental treatment. We did use LAI measurements from both control and treatment plots where available (SETRES).

We also included observations from the Duke Free-Air Carbon Enrichment (FACE) study where the atmospheric CO₂ was increased by 200 ppm above ambient concentrations. Based on the data presented in McCarthy et al. (2010), the study included six control plots, four CO₂ fumigated rings (including the unfertilized half of the prototype), two nitrogen fertilization treatments (115 kg N ha⁻¹ yr⁻¹ applied annually), and one CO₂ by nitrogen addition treatment (fertilized half of prototype). The Duke FACE study included observations of stem biomass (loblolly pine and hardwood), coarse root biomass (loblolly pine and hardwood), fine root biomass (combined loblolly pine and hardwood), stem density (loblolly pine only), leaf turnover (combined loblolly pine and hardwood), fine root production (combined loblolly pine and hardwood), and monthly LAI (loblolly pine and hardwood).

Finally, we included two AmeriFlux sites with eddy-covariance towers in loblolly pine stands. The US-DK3 site was located in the same forest as the Duke FACE site described above (Novick et al., 2015). The US-NC2 site was located in coastal North Carolina (Noormets et al., 2010). We used monthly gross ecosystem production (GEP; modeled gross primary productivity from net ecosystem exchange measured at an eddy-covariance tower) and evapotranspiration (ET) estimates from the sites. The monthly GEP and ET were gap-filled by the site principal investigator. The GEP was a flux-partitioned product created by the site principal investigator. The biometric data from the US-DK3 site were assumed to be the same as the first control ring. The biometric data from the US-NC2 site included observations of stem biomass (loblolly pine and hardwood), coarse root biomass (loblolly pine and hardwood), fine root biomass (combined loblolly pine and hardwood), stem density (loblolly pine only), leaf turnover (combined loblolly pine and hardwood), and fine root production (combined loblolly pine and hardwood).

### 2.2 Ecosystem model

We used a modified version of the Physiological Principles Predicting Growth (3-PG) model to simulate vegetation dynamics in loblolly pine stands (Bryars et al., 2013; Gonzalez-Benecke et al., 2016; Landsberg and Waring, 1997). 3-PG is a stand-level vegetation model that runs at a monthly time step and includes vegetation carbon dynamics and a simple soil
<table><caption>**Table 1.** Regional observational data streams used in data assimilation.</caption><thead><tr><th>Data stream</th><th>Measurement frequency</th><th>Measurement or estimation technique</th><th>Uncertainty</th><th>Stream ID for Table 3</th></tr></thead><tbody><tr><td>Foliage biomass (Pine)</td><td>Annual or less</td><td>Allometric relationship</td><td>Based on propagating the allometric model uncertainty in Gonzalez-Benecke et al. (2014). Varied by observation.</td><td>1</td></tr><tr><td>Foliage biomass (hardwood)</td><td>Annual or less</td><td>Allometric relationship</td><td>Assumed zero</td><td>2</td></tr><tr><td>Stem biomass (pine)</td><td>Annual or less</td><td>Allometric relationship</td><td>Based on propagating the allometric model uncertainty in Gonzalez-Benecke et al. (2014). Varied by observation.</td><td>3</td></tr><tr><td>Stem biomass (hardwood)</td><td>Annual or less</td><td>Allometric relationship</td><td>Assumed zero</td><td>4</td></tr><tr><td>Coarse root biomass (combined)</td><td>Annual or less</td><td>Allometric relationship</td><td>Assumed zero*</td><td>5</td></tr><tr><td>Fine root biomass (combined)</td><td>Annual or less</td><td>Allometric relationship</td><td>SD: 10% of observation</td><td>6</td></tr><tr><td>Foliage biomass production (combined)</td><td>Annual</td><td>Litterfall traps</td><td>SD: 10% of observation</td><td>7</td></tr><tr><td>Fine root biomass production (combined)</td><td>Annual</td><td>Mini-rhizotrons</td><td>SD: 10% of observation</td><td>8</td></tr><tr><td>Pine stem density</td><td>Annual or less</td><td>Counting individuals</td><td>1% (assumed small)</td><td>9</td></tr><tr><td>Leaf area index (pine)</td><td>Monthly to annual</td><td>Litter traps or LI 2000</td><td>SD: 10% of observation</td><td>10</td></tr><tr><td>Leaf area index (hardwood)</td><td>Monthly to annual</td><td>Litter traps or LI 2000</td><td>SD: 10% of observation</td><td>11</td></tr><tr><td>Leaf area index (combined)</td><td>Only used if not separated into pine and hardwood</td><td>Litter traps or LI 2000</td><td>SD: 10% of observation</td><td>12</td></tr><tr><td>Gross ecosystem production</td><td>Monthly</td><td>Modeled from flux eddy-covariance net ecosystem exchange</td><td>SD: 10% of observation</td><td>13</td></tr><tr><td>Evapotranspiration</td><td>Monthly</td><td>Eddy covariance</td><td>SD: 10% of observation</td><td>14</td></tr></tbody></table>

water bucket model (Fig. 2). While a complete description of the 3-PG model and our modifications can be found in the Supplement Sect. 1, the key concept for interpreting the results is that gross primary productivity (GPP) was simulated using a light-use efficiency approach where the absorbed photosynthetically active radiation (APAR) was converted to carbon based on a quantum yield (Supplement Sect. 1.1). Quantum yield was simulated using a parameterized maxi-

mum quantum yield (alpha) that was modified by environmental conditions including atmospheric CO₂, available soil water (ASW), and soil fertility (Supplement Sect. 1.2–1.3). The ASW and soil fertility modifiers were values between 0 and 1, while the atmospheric CO₂ modifier had a value of 1 at 350 ppm (thus values greater than 1 at higher CO₂ concentrations).

<footnote>* The relatively low number of observations prevented convergence when using the observational uncertainty model, so observational uncertainty was assumed to be zero to allow convergence.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 1.** Map of loblolly pine distribution, plot locations used in data assimilation, and the experiment type associated with each plot. The control-only treatments were plots without any associated experimental treatment or flux measurements. Fertilized treatments were plots with nutrient additions. CO₂ treatments were plots with free-air concentration enrichment treatments. The flux treatments were plots with eddy-covariance measurements of ecosystem-scale carbon and water exchange. The water treatments included throughfall exclusion and irrigation experiments.

<figure><img src="image_5.png" /></figure>

**Figure 2.** A diagram of the monthly time-step 3-PG model used in this study. The stocks are represented by the boxes and the fluxes by the arrows. An influence of a stock on a flux that is not directly related to that stock is represented by the dotted lines. The environmental influences on a flux are described using italics. A description of the model can be found in the Supplement.

Elevated CO₂ modified tree physiology by increasing quantum yield, based on an increasing but saturating relationship with atmospheric CO₂ (Supplement Sect. 1.2). Based on initial results from the data assimilation, we also added a function where the allocation to foliage relative to stem biomass decreased as atmospheric CO₂ increased (Sup-

plement Sect. 1.2). ASW and quantum yield were positively related through a logistic relationship between relative ASW and the quantum yield modifier, where relative ASW was the ratio of simulated ASW to a plot-level maximum ASW. Soil fertility and quantum yield were proportionally related, where quantum yield was scaled by an estimate of relative
stand-level fertility (a value of 1 was the maximum fertility). The fertility modifier (or soil fertility rating, FR) was constant throughout a simulation of a plot and was either based on site characteristics or directly optimized as a stand-level parameter (Supplement Sect. 1.3). For plots with nutrient fertilization, FR was a directly optimized parameter or set to 1, depending on the level of fertilization (see below). For unfertilized plots, we used site index (SI), a measure of the height of a stand at a specified age (25 years), to estimate FR. This approach is in keeping with previous efforts (Gonzalez-Benecke et al., 2016; Subedi et al., 2015); however, SI does not solely represent the nutrient availability of an ecosystem. For a given climate SI captures differences in soil fertility, where a lower SI corresponded to a site with lower fertility, but regional variation in SI also included the influence of climate on growth rates that were already accounted for in the other environmental modifiers in the 3-PG model. When a climate term is not used in the empirical FR model, FR is relative to the highest SI in the region, which does not occur in the northern extent of the region even in fertilized plots due to climatic constraints. Thus, we also included the historical (1970–2011) 35-year mean annual temperature (MAT) as an additional predictor, resulting in an empirical relationship that predicted FR as an increasing, but saturating, function of SI within areas of similar long-term temperature. For our application of the 3-PG model using DA, we removed the previously simulated dependence of total root allocation on FR (Bryars et al., 2013; Gonzalez-Benecke et al., 2016) because we separated coarse and fine roots. Other environmental conditions influenced GPP, including temperature, frost days, and vapor pressure deficit (VPD). A description of these modifiers can be found in Supplement Sect. 1.2.

Each month, net primary production (a parameterized and constant proportion of GPP) was allocated to foliage, stem (stem wood, stem bark, and branches), coarse roots, and fine roots (Supplement Sect. 1.4). Differing from previous applications of 3-PG to loblolly pine ecosystems, we modified the model to simulate fine roots and coarse roots separately. 3-PG also simulated simple population dynamics by including stem density as a state variable. Stem density and stem biomass pools were reduced by both density-dependent mortality, based on the concept of self-thinning (Landsberg and Waring, 1997), and density-independent mortality, a new modification where a constant proportion of individuals die each month (Supplement Sect. 1.5). Finally, we added a simple model of hardwood understory vegetation to enable the assimilation of GEP and ET observations from eddy-covariance tower studies with significant understories (Supplement Sect. 1.7).

The water cycle was a simple bucket model with transpiration predicted using a Penman-Monteith approach (Bryars et al., 2013; Gonzalez-Benecke et al., 2016; Landsberg and Waring, 1997) (Supplement Sect. 1.6). The canopy conductance used in the Penman-Monteith subroutine was modified by environmental conditions. The modifiers included the

same ASW and VPD modifier as used in the GPP calculation. Maximum canopy conductance occurred when simulated LAI exceeded a parameterized value of LAI (LAIgcx). Evaporation was equal to the precipitation intercepted by the canopy. Runoff occurred when the ASW exceeded a plot-specific maximum ASW. As in prior applications of 3-PG, ASW was not allowed to take a value below a minimum ASW, resulting in an implicit irrigation in very dry conditions. This assumption may cause the model to be less sensitive to low ASW, but the optimized parameterization may compensate for this.

The 3-PG model used in this study simulated the monthly change in 11 state variables per plot: four stocks for loblolly pines, five stocks for understory hardwoods, loblolly pine stem density (stems ha⁻¹), and ASW. The key fluxes that were used for DA included monthly GEP, monthly ET, annual root turnover, and annual foliage turnover. In total, 46 parameters were required by 3-PG. The model required mean daily maximum temperature, mean daily minimum temperature, mean daily PAR, total frost days per month, total rain per month, annual atmospheric CO₂, and latitude. Each plot also required maximum ASW, SI, MAT, and the initial condition of the 11 state variables as model inputs (Fig. 3).

We used the first observation at the plot as the initial conditions for the loblolly pine vegetation states (foliage biomass, stem biomass, coarse root biomass, fine root biomass, and stem number). When observations of coarse biomass and fine root biomass were not available, these stocks were initialized as a mean region-wide proportion of the observed stem biomass. However, the value of initial root biomass in plots without observations was not important because root biomass did not influence any other functions in the model. The hardwood understory stocks at US-DK3 and US-NC2 were also initialized using the first set of observations. Initial fine root and coarse biomass were distributed between loblolly pine and hardwoods based on their relative contribution of total initial foliage biomass. The initialized ASW was assumed to be equal to the maximum ASW because most plots were initialized in winter months when plant demand for water was minimal. The maximum ASW in each plot was extracted from the Soil Survey Geographic Database (SSURGO) soils dataset (Soil Survey Staff, 2013). The value we used corresponded to the maximum ASW for the top 1.5 m of the soil. We assumed that the minimum ASW was zero. Because we focused on a region-wide optimization, we used region-wide 4 km estimates of observed monthly meteorology as inputs and to calculate the 35-year MAT for each plot (Abatzoglou, 2013). SI was based on height measurements at age 25 in each plot or calculated by combining observations of height at younger ages with an empirical model (Dieguez-Aranda et al., 2006).

We simulated ecosystem manipulation experiments in the 3-PG model by altering the environmental modifiers or by modifying the environmental inputs. Nutrient addition experiments were simulated by setting FR equal to 1 for the studies
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Key climatic and stand characteristic inputs to the regional 3-PG simulations: (a) mean annual temperature (1979–2011) as a summary of the gradient in monthly temperature inputs used in simulations, (b) maximum available soil water for the top 1.5 m of soil from SSURGO, (c) mean annual precipitation (1979–2011) as a summary of the gradient in monthly precipitation inputs used in simulations, and (d) site index. The area shown is the natural range of loblolly pine (*Pinus taeda* L.).</figcaption></figure>

that applied nutrients at regular intervals to remove nutrient deficiencies (RW18, SETRES, Waycross). FR was directly estimated for fertilized plots in two of the studies either because nutrients were only added once at the beginning of the study (PINEMAP), thus potentially not removing nutrient limitation, or because nitrogen was the only element added (Duke FACE), thus allowing the potential for nutrient limitation by other elements. For these plots, we also assumed that the FR of the fertilized plot was equal to or larger than the control plot. Throughfall exclusion experiments were simulated by decreasing the throughfall by 30% in the treatment plots. The SETRES irrigation experiments were simulated by adding 650 mm to ASW between April and October. CO₂ enrichment experiments were simulated by setting the atmospheric CO₂ input equal to the treatment mean from the elevated CO₂ rings (570 ppm). One plot (US-NC2) included a thinning treatment during the period of observation. We simulated the thinning by specifying a decrease in the stem count that matched the proportion removed at the site, with the biomass of each tree equivalent to the average of trees in the plot.

## 2.3 Data assimilation method

We used a hierarchical Bayesian framework to estimate the posterior distributions of parameters, latent states of stocks and fluxes, and process uncertainty parameters. The latent states represented a value of the stock or flux before uncertainty was added through measurement. The approach was as follows.

Consider a stock or flux (m) for a single plot (p) at time t ($q_{p,m,t}$). $q_{p,m,t}$ is influenced by the processes represented in the 3-PG model and a normally distributed model process error term,

$$q_{p,m,t} \sim N(f(\theta, \text{FR}_p), \sigma_m), \quad (1)$$

where $\theta$ is a vector of parameters that are optimized, $\text{FR}_p$ is the site fertility, and $\sigma_m$ is the model process error. Not shown are the vector of parameters that were not optimized (Supplementary Table S1), the plot ASW, an array of climate inputs, and the initial conditions because these were assumed known and not estimated in the hierarchical model. The process error assumed that the error linearly scales with the magnitude
of the prediction:

$$ \sigma_m^2 = \gamma_m + \rho_m f(\theta, \text{FR}_p). \qquad (2) $$

While the structure of the Bayesian model allowed for all data streams to have process uncertainty that scales with the prediction, in this application we only allowed stem biomass, GEP, and ET process uncertainty to scale because they had large variation across space (stem biomass) and through time (i.e., there should be lower process uncertainty in the winter when GEP is lower). For the other data streams, the linear scaling term was removed by fixing $\rho_m$ at 0.

FRp did not have an explicit probability distribution. Rather the probability density was evaluated as 1 if the plot was not fertilized, thus causing FRp to be estimated from SI and MAT (Supplement Eq. 15), or if it was a fertilized plot and had an FRp equal or higher than that of its non-fertilized control plot. The probability density was evaluated as 0 if the estimated FRp in a fertilized plot was less than the FRp in the control plot or if FRp was not contained in the interval between 0 and 1.

$$ \text{FR}_p \sim \begin{cases} 1 & \text{if non-fertilized, } \text{FR}_p \ge 0, \text{ and } \text{FR}_p \le 1 \\ 1 & \text{if } \text{FR}_p = 1 \text{ and fertilization levels are assumed to remove nutrient deficiencies} \\ 0 & \text{if } \text{FR}_p < 1 \text{ and fertilization levels are assumed to remove nutrient deficiencies} \\ 1 & \text{if fertilized but levels are not assumed to remove deficiencies and } \text{FR}_p \ge \text{FR of control plot} \\ 0 & \text{if fertilized but levels are not assumed to remove deficiencies and } \text{FR}_p < \text{FR of control plot} \\ 0 & \text{if } \text{FR}_p < 0 \text{ or } \text{FR}_p > 1 \end{cases} \qquad (3) $$

Our model included the effect of observational errors for measurements of stocks and fluxes. For a single stock or flux for a plot at time t there was an observation ($y_{p,m,t}$). The normally distributed observation error model was

$$ y_{p,m,t} \sim N(q_{p,m,t}, \tau_{p,m,t}^2), \qquad (4) $$

where $\tau_{p,m,t}^2$ represented the measurement error of the observed state or flux. By including the observational error model, $q_{p,m,t}$ represented the latent, or unobserved, stock or flux. The variance was unique to each observation because it was represented as a proportion of the observed value. The $\tau_{p,m,t}^2$ was assumed known (Table 1) and not estimated in the hierarchical model.

The hierarchical model required prior distributions for all optimized parameters, including the parameters for the 3-PG model ($\theta$), FRp, and the process error parameters. The prior distributions for ($p(\theta)$) are specified in Table 3. Some parameters were informed by previous research in loblolly pine ecosystems, while other parameters were “uninformative” with flat distributions that had broad, but physically reasonable, bounds. The prior distributions for the process error parameters were non-informative and had a uniform distribution with upper and lower bounds that spanned the range of

reasonable error terms.

$$ \gamma_m \sim U(0.001, 100) \qquad (5) $$

$$ \rho_m \sim U(0, 10) \qquad (6) $$

By combining the data, process, and prior models, our joint posterior that includes all 13 data streams, plots, months with observations, and fitted parameters was

$$ \begin{align} p(\theta, y, \gamma, q | y, \tau, \text{priors}) &\propto \nonumber \\ & \prod_{p=1}^{P} \prod_{m=1}^{M} \prod_{t=1}^{T} N(q_{p,m,t} | f(\theta, \text{FR}_p), \gamma_m + \rho_m f(\theta, \text{FR}_p)), \nonumber \\ & \prod_{p=1}^{P} \prod_{m=1}^{M} \prod_{t=1}^{T} N(y_{p,m,t} | q_{p,m,t}, \tau_{p,m,t}^2), \tag{7} \\ & \prod_{p=1}^{P} p(\text{FR}_p) \prod_{f=1}^{F} p(\theta_f) \prod_{m=1}^{M} p(\gamma_m) \prod_{m=1}^{M} p(\rho_m), \nonumber \end{align} $$

where bolded components represent vectors, P is the total number of plots, M is the total number of data streams, T is the total months with observations, and F is the total number of 3-PG parameters that are optimized.

We numerically estimated the joint posterior distribution using the Monte Carlo Markov Chain–Metropolis Hasting (MCMC-MH) algorithm (Zobitz et al., 2011). This approach has been widely used to approximate parameter distributions in ecosystem DA research (Fox et al., 2009; Trudinger et al., 2007; Williams et al., 2005; Zobitz et al., 2011). Briefly, the algorithm proposed new values for the model parameters, uncertainty parameters, latent states, and FR. The proposed values were generated using a random draw from a normal distribution with a mean equal to the previously accepted value for that parameter and standard deviation equal to the parameter-specific jumping size. The ratio of the proposed calculation of Eq. (7) to the previously accepted calculation of Eq. (7) was used to determine if the proposed parameter was accepted. If the ratio was greater than or equal to 1, the proposed value was always accepted. If the ratio was less than 1, a random number between 0 and 1 was drawn and the proposed value was accepted if the ratio was greater than the random number. This allowed less probable parameter sets to be accepted, thus sampling the posterior distribution. We adapted the size of the jump size for each parameter to ensure the acceptance rate of the parameter set was between 22 and 43% (Ziehn et al., 2012) by adjusting the jump size if the acceptance rate for a parameter was outside the 22–43% range. All MCMC-MH chains were run for 30 million iterations with the first 15 million iterations discarded as the burn-in. Four chains were run and tested for convergence using the Gelman–Rubin convergence criterion, where a value for the criterion less than 1.1 indicated an acceptable level of convergence. We sampled every 1000th parameter in the final 15 million iterations of the MCMC-MH chain and used this thinned chain in the analysis described below. The 3-PG
model and MCMC-MH algorithm were programmed in Fortran 90 and used OpenMP to parallelize the simulation of each plot within an iteration of the MCMC-MH algorithm.

## 2.4 Data assimilation evaluation

Using the observations, model, and hierarchical Bayesian method described above, we assimilated both the non-manipulated and manipulated plots (Base assimilation; Table 4). We assessed model performance first by calculating the RMSE and bias of stem biomass predictions (the most common data stream). In the evaluation, we only used the most recent observed values to increase the time length between initialization and validation. Second, we assessed the predictive capacity by comparing model predictions to data not used in the parameter optimization in a cross-validation study. In this evaluation, we repeated the Base assimilation without 160 FMRC thinning study plots (Table 2), predicted the 160 plots using the median parameter values, and calculated the RMSE and bias stem biomass of the independent set of plots. Rather than holding out all 160 plots from a single assimilation and not generating a converged chain, we divided the 160 plots into four unique sets of 40 plots and repeated the assimilation for each set. Finally, we compared the predicted responses to experimental manipulation to the observed responses. We focused the comparison on the percentage difference in stem biomass between the control and treatment plots. We used a paired *t* test to test for differences between the predicted and observed responses within an experimental type (irrigated, drought, nutrient addition, and elevated CO₂). We combined the single and multi-factor treatments for analysis. For the analysis of the nutrient addition studies, we only used plots where FR was assumed to be 1 so that we were able to simulate the treatments without requiring the optimization of a site-specific FR parameter.

During preliminary analysis, we found that the Base assimilation predicted lower stem biomass than observed in the elevated CO₂ plots in the Duke FACE study. Further analysis investigating the cause of the bias in the CO₂ plots showed that three parameters (wSx1000, ThinPower, and pCRS) were required to be unique to the Duke FACE study in order to reduce the bias. Therefore, the Base assimilation included unique parameters for wSx1000, ThinPower, and pCRS parameters in all plots in the Duke FACE and US-DK3 studies. To highlight the need for the site-specific parameters, we repeated the Base assimilation approach without the three additional parameters for the Duke studies (NoDkPars assimilation).

## 2.5 Sensitivity to the inclusion of ecosystem experiments

We also evaluated how parameter distributions and the associated environmental sensitivity of model predictions depended on the inclusion of ecosystem experiments in data as-

similation. First, we repeated the Base assimilation, this time excluding the plots that included the manipulated treatments (NoExp). We removed all manipulation types at once, rather than individual experimental types, because all experimental types involved multi-factor studies. The NoExp assimilation had the same number of data streams as the Base assimilation because it included the control treatments from the experimental studies. The NoExp assimilation represented the situation where only observations across environmental gradients were available. Second, we compared the parameterization of the ASW, soil fertility, and atmospheric CO₂ environmental modifiers from the Base to the NoExp assimilation. The modifier equations are described in Supplement Sects. 1.2 and 1.3. Third, we repeated the same independent validation exercise for the 160 FMRC plots as described above for the Base assimilation. Fourth, we predicted the treatment plots in the irrigated, drought, nutrient addition (only plots where FR was assumed to be 1), and elevated CO₂ plots. As for the Base assimilation, we used a *t* test to compare the experimental response between the NoExp assimilation and observed values and between the NoExp and Base assimilations. Since the experimental treatments were not used in the optimization, this was an independent evaluation of predictive capacity.

## 2.6 Regional predictions with uncertainty

To demonstrate the capacity of the data assimilation system to create regional predictions with uncertainty, we simulated the regional response to a decrease in precipitation, an increase in nutrient availability, and an increase in atmospheric CO₂ concentration, each as a single factor change from a 1985–2011 baseline. Each prediction included uncertainty by integrating across the parameter posterior distributions using a Monte Carlo sample of the parameter chains. Our region corresponded to the native range of loblolly pine and used the HUC12 (USGS 12-digit Hydrological Unit Code) watershed as the scale of simulation. For each HUC12 in the region, we used the mean SI, 30-year mean annual temperature, ASW aggregated to the HUC12 level, and monthly meteorology from Abatzoglou (2013) as inputs (Fig. 3). The SI of each HUC12 was estimated from biophysical variables in the HUC12 using the method described in Sabatia and Burkhart (2014). This SI corresponded to an estimated SI for stands without intensive silvicultural treatments or advanced genetics of planted stock.

To sample parameter uncertainty, we randomly drew 500 samples, with replacement, from the Base assimilation MCMC chain and simulated forest development from a 1985 planting to age 25 in 2011 in each HUC. We chose age 25 as the final age because it is a typical age of harvest in the region. For each sample, we repeated the regional simulation with (1) a 30% reduction in precipitation, (2) FR set to 1, and (3) atmospheric CO₂ increased by 200 ppm. Within a parameter sample, we calculated the percent change in stem
<table><caption>**Table 2.** Descriptions of the studies used in data assimilation.</caption><thead><tr><th>Study name</th><th>Number of locations</th><th>Number of plots per site</th><th>Experimental treatments (plots)</th><th>Data streams (Table 2)</th><th>Measurement years</th><th>Measurement stand ages (years)</th><th>Reference</th></tr></thead><tbody><tr><td>FMRC<sup>a</sup><br> thinning study</td><td>163</td><td>1</td><td>None</td><td>1, 3, 9</td><td>1981–2003</td><td>8–30</td><td>Burkhart et al. (1985)</td></tr><tr><td>FPC<sup>b</sup><br> Region-wide 18</td><td>18</td><td>2</td><td>Nutrient addition</td><td>1, 3, 9</td><td>2011–2014</td><td>12–21</td><td>Albaugh et al. (2015)</td></tr><tr><td>PINEMAP<sup>c</sup></td><td>4</td><td>16</td><td>Nutrient addition, 30% throughfall, nutrient × throughfall</td><td>1, 3, 9</td><td>2011–2015</td><td>3–13</td><td>Will et al. (2015)</td></tr><tr><td>Waycross</td><td>1</td><td>2</td><td>Nutrient addition</td><td>3, 9, 10</td><td>1991–2010</td><td>4–23</td><td>Bryars et al. (2013)</td></tr><tr><td>SETRES<sup>d</sup></td><td>1</td><td>16</td><td>Nutrient addition, irrigation, nutrient × irrigation</td><td>1, 3, 5, 6, 9, 10</td><td>1991–2006</td><td>8–23</td><td>Albaugh et al. (2004)</td></tr><tr><td>Duke FACE<sup>e</sup><br> and US-DK3 flux</td><td>1</td><td>12</td><td>CO<sub>2</sub>, nutrient addition, CO<sub>2</sub> × nutrient addition</td><td>2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14</td><td>1996–2004</td><td>13–22</td><td>McCarthy et al. (2010); Novick et al. (2015)</td></tr><tr><td>NC2 flux</td><td>1</td><td>1</td><td>None</td><td>2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14</td><td>2005–2014</td><td>12–22</td><td>Noormets et al. (2010)</td></tr><tr><td>Total</td><td>187</td><td>294</td><td></td><td></td><td>1981–2014</td><td>4–30</td><td></td></tr></tbody></table>

biomass at age 25 between the control simulation and the three simulations with the environmental changes. We focused our regional analysis on the distribution of the percent change in stem biomass.

# 3 Results

## 3.1 Data assimilation evaluation

Our multisite, multi-experiment, multi-data stream DA approach (Base assimilation) increased confidence in the model parameters (Table 5). Averaged across parameters, the posterior 99% quantile range from the Base assimilation was 60% less than the prior range. The largest reduction in parameter uncertainty was for the parameters associated with light-use efficiency (alpha) and the conversion of GPP to net primary productivity (NPP) (y), which on average had ranges that were 85% lower in the posterior than the prior. Parameters associated with allocation and allometry had a 63% reduc-

tion in the range while parameters associated with mortality processes had a 70% reduction in the range. Parameters associated with environmental modifiers had the least reduction in the range with a 40% decrease. In addition to the parameters associated with the 3-PG model, the model process error parameters for each data stream were well constrained with large reductions in the range (>99% decrease; Supplemental Table S2)

The Base assimilation reliably predicted data from the regionally distributed non-manipulated plots that were not used in the optimization. The mean bias in stem biomass of the cross-validation was -3.7% and the RMSE was 21.8 Mg ha⁻¹ (Fig. 4a). Furthermore, the response of stem biomass to irrigation (d f = 7, p = 0.18), nutrient addition (d f = 26, p = 0.29), and elevated CO₂ (d f = 4, p = 0.43) was not significantly different between the observed and the Base assimilation (Fig. 5). The Base assimilation was significantly more sensitive to drought than observed (n = 31, p < 0.001; Fig. 5).

<footnote><sup>a</sup> Forest Modeling Research Cooperative. <sup>b</sup> Forest Productivity Cooperative. <sup>c</sup> PINEMAP. <sup>d</sup> Southeast Tree Research and Education Site. <sup>e</sup> Free-Air Carbon Enrichment.</footnote>
<table><caption>**Table 3.** The prior distributions of all 3-PG model parameters optimized using data assimilation. NPP: net primary production.</caption><thead><tr><th>Parameter</th><th>Parameter description</th><th>Units</th><th>Prior distribution</th><th>Prior parameters</th><th>Reference for prior (see footnote)</th></tr></thead><tbody><tr><td colspan="6">Allocation and structure</td></tr><tr><td>pFS2</td><td>Ratio of foliage to stem allocation at stem diameter: 2 cm</td><td>–</td><td>Uniform</td><td>Min: 0.08<br>Max: 1.00</td><td>Uninformed</td></tr><tr><td>pFS20</td><td>Ratio of foliage to stem allocation at stem diameter: 20 cm</td><td>–</td><td>Uniform</td><td>Min: 0.10<br>Max: 1.00</td><td>Uninformed</td></tr><tr><td>pRF</td><td>Ratio of fine roots to foliage allocation</td><td>–</td><td>Uniform</td><td>Min: 0.05<br>Max: 2.00</td><td>Uninformed</td></tr><tr><td>pCRS</td><td>Ratio of coarse roots to stem allocation</td><td>–</td><td>Uniform</td><td>Min: 0.15<br>Max: 0.35</td><td>1</td></tr><tr><td>SLA0</td><td>Specific leaf area at stand age 0</td><td>m<sup>2</sup> kg<sup>-1</sup></td><td>–</td><td>mean: 5.53<br>SD: 0.44</td><td>2</td></tr><tr><td>SLA1</td><td>Specific leaf area for mature aged stands</td><td>m<sup>2</sup> kg<sup>-1</sup></td><td>Normal</td><td>mean: 3.58<br>SD: 0.11</td><td>2</td></tr><tr><td>tSLA</td><td>Age at which specific leaf area is 0.5 (SLA0 + SLA1)</td><td>Years</td><td>Normal</td><td>mean: 5.97<br>SD: 2.15</td><td>2</td></tr><tr><td>fCpFS700</td><td>Proportional decrease in allocation to foliage between 350 and 700 ppm CO<sub>2</sub></td><td>–</td><td>Uniform</td><td>Min: 0.50<br>Max: 1.00</td><td>Uninformed</td></tr><tr><td>StemConst</td><td>Constant in stem mass vs. diameter relationship</td><td>–</td><td>Normal</td><td>mean: 0.022<br>SD: 0.005</td><td>3</td></tr><tr><td>StemPower</td><td>Power in stem mass vs. diameter relationship</td><td>–</td><td>Normal</td><td>mean: 2.77<br>SD: 0.2</td><td>3</td></tr><tr><td colspan="6">Canopy photosynthesis, autotrophic respiration, and transpiration</td></tr><tr><td>alpha</td><td>Canopy quantum efficiency (pines)</td><td>mol C mol PAR<sup>-1</sup></td><td>Uniform</td><td>Min: 0.02<br>Max: 0.06</td><td>Uninformed</td></tr><tr><td>y</td><td>Ratio NPP / GPP</td><td>–</td><td>Uniform</td><td>Min: 0.30<br>Max: 0.65</td><td>4</td></tr><tr><td>MaxCond</td><td>Maximum canopy conductance</td><td>m s<sup>-1</sup></td><td>Uniform</td><td>Min: 0.005<br>Max: 0.03</td><td>2</td></tr><tr><td>LAIgcx</td><td>Canopy LAI for maximum canopy conductance</td><td>–</td><td>Uniform</td><td>Min: 2<br>Max: 5</td><td>2, 5, 6</td></tr><tr><td colspan="6">Environmental modifiers of photosynthesis and transpiration</td></tr><tr><td>kF</td><td>Reduction rate of production per °C below zero</td><td>–</td><td>Normal</td><td>mean: 0.18<br>SD: 0.016</td><td>2</td></tr><tr><td>T<sub>min</sub></td><td>Minimum monthly mean temperature for photosynthesis</td><td>°C</td><td>Normal</td><td>mean: 4.0<br>SD: 2.0</td><td>2, 5, 6</td></tr><tr><td>T<sub>opt</sub></td><td>Optimum monthly mean temperature for photosynthesis</td><td>°C</td><td>Normal</td><td>mean: 25.0<br>SD: 2.0</td><td>2, 5, 6</td></tr><tr><td>T<sub>max</sub></td><td>Maximum monthly mean temperature for photosynthesis</td><td>°C</td><td>Normal</td><td>mean: 38.0<br>SD: 2.0</td><td>2, 5, 6</td></tr></tbody></table>

The plots at the Duke Forest study had a higher carrying capacity of stem biomass before self-thinning (WSx1000), lower self-thinning rate (ThinPower), and smaller allocation to coarse root (pCRS) than values optimized from the other plots across the region (Table 6). The DA approach without these three study-specific parameters (NoDkPars) predicted significantly lower accumulation of stem biomass in response to elevated CO₂ than observed (d f = 4, p = 0.002;

Fig. 5). The NoDKPars assimilation optimized the CO₂ fertilization parameter (fAlpha700) to a value that predicted 45% less light-use efficiency at 700 ppm (1.13 in NoDKPar vs. 1.33 in Base; Table 6) than the Base assimilation.
<table><caption>**Table 3. Continued.**</caption><thead><tr><th>Parameter</th><th>Parameter description</th><th>Units</th><th>Prior distribution</th><th>Prior parameters</th><th>Reference for prior (see footnote)</th></tr></thead><tbody><tr><td>SWconst</td><td>Moisture ratio deficit when downregulation is 0.5</td><td>–</td><td>Uniform</td><td>Min: 0.01<br>Max: 1.8</td><td>Uninformed</td></tr><tr><td>SWpower</td><td>Power of moisture ratio deficit</td><td>–</td><td>Uniform</td><td>Min: 1<br>Max: 13</td><td>Uninformed</td></tr><tr><td>CoeffCond</td><td>Defines stomatal response to VPD</td><td>mbar<sup>-1</sup></td><td>Normal</td><td>mean: 0.041<br>SD: 0.003</td><td>2</td></tr><tr><td>fAlpha700</td><td>Proportional increase in canopy quantum efficiency between 350 and 700 ppm CO<sub>2</sub></td><td>–</td><td>Uniform</td><td>Min: 1.00<br>Max: 1.8</td><td>Uninformed</td></tr><tr><td>MaxAge</td><td>Maximum stand age used to compute relative age</td><td>Years</td><td>Uniform</td><td>Min: 16<br>Max: 200</td><td>Uninformed</td></tr><tr><td>nAge</td><td>Power of relative age in the age modifier</td><td>–</td><td>Uniform</td><td>Min: 0.2<br>Max: 4.0</td><td>Uninformed</td></tr><tr><td>rAge</td><td>Relative age to where age modifier was 0.5</td><td>–</td><td>Uniform</td><td>Min: 0.01<br>Max: 3.00</td><td>Uninformed</td></tr><tr><td>FR1</td><td>Fertility rating parameter 1 (mean annual temperature coefficient)</td><td>–</td><td>Uniform</td><td>Min: 0.0<br>Max: 1.0</td><td>Uninformed</td></tr><tr><td>FR2</td><td>Fertility rating parameter 2 (site index age 25 coefficient)</td><td>–</td><td>Uniform</td><td>Min: 0.0<br>Max: 1.0</td><td>Uninformed</td></tr><tr><td colspan="6"><strong>Mortality</strong></td></tr><tr><td>wSx1000</td><td>Maximum stem mass per tree at 1000 trees ha<sup>-1</sup></td><td>kg tree<sup>-1</sup></td><td>Normal</td><td>mean: 235<br>SD: 25</td><td>2, 5, 6</td></tr><tr><td>ThinPower</td><td>Power in self-thinning law</td><td>–</td><td>Uniform</td><td>Min: 1.0<br>Max: 2.5</td><td>2, 5, 6</td></tr><tr><td>mS</td><td>Fraction of mean stem biomass per tree on dying trees</td><td>–</td><td>Uniform</td><td>Min: 0.1<br>Max: 1.0</td><td>Uninformed</td></tr><tr><td>Rttover</td><td>Average monthly root turnover rate</td><td>month<sup>-1</sup></td><td>Uniform</td><td>Min: 0.017<br>Max: 0.042</td><td>7</td></tr><tr><td>MortRate</td><td>Density-independent mortality rate (pines)</td><td>month<sup>-1</sup></td><td>Uniform</td><td>Min: 0.0002<br>Max: 0.004</td><td>Uninformed</td></tr><tr><td colspan="6"><strong>Understory hardwoods</strong></td></tr><tr><td>alpha_h</td><td>Canopy quantum efficiency (understory hardwoods)</td><td>mol C mol PAR<sup>-1</sup></td><td>Uniform</td><td>Min: 0.005<br>Max: 0.07</td><td>Uninformed</td></tr><tr><td>pFS_h</td><td>Ratio of foliage to stem partitioning (understory hardwoods)</td><td>–</td><td>Uniform</td><td>Min: 0.2<br>Max: 3.0</td><td>Uninformed</td></tr><tr><td>pR_h</td><td>Ratio of foliage to fine roots (understory hardwoods)</td><td>–</td><td>Uniform</td><td>Min: 0.05<br>Max: 2</td><td>Uninformed</td></tr><tr><td>SLA_h</td><td>Specific leaf area (understory hardwoods)</td><td>m<sup>2</sup> kg<sup>-1</sup></td><td>Normal</td><td>mean: 16<br>SD: 3.8</td><td>8</td></tr><tr><td>fAlpha700_h</td><td>Proportional increase in canopy quantum efficiency between 350 and 700 ppm CO<sub>2</sub> (understory hardwood)</td><td>–</td><td>Uniform</td><td>Min: 1.00<br>Max: 2.5</td><td>Uninformed</td></tr></tbody></table>

## 3.2 Sensitivity to the inclusion of ecosystem experiments

Excluding the experimental treatments from the data assimilation did not strongly influence the predictive capacity of

the model. The RMSE validation plots in NoExp assimilation decreased slightly compared to Base assimilation (21.8 to 18.0 Mg ha<sup>-1</sup>), while the bias slightly increased (–3.7 to –4.1 %) (Fig. 4b). Excluding the experimental treatments resulted in a significantly lower response of stem biomass

<footnote>1: Albaugh et al., 2005. 2: Gonzalez-Benecke et al., 2016. 3: Gonzalez-Benecke et al., 2014. 4: DeLucia et al., 2007. 5: Bryars et al., 2013. 6: Subedi et al., 2015.</footnote>

<footnote>7: Matamala et al., 2003. 8: LeBauer et al., 2010. Uninformed priors had large, ecologically reasonable bounds.</footnote>
<figure><img src="image_7.png" /></figure>

**Figure 4.** Model evaluation of stem biomass when assimilating (a) observations across environmental gradients and ecosystem manipulation experiments (Base; Table 4) and (b) only observations across environmental gradients (NoExp; Table 4). The gray circles correspond to predictions where all plots were used in data assimilation. The black triangles correspond to predictions where 160 plots were not included in data assimilation and represent an independent evaluation of model predictions (cross-validation). For each plot, we used the measurement with the longest interval between initialization and measurement for evaluation.

<figure><img src="image_8.png" /></figure>

**Figure 5.** The mean response, expressed as a percentage change in stem biomass from the control treatment, for irrigation, drought (as a reduction in throughfall), nutrient addition, and elevated CO₂ experiments. The observed response and the response simulated by the Base, NoExp, and NoDkPars assimilation approaches are shown. The # sign signifies that the value below the marker was significantly different from the observed response (p < 0.05). The * sign signifies that the value below the marker was significantly different from the response in the Base assimilation (p < 0.05). Error bars are ±1 standard deviation.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Optimized environmental response functions in the 3-PG model for the (a) soil fertility influence on photosynthesis, (b) available soil water influence on photosynthesis and conductance, and (c) atmospheric CO₂ influence on photosynthesis. The function shapes were derived from the parameters in the Base, NoExp, and NoDkPars assimilations (Table 4).</figcaption></figure>

to elevated CO₂ than observed (df = 4, p < 0.001; Fig. 5). Furthermore, there was a slight negative response of stem biomass to CO₂ in the NoExp assimilation because the parameter governing the change in foliage allocation at elevated CO₂ (fCpFS700) was unconstrained by observations (Table 6). This led to convergence on the lower bound of the prior distribution (0.5) where foliage allocation decreased with increased atmospheric CO₂. The predictions of irrigation, drought, and nutrient addition experiments were not significantly different between the Base and NoExp assimilations (Fig. 5).

The parameters and associated response functions in the 3-PG for nutrients, ASW, and atmospheric CO₂ differed between the Base and NoExp assimilations (Fig. 6). First, the parameterization of the soil fertility (FR) showed a stronger dependence on SI in the NoExp assimilation than in the Base assimilation (Fig. 6a). For a given SI there was a lower FR and thus stronger nutrient limitation, when experimental treatments were excluded from assimilation. Second, the parameterization of the function relating photosynthesis and canopy conductance to ASW resulted in lower photosynthesis

and maximum conductance when ASW was less than 50% of the maximum ASW in the NoExp than in the Base assimilations (Fig. 6b). Finally, the response of photosynthesis to atmospheric CO₂ was functionally zero in the NoExp assimilation, thus highlighting the importance of the elevated CO₂ treatments in the Duke FACE study for constraining the parameterization of the CO₂ response function (Fig. 6c).

## 3.3 Regional predictions with uncertainty

Regionally (i.e., the native range of loblolly pines), stem biomass at age 25 ranged from 52 to 292 Mg ha⁻¹ with the most productive areas located in the coastal plains and the interior of Mississippi and Alabama (Fig. 7a). The least productive locations were the western and northern extents of the native range. The width of the 95% quantile interval for each HUC12 unit ranged from 6.2 to 29.8 Mg ha⁻¹ with the largest uncertainty located in the most productive HUC12 units and in the far western extent of the region (Fig. 7b).

The predicted change in stem biomass at age 25 from an additional 200 ppm of atmospheric CO₂ (over the 1985–2011 concentrations) was similar to the change associated with a
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Regional predictions of stem biomass stocks for a 25-year-old stand planted in 1985. Parameters used in the predictions were from the Base assimilation approach described in Table 5. (b) The width of the 95% quantile interval associated with uncertainty in model parameters.</figcaption></figure>

removal of nutrient limitation (by setting FR to 1) (Fig. 8a, c). The median change associated with elevated CO₂ for a given HUC12 unit ranged from 19.2 to 55.7% with a regional median of 21.7% (Fig. 8a). The change associated with the removal of nutrient limitation ranged from 6.9 to 303.7% for a given HUC12 unit, with a regional median of 24.1% (Fig. 8b). The response to elevated CO₂ was more consistent across space than the response to nutrient addition. The largest potential gains in productivity from nutrient addition were predicted in central Georgia, the northern extent of the region, and the western extents, areas with the lowest SI (Fig. 3).

Stem biomass was considerably less responsive to a 30% decrease in precipitation than to nutrient addition and an increase in atmospheric CO₂. The median change in stem biomass when precipitation was reduced from the 1985–2011 levels ranged from -11.6 to -0.1% for a given HUC12 unit with a regional median of -5.1% (Fig. 8c). Central Georgia was the most responsive to precipitation reduction, reflecting the relatively low annual precipitation and warm temperatures (Fig. 3).

For a given location, the predicted response to elevated CO₂ had larger uncertainty than the predicted response to precipitation reduction and nutrient limitation removal (Fig. 8c, d, f). The uncertainty, defined as the width of the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Predictions of the percentage change in stem biomass at age 25 in response to (a, b) a 200 ppm increase in atmospheric CO₂ over 1985–2011 concentrations, (c, d) a 30% reduction in precipitation from 1985–2011 levels, and (e, f) a removal of nutrient limitation by setting the soil fertility rating in the model equal to 1. The left column is the median prediction and the right column is the width of the 95% quantile interval (C.I.: credible interval) associated with parameter uncertainty. The predictions used the Base assimilation.</figcaption></figure>

95% quantile interval, was consistent across the region for the response to elevated CO₂ (Fig. 8b). The uncertainty in the response to precipitation reduction and nutrient limitation removal was largest in the regions with the largest predicted change (Fig. 8d, f).

## 4 Discussion

Using DA to parameterize models for predicting ecosystem change requires disentangling the vegetation responses to temperature, precipitation, nutrients, and elevated CO₂. To address this challenge, we introduced a regional-scale hierarchical Bayesian approach (Data Assimilation to Predict Productivity for Ecosystems and Regions, Dapper) that assimilated data across environmental gradients and ecosystem
<table><caption>**Table 4.** Description of the different data assimilation approaches used.</caption><thead><tr><td>Simulation name</td><td>Treatments included in assimilation</td><td>Number of plots</td></tr></thead><tbody><tr><td>Base</td><td>All plots and experiments in the region were used simultaneously. Includes unique pCRS, wSx1000, and ThinPower parameters for plots in the Duke FACE study.</td><td>294</td></tr><tr><td>NoExp</td><td>Same as Base assimilation but excluding all plots with experimental manipulations. Includes control plots that are part of experimental studies.</td><td>208</td></tr><tr><td>NoDkPars</td><td>Same as Base assimilation but without pCRS, wSx1000, and ThinPower parameter for plots in the Duke FACE and US-DK3 studies.</td><td>294</td></tr></tbody></table>

manipulation experiments into a modified version of the 3-PG model. Furthermore, we synthesized observations of carbon stocks, carbon fluxes, water fluxes, vegetation structure, and vegetation dynamics that spanned 35 years of forest research in a region (Table 1, Fig. 1) with large and dynamic carbon fluxes (Lu et al., 2015). By combining the DAPPER system with the regional set of observations, we were able to estimate parameters in a model with high predictive capacity (Fig. 4) and with quantified uncertainty on parameters (Table 5) and regional simulations (Figs. 7 and 8).

Our hierarchical approach (Eq. 7) was designed to partition uncertainty among parameters, model process, and measurements (Hobbs and Hooten, 2015). Separating the parameter and process uncertainty is required to estimate prediction intervals, as prediction intervals only include parameter and process errors (Dietze et al., 2013; Hobbs and Hooten, 2015). Previous forest ecosystem DA efforts have either focused on parameter uncertainty, by using measurement uncertainty as the variance term in a Gaussian cost function (Bloom and Williams, 2015; Keenan et al., 2012; Richardson et al., 2010) or on total uncertainty by directly estimating the Gaussian variance term (Ricciuto et al., 2008). Our approach allowed the estimation of the probability distribution of forest biomass before uncertainty is added through measurement. Considering that the method of DA can potentially have a large influence on posterior parameter distributions (Trudinger et al., 2007), future research should focus on comparing the hierarchical approach presented here to other approaches by using the same data constraints with alternative cost functions.

## 4.1 Sensitivity to the inclusion of ecosystem experiments

The most important experimental manipulation for constraining model parameters was the Duke FACE CO₂ fertilization study because the CO₂ fertilization parameters (fCalpha700 and fCpFS700) converged on the lower bounds of their prior distributions when the experiments were excluded from the assimilation. In contrast, excluding the nutrient fertilization, drought, and irrigation studies did not substantially alter the predictive capacity of the model. This finding suggests that data assimilation using plots across environmental gradients alone can constrain parameters associated with water and nutrient sensitivity. However, regardless of whether the experiments were included in the assimilation, the optimized model predicted higher sensitivity to drought than observed, highlighting that future studies should focus on improving the sensitivity to drought.

The 3-PG model included a highly simplified representation of interactions between the water and carbon cycles that resulted in parameterizations that may contain assumptions that require additional investigation. First, transpiration was modeled as a function of a potential canopy transpiration that occurred if leaf area was not limiting transpiration. The LAI at which leaf area was no longer limiting was a parameter that was optimized (LAIgcx in Table 5), resulting in a value of 2.2. Interestingly, this optimized value is consistent with the scant literature on this topic. In their analysis of multiyear measurements of transpiration in loblolly pine, Phillips and Oren (2001) observed that transpiration per unit leaf area was relatively insensitive to increases in leaf area above an LAI of approximately 2.5. Iritz and Lindroth (1996) reviewed transpiration data from a range of crop species and found only small increases in transpiration above LAI of 3–4. These authors suggest that the threshold-type responses observed were related to the range of LAI at which self-shading increases most rapidly, therefore limiting increases in transpiration. The resulting model behavior of “flat” transpiration above 2.2 LAI, with gradually decreasing photosynthesis above that value, results in increasing water use efficiency at higher LAI values. Second, the relationship between relative ASW and the modifier of photosynthesis and transpiration predicted a modifier value greater than zero when the relative ASW was zero. This resulted in positive values from photosynthesis and transpiration when the average ASW during the month was zero. In practice, the monthly ASW was rarely zero during simulations, which presents a challenge constraining the shape of the ASW modifier. The priors for the two ASW modifiers (SWconst and SWpower) had ranges that permitted the modifier to be zero. Therefore, additional data are likely needed during very dry conditions to develop a more physically based parameterization. Alternatively, the parameterization of a non-zero soil moisture modifier at zero ASW may be due to trees having access to water at soil depths deeper than the top 1.5 m of soil represented by the
<table><caption>**Table 5.** The optimized medians, range of the 99 % quantile intervals of the posterior distributions and the 99 % quantile range for priors with normally distributed priors or the range of the upper and lower bounds for priors with uniform distributions. C.I.: credible interval.</caption><thead><tr><th>Parameter</th><th>Posterior median</th><th>Posterior 99% C.I. range</th><th>Prior range</th><th>Posterior/ prior range</th></tr></thead><tbody><tr><td>Allocation and structure</td><td></td><td></td><td></td><td>Parameter group mean: 0.38</td></tr><tr><td>pFS2</td><td>0.58</td><td>0.55-0.61</td><td>0.08-1.00</td><td>0.06</td></tr><tr><td>pFS20</td><td>0.57</td><td>0.55-0.59</td><td>0.10-1.00</td><td>0.05</td></tr><tr><td>pR</td><td>0.11</td><td>0.07-0.15</td><td>0.05-2.00</td><td>0.04</td></tr><tr><td>pCRS</td><td>0.26</td><td>0.25-0.27</td><td>0.15-0.35</td><td>0.11</td></tr><tr><td>pCRS (Duke)</td><td>0.21</td><td>0.18-0.23</td><td>0.15-0.35</td><td>0.20</td></tr><tr><td>SLA0</td><td>8.44</td><td>7.67-9.25</td><td>4.4-6.66</td><td>0.70</td></tr><tr><td>SLA1</td><td>2.84</td><td>2.72-2.96</td><td>3.59-4.16</td><td>0.43</td></tr><tr><td>tSLA</td><td>4.13</td><td>3.88-4.41</td><td>0.43-11.51</td><td>0.05</td></tr><tr><td>fCpFS700</td><td>0.74</td><td>0.60-0.90</td><td>0.50-1.00</td><td>0.60</td></tr><tr><td>StemConst</td><td>0.022</td><td>0.009-0.035</td><td>0.009-0.035</td><td>1.00</td></tr><tr><td>StemPower</td><td>2.78</td><td>2.29-3.27</td><td>2.25-3.29</td><td>0.95</td></tr><tr><td>Canopy photosynthesis, autotrophic respiration, and transpiration</td><td></td><td></td><td></td><td>Parameter group mean: 0.14</td></tr><tr><td>alpha</td><td>0.029</td><td>0.026-0.031</td><td>0.02-0.06</td><td>0.14</td></tr><tr><td>y</td><td>0.50</td><td>0.47-0.53</td><td>0.30-0.65</td><td>0.15</td></tr><tr><td>MaxCond</td><td>0.011</td><td>0.01-0.012</td><td>0.005-0.03</td><td>0.09</td></tr><tr><td>LAIgcx</td><td>2.2</td><td>2.0-2.48</td><td>2.0-5.0</td><td>0.16</td></tr><tr><td>Environmental modifiers of photosynthesis and transpiration</td><td></td><td></td><td></td><td>Parameter group mean: 0.61</td></tr><tr><td>kF</td><td>0.16</td><td>0.12-0.2</td><td>0.14-0.22</td><td>1.04</td></tr><tr><td>T<sub>min</sub></td><td>-5.56</td><td>-8.88 to -2.69</td><td>-1.15 to 9.15</td><td>0.60</td></tr><tr><td>T<sub>opt</sub></td><td>23.42</td><td>21.1-26.31</td><td>19.85-30.15</td><td>0.51</td></tr><tr><td>T<sub>max</sub></td><td>39.56</td><td>34.71-44.39</td><td>32.85-43.15</td><td>0.94</td></tr><tr><td>SWconst</td><td>1.09</td><td>0.91-1.56</td><td>0.01-1.8</td><td>0.36</td></tr><tr><td>SWpower</td><td>8.86</td><td>3.39-12.98</td><td>1.00-13.00</td><td>0.80</td></tr><tr><td>CoeffCond</td><td>0.036</td><td>0.029-0.043</td><td>0.034-0.048</td><td>0.91</td></tr><tr><td>fAlpha700</td><td>1.33</td><td>1.18-1.52</td><td>1.0-1.80</td><td>0.43</td></tr><tr><td>MaxAge</td><td>151.5</td><td>54.4-199.6</td><td>16.0-200.0</td><td>0.79</td></tr><tr><td>nAge</td><td>3.35</td><td>1.77-3.99</td><td>1.00-4.00</td><td>0.74</td></tr><tr><td>rAge</td><td>2.25</td><td>0.81-2.99</td><td>0.01-3.00</td><td>0.73</td></tr><tr><td>FR1</td><td>0.073</td><td>0.061-0.086</td><td>0.00-1.00</td><td>0.03</td></tr><tr><td>FR2</td><td>0.17</td><td>0.15-0.19</td><td>0.0-1.0</td><td>0.04</td></tr><tr><td>Mortality</td><td></td><td></td><td></td><td>Parameter group mean: 0.37</td></tr><tr><td>wSx1000</td><td>176.9</td><td>169.6-184.4</td><td>165.6-294.4</td><td>0.15</td></tr><tr><td>wSx1000 (Duke)</td><td>243.3</td><td>196.89-305.02</td><td>165.6-294.4</td><td>0.76</td></tr><tr><td>ThinPower</td><td>1.68</td><td>1.60-1.78</td><td>1.00-2.5</td><td>0.12</td></tr><tr><td>ThinPower (Duke)</td><td>1.26</td><td>1.00-1.85</td><td>1.00-2.5</td><td>0.56</td></tr><tr><td>mS</td><td>0.52</td><td>0.37-0.71</td><td>0.10-1.00</td><td>0.38</td></tr><tr><td>Rttover</td><td>0.023</td><td>0.017-0.031</td><td>0.017-0.042</td><td>0.55</td></tr><tr><td>MortRate</td><td>0.001</td><td>9e-04-0.0011</td><td>2e-04-0.004</td><td>0.06</td></tr><tr><td>Understory hardwoods</td><td></td><td></td><td></td><td>Parameter group mean: 0.28</td></tr><tr><td>alpha_h</td><td>0.02</td><td>0.02-0.02</td><td>0.005-0.07</td><td>0.01</td></tr><tr><td>pFS_h</td><td>1.78</td><td>1.54-2.06</td><td>0.2-3.0</td><td>0.19</td></tr><tr><td>pR_h</td><td>0.21</td><td>0.06-0.43</td><td>0.05-2.00</td><td>0.19</td></tr><tr><td>SLA_h</td><td>16.3</td><td>14.1-19.0</td><td>6.2-25.8</td><td>0.25</td></tr><tr><td>fAlpha700_h</td><td>1.84</td><td>1.58-2.17</td><td>1.0-2.50</td><td>0.74</td></tr></tbody></table>
<table><caption>**Table 6.** Median and range of the 99% quantile intervals of the posterior distributions for the parameters in the NoExp and NoDkPars assimilations</caption><thead><tr><th>Parameter</th><th>NoExp median</th><th>NoExp 99% range</th><th>NoDkPars median</th><th>NoDkPar 99%</th></tr></thead><tbody><tr><td colspan="5">Allocation and structure</td></tr><tr><td>pFS2</td><td>0.63</td><td>0.61-0.68</td><td>0.57</td><td>0.55-0.60</td></tr><tr><td>pFS20</td><td>0.63</td><td>0.60-0.65</td><td>0.57</td><td>0.55-0.59</td></tr><tr><td>pR</td><td>0.11</td><td>0.06-0.16</td><td>0.11</td><td>0.08-0.15</td></tr><tr><td>pCRS</td><td>0.29</td><td>0.27-0.30</td><td>0.26</td><td>0.25-0.27</td></tr><tr><td>pCRS (Duke)</td><td>0.25</td><td>0.23-0.28</td><td>n/a</td><td>n/a</td></tr><tr><td>SLA0</td><td>7.47</td><td>6.57-8.41</td><td>8.56</td><td>7.73-9.32</td></tr><tr><td>SLA1</td><td>3.00</td><td>2.88-3.12</td><td>2.89</td><td>2.79-2.99</td></tr><tr><td>tSLA</td><td>4.75</td><td>4.30-5.26</td><td>4.12</td><td>3.90-4.38</td></tr><tr><td>fCpFS700</td><td>0.50</td><td>0.50-0.53</td><td>0.94</td><td>0.83-1.00</td></tr><tr><td>StemConst</td><td>0.022</td><td>0.01-0.04</td><td>0.02</td><td>0.01-0.04</td></tr><tr><td>StemPower</td><td>2.79</td><td>2.27-3.26</td><td>2.77</td><td>2.28-3.30</td></tr><tr><td colspan="5">Canopy photosynthesis, autotrophic respiration, and transpiration</td></tr><tr><td>alpha</td><td>0.030</td><td>0.028-0.033</td><td>0.029</td><td>0.026-0.031</td></tr><tr><td>y</td><td>0.48</td><td>0.45-0.51</td><td>0.49</td><td>0.46-0.52</td></tr><tr><td>MaxCond</td><td>0.017</td><td>0.015-0.021</td><td>0.011</td><td>0.011-0.012</td></tr><tr><td>LAIgcx</td><td>4.4</td><td>3.9-5.0</td><td>2.1</td><td>2.0-2.5</td></tr><tr><td colspan="5">Environmental modifiers of photosynthesis and transpiration</td></tr><tr><td>kF</td><td>0.15</td><td>0.11-0.20</td><td>0.16</td><td>0.11-0.20</td></tr><tr><td>T<sub>min</sub></td><td>-7.8</td><td>-10.97 to -4.95</td><td>-6.04</td><td>-9.06 to -3.03</td></tr><tr><td>T<sub>opt</sub></td><td>21.55</td><td>19.15-24.39</td><td>22.71</td><td>20.54-25.42</td></tr><tr><td>T<sub>max</sub></td><td>40.56</td><td>36.51-45.62</td><td>39.82</td><td>35.62-44.56</td></tr><tr><td>SWconst</td><td>0.93</td><td>0.8-1.1</td><td>1.14</td><td>0.91-1.62</td></tr><tr><td>SWpower</td><td>6.27</td><td>2.98-11.49</td><td>7.99</td><td>3.29-12.95</td></tr><tr><td>CoeffCond</td><td>0.041</td><td>0.034-0.047</td><td>0.036</td><td>0.030-0.042</td></tr><tr><td>fAlpha700</td><td>1.01</td><td>1.0-1.06</td><td>1.15</td><td>1.10-1.25</td></tr><tr><td>MaxAge</td><td>152.84</td><td>54.18-199.5</td><td>152.0</td><td>49.2-199.3</td></tr><tr><td>nAge</td><td>3.36</td><td>1.93-3.99</td><td>3.36</td><td>1.89-3.99</td></tr><tr><td>rAge</td><td>2.26</td><td>0.80-2.99</td><td>2.24</td><td>0.83-2.99</td></tr><tr><td>FR1</td><td>0.12</td><td>0.09-0.14</td><td>0.08</td><td>0.07-0.09</td></tr><tr><td>FR2</td><td>0.20</td><td>0.16-0.24</td><td>0.17</td><td>0.15-0.19</td></tr><tr><td colspan="5">Mortality</td></tr><tr><td>wSx1000</td><td>191.6</td><td>180.2-210.2</td><td>181.32</td><td>173.26-196.32</td></tr><tr><td>wSx1000 (Duke)</td><td>235.1</td><td>175.0-297.5</td><td>n/a</td><td>n/a</td></tr><tr><td>ThinPower</td><td>1.76</td><td>1.61-1.92</td><td>1.59</td><td>1.46-1.72</td></tr><tr><td>ThinPower (Duke)</td><td>1.42</td><td>1.01-2.02</td><td>n/a</td><td>n/a</td></tr><tr><td>mS</td><td>0.54</td><td>0.33-0.80</td><td>0.50</td><td>0.25-0.71</td></tr><tr><td>Rttover</td><td>0.019</td><td>0.02-0.03</td><td>0.022</td><td>0.017-0.030</td></tr><tr><td>MortRate</td><td>0.0013</td><td>0.0011-0.0014</td><td>0.0011</td><td>9e-04-0.0013</td></tr><tr><td colspan="5">Understory hardwoods</td></tr><tr><td>alpha_h</td><td>0.031</td><td>0.025-0.040</td><td>0.02</td><td>0.017-0.023</td></tr><tr><td>pFS_h</td><td>2.39</td><td>1.86-2.96</td><td>1.79</td><td>1.59-2.09</td></tr><tr><td>pR_h</td><td>0.25</td><td>0.05-0.67</td><td>0.21</td><td>0.06-0.41</td></tr><tr><td>SLA_h</td><td>12.37</td><td>9.96-15.07</td><td>16.42</td><td>14.37-18.55</td></tr><tr><td>fAlpha700_h</td><td>1.08</td><td>1.00-1.83</td><td>1.83</td><td>1.56-2.15</td></tr></tbody></table>

n/a: not applicable; NoDkPars assimilation did not include Duke-specific parameters.
bucket in 3-PG. Overall, it is important to view the parame-
terization presented here as a phenomenological relationship
that is consistent with observations from drought and irri-
gation experiments as well as observations across regional
gradients in precipitation.

Constraining the sensitivity to atmospheric CO2 differs
from constraining the sensitivity to ASW because, unlike the
multiple constraints on water sensitivity (drought, irrigation,
and gradient studies), environmental conditions created by
the few elevated CO2 plots provided unique constraint on
parameters. Our finding demonstrated that DA efforts should
test for bias in unique ecosystem experiments before finaliz-
ing a set of model parameters used in optimization. In par-
ticular, we found that the parameter governing the photosyn-
thetic response to elevated CO2 (fCalpha700) was substan-
tially lower when all parameters were assumed to be shared
across all plots than when the CO2 fertilization experiment
was allowed to have unique parameters. The need for the
three unique parameters at the Duke FACE study parameters
can be explained by the constraint provided by multiple data
streams and multiple plots. An assumption of the model was
that an increase in stem biomass caused a decrease in stem
density through self-thinning, unless the average tree stem
biomass was below a parameterized threshold (WSx1000).
Therefore, an increase in photosynthesis and stem biomass
through CO2 fertilization could cause a decrease in stem
density. For a single study, it is straightforward to simulta-
neously fit the CO2 fertilization and self-thinning parame-
ters to fit stem biomass and stem density observations for
the site. However, regional DA presents a challenge because
the self-thinning parameters are well constrained by the stem
biomass and stem density observations across the region but
the CO2 fertilization parameters are not. As a result of the
regional DA, the self-thinning parameters caused a stronger
decrease in stem density than observed in the Duke FACE
study. Therefore, the optimization favored a solution where
there was a lower response to CO2 and thus a smaller de-
crease in stem density. Allowing the Duke FACE study to
have unique self-thinning parameters resulted in lower rates
of self-thinning and allowed for simulated stem biomass to
respond to CO2 in a way that matched the observations with-
out penalizing the optimization by degrading the fit to the
stem density.

Our finding that the Duke FACE study required unique
self-thinning parameters to reduce bias in the simulated stem
biomass suggests that when using DA to optimize parame-
ters that are shared across plots, careful examination of pre-
diction bias in key sites that provide a unique constraint on
certain parameters (like the Duke FACE) is critical. Based
on this example, we suggest that DA efforts using multiple
studies and multiple experiment types identify whether par-
ticular experiments at a limited number of sites have the po-
tential to uniquely constrain specific parameters. In this case,
additional weight or site-specific parameters may be needed
to avoid having the signal of the unique experiment over-

whelmed by the large amount of data from the other sites and
experiments. Additionally, the finding suggests that multisite
DA should consider using hierarchical approaches to predict-
ing mortality, particularly because mortality is often not sim-
ulated as mechanistically as growth. A hierarchical approach,
where each plot has a set of mortality parameters that are
drawn from a regional distribution, could avoid having unex-
plained variation in mortality rates leading to bias in the pa-
rameterization of growth-related processes (i.e., growth re-
sponses to CO2, drought, nutrient fertilization). The hierar-
chical approach to mortality could also highlight patterns in
mortality rates across a region and allow for additional inves-
tigations into the mechanisms driving the patterns.

## 4.2 Regional predictions with uncertainty

Our predictions of how stem biomass responds to elevated CO2, nutrient addition, and drought were designed to illustrate the capacity of the DAPPER approach to simulate the uncertainty in future predictions. By using DA, our regional predictions and the uncertainty are consistent with observations but are associated with key caveats. First, only parameter uncertainty was presented in the regional simulations. There is additional uncertainty associated with model process error. We showed the parameter uncertainty because it isolated the capacity to parameterize the individual environmental response functions in the model. Second, the response to drought may be too strong because of the bias in the model predictions of the drought studies. However, there is potential that the drought studies underestimated the sensitivity to ASW since they are relatively short term (<5 years) and manipulate local ASW without manipulating large-scale ASW (i.e., regional water tables). Third, the large responses to nutrient fertilization at the western and northern extents of the study region may be too high. The large responses are attributed to the low SI and the low predicted site fertility rating (FRp). The low SI may be attributable to water limitation and temperature limitation that is not fully accounted for in the parameterization. Additional nutrient addition experiments in the northern and western extent along with further development of the representation of nutrient availability in the 3-PG model may allow for a more robust representation of soil fertility. Finally, the baseline fertility used in our regional analysis was derived from an empirical model of SI that was developed using field plots with minimal management (Sabatia and Burkhart, 2014). Subsequently our estimate of baseline fertility is likely on the low end of forest stands currently in production and the response to nutrient addition may be higher than a typical stand under active management.
## 5 Conclusions

DA is increasingly used for developing predictions from ecosystem models that include uncertainty estimation, due to its ability to represent prior knowledge, integrate observations into the parameterization, and estimate multiple components of uncertainty, including observation, parameter, and process representation uncertainty (Dietze et al., 2013; Luo et al., 2011; Niu et al., 2014). Our application of DA to loblolly pine plantations of the southeastern US demonstrated that these ecosystems are well suited as a test-bed for the development of DA techniques, particularly techniques for assimilating ecosystem experiments. We found that assimilating observations across environmental gradients can provide substantial constraint on many model parameters but that ecosystem manipulative experiments, particularly elevated CO₂ studies, were critical for constraining parameters associated with forest productivity in a more CO₂-enriched atmosphere. This highlights the importance of whole-ecosystem manipulation CO₂ experiments for helping to parameterize and evaluate ecosystem models. Finally, we present an approach for the development of future predictions of forest productivity for natural resource managers that leverage a rich dataset of integrated ecosystem observations across a region.

*Data availability.* Observations used in the DA can be found in the following: the Duke FACE study can be found in McCarthy et al. (2010), the PINEMAP studies are available through the TerraC database (http://terrac.ifas.ufl.edu), the US-DK3 eddy-flux tower data are available through the AmeriFlux database (http://ameriflux.lbl.gov), the Waycross data can be found in Bryars et al. (2013), the US-NC2 data are available upon request from Asko Noormets, and the FMRC and FPC are available through membership with the cooperatives. The parameter chains and 3-PG model code are available upon request from R. Quinn Thomas (rqthomas@vt.edu). SSURGO soils database can be found at https://sdmdataaccess.sc.egov.usda.gov.

**The Supplement related to this article is available online at https://doi.org/10.5194/bg-14-3525-2017-supplement.**

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* Funding support came from USDA-NIFA Project 2015-67003-23485 and the Pine Integrated Network: Education, Mitigation, and Adaptation project (PINEMAP), a Coordinated Agricultural Project funded by the USDA National Institute of Food and Agriculture, Award no. 2011-68002-30185. Additional funding support came from USDA-NIFA McIntire-Stennis Program. The Virginia Space Grant Consortium Graduate

STEM Research Fellowship Program provided partial support for Annika L. Jersild. Computational support was provided by Virginia Tech Advanced Research Computing. This research was also supported by grants from the French Research Agency (MACACC ANR-13-AGRO-0005 and MARIS ANR-14-CE03-0007). We thank Luke Smallman and Mat Williams for helpful discussions about data assimilation, the corporate and government agency members of the FPC and FMRC research cooperatives for supporting the extensive long-term experimental and observational plots in those datasets, and two anonymous reviewers for helpful feedback on the paper. This material is based upon work supported by the US Department of Energy, Office of Science, Office of Biological and Environmental Research, under contract number DE-AC05-00OR22725.

Edited by: Sönke Zaehle

Reviewed by: two anonymous referees

## References

Abatzoglou, J. T.: Development of gridded surface meteorological data for ecological applications and modelling, Int. J. Climatol., 33, 121–131, https://doi.org/10.1002/joc.3413, 2013.

Albaugh, T., Fox, T., Allen, H., and Rubilar, R.: Juvenile southern pine response to fertilization is influenced by soil drainage and texture, Forests, 6, 2799–2819, https://doi.org/10.3390/f6082799, 2015.

Albaugh, T. J., Lee Allen, H., Dougherty, P. M., and Johnsen, K. H.: Long term growth responses of loblolly pine to optimal nutrient and water resource availability, Forest Ecol. Manag., 192, 3–19, https://doi.org/10.1016/j.foreco.2004.01.002, 2004.

Albaugh, T. J., Allen, H. L., and Kress, L. W.: Root and stem partitioning of *Pinus taeda*, Trees, 20, 176–185, https://doi.org/10.1007/s00468-005-0024-4, 2005.

Albaugh, T. J., Albaugh, J. M., Fox, T. R., Allen, H. L., Rubilar, R. A., Trichet, P., Loustau, D., and Linder, S.: Tamm Review: Light use efficiency and carbon storage in nutrient and water experiments on major forest plantation species, Forest Ecol. Manag., 376, 333–342, https://doi.org/10.1016/j.foreco.2016.05.031, 2016.

Allen, C. B., Will, R. E., and Jacobson, M. A.: Production efficiency and radiation use efficiency of four tree species receiving irrigation and fertilization, Forest Sci., 51, 556–569, 2005.

Bartkowiak, S. M., Samuelson, L. J., McGuire, M. A., and Teskey, R. O.: Fertilization increases sensitivity of canopy stomatal conductance and transpiration to throughfall reduction in an 8-year-old loblolly pine plantation, Forest Ecol. Manag., 354, 87–96, https://doi.org/10.1016/j.foreco.2015.06.033, 2015.

Bloom, A. A. and Williams, M.: Constraining ecosystem carbon dynamics in a data-limited world: integrating ecological “common sense” in a model-data fusion framework, Biogeosciences, 12, 1299–1315, https://doi.org/10.5194 bg-12-1299-2015, 2015.

Bryars, C., Maier, C., Zhao, D., Kane, M., Borders, B., Will, R., and Teskey, R.: Fixed physiological parameters in the 3-PG model produced accurate estimates of loblolly pine growth on sites in different geographic regions, Forest Ecol. Manag., 289, 501–514, https://doi.org/10.1016/j.foreco.2012.09.031, 2013.
Burkhart, H. E., Cloeren, D. C., and Amateis, R. L.: Yield relationships in unthinned loblolly pine plantations on cutover, site-prepared lands, South. J. Appl. For., 9, 84–91, 1985.

Carlson, C. A., Fox, T. R., Allen, H. L., Albaugh, T. J., Rubilar, R. A., and Stape, J. L.: Growth responses of loblolly pine in the Southeast United States to midrotation applications of nitrogen, phosphorus, potassium, and micronutrients, Forest Sci., 60, 157–169, https://doi.org/10.5849/forsci.12-158, 2014.

DeLucia, E. H., Drake, J. E., Thomas, R. B., and Gonzalez-Meler, M.: Forest carbon use efficiency: is respiration a constant fraction of gross primary production?, Glob. Change Biol., 13, 1157–1167, https://doi.org/10.1111/j.1365-2486.2007.01365.x, 2007.

Dieguez-Aranda, U., Burkhart, H. E., and Amateis, R. L.: Dynamic site model for loblolly pine (*Pinus taeda* L.) plantations in the United States, Forest Sci., 52, 262–272, 2006.

Dietze, M. C., LeBauer, D. S., and Kooper, R.: On improving the communication between models and data, Plant Cell Environ., 36, 1575–1585, https://doi.org/10.1111/pce.12043, 2013.

Ewers, B. E., Oren, R., Phillips, N., Stromgren, M., and Linder, S.: Mean canopy stomatal conductance responses to water and nutrient availabilities in Picea abies and Pinus taeda, Tree Physiol., 21, 841–850, 2001.

Fox, A., Williams, M., Richardson, A. D., Cameron, D., Gove, J. H., Quaife, T., Ricciuto, D., Reichstein, M., Tomelleri, E., Trudinger, C. M., and Van Wijk, M. T.: The REFLEX project: Comparing different algorithms and implementations for the inversion of a terrestrial ecosystem model against eddy covariance data, Agr. Forest Meteorol., 149, 1597–1615, https://doi.org/10.1016/j.agrformet.2009.05.002, 2009.

Fox, T. R., Jokela, E. J., and Allen, H. L.: The development of pine plantation silviculture in the Southern United States, J. Forest, 105, 337–347, 2007.

Gonzalez-Benecke, C. A. and Martin, T. A.: Water availability and genetic effects on water relations of loblolly pine (*Pinus taeda*) stands, Tree Physiol., 30, 376–392, https://doi.org/10.1093/treephys/tpp118, 2010.

Gonzalez-Benecke, C. A., Gezan, S. A., Albaugh, T. J., Allen, H. L., Burkhart, H. E., Fox, T. R., Jokela, E. J., Maier, C. A., Martin, T. A., Rubilar, R. A., and Samuelson, L. J.: Local and general above-stump biomass functions for loblolly pine and slash pine trees, Forest Ecol. Manag., 334, 254–276, https://doi.org/10.1016/j.foreco.2014.09.002, 2014.

Gonzalez-Benecke, C. A., Teskey, R. O., Martin, T. A., Jokela, E. J., Fox, T. R., Kane, M. B., and Noormets, A.: Regional validation and improved parameterization of the 3-PG model for *Pinus taeda* stands, Forest Ecol. Manag., 361, 237–256, https://doi.org/10.1016/j.foreco.2015.11.025, 2016.

Hobbs, N. T. and Hooten, M. B.: Bayesian Models: A Statistical Primer for Ecologists, Princeton University Press, Princeton, NJ, USA, 2015.

Iritz, Z. and Lindroth, A.: Energy partitioning in relation to leaf area development of short-rotation willow coppice, Agr. Forest Meteorol., 81, 119–130, https://doi.org/10.1016/0168-1923(95)02306-2, 1996.

Keenan, T. F., Carbone, M. S., Reichstein, M., and Richardson, A. D.: The model-data fusion pitfall: assuming certainty in an uncertain world, Oecologia, 167, 587–597, https://doi.org/10.1007/s00442-011-2106-x, 2011.

Keenan, T. F., Davidson, E., Moffat, A. M., Munger, W., and Richardson, A. D.: Using model-data fusion to interpret past trends, and quantify uncertainties in future projections, of terrestrial ecosystem carbon cycling, Glob. Change Biol., 18, 2555–2569, https://doi.org/10.1111/j.1365-2486.2012.02684.x, 2012.

Landsberg, J. and Waring, R.: A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning, Forest Ecol. Manag., 95, 209–228, https://doi.org/10.1016/S0378-1127(97)00026-1, 1997.

Le Quéré, C., Moriarty, R., Andrew, R. M., Canadell, J. G., Sitch, S., Korsbakken, J. I., Friedlingstein, P., Peters, G. P., Andres, R. J., Boden, T. A., Houghton, R. A., House, J. I., Keeling, R. F., Tans, P., Arneth, A., Bakker, D. C. E., Barbero, L., Bopp, L., Chang, J., Chevallier, F., Chini, L. P., Ciais, P., Fader, M., Feely, R. A., Gkritzalis, T., Harris, I., Hauck, J., Ilyina, T., Jain, A. K., Kato, E., Kitidis, V., Klein Goldewijk, K., Koven, C., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lenton, A., Lima, I. D., Metzl, N., Millero, F., Munro, D. R., Murata, A., Nabel, J. E. M. S., Nakaoka, S., Nojiri, Y., O'Brien, K., Olsen, A., Ono, T., Pérez, F. F., Pfeil, B., Pierrot, D., Poulter, B., Rehder, G., Rödenbeck, C., Saito, S., Schuster, U., Schwinger, J., Séférían, R., Steinhoff, T., Stocker, B. D., Sutton, A. J., Takahashi, T., Tilbrook, B., van der Laan-Luijkx, I. T., van der Werf, G. R., van Heuven, S., Vandenmark, D., Viovy, N., Wiltshire, A., Zaehle, S., and Zeng, N.: Global Carbon Budget 2015, Earth Syst. Sci. Data, 7, 349–396, https://doi.org/10.5194/essd-7-349-2015, 2015.

LeBauer, D. S., Dietze, M., Kooper, R., Long, S., Mulrooney, P., Rohde, G. S., and Wang, D.: Energy Biosciences Institute, University of Illinois at Urbana-Champaign, Biofuel Ecophysiological Traits and Yields Database (BETYdb), available at: https://www.betydb.org (last access: 16 May 2016), 2010.

Lu, X., Lu, X., Kicklighter, D. W., Kicklighter, D., Melillo, J. M., Melillo, J. M., Reilly, J. M., Reilly, J. M., Xu, L., and Wu, L.: Land carbon sequestration within the conterminous United States: Regional- and state-level analyses, J. Geophys. Res.-Biogeo., 120, 379–398, https://doi.org/10.1002/2014JG002818, 2015.

Luo, Y., Weng, E., Wu, X., Gao, C., Zhou, X., and Zhang, L.: Parameter identifiability, constraint, and equifinality in data assimilation with ecosystem models, Ecol. Appl., 19, 571–574, https://doi.org/10.1890/08-0561.1, 2009.

Luo, Y., Ogle, K., Tucker, C., Fei, S., Gao, C., LaDeau, S., Clark, J. S., and Schimel, D. S.: Ecological forecasting and data assimilation in a data-rich era, Ecol. Appl., 21, 1429–1442, https://doi.org/10.1890/09-1275.1, 2011.

MacBean, N., Peylin, P., Chevallier, F., Scholze, M., and Schürmann, G.: Consistent assimilation of multiple data streams in a carbon cycle data assimilation system, Geosci. Model Dev., 9, 3569–3588, https://doi.org/10.5194/gmd-9-3569-2016, 2016.

Matamala, R., González-Meler, M. A., Jastrow, J. D., Norby, R. J., and Schlesinger, W. H.: Impacts of fine root turnover on forest NPP and soil C sequestration potential, Science, 302, 1385–1387, https://doi.org/10.1126/science.1089543, 2003.

McCarthy, H. R., Oren, R., Johnsen, K. H., Gallet-Budynek, A., Pritchard, S. G., Cook, C. W., LaDeau, S. L., Jackson, R. B., and Finzi, A. C.: Re-assessment of plant carbon dynamics at the Duke free-air CO₂ enrichment site: interactions of atmospheric [CO₂] with nitrogen and water avail-
ability over stand development, New Phytol., 185, 514–528,
https://doi.org/10.1111/j.1469-8137.2009.03078.x, 2010.

McKeand, S., Mullin, T., Byram, T., and White, T.: Deployment of genetically improved loblolly and slash pines in the south, J. Forest., 101, 32–37, 2003.

Medlyn, B. E., Zaehle, S., De Kauwe, M. G., Walker, A. P., Dietze, M. C., Hanson, P. J., Hickler, T., Jain, A. K., Luo, Y., Parton, W., Prentice, I. C., Thornton, P. E., Wang, S., Wang, Y.-P., Weng, E., Iversen, C. M., McCarthy, H. R., Warren, J. M., Oren, R., and Norby, R. J.: Using ecosystem experiments to improve vegetation models, Nature Climate Change, 5, 528–534,
https://doi.org/10.1038/nclimate2621, 2015.

Niu, S., Luo, Y., Dietze, M. C., Keenan, T. F., Shi, Z., Li, J., and Chapin III, F. S.: The role of data assimilation in predictive ecology, Ecosphere, 5, 65, https://doi.org/10.1890/ES13-00273.1, 2014.

Noormets, A., Gavazzi, M. J., McNulty, S. G., Domec, J.-C., Sun, G., King, J. S., and Chen, J.: Response of carbon fluxes to drought in a coastal plain loblolly pine forest, Glob. Change Biol., 16, 272–287, https://doi.org/10.1111/j.1365-2486.2009.01928.x, 2010.

Novick, K. A., Oishi, A. C., Ward, E. J., Siqueira, M. B. S., Juang, J.-Y., and Stoy, P. C.: On the difference in the net ecosystem exchange of CO₂ between deciduous and evergreen forests in the southeastern United States, Glob. Change Biol., 21, 827–842,
https://doi.org/10.1111/gcb.12723, 2015.

Oren, R., Ellsworth, D., Johnsen, K., Phillips, N., Ewers, B., Maier, C., Schafer, K., McCarthy, H., Hendrey, G., McNulty, S. G., and Katul, G.: Soil fertility limits carbon sequestration by forest ecosystems in a CO₂-enriched atmosphere, Nature, 411, 469–472, https://doi.org/10.1038/35078064, 2001.

Pan, Y., Birdsey, R. A., Fang, J., Houghton, R., Kauppi, P. E., Kurz, W. A., Phillips, O. L., Shvidenko, A., Lewis, S. L., Canadell, J. G., Ciais, P., Jackson, R. B., Pacala, S. W., McGuire, A. D., Piao, S. L., Rautiainen, A., Sitch, S., and Hayes, D.: A large and persistent carbon sink in the world's forests, Science, 333, 988–993, https://doi.org/10.1126/science.1201609, 2011.

Phillips, N. and Oren, R.: Intra- and inter-annual variation in transpiration of a pine forest, Ecol. Appl., 11, 385–396, 2001.

Raymond, J. E., Fox, T. R., Strahm, B. D., and Zerpa, J.: Differences in the recovery of four different nitrogen containing fertilizers after two application seasons in pine plantations across the southeastern United States, Forest Ecol. Manag., 380, 161–171,
https://doi.org/10.1016/j.foreco.2016.08.044, 2016.

Ricciuto, D. M., Davis, K. J., and Keller, K.: A Bayesian calibration of a simple carbon cycle model: The role of observations in estimating and reducing uncertainty, Glob. Biogeochem. Cy., 22, GB2030, https://doi.org/10.1029/2006GB002908, 2008.

Richardson, A. D., Williams, M., Hollinger, D. Y., Moore, D. J. P., Dail, D. B., Davidson, E. A., Scott, N. A., Evans, R. S., Hughes, H., Lee, J. T., Rodrigues, C., and Savage, K.: Estimating parameters of a forest ecosystem C model with measurements of stocks and fluxes as joint constraints, Oecologia, 164, 25–40,
https://doi.org/10.1007/s00442-010-1628-y, 2010.

Sabatia, C. O. and Burkhart, H. E.: Predicting site index of plantation loblolly pine from biophysical variables, Forest Ecol. Manag., 326, 142–156,
https://doi.org/10.1016/j.foreco.2014.04.019, 2014.

Samuelson, L. J., Butnor, J., Maier, C., Stokes, T. A., Johnsen, K., and Kane, M.: Growth and physiology of loblolly pine in response to long-term resource management: defining growth potential in the southern United States, Can. J. Forest Res., 38, 721–732, https://doi.org/10.1139/X07-191, 2008.

Shvidenko, A., Barber, C. V., and Persson, R.: Forest and Woodland Systems, in: Ecosystems and Human Well-being: Current State and Trends, Volume, edited by: Hassan, R., Scholes, R., and Ash, N., 585–621, Island Press, Washington, USA, 2005.

Soil Survey Staff: Natural Resources Conservation Service, United States Department of Agriculture, Soil Survey Geographic (SSURGO) Database, available online at: https://sdmdataaccess.sc.egov.usda.gov, last access: 12 November 2013.

Subedi, S., Fox, T., and Wynne, R.: Determination of fertility rating (FR) in the 3-PG model for loblolly pine plantations in the Southeastern United States based on site index, Forests, 6, 3002–3027,
https://doi.org/10.3390/f6093002, 2015.

Tang, Z., Sayer, M. A. S., Chambers, J. L., and Barnett, J. P.: Interactive effects of fertilization and throughfall exclusion on the physiological responses and whole-tree carbon uptake of mature loblolly pine, Can. J. Botany, 82, 850–861,
https://doi.org/10.1139/b04-064, 2004.

Trudinger, C. M., Raupach, M. R., Rayner, P. J., Kattge, J., Liu, Q., Pak, B., Reichstein, M., Renzullo, L., Richardson, A. D., Roxburgh, S. H., Styles, J., Wang, Y.-P., Briggs, P., Barrett, D., and Nikolova, S.: OptIC project: An intercomparison of optimization techniques for parameter estimation in terrestrial biogeochemical models, J. Geophys. Res., 112, G02027–17,
https://doi.org/10.1029/2006JG000367, 2007.

Ward, E. J., Domec, J.-C., Laviner, M. A., Fox, T. R., Sun, G., McNulty, S., King, J., and Noormets, A.: Fertilization intensifies drought stress: Water use and stomatal conductance of Pinus taeda in a midrotation fertilization and throughfall reduction experiment, Forest Ecol. Manag., 355, 72–82,
https://doi.org/10.1016/j.foreco.2015.04.009, 2015.

Weng, E. and Luo, Y.: Relative information contributions of model vs. data to short- and long-term forecasts of forest carbon dynamics, Ecol. Appl., 21, 1490–1505, https://doi.org/10.1890/09-1394.1, 2011.

Will, R., Fox, T., Akers, M., Domec, J.-C., González-Benecke, C., Jokela, E., Kane, M., Laviner, M., Lokuta, G., Markewitz, D., McGuire, M., Meek, C., Noormets, A., Samuelson, L., Seiler, J., Strahm, B., Teskey, R., Vogel, J., Ward, E., West, J., Wilson, D., and Martin, T.: A range-wide experiment to investigate nutrient and soil moisture interactions in loblolly pine plantations, Forests, 6, 2014–2028, https://doi.org/10.3390/f6062014, 2015.

Williams, M., Schwarz, P., Law, B. E., Irvine, J., and Kurpius, M.: An improved analysis of forest carbon dynamics using data assimilation, Glob. Change Biol., 11, 89–105,
https://doi.org/10.1111/j.1365-2486.2004.00891.x, 2005.

Zobitz, J. M., Desai, A. R., Moore, D. J. P., and Chadwick, M. A.: A primer for data assimilation with ecological models using Markov Chain Monte Carlo (MCMC), Oecologia, 167, 599–611,
https://doi.org/10.1007/s00442-011-2107-9, 2011.

Ziehn, T., Scholze, M., and Knorr, W.: On the capability of monte carlo and adjoint inversion techniques to derive posterior parameter uncertainties in terrestrial ecosystem models, Global Bio.-Geochem. Cy., 26, GB3025, doi:10.1029/2011GB004185, 2012.