Biophysical and Economic Analysis of Black Spruce Regeneration in Eastern Canada Using Global Climate Model Productivity Outputs

We explore the biophysical potential and economic attractiveness of black spruce (Picea mariana) regeneration in eastern Canada under the high greenhouse gas emission scenario (RCP 8.5) of the Intergovernmental Panel on Climate Change (IPCC). The study integrates net primary productivity and net ecosystem productivity estimates from three major global climate models (GCMs), growth and yield equations specific to black spruce, and economic analyses to determine spatially varying investment values of black spruce regeneration—both including and excluding carbon sequestration benefits. Net present value (NPV) was used to represent financial attractiveness. It was assumed that stands would not be harvested at volumes less than 80 m3·ha−1. A baseline case with the stumpage price set to $20 m−3, stand establishment cost $500 ha−1, and the discount rate 4%, was used with wide-ranging sensitivity analyses conducted around these assumptions. These values represent the wide range of choices and outcomes possible with black spruce regeneration investments. The results indicated a latitudinal gradient in economic attractiveness, with higher forest productivity and NPVs in the southern portion of the study area; however, black spruce regeneration was not economically attractive unless regeneration costs were very low (representing something closer to a natural regeneration type scenario) and/or carbon sequestration benefits of at least $5 ton−1 CO2 were realized. In general, the optimal harvest rotation age increased with increasing carbon price by approximately 9 to 18 years. Results of this study highlight the importance of future price and yield expectations and establishment costs in evaluating forest investments.


Introduction
Boreal forests cover approximately 1.1 billion hectares (ha) worldwide and store about 53.9 Pg C as live biomass and 217.5 Pg C as organic matter in litter, deadwood and soil [1,2].They are a significant carbon sink with annual net primary productivity (NPP) of 2.6 Pg•C•year −1 [2,3].The productivity and carbon storage capacity of boreal forests is impacted by climate change-related stresses, which are more severe in boreal (i.e., northern) regions, and include extreme temperatures, droughts, wildfires and insect outbreaks [2][3][4][5].Therefore, it is important to better understand responses of boreal forests to climate change pressures over the course of this century, from both timber production and carbon sequestration perspectives.Several studies have evaluated the economic feasibility of afforestation and reforestation in North America, but assessments that attempt to account for carbon sequestration through direct linkages to growth processes and their dependence on seasonal and inter-annual climate variations are lacking.Past studies (e.g., [6][7][8][9][10]) tend to utilize generalized forest growth and yield curves, which do not account for climate variability.Excluding climate variability could be problematic because seasonal and inter-annual variations, as well as long-term trends due to climate change, are expected to impact future forest productivity [11].
Several studies have examined the effect of carbon revenues on harvest timing and stand values, when both timber and carbon sequestration benefits are considered [6,[12][13][14][15].These studies show that additional revenue flows (when carbon has an economic value) generally result in slightly delayed forest harvests, thus increasing the optimal harvest age [16][17][18].The permanence of carbon sequestered (or lack thereof) ultimately affects the price a forest owner would receive for sequestered carbon.Given that carbon markets are not well established in many jurisdictions across the world including Canada, there is significant interest in better understanding how price levels may influence carbon-related forest investment decisions.However, there is also large uncertainty in future forest productivity, hence there is a need for tools and/or models to integrate biophysical potentials and economic values.
In this study, we evaluate the economic feasibility of afforestation and reforestation projects under a changing climate.We focus on black spruce (Picea mariana) regeneration in the province of Ontario, Canada.Black spruce is a dominant species in the boreal region of this province (and in Canada as a whole)-covering approximately 41 percent of the province's forested area and accounting for up to 80% of the annual allowable forest cut [19].The specific objectives of this study are to (i) generate future merchantable timber yield and carbon sequestration estimates for black spruce, utilizing GCM-based NPP and net ecosystem productivity (NEP) estimates that account for projections of climate variability; (ii) examine the economic potential of black spruce regeneration, estimated using net present value (NPV) analyses both inclusive and exclusive of carbon sequestration benefits; and (iii) conduct sensitivity analyses to determine how uncertainty in the economic inputs affects economic attractiveness of black spruce forests.

Overview
The study focused on the province of Ontario in eastern Canada, covering the area between 40 • and 52 • N and 75 • and 92 • W (Figure 1), with a northern boundary that roughly coincides with the northern limit of the province's active forest management zone [19].The study area has been divided into three regions-northwest, northeast, and southern Ontario as shown in Figure 1.Although far southern Ontario is not part of the active forest management zone, or part of black spruce's natural range, the species can in principle grow in this region; thus, it is included to help illustrate variation in economic values across a broad spatial domain.
The analysis required the development and integration of data from several observation and modelling sources.Annual forest growth (i.e., NPP) and carbon sequestration (i.e., NEP) estimates were obtained over the 2006 to 2100 period for three widely used GCMs and one emissions pathway from the Coupled Model Intercomparison Project phase 5 (CMIP5) initiative as described in the Intergovernmental Panel on Climate Change (IPCC) AR5 report [20].The NPP values were used to modify merchantable timber volume estimates obtained from a black spruce growth and yield model, while the NEP values provided annual estimates of carbon sequestered by the growing forest.A baseline scenario of costs, wood and carbon values (prices) and a discount rate to reflect the time value of money were used to calculate NPV (in dollars/hectare) of regeneration investments.NPV calculations followed the Faustmann model [21] for the wood production value calculations and a modified Hartman model for the carbon value calculations [6,22].A broad set of economic

GCM-Simulated Forest Productivity Estimates
Besides commonly used projections of temperature and precipitation, GCMs provide a number of carbon-related outputs.Photosynthetic carbon assimilation is expressed as gross primary production (GPP) and this assimilated carbon is used as a source of energy for trees to grow [23].Plant biomass accumulation is represented by NPP, which is calculated as GPP minus carbon losses through respiration from live plant components or autotrophic respiration [23,24].The net accumulation of carbon by an ecosystem is represented by NEP, which is the difference between NPP and carbon loss through decomposition of dead organic matter or heterotrophic respiration [20].NEP is often used to estimate carbon sequestration for a given area and time period, assuming disturbance or dissolved organic carbon related losses are negligible.Note that these carbon-related GCM outputs incorporate the projected impact of increased levels of atmospheric carbon (i.e., CO2 fertilization) on plant growth and respiration.
As previously noted, estimates of NPP and NEP were downloaded for three GCMs (i.e., CanESM2, MIROC, and MPI-ESM).(i) The Canadian Earth System Model (CanESM2) includes a fourth generation atmospheric general circulation model (CanCM4), a physical ocean component (OGCM4), the Canadian Model of Ocean Carbon (CMOC) and a process-based dynamic vegetation model, known as the Canada Terrestrial Ecosystem Model (CTEM) [25,26].Model simulations were performed at a rather coarse resolution of 2.8° × 2.8° (Figure 1a); (ii) The Japanese Model for

GCM-Simulated Forest Productivity Estimates
Besides commonly used projections of temperature and precipitation, GCMs provide a number of carbon-related outputs.Photosynthetic carbon assimilation is expressed as gross primary production (GPP) and this assimilated carbon is used as a source of energy for trees to grow [23].Plant biomass accumulation is represented by NPP, which is calculated as GPP minus carbon losses through respiration from live plant components or autotrophic respiration [23,24].The net accumulation of carbon by an ecosystem is represented by NEP, which is the difference between NPP and carbon loss through decomposition of dead organic matter or heterotrophic respiration [20].NEP is often used to estimate carbon sequestration for a given area and time period, assuming disturbance or dissolved organic carbon related losses are negligible.Note that these carbon-related GCM outputs incorporate the projected impact of increased levels of atmospheric carbon (i.e., CO 2 fertilization) on plant growth and respiration.
As previously noted, estimates of NPP and NEP were downloaded for three GCMs (i.e., CanESM2, MIROC, and MPI-ESM).(i) The Canadian Earth System Model (CanESM2) includes a fourth generation atmospheric general circulation model (CanCM4), a physical ocean component (OGCM4), the Canadian Model of Ocean Carbon (CMOC) and a process-based dynamic vegetation model, known as the Canada Terrestrial Ecosystem Model (CTEM) [25,26].Model simulations were performed at a rather  1a); (ii) The Japanese Model for Interdisciplinary Research on Climate (MIROC) consists of an atmospheric general circulation model (MIROC-AGCM), an aerosol component (SPRINTARS v5), an ocean GCM with sea-ice components (COCO v3.4), and a land surface model (MATSIRO) [27].MIROC shares the same spatial resolution as CanESM2 (Figure 1a); (iii) The German Max-Planck-Institute-Earth System Model (MPI-ESM) consists of an atmospheric general circulation model (ECHAM6), an ocean model (MPIOM), marine geochemistry model (HAMOCC5), and a dynamic vegetation model describing physical and biogeochemical aspects of soil and vegetation (JSBACH) [28,29].The spatial resolution of the MPI outputs is 2.5 • × 2.5 • (Figure 1b).
Simulations performed by all three models used a common experimental protocol prepared for CMIP5 participating models [30].This protocol included carbon fluxes from natural vegetation and soil as well as some measure of disturbance such as land use change.Performance of these models in simulating carbon fluxes and related components for historical (1850-2005) and future (2006-2100) time periods has been evaluated previously using a moderate greenhouse gas emission scenario, i.e., Representative Concentration Pathway 4.5 (RCP 4.5) [31].Here, we employ a high emission scenario, RCP 8.5, as emissions are currently tracking at or above this level [32].

Growth and Yield Estimates and NPP Adjustments
Yield equations for black spruce, that incorporated broad-scale climate influences, were obtained from Ung et al. [33]: where V represents gross total volume of live merchantable trees (m 3 •ha −1 ), t represents stand age (years), D represents mean annual temperature ( • C), P represents mean annual precipitation (mm), C d represents the Duan correction factor [34] and β 0 and β 1-5 are fitted parameters.Despite including terms for temperature and precipitation in the model, Ung et al. [33] warn that their equations are not intended for use under climate change as the models could be pushed outside the range of values for which they were originally developed.As evidence of this, preliminary work with Equation (1) produced unrealistically high black spruce yields (i.e., greater than 1000 m 3 •ha −1 ) as mean annual temperatures were pushed beyond typical historical limits.Thus, while the Ung et al. [33] yield equation provided reasonable estimates under current climate conditions (i.e., ~200-250 m 3 •ha −1 ), another approach was needed to incorporate climate change.
We developed an approach that incorporated both GCM-simulated NPP and the Ung et al. [33] yield equation for black spruce.Because NPP is a direct estimate of plant productivity, it follows that NPP will reflect black spruce growth (and hence yields) through time.We first calculated an NPP modifier (R NPP ) for each grid cell and year over the 2011-2100 period as a ratio of future to current NPP: where NPP baseline is the GCM-simulated annual average NPP estimate for the 1971 to 2000 period, and NPP predicted is the GCM-simulated NPP estimate for each year over the 2011 to 2100 period.R NPP ranges from 0 to ∞ where a value between 0 and 1 indicates a decrease in NPP and any value higher than 1 indicates an increase in NPP.R NPP is multiplied by the change in volume over each year and then added onto the total volume accumulated in previous years: where V t is total yield (m 3 •ha −1 ) at year t, V Ung is the volume calculated using Equation (1) for stand age at year t and year t − 1, and V t−1 is accumulated volume in previous years.Note that the V Ung values were calculated using temperature and precipitation estimates for the 1971-2000 period; these Forests 2017, 8, 106 5 of 18 estimates were obtained from previously developed spatial climate models (see McKenney et al. [35] for details).Annual yields were estimated for all grid cells (~10 km resolution) across the study area, based on a plantation initiation date of 2016, but modified through time using the NPP ratio modifier.

Economic Calculations
NPV was calculated using two different harvest models-the basic Faustmann model, which considers only revenues from wood production [21,36], and a modified version of Faustmann, the so-called Hartman model, which can in principle include both non-wood value flows (e.g., carbon sequestration revenues) and wood production revenues [22].Because the exact price of carbon is still unknown or yet to be defined in Canada, we used a range of values from Canadian $5 to $20 ton −1 •CO 2 .These values reflect fledgling carbon prices in North America, which have ranged from $3.36 (US $3) ton −1 •CO 2 at the Chicago Climate Exchange in 2006 [9] to $14 (US $12.5) ton −1 •CO 2 at the California Carbon Dashboard in 2016 [37].Note that, for summary purposes, the NPV results were organized into three broad regions-northwest, northeast, and southern Ontario (Table 2).

Wood-Only Calculations
The Faustmann model provides the theoretical basis for stand-level optimal harvest age decisions when only wood fibre values are considered [38,39]: where NPV is the net present value calculated on a representative per hectare basis ($•ha −1 ), P is the stumpage or standing timber value ($•m −3 ), V(t) is the yield (m 3 •ha −1 ) obtained from Equations ( 1)-( 4), C is the establishment cost ($•ha −1 ), r is the discount rate or time value of money (%), and t is the age of the forest.The optimal economic harvest age is the age at which the NPV is maximized.
For the baseline scenario, the price of stumpage is set to $20 m −3 , stand establishment cost is set to $500 ha −1 , and the discount rate is set to 4% [8,9].Furthermore, the minimum harvest volume was set at 80 m 3 •ha −1 for the simulations, which reflects the approximate lower limit of commercial viability in the study region.

Wood and Carbon Sequestration Calculations
Hartman [22] provided a formulation for economic assessments of the optimal harvest age of a forest stand when both timber value and flows of other services (e.g., carbon sequestration) are considered: where PV C is the present value of sequestered carbon, P cs is the carbon price ($•ton −1 •CO 2 ), F(t) is the carbon sequestration function, F (t) is the integral of F(t), β c is the conversion factor (set to 0.3 for the current work) to reduce the total carbon mass to biomass volume which is consistent with an estimated carbon content of wood of approximately 200 kg•m −3 [40], r is the discount rate (%), and t is the age of the forest stand.The term to the left of the minus sign in Equation ( 6) represents the NPV of carbon stored over the period of the standing trees, whereas the term to the right represents the carbon value loss when the stand is harvested, as forest biomass is expected to be set to zero at the time of harvest [15].The mass of carbon in standing trees (biomass: ton•ha −1 ) as a function of stand age, is given by: where v 1 , v 2 , and v 3 are fitted parameters with values of 189.6, 0.0268, and 2.56, respectively.This equation was parameterized by fitting GCM-simulated NEP data over the 1901 to 2100 period.
Inserting Equation ( 6) into the Faustmann model (Equation ( 5)) gives the NPV for both timber production and carbon sequestration: We reiterate that the optimal rotation age is that at which NPV is maximized.Importantly, Hartman showed that the "optimal" decision may be to leave a stand unharvested, depending on details such as growth rates, prices of the values under consideration, and the relevant discount rate.

Sensitivity Analysis
Sensitivity analyses are a useful tool to examine possible outcomes when uncertainty is inherent in the problem formulation.In this case, we have used them to explore the response of NPV to parameter uncertainty.For each model parameter, simulations were repeated with the parameter adjusted to represent known uncertainty in these values, while all other parameters were maintained at their original values (Table 1).The zero establishment cost is intended to represent the possible example of natural regeneration occurring with no additional stand management activities after an area is harvested.Ultimately, regeneration costs will vary according to a managers objectives and a stand's condition after harvest.

Climate, NPP, and Growth Projections
Figure 2 provides temporal trends in GCM-simulated temperature, precipitation and NPP for each grid cell across the study area over the course of the current century.An increasing temperature trend was projected by all three models, although the rate of temperature increase varied slightly among the GCMs.Over the study area, the average rate of temperature increase was 0.76, 1.03 and 1.16 • C per decade as simulated by CanESM2, MIROC and MPI models, respectively (p-values < 0.05).Temperature projections by the MPI model were slightly lower than those of the other two models.
There was also a slight increasing trend in the precipitation projections, though this is largely masked by significant spatial and temporal variability in the precipitation estimates.The average rate of precipitation increase was 0.33, 1.20 and 0.63 mm per decade as simulated by CanESM2, MIROC and MPI models.
CanESM2 and MPI exhibited gradual increases in NPP over time of 1.53 and 2.17 g•C•m −2 per decade, respectively (p-values = 0.18 and 0.05), while MIROC-simulated NPP showed a curvilinear relationship over the course of the century (Figure 2f).Liu et al. [41] reported an average NPP in eastern Canada of approximately 259 g•C•m −2 year −1 ; the GCM-simulated NPP values reported here are smaller in comparison, reaching a maximum of ~170 g•C•m −2 •year −1 .This issue may arise because the GCMs use single class representations of land cover to describe large grid cells, which can lead to prediction errors [42].3).Pearson's linear correlation analysis indicated a positive correlation between temperature and NPP as well as between precipitation and NPP as reported in previous studies [43][44][45].ANOVA indicated that 43% of the variability in NPP is explained by temperature and precipitation (29.7% and 12.7% respectively).Experimental studies have also reported that carbon assimilation by black spruce may increase due to climate warming [46,47], although it is possible that such a response may be inhibited by other warming-induced processes affecting plant growth, such as drought-related impacts on photosynthesis [48][49][50].
Similar to NPP, projections of merchantable wood volume were highest in the southern part of the study area and lowest in the north (Figure 4).The projected average yield volumes in 2100 for CanESM2, MIROC, and MPI in southern Ontario were 270, 277, and 253 m 3 •ha −1 , respectively.In contrast, the projected average yield volumes in 2100 for CanESM2, MIROC, and MPI in northern Ontario were 170, 170, and 156 m 3 •ha −1 , respectively.Again, this north-south productivity gradient is likely driven by a combination of higher NPP, mean annual temperature, and annual precipitation in  3).Pearson's linear correlation analysis indicated a positive correlation between temperature and NPP as well as between precipitation and NPP as reported in previous studies [43][44][45].ANOVA indicated that 43% of the variability in NPP is explained by temperature and precipitation (29.7% and 12.7% respectively).Experimental studies have also reported that carbon assimilation by black spruce may increase due to climate warming [46,47], although it is possible that such a response may be inhibited by other warming-induced processes affecting plant growth, such as drought-related impacts on photosynthesis [48][49][50].
Similar to NPP, projections of merchantable wood volume were highest in the southern part of the study area and lowest in the north (Figure 4).The projected average yield volumes in 2100 for CanESM2, MIROC, and MPI in southern Ontario were 270, 277, and 253 m 3 •ha −1 , respectively.In contrast, the projected average yield volumes in 2100 for CanESM2, MIROC, and MPI in northern Ontario were 170, 170, and 156 m 3 •ha −1 , respectively.Again, this north-south productivity gradient is likely driven by a combination of higher NPP, mean annual temperature, and annual precipitation in the south [51][52][53].Our yield estimates, are comparable to (though slightly higher than) published yield values for Ontario, which generally range between 120 and 250 m 3 •ha −1 at 80-100 years of age [54][55][56].

Economic Benefits
Based on wood value alone, black spruce regeneration was not economically attractive across the study area under our baseline assumptions (Figure 5, Table 2).At the optimal economic harvest ages, the average net present values for CanESM2 in the northwestern, northeastern, and southern regions of the study were −$136, −$265, and −$51 ha −1 , respectively.Similarly, NPVs for these regions using MIROC were −$137, −$269, and −$39 ha −1 ; and using MPI were −$226, −$263, and −$90 ha −1 , respectively.The optimal harvest ages for CanESM2, MIROC, and MPI were 39, 38, and 38 years, respectively, for the baseline scenario.These rotation ages are well below the traditional rotation ages for this species, which are typically determined using the maximum sustained yield (MSY) approach.This disparity reflects the fact that the time value of money is not considered in the MSY approach, only biophysical yields (see discussion in Yang et al. [57]).Interestingly, despite the decline in projected NPP in the second half of the century by the MIROC model (see Figure 3) both optimal harvest age and NPV did not vary much from the other two GCMs.This is due to the discount rate,

Economic Benefits
Based on wood value alone, black spruce regeneration was not economically attractive across the study area under our baseline assumptions (Figure 5, Table 2).At the optimal economic harvest ages, the average net present values for CanESM2 in the northwestern, northeastern, and southern regions of the study were −$136, −$265, and −$51 ha −1 , respectively.Similarly, NPVs for these regions using MIROC were −$137, −$269, and −$39 ha −1 ; and using MPI were −$226, −$263, and −$90 ha −1 , respectively.The optimal harvest ages for CanESM2, MIROC, and MPI were 39, 38, and 38 years, respectively, for the baseline scenario.These rotation ages are well below the traditional rotation ages for this species, which are typically determined using the maximum sustained yield (MSY) approach.This disparity reflects the fact that the time value of money is not considered in the MSY approach, only biophysical yields (see discussion in Yang et al. [57]).Interestingly, despite the decline in projected NPP in the second half of the century by the MIROC model (see Figure 3) both optimal harvest age and NPV did not vary much from the other two GCMs.This is due to the discount rate, which reduces    When wood and carbon values were considered, NPVs were positive across the study area (Figure 6; Table 2).Average NPVs at their optimal harvest age for the baseline scenario for CanESM2 in northwestern, northeastern, and southern Ontario were $430, $305, and $522 ha −1 , respectively.NPVs for these regions using MIROC were $424, $304, and $529 ha −1 , respectively, and using MPI were $341, $306, and $485 ha −1 , respectively.A price of $5 ton −1 CO 2 appears to be a minimum threshold at which black spruce regeneration investments start to become attractive (see also Yemshanov et al. [8]).This suggests that the overall investment results rely heavily on the carbon price assumptions/expectations.We note the maps should only be used as general indicators of relative potential value due to the coarse scale of the spatial model inputs used in the study.Optimal harvest ages ranged from 42-44 years for the three GCMs, which were slightly longer than those identified in the wood-only analysis.In comparison, Yemshanov et al. [8] reported an optimal rotation age of 49 years for coniferous forests in eastern Canada inclusive of carbon benefits.When wood and carbon values were considered, NPVs were positive across the study area (Figure 6; Table 2).Average NPVs at their optimal harvest age for the baseline scenario for CanESM2 in northwestern, northeastern, and southern Ontario were $430, $305, and $522 ha −1 , respectively.NPVs for these regions using MIROC were $424, $304, and $529 ha −1 , respectively, and using MPI were $341, $306, and $485 ha −1 , respectively.A price of $5 ton −1 CO2 appears to be a minimum threshold at which black spruce regeneration investments start to become attractive (see also Yemshanov et al. [8]).This suggests that the overall investment results rely heavily on the carbon price assumptions/expectations.We note the maps should only be used as general indicators of relative potential value due to the coarse scale of the spatial model inputs used in the study.Optimal harvest ages ranged from 42-44 years for the three GCMs, which were slightly longer than those identified in the wood-only analysis.In comparison, Yemshanov et al. [8] reported an optimal rotation age of 49 years for coniferous forests in eastern Canada inclusive of carbon benefits.We note that a number of studies have indicated that water limitations and heat stress may limit the ability of black spruce to grow and sequester carbon, but there is no general consensus on the We note that a number of studies have indicated that water limitations and heat stress may limit the ability of black spruce to grow and sequester carbon, but there is no general consensus on the impacts of climate warming on the productivity of black spruce and forests in general [49,[58][59][60].Girardin et al. [45,50] looked at the impacts of climate warming, drying, and increasing CO 2 concentrations on the productivity of black spruce forests and reported that inter-annual variability in black spruce productivity is significantly driven by soil water availability across broad areas of the western to eastern Canadian boreal forest and by autotrophic respiration in warm southern boreal regions.Thus, it is plausible that pending climate change could result in a decrease in carbon uptake and affect actual outcomes in a manner not captured by the GCM-based NPP adjustments utilized here.The general issue of predicting forest productivity under a rapidly changing climate remains an important research topic.

Sensitivity Analyses
Sensitivity analyses were performed for both the wood-only and wood + carbon scenarios to determine the effect of economic uncertainty on model results.Table 3 summarizes the sensitivity analyses of the wood-only scenario for a range of r, C, and Stumpage Price (P) values.When r was increased from 4% to 8%, the optimal rotation period was shortened by about 6 years for all three GCMs.Furthermore, because of this, NPV was significantly lower at the new adjusted optimal harvest ages.The largest relative drop in NPV was observed in the southern region of the study area, whereas the smallest reduction was seen in the northeastern region.Conversely, when r was decreased from 4% to 2%, rotation periods lengthened by 7-8 years and NPVs increased by up to $1086/ha across the three GCMs.Again, NPV changes were greatest in southern Ontario, reflecting the higher growth potential in this region.
Establishment costs (C) were adjusted from $500 to $0, $200 and $1000 per hectare.When C was reduced to $0, the rotation periods were shorter: from 39 years to 33 years for CanESM2, from 38 years to 34 years for MIROC, and from 38 years to 33 years for MPI.Furthermore, NPVs increased by approximately $300.When C was reduced to $200, a similar pattern emerged, with rotations periods shortening by two to four years depending on the GCM and increased NPVs.Conversely, when C was increased to $1000, the rotation periods lengthened by 7 to 8 years and NPVs declined making the overall investment questionable for all three GCMs.This has been shown in other studies as well, where the economic attractiveness of plantations increases in areas where opportunity costs (e.g., land values) are low [9].Finally, when P was adjusted from $20 to $50, the result was a shorter rotation period by 2 to 5 years depending on the GCM.Increasing P turns NPV from negative to positive, making the investment more attractive.Although the increases in NPV vary by region, the highest values occur primarily in the southern portion of the study area, with the lowest values in northeastern Ontario.While perhaps counterintuitive, these shortened rotation periods are due to the increased value of future rotations, hence causing the rotation ages to decline to capture these increases in value.
Table 4 illustrates the results of the sensitivity analyses of the wood + carbon scenario, where r, C, P and P cs were adjusted to observe the effects of parameter uncertainty on the resulting NPVs.Decreasing r from 4% to 2% resulted in longer rotation periods by 3 to 7 years depending on the GCM; NPV also increased by as much as $1,000 in the southern region of the study area.When r was increased from 4% to 8%, shorter rotation periods were identified for all three GCMs, and NPVs were reduced.When C was adjusted, the results were similar to those described above for the wood-only scenario: as C increases, the rotation periods become longer and NPVs decrease.When P was adjusted from $20 to $50, the rotation periods shortened by 4 to 9 years depending on the GCM; however, NPVs increased by almost 3 fold in southern Ontario, making the investment significantly more attractive.The results highlight the importance of future price and yield expectations and establishment costs in evaluating forest investments.Forest managers will have to make judgements about their ability to adjust costs and still attain desirable future yields.In this study, we have not adjusted yields based on establishment costs given the lack of empirical information on this issue.This will be a subject of future research.
The optimal harvest age and related NPVs increased with increasing carbon price.In comparison to the wood-only results, the positive NPVs associated with the range of P cs values, indicate the significance of including C uptake benefits in this investment.These results also provide some evidence for the potential cost-effectiveness of black spruce regeneration efforts, as these values could be compared to other possible carbon sequestration activities.The results from the sensitivity analyses also suggest that the discount rate is a critical factor in determining the optimal harvest age and value.To properly compare carbon sequestration costs, the discount rate should be consistently applied across the options being examined.Because northern temperate forest regeneration efforts require long-term investments, perceptions of the time value of money are very important [61,62].

Conclusions and Summary
Our findings suggest that investments in black spruce regeneration for timber production in Ontario are relatively unattractive without introducing financial benefits related to carbon sequestration.This underlines the importance of appropriate carbon pricing policies for making slow-growing tree species a more attractive economically renewable resource [63].The current work contributes to the ongoing effort to identify the effectiveness of particular carbon sequestration prices and projects in forest management.
Our results are portrayed spatially, which is often overlooked in economic analyses.Non-spatial analyses often apply average values over large areas, thus ignoring significant geographic variation of key biological and financial factors.This can be problematic because climate change is likely to affect ecosystems at multiple spatial scales-from local to regional to global.Our use of GCM-derived NPP values and climate-driven yield equations allowed us to present results spatially-albeit at a relatively low resolution-and thus we were able to distinguish regions where black spruce regeneration investment should be most attractive.
A significant aspect of this study is that it incorporates inter-annual future climate variability into the growth and yield model, which in principle allows for more realistic growth estimates.Past studies have utilized generalized forest growth and yield curves where seasonal and inter-annual climate variability are ignored [6][7][8].Interestingly, our findings suggest that the range of black spruce growth and yield estimates for the end of the current century using annual future climate inputs do not vary substantially from those obtained using standard curves [64].We do note, however, the issue of predicting future forest productivity under a changing climate is complex and requires further research, and is the subject of significant ongoing research.
Further research is required to quantify additional risk factors associated with investments in forest regeneration.In our study, risk factors such as wildfire, disease outbreaks, or drought were not explicitly included in the economic calculus (which was done at the stand level).However, while the inclusion of such factors would help to improve future estimates of the economic values associated with forest regeneration, the general result of shortening rotations with increased risk is well known.Martell [65] reported that, when probabilistic fire occurrence was incorporated into a stochastic forest stand rotation model, it led to shorter rotation intervals.Daigneault et al. [66] also reported a reduction in rotation length with increasing fire probability, but noted that this reduction could be offset by the introduction of a carbon pricing system.Inclusion of climate change-modified forest fire risk would be a useful addition, because fire represents a major disturbance in the Canadian boreal forest that could impact large-scale patterns in biodiversity, carbon, vegetation, and forest management strategies [66,67].
There are many uncertainties involved in this work.Our sensitivity analyses have attempted to explore some of the implications of this uncertainty as it relates to the underlying economic parameters; however, there is likely to remain significant uncertainty regarding the future growth of forests in general-with black spruce being a particularly important boreal species.Our use of GCM-based NPP-modifiers in combination with empirical yield equations represents a novel approach for incorporating climate change into future estimates of forest productivity.Ongoing efforts should continue to incorporate new insights into modelling future growth and yield and related economic values as they become available.

Figure 1 .
Figure 1.Location of study are in North America and the Province of Ontario in Canada.Grid sizes are based on the resolution of global climate model outputs for (a) CanESM2 and MIROC (~2.8°) and (b) MPI-ESM (~2.5°).

Figure 1 .
Figure 1.Location of study are in North America and the Province of Ontario in Canada.Grid sizes are based on the resolution of global climate model outputs for (a) CanESM2 and MIROC (~2.8 • ) and (b) MPI-ESM (~2.5 • ).

ForestsFigure 2 .
Figure 2. Time series of simulated temperature, precipitation and Net Primary Production (NPP) from three global climate models; (a-c) CanESM2 (Canadian); (d-f) MIROC (Japanese); and (g-i) MPI (European) for IPCC future climate scenario RCP 8.5 for each grid (~10 km × 10 km) over the study area.Each line represents a grid-specific value from 2006 to 2100.Fitted average lines are also shown.

Figure 2 .
Figure 2. Time series of simulated temperature, precipitation and Net Primary Production (NPP) from three global climate models; (a-c) CanESM2 (Canadian); (d-f) MIROC (Japanese); and (g-i) MPI (European) for IPCC future climate scenario RCP 8.5 for each grid (~10 km × 10 km) over the study area.Each line represents a grid-specific value from 2006 to 2100.Fitted average lines are also shown.

Figure 3 .
Figure 3. Spatial distribution of projections of annual Net Primary Production (NPP) over time as simulated by the Global Climate Models: (a) CanESM2; (b) MIROC; and (c) MPI.

Figure 3 .
Figure 3. Spatial distribution of projections of annual Net Primary Production (NPP) over time as simulated by the Global Climate Models: (a) CanESM2; (b) MIROC; and (c) MPI.

Forests
the NPP modifiers-particularly toward the end of the century when the modifiers are the largest.Forests 2017, 8, 106 10 of 18 which reduces the impact of the NPP modifiers-particularly toward the end of the century when the modifiers are the largest.

Figure 5 .
Figure 5. Spatial distribution of Net Present Values (NPVs) based on wood value only for the (a) CanESM2; (b) MIROC; and (c) MPI models.

Figure 5 .
Figure 5. Spatial distribution of Net Present Values (NPVs) based on wood value only for the (a) CanESM2; (b) MIROC; and (c) MPI models.

Figure 6 .
Figure 6.Spatial distribution of Net Present Values (NPVs) based on wood and carbon sequestration benefits for (a) CanEMS2; (b) MIROC; and (c) MPI model.

Figure 6 .
Figure 6.Spatial distribution of Net Present Values (NPVs) based on wood and carbon sequestration benefits for (a) CanEMS2; (b) MIROC; and (c) MPI model.

Table 2 .
Net Present Values (NPV, $•ha −1 ) of timber and carbon by region for the baseline scenarios.
* OHA represents optimal harvest year; W represents wood only scenario; W + C represents both wood and carbon value scenario.Forests 2017, 8, 106