Optimal Carbon Storage in Mixed-Species Size-Structured Forests

We extend the study of economically optimal carbon storage to a previously unexplored forest type, mixed-species size-structured stands. The ecological model applied in the study is a transition matrix model with growth functions for boreal Norway spruce (Picea abies (L.) Karst.), birch (Betula pendula Roth and B. pubescens Ehrh.), and other broadleaves. The other broadleaved trees are assumed to have no commercial value. We maximize the sum of timber revenues and the value of carbon storage by optimizing the timing and intensity of thinnings and the potentially infinite rotation age. The optimization problem is solved in its general dynamic form using gradient-based interior point methods and a genetic algorithm. We present results for a mixed stand of Norway spruce and birch, and a mixed stand of Norway spruce, birch, and other broadleaves, and compare these to a pure Norway spruce stand. We show that carbon pricing increases stand volume by postponing harvests and limiting them to larger trees, and changes the optimal species composition by increasing the share of Norway spruce relative to birch. Further, carbon pricing incentivizes maintaining other broadleaves in the stand despite their lack of commercial value, thus increasing tree species diversity. We find that sawlog and total yields increase with carbon price. We show that the higher the number of tree species in a stand, the lower the marginal cost of carbon storage.


Introduction
Forest sinks, currently sequestering approximately 30% of global emissions, are likely to play a crucial role in limiting global warming to 1.5 °C by the end of the century (Pan et al. 2011;Walsh et al. 2017). Per hectare, the value of carbon storage may well exceed that of timber production and other provisioning services (Chiabai et al. 2011). These and all other forest ecosystem services rely on forest health, which is under significant pressure due to climate change, biodiversity loss, and the increased transmission of pests (Trumbore et al. 2015). An emerging body of ecological literature highlights the importance of heterogeneous age, size, and species structures for maintaining forest resilience under various threats (Gauthier et al. 2015;Jactel et al. 2017). Forests with more than one tree species are found throughout the world and account for 68% of Europe's forested area (San-Miguel-Ayanz et al. 2015). Due to climate change, management that supports structurally diverse stands may become necessary even in areas currently dominated by single-species even-aged stands (Dymond et al. 2014). However, economic research on carbon storage in structurally diverse stands is very scarce. This study presents a dynamic bioeconomic model, where the combined production of timber and carbon storage in boreal mixed-species size-structured stands is optimized using an empirical growth model and a detailed economic setup.
Changing stand-level forest management practices to enhance carbon storage may be a fast and cost-efficient mitigation option in areas with high forest cover, e.g. in the boreal region (Pihlainen et al. 2014). The stand-level economics of carbon storage has been studied since Plantinga and Birdsey (1994) and van Kooten et al. (1995), who show that carbon pricing tends to moderately increase the rotation period (i.e. time period from regeneration to clearcut). The model used by van Kooten et al. (1995) and extended in numerous later studies assumes a homogenous forest stand that is artificially regenerated and may only be harvested by clearcutting-i.e. the generic version of the Faustmann (1849) optimal rotation model as presented by Samuelson (1976). Using a similar model, Akao (2011) shows analytically how the effects of carbon storage on optimal rotation age depend on assumptions concerning the carbon release from wood products. While these and similar models offer valuable insight on optimal carbon storage in even-aged plantation-type forests, they are unequipped to address important questions related to more heterogeneous forests.
Certain recent even-aged single-species models with carbon storage have been extended to include thinnings, i.e. partial harvesting (Pohjola and Valsta 2007;Daigneault et al. 2010;Niinimäki et al. 2013;Pihlainen et al. 2014). The two latter studies, using very detailed process-based growth models, show that carbon pricing changes (in addition to the rotation age) the optimal internal structure of the stand along the rotation period. However, as the models used in these studies do not include natural regeneration, i.e. the emergence of new saplings to the stand without planting, they cannot be used to analyse uneven-aged forests.
Uneven-aged forestry, also called continuous cover forestry, is an alternative to the clearcutting-based rotation regime. In uneven-aged forestry, the stand is managed by thinnings only, and natural regeneration leads to a heterogeneous age and size distribution. Compared to rotation forestry, management that continuously maintains forest cover is likely to support more ecosystem services (Calladine et al. 2015;Peura et al. 2018) and to improve resilience against climate change (Thompson et al. 2009;Dymond et al. 2014). The first attempts at optimizing uneven-aged management include de Liocourt (1898) and Adams and Ek (1974). As optimizing uneven-aged forestry involves dynamic complexities and related computational challenges, many studies have resorted to simplified static models and models without sound economic bases (see discussion in Haight (1989), p. 287-295 and).
Related to the question of age-and size structure is the question of species composition. A body of forest ecological and silvicultural research exists on tree species interactions and their effects on growth and yield (Rothe and Binkley 2001;Forrester 2014). Recent research has highlighted the importance of mixed-species forests for resilience (Jactel et al. 2017;Brandl et al. 2020) and biodiversity (Felton et al. 2010;Cavard et al. 2011). However, economically oriented research on mixed-species stands remains rather scarce. Knoke et al. (2005) study even-aged mixtures of Norway spruce (Picea abies (L.) Karst.) and European beech (Fagus sylvatica L.) using Monte Carlo simulation, showing that while mixed-species management reduces profitability compared to single-species management, it enables portfolio diversification and thus risk attenuation.
The complex problem of managing mixed-species uneven-aged stands has been studied in a dynamic setting by Haight and Getz (1987), Haight and Monserud (1990a, b), and Rämö and Tahvonen (2015). Tahvonen et al. (2019) study timber production and ecosystem services in boreal mixed-species stands using a generalized model that includes both continuous cover forestry and rotation forestry and optimizes the choice between them. Parkatti and Tahvonen (2020) explore various boreal species mixtures and management objectives, demonstrating that overyielding (i.e. higher productivity of mixed stands compared to monocultures) in terms of physical output does not reveal the economically preferable species combination. Knoke et al. (2020) show that multi-criteria optimization integrating uncertainty and biological interactions supports the application of both mixed-species management and continuous cover forestry.
While recent years have seen a marked increase in interest in both uneven-aged and mixed-species forests, economic studies on carbon storage in such non-homogeneous stands are very few, and most of the existing contributions have limited their scope to static settings (e.g. Trasobares and Pukkala 2004;Buongiorno et al. 2012). Exceptions include Boscolo and Vincent (2003), studying production non-convexities in tropical rainforests, and Goetz et al. (2010), applying a dynamic model of uneven-aged stands of Scots pine (Pinus sylvestris L.) in Spain. A model for studying optimal carbon storage with endogenous management regimes-continuous cover forestry or rotation forestry-is presented and applied to boreal Norway spruce in  and . The results of these two studies suggest that carbon pricing changes thinning schedules, increases the rotation age, and eventually causes a regime shift from rotation forestry to continuous cover forestry. The present study extends that approach to boreal mixed-species stands, and as far as we know, is the first one to study carbon storage in mixed-species size-structured stands without a priori assuming the management regime.
Using a statistical-empirical size-structured growth model and including fixed and variable harvesting costs, we determine the economically efficient methods for increasing carbon storage in boreal mixed stands. We explore the cases of a mixed stand of Norway spruce and birch, and a mixed stand consisting of Norway spruce, birch (silver birch Betula pendula Roth and downy birch B. pubescens Ehrh.), and commercially non-valuable other broadleaves, and compare these to a pure Norway spruce stand. We show that carbon pricing changes optimal harvests by postponing thinnings and by increasing the size of harvested trees and trees left standing. While carbon pricing increases overall stand density, it also changes the species composition by increasing the share of Norway spruce relative to birch, and by incentivizing maintaining other broadleaves in the stand despite their lack of commercial value. Additionally, we show that moderate carbon price levels increase the steady-state yields of both Norway spruce and birch. Finally, we show that the marginal cost of carbon storage is the lower, the higher the tree species diversity.
We continue by introducing the growth model and the optimization problem. Thereafter we present the empirical parameter values and the computational methods. This is followed by results on optimal stand management, on timber production and carbon storage, and on forestry revenues and carbon storage costs. Finally, we discuss our results in light of earlier studies and draw conclusions.

3 2 The Growth Model and the Optimization Problem
We denote the number of trees of species i in size class s at the beginning of period t by x ist , i = 1, 2, ..., m, s = 1, 2, ..., n, t = t 1 , t 1 + 1, ..., T . Accordingly, the stand state at the beginning of period t can be given as Let us denote the fraction of trees of species i moving to size class s + 1 in period t by is t , i = 1, 2, ..., m, s = 1, 2, ..., n , where in t ≡ 0. The natural mortality of trees of species i in size class s in period t is is t , i = 1, 2, ..., m, s = 1, 2, ..., n. Thus the fraction of trees remaining in the same size class equals 1 − is t − is t , i = 1, 2, ..., m, s = 1, 2, ..., n . Natural regeneration is described by ingrowth, i.e. trees entering the smallest size class. Tree ingrowth of species i at the beginning of period t is denoted by i t , i = 1, 2, ..., m . Additionally, we denote the number of trees harvested from species i and size class s at the end of period t by h ist , and trees felled but left in the forest by f ist , i = 1, 2, ..., m, s = 1, 2, ..., n, t = t 1 , t 1 + 1, ..., T . Hence, stand development may be described by the difference equations where i = 1, 2, ..., m, s = 1, 2, ..., n, t = t 1 , t 1 + 1, ..., T.
We assume that the stand is artificially regenerated immediately after a clearcut at t = 0 . At the beginning of period t = t 1 , we have an initial stand state t 1 composed of trees that have reached the smallest size class. In addition to trees grown from planted saplings, trees of various species may have naturally emerged in the stand. The length of the rotation period, i.e. the time between clearcuts, is denoted by T ∈ t 1 , ∞ .
Let w ≥ 0 (€ ha −1 ) denote the cost of artificial regeneration. We denote the discount factor by b = 1∕(1 + r) , where r refers to the annual interest rate. The length (in years) of a period is denoted by Δ . Denoting the timber product assortments of sawlog and pulpwood by k = 1, 2 , respectively, the timber assortment volumes of species i in a tree of size class s are given by v isk . The total stem volume of a tree of species i in size class s is v is = v is1 + v is2 . The (roadside) prices (€ m −3 ) for sawlog and pulpwood for species i are p i1 and p i2 . Gross harvesting revenues per period, R t , depend on the number, size and species of trees harvested, and are given by We take into account that certain tree species may have no commercial value, i.e. for some l ∈ {1, 2, ..., m} , p l1 = p l2 = 0 . This implies that when such trees emerge in the stand through natural regeneration, it may be optimal to prevent them from competing with valuable trees by felling them but leaving them in the stand. In this case, the trees are neither pre-processed nor hauled out of the forest, implying that the variable costs of felling consist of only cutting costs. The sum of variable harvesting and felling costs in period t is defined by C u t , t , u = th, cl , where u stands for thinnings and clearcuts, respectively. This takes into account the higher time consumption in thinnings compared to clearcuts.
Additionally, a fixed cost C incurs whenever harvesting or felling takes place in the stand. The fixed cost C covers e.g. planning and the transportation of machinery to the stand site. Because of the fixed cost, it may not be optimal to harvest the stand during every period. This is taken into account by the binary variables t ∈ {0, 1}, t = t 1 , t 1 + 1, ..., T and by the Boolean operators h ist = t h ist and f ist = t f ist . When the choice is t = 1 , the levels of h ist ≥ 0, f ist ≥ 0, i = 1, 2, ..., m, s = 1, 2, ..., n can be freely optimized. When t = 0 , the only admissible choice is h ist = f ist = 0, i = 1, 2, ..., m, s = 1, 2, ..., n.
We assume that society pays the forest owner for the carbon dioxide (CO 2 ) that is absorbed as the stand grows and charges for the CO 2 released as a consequence of harvesting and the decomposition of deadwood. When pricing carbon storage, we consider carbon storage in living trees and deadwood but exclude harvested wood products. 1 We x ist v is denote the total stem volume of the stand at the beginning of period t. Deadwood is created both through natural mortality and from trees felled but left in the forest. The deadwood volume formed in period t through natural mortality equals We denote the annual decay rate of deadwood by g . Per unit of deadwood, the present value of future emissions due to decay equals p c (r) , where (cf. . We let denote the quantity of CO 2 in one wood volume unit, and let p c ≥ 0 (€ tCO 2 −1 ) denote the economic value of CO 2 . We assume that the carbon price is constant, reflecting constant damages from carbon emissions over time (Richards and Stokes 2004). Thus the economic value of net carbon sequestration (or net negative emissions) in period t can be given as for t = 0, 1, ..., T , where t+1 − t denotes net volume growth in period t and the term [1 − (r)] d ,t + d f ,t takes into account carbon storage in and release from deadwood.
Following Faustmann (1849), the infinite series of identical rotations can be written in compact form, and the problem of optimizing stand management over an infinite horizon can be given as subject to (1)-(3), (4), (6) and The optimized variables are the number of harvested and felled trees per size class and species in thinnings, the timing of thinnings, and the rotation age T. The choice of T determines the optimal management regime: clearcutting is optimal if the objective functional is maximized by a finite rotation age. If no maximum with finite T exists and the net present value converges toward the continuous cover forestry net present value from below as T → ∞ , continuing the thinnings while maintaining forest cover is optimal.

Empirical Specifications and Data
We apply an empirical transition matrix growth model by Bollandsås et al. (2008) at latitude 61.9 ºN. The model has been estimated using the National Forest Inventory of Norway and includes functions for ingrowth, mortality, and diameter increment in single-and mixed-species stands with both intra-and interspecies competition. The model includes species-specific growth parameters for Norway spruce and group-specific growth parameters for birch (silver and downy birch) and other broadleaves (e.g. oaks [Quercus sp.], maples [Acer sp.], European beech and Eurasian aspen [Populus tremula L.]). We study an average productivity site ( SI = 15 ), implying that the height of the dominant trees at the age of 40 (100) years is 15 (24) metres. We use 12 size classes with diameters (midpoints) ranging from 7.5 to 62.5 cm in 5.0-cm intervals. The length of a period ( Δ ) is five years.
The empirical setting approximates a typical Nordic situation where a cohort of an economically preferred native coniferous tree species is artificially established on bare forest land, followed by the natural regeneration of not only this species but possibly also certain economically less-valuable species such as broadleaf trees. The broadleaves may enter the stand through air-borne seeding, originating from forests in the vicinity. On all studied stand types, Norway spruce is planted at the beginning of the time horizon, and 1750 Norway spruce trees emerge in the smallest size class 20 years after artificial regeneration (i.e. t 1 = 4 ). We study three different cases: the first case, included for reference, is a pure Norway spruce stand. In the second case, the initial stand structure (i.e. the stand 20 years after regeneration) consists of 1750 small spruce trees and 1000 birch trees (silver and downy birch) that have entered the smallest size class. In the third case, in addition to the 1750 spruce trees and 1000 birch trees, 500 other broadleaf trees have naturally emerged to the smallest size classes by time t 1 . Such variation in initial stand composition may result from either varying tree species composition on neighbouring forest sites or from seedlings of broadleaved trees being removed by the forest manager (as is currently typical in e.g. Finland).
Following Bollandsås et al. (2008), the estimated number of trees of species i entering the smallest size class (i.e. natural regeneration) during the five-year period t is given as.
where B t is the total stand basal area (m 2 ha −1 ) at the beginning of period t, PB it the proportional basal area of species i (expressed as percentage of total basal area) at the beginning of period t, and S the H 40 site index.
The fraction of trees of species i that move to size class s + 1 during period t is given as where M s is the diameter (midpoint) of size class s, A st is the total basal area (m 2 ha −1 ) of size classes s + 1, ..., n at the beginning of period t, L is the latitude, and W is the size class interval (5.0 cm). Finally, the estimated natural mortality (i.e. fraction of trees that die) of tree species i in size class s during the 5-year period t is given as Coefficient values for the growth functions given above are presented in the Appendix, Table 4. Table 5 in the Appendix presents the species and size class -specific sawlog and pulpwood volumes, v is1 and v is2 (Rämö and Tahvonen 2015), assuming that other broadleaves have similar volumes as birch. The roadside prices for Norway spruce sawlog and pulpwood are €58.44 m −3 and €34.07 m −3 , respectively. The corresponding prices for birch are €49.73 m −3 and €30.50 m −3 . Other broadleaves emerge naturally onto many Fennoscandian forest sites, but based on the current market situation in Finland, we assume they have no commercial value.
The specification of harvesting cost functions follows from empirical logging experiments by Nurminen et al. (2006), based on the performance of modern harvesters and specified separately for thinnings and clearcuts. The per-period variable harvesting and felling costs are given as where v is is the total stem volume of species i in size class s, and C iuz ; i = 1, ..., m, u = th, cl, z = 0, ..., 8 are parameters ( Table 6 in the Appendix). The variable costs include the cost of cutting (16a), hauling (16b), and non-commercial felling (16c). The per-tree cost of felling is lower than that of commercial harvesting, as no preprocessing or hauling is required. Further, cutting is more expensive in thinnings than in clearcuts, which is taken into account by coefficient C iu1 . The fixed harvesting cost C equals €500 ha −1 . The cost of artificial regeneration carried out at the beginning of the time horizon is assumed to equal €1000 ha −1 (Niinimäki et al. 2012). We assume a 3% annual interest rate.
Based on Lehtonen et al. (2004), the stemwood density factor is 0.38 tonnes of dry matter per cubic metre of fresh volume. Following Niinimäki et al. (2013), the CO 2 content of a wood volume unit ( ) is obtained by multiplying the density factor with the share of carbon in biomass dry weight (0.5) and the coefficient that converts tonnes of carbon to tonnes of CO 2 (44/12). Thus, we set = 0.697 tCO 2 m −3 . For the decay rate of deadwood, we use g = 0.055 based on the average of relevant stem sizes in Hyvönen and Ågren (2001).

Computational Methods
Because harvest timing variables are integers, but harvest and felling intensities are continuous, the task is to solve a mixed-integer nonlinear programming problem. To do this, we apply a bi-level optimization method (Colson et al. 2007) previously applied to a forestry setting e.g. in Tahvonen and Rämö (2016). The lower-level problem is computed using version 10.3 of the Knitro optimization software, which applies advanced gradient-based interior point algorithms (Byrd et al. 2006). Given any vector of harvest timing binaries, the maximized objective value of the lower-level problem forms the objective value. The harvest timing vector is optimized using a genetic algorithm (Deb and Sinha 2010;Sinha et al. 2017), and the optimal harvest schedules are solved for a series of rotation lengths. The optimal rotation is finite if the objective function obtains a maximum with some T ∈ [60, 180) years. Examples of such solutions are shown in Fig. 7 in the Appendix. If the value of the objective function continues to increase as the rotation period is lengthened (see Fig. 8 in the Appendix for illustration), the optimal rotation is infinite, and the optimal continuous cover solution is obtained by lengthening the horizon to reach a close approximation of the infinite horizon solution. To handle potential non-convexities, we apply multiple randomly chosen initial points in the optimization.

The Effects of Carbon Pricing on Optimal Stand Management
We present results for three cases: a pure Norway spruce stand, a mixed stand of Norway spruce and birch, and a stand containing Norway spruce, birch, and the commercially non-valuable other broadleaves. The optimal rotation period is infinite, i.e. it is optimal to manage the stand by applying continuous cover thinnings, in all the studied cases given a 3% annual interest rate and carbon prices €0-€50 tCO 2 −1 ( Table 1). The first thinning after stand regeneration is carried out at stand ages 40-55 years depending on carbon price and species combination. The beginning of thinnings is postponed by carbon pricing and hastened by a larger number of tree species naturally regenerating on the site (Table 1). The latter effect is due to the need to control stand density and (when relevant) to enable the eradication of the other broadleaves. The large initial tree cohorts are intensively utilized in a series of thinnings before approaching the steady state ( Figs. 1 and 3).
The steady-state harvest cycles entail a harvesting operation every 20 or 25 years where all trees above a species-specific diameter are harvested and the smaller trees are left standing (i.e. thinning from above). Targeting harvests to the largest trees of each species is a universal feature in our optimal thinning solutions. As relative value growth is high in small trees, it is optimal to postpone harvesting until trees have grown to a sawlog-yielding size. However, the size to which trees are allowed to grow before they are harvested depends on carbon price. Specifically, the higher the carbon price, the larger the steady-state mean size   (Table 1). Thus, the timing of the harvests and the size of the harvested trees adjust to the carbon price to maintain optimal average stand density and economic return along the steady-state cycle.
In the case of the pure Norway spruce stand, zero carbon price implies that trees in diameter classes (midpoints) 22.5-42.5 cm are cut every 25 years, leaving only trees in diameter classes 7.5-17.5 cm standing, whereas given a €50 tCO 2 −1 carbon price trees with diameters 32.5-47.5 cm are cut every 20 years (Table 1). The resulting mean stand volumes over the steady state harvest cycle are 79 and 169 m 3 ha −1 .
In the case of a mixed stand of Norway spruce and birch and given zero carbon price, steady-state harvesting occurs every 25 years and targets spruces in diameter classes 22.5-42.5 cm and birches in diameter classes 17.5-32.5 cm (Table 1, Fig. 2). As shown in Figs. 2 and 4 depicting optimal steady-state structures, the number of trees decreases with tree size class. 2 Increasing the carbon price from €0 to €25 tCO 2 −1 implies that the steadystate harvesting size of both spruce and birch increases by one diameter class. Increasing the carbon price to €50 tCO 2 −1 shifts the steady-state harvest of Norway spruce to diameter classes 32.5-47.5 cm. In contrast, the average size of the harvested birch decreases somewhat compared to management with €25 tCO 2 −1 carbon price. Silver and downy birches are pioneer species with limited shade tolerance, implying that their growth is affected by the high stand density. Thus, when the carbon price is high, the share of spruce relative to birch increases ( Figs. 1 and 2).
In the most complex case, i.e. the stand containing Norway spruce, birch, and the commercially non-valuable other broadleaves, steady-state harvests given zero carbon price target Norway spruce and birch similarly as in the case of a spruce-birch stand. During each harvesting operation, all other broadleaves are felled and left in the forest (Fig. 4), as this is the most inexpensive way to prevent them from competing for resources with commercially valuable trees. With a carbon price of €25 tCO 2 −1 , it is optimal to limit the steady-state harvesting of both spruce and birch to larger size classes and to allow a small number of other broadleaves to grow in the stand, felling those from diameter classes 17.5-32.5 cm (Figs. 3 and 4). With a carbon price of €50 tCO 2 −1 , a high mean stand volume is maintained by harvesting large spruce and birch trees and by allowing a number of other broadleaves to reach a notably large size (42.5-52.5 cm) by the time they are felled. At the steady state, the standing volume of other broadleaves is approximately the same as that of birch (Fig. 3). As other broadleaves can subsist even in a dense stand, it is optimal to increase carbon storage through maintaining some of the other broadleaves in the stand, despite their lack of commercial value.

The Effects of Carbon Pricing on Timber Production and Carbon Storage
Mean annual sawlog yields at the steady state range from 3.6 to 5.7 m 3 ha −1 while mean annual total yield (sum of sawlog and pulpwood) range from 4.6 to 6.1 m 3 ha −1 (Table 2). For all studied species mixtures, both sawlog yield and total yield at the steady state increase with carbon price given the carbon price range €0-€50 tCO 2 −1 . The sawlog-total ratio of the yield, ranging from 0.78 to 0.93, is increasing in carbon price because the sawlog-pulp ratio is higher in larger size classes (see Table 5 in the Appendix).

3
The pure Norway spruce stand has the highest steady-state yields followed by the spruce-birch mixture stand. In the latter stand type, spruce yield increases from 3.9 to 5.5 m 3 ha −1 as the carbon price increases from zero to €50 tCO 2 −1 , while birch yield (0.8 m 3 ha −1 ) is higher with a carbon price of €25 tCO 2 −1 than with either zero carbon price or €50 tCO 2 −1 carbon price. This difference is explained by the high shade tolerance of spruce relative to birch. Additionally, to optimally utilize both species simultaneously, the growing volume of Norway spruce must be kept slightly lower than in the pure spruce stand.
Given zero carbon price, the mixed stand containing Norway spruce, birch, and other broadleaves produces similar yields as the mixed stand without other broadleaves (Table 2).

Fig. 2
Optimal steady-state structures in stands of Norway spruce and birch, with carbon prices €0, €25, and €50 tCO 2 −1 . Size classes begin from a diameter of 7.5 cm and increase in 5-cm intervals. Note: 3% interest rate This is because the other broadleaves are felled at the earliest convenience. However, when carbon price increases to €25 tCO 2 −1 , certain other broadleaves are allowed to grow to a medium size. The competition caused by the small number of other broadleaves is not sufficient to notably lower the spruce or birch yields compared to the spruce-birch stand. Given a carbon price of €50 tCO 2 −1 , carbon storage is increased not only by increasing the growing volume of spruce but also by allowing the other broadleaves to take up more space in the stand. Hence the addition of the other broadleaves lowers optimal yields somewhat compared to the spruce-birch stand. The development of carbon stocks over time shows that each harvest decreases the quantity of carbon stored in standing trees (Fig. 5). In addition to living trees, carbon is stored in deadwood. In the stand consisting of Norway spruce and birch (Fig. 5a), the relatively high stand density caused by the large initial tree cohorts results in increased   natural mortality during stand ages of approximately 40-80 years. Thus, the carbon stock in deadwood first increases and then decreases towards the low steady-state level. Due to higher stand density, this effect is stronger with carbon pricing. In the case of a stand containing Norway spruce, birch, and other broadleaves (Fig. 5b), zero carbon price implies that the carbon stock in deadwood increases sharply after the first harvesting operation, during which the initial cohort of other broadleaves is felled. However, given a €50 tCO 2 −1 carbon price, the carbon stock in deadwood remains relatively large throughout time, because while other broadleaves are still controlled by fellings, some are first allowed to reach a large diameter, producing a considerable pulse of deadwood after each harvesting and felling operation.
Steady-state mean carbon storage in standing trees and in deadwood increase with carbon price irrespective of species composition (Table 2). In the case of the Norway spruce and birch stand, the mean carbon storage in standing trees (deadwood) equals 55.2 tCO 2 ha −1 (3.0 tCO 2 ha −1 ) with zero carbon price, and increasing the carbon price to €50 tCO 2 −1 increases mean storage to 114.9 tCO 2 ha −1 (5.7 tCO 2 ha −1 ). In the case of the stand also containing other broadleaves, mean carbon stocks in both standing trees and deadwood are slightly higher, especially with a carbon price of €50 tCO 2 −1 . This is due to the other broadleaves that are maintained in the stand exclusively for the purpose of carbon storage. The steady state is typically reached as late as approximately 200-250 years after stand regeneration. This implies that the (value of) carbon storage occurring during the long transition phase towards the steady state is very important for the economic outcome. Discounted CO 2 sequestration (tCO 2 ha −1 ) is the sum of all periodic net carbon fluxes within the infinite planning horizon, each discounted to the present (stand regeneration) moment. For example, given the mixed stand containing Norway spruce and birch, and a carbon price of €25 tCO 2 −1 , the net carbon sequestration over the infinite horizon is equivalent to 64.9 tonnes of CO 2 emissions abated immediately (Table 2). Discounted CO 2 sequestration increases monotonously with carbon price, and is the higher, the greater the number of tree species naturally regenerating in the stand.

Forestry Revenues and the Cost of Carbon Storage
Discounted timber income decreases with carbon price irrespective of species composition (Table 3), despite the fact that steady-state timber yields increase with carbon price (Table 2). This is partly because harvests are carried out later when carbon is valued, which affects the net present values due to discounting. Another reason is that diverging from the optimal timber-only solution may imply higher harvesting costs per timber unit. Given any carbon price, the highest discounted timber income is produced by the mixture stand of Norway spruce and birch. While the roadside prices and commercial volumes of birch are somewhat lower than those of Norway spruce, the naturally emerging birch is a valuable addition to the stand as it is possible to manage them over time in a manner that yields higher net benefits than a pure Norway spruce stand. However, the stand consisting of other broadleaves in addition to spruce and birch yields lower discounted timber incomes than either the pure Norway spruce stand or the spruce-birch mixture, because controlling the other broadleaves by fellings incurs costs.
The economic value of carbon storage is considered in its full extent in our model (i.e. we do not consider only carbon storage that is additional compared to baseline management) and is represented by carbon subsidies. Net present values improve considerably when carbon storage is priced, because the earned carbon subsidies more than compensate the decrease in timber income caused by adapting management (Table 3). In the spruce-birch mixture case, a carbon price of €25 tCO 2 −1 increases net present value by 48% compared to zero carbon price. The corresponding figure for the stand with Norway spruce, birch, and other broadleaves is 56%. Thus, with higher (but still plausible) carbon prices, the economic value of carbon storage may overweigh the economic value of timber production.
The mixed stand with other broadleaves has the lowest net present value based on timber income only. However, given a €25 tCO 2 −1 carbon price, this stand type fares better than the pure Norway spruce stand; given a €50 tCO 2 −1 carbon price, it also outperforms the spruce-birch mixture. This implies that taking carbon benefits into account may improve the economic attractiveness of diverse species mixtures in forestry.
The economic cost of carbon abatement in forestry, i.e. the cost of increasing carbon storage, can be measured by foregone timber income. To obtain marginal abatement costs, we divide the incremental decrease (that is, the decrease compared to the solution given a lower carbon price) in timber income for each optimal solution by the corresponding incremental increase in discounted CO 2 sequestration. Figure 6 presents marginal costs in Norway spruce stands, in mixed stands of Norway spruce and birch, and in stands consisting of Norway spruce, birch, and other broadleaves. Regardless of species composition, marginal abatement costs increase with the quantity of carbon abatement. In the pure Norway spruce stand, marginal costs range from €6 to €35 tCO 2 −1 for 5-19 tonnes of carbon abatement per hectare, while the mixed stand of Norway spruce and birch has slightly lower marginal costs. Increasing carbon storage is least expensive in the stand consisting of Norway spruce, birch, and other broadleaves.

Discussion and Conclusions
The stand-level economics of carbon storage has heavily focused on single-species evenaged forestry. However, studies analysing carbon storage with a single-species model that optimizes the choice between an even-aged rotation regime and an uneven-aged continuous Fig. 6 Marginal abatement costs in stands of Norway spruce, Norway spruce and birch, and Norway spruce, birch, and other broadleaves. Note: 3% interest rate cover regime show that carbon pricing often causes a regime shift: instead of clearcutting, the stand is utilized by continuing thinnings indefinitely ). The present study applies a similar generalized model framework to mixed-species stands. With a 3% interest rate and a regeneration cost of €1000 (in the low range, see Niinimäki et al. 2012) continuous cover forestry is superior to rotation forestry regardless of carbon pricing. As shown in Fig. 7 in the Appendix, with an interest rate of 1% rotation forestry is the economically dominating management regime given the studied range of carbon prices (€0-€50 tCO 2 −1 ). Besides a lower interest rate, factors supporting the optimality of rotation forestry within the studied setting include low regeneration costs, high site productivity, and high number of viable young spruce trees resulting from the artificial regeneration efforts (see Tahvonen and Rämö 2016 on the effects of these parameters).
The results of our present study show that carbon pricing changes thinnings by postponing them and by increasing the average size of harvested trees and trees left standing, i.e. it is optimal to allow trees to grow larger before they are harvested. This is in line with . However, we find that carbon pricing has a non-uniform effect on various tree species depending on their ecological characteristics. Carbon pricing monotonously increases the optimal standing volume of Norway spruce, a shade-tolerant species, but the standing volume of silver and downy birches, relatively shade-intolerant pioneer species, are higher under a moderate carbon price than under either zero or high carbon price.
Our study examines cases where both commercially valuable and non-valuable tree species exist in the stand. Such a situation in a boreal context was first studied by Rämö and Tahvonen (2015), who show that if the non-valuable other broadleaves are left undisturbed (e.g. for aesthetic reasons), they begin to dominate the stand rather quickly, impairing the growth of valuable tree species and eroding the economic profitability of forestry. Our present study, similarly to Tahvonen et al. (2019), applies empirically estimated cost functions for harvesting and for merely felling without pre-processing or hauling. Tahvonen et al. (2019) find that if biodiversity benefits are not taken into account, it is optimal to fell all the other broadleaves to prevent them from competing with commercially valuable species. Our results with zero carbon price confirm this finding. While the cost of felling is lower than that of harvesting, controlling the commercially worthless trees is still costly. This underlines that continuous cover forestry in mixed-species size-structured stands does not constitute 'high grading' when management is optimized over the infinite time horizon assuming well-defined property rights.
The results of our present study show that while carbon pricing slightly lowers the optimal share of the shade-intolerant birch relative to spruce, it also incentivizes maintaining some of the commercially non-valuable other broadleaves in the stand, as they contribute to the carbon stocks both while standing and when eventually felled and decaying. Similar effects are found by Tahvonen et al. (2019) when utility is derived from ecosystem services. Additionally, while the other broadleaves are treated as one "species" in the growth model, in reality they may consist of several tree species, and hence their presence may considerably increase the number of tree species in the stand. The increase in tree species diversity as a side product of the internalization of carbon storage benefits is a completely new result in the economic literature on carbon sequestration in forests. Moreover, we find that carbon pricing increases the number of large trees and the deadwood quantity. As such features, along with tree species diversity, support forest biodiversity (McElhinny et al. 2005;Gao et al. 2015), our results suggest that incentives for increasing carbon storage may yield considerable ecological side benefits.
The 3% interest rate assumed in the present study is in the mid-range of interest rates typically applied in stand-level economic studies in the Fennoscandian context. The role of the interest rate in mixed-species forestry has been explored in Rämö and Tahvonen (2015). The authors show that in boreal uneven-aged stands, a higher interest rate promotes species diversity by favouring lower stand volume, which in turn allows birch to grow at an economically sufficient rate. Moreover, the effect of carbon pricing on stand management is more pronounced, the higher the interest rate . This follows from the greater incentive to shift net emissions forward in time through carbon storage. Hence, it is reasonable to assume that within our mixed-species optimization framework, interest rates higher than 3% would imply a greater relative share of birch in both standing volume and yield and stronger changes in harvest schedules in reaction to carbon pricing.
As the model presented in our study is a stand-and forest owner-level setup, it does not include market interactions explicitly. According to our results, carbon pricing increases mean total yield and the sawlog share of the total yield regardless of the species composition of the stand. This is because carbon pricing motivates increasing the stand density, and a higher level of growing volume implies higher annual growth as long as stand density is below the growth-maximizing level. The latter is clearly the case given 3% interest rate and zero carbon price. However, very high levels of carbon price would likely begin to decrease yields because only notably large trees would be harvested from a very dense stand.
According to our results, marginal costs range from below €3 to €34 tCO 2 −1 for 5-25 tonnes of carbon abatement per hectare in the stand consisting of Norway spruce, birch, and other broadleaves, and is more expensive in stands of Norway spruce and birch. Highest marginal costs are found for the pure Norway spruce stand. The potential carbon abatement that can be achieved within a similar cost range is somewhat higher in Niinimäki et al. (2013) and , both studying boreal Norway spruce stands. This is likely due to their considering a wider range of carbon pools (i.e. living and decaying nonmerchantable tree compartments such as branches). Compared to carbon storage cost estimations presented for the boreal region in the meta-regression analysis by van Kooten et al.
The result that marginal costs for carbon abatement are lower, the higher the number of tree species in the stand, is a new finding within the economic literature. However, it is in line with the ecological literature stating that tree species diversity supports carbon storage (e.g. Ruiz-Benito et al. 2014;Mensah et al. 2016). These findings are explained by a higher productivity of diverse ecosystems due to more efficient resource utilization, and related to the concept of overyielding in mixed-species stands (see e.g. Lu et al. 2016).
However, Parkatti and Tahvonen (2020) show that the relative productivities of various species mixtures are highly dependent on stand management. This is especially true when management is optimized to maximize economic profit instead of physical output, and hence timber prices and the costs of both harvesting commercial tree species and controlling non-valuable ones play a large role. In the case of the present study, we have tree species with varying commercial value (Norway spruce is somewhat more valuable than birch, and other broadleaves are non-valuable) and shade-tolerance (spruce tolerates shading very well, other broadleaves somewhat, and birch poorly). From an economic viewpoint, it seems that these differing features of various tree species may play a role in combining timber production and carbon storage in a cost-efficient manner. Hence, our most species-diverse stand type that also has the lowest net present value given zero carbon price becomes the one with the best economic performance when carbon storage benefits are internalized. This implies that the currently common practice of eradicating lesser-valued broadleaves during the tending of the seedlings as well as in thinnings may be suboptimal when carbon storage is taken into account.
In our model, the carbon storage externality is valued to its full extent. In practice, policy makers wishing to curb public expenditure are likely to prefer a subsidy system based on carbon storage that is additional to the business-as-usual management solution. At the stand level, both approaches create the same incentives for adapting management (Tahvonen and Rautiainen 2017). The present study omits carbon storage in harvested wood products, which is the current approach in the New Zealand forest carbon scheme (Adams and Turner 2012). It is interesting to note that our results show no qualitative discrepancies to those obtained in , where carbon storage in (and release from) sawlog and pulpwood products is included in the model. Adding this feature to the present model would slightly moderate the effects of carbon pricing, as the carbon release caused by harvesting would be slower.
The approach of the present study is fully deterministic, and hence uncertainties related to prices and environmental factors are omitted. Stochastic phenomena with high relevance for timber production and carbon storage include fire hazard and pest damage, and as both of these are linked to the structure and quantity of living and dead trees in the stand, a model extension towards that direction would be interesting. These questions are recently studied in Malo et al. (2021) but without carbon pricing. Also beyond the scope of this study are questions related to transaction costs e.g. from carbon measurement, reporting, and auditing (see e.g. Köhl et al. 2020 in the context of REDD + and Kelly et al. 2017 regarding a voluntary-based forest offset market within the Californian cap-and-trade protocol). Finally, we feel that issues related to policy design and implementation, as well as the role of forest management within the wider market-level context of climate policy (issues such as substitution), are among top-priority topics for future research.
In this paper, we have extended the study of economically optimal carbon storage to mixed-species size-structured stands. We show that carbon pricing postpones thinnings and limits them to larger tree diameter classes. As high stand densities affect Norway spruce less than birch, carbon pricing also changes the optimal species composition by increasing the volume share of the former. If saplings of commercially non-valuable other broadleaves emerge in the stand, they are allowed to grow larger only if carbon storage has economic value. Finally, we find that the marginal cost of increasing carbon storage is lower, the higher the number of naturally regenerating tree species at the site. As diverse age, size, and species structures within forest stands promote resilience under climate change (Gamfeldt et al. 2013), our results suggest that cost-efficient mitigation and adaptation may go hand in hand.

Fig. 7
Bare land value as a function of rotation age in a stand of Norway spruce and birch, with 1% interest rate and carbon prices €0, €25, and €50 tCO 2 −1 . Optimal rotation ages denoted by circle symbols Fig. 8 Bare land value as a function of rotation age in a stand of Norway spruce and birch, with 3% interest rate and carbon prices €25 and €50 tCO 2 −1 . The bare land value of the rotation solutions approach the bare land value of the continuous cover solutions from below as the rotation age approaches infinity  Table 6 Parameter values for the harvesting cost function (Nurminen et al. 2006) th denotes thinning, cl denotes clearcuts. See Eqs. 16a-16c for symbols