Comparing models can be seen as an alternative way to validate the credibility of model outcomes and results. These two models belong to empirical forest-inventory–based bookkeeping models. We showed that when maximum harmonization is strived for in the input and scenarios, these models produce comparable results and, to a large extent, the same trends. However, when the models employ different approaches (such as for the soil compartment), the output dynamics differ. Our results highlighted those small differences in modelling principles of CBM and EFISCEN combined with the available data can yield differences in output estimates on various components. Despite the efforts invested in harmonization of the input parameters, differences remained between the estimates of the two models. Here, we discuss in detail the potential causes of these differences.
Uncertainty of input data
Developing an uncertainty or sensitivity analysis for CBM and EFISCEN was beyond the aim of our study. However, it is obvious that the estimates of both models are affected by uncertainty. According to the NFI [47, 51], the sampling error estimated for each national-scale valid indicator was typically low (e.g., 6% for fellings and 1.7% for the merchantable standing stock). Since the input data were disaggregated, it is expected that the standard error will increase (e.g., up to 12% for highly represented, e.g., Fagus sylvatica, Norway spruce and other broadleaved forests, and up to 100% for low represented, e.g., other coniferous and silver fir forests).
For CBM, an important source of uncertainty originates in the fit of yield, increment and Boudewyn models. These models were fitted on nationally available data. For the yield models, the residual standard error of the yield models, also relative to the mean predicted value, varied by forest type between 49 m3 ha−1, or− 20.8%, for F. sylvatica and 172 m3 ha−1, or 94.8%, for Abies alba forest type. For NAI models, the residual standard error varied between 0.68 m3 ha−1 yr−1 for Robinia pseudoacacia and 3.2 m3 ha−1 yr−1 for the “Other Broadleaved” forest type. The Boudewyn models, used to predict the proportion of biomass compartments, showed residual standard errors between 1.1% for the bark of “PredBroad” forest type and 11.5% for the proportion of R. pseudoacacia stemwood. EFISCEN uses BEF, developed for each forest type and each biomass compartment, by age class. The standard deviation of BEF varied between 6.3% for stem biomass of A. alba trees and 0.2% for branch biomass of the “Other Broadleaved” forest type. These standard deviations include, however, the variation by age class. Therefore, the standard deviation of allocation in biomass compartments within age class is expected to be smaller than the range presented.
To reduce the uncertainty sourced in these models (which act as input data in either EFISCEN or CBM), additional observations at tree and stand level should be collected from Romanian forests based on a robust sampling design. Nevertheless, data analysis should engage appropriate model selection criteria to minimize the uncertainty related to model selection.
Merchantable standing stock
Compared with the NFI estimate, in 2010, the simulated merchantable stock volume was 7% larger in CBM and only 1.5% larger in EFISCEN (see Fig. 2). However, for C stock, the difference between models and the NFI was 9.9%. Furthermore, between the models, the differences in merchantable standing stock volume increased in time (from 6.8% to 8.1%), whereas the differences in C stock reduced (from 13.7% to 11.7%), mainly because of the change in the share among forest types with different wood densities. These differences between merchantable standing stock volume and carbon stock may generate confusion when making climate decisions based on forestry indicators. For example, the forest can be assumed sustainable from a forestry perspective (when the harvest volume is slightly lower than the increment volume); however, it may be assumed not sustainable from a corresponding climatic perspective (C from harvest is larger than C from forest growth). This may happen, for example, when harvesting more broadleaved species with higher wood density (therefore the unit of volume contains more biomass and more carbon), while the growth of living biomass is ensured mainly by coniferous species (with generally low wood density).
These differences originated in the way CBM reconstructed initial standing volume from user-defined yield tables, while EFISCEN used the exact values reported by the NFI. On forest types, the actual deviations of the initial standing stock were up to ± 20% for the values initialized by CBM. We identified two limitations that affected the fit of merchantable standing stock volume models: (i) eight out of 10 forests types were subject to shelterwood systems, therefore, with lower standing stocks for stands older than 100 years, for which reason the developed yield curves predicted systematically larger stock (assuming that standing stock was not reduced by the shelterwood system); (ii) the yield curves were derived from age class–dependent standing stock volumes per forest type and per owner type, available as region averages (at NUTS-2 level) as no detailed information at the NFI plot level was available. The NFI estimate resulted from much more detailed stratification on approximately 25 forest types [47, 51] compared with our structure of 10 forest types. A further simplification has to do with omitting the specific models for mixed forests, despite the existing knowledge on the growth of mixed tree species . Moreover, we did not include in the simulation the effect of specific regional  or local growth conditions as the site index. This omission may further affect to some extent the accuracy of growth and yield projections in both models at the more local scale. In fact, both models can deal with growth functions that would represent the share of each species in a mixed forest, but that is not seen as a mixed-forest model.
The C stocks reported here are within the typical range of standing C stock for most of the European forests . Comparatively, Bouriaud at el.  found that above-ground biomass in Romanian beech forests increased with stand age across all management types and treatments, reaching about 300 t dry biomass ha−1 equivalent to approximately 150 tC ha−1 at the age of 100 years, similar to our estimate. Similar estimates of standing stock were reported by the National Forestry Accounting Plan of Romania , where the yield table was based on models developed by Giurgiu and Draghiciu . We question whether the initial overestimation of the merchantable stock has any marginal impact on the accuracy of the CBM-based estimates. We speculate that the overestimated initialized standing stocks propagate until the end of the ongoing production cycle, while applying an inaccurate increment would affect both the first and the following cycles.
Furthermore, we could not detect any substantial effect of CBM annualization of the area in the initial simulation year (i.e., CBM divides the input area for each age class in 10 equal areas corresponding to a one-year time step).
To keep the required initialization data to a minimum, only the area and the mean growing stock volume per age class as reported by the NFI were retained in EFISCEN for the initial year of simulation. For the simulation period, the volume distribution over age classes (matrix columns) was generated by an empirically based function . The aggregation of all individual volumes to a nationally aggregated volume in a different way than the NFI may have caused the 1.5% overestimation in EFISCEN (compared with the NFI estimate).
Through the Results and Discussion sections, all elements related to merchantable stock in the initial year of the simulation pertain as validation effort. Using an additional evaluation, we assessed that the yield curves used in CBM practically matched the Romanian yield curves corresponding to site productivity between the third and fourth class .
Net annual increment
Because of the particular type of indicators in Romanian forestry addressing the volume of the entire above-ground woody biomass, it is practically impossible to validate NAI of merchantable part only as used in this modelling exercise, without making assumptions.
Despite differences regarding the initial merchantable standing stock, its relative rate of change was simulated similarly by both models, although EFISCEN showed a slightly higher increment rate. This may be related to the enhanced regeneration that takes place in EFISCEN leading to a younger age class distribution. Also, less area harvested by CBM, given higher merchantable standing stock, leads to overall older forests. It is worth noting that EFISCEN implements the volume-increment dynamics as a percentage of the growing stock . On the other hand, CBM requires user-defined volume-increment functions, so the accuracy as well as the consistency of yield and increment curves depend on the user-selected pre-processing operations, i.e., independent of the modelling frame of the CBM model.
The harvesting allocation algorithm had a noticeable effect on the estimates of all forest indicators because it affects age classes, thus also increment, thus also the sink, etc. Both models satisfied a demand of 23 million m3 yr−1. Harvest allocation also had an effect on merchantable stocks per species. In CBM, a low harvest rate of broadleaved forests, according to historical rate, resulted in the increasing proportion of broadleaves and the decreased proportion of coniferous forests in standing stock. After 50 simulated years of forest management, the standing stock was composed of relatively more broadleaved forests (i.e., of species with higher wood density) and fewer coniferous forests, contrary to what is simulated by EFISCEN. This also explains the higher C stock in 2060 in CBM compared with EFISCEN.
Conceptually, the implementation of harvest is very different in CBM compared with EFISCEN. EFISCEN applies a so called “free allocation”, while CBM uses a “detailed allocation”. The free-allocation approach allows EFISCEN to manage internally the split of total harvest demand on thinning and final fellings by species based on pre-defined management regimes per species and age class. The EFISCEN model then searches for its total requested demand, given the state of the forest, and allocates it across available forests types. In CBM, for an optimum projection of the forest-management interventions, the harvest needs to be projected by an external effort, e.g., allocation of harvest based on simulated future dynamics of standing stock, as in ref. , or through specific tools, e.g., the Remsoft Spatial Planning System as used in Canada . Most likely, the constant harvest scenario that we applied was satisfied given the relatively low share of the harvest rate from the total increment rate (some 66%) or from available standing stock subject to applicable silvicultural interventions (< 90%); therefore, enough wood resource was available across all forest types.
The way the model defines the harvest does have a significant effect on age classes and thus on simulation results (CBM defines the harvest as the amount of carbon targeted, i.e., having a “climate reporting”-oriented approach, whereas EFISCEN defines the harvest in terms of volume, having therefore a “forestry-centred” approach).
Both models showed an increase in forest age over time, whereby in CBM the ageing was much stronger and had less area shift in the regeneration class. Such a dynamic highlights the models’ internalized concept of allocation of harvest demand to forest-management practices, and especially how stand replacements, i.e., final fellings, are distributed according the availability of standing biomass by each of the two models (see Fig. 1).
Allocation of biomass in other compartments of the stands
Of importance in simulations of standing C stocks and CO2 removals were the non-merchantable biomass compartments (e.g., branches, bark, roots, foliage). Despite input-harmonization efforts, CBM simulated an annual average of 46–52% more biomass in these compartments compared with EFISCEN (Table 1). Still, the accuracy of biomass proportions may well be questioned in both cases. EFISCEN uses a straightforward approach in which a BEF value (specific to the forest age and type of non-stemwood biomass compartment) is applied directly to the standing volume parameterized by age and species to estimate the biomass of that compartment. CBM requires as input the proportions of the biomass compartments (i.e., stemwood, bark, branches and foliage) estimated as a function of merchantable volume using a simultaneous fit of all biomass compartments. It is noticeable that employing a simultaneous fitting approach ensures that the sum of compartment predictions equals the prediction of total biomass (so giving due consideration of compartments size within the whole architecture of the standing living biomass). However, both models add other compartments using proportion of the merchantable standing stock, so only the total standing biomass estimate is affected by these proportions.
The quality of the data used to fit these models is also essential; therefore, it is important that the users have access to appropriate data (measured based on robust sampling designs). Instead of using country-specific data to estimate the parameters of Boudewyn equations , another alternative is to use the CBM’s Canadian library to select the most suitable parameters , although this may result in strange correspondence mapping, e.g., coniferous to broadleaved species. Allocation in compartments could not be reasonably calibrated against any measured data from Romania, but we expect each model is self-consistent (proportions of other biomass compartments stayed rather constant in time).
The sink in standing living biomass
Overall, there was an enhanced, mutually related model effect on CO2 fluxes for the living biomass pool. For example, the annual sinks showed a 15% difference assuming averages during 2010 to 2060, i.e., 17.5 million tCO2 in EFISCEN versus 19.3 million tCO2 in CBM. Despite different but equally justifiable procedures, there is a cumulative effect when the small, apparently insignificant differences or percentages are applied to the relatively lower carbon stocks in EFISCEN versus the relatively higher carbon stocks in CBM. Such a combined arithmetical effect applies to NAI (e.g., 5% difference), the harvest level achieved (e.g., 2% underachievement of the harvest target by CBM, which means approximately 0.2–0.5 Mt CO2 yr−1), shares of other biomass compartments and changes in the contribution of forest types with different wood density in the total standing stock (e.g., impact ranges between 1.2 and 2.5 Mt CO2 yr−1). Apparently, this is in the range of uncertainty of the forest-sink estimates (usually around 20%) .
Mortality and dead wood carbon stock
These two parameters were used to validate the loss from living biomass against the NFI measured data. With regard to carbon stocks in dead wood pools in the initial year, Yasso07 relies on an equilibrium approach (steady state amount), while CBM performs a semi-equilibrium procedure (iteration of saturation levels through repeated disturbances and taking into effect the latest major silvicultural intervention). Default decomposition parametrization of the models was used, with the exception of mortality rate, which was harmonized between the two models. Despite an estimated larger carbon stock in living biomass in CBM (a larger amount of matter transferred to DOM expected given the higher turnover values), the initialization of dead wood pools was lower likely because of a faster decomposition rate in CBM. Uncertainty of C stock in dead wood is expected to increase substantially from standing to lying dead wood, as the NFI only provides volume information, e.g., no decomposition or density information. Obviously, using the same wood density of live tree species would result in overestimation of C stock.
Apparently, both models show an initialization issue with regard to mortality and standing dead wood. Mortality seems higher before the initialization than in simulations because of missing forestry operations and salvage of trees otherwise transferred to dead wood. This seems contrary to the fact that both dead wood stock and merchantable standing stock increased between the NFI-1 and the NFI-2.
The forest-management practices as sampled by the NFI were also simplified for the purpose of simulations. We used a business-as-usual interventions intensity and a regeneration delay after final cut (see Additional file 1: Appendix S2), capturing the patterns of forest types throughout age classes. The two models have a very similar solution to specify silvicultural practices; therefore, these were well harmonized. Nevertheless, our study excluded three types of potential changes that might have influenced the projections: (i) natural disturbances, (ii) wide application of shelter systems and (iii) change of forest composition. The most significant natural disturbance was the windthrow (a stand-replacing disturbance). Given the high rate of salvage logging following windthrow, the impact on age class structure is expected to be low. Therefore, the resulting living biomass pool following windthrow may not be significantly different from typical harvesting. However, we expect significant differences for dead organic matter input, but this is beyond the scope of our study.
The current versions of the models are not spatially explicit, so it is not possible to apply natural disturbances in places where they most likely occur (e.g., in forests affected by wrong past management decisions or forests prone to insect outbreaks). Instead, they are applied randomly across largely defined criteria (forest types and age classes). Furthermore, for comparability between the two models, no shelter systems (which is largely applied in Romanian forestry) and no changes in forest composition were applied (although available from the NFI), as they cannot be easily simulated by either CBM or EFISCEN.
Specification of deforestation
The two models apply different solutions that may result in either negligible or significant impact on projections depending on the magnitude of deforestation. CBM applies user-defined criteria for deforestation and accounts explicitly for losses in all carbon pools during deforestation, following the IPCC guidance for national GHG inventories . EFISCEN assumes that deforestation takes place after a final felling, so the impact on forest indicators may be significant if there is a different magnitude of deforested area. Nevertheless, this analysis excludes the carbon loss by deforestation in both models. The deforestation rate applied in this simulation seems to have a negligible impact on Romanian forests, as demonstrated by both models.
Further research, data needs and improvements
The availability of NFI data as regionally aggregated instead of plot-level might have affected the simulation effort. Thus, data and related metadata should be made available freely and openly to the scientific community to ensure good governance of forest and forestry information. Repeated measurements as well as improved modelling tools are further required to allow assimilation of most recent data in the modelling exercises .
Both models require pre-processing, and despite they are both open source, the pre-processing may limit the full reproducibility, i.e., if different models are used to fit the yield and increment. For this reason, more transparency is needed when models are improved or used.
A harvest-optimization tool is needed for CBM, which requires explicit harvest demand according to forest types for the simulated period. While this is a straightforward task for non-stand–replacing disturbances based on intensity of thinning interventions, it is a difficult decision for old stands, where forestry principles may result in various responses on forest dynamics: normalization/optimization of forest age structure (e.g., either on standing volume or area), or optimization of transfers to older protected forests (e.g., proportions of area excluded from wood production).
Improvements are already in progress for both models. The new EFISCEN-Space will be based on a modelling approach running on each NFI plot and its individual trees, accounting for information (such as tree densities) and individual tree data (such as diameter and height and running on climate sensitive growth functions). These NFI plot data will allow for better representation of mixed forests, uneven-aged forests, actual forest management and site-specific growth conditions, thereby making a climate-sensitive modelling approach possible. Refining the representation of climate-change impacts is the subject of ongoing research on both models. CBM moves to open-source versions while it strives to implement geospatial simulations [61, 62].