Quantifying forest growth uncertainty on carbon payback times in a simple biomass carbon model

In 2018 Sterman et al (2018a) published a simple dynamic lifecycle analysis (DLCA) model for forest-sourced bioenergy. The model has been widely cited since its publication, including widespread reporting of the model’s headline results within the media. In adapting a successful replication of the Sterman et al (2018a) model with open-source software, we identified a number of changes to input parameters which improved the fit of the model’s forest site growth function with its training data. These relatively small changes to the input parameters result in relatively large changes to the model predictions of forest site carbon uptake: up to 92 tC.ha−1 or 18% of total site carbon at year 500. This change in estimated site carbon resulted in calculated payback periods (carbon sequestration parity) which differed by up to 54 years in a clear-fell scenario when compared with results obtained using previously published parameters. Notably, this uncertainty was confined to forests which were slower growing and where the model’s training dataset was not sufficiently long for forests to reach maturity. We provide improved parameterisations for all forest types used within the original Sterman et al (2018a) paper, and propose that these provide better fits to the underlying data. We also provide margins of error for the generated growth curves to indicate the wide range of possible results possible with the model for some forest types. We conclude that, while the revised model is able to reproduce the earlier Sterman et al (2018a) results, the headline figures from that paper depend heavily on how the forest growth curve is fitted to the training data. The resulting uncertainty in payback periods could be reduced by either obtaining more extensive training data (including mature forests of all types) or by modification of the forest growth function.


Introduction
The use of biofuels is a significant element in global climate change mitigation strategies (Smith et al 2014, Rogelj et al 2018. These fuels have been attractive because, while biofuels do release CO 2 on combustion, the regeneration of biofuel feedstocks can result in the absorption of an equal quantity of CO 2 . This regenerative ability in biofuels (unlike fossil fuels which do not regenerate during any meaningful timeframe: Rose et al 2014) has historically led to an assumption of 'carbon neutrality' of energy from biofuels (Buchholz et al 2016, Smith et al 2014. While this assumption can be broadly true (assuming that the feedstock regenerates as expected, omitting supply chain emissions, and depending on the definition of 'carbon neutral': Malmsheimer et al 2011) a period of elevated atmospheric CO 2 may still exist because of the disparity between the rates of carbon emission and uptake relative to a counterfactual scenario. This has been described as a 'carbon debt' (Fargione et al 2008) site carbon is 'spent' on combustion, and takes a period of time to be 'repaid'. The time required to pay back the debt varies depending on the carbon emissions of the relevant supply chains, and recovery of the forest where the feedstock originated.
An extensive array of possible methodological choices (Buchholz et al 2016, Breton et al 2018 terminology (Malmsheimer et al 2011, Mitchell et al 2012, Bentsen 2017 and accounting assumptions (Lamers and Junginger 2013, Ter-Mikaelian et al 2015, Bentsen 2017 as well as the possibility of long periods of carbon debt, has led to a wide range of estimated payback periods (the time taken for a biomass scenario to outperform a counterfactual in terms of emissions over time: described as Carbon Sequestration Parity in Mitchell et al 2012) and conflicting opinions about the advisability of using forest biomass for energy. Variability in the results available from the literature, debate about the inclusion of different factors, and over-simplified communication of the benefits and costs of biomass use (Searchinger 2012, Searchinger and Heimlich 2015, Moomaw 2018a, Isaacs 2019) has been interpreted more widely as a lack of consensus. This has resulted in a range of responses in the policy sphere (e.g. FERN 2011, RSPB 2012, Brack 2017 which are not necessarily fully informed (Slade et al 2018).
In 2018  In testing our replication of the model, we identified uncertainties in the modelled forest growth curves, which have a profound impact on the reliability of results from five of the eight region/species combinations previously published by Sterman et al (2018a).

Research objectives
The work was carried out with reference to two key research objectives: Firstly, the paper by Sterman et al (2018a) has been widely cited in the academic literature, policy briefings/ recommendations and the press (as described above). We produced a model based on a similar framework which is simple/easy to modify; fully open source; and easy to calibrate to different assumptions. This tool was intended to facilitate comparisons between the wide range of different techniques, experimental decisions and assumptions which already exist in the literature.
Secondly, our work identified a level of uncertainty in the results as obtained by Sterman et al (2018a) that they did not discuss. We draw attention to this discrepancy, identify possible solutions, and caution against undue reliance on the affected results.

Method
Sterman et al (2018a) model analysis The Sterman et al (2018a) model is available under an open access licence, although it operates inside a closedsource framework which requires a full, 'professional', licenced copy of the VENSIM software (Venata Systems 2017) to work. The model is formed of two key components: (1) equations describing the supply chain, which are used to calculate the emissions associated with fuel production and use (a life cycle assessment: LCA calculation) and (2) a site model which calculates ongoing changes in carbon on a forest site and in the atmosphere after felling has taken place (a dynamic component to the LCA). While closer analysis of the supply chain variables and calculation may be warranted, here we focus on the forest site model.
The site model makes use of the flexible Chapman-Richards growth function (Richards 1959, Pienaar and Turnbull 1973, Zhao-gang and Feng-ri 2003 parameterised to represent aboveground carbon in forest sites. This function is used to represent gross primary productivity (GPP) of a forest site, while simple multiplication by dimensionless constants is used to approximate forest respiration; organic carbon deposition; and heterotrophic respiration.
The model was parameterised for eight different forest types in the USA (listed in table 1) by using a leastsquares non-linear regression method, which is common ('invariably' used: Burkhart and Tomé 2012, p 239) within forest model development. The method uses a solution-finding algorithm to fit a mathematical representation of the forest growth curve to representative data-points with the minimum degree of error. This provides an approximation of actual forest growth behaviour which can be used to predict forest growth values where actual data is not present (e.g. in different time-steps, or over extended time periods). In this case, the model was trained using average values for forest and soil carbon produced by the USDA (Smith et al 2006). This dataset has been constructed using a combination of sample plots, and interpolation of data using the FORCARB2 model (Smith et al 2006, p 13). While uncertainties exist, and the data is unsuitable for smaller scale (stand-level) modelling, (Smith et al 2006, p 17)  This approach has limitations, specifically resolution and predictive capacity. Generalisation implicit in the training data (average forest growth over a wide range of local site types) limits the applicability of the model to individual stands or spatially explicit areas -providing only high-level estimates. The lack of process-based calculation limits the model's ability to take account of varying site conditions (climate, soil, exposure), complex silvicultural systems (thinning, or selection fellings), and stochastic factors (e.g. fire or pest outbreak). This method does, however, have value as a tool to look at the approximate effect of management of large areas of forest to simple silvicultural systems in a more abstract/theoretical setting. This is particularly useful when looking at the carbon dynamics of biomass where a number of misconceptions persist (Ter-Mikaelian et al 2015, Bentsen 2017).

Replicating the Sterman et al (2018a) model
We recreated the model developed in Sterman et al (2018a) in a general-purpose, object-oriented, high-level programming language: Python. This was carried out using the Anaconda Python distribution (available from www.anaconda.com) which is already widely used in the research community and is distributed with an extensive ecosystem of open-source libraries and modules.
Our implementation: the 'Simple Biomass Comparison Model' (SBCM: included as supplementary material, most recent version downloadable from github.com/Priestley-Centre/SBCM) follows the same pattern as described above. SBCM compares emissions over time resulting from 1) fossil fuel production and use, and 2) biomass production, use, and subsequent forest site recovery.
The model was initially compared to the results published by Sterman et al (2018a) to check the success of the replication, this was carried out using a scenario based on a clear-fell of forest in order to supply biomass (equivalent to scenario S3 in Sterman et al 2018a). The biomass supply scenario was compared with a counterfactual (fossil fuel) scenario with energy derived from coal with an implicit assumption that there would be no emissions from the forest site in the absence of biomass production. While this may not always be an appropriate comparison  Agreement between SBCM and the Sterman et al (2018a) model (using the same parameters) in terms of forest growth over the first hundred years is very good (<1tC.ha −1 ) but not exact in terms of maximum values for carbon over extended time periods (varying by up to 36tC.ha −1 for above-ground carbon at year 500). In view of this uncertainty, a curve fitting approach was developed in Python software to see whether improvements to the model parameterisation could be made, and whether a better fit of the training data would result in a closer agreement with Sterman et al (2018a).
Sterman et al (2018a) use the [optimizer] function in Vensim to apply least squares non-linear regression and Markov Chain Monte Carlo methods to determine model parameters. This is described in detail in their supplementary material (Sterman et al 2018a available online at stacks.iop.org/ERC/2/045001/mmedia) but, in summary, they restricted the matching algorithm to parameters which generate values which satisfy two set conditions: We used a similar approach in our replication of the model, using the [scipy] Python library. The [scipy. optimize] function for Python contains a number of different algorithms for curve fitting (table 2) and for the weighting given to outlying data (table 3). The training data produced by the USDA (Smith et al 2006) which in turn were based on a series of estimates from the Forestry Inventory and Analysis Database (USDA Forestry Service 2005) contains values for forest and soil carbon for forests over the first 90 or 125 years of growth (depending on forest type). This contains the full growth curve (planting to maturity) of plantation forests; however, natural forests (and forest soils) take far longer to reach equilibrium and data over this timescale was not present. This projection of the growth curve beyond the length of the fitting data results in additional uncertainties (as discussed below).
A full range of possible combinations of algorithm and loss function were attempted in two permutations: firstly, with the model parameters unconstrained (simply looking for the best fit possible, with very loose limits on possible parameter values) and, secondly with constraints applied-requiring the first value in the results to equal the first value in the USDA training data (±1 tC.ha −1 ). In both cases the fit of the modelled output to the training data was assessed by calculating the RMSE.
Using a subset of possible parameterisations (40 of 240) which resulted in an improved RMSE score over Sterman et al (2018a) SBCM was run to determine: the changes over time in forest and soil carbon storage, the carbon storage when mature (at 'equilibrium' as discussed by Sterman et al 2018a-assumed to be reached when forest and soil carbon is 99% of the potential maximum value) and the effect that parameter changes have on payback period.
In each case, the supply-chain model used the original parameterisation for supply chain efficiencies and emissions (from Sterman et al 2018a) without incorporating the modification of the efficiency parameter used by Dwivedi et al (2019) to allow a comparison with the original model results. Trust region reflective (Branch et al 1999). dogbox A trf implementation using a rectangular trust region (Voglis andLagaris 2004, Nocedal andWright 2006). Table 3. Loss functions providing different weight to outliers (Scipy.org 2019).

SBCM evaluation
Results for forest growth obtained by running SBCM using the standard parameters in Sterman et al (  and soil carbon predicted by SBCM (which is unsurprising due to the numerical instability of the Chapman -Richards growth function: Ratkowsky 1983, Burkhart andTomé 2012). In five of the eight forest types assessed, revised parameterisation resulted in substantial expansions of the uncertainty associated with site carbon, particularly in levels of soil carbon (see figure 5, which shows the range of possible outcomes using improved parameterisations). It is notable that this increase in uncertainty is confined to slower growing species mixtures  where the entire growth curve was not present within the training dataset. In contrast, all three plantation forests show very similar outcomes regardless of parameterisation ( figure 5 below). This supports the assertion in Burkhart and Tomé (2012) that the Chapman-Richards growth curve can result in more accurate outcomes when an estimated asymptote is added, as these forest types contain values close to the asymptote within the training data.

Changes in payback periods
Taking the range of growth curves into account and running the model for the full range of improved parameterisations (again using a 95% clear-fell scenario: equivalent to Sterman et al 2018a scenario S3) we found that the range of possible payback periods expands (figure 6). In each case, the values obtained using the Sterman et al (2018a) parameterisation fall at the top of the range of possible outcomes, and multiple estimates of shorter payback periods also fit the available data to a comparable degree. This is due primarily to the rate of emission to the atmosphere from the soil carbon pool. The rate of emission is directly proportional to the carbon pool size and, since this size and the time required for the pool to saturate are uncertain, this has a strong effect on the net emissions of the system shortly after felling.
There is a considerable difference between forests labelled by Sterman et al (2018a) as 'natural' (predominantly naturally regenerated hardwoods) and plantations (planted pine forests in the southern USA). Natural forests showed greater variation in terms of the time taken to reach maturity, the embodied carbon in a mature forest site, and the payback times associated with their use for biomass fuel. In contrast, plantation forests exhibited very low levels of uncertainty throughout. We argue that this pattern is due to extension of growth curves beyond the end of the fitting data. Of the forest types studied, all of the natural forests required long periods of time to reach maturity (equilibrium). The length of the datasets available to train the model were either 90 or 125 years, and this often fell short of capturing the entire growth curve. The faster growing plantation forests reached maturity comfortably within the timeframe available. A summary of the dataset length and years required to reach maturity is shown in table 4.
A further illustration of this effect can be seen in figure 7. There is a statistically significant correlation (p=6.26×10 −6 ) between the degree to which data has been extrapolated and standard deviation of the site carbon at equilibrium.

Discussion and conclusions
While there are small discrepancies between the forest site model results from the newly developed SBCM and Sterman et al (2018a) the overall match is good (figure 1), and we conclude that SBCM replicates the earlier model well.
In investigating disagreements between previously published results (Sterman et al 2018a) and the output of SBCM, we were not able to adjust the model parameterisations to achieve a large improvement in the fit with the training data ( figure 4). It is significant, however; that the small improvements that were achieved, resulted in substantial changes to carbon uptake estimates generated by the model (figure 5). New parameterisations revealed large uncertainties in terms of the time carbon stocks take to reach equilibrium as well as the quantity of carbon stored on forest sites. This is believed to be due to the sensitivity of the growth function used within the model to minor changes in input parameters (numerical instability: discussed in Ratkowsky 1983). It is notable that the wide variation in possible growth rates exists only in the forests described by Sterman et al (2018a) as 'natural'. The common feature in these forests is that the full growth curve is not fully represented within the training data, as forest growth to maturity takes longer than the 125 years available from Smith et al (2006). This failure to specify a value on or near the asymptote of the curve is known to provide less accurate results when using the Chapman-Richards growth function (Burkhart and Tomé 2012). Faster growing 'plantation' forests do not suffer from this problem: the entire growth curve (including the asymptote) is contained within the training dataset; and consequently, the variation in output is very low. We conclude that the modelled outputs for these site types are more reliable and less prone to error than those from slower  growing forests. In all forest types, the degree of variation is well correlated with the degree to which the growth curve has been extrapolated beyond the available training data (as shown in figure 7). The variability exhibited by the 'natural' forest types is significant given the reliance of the original model on the assumption that forests are biologically mature (at equilibrium) before felling and that this has occurred by year 500. In the absence of an extended input dataset, we conclude that it is not possible to determine which parameterisation for 'natural' forests is most appropriate when using the current growth function (indeed, the values preferred by Sterman et al (2018a) may be the best fit with real-world situations).
We acknowledge the argument by Prisley et al (2018) that the lack of more nuanced silvicultural systems is a weakness in the model. We suggest that this, when combined with the uncertainty inherent in the slower growing 'natural' forests, raises important questions about the ability of the model to accurately predict payback times in these forest types. Based on our range of possible parameterisations, the range of possible carbon payback periods when modelling 'natural' forests expanded by between 21 and 54 years. This is in contrast to range of possible payback periods when modelling 'plantation' forests, which in each case was within one year (figure 6). We found that in every case the values published by Sterman et al (2018a) represented the maximum payback period achieved, suggesting a potential bias in their parameterisation. We conclude that, while the paper produced by Sterman et al (2018a) has been widely cited and discussed, its headline figures contain a level of uncertainty which is not apparent. This uncertainty can be limited by use of a growth function which does not suffer from the same instability, or collection of an extended dataset for very old forests (to provide an asymptote, limiting the range of possible growth curves produced).
We recommend the inclusion of robust error/uncertainty reporting (included in tables 5a and 5b) when discussing the climate effect of biomass fuel production from these site types when using this model. We have provided a summary of our best parameterisations below (table 6) and these are also available in the supplementary information.
This work sits within a wider programme of research to address the large uncertainties in the literature when calculating the carbon dynamics of forest growth and biomass fuel use. This work is ongoing, and is being carried out through expansion of the SBCM model framework. In the short term we intend to mitigate the uncertainty addressed in this paper by restricting assessments to either time horizons which do not extrapolate beyond the training data (max 125 years) or are based on forest types which are not affected (the 'plantation' forests). Future work will be based around incorporation of an alternative growth function, development of a wider range of silvicultural options (   Recommended margins of error and parameters.