A forest biogeochemistry model intercomparison on the East Coast of the United States

Recent advances in forest ecosystem modeling allow the simulation of a suite of dynamics from site- to landscape-scale. In order to scale models efficiently from trees to landscapes, different model reduction strategies are employed. Yet, the results of these strategies and the assumptions they entail are rarely compared. Here, we conducted a model intercomparison exercise using two such forest biogeochemistry models, PPA-SiBGC and LANDIS-II NECN. We simulated past-decade conditions at flux tower sites in Harvard Forest, MA, USA and Jones Ecological Research Center, GA, USA. We mined the wealth of field data available for both sites to perform model parameterization, validation, and intercomparison. We assessed model performance using the following time-series metrics: net ecosystem exchange, aboveground net primary production, aboveground biomass, C, and N, belowground biomass, C, and N, soil respiration, and, species total biomass and relative abundance. We also assessed static observations of soil organic C and N, and concluded with an assessment of general model usability, performance, and transferability. Despite substantial differences in design, both models achieved good accuracy across the range of metrics. While LANDIS-II NECN performed better for interannual NEE fluxes due to its basis in the Century model, the PPA-SiBGC model indicated better overall correspondence to observational data for both sites across the 11 temporal and 2 static metrics tested .


Introduction
Century model [19,20]. The original ten soil layers in Century have been replaced by a single soil layer, 180 with functions for growth and decay borrowed directly from Century v4.5. The NECN succession 181 model Figure 1 is thus a process-based model that simulates C and N dynamics along the plant-soil 182 continuum at a native monthly timestep. gap model [63,64]. The Perfect Plasticity Approximation, or PPA [33,34], was derived from the dual 196 assumptions of perfect crown plasticity (e.g., space-filling) and phototropism (e.g., stem-leaning), both 197 of which were supported in empirical and modeling studies [36]. The discovery of the PPA was rooted 198 in extensive observational and in silico research [33]. The PPA model was designed to overcome the most computationally challenging aspects of gap models in order to facilitate model scaling from the 200 landscape to global scale.

201
The PPA and its predecessor, the size-and-age structured (SAS) equations [30,65], are popular 202 model reduction techniques employed in current state-of-the-art terrestrial biosphere models [13]. The

203
PPA model can be thought of metaphorically as Navier-Stokes equations of forest dynamics, capable of 204 modeling individual tree population dynamics with a one-dimensional von Foerster partial differential 205 equation [33]. The simple mathematical foundation of the PPA model is provided in Equation 1.
where k is the number of species, j is the species index, N j (z) is the density of species j at height z, 207 A j (a ⇤ , z) is the projected crown area of species j at height z, and dz is the derivative of height. In other 208 words, we discard the spatial location of individual trees and calculate the height at which the integral 209 of tree crown area is equal to the ground area of the stand. This height is known as the theoretical z ⇤ 210 height, which segments trees into overstory and understory classes [33].

211
The segmentation of the forest canopy into understory and overstory layers allows for separate 212 coefficients or functions for growth, mortality, and fecundity to be applied across strata, whose first Empirical observations were relied upon for the C and N content of tree species compartments.

219
Stoichiometric relations were used to estimate N from C, based on empirical measurements provided 220 for both sites. All values were calculated directly from observations. Previously published equations 221 [71] and parameters [72] were used to model crown allometry. Together with inventory data, general 222 biomass equations were used to estimated dry weight mass (kg) for tree stems, branches, leaves, and, 223 fine and coarse roots [73]. Carbon content is assumed to be 50% of dry mass, supported by data.

224
Monthly soil respiration is modeled using the approach of Raich et al.   A table of observed species abundances for the 2009-2013 period are provided in Table 4. The selection of simulation years was based on the availability of EC flux tower data used in

298
The areal extent of the single-site model intercomparisons were designed to correspond to 299 available field measurements. At both sites, tree inventories were conducted in 10,000 m 2 , or 300 one-hectare, areas. All target metrics were converted to an annual areal basis to ease interpretation, 301 comparison, and transferability of results. Importantly, an areal conversion will allow comparison to 302 other sites around the world. While flux tower measurements for both sites were already provided 303 on an areal (m 2 ) basis, many other variables were converted to harmonize metrics between models 304 and study sites. For example, moles CO 2 measurements were converted to moles C through  Table 5.  .
where n is the sample size, y i is the ith observed value,ŷ i is the ith predicted value, y is the mean 324 observed value, andŷ is the mean predicted value. The calculation of RMSE follows the standard 325 formulation, as shown in Equation 3.
where n is the sample size and e t is the error for the tth value, or the difference between observed where again n is the sample size and e t is the error for the tth value. Our calculation of mean 329 error (ME) or bias is the same as MAE, but without taking the absolute value. quantify the direction (i.e., bias) and/or magnitude of error. We rely on RMSE and MAE to provide 340 information on error or residual magnitude, and ME to provide information on bias. We utilize a 341 visual analysis to assess error directionality over time, as this can be poorly characterized by a single 342 coefficient, masking periodicity. 343 We compute R 2 , RMSE, MAE, and ME for time-series of the metrics described in Table 5    where this is the case, while the correlation of this metric to observed values was also lower than that 388 of PPA-SiBGC. Overall results for the HF-EMS site model intercomparison are shown in Table 7. Time-series figures allow a visual analysis of the temporal dynamics between observations and 390 model predictions in order to assess the ability of models to capture interannual variability. Both 391 models effectively captured temporal dynamics in biomass, C, and, species biomass and abundance.

392
In Figure 4, the temporal differences in modeled NEE and aboveground C are shown for the two 393 models in comparison to observations for the HF-EMS site. While LANDIS-II NECN predicted NEE 394 showed a higher correlation with observations, the magnitude of error and bias were also higher.

395
Furthermore, LANDIS-II NECN predicted that the HF-EMS site is a net C source, rather than sink, in 396 contrary to observations. Meanwhile, PPA-SiBGC outperformed LANDIS-II NECN in aboveground C 397 per both R 2 and RMSE. Both models overpredicted species cohort biomass, while LANDIS-II NECN 398 underpredicted total aboveground C.
[a] For the JERC-RD site, both models showed stronger fidelity to data than for the HF-EMS site. While both models showed higher performance at the JERC-RD site, an analysis of simulated 417 species biomass and abundance again indicates greater fidelity of the PPA-SiBGC model to data, as shown in Figure 6. While LANDIS-II NECN overpredicts the rate of longleaf pine growth, PPA-SiBGC 419 nearly perfectly matches observed species abundance and biomass trajectories for all species present.

420
While the correlations are high, PPA-SiBGC overpredicts the magnitude of biomass here.

421
[a]         for this validation and intercomparison exercise, fully mechanistic processes were not required; both 513 models required a host of empirical parameters that limit their prognostic abilities in their current form.

514
Future developments with both models should improve upon this by adopting more mechanistic 515 processes.

516
Replacing the simple biogeochemistry approach of PPA-SiBGC with mechanistic processes would      The following abbreviations are used in this manuscript:

ANPP
Aboveground net primary production API      Figure A5. HF-EMS flux tower and landcover classes Next is the JERC-RD flux tower with landcover classes Figure A6.  Figure A6. JERC-RD flux tower and landcover classes