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

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 Predicting Growth (3-PG) forest growth model. The observations included major experiments that manipulated atmospheric carbon dioxide (CO2) concentration, water, and nutrients, along with nonexperimental surveys that spanned environmental gradients across an 8.6× 105 km2 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 CO2 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 CO2 strongly depended on whether ecosystem experiments were assimilated and whether the assimilated field plots in the CO2 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 CO2, decreased precipPublished by Copernicus Publications on behalf of the European Geosciences Union. 3526 R. Q. Thomas et al.: Leveraging 35 years of Pinus taeda research itation, and increased nutrient availability that include estimates of uncertainty for the southeastern US. Overall, we (1) demonstrated how three decades of research in southeastern US planted pine forests can be used to develop DA techniques 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 predictions of forest productivity for natural resource managers that leverage a rich dataset of integrated ecosystem observations across a region.

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 2 ) concentration, water, and nutrients, along with nonexperimental surveys that spanned environmental gradients across an 8.6 × 10 5 km 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 2 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 2 strongly depended on whether ecosystem experiments were assimilated and whether the assimilated field plots in the CO 2 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 2 , decreased precip-

Introduction
Forest ecosystems absorb and store a large fraction of anthropogenic carbon dioxide (CO 2 ) emissions (Le Quéré et al., 2015;Pan et al., 2011) and supply wood products to a growing human population (Shvidenko et al., 2005).Therefore, predicting future carbon sequestration and timber supply is critical for adapting forest management practices to future environmental conditions and for using forests to assist with the reduction in atmospheric CO 2 concentrations.The key sources of information for developing these predictions are results from global change ecosystem manipulation experiments, observations of forest dynamics across environmental 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 increasingly used framework for integrating ecosystem observations 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, methods that can generate probability distributions for ecosystem model parameters and initial states.DA allows for the explicit accounting of observational uncertainty (Keenan et al., 2011), the incorporation of multiple types of observations with different timescales of collection (MacBean et al., 2016;Richardson et al., 2010), and the representation of prior knowledge through informed parameter prior distributions or specific relationships among parameters (Bloom and Williams, 2015).
Using DA to parameterize ecosystem models with observations from multiple locations that leverage ecosystem manipulation experiments and environmental gradients will allow for predictions to be consistent with the rich history of global change research in forest ecosystems.Ecosystem manipulation experiments provide a controlled environment in which data collected can be used to describe how forests acclimate and operate under altered environmental conditions (Medlyn et al., 2015) and can potentially allow for the optimization of model parameters associated with the altered environmental factor in the experiment.Furthermore, the assimilation 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 availabilities, and soil water dynamics.These studies leverage time and space to quantify the sensitivity of forest dynamics to environmental variation.However, covariation of environmental variation can pose challenges separating the responses to individual environmental factors.Overall, assimilating 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 multiple 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 CO 2 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 CO 2 by nutrient addition treatments (Mc-Carthy 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 multidecadal 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 2 across the southeastern US.

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.4kg ha −1 N + 13.44 kg ha −1 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 −1 N + 27 kg ha −1 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 (SE-TRES) 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 −1 + 17 kg P ha −1 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 2 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 2 fumigated rings (including the unfertilized half of the prototype), two nitrogen fertilization treatments (115 kg N ha −1 yr −1 applied annually), and one CO 2 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 eddycovariance 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).

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 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 2 , 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 2 modifier had a value of 1 at 350 ppm (thus values greater than 1 at higher CO 2 concentrations).Elevated CO 2 modified tree physiology by increasing quantum yield, based on an increasing but saturating relationship with atmospheric CO 2 (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 2 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 www.biogeosciences.net/14/3525/2017/Biogeosciences, 14, 3525-3547, 2017 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  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 eddycovariance 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 plotspecific 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 −1 ), 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 2 , 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 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 2 enrichment experiments were simulated by setting the atmospheric CO 2 input equal to the treatment mean from the elevated CO 2 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.

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, where θ is a vector of parameters that are optimized, FR p is the site fertility, and σ m is the model process error.Not shown are the vector of parameters that were not optimized (Supplemental 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 www.biogeosciences.net/14/3525/2017/Biogeosciences, 14, 3525-3547, 2017 of the prediction: 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 ρ m at 0. FR p did not have an explicit probability distribution.Rather the probability density was evaluated as 1 if the plot was not fertilized, thus causing FR p to be estimated from SI and MAT (Supplement Eq. 15), or if it was a fertilized plot and had an FR p equal or higher than that of its non-fertilized control plot.The probability density was evaluated as 0 if the estimated FR p in a fertilized plot was less than the FR p in the control plot or if FR p was not contained in the interval between 0 and 1.
1 if non-fertilized, FR p ≥ 0, and FR p ≤ 1 1 if FR p = 1 and fertilization levels are assumed to remove nutrient deficiencies 0 if FR p < 1 and fertilization levels are assumed to remove nutrient deficiencies 1 if fertilized but levels are not assumed to remove deficiencies and FR p ≥ FR of control plot 0 if fertilized but levels are not assumed to remove deficiencies and FR p < FR of control plot 0 if FR p < 0 or FR p > 1 (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 where τ 2 p,m,t 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 τ 2 p,m,t 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 (θ ), FR p , and the process error parameters.The prior distributions for (p(θ )) 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.
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 p(θ , y, γ , q|y, τ , priors) ∝, 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.

Data assimilation evaluation
Using the observations, model, and hierarchical Bayesian method described above, we assimilated both the nonmanipulated 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 2 ).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 2 plots in the Duke FACE study.Further analysis investigating the cause of the bias in the CO 2 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).

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 2 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 2 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.

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 2 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 2 increased by 200 ppm.Within a parameter sample, we calculated the percent change in stem www.biogeosciences.net/14/3525/2017/Biogeosciences, 14, 3525-3547, 2017 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.

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 −1 (Fig. 4a).Furthermore, the response of stem biomass to irrigation (df = 7, p = 0.18), nutrient addition (df = 26, p = 0.29), and elevated CO 2 (df = 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).q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Observed (Mg biomass ha q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Observed (Mg biomass ha  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.4).
to elevated CO 2 than observed (df = 4, p < 0.001; Fig. 5).Furthermore, there was a slight negative response of stem biomass to CO 2 in the NoExp assimilation because the parameter governing the change in foliage allocation at elevated CO 2 (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 2 .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 2 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 photosynthe-sis 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 2 was functionally zero in the NoExp assimilation, thus highlighting the importance of the elevated CO 2 treatments in the Duke FACE study for constraining the parameterization of the CO 2 response function (Fig. 6c).

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 −1 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 ha −1 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 2 (over the 1985-2011 concentrations) was similar to the change associated with a removal of nutrient limitation (by setting FR to 1) (Fig. 8a, c).The median change associated with elevated CO 2 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 2 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 2 .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 2 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 95 % quantile interval, was consistent across the region for the response to elevated CO 2 (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).

Discussion
Using DA to parameterize models for predicting ecosystem change requires disentangling the vegetation responses to temperature, precipitation, nutrients, and elevated CO 2 .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  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.

Sensitivity to the inclusion of ecosystem experiments
The most important experimental manipulation for constraining model parameters was the Duke FACE CO 2 fertilization study because the CO 2 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 bucket in 3-PG.Overall, it is important to view the parameterization presented here as a phenomenological relationship that is consistent with observations from drought and irrigation experiments as well as observations across regional gradients in precipitation.
Constraining the sensitivity to atmospheric CO 2 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 CO 2 plots provided unique constraint on parameters.Our finding demonstrated that DA efforts should test for bias in unique ecosystem experiments before finalizing a set of model parameters used in optimization.In particular, we found that the parameter governing the photosynthetic response to elevated CO 2 (fCalpha700) was substantially lower when all parameters were assumed to be shared across all plots than when the CO 2 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 CO 2 fertilization could cause a decrease in stem density.For a single study, it is straightforward to simultaneously fit the CO 2 fertilization and self-thinning parameters 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 CO 2 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 CO 2 and thus a smaller decrease 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 CO 2 in a way that matched the observations without 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 parameters that are shared across plots, careful examination of prediction 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 particular experiments at a limited number of sites have the potential 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 predicting mortality, particularly because mortality is often not simulated 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 unexplained variation in mortality rates leading to bias in the parameterization of growth-related processes (i.e., growth responses to CO 2 , drought, nutrient fertilization).The hierarchical approach to mortality could also highlight patterns in mortality rates across a region and allow for additional investigations into the mechanisms driving the patterns.

Regional predictions with uncertainty
Our predictions of how stem biomass responds to elevated CO 2 , 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 (FR p ).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.

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 2 studies, were critical for constraining parameters associated with forest productivity in a more CO 2 -enriched atmosphere.This highlights the importance of whole-ecosystem manipulation CO 2 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.

Figure 1 .
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 2 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 2 .
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.

Figure 3 .
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.).

Figure 4 .
Figure 4. Model evaluation of stem biomass when assimilating (a) observations across environmental gradients and ecosystem manipulation experiments (Base; Table4) and (b) only observations across environmental gradients (NoExp; Table4).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 5 .
Figure5.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 2 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 6 .
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 2 influence on photosynthesis.The function shapes were derived from the parameters in the Base, NoExp, and NoDkPars assimilations (Table4).

Figure 7 .
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.

Figure 8 .
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 2 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.

Table 1 .
Regional observational data streams used in data assimilation.

Table 2 .
Descriptions of the studies used in data assimilation.
a Forest Modeling Research Cooperative.b Forest Productivity Cooperative.c PINEMAP. d Southeast Tree Research and Education Site.e Free-Air Carbon Enrichment.

Table 4 .
Description of the different data assimilation approaches used.