A forest model intercomparison framework and application at two temperate forests along the East Coast of the United States

: Forest models often reﬂect the dominant management paradigm of their time. Until the 1 late 1970s, this meant sustaining yields. Following landmark work in forest ecology, physiology, and 2 biogeochemistry, the current generation of models is further intended to inform ecological and climatic 3 forest management in alignment with national biodiversity and climate mitigation targets. This has 4 greatly increased the complexity of models used to inform management, making them difﬁcult to 5 diagnose and understand. State-of-the-art forest models are often complex, analytically intractable, 6 and computationally-expensive, due to the explicit representation of detailed biogeochemical and 7 ecological processes. Different models often produce distinct results while predictions from the same 8 model vary with parameter values. In this project, we developed a rigorous quantitative approach for 9 conducting model intercomparisons and assessing model performance. We have applied our original 10 methodology to compare two forest biogeochemistry models, the Perfect Plasticity Approximation 11 with Simple Biogeochemistry (PPA-SiBGC) and Landscape Disturbance and Succession with Net 12 Ecosystem Carbon and Nitrogen (LANDIS-II NECN). We simulated past-decade conditions at ﬂux 13 tower sites located within Harvard Forest, MA, USA (HF-EMS) and Jones Ecological Research Center, 14 GA, USA (JERC-RD). We mined ﬁeld data available for both sites to perform model parameterization, 15 validation, and intercomparison. We assessed model performance using the following time-series 16 metrics: net ecosystem exchange, aboveground net primary production, aboveground biomass, C, 17 and N, belowground biomass, C, and N, soil respiration, and, species total biomass and relative 18 abundance. We also assessed static observations of soil organic C and N, and concluded with 19 an assessment of general model usability, performance, and transferability. Despite substantial 20 differences in design, both models achieved good accuracy across the range of pool metrics. While 21 LANDIS-II NECN showed better ﬁdelity to interannual NEE ﬂuxes, PPA-SiBGC indicated better 22 overall performance for both sites across the 11 temporal and 2 static metrics tested (HF-EMS R 2 = 23 0.73, + 0.07, RMSE = 4.84, − 10.02; JERC-RD R 2 = 0.76, + 0.04, RMSE = 2.69, − 1.86). To facilitate 24 further testing of forest models at the two sites, we provide pre-processed datasets and original 25 software written in the R language of statistical computing. In addition to model intercomparisons, 26 our approach may be employed to test modiﬁcations to forest models and their sensitivity to different 27 parameterizations. 28 Time-series ﬁgures allow a visual analysis of the temporal dynamics between observations and model predictions in order to assess the ability of models to capture interannual variability. Both models effectively captured temporal dynamics in biomass, C, and, species biomass and abundance.


Introduction
where k is the number of species, j is the species index, N j (z) is the density of species j at height z, 243 A j (a * , z) is the projected crown area of species j at height z, and dz is the derivative of height. In other 244 words, we discard the spatial location of individual trees and calculate the height at which the integral 245 of tree crown area is equal to the ground area of the stand. This height is known as the theoretical z * 246 height, which segments trees into overstory and understory classes [48].

247
The segmentation of the forest canopy into understory and overstory layers allows for separate 248 coefficients or functions for growth, mortality, and fecundity to be applied across strata, whose first 249 moment accurately approximates the dynamics of individual-based forest models. Recent studies have 250 shown that the PPA model faithfully reduces the dynamics of the more recent neighborhood dynamics 251 (ND) SORTIE-ND gap model [81] and is capable of accurately capturing forest dynamics [82,83]. Empirical observations were relied upon for the C and N content of tree species compartments.

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

260
Monthly soil respiration is modeled using the approach of Raich et al. [84], while soil organic C is 261 All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; modeled using the simple generalized approach of Domke et al. [85]. Species-and stratum-specific 262 parameters for growth, mortality, and fecundity were calculated from observational data available for 263 both sites. In the following sections, we describe the two forested sites on the East Coast of the United States: 266 HF-EMS and the JERC-RD. A critical factor in the selection of the sites was the availability of eddy 267 covariance flux tower data needed to validate NEE in the models.  Table 1. 287 Table 1. HF-EMS forest inventory summary for the 34 tower plots in 2002; DBH = depth at breast height (cm); BA = basal area per hectare (m 2 ); Stocking = n trees per hectare; QMD = quadratic mean diameter (cm); SDI = Reineke's stand density index [90]  All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; period. Summary statistics for the RD tower site for the year 2008 are provided in Table 3.

315
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; The areal extent of the single-site model intercomparisons were designed to correspond to 336 available field measurements. At both sites, tree inventories were conducted in 10,000 m 2 , or 337 one-hectare, areas. All target metrics were converted to an annual areal basis to ease interpretation, 338 comparison, and transferability of results. Importantly, an areal conversion will allow comparison to 339 other sites around the world. While flux tower measurements for both sites were already provided 340 on an areal (m −2 ) basis, many other variables were converted to harmonize metrics between models 341 All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; and study sites. For example, moles CO 2 measurements were converted to moles C through 342 well-described molecular weights, all other measures of mass were converted to kg, and all areal and 343 flux measurements were harmonized to m −2 . A table of metrics and units used in the intercomparison 344 of LANDIS-II and PPA-SiBGC is provided in Table 5.

345
Aboveground net primary production kg mass m −2 year −1 B Sp Species aboveground biomass kg mass m −2 n Sp Species relative abundance % In the subsequent section, we describe the model intercomparison methodology. .
where n is the sample size, y i is the ith observed value,ŷ i is the ith predicted value, y is the mean 361 observed value, andŷ is the mean predicted value. The calculation of RMSE follows the standard 362 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 364 and predicted values. The calculation of MAE is similarly unexceptional, per Equation 4.

365
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; where again n is the sample size and e t is the error for the tth value. Our calculation of mean 366 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 377 information on error or residual magnitude, and ME to provide information on bias. We utilize a 378 visual analysis to assess error directionality over time, as this can be poorly characterized by a single 379 coefficient, masking periodicity.

380
We compute R 2 , RMSE, MAE, and ME for time-series of the metrics described in Table 5 on 381 page 11. These include NEE, above-and below-ground biomass, C, and N, soil organic C and N, soil 382 respiration (r soil ), aboveground net primary production (ANPP), and, species aboveground biomass 383 and relative abundance. All of these metrics are pools with the exception of NEE, r soil , and ANPP  (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018;  where this is the case, while the correlation of this metric to observed values was also lower than that 425 of PPA-SiBGC. Overall results for the HF-EMS site model intercomparison are shown in Table 7.

426
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; Time-series figures allow a visual analysis of the temporal dynamics between observations and 427 model predictions in order to assess the ability of models to capture interannual variability. Both 428 models effectively captured temporal dynamics in biomass, C, and, species biomass and abundance.

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

432
Furthermore, LANDIS-II NECN predicted that the HF-EMS site is a net C source, rather than sink, in 433 contrary to observations. Meanwhile, PPA-SiBGC outperformed LANDIS-II NECN in aboveground C 434 per both R 2 and RMSE. Both models overpredicted species cohort biomass, while LANDIS-II NECN 435 underpredicted total aboveground C.

436
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; 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 454 species biomass and abundance again indicates greater fidelity of the PPA-SiBGC model to data, as 455 All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; shown in Figure 6. While LANDIS-II NECN overpredicts the rate of longleaf pine growth, PPA-SiBGC 456 nearly perfectly matches observed species abundance and biomass trajectories for all species present.

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

458
[a]   (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.     Table 9, PPA-SiBGC was between 1,200 and 2,800% faster than LANDIS-II NECN in our timing tests.

505
This was surprising given that PPA-SiBGC models true cohorts (i.e., individual trees) in an interpreted 506 language while LANDIS-II models theoretical cohorts (i.e., cohorts without a physical basis) in a 507 compiled language. The difference in speed is likely attributable to the parsimony of the PPA-SiBGC 508 model. All rights reserved. No reuse allowed without permission.

509
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; the growth, mortality, and fecundity parameters used by PPA-SiBGC are easy to calculate using 522 common field inventory data. PPA-SiBGC is simpler to transfer in this regard given the wide 523 availability of forest inventory data.   (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; programming languages (e.g., [101]). Thus, SORTIE-NG is intended to be the first forest model to 571 bridge the divide between big-leaf, gap, and landscape models, and to be designed from the outset as a 572 probabilistic modeling framework [100]. Future model extensions are in the planning stages, including 573 the first machine-learning modeling interface included in an ecosystem model.

574
While implemented in a 'close-to-metal' language (i.e., C++17) and designed for efficiency, 575 SORTIE-NG is vastly more computationally demanding than the PPA-SiBGC model used in this 576 paper. Yet, we anticipate that SORTIE-NG will be able to improved the fidelity to observed fluxes, 577 which is the major shortcoming of both models considered in this paper. Similarly, there is a new      (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

6.
Vehkamäki, S., The concept of sustainability in modern times. In Sustainable use of renewable natural resources;   All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
1006 All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; Patterns of daytime and nighttime NEE are shown in Figure A4. Again, this was calculated by 1007 taking daily mean NEE values for three-hour windows surrounding noon and midnight, respectively 1008 (1100-1300 and 2300-0100 hours).

1009
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
1012 All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.  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 All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018;  (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/464578 doi: bioRxiv preprint first posted online Nov. 7, 2018; (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.