Pro ﬁ tability of continuous-cover forestry in Norway spruce dominated peatland forest and the role of water table

: Continuous-cover forestry (CCF) is expected to reduce the negative environmental impacts of peatland forestry in comparison with rotation forestry (RF), but the unknown pro ﬁ tability of CCF on peatlands limits its application in practice. The pro ﬁ tability of CCF was analyzed by simulating management scenarios with a process-based ecosystem model, EFIMOD, which was complemented to describe the interplay between tree growth and water table depth, which is typical of peatland forests. A variety of harvest intervals and post-harvest basal areas for a mature Norway spruce ( Picea abies (L.) Karst.) dominated stand was simulated on a nutrient-rich peatland site. Conventional RF was simulated for comparison. CCF provided a higher pro ﬁ t than RF. The best ﬁ nancial performance was obtained with a 15-year harvest interval regardless of interest rate, although the overall pro ﬁ tability of CCF depended on the interest rate used. Ditch network maintenance was needed to maintain the stand growth only when the post-harvest basal area was smaller than 10 m 2 ·ha (cid:1) 1 . There were many CCF scenarios in which the difference in the net present value of harvest revenues was within 10% compared with the best CCF scenario. Hence, there are many relatively pro ﬁ table CCF harvesting alternatives for forest management in boreal spruce-dominated peatland forests.


Introduction
Globally, there are about 15 Mha of forests growing on peatland (Joosten and Clarke 2002), as a result of either draining and afforestation (e.g., Sloan et al. 2018) or, as is more common in northern Europe, draining naturally treed peatlands to improve forest growth (e.g., Sarkkola et al. 2004;Remm et al. 2013).In Finland, for instance, drained peatland forests cover more than 20% of the total forest area and comprise almost one-quarter of the total annual forest growth and timber reserves (Natural Resources Institute Finland 2018).The growing bioeconomy is increasing pressure to utilize the timber reserves of peatland forests, but simultaneously, demands to reduce the significant environmental impacts of peatland forestry have increased, particularly on the climate and water environment.
In peatland forests, both tree growth and the environmental challenges generally depend on the depth to the ground water table (WT) in peat, as reviewed by, for example, Nieminen et al. (2018).Even when drained, the WT in boreal peatland forests is generally relatively close to the peat surface (Ahti 1987;Ojanen et al. 2013), and there is a strong bidirectional correlation between the WT and the size of the tree stand (Hökkä et al. 2008;Sarkkola et al. 2010).Negative environmental impacts associated with conventional rotation forestry (RF) on peatlands are generally highest at the extremes of the range of WT variation during the rotation and following management to control the WT.These include the following: the clearcut phase, when the WT generally rises due to greatly decreased evapotranspiration (e.g., Marcotte et al. 2008); the phase when the stand is mature and lowers the WT through high evapotranspiration (Hökkä et al. 2008); and the conventional practice to manage the WT by regular ditch network maintenance (DNM) operations (Sikström and Hökkä 2016).Export load of suspended solids, nutrients, and organic carbon (C) to watercourses related to DNM and clearcuts is high (Nieminen et al. 2010;Kaila et al. 2014).The most productive, nutrient-rich drained peatland forests typically emit more carbon dioxide (CO 2 ) from the peat soil than they sequester from new litter inputs, showing a negative soil C balance (Ojanen et al. 2013).The higher the net CO 2 emission is, the deeper the WT lies (Ojanen and Minkkinen 2019).With increasing utilization pressure, it is critical to find methods to manage peatland forests so that both greenhouse gas emissions and element loading to watercourses are minimized.
Continuous-cover forestry (CCF) could provide a management regime to avoid some environmental detriments, while maintaining or even improving economic performance (Nieminen et al. 2018).It has been suggested that DNM is not necessary in CCF up to certain limits of removal intensity, as the evapotranspiration of the remaining stand maintains the WT sufficiently deep below the peat surface to support economically satisfactory tree growth.With CCF, the adverse effects caused by the excessively high WT following clearcuts, which are largely related to the formation of reductive conditions close to the soil surface (Kaila et al. 2014), could be avoided.Avoiding WTs closer to the soil surface than about 30 cm also keeps the emissions of the potent greenhouse gas methane minimal (Ojanen et al. 2010).Additionally, the WT remaining higher than under the efficient DNM regime would help protect the high soil C stores of peatland forests (Ojanen et al. 2013).
There are different cutting regime options for CCF (Nieminen et al. 2018).For those nutrient-rich peatland forests usually dominated by Norway spruce (Picea abies (L.) Karst.) that are in northern Europe, uneven-aged management with selection cuttings is one option.In selection cuttings, the trees to be felled are selected according to size classes.Typically, the largest trees are selected when the aim is to maximize profitability (e.g., Rämö and Tahvonen 2017).Tree stands in Norway spruce dominated peatland forests generally have a wide diameter distribution (Sarkkola et al. 2003(Sarkkola et al. , 2004)), which assists the transition from RF to uneven-aged management (Juutinen et al. 2018).Yet, currently, there is no information on what would be considered feasible harvest intervals and intensities for CCF in peatland forests.
For any management regime, information on the prospective profitability is critical.Previous studies on the financial performance of peatland forest management have dealt with RF (e.g., Ahtikoski et al. 2012;Hökkä et al. 2017).A recent study by Ahtikoski and Hökkä (2019), for example, analyzed the profitability of four management regimes varying from passive to highly intensive management of Scots pine (Pinus sylvestris L.).The study found that intensive management with DNM and fertilization outperformed the other management options; fertilization was relatively more important for profitability than DNM.Hökkä et al. (2017) emphasized the need to carefully consider the DNM strategy, because DNM may double the negative environmental effects on surface waters in terms of suspended solids and nutrients without considerably improving either financial performance or roundwood supply.
For uneven-aged mineral-soil forests managed with CCF in northern Europe, previous studies have found that an economically optimal steady-state management regime typically involves harvest intervals of 10-20 years with a post-harvest basal area (BA) of 4-10 m 2 •ha À1 .This depends on site characteristics, including temperature sum and fertility, as well as on economic factors such as interest rate level (Tahvonen 2011(Tahvonen , 2016;;Tahvonen and Rämö 2016;Rämö and Tahvonen 2017;Juutinen et al. 2018).These studies also found that in mineral-soil forests, CCF can, in certain cases, outperform RF; however, in peatland forests, the situation is inherently different due to the correlation between the WT and the retained stand that needs to be considered to reach both the economic and environmental gains.
Here, the profitability of managing southern boreal, productive mature spruce forests on peatland are investigated by a wide spectrum of different CCF harvesting alternatives.Selection cuttings were simulated with varying intervals and intensities and compared with a standard RF management option.A matrix of the financial performance against the harvest interval and intensity for each CCF alternative is created following Juutinen et al. (2018).The modelbased analysis yields information on the profitability of DNM in different cases and how the harvest interval and intensity affect the profitability of CCF.While developing management guidance for forest owners, it is important to identify the most profitable regime for timber production, but as forest owners may also have other objectives than timber production, it is equally important to reveal how the harvest interval and intensity affect the financial performance of CCF when wide fluctuations in both are allowed.Such an analysis enables the investigation of the limits of applying CCF in practical forestry.In addition, the economic performance of RF is analyzed in comparison with CCF.The analysis is conducted by using the process-based ecosystem model EFIMOD (Komarov et al. 2003;Shanin et al. 2015).To the researchers' knowledge, this is the first study to analyze the financial performance of CCF in peatland forests.

Model of stand dynamics
The effects of CCF harvestings (here carried out as selection cuttings) on tree stand dynamics were quantified with the individualtree based and spatially explicit forest ecosystem process model EFIMOD (Komarov et al. 2003).The model has been developed for the simulation of C and nitrogen (N) pool dynamics in a tree-soil system, along with population dynamics of forest stands in a wide range of boreal forests.It has previously been validated with empirical data for simulations of CCF in boreal forests on mineral soils (Shanin et al. 2016).A detailed description of EFIMOD can be found in Shanin et al. (2015Shanin et al. ( , 2016)).To apply the model to the simulation of forest dynamics on peatland sites, several additional procedures for linking the ground water table (WT) to stand growth and peat nutrient availability were implemented in this study.The structure of EFIMOD (Fig. 1) and the basic calculation procedures, as well as the following additional procedures that were implemented to apply the original version of the model to simulations of forest dynamics on peatland sites, have been described in more detail (including limitations and uncertainties) in the Supplementary data 1 .

Effect of tree stand on WT
The stand dynamics were linked to the WT level in late summer conditions by adopting the empirical equation by Hökkä et al. (2008) in which the response of the WT to summer precipitation (Sarkkola et al. 2012) was added and the plot random effects of the model were excluded.The final form of the equation, including all empirical coefficients, is where DDitch is the distance to the nearest ditch (m); DDepth is the ditch depth (cm), which is assumed to be equal for all ditches on a simulation plot; V is the local stand volume calculated as the total volume of trees closer than 6 m from the given cell and then recalculated to m 3 •ha À1 ; and MPrec is the mean monthly precipitation during the summer period (June-August; mm).The model was tested with the experimental data presented in Sarkkola et al. (2010) and showed a generally good correspondence (NRMSD = 0.0135, for more details, see Supplementary data 1 ).The value of the WT is calculated for each 0.5 Â 0.5 m 2 cell on the simulation plot.
2.1.3.The response of stand net primary production and regeneration to the WT The regeneration is calculated according to the empirical equations suggested by Pukkala et al. (2012) in which the amount of seedlings and species composition is dependent on total stand density and composition of the canopy layer.These equations were originally proposed for mineral soil sites, but the ingrowth rate estimated for the study sites (depending on stand BA, 40-80 trees•ha À1 reaching the diameter at breast height (1.3 m; DBH) of 0.5 cm annually) fitted well to the experimental data (Sarkkola et al. 2003) on the amount of seedlings of the corresponding size class.The influence of the WT on the net primary production and stand regeneration was defined by the coefficient D, which gains values between 0 and 1 and is calculated based on the threshold values of the WT: where tWT 0 denotes the depth of the WT above which the production process is completely stopped (here set to 0 cm, i.e., the WT at the soil surface); tWT 1 (the value set at 35 cm below the soil surface) denotes the threshold value of the WT below, which does not have a significant effect on stand growth (Sarkkola et al. 2012); and between these two values, the value of D changes linearly.The value of D is first calculated for each cell of the simulation plot, and then D for each individual tree is calculated as a mean value among all cells included in its nutrition zone (Komarov et al. 2003;Shanin et al. 2016).The value of D was used as a reducing multiplier for the biomass production of an individual tree, which is calculated in the original production subroutine in EFIMOD.The same modifier (the whole-plot mean value) was used at the stand level to modify the ingrowth rate.It was also assumed that trees cannot regenerate on ditches, and the roots of existing trees cannot proliferate into ditches.

The response of stand growth to soil nutrient supply
For calculating the response of stand growth to soil nutrient supply, the existing subroutine in EFIMODwhich calculates the biomass increment in relation to N supply using speciesspecific values of relative N consumption (in grams of nutrient per kilogram of biomass production)was extended by the simulation of the growth response to the supply of phosphorus (P) and potassium (K) using the same approach.The hydrological simulation model SUSI was used to generate the input for N, P, and K supply (Laurén et al. 2020; for details, see Supplementary data 1 ).

Proportional timber assortment distributions
To assess the timber assortment distribution and the total volume of harvested wood by assortments, stem taper curve functions were applied using the TapeR package (Kublin et al. 2013) for the R statistical environment (R Core Team 2014), and the curves were calibrated with empirical data on the stem shape (Laasasenaho 1982; for details, see Supplementary data 1 ).

Stand data used in the simulations
A virtual stand was used as input data in the model calibration that was compiled by aggregating the empirical plot-wise tree stand data from two experimental sites in drained Norway spruce dominated peatland forests located in the southern boreal zone of Finland: Janakkala (61°0 0 N, 24°45 0 E, 120 m a.s.l.) and Heinävesi (62°27 0 N, 29°03 0 E, 115 m a.s.l.).Both sites were characterized by mature Norway spruce stands with some admixture of downy birch (Betula pubescens Ehrh.).
The Heinävesi site (2.6 ha) was initially drained in 1963.The site was classified as Vaccinium myrtillus type drained peatland forest (Mtkg I, according to Laine et al. 2018), and the approximately 1 m thick peat layer was Sphagnum peat with woody remains.The tree stand was of a pre-drainage origin with widely varying ages of individual trees.In 1994, the stand BA was 34 m 2 •ha À1 , the dominant height was 20.5 m, and the standing stem volume was $290 m 3 •ha À1 (Repola et al. 2006;site 5018).A thinning experiment (randomized blocks with three replicates) was established in 1993-1994 involving a control (no cuttings) and light, moderate, and heavy thinning, mainly from below, to retained BAs of 26, 21, and 16 m 2 •ha À1 , respectively.Following harvesting, the ditches were cleaned to a standard depth of $0.8 m.The stand parameters were monitored at 5-year intervals until October 2008.During the first 10 years after the establishment of the experiment, the mean annual volume increment of the growing stock varied from 8.2 to 10.4 m 3 •ha À1 annum À1 , depending on treatment (Repola et al. 2006).
The Janakkala site (6 ha) was initially drained in the 1940s when the main ditches were dug.The ditch network was complemented by digging additional drains in the beginning of the 1960s.The site was classified as herb-rich type peatland forest (Rhtkg II), and the >1 m thick peat layer was mostly sedge (Carex spp.) peat.The tree stand originated from advance-growth spruce seedlings evolved after drainage under a birch cover stand that had been harvested earlier.The stand has been managed by successive commercial thinnings from below, the last of which was carried out in the 1980s.Consequently, the stand structure was distinctively even sized, and the stand diameter distribution was close to normal or bimodal with the second peak in the small diameter classes.Spruce and birch seedlings occurred as undergrowth with varying sizes and densities.Monitoring data were not available, with only a one-time measurement describing the stand structure.
For both sites, the following tree characteristics were available: the stem DBH (cm) from all trees and tree height (m) and crown length (m) from 7-15 systematically selected sample trees on each plot.

Calibration and validation of the EFIMOD model
According to the availability of data, the calibration of the model was aimed at adequately representing (i) the stand growth and (ii) the diameter distribution of trees.The calibration involved the parameters of the initial distribution and self-thinning of trees among DBH classes (for details, see Supplementary data 1 ).

Input data and simulation scenarios
Based on the stand characteristics of the two sites (the distribution of trees among DBH classes; stand density, i.e., stem number per hectare; mean distance between trees of certain DBH classes; and the presence or absence of clustering of trees of certain DBH classes), the initial dataset for the management scenarios was generated representing a mixed uneven-aged stand with 1212 trees•ha À1 and a BA of 24.6 m 2 •ha À1 .The DBH of trees varied within 0.25-38.45cm (mean of 13.47 cm), and the tree height was in the range of 1.30-27.24m (mean of 12.88 m).The resulting distribution of trees among DBH classes, being a result of fitting the model to the experimental data from study sites, is bimodal, peaked at middle (15-20 cm) and smaller (<5 cm) size classes.Eighty-nine percent of the BA was spruce and 11% was birch.In both tree species, the size-class distribution was similar.The spatial allocation of individual trees of both species showed no clustering.The distribution of trees among diameter classes was set the same as in the calibration but with setting the number of the smallest trees (DBH less than 4.6 cm) based on the data from the Janakkala site.As an average of the two sites, this "virtual" stand represented the following long-term (1981-2010) climatic conditions: a mean annual temperature of +4 °C, a growing season temperature sum of 1280 dd (>5 °C), and an annual precipitation of 540 mm, of which 200 mm fell during the summer months (June-August).To clearly estimate the effect of thinning parameters, we assumed the climatic conditions to be stable during the whole simulation period.
The age of individual trees was not a driving variable in the scenarios, but it affected some age-dependent functions and parameters in EFIMOD.Therefore, as the age of individual trees was not provided in the datasets, the trees were divided into four age cohorts: 22, 50, 75, and 125 years, according to their DBH, using the threshold values of DBH for attributing a tree to a certain age cohort taken from yield tables (Shvidenko et al. 2008;Usoltsev 2016).
The size of the simulation plot was set at 100 Â 100 m (1 ha).Ditches (1 m wide) were assumed to be parallel, with a distance of 50 m between their mid-lines.To simulate ditch shallowing with time due to ditch bank erosion and ingrowth of vegetation, the following equation was used to determine the ditch depth (Lauhanen et al. 1998): where D age is the age of the ditch (years).The initial ditch depth was set at 100 cm.
The simulation scenarios for the CCF were the same as in Shanin et al. (2016) and Juutinen et al. (2018) and consisted of a series of selection cuttings initiated at the first step of simulation.The period between the two consecutive cuttings was defined by the harvest interval R, set at 10, 15, 20, and 30 years, and the intensity of removal was defined by the post-harvest stand BA, set at 4, 6, 8, 10, 12, 14, and 16 m 2 •ha À1 .Finally, a set of simulation scenarios was obtained by varying both the harvest interval and intensity.In general, cutting assumed the removal of the largest trees, but as the silvicultural recommendations require that some dominant trees need to be retained to ensure sufficient seed production (Äijälä et al. 2019), the retention of 5% of trees was simulated in the largest size class (DBH > 25 cm).The minimum DBH for removed trees was generally set at 16 cm, but the required intensity of removal was checked to ensure that it was reached in all cases to produce comparable results for the different scenarios.The selection of trees for removal was thus primarily based on tree size, but if there were no trees with DBH > 16 cm in the radius of 6 m from the tree selected for removal, this tree was retained and substituted by a tree that was further in the list, if available.This procedure was used to ensure that trees of different size classes occurred in all parts of the stand as homogeneously as possible.Harvesting in EFIMOD is simulated as removal of stems, while branches, foliage, and belowground parts are left on site as felling residues.The total duration of the simulations was set at 240 years, which is divisible by any of the chosen harvest intervals.
For RF, the management simulation was done according to the current Finnish forest management recommendations (Äijälä et al. 2019).The initial stand was the same as in the CCF simulations (1212 trees•ha À1 ; BA, 24.6 m 2 •ha À1 ).Following a clearcut conducted in the beginning of the simulation, a new stand was established by planting 1800 spruce seedlings with spot mounding as a soil preparation method.Then, tending of the sapling stand (time period of 5 years) and pre-commercial thinning (time period of 16 years) were conducted following the recommendations, with removal of the smallest trees (thinning from below).The first commercial thinning occurred when the BA exceeded 26 m 2 •ha À1 at a dominant height of 13 m (this took place at a time period of 60 years) The second thinning followed the same BAdominant height principle (time period of 75 years).Both thinnings were conducted from above when the largest trees were removed.Finally, the stand was clearcut (time period of 89 years) a second time.
Due to the shallowing of ditches, the WT may rise along with the increasing ditch age, and in combination with heavy removals, this may result in WT levels that decrease biomass production and suppress the regeneration of trees.As low-density stands are unable to lower the WT through evapotranspiration, such a situation may lead to the further paludification of the site and a negative feedback loop for stand development.To prevent this, the scenarios assumed that DNM is conducted when necessary.In the simulated CCF, the interval between DNM actions was set to 60 years in the cases of post-harvest BA of 4 and 6 m 2 •ha À1 and to 120 years in the case of post-harvest BA of 8 m 2 •ha À1 .In scenarios with a lower intensity of removal, DNM was not simulated, as in these cases, the stand density was high enough to assumedly maintain the WT at an acceptable level.For comparison, the scenarios with BA of 4, 6, and 8 m 2 •ha À1 were also simulated without DNM.Regarding RF, the DNM was conducted at the same time as the second clearcutting (once per rotation period).

The analysis of profitability
The profitability of management was defined as the net present value (NPV) of harvest revenues.Let the harvest interval (years) be denoted in scenario s by ts ( ts = 10, 15, 20, and 30) and the duration of the simulation period (years) by T, fixed in each scenario (T = 240).The volumes of harvest (m 3 •ha À1 ) at time period t in scenario s are h lts and h pts for sawlog wood and pulpwood (same for all tree species, for clarity).Denote the roadside prices (e•m À3 ) by p l and p p , respectively.The harvest revenues (e•ha À1 ) at period t in scenario s are R ts = p l h lts + p p h pts .Let C ts refer to the variable harvest costs (e•ha À1 ) at period t in scenario s, including cutting and haulage costs: C ts = c ts (h lts + h pts + h wts ), where c ts refers to unit harvest costs (e•m À3 ) and h wts refers to waste wood volume.The fixed harvest costs (e•ha À1 ) at period t are denoted by Ĉt and the cost of DNM is denoted by ĈDNMt .The time interval in which DNM is conducted in scenario s is tDNMs .The interest rate is r.It is assumed that by T, the stand has reached a steady state, and thus the final harvest is repeated every ts years until infinity.Given this notation, the NPV was calculated for CCF scenario s using the following equation: À ĈDNMT e Àr tDNMs 1 À e Àr tDNMs À Á À1 h i Regarding RF, both the ongoing rotation of the initial stand and a bare land value, which captures the NPV of all future rotations, were considered.For bare land, the NPV was calculated as where T* denotes the clearcutting stand age for bare land (T* = 89).In addition to clearcutting, the net revenues include the revenues and costs of thinnings (t = 60, t = 75).The variable harvest costs (e•ha À1 ) are denoted by C tEM .The regeneration costs (e•ha À1 ) at period t are W t , including the costs of site preparation (t = 0), planting (t = 0), tending of a sapling stand (t = 5), and pre-commercial thinning (t = 16).DNM operations are conducted after clearcutting (t = 89).Because the initial stand was mature and was clearcut immediately (t = 0), the NPV for the ongoing rotation was calculated using the following equation:

EM
The variable unit harvest costs for CCF and RF were calculated by using effective time consumption models adjusted for peatland forests from Väätäinen et al. (2010).Unit costs and roadside prices were based on time-series data collected between 2013 and 2017 forming the most recent continuous time series for both costs and roadside prices at the time of this study.After converting the original nominal costs and prices into real terms through deflation (Statistics Finland 2019; cost-of-living index, 1951:10 = 100), the arithmetic average was calculated for each time series.The arithmetic averages were e382.8•haÀ1 for site preparation costs, e633.2•haÀ1 for planting costs (including labour and material costs of seedlings), e374.8•haÀ1 for costs of tending a sapling stand, e429.6•haÀ1 for pre-commercial thinning costs, e391.9•haÀ1 for fertilization, and e210.6•haÀ1 for DNM costs.Then, e57.73•mÀ3 as the roadside price of spruce sawlogs (e30.40•mÀ3 for spruce pulp`wood) and e47.39•mÀ3 as the roadside price of birch sawlogs (e29.58•mÀ3 ) were applied in the analyses.The unit costs for the harvest and haulage were e95•h À1 and e80•h À1 , respectively.The fixed harvest costs were e100•ha À1 .

Model performance
The predicted growth rate generally agreed well with the yield tables (Shvidenko et al. 2008; Usoltsev 2016) for all three control variables including tree height, DBH, and stand density (Fig. 2).The normalized root-mean-square deviation (NRMSD), calculated as RMSD/(max(y) À min(y)), was based on the difference between values predicted by the models and those from the yield tables and comprised 0.0534 for the mean height, 0.0724 for the mean DBH, and 0.0179 for the stand density.It may seem that the model reproduced the stand density remarkably better than the mean DBH and the height of the trees; however, the important fact is that both stand growth series used for validation (Fig. 2) were very close to the model predictions in terms of stand density and slightly differed in mean height and DBH.
The simulated mean DBH growth with EFIMOD was 0.187 cm•annum À1 for the stand age of 50-110 years (corresponding to the stand characteristics of the experimental sites), while the measured annual DBH increment was 0.202 cm•annum À1 for the Heinävesi site.The growing stock mean annual increment (MAI) calculated, taking into account stand development from a young stage to a mature stage without cuttings, also fitted the range of 4.5-6 m 3 •ha À1 reported for the medium yield classes of drained peatlands (Gustavsen et al. 1998).
The model represented the desired diameter distribution quite properly (Fig. 3).The Pearson's x 2 test also confirmed the goodnessof-fit (p value < 0.05).

Stand development and timber production
Regardless of harvesting interval, when the post-harvest BA was 10 m 2 •ha À1 or higher, the WT never rose to levels affecting stand growth in CCF scenarios (i.e., 35 cm below soil surface, Fig. 4).Therefore, in these CCF scenarios, DNM had no effect on the stand growth, and the MAI was between 6.8 and 7.6 m 3 •ha À1 •annum À1 (Fig. 5).However, with BAs of 6 m 2 •ha À1 or lower in CCF scenarios, the WT rose above the critical level.Consequently, stand growth and timber production decreased (Fig. 6), and the MAI remained between 2.1 and 3.6 m 3 •ha À1 •annum À1 (Fig. 5).In these CCF cases, applying DNM maintained the WT below the critical level (Fig. 4) and increased the MAI by up to 100% (Fig. 5).While in CCF cases in which the BA was 8 m 2 •ha À1 , the WT did rise briefly above the level affecting the stand growth (Fig. 4), and the effect of DNM on the MAI was only small (Fig. 5).In the case of RF with DNM, the MAI was 6.3 m 3 •ha À1 •annum À1 .
Generally, in CCF scenarios, the longer the harvesting interval was, the larger the average size of harvested trees was.Similarly, when the post-harvest BA was increased, the trees were harvested at larger sizes.Consequently, the sawlog ratio increased from 61%, with a post-harvest BA of 4 m 2 •ha À1 and an interval of 10 years (R10BA04), up to 82% when the post-harvest BA was 16 m 2 •ha À1 and the interval was 30 years (R30BA16; Fig. 7).For comparison, the sawlog ratio in the RF scenario was 75%.

Profitability
Regarding a 3% interest rate, the NPVs of the harvest net revenues for the scenarios with and without DNM show that DNM was profitable with a post-harvest BA of between 4 and 6 m 2 •ha À1 , regardless of the applied length of the harvest interval (Fig. 8).
The scenario with a BA of 8 m 2 •ha À1 after harvest was a border case, as DNM was typically profitable with harvest intervals of 10-20 years, but not 30 years.In scenarios with a BA of 10 m 2 •ha À1 or higher, the WT never rose over the critical level, thus DNM had no effect on stand growth and therefore was not profitable regardless of the interest rate level.
The profitability of CCF depended strongly on the interest rate.Interestingly, with a 1% interest rate, the highest NPV was associated with the 15-year harvest interval and the 16 m 2 •ha À1 BA scenario (R15BA16), and with a 3% interest rate, the highest NPV was associated with the 15-year harvest interval and the 10 m 2 •ha À1 BA (R15BA10).In these two scenarios, DNM was not profitable.In contrast, with a 5% interest rate, the 15-year harvest interval and the 4 m 2 •ha À1 scenario (R15BA06) with DNM was the most profitable scenario.Further, DNM was profitable in all scenarios with the BA below 8 m 2 •ha À1 with a 5% interest rate.
The profitability of CCF varied between the management scenarios.With a 3% interest rate, for instance, the NPV of the worst scenario with DNM (R30BA04, NPV = e9632•ha À1 ) was 28.4% lower than that of the best scenario (R15BA10, NPV = e11 607•ha À1 ).There were, however, many scenarios in which profitability was relatively close to the best scenario.In the second (R15BA12) and third (R15BA06) best scenarios, the profitability was less than 3% lower than in the best scenario.In fact, in all 15-year harvest interval scenarios, the difference in profitability compared with the best scenario was less than 6%.It is also worth pointing out that with a 3% (and a 5%) interest rate, profitabilities of 10-year harvest interval scenarios were better than 20-year (or 30-year)  Growth data from the Heinävesi site were used in the calibration and compared against yield tables (Shvidenko et al. 2008;Usoltsev 2016; for further information, see Supplementary data 1 ).Grey dots correspond to the data from the yield tables, the solid line represents the model prediction, and the dashed lines represent the prediction interval.harvest interval scenarios, but the opposite was true with a 1% interest rate.
The NPV of RF was e6193•ha À1 with a 3% interest rate.With 1% and 5% interest rates, the figures were e19 841•ha À1 and e4647•ha À1 , respectively.Hence the profitability of RF was generally lower than that of CCF scenarios.Only those two scenarios with the lowest BA after harvest and the shortest (R04BA10) or longest (R04BA30) harvest interval were less profitable than RF when the interest rate was 1%.

Discussion
According to this study, CCF was clearly more profitable than RF for a mature Norway spruce dominated stand on a nutrientrich peatland site.The NPV of the harvest revenues was higher in each CCF scenario than in RF.The outcome was the same with all the interest rate levels except one scenario with a 1% interest rate.Importantly, however, the RF was simulated according to the current forest management recommendations for Norway spruce (Äijälä et al. 2019).The optimal management regime may provide a higher NPV than the recommended regime for RF (e.g., Ahtikoski and Hökkä 2019).While, on mineral soils, WT has no effect on stand growth, some comparisons can be made on a general level.Previous studies on optimized CCF management on mineral soils have found that the profitability of RF is better than that of CCF with low interest rate levels (Tahvonen et al. 2010;Tahvonen 2011).Regarding the mature stand considered in this study, the optimal management would likely include clearcutting in the beginning of simulation similarly as was included in the recommended RF.Hence, a large share of net revenues would be obtained in the beginning of simulation also in the optimal regime, and therefore, it would provide rather similar NPV than Fig. 4. Water table depth (metres below soil surface) for the continuous-cover forestry scenarios of the simulated Norway spruce stand.Cases with ditch network maintenance (DNM) are shown in grey, only cases where DNM had an effect are shown.Dotted horizontal line represents the depth above which the water table is assumed to affect stand growth (35 cm).Thick horizontal lines denote median values, boxes delineate first and third quartiles, and whiskers delineate the range of variation.In the scenario codes, R10 to R30 refer to a harvesting interval from 10 to 30 years and BA04 to BA16 refer to a post-harvest basal area (BA) from 4 m 2 •ha À1 to 16 m 2 •ha À1 , e.g., R15BA10 denotes harvesting interval of 15 years with a post-harvest BA of 10 m 2 •ha À1 .The simulated period was 240 years in each scenario.
Fig. 5. Mean annual increment with and without ditch network maintenance (DNM) for the continuous-cover forestry scenarios and rotation forestry (RF).Codes are as in Fig. 3. the recommended regime for the studied stand that had trees of different ages and sizes initially.In addition, previous studies have found that the financial performance of CCF is relatively high when an initial mature stand is considered (Tahvonen et al. 2010;Juutinen et al. 2018).When the initial stand has an evenaged forest structure, it may take a long period to achieve an uneven-aged stand structure (Wikström 2000;Rämö and Tahvonen 2017), and it may even be optimal to regenerate the stand by clearcutting and not start to apply CCF until in the second tree generation (Tahvonen and Rämö 2016).
The profitability of CCF scenarios depended on the interest rate used also in terms of harvest interval and intensity.Regarding the most profitable scenarios, the harvest intensity increased with interest rate, as has also been found in studies on CCF on mineral soils (Yousefpour and Hanewinkel 2009;Tahvonen and Rämö 2016;Juutinen et al. 2018).However, in this study the harvest interval was not that sensitive to the interest rate as the best financial performance was obtained with a 15-year harvest interval regardless of the interest rate used.With a 3% interest rate, the best scenario included a 15-year harvest interval and a 10 m 2 •ha À1 post-harvest BA, which is in line with previous studies on CCF on mineral soils (Rämö and Tahvonen 2017;Parkatti et al. 2019).The analyses revealed that there were many almost equally profitable CCF scenarios in which the differences in the NPV of harvest revenues were within 10% compared with the best scenario.Schou et al. (2012) found similar results when evaluating different transition scenarios from even-aged stand to CCF.Hence, forest owners have many relatively profitable CCF harvesting alternatives when making decisions about managing their peatland forests.
In this study, the WT dynamics was included in the simulations, and the DNM was subjected to thinning regimes that necessitated the DNM operations due to raised WT.Without such consideration, it is likely that stand projections for the management regimes were biased.The WT influenced the profitability in two ways.First, the WT affected growth (primary production and stand regeneration), particularly when heavy thinning operations were conducted.Heavy thinning operations decreased stand water use potential and caused watering up to such an extent that growth of the retained trees was also reduced.Second, the WT also affected the need of DNM operations.The results indicated that DNM was needed only in CCF scenarios in which the post-harvest BA was less than 10 m 2 •ha À1 , which supports the arguments by Nieminen et al. (2018) and is consistent with findings of Leppä et al. (2020), who showed that selection harvesting to post-harvest BA of 12 m 2 •ha À1 or higher only slightly increased WT.From a viewpoint of profitability, however, the outcome was more complex.With a 3% interest rate, the most profitable CCF scenario did not include DNM, but with a 5% rate, the best scenario included DNM.Hence, the higher the interest rate was, the more intensive the management should be to maximize the NPV of peatland forestry.This is in line with studies on CCF on mineral soils (Parkatti et al. 2019).Similarly, Ahtikoski and Hökkä (2019) found that in RF in pine-dominated peatland forests, the financial performance of intensive management (with Fig. 6.Simulated dynamics of growing stock, volumes of harvested timber, and water table depth (metres below soil surface) for the R15BA04 scenario with and without ditch network maintenance (DNM).R15BA04 means harvesting interval of 15 years with a post-harvest basal area (BA) of 4 m 2 •ha À1 .Critical level represents the depth above which the water table is assumed to affect stand growth (35 cm).
[Colour online.]DNM and fertilization) relative to passive management (without any silvicultural measures) increased with interest rate.
Regarding DNM in CCF scenarios, it should be noted that the stand must reach the minimum post-harvest requirement of 8 m 2 •ha À1 BA as specified in the Forest Act 1308/2013 (Finlex, https://www.finlex.fi/fi/laki/alkup/2013/20131308,accessed 16 March 2020) for the type of examined stand in this study.A stand must be artificially regenerated if the legal minimum requirement is violated.The profitability of CCF would, of course, be lower in this case.Therefore, given the regeneration costs, the results indicate that the best CCF management alternative would be to apply a 15-year harvest interval and a 10 m 2 •ha À1 BA without DNM with 5% (and 3%) interest rates.
Previous studies on RF have found that it is optimal to conduct two DNM operations within a rotation when maximizing NPV in drained peatland Scots pine stands (Ahtikoski and Hökkä 2019).The results showed that only one DNM within a rotation was needed to keep WTD below the critical level in RF in peatland spruce stands.There are probably two equally important reasons for the difference: (i) in the simulation study of Ahtikoski and Hökkä (2019), ditch shallowing and the WT were not explicitly included in the analysis as in this study; and (ii) the higher leaf area of a spruce stand of a specific volume compared with a pine stand, and the better ability of spruce stands to sufficiently sustain a deep WT.
In the timber productivity of CCF compared with RF, earlier research on mineral soils has found mixed results (Kuuluvainen et al. 2012).In the results, apart from scenarios with a post-harvest BA of 4 m 2 •ha À1 , the productivity of CCF exceeded that of RF; however, the difference in the MAI between CCF scenarios and RF was less than 1 m 3 •ha À1 •annum À1 .Compared with earlier studies, the MAI of CCF was on the higher end (e.g., Pukkala et al. 2010;Parkatti et al. 2019).The sawlog ratios of all scenarios, including RF, were in line with those of previous studies (Pukkala et al. 2012;Parkatti et al. 2019).
There is currently little practical experience from CCF in peatland forests.Consequently, there is still some uncertainty related to the regeneration and growth rate of spruce during several consecutive selective cutting cycles with varying post-harvest BAs; however, previous studies on the regeneration of small openings (78-314 m 2 ) suggest that a high number (>10 000) of seedlings can be established on openings without site preparation (Hökkä et al. 2012).Furthermore, they showed that the number of >0.1 m tall spruce seedlings was higher in the small openings than in the clearcuts, where intensive competition with understory vegetation may reduce the growth of seedlings.The understanding of the dynamics of the WT is largely derived from pine-dominated peatland forests (Sarkkola et al. 2010(Sarkkola et al. , 2012)), and only recent studies on Norway spruce dominated mire stands have shown how the standing trees with a higher leaf per volume or BA unit than pine stands impact WT (Leppä et al. 2020).
In this study, the 240-year simulation period was used to reduce the end-point effect in the NPV calculation.The end-point effect emerged because a steady-state harvest was assumed to recur infinitely after the last simulation period in CCF scenarios (eq.4).The longer the simulation period is, the smaller is the influence of the end-point effect due to discounting.Obviously, forest growth and management simulations involve many uncertainties when such a long period is simulated.The applied process-based model can predict tree growth responses to variation of the resource availability over the long simulation periods, including stand development phases in changing climate (e.g., Komarov et al. 2003;Shanin et al. 2013), but in our current model version, stand dynamics were linked to the WT with an empirical equation, which does not support simulations in the changing climate.Leppä et al. (2020) applied a process-based ecohydrological model for simulations of the responses of WT to selection harvests and their findings support our result that WT can be maintained below the critical level by standing trees without DNM operations, and the importance of the tree stand in controlling the WT is likely to increase in the future climate.
In this study, the NPV was used as an indicator of profitability.The NPV provides an accurate basis for comparing different forest management scenarios; however, it requires a set of assumptions.For example, it is assumed that the considered project is marginal, and therefore, the NPV is more suitable for large-scale (risk-neutral) forest owners.The stand-level approach used in this study may be regarded as marginal because Finish forest owners typically own many stands (Juutinen et al. 2020).There are also many different options and assumptions involved in the determination of the interest rate used in calculating the NPV.Several interest rates were used in this study to elaborate the sensitivity of results.In addition, the CCF and RF management scenarios may involve different risks in terms of catastrophic events such as fires and damages due forest deceases that could be also considered when analyzing their profitability.In this case, a utility Fig. 7. Sawnwood-pulpwood ratios (waste wood is not accounted for) for the continuous-cover forestry scenarios, ditch network maintenance was applied when needed.Numbers along the upper axis denote the mean ratio for the whole scenario.Codes are as in Fig. 3. function approach could be used instead of NPV (e.g., Brunette and Couture 2008).

Conclusions
In peatland forests, tree growth decreases if the WT level is too high.The conventional solution to this problem is to maintain adequate drainage artificially by means of DNM.This, however, entails costs and causes negative environmental impacts such as the deterioration of water quality.The findings of this study suggest that CCF is an economically feasible forest management alternative for Norway spruce.Regarding the most profitable CCF alternatives, the stand volume after harvest is high enough to keep the WT below the critical level without DNM operations.In addition, there are many relatively profitable CCF harvesting alternatives for forest owners when they are making decisions about managing their forests on boreal spruce peatland sites.
Future studies are, however, required to examine the regeneration and growth rate of spruce under CCF management on peatland and how the initial peatland stand structure influences the profitability of CCF.Further studies are also needed to investigate what the impacts of CCF are on peatland, GHG emissions, and loading to watercourses.

Fig. 1 .
Fig. 1.Flowchart of the EFIMOD system of models used in the analysis.

Fig. 3 .
Fig. 3.The measured distributions of stem diameter at breast height (1.3 m; DBH) of the Norway spruce stands in the Janakkala and Heinävesi sites (histograms), and the simulated DBH distribution (solid line).The midpoints of each 5 cm DBH class are shown.The smallest size class (DBH of 0-5 cm, midpoint at 2.5 cm) includes only trees with DBH ranging between 4.6 and 5.0 cm, because smaller trees were not inventoried in Heinävesi.During the simulation experiments, the model was initialized using the number of smallest trees available from the Janakkala site.

Fig. 2 .
Fig.2.Calibration of the growth rate used in the simulations of a Norway spruce stand on drained nutrient-rich peatland in southern Finland.Growth data from the Heinävesi site were used in the calibration and compared against yield tables(Shvidenko et al. 2008;Usoltsev 2016; for further information, see Supplementary data 1 ).Grey dots correspond to the data from the yield tables, the solid line represents the model prediction, and the dashed lines represent the prediction interval.

Fig. 8 .
Fig. 8. Net present values (NPV) of the continuous-cover forestry scenarios by interest rate.In the left panel, the scenarios with 4, 6, and 8 m 2 •ha À1 post-harvest basal area (BA) include ditch network maintenance (DNM).