Modeling forest carbon cycle using long-term carbon stock field measurement in the Delaware River Basin

Process-based models are a powerful approach to test our understanding of biogeochemical processes, to extrapolate ground survey data from limited plots to the landscape scale, and to simulate the effects of climate change, nitrogen deposition, elevated atmospheric CO2, increasing natural disturbances, and land-use change on ecological processes. However, in most studies, the models are calibrated using ground measurements from only a few sites, though they may be extrapolated to much larger areas. Estimation accuracy can be improved if the models are parameterized using long-term carbon (C) stock data from multiple sites representative of the simulated region. In this study, forest biomass C stocks measured in 61 forested plots located in three research sites in the Delaware River Basin (DRB) were used to modify the PnET-CN model in three ways: (1) Field-measured mortality rates in each forest type were used to parameterize the wood turnover rate; (2) a numerical approach was used to calibrate the relationship between foliage N and maximum photosynthesis rate; and (3) stand age was incorporated into the model as an input variable, which determines the year of the last disturbance. The results showed that these model modifications improved model performance in capturing the spatial variation of forest C dynamics in the DRB forests. The spatial distribution of forest C pools and fluxes in the three sites was mapped using the modified model. The modified model was also used in experimental scenarios, which predicted that 39% of forest C sequestered over the past decade could be attributed to the combined effects of elevated CO2 and N deposition. This study demonstrated an effective method for using long-term biometric measurements of forest biomass C stocks to constrain and improve a process-based ecosystem model at a regional scale. Further research should target improving model parameters that are sensitive to the spatial variation of forest C dynamics.


INTRODUCTION
Process-based ecosystem models are a powerful approach to test our understanding of biogeochemical processes, to extrapolate ground survey data from limited plots to the landscape scale, and to simulate the effects of climate change, N deposition, elevated atmospheric CO 2 , increasing natural disturbances, and land-use change on ecological processes (Campbell et al. 2009, Felzer 2012. They can also be used for projecting future ecosystem changes and dynamics based on our current knowledge, which is important for future resource management and actions for climate change adaptation (Chen et al. 2000, Mickler et al. 2002, Hudiburg et al. 2013. In most modeling studies of forest C cycles, models are parameterized or calibrated using carbon fluxes measurements from eddy covariance techniques from one or several flux towers (Carvalhais et al. 2010, Liu et al. 2015. The eddy covariance data sets consist of fine time resolution (hourly observations), but are very low in their spatial resolution (i.e., few sites covering a large area). Therefore, the parameters associated with rapid processes such as photosynthesis are well constrained, but the parameters associated with long-term processes such as the turnover rates of biomass and soil C are poorly constrained (Braswell et al. 2005, Richardson et al. 2010. Moreover, the ability of the model to capture heterogeneity of the forest C and N cycles over complex landscapes is difficult to evaluate using eddy covariance data because the footprint area of each flux tower is normally smaller than 10 km 2 , and a limited number of flux towers are used to represent large areas of forests (Owen et al. 2007).
Over the past decade, a large number of longterm biometric measurements have accumulated data for forest biomass and soil C stocks (Xu et al. 2014). Forest Inventory Analysis (FIA) and remote sensing techniques have provided useful information on the spatial distribution of forest C stock changes (Woodbury et al. 2007, Lichstein et al. 2014. Several studies have indicated that flux and biometric observations can play complementary roles in constraining model parameters (Kondo et al. 2013, Du et al. 2015. The performance of a process-based simulation model could be improved if it can be parameterized using better data sources from long-term biometric C stock data from multiple sites representative of the simulated region (Carvalhais et al. 2010, Molina-Herrera et al. 2015. The PnET models (Aber and Federer 1992) were based on general principles of ecophysiology and biogeochemistry and have been validated in various vegetation types (Goodale et al. 2002, Chen et al. 2004, Thorn et al. 2015. PnET-CN is the model version simulating closed C and N cycles . Applications of the PnET-CN model have provided important insights into the interactive effects of climate change, N deposition, increasing atmospheric CO 2 and ozone, and natural and anthropogenic disturbance on forest C and N cycles , Pan et al. 2004a, Chiang et al. 2008. The PnET-CN model was originally developed for single-site simulations in the northern hardwood forests in the northeast United States . The Delaware River Basin (DRB), situated in the southern edge of the northern hardwood forest, features diverse forest types and land-use histories, representing a transition zone from typical northern hardwood forests to Mid-Atlantic oak-dominated deciduous forests. To use the PnET-CN model effectively in the DRB forest, the model needs to be localized. A multi-agency program known as the Collaborative Environmental Monitoring and Research Initiative (CEMRI) has collected long-term biometric measurement of forest C stocks in the DRB (Murdoch et al. 2008). Recent re-measurements and data analysis of the CEMRI sites (Xu et al. 2016b) provided a good opportunity to examine the PnET-CN model performance outside its original range in the DRB landscapes along a gradient, and to explore approaches for parameterizing and improving the model using long-term field measurements from multiple sites.
In this study, forest biomass C stocks and their change over the recent decade, measured in 61 plots from three sites, were used to constrain the PnET-CN model applying for the DRB forests. The major objectives of this study were (1) to explore multiple parameterization approaches for improving model performance using longterm biometric measurements; (2) to verify the model specifically improved for our study sites and compare model performance between the original and modified models; (3) to apply the modified model to simulate the spatial distribution of C fluxes and pools in the DRB forests; and (4) to simulate the effects of environmental changes such as N deposition and elevated atmospheric CO 2 concentration on the spatial distributions and temporal changes in the DRB forest C cycle. The newly adapted model for DRB is expected to provide more accurate modeling results which are useful to our understanding of regional forest C dynamics and to forest management.

Model description
The PnET-CN model ) is a process-based ecosystem model that uses information on vegetation, climate, and soil to simulate the carbon, nitrogen, and water dynamics of forest ecosystems in the northeastern United States. After years of development, three versions of the model have formed a nested sequence in the PnET family of models: PnET-Day, a daily time step model of forest canopy C balance (Aber et al. 1996); PnET-II, a monthly time step model which incorporates soil respiration and full water balance (Aber et al. 1995); and PnET-CN, a monthly time step model that includes closed carbon and nitrogen cycles and the effects of elevated atmospheric CO 2 and ozone . The PnET-CN model contains a set of modules to simulate major processes of forest carbon, nitrogen, and water cycles. The input data include climate variables (temperature, precipitation, radiation, etc.), site parameters (latitude, elevation, soil water-holding capacity [WHC], etc.), and vegetation parameters (forest type, canopy trait, C allocation, N content, etc.). The major outputs include above-and belowground biomass, net primary productivity (NPP), net ecosystem productivity (NEP), soil C content, N leaching, and water yield.
The PnET-CN model was built on the principle that the maximum photosynthetic rate (Amax) is a function of foliar N concentration. This represents the interaction between C and N cycles and links leaf traits with physiological processes at the ecosystem level. Photosynthesis is further constrained by temperature and water availability, and some photosynthetic product is consumed by respiration while the rest is allocated to living C pools (foliage, wood, and fine root). Living biomass is transferred to dead C pools (e.g., woody debris and soil organic matter) through litterfall, mortality, or disturbance. The model also incorporates the effect of atmospheric N deposition, which will increase the available soil N pool, and elevated CO 2 concentration, which will increase the potential photosynthesis rate and decrease stomatal water conductance.
Unlike most process-based ecosystem models and Dynamic Global Vegetation Models, PnET-CN uses finer vegetation classes down to forest type groups rather than coarse plant functional types (Pan et al. 1996), which gives the model an advantage in making it more suitable for forest ecosystem studies at a regional scale. The PnET-CN model originally included parameters for four vegetation types (northern hardwood, spruce-fir, red oak-red maple, and pine forests). A new forest type (mixed forests) has been created by numerically combining the results of red oak-red maple and pine forest types to represent the mixed forest in the Mid-Atlantic region (Pan et al. 2004b(Pan et al. , 2009. Regional spatially explicit versions of the model based on geographic information systems (GIS version of the PnET-CN model) were developed to extrapolate the simulation from single-plot to the regional level (Pan et al. 2009). More details about the model structure and processes have been described elsewhere , Ollinger et al. 2002, Pan et al. 2004a).

Study sites and field measurement
The Delaware River is one of the major rivers in the Mid-Atlantic region of the United States, draining an area of about 33,000 km 2 in Pennsylvania, New Jersey, New York, Delaware, and Maryland. The Delaware River Basin is characterized by a humid continental climate, with mean annual temperature of 9-12°C and mean annual precipitation of 1143 mm (Kauffman et al. 2008). The DRB is located in the eco-zone of deciduous forests and is ecologically diverse, comprised of five physiographic provinces and multiple species assemblages that represent most of the major eastern U.S. forest types (Murdoch et al. 2008). Three areas in the DRB were selected as intensive monitoring and research sites for process-level studies in forested landscapes: the Neversink River Basin (NS) in the northern, mostly forested region of the Appalachian Plateau Province; the Delaware Water Gap Area (DEWA) with three small watersheds (Adams Creek, Dingmans Falls, and Little Bushkill) lying in the central Appalachian Plateau Province; and the French Creek Watershed (FC) in the midbasin Piedmont Province (Table 1, Fig. 1).
During 2001-2003, 61 intensive inventory plots were randomly located across the three sites.
In each plot, all trees with diameter at breast height (dbh) greater than 5 inches (12.7 cm) were Notes: All data were extracted from the model input GIS database, and mean values for each site are shown. Annual temperature and precipitation are 30-yr means from 1981 to 2010 (Thornton et al. 2014). Wet deposition is inorganic N deposition from 1983 to 2007 (Grimm and Lynch 2004).  measured and marked, and the specific locations of the plots were mapped. In 2012-2014, these plots were revisited, and biomass parameters were re-measured using the same protocols (Xu et al. 2016b). The plot design and sampling method followed the FIA protocols in the two measurements ( Fig. 2; U.S. Department of Agriculture 2015). Each plot had four round subplots, covering a total area of 672.44 m 2 . Live and dead trees were measured in each subplot. dbh, total and bole height, tree species, and status change of each tree were recorded. Each subplot also had one microplot (area, 13.49 m 2 ). Live and dead saplings (1 inch < dbh < 5 inch) and seedlings (dbh < 1 inch) were measured in the microplots. Within each plot, two trees close to the subplots that represented the dominant species and growing condition of the forest stand were selected as site trees. The ages of the site trees were measured by counting rings in a tree core. The stand ages of plots were determined as the mean age of the two site trees. Field measurement data from the two measurement visits (2001-2003 and 2012-2014) were compiled into a single database for biomass C and biomass C change calculations (Cole et al. 2013). Species-specific allometric equations (Jenkins et al. 2004) were used to calculate aboveground tree biomass, and the general equations were used to calculate coarse roots biomass. To compare with model output, each plot was assigned to a vegetation type. The total biomass of living trees, saplings, and seedlings was summed to calculate the observed live biomass in the two measurements. A conversion factor of 0.5 was used to convert biomass to C stock.

Model input data
In this study, the GIS version of the PnET-CN model (Pan et al. 2009), which is spatially explicit, was used in the three intensive monitoring and research sites in the DRB with spatial resolution of 1 9 1 km. The required model input data include elevation, vegetation type, and soil WHC for each pixel, and monthly maximum and minimum air temperature, precipitation, photosynthetically active radiation (PAR), wet NH 4 and NO 3 deposition, atmospheric CO 2 , and ozone concentrations for the entire simulation period. The climatic data used include monthly averages and historical data. The historical data for precipitation and maximum and minimum air temperature for each pixel and each month were extracted from the database developed by Prism Climate Group (http://www.prism.oregon state.edu/). The monthly historical data were available from 1895 to 2010. Photosynthetically active radiation data were derived from Daymet data (http://daymet.ornl.gov/) which were available from 1980 to 2010. The climate data from 1980 to 2010 were used to calculate the 30-yr monthly averages, and the monthly average data were used in the years when the historical data were not available. The ozone data layer was developed by Teague Prichard and US Forest Service based on the SUM60 data from May to September monitored by EPA (personal . The historical data are available from 1993 to 2010. The monthly average ozone deposition was calculated based on all the historical data, and then used when historical data were not available. Ammonium (NH 4 ) and nitrate (NO 3 ) wet deposition data from December 1982 to November 2007 were created by Jeff Grimm based on standardized precipitation grids (Grimm and Lynch 2004). We assumed that N deposition before 1931 was approximately 20% of the average deposition from 1993 to 2007, and the deposition rate increased linearly from 1931 to 1993, and after 2007. This linear relationship was interpolated for each pixel.
In this study, the vegetation type data were reclassified from USFS forest type groups (http://data.fs.usda.gov/geodata/rastergateway/ forest_type/conus_forest_type_group_metadata. php). Based on the map of forest type groups, the NS area is dominated by northern hardwood forest; FC area is dominated by red oak-red maple forest; DEWA is located in a transition area between these two vegetation types, with a small proportion of the area covered by mixed forest. Soil WHC data were derived from STATSGO soil database (http://websoilsurvey.nrcs.usda.gov/; Soil Survey Staff 1995). To run the GIS version of the PnET-CN model with a standardized pixel size, all data layers of different spatial resolutions were adjusted to 1 9 1 km.

Model modification
To apply the PnET-CN model to the DRB forests, we modified the model in three ways.
Parameterization. -The wood turnover rate (Woodtrn) is a parameter representing the rate of biomass loss from mortality. It controls what proportion of live wood biomass can be transferred to dead biomass on a yearly basis, and was set to 0.025 for all the vegetation types in the original model. The wood turnover rate can also be calculated using field-measured biomass data for each vegetation type. Tree mortality rate (M) between the two measurements was calculated as the biomass of trees that were live in the first measurement and died before the second measurement, divided by the total live tree biomass. The relationship between the mortality rate between measurements and the yearly wood turnover rate of trees (Treetrn) is expressed as: where n is the number of years between the two measurements, M is the measured tree biomass mortality rate between the two measurements, Treetrn is the yearly turnover rate of trees, and a is the proportion of biomass loss from branch and limbs according to tree mortality. The average n and M were calculated for each vegetation type separately, and Treetrn values were determined algebraically using Eq. 1. To include the biomass turnover of branch and limbs from living trees, the biomass ratio between small pieces of coarse woody debris (CWD, diameter <5 in.) and large pieces of CWD (diameter >5 in.) in the 2012-2014 measurements was used to estimate the proportion of biomass loss from branch and limbs according to tree mortality. In all the field measurement sites, the biomass of dead branch and limbs was equal to 7.25% of the dead tree biomass (a = 7.25%). Therefore, the Treetrn for each vegetation type was increased by 7.25% to calculate Woodtrn of the ecosystem (Eq. 2).
Calibration.-PnET model results are sensitive to the parameters AmaxA and AmaxB, which are the intercept and slope of the linear relationship between maximum photosynthesis rate and foliage N concentration. However, this relationship can vary with forest conditions such as species composition and nutrient availability (Schulze et al. 1994. In this study, a numerical approximation process was used to find the best AmaxA and AmaxB for each plot. The model ran multiple times with self-adjusted AmaxA and AmaxB on a step of 0.02. The best sets of parameters were reported when the simulated biomass in the last year (2010) matched with the observed value in the 2012-2014 measurement (within a range of AE200 g/m 2 ). The optimized parameters of each vegetation type were calculated as the mean of the best-fit AmaxA and AmaxB of all the plots. In total, 11 northern hardwood forest type plots, 39 red oakred maple forest type plots, and 3 pine forest type plots passed the approximation process. Their best-fit AmaxA and AmaxB values were used to calculate the optimized parameters. We had only one plot representing the spruce-fir forest type in the DRB sites, and as a result, the ❖ www.esajournals.org AmaxA and AmaxB values for this vegetation type were not calibrated.
Modification of input data.-In the original settings of the PnET-CN model, it is assumed that no disturbance occurred over the past 210 yr, and as a result, modeled forest ages are over 200 yr. However, forest plots in the DRB have very diverse land-use history and ages of the forests differ among plots. The PnET-CN model has a harvesting section in the scenario module, in which a constant proportion of biomass can be removed at the harvest year. To improve the model, we created a new input variable for stand age, and used it to calculate the year of the last disturbance. In that year, 80% of the biomass was disturbed and 50% of the disturbed biomass was removed from the biomass C pool ). Therefore, the modified PnET-CN model was able to reflect forest recovery processes after disturbances and better represent the different successional stages of each plot. Field-measured stand age for each plot was used for the model calibration and verification process. To run the GIS version of the model in the three sites, the stand age for each pixel was derived from a national stand age map based on FIA plot measurement (Pan et al. 2011).

Model verification
Live biomass data from the 2012-2014 measurement were compared to model output during the calibration step, while live biomass data from the first field measurement (2001)(2002)(2003) and the change in biomass between the two measurements were used to verify the modified model. Biomass was not used directly during calibration, but used as a benchmark for calibrating maximum photosynthetic rates in the model (AmaxA and AmaxB) until simulated the biomass matching with the observations at the end of the study period. Our assumption is that if these most important and sensitive parameters (which affect all the carbon dynamics processes) are correctly set up, we could expect that the modified model would be able to appropriately simulate annual photosynthesis, annual net primary production, annual biomass accumulation over simulation years, and finally total biomass in the study sites. As such, the past measurements of biomass and biomass change between measurements were used as independent data sets to verify the modified model with new maximum photosynthetic rates.
The average yearly biomass change of each plot was calculated as the difference of measured live biomass divided by the number of years between the 2001-2003 and 2012-2014 measurements. Since the time interval between the two measurements was not exactly ten years in most plots, and since we assume live biomass in the second measurement represents the observed biomass in 2010, the observed biomass in the first measurement was adjusted to the biomass in 2000 based on the average yearly change of biomass for each plot.
The normalized root mean square error (NRMSE) was used to evaluate overall model performance. The NRMSE assesses the difference between predicted (P) and observed (O) variables and can be expressed as: ( 3) where O max and O min are the maximum and minimum observed values, respectively; i is the ith observation plot; and n is the total number of plots. A NRMSE value close to 0 indicates perfect agreement, while a value of 100% indicates poor agreement. The NRMSE for the original PnET model and the modified model was compared, such that an observed decrease in NRMSE would indicate that the modification processes improved model performance.

Simulations
The original and modified models were used to simulate the C cycling in forested areas in each of the three sites. The input data at 1 9 1 km resolution were used to run the model for 210 yr (from 1801 to 2010). From the model output, four variables were mapped and compared between the original and modified PnET models. Vegetation biomass (VegM) and soil mass (SoilM) represent the major C pools, and NPP and NEP represent the major C fluxes in the forest C cycle. The NPP and NEP of the last ten years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) were averaged and reported.
To evaluate the effect of N deposition and elevated CO 2 , and their combined effects, four scenarios were created under which to run the modified PnET model ( Table 2). The full scenario represented the combined effects of elevated CO 2 ❖ www.esajournals.org and N deposition. In this scenario, atmospheric CO 2 concentration was ramped according to the Mauna Loa record. Atmospheric N inputs were calculated using historic data for local N deposition during 1982-2007 and ramped in other years. In the CO 2 removed scenario, the atmospheric CO 2 concentration was set to a fixed pre-industrial level (280 ppm). In the N removed scenario, the N input from N deposition was set to 0 g/m 2 . In the CO 2 +N removed scenario, both the effects of elevated CO 2 and N deposition were removed. Net primary productivity over the past 100 yr in each of the scenarios was plotted for each site (e.g., NS, DEWA, and FC) to illustrate the effects of elevated CO 2 and N deposition over time.

Model performance
The wood turnover rates calculated using field measurements were 0.0147 yr À1 in northern hardwood forest, 0.0208 yr À1 in red oak-red maple forest, and 0.0129 yr À1 in pine forest (Table 3). For all forest types, the measured wood turnover rate was smaller than the value used in the original PnET-CN model (0.025 yr À1 ). In the northern hardwood forest, the average best-fit AmaxA value was À48.6 lmol CO 2 Ám À2 ÁleafÁs À1 and AmaxB was 69.3 lmol CO 2 Ám À2 ÁleafÁs À1 , both close to the values used in the original model. In red oak-red maple forest, the average best-fit AmaxA value was À53.6 lmol CO 2 Ám À2 ÁleafÁs À1 and AmaxB was 64.3 lmol CO 2 Ám À2 ÁleafÁs À1 , smaller than the original values. In the pine forest, the average best-fit AmaxA and AmaxB values were 7.8 lmol CO 2 Ám À2 ÁleafÁs À1 and 24.0 lmol CO 2 Ám À2 ÁleafÁs À1 , respectively, both greater than the original values (Table 3).
Compared with the observed live biomass in the 61 plots during the two measurements (i.e., 2000 and 2010), the original model overestimated the mean biomass in FC and DEWA (above the 1:1 line in Fig. 2a), but underestimated the mean biomass in NS (below the 1:1 line in Fig. 2a). In   Fig. 2a). The original model was also unable to adequately simulate spatial variations of biomass on the landscapes and resulted in a narrow range of predicted biomass values, with the SD of spatial variations about 70% smaller than the observed SD in each site (Table 4). The modified model improved the prediction of mean biomass of the 61 plots and also increased the range of spatial variation, though the variation was still smaller than the observed value. In particular, the modified model accurately predicted the mean biomass in 2000, which was not affected by the calibration procedure but matched very well with the observed values for model verification (the blue dots falling on the 1:1 line in Fig. 2b). The model-predicted and observed biomass in 2010 was comparable in FC and DEWA. In NS, the model underestimated biomass, mostly because some plots in this area contained a small amount of conifer trees such as hemlocks (Tsuga canadensis), which tend to have larger biomass than typical northern hardwood forests. The biomass change between 2000 and 2010 predicted by the modified model was greater than the biomass change predicted by the original model in FC and NS, which was closer to the large biomass change observed in the field measurements. However, the biomass change in DEWA was smaller in the modified model because of the greater stand age in this site. Although the modified model significantly improved the predictions of mean biomass, it did not necessarily reduce prediction error on a plot-by-plot basis. The modified model slightly reduced NRMSE for FC and DEWA, but slightly increased NRMSE for NS (Table 4). The overall NRMSE of predictions for all plots in three sites by the modified model were similar to those of the original model.

C pools and fluxes in DRB
After the modification process, modelpredicted vegetation biomass over the entire areas of the three sites was greater in some pixels but smaller in other pixels in each of the three sites (Figs. 3a, b, 4a, b and 5a, b), with increased spatial variations due to the parameter modifications in mortality rate, maximum photosynthesis rate, and stand age. Vegetation biomass in FC and DEWA was slightly smaller in the modified model compared to the original model, but greater in NS (Table 5). Soil C was smaller in the modified model in all sites ( Table 5). The model modification process also resulted in a smaller NPP in all sites. Although NEP was either greater or smaller in different pixels compared to the original model (Figs. 3g, h, 4g, h and 5g, h), on average, the NEP predicted by the modified model was greater than the NEP predicted by the original model in all three sites. Compared to the original model, the modified model performed better in capturing the Notes: The mean, standard deviation (SD) of spatial variations, and normalized root mean square error (NRMSE) for estimating predicted errors by the original model and modified PnET-CN models are shown. The observed forest biomass data for 2000 and biomass change over decade were not used in the modification process and therefore could be considered as verification of model results. spatial complexity of the C pools and fluxes over the whole sites (Figs. 3-5). The SD of spatial variations among all the pixels was larger in the modified model, especially for NEP, which is the net C sink of the ecosystem ( Table 5).
Results of the modified model showed that, among the three sites, NS had the largest biomass C pool (168 Mg C/ha), NPP (41 Mg CÁha À1 Áyr À1 ), and NEP (0.79 Mg CÁha À1 Áyr À1 ), while DEWA had the largest soil C pool (42 Mg C/ha; Table 5). Within each site, large C pools were observed in the northern part of the FC (Fig. 3b, d), northerncentral part of the DEWA (Fig. 4b, d), and northern part of the NS (Fig. 5b, d). Large C sinks were observed in the southeast part of FC (Fig. 3h), center part of DEWA (Fig. 4h), and center and east part of NS (Fig. 5h).

Effects of environmental changes
In scenario simulations with the modified model, NPP values changed dramatically over time (Fig. 6). Net primary productivity decreased substantially in the year of major disturbances and increased gradually in the following years. Removing N deposition (N removed scenario) led to a decrease in NPP over last 100 yr of modeling. Setting the atmospheric CO 2 concentration to a fixed value (CO 2 removed scenario) resulted in even smaller modeled NPP in DEWA and NS. Removing both N deposition and elevated CO 2 (CO 2 +N removal scenario) led to the smallest NPP in all sites. The effect size between the CO 2 removed and CO 2 +N removed scenarios represented the effect of N deposition, while the effect size between the N removed and CO 2 +N removed scenarios represented the effect of elevated CO 2 . The effect of N deposition appeared to be greater than that of elevated CO 2 during 1970-1990in FC, and during 1970-1980 in DEWA and NS (Fig. 6b, c). However, the effect of elevated CO 2 was greater than the effect of N deposition after 1990 in DEWA and NS (Fig. 6b, c). The combined effect of N deposition and elevated CO 2 was also larger than sum of the two separate effects. The difference between the full scenario and other scenarios also increased over time, especially in the last 30 yr of the model simulations (Fig. 6). The forest biomass predicted by the modified model under the CO 2 +N removed scenario was 17-20% smaller than the biomass predicted in the full scenario, indicating that elevated CO 2 and N deposition increased forest biomass C stocks in all three sites (Table 5). The removal of both elevated CO 2 and N deposition decreased NPP in the last ten years by 17-21% and decreased NEP by 29-48% compared to the full scenario.

DISCUSSION
To improve the performance of a process-based model at the regional scale, model calibration or parameterization should target the parameters that the major output variables are sensitive to, and represent the spatial variation of forest C dynamics (Thorn et al. 2015). In this study, we modified the wood turnover rate, the maximum Fig. 4. Spatial distributions of forest biomass (a, b), soil C (c, d), net primary productivity (NPP) (e, f), and net ecosystem productivity (NEP) (g, h) as simulated by the original (a, c, e, g) and modified (b, d, f, h) PnET models in the Delaware Water Gap study site in the Delaware River Basin. photosynthesis rate, and the disturbance years in the PnET-CN model. The main objective of these modifications was not only to improve the accuracy of the model prediction compared to field observation, but also to increase the ability of the model to capture the spatial complexity of forest C cycling at the landscape scale.
Importance of wood turnover rate, maximum photosynthesis rate, and stand age in C cycle modeling Mortality is a key driver of the spatial distribution of forest biomass and also affects soil C pools and sinks. A better estimation of mortality rate can strongly improve modeling of biomass Fig. 5. Spatial distributions of forest biomass (a, b), soil C (c, d), net primary productivity (NPP) (e, f), and net ecosystem productivity (NEP) (g, h) as simulated by the original (a, c, e, g) and modified (b, d, f, h) PnET model in the Neversink study site of the Delaware River Basin. and its spatial variation (Delbart et al. 2010). Comparisons among multiple process-based models often result in divergent predictions in C sinks, which were caused by different formulations for turnover rates (Warszawski et al. 2013). Compared with C flux data, long-term biomass change data were more effective in constraining C turnover rates (Du et al. 2015). In this study, the wood turnover rates calculated from the observed mortality rate were directly used in the modified model. The smaller observed wood turnover rates in all three vegetation types are very likely due to the forests in the DRB generally being in their mid-ages, and subject to smaller mortality rates (Coomes andAllen 2007, Xu et al. 2016b). Our results showed that smaller wood turnover rates resulted in larger C accumulation in the live biomass C pool, but reduced the C stocks in soil. Increasing autotrophic respiration from larger biomass leads to a decreasing NPP. However, heterotrophic respiration is likely to be limited by the smaller C input to the soil C pool. As the balance between NPP and heterotrophic respiration, NEP did not have a consistent trend in its response to a smaller wood turnover rates. The effect of modified wood turnover can be detected in northern part of NS, where Amax did not change substantially and the forests are relatively older (Fig. 5).
Lower optimized AmaxA and AmaxB values for the red oak-red maple forest type, but not for the northern hardwood forest type, indicated that the original model parameters (which were developed in the northern hardwood forests in the New England area) could be used in the northern hardwood forests in the DRB, but may substantially overestimate forest biomass and productivity in other forest types. To adapt their climate conditions, northern hardwood trees usually exhibit greater photosynthetic rates in the relatively shorter growing season (Pan et al. 2006). When these photosynthetic rates were used in the southern oak-dominated deciduous forests, it caused overestimates in productivity and biomass. Because AmaxA and AmaxB are the parameters of a linear function of foliage N concentration and photosynthesis rate, consistently lower photosynthetic rates related to foliage N concentration in southern trees (Prentice et al. 2014) could reflect relatively lower N use efficiency, while northern hardwoods are likely to have greater N use efficiency under more N limited conditions (Aber and Melillo 2002). In the red oak-red maple forest, the smaller AmaxA and AmaxB values used in the modified model led to smaller living biomass and soil C pools and smaller NPP and NEP.
Successional stage is a critical factor determining forest C dynamics. Many process models have been used to simulate the effect of forest harvesting and forest regrowth processes (Chiang et al. 2008, Potter et al. 2008  Notes: Values are reported as mean AE standard deviation. The effects of elevated CO 2 and N deposition were calculated as the percentage change of each variable from the full scenario to the CO 2 +N removed scenario using the modified model. Negative effect values indicate that elevated CO 2 and N deposition resulted in increased values of the variables in the full scenario. 2014a). In the PnET model, after a destructive disturbance and losses of C stocks, forest biomass and soil C pools increase during recovery, and NPP and NEP increase in earlier development years or decades and then decrease after the canopy is closed and with increasing respiration (Wang et al. 2014b). As a result, younger forests are expected to have smaller C pools and larger productivity compared to older forests approaching an equilibrium stage. This trend has been tested and verified by long-term field measurements and chronosequence studies (Law et al. 2003. Our results show that incorporating stand age data (Pan et al. 2011) into the PnET model can dramatically change the simulated distribution of C pool and fluxes in the DRB forests . For example, the center area of DEWA experienced a clear cut in the late 1960s, but the impact of this disturbance was not detected by the original model. In the modified PnET model results, this area had relatively smaller biomass and soil C stocks and became a hotspot of C sequestration with the largest NEP.

Performance of the improved model
The goal of the model modifications was twofold: to account for the real and substantial spatial variability in forest C stocks and fluxes, which is expressed by the SD around the mean model output for each site, and to improve the accuracy of the model output, which is expressed in the differences in modeled vs. measured data and the NRMSE. In the end, the hoped-for goal was simultaneous increases in SD and decreases in NRMSE.
The modification process generally improved the performance of the PnET-CN model in the DRB forests. The observed and predicted mean biomass in each site matched well. Our modifications also largely increase the SD of all the model outputs (e.g., error bars in Fig. 2, and SD in Tables 4 and 5), which suggests the model can better capture the spatial variation in these small regions. However, the modified model still predicted SD much smaller than the observed SD. The highly scattered data points in Fig. 2b highlight the difficulty of simulating the heterogeneity of forest C dynamics. One of the reasons is that the calibrated photosynthetic rates had to be averaged by forest types while the calibrations were performed at the individual plot scale. As a result, some spatial variation was lost, which is an inevitable limit in the application of a spatially explicit process model and comparing the landscape-scale results with plot-level measurements.
Although the modified model significantly improved the mean values of predicted biomass, prediction errors based on NRMSE did not always decrease, and in the case of NS, actually increased. As mentioned earlier, spatial heterogeneity in biomass observed in plots was the result of complex biotic and abiotic factors. Even after modification of a few key parameters, the model still cannot fully account for the effects of diverse factors when the model must use averaged parameters and driving climate/soil data derived from coarse resolution data sets (Lichstein et al. 2014). As an example, the low-elevation areas in NS commonly contain patches of hemlock trees, which tend to have larger growth rates and hence can accumulate more biomass C stock at young stand ages. However, because hemlock is not identified as a forest type in the 1-km-resolution forest type map, the model-predicted biomass in the low-elevation NS areas was limited by the observed young stand age. Although the modified model could be improved by better observations of stand age, this would cause even greater NRMSE from the observed high biomass in hemlock forests. This case also indicated the potential for further model improvement by calibrating additional forest types and better accounting for the effect of species composition (Porte and Bartelink 2002).
Soil C pools estimated by the original model were smaller than the field-measured soil C stocks from the same research sites (Xu et al. 2016a). Reduced wood turnover rate in the modified model resulted in even smaller soil C pools. A possible reason for this model-data inconsistency is that the soil profile depth in the PnET-CN model is not clearly defined. Soil C pools estimated by the model are therefore not easily comparable with field measurements that are sampled by depth increments. The ability of the PnET model to represent the soil C dynamics is also limited by the simplified belowground processes and the lack of spatialized input information of soil characteristics.

Effects of N deposition and elevated CO 2 in the DRB forest
In field measurements, the observed biomass C change is a combined effect of natural, baseline forest growth and growth enhancement due to environmental changes such as elevated CO 2 and N deposition (Caspersen et al. 2000). Process-based models provide a useful tool to quantify the contribution of these environmental changes separately, which cannot be identified by observations alone , Pan et al. 2009).
Results of model scenarios (Fig. 6) demonstrated that elevated CO 2 and N deposition had a large impact on NPP in recent decades. The effect of N deposition was larger than that of elevated CO 2 before 1980. But after 1980, the effect of elevated CO 2 was greater in two of the three sites. The effect of N deposition estimated in this study was smaller than the effect estimated by Pan et al. (2009). Possible reasons include the decreasing N deposition in the historical input data after 2000, and the smaller AmaxA and AmaxB values used in the modified model. In addition, young forests appeared to have smaller response to N deposition than mature forests as we improved representation of land-disturbance consequence with spatial forest age data in this study. It is also implied that the effect of elevated CO 2 will likely increase continuously in the future even in forests, which was considered to be limited by N supply.
The combined effect of the two factors was greater than the sum of the two separated effects, which demonstrated a positive interaction between N deposition and elevated CO 2 and is consistent with the results from controlled experiments and with other model simulations (Nowak et al. 2004, Richardson et al. 2010. By removing these factors from the modified PnET-CN model, we found that 19% of the current biomass C stock and 39% of the NEP were attributable to the effects of elevated CO 2 and N deposition, suggesting that environmental change had substantially changed the C dynamics in the DRB forest in recent decades. Because of the effects, these forests remained vigorous with high productivity when they would otherwise have lower productivity as they grow older. Using similar scenarios designs, the effect of other factors such as ozone damage, climate change, and landscape change induced by disturbance and their interactions can also be addressed using the PnET-CN model. The modified model would be able to provide more precise evaluation of these effects, and their spatial distribution within the study sites. Implications of using carbon stock field measurements for improving a spatially explicit model Model simulations extended our understanding of the DRB forest C cycle spatially and temporally (from 61 plots to the area of three sites, and from two measurements ten years apart to a 210-yr time sequence) and generated additional information about C stocks and fluxes (from measured live biomass to modeled biomass and soil C pools, as well as NPP and NEP). In this study, forest biomass C stock measurements were used as the key benchmark in the model verification processes. Compared to C fluxes variables, which are frequently used in model calibrations based on eddy covariance data, forest biomass better represents the spatial variation of forest C cycling and synthesizes more information about the long-term processes of forest dynamics since it is an accumulative variable (Miehle et al. 2006, Kondo et al. 2013. Moreover, field measurement of biomass is feasible in all the forested area, whereas eddy covariance towers are far scarcer and cannot be installed in highly complex topographies. The FIA database has collected periodic measurements of forest biomass all over the country, and the structure of the database is similar to our observational data from the three sites in the DRB forest (Potter et al. 2008). Biomass data as perhaps the most available and robust field measurements should be considered more frequently as the benchmark data for constraining and improving process-based models. The methods of model modification used in this study could therefore be easily used at a larger scale and in different locales to improve model predictions of forest C and N cycling.
Because spatially explicit process-based models have become important tools for exploring heterogeneous carbon dynamics on landscapes, the PnET models have been developed or coupled with other models for spatial applications (Pan et al. 2006, De Bruijn et al. 2014, Gustafson et al. 2015. However, the reliability of those models is often limited by the spatial extent of the observation data that were used to parameterize the models. Our study demonstrated a feasible method for using forest biomass measurements from multiple plots to improve representations of a process-based model on landscapes, which can provide more relevant and useful information for local forest management, such as the distribution of carbon sequestration hotspots and the impacts of N deposition and elevated CO 2 on specific forests in the study area. Our study shows that the relationship between maximum photosynthesis rate and foliage N concentration is one of the key principles in the PnET model that links the C and N cycles and determines forest productivity on a monthly basis (Aber and Federer 1992). However, few studies have measured this relationship in different forest types, and the spatial variations of the parameters are not well constrained. In fact, the lack of information on these parameters is a major source of uncertainty in using the PnET model in different vegetation types (Thorn et al. 2015). Although remote sensing involving with hyperspectral sensors and expensive field probes tried to fill this data gap (Martin et al. 2008), the numerical approach used in our model calibration provided an effective method to quantify this relationship and to test the variability of these parameters across different forest compositions, ages, and environmental conditions.
To further improve the PnET-CN model, which is particularly important and necessary for regional or local applications, input data of higher spatial resolution should be used and more long-term field measurement data are needed (Campbell et al. 2009). Foliage N concentration is an important and sensitive variable that determines maximum photosynthetic rates in the PnET model (Pan et al. 2004b). If field measurements of foliage N concentrations are available, they should be used to parameterize the model for each forest type. Soil WHC is one of the most important input data and is the only soil property that is required in the PnET-CN model. The NRCS Web Soil Survey, a new soil database with high spatial resolution, has recently been published (Soil Survey Staff 2014) and could be used as a better source of soil WHC data to improve the model. Large areas of the northeastern United States are covered by mixed forests, and species composition can change over time. More accurate model predictions will result if the model applied to this region uses input data that retain fine-resolution forest type variability, and with more forest types parameterized in the model (Jenkins et al. 1999). Incorporating forest species dynamics in the model may also greatly improve the ability of the PnET model to capture spatial variability. Improving model simulations to generate more robust output and spatial variations that reflect actual spatial heterogeneity is necessary to extrapolating results of the field survey to the entire DRB forest and to test the effects of additional environmental changes (e.g., ozone) or disturbances (e.g., invasive pests). Ultimately, robust model output is essential to improving our understanding forest carbon dynamics and for improved forest carbon management.