## 2 Materials and Methods
### Model Intercomparison  
Intercomparison of the PPA-SiBGC and LANDIS-II models at the HF-EMS and JERC-RD EC flux tower sites was conducted using a collection of object-oriented functional programming scripts written in the R language for statistical computing [83]. These scripts were designed to simplify model configuration, parameterization, operation, calibration/validation, plotting, and error calculation. The scripts and our parameters are available on GitHub ([[https://github.com/adam-erickson/ecosystem-model-comparison](https://github.com/adam-erickson/ecosystem-model-comparison)]([https://github.com/adam-erickson/ecosystem-model-comparison](https://github.com/adam-erickson/ecosystem-model-comparison))), making our results fully and efficiently reproducible. The directory structure of the repository is shown in Figure S1 in the Supplementary Materials. The R scripts are also designed to automatically load and parse the results from previous model simulations, in order to avoid reproducibility issues stemming from model stochasticity. We use standard regression metrics applied to the time-series of observation and simulation data to assess model fitness. The metrics used include the coefficient of determination (\(R^{2}\)), root mean squared error (RMSE), mean absolute error (MAE), and mean error (ME) or bias, calculated using simulated and observed values. Our implementation of \(R^{2}\) follows the Bravais-Pearson interpretation as the squared correlation coefficient between observed and predicted values [130]. This implementation is provided in Equation (2). \[R^{2}=r^{2}=\left(\frac{\sum_{i=1}^{n}(y_{i}-\overline{y})(\hat{y}_{i}- \overline{y})}{\sqrt{\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}(\hat{y}_{i}- \overline{y})^{2}}}\right)^{2} \tag{2}\] \begin{table}
\begin{tabular}{c c c} \hline \hline
**Abbreviation** & **Metric** & **Units** \\ \hline NEE & Net ecosystem exchange & kg C m\({}^{-2}\) year\({}^{-1}\) \\ \(B_{AG}\) & Aboveground biomass & kg mass m\({}^{-2}\) \\ \(C_{AG}\) & Aboveground C & kg C m\({}^{-2}\) \\ \(N_{AG}\) & Aboveground N & kg N m\({}^{-2}\) \\ \(B_{BG}\) & Belowground biomass & kg mass m\({}^{-2}\) \\ \(C_{BG}\) & Belowground C & kg C m\({}^{-2}\) \\ \(N_{BG}\) & Belowground N & kg N m\({}^{-2}\) \\ \(C_{SO}\) & Soil organic C & kg C m\({}^{-2}\) \\ \(N_{SO}\) & Soil organic N & kg N m\({}^{-2}\) \\ \(r_{soil}\) & Soil respiration C & kg C m\({}^{-2}\) year\({}^{-1}\) \\ ANPP & Aboveground net primary production & kg mass m\({}^{-2}\) year\({}^{-1}\) \\ \(B_{Sp}\) & Species aboveground biomass & kg mass m\({}^{-2}\) \\ \(n_{Sp}\) & Species relative abundance & \% \\ \hline \hline \end{tabular}
\end{table} 
Table 5: Model intercomparison abbreviations, metrics, and units. where \(n\) is the sample size, \(y_{i}\) is the \(i\)th observed value, \(\hat{y_{i}}\) is the \(i\)th predicted value, \(\overline{y}\) is the mean observed value, and \(\overline{y}\) is the mean predicted value. The calculation of RMSE follows the standard formulation, as shown in Equation (3). \[\text{RMSE}=\sqrt{\frac{1}{n}\sum_{t=1}^{n}e_{t}^{2}} \tag{3}\]  
where \(n\) is the sample size and \(e_{t}\) is the error for the \(t\)th value, or the difference between observed and predicted values. The calculation of MAE is similarly unexceptional, per Equation (4). \[\text{MAE}=\frac{1}{n}\sum_{t=1}^{n}|e_{t}| \tag{4}\]  
where again \(n\) is the sample size and \(e_{t}\) is the error for the \(t\)th value. Our calculation of mean error (ME) or bias is the same as MAE, but without taking the absolute value. While Nash-Sutcliffe efficiency (NSE) is often used in a simulation model context, we selected the Bravais-Pearson interpretation of \(R^{2}\) over NSE to simplify the interpretation of results. The NSE metric replaces \(1-(SS_{predictions}/SS_{observations})\) with \((SS_{observations}-SS_{predictions})/SS_{observations}\), where \(SS\) is the sum of squares. Thus, NSE is analogous to the standard \(R^{2}\) coefficient of determination used in regression analysis [131].