Incorporating a model for ground lichens into multi-functional forest planning for boreal forests in Finland

Forest Ecology and Management 460 (2020) 117912

Contents lists available at ScienceDirect

Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco

Incorporating a model for ground lichens into multi-functional forest planning for boreal forests in Finland

T

Jari Miinaa, , Ville Hallikainenb, Kari Härkönenc, Päivi Meriläd, Tuula Packalena, Pasi Rautiob, Maija Salemaac, Tiina Tonteric, Anne Tolvanend ⁎

a

Natural Resources Institute Finland (Luke), Yliopistokatu 6 B, FI-80100 Joensuu, Finland Natural Resources Institute Finland (Luke), Ounasjoentie 6, FI-96200 Rovaniemi, Finland c Natural Resources Institute Finland (Luke), Latokartanonkaari 9, FI-00790 Helsinki, Finland d Natural Resources Institute Finland (Luke), Paavo Havaksentie 3, FI-90014 University of Oulu, Finland b

ARTICLE INFO

ABSTRACT

Keywords: Model validation Forest management Simulation Epigeic lichens Terricolous lichens Reindeer pastures

The quantitative assessment of forest management impacts on ecosystem services calls for ecological models to be linked to forest planning systems. In these systems, such models are expected to produce valid and accurate predictions to quantify the effects of alternative forest management strategies on ecosystem services. As an example, we present an evaluation of the model for the percent cover of ground lichens in forests on mineral soils in Finland. First, the performance of the lichen cover model fitted to the data of forest vegetation surveyed in 1985–1986 was comprehensively validated with the data resurveyed on the same plots in 1995. The differences between the observed and predicted values in lichen cover were calculated and analysed. Second, stand-level analyses were conducted to validate the predicted development of lichen cover in accordance with the rotation length development and management of forest stands. Third, the impacts of forest management scenarios on the average lichen cover at regional level were calculated to validate the results over larger areas. Lichen cover was predicted by site fertility, growing stock, forest harvests, soil preparation and geographic location. Our evaluation suggested that the model logically predicted the 10-year change in lichen cover at the stand level. A negative bias was found in predictions especially on xeric sites in northern Finland, although this was not related to stand management, but probably to global greening and reindeer grazing. The model based on the crosssectional vegetation survey data showed its ability to describe the long-term changes in lichen cover along with the development and management of stands in different site fertility classes. On xeric and barren sites, the analyses showed decreasing trends in the average lichen cover at the regional level during the 30-year simulation period. However, the development of the average lichen cover varied only slightly among the alternative forest management scenarios. This study provides a modelling approach to evaluate the impacts of forest management on ecosystem services, here the cover of ground lichens, as an indicator of forest biodiversity and quality of reindeer pastures.

1. Introduction Even-aged forest management practices are aimed at sustainable timber production. In boreal forests, from the point of view of ecological sustainability there are great challenges in developing methods and tools for decision making that guarantee preservation of biodiversity and other ecosystem services (Pohjanmies et al., 2017), as forest management causes rapid environmental changes and disrupts natural succession of an ecosystem, including understorey succession (Økland et al., 2003; Hart and Chen, 2006; Tonteri et al., 2016). In the understorey of boreal forests, ground lichens represent a ⁎

distinct species group with special characteristics and ecosystem functions. Lichens contribute to several ecosystem functions e.g. in biogeochemical nutrient cycling and in nutrient webs of northern forests (see Asplund and Wardle, 2017). They also provide insulation and moisture retention on forest floor (see Odland et al., 2018). Lichens may also serve as ecological indicators of forest biodiversity and of the presence of old-growth forests (McMullin and Wiersma, 2019). In addition, lichens provide critical winter forage to semi-domesticated reindeer (Rangifer tarandus tarandus L.), since more than 80% of reindeer winter time diet consists of ground and arboreal lichens when lichens are abundant (Heggberget et al., 2002).

Corresponding author. E-mail address: [email protected] (J. Miina).

https://doi.org/10.1016/j.foreco.2020.117912 Received 19 September 2019; Received in revised form 10 January 2020; Accepted 15 January 2020 0378-1127/ © 2020 Elsevier B.V. All rights reserved.

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

In Finland, lichen-rich forests are found on nutrient-poor, dry mineral soils, which are more common in the northern than in the southern parts of the country. The reindeer lichens (Cladonia spp.) are the most common lichen genus in these forests, other common genera include Cetraria spp., Nephroma spp., Peltigera spp. and Stereocaulon spp. (Reinikainen et al., 2000). The reindeer lichens thrive in old, sparse forests and grow extremely slowly, making them sensitive to mechanical damage and to competition by vegetation. Mechanical disturbance associated with forest management and associated environmental change, along with other environmental changes such as climate warming or air pollution (see Odland et al., 2018), may have negative impacts on reindeer lichens. Recently, forest ground lichens have been reported to have declined in various parts of the boreal and temperate zones (Berg et al., 2008; McMullin et al., 2013; Sandström et al., 2016; Reinecke et al., 2014; Stefanska-Krzaczek et al., 2018). In northern Europe, the decline of lichens has been linked to reindeer grazing and trampling, but also to the effects of forestry practices (Berg et al., 2008; Kivinen et al., 2010; McMullin et al., 2013; Akujärvi et al., 2014; Korosuo et al., 2014; Kumpula et al., 2014; Sandström et al., 2016), a fact that requires further consideration when striving towards sustainable forest management. Due to the many functions and high indicator value of lichens, modelling the response of lichens to forest management and further linkage of the model with stand simulators is highly pertinent, as it enables us to predict their abundance for multi-functional forest management. Lichen cover will provide decision makers with useful datadriven information about the effects of forest management strategies on ecosystem services, thus facilitating decisions in terms of alternative forest management scenarios. Previously, Dettki and Esseen (2003) applied a model for epiphytic lichens to predict the effects of alternative forest management scenarios on such lichens in northern Sweden. Moreover, Korosuo et al. (2014) analysed different forest management scenarios with respect to their suitability for providing reindeer pasture areas (abundance of arboreal and ground lichens). The need for this type of information is growing, along with the increasing awareness about the non-timber benefits of forests (Haara et al., 2018). On sites with similar levels of fertility, the plant and lichen cover may vary over wide ranges depending on e.g. tree species, stand age, stand density and type of harvests (Tonteri et al., 2016; Eldegard et al., 2019). Both intermediate (mostly thinning) and regeneration (mostly clearcutting) harvests open the tree canopy and thus modify the environmental conditions experienced by understorey plants and lichens. In ecological models, stand basal area and categorical variables describing of previous harvest activities have been used as proxies of forest harvest effects on plants or lichens (Miina et al., 2009; Turtiainen et al., 2013, 2016; Tonteri et al., 2016). High spatial variability – both determinate and stochastic – in the abundance of lichens calls for systematic, long-term data from various growing site conditions. Ecological models are frequently based on datasets collected from forests with management histories different from those of tomorrow’s forests. Changing forest management regimes, for example more intensive soil preparation during forest regeneration, may cause delay in the recovery of lichens as the stand ages. Also changes in climate (Hughes, 2000), atmospheric deposition of nitrogen (Dirnböck et al., 2014) and ungulate grazing (Akujärvi et al., 2014; Kumpula et al., 2014) may create new types of relationships between lichens and site and stand characteristics in the future. Because only those site and stand variables that are included in forest planning systems can be used as model predictors, many influencing processes are unmeasured or unknown, and hence, only a part of the total variation in lichen cover can be explained by models. Despite these difficulties, ecological models could be used to support multi-functional forest management decision making in order to maintain the significant roles of understorey vegetation in biodiversity, nutrient cycling and carbon sequestration. However, prior to applying the ecological model,

the model performance should be carefully evaluated. In this context, the aims of the study were (1) to fit a model for ground lichen cover observed at the stand level to be implementable in the decision support systems for forestry; (2) to validate the model by comparing the observed and predicted values of the cover of ground lichens; (3) to validate the performance of model predictions in longterm simulations, and (4) to predict the effects of alternative forest management scenarios on the average cover of ground lichens at the regional level in Finland. 2. Material and methods 2.1. Data on the percent cover of ground lichens The lichen cover model was fitted using the site and stand characteristics and the cover of plant and lichen species in the understorey measured on the permanent sample plots in forests on mineral soils in Finland in 1985–1986 (Reinikainen et al., 2000). These sample plots are a part of a systematic sampling network of 3000 permanent plots established in forests on both mineral soils and peatlands in connection with the 8th Finnish National Forest Inventory (NFI). Plant and ground lichen species, excluding crustose lichens, were identified and their percentage cover was visually estimated on four 2-m2 permanent quadrats within circular 300-m2 sample plots. In the data analysis, species abundances on the four quadrats were averaged for each sample plot. The lichen cover variable was constructed by calculating the total sum of all ground lichen species. The site and stand characteristics on the sample plots were measured and recorded following the field instructions of the NFI. In our analyses, the modelling data (the 8th NFI data, below referred to as NFI85 data) consisted of a total of 1721 sample plots on 822 sample plot clusters in forests on mineral soils in Finland (Table 1). For more details on systematic cluster sampling applied in the vegetation survey, see Reinikainen et al. (2000). It should be noted that the data set of 443 sample plots used by Tonteri et al. Table 1 Main characteristics of the sample plots used in modelling and validating the lichen cover model. Variable

Lichen cover, % Stand age, years Stand basal area, m2ha−1 Proportion of spruce, % Temperature sum, dd

Validating data (NFI95) (1692 sample plots)

Mean

SD

Range

Mean

SD

Range

3.4 63.3 14.4 32.1 1098

9.2 49.3 10.4 35.9 158

0.0–92.8 0–325 0.0–53.8 0–100 651–1360

2.3 62.6 15.8 30.9 1097

6.9 49.4 10.7 35.2 160

0.0–70.4 0–334 0.0–55.9 0–100 651–1360

Variable

N

%

N

%

No harvests in 10 years Regeneration harvest in 0–5 years Regeneration harvest in 6–10 years Intermediate harvest in 0–5 years Intermediate harvest in 6–10 years Drainage in 50 years Paludified sites OMaT (herb-rich sites on brown soils) OMT (herb-rich sites) MT (mesic sites) VT (sub-xeric sites) CT (xeric sites) ClT (barren sites) No soil preparation in 25 years Ploughing in 25 years DNMa in 25 years Other soil preparation in 25 years Reindeer management area

1056 76 81 302 206 98 110 40 270 801 528 76 6 1463 74 92 92 273

61.3 4.4 4.7 17.5 12.0 5.7 6.4 2.3 15.7 46.5 30.7 4.4 0.3 85.0 4.3 5.3 5.3 15.9

1054 105 75 224 234 144 65 51 305 798 443 76 19 1304 91 101 196 272

62.2 6.2 4.4 13.2 13.8 8.5 3.8 3.0 18.0 47.2 26.2 4.5 1.1 77.1 5.4 6.0 11.6 16.1

a

2

Modelling data (NFI85) (1721 sample plots)

DNM = ditch network maintenance.

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

(2016) to analyse the change in the cover of reindeer lichens was a subset of the data used in this study.

quantitatively and qualitatively: data, model design and implementation, as well as model output verification, sensitivity and corroboration (Augusiak et al., 2014). Under the most optimal situation, models are validated with independent data that have not been used for model fitting (Yang et al., 2004). Unfortunately, in many cases, independent data are not available. The lichen cover model was validated with the data of the sample plots resurveyed in 1995. In the model validation, measurements on a total of 1692 sample plots on 816 sample plot clusters in forests on mineral soils in Finland were used (referred to as the NFI95 data) (Table 1). In the model validation, the percent cover of lichens in 1995 was predicted with the fixed part of the lichen cover model (i.e. marginal predictions), using the site and stand variables measured in 1995 as predictors. However, in 1995, site fertility of a total of 494 sample plots (29% of the remeasured sample plots) was reclassified. Because the reclassification of site fertility violated the model validation, we used the site fertility determined in 1985–1986 as a predictor in model validation for these 494 sample plots. The reclassification of site fertility is probably partly due to global greening, i.e. increased site productivity as a response to warming and changes in other drivers such as forest management, nitrogen deposition and/or CO2 concentration (Kauppi et al., 2014). Predicted and measured percent covers of lichens as well as changes in percent cover between the inventories in 1985–1986 and 1995 were plotted against each other, and Pearson's correlation coefficient was calculated. The performance of the models was evaluated by examining the magnitude and distribution of the prediction errors as a function of the predictions as well as the predictors of the model. The aim was to detect any obvious dependences or patterns that indicate systematic discrepancies. The predictive accuracy of the model was determined by characterising the differences between predicted and measured values (e.g. Gauch et al., 2003). The mean of the prediction errors (err = measured - predicted) across all observations was computed as an estimate of the model’s prediction bias (accuracy), and the standard deviation of the prediction errors (SD(err)) was considered as the model’s precision. The statistical significance of the prediction bias was assessed using a twotailed t test. Finally, the mean square error, combining both prediction bias and precision, was obtained as MSE = bias2 + SD(err)2 and RMSE, referring to the root mean square error. The relative error statistics were calculated by dividing the bias, SD(err) and RMSE by the mean value of the measured lichen cover. Also, the proportion of explained variance (R2 = 1 − (VAR(err)/VAR(measured))) was calculated to validate the model predictions.

2.2. Modelling the lichen cover Due to a highly skewed distribution and non-integer values of the response variable, a generalised linear mixed-effects model using a loglink function and a quasi-Poisson distribution assumption was used in modelling the lichen cover. In addition, the response variable had a substantial proportion of its values equal or close to zero. The quasiPoisson distribution was used to account for overdispersion, since the expected value of variance was greater than the expected value of mean. The general multi-level quasi-Poisson model for lichen cover was as follows:

yij ~ Poisson(

ij )

ln( ij ) = f (Xij , ) + ui

(1)

where y is the observed cover of ground lichens on the sample plot (i.e. the mean percent cover of four 2-m2 quadrats on the sample plot, see Tonteri et al. (2016) for details); the conditional distribution of y, given the expected value π, is the Poisson distribution; ln(π) is a log-link function; f(·) is a linear function with fixed predictors Xij and fixed parameters β; subscripts i and j refer to sample plot cluster and sample plot (Reinikainen et al., 2000), respectively; and ui is the normally distributed, between-sample plot cluster random effect with a mean of 0 and constant variance. Variables available in the NFI data were considered as predictors in the model. The predictors included in the model had to be logical and significant at the 0.05 level, and with no systematic errors observed in the residuals. Also non-significant effects of the categorical variable were included in the final model, if the joined Wald χ2 test of the effects was significant. The categorical predictors were formed expressing regeneration and intermediate harvests 0–5 and 6–10 years since harvest (the reference is no harvests in the last 10 years); site fertility classes MT (mesic), VT (sub-xeric), CT (xeric) and ClT (barren sites) (ref. OMaT or OMT, herbrich sites); water management: drainage in the last 50 years and paludified mineral soil sites (ref. no drainage in the last 50 years); soil preparation method: ploughing, ditch network maintenance (DNM) and other (i.e. disc trenching, patching or mounding) (ref. no soil preparation in the last 25 years); and reindeer management area (RMA) if the sample plot is inside the RMA (ref. outside the RMA). The continuous predictors were stand age (years); stand basal area (m2ha−1); the proportion of Norway spruce (Picea abies [L.] H. Karst.) of stand basal area (%); and temperature sum, i.e. the sum of the positive differences between daily mean temperatures and +5 °C (dd). To include the effect of reindeer grazing, the reindeer management area (RMA) (see e.g. Kumpula et al., 2014; Tonteri et al., 2016) was used as a predictor. The RMA consists of almost the entire region of Lapland (excluding Kemi, Tornio and Keminmaa), but also includes several municipalities in the Northern Ostrobothnia region. The RMA is not indicated in the Finnish NFI data or forest simulators. In this study, the RMA was determined as the area north of the line defined by the geographic coordinates of the national uniform coordinate system as follows: latitude = 9992000–0.7841 * longitude. In Eq. (1), the effects of reindeer grazing and trampling on the lichens are described especially by the interaction of RMA and site fertility. The above generalised linear mixed-effects model was fitted by using the function glmmPQL in the R package MASS (Venables and Ripley, 2002).

2.4. Forest planning systems used in model evaluation The calculation of error statistics (e.g. bias and precision) is only a part of model evaluation (Vanclay and Skovsgaard, 1997). In our case, a model for ground lichens should also be able to logically predict the abundance and cover of lichens for any prevailing condition and time horizon along with the development and management of a stand. This kind of evaluation would provide valuable information on the accuracy and predictive ability of the model for model appliers (cf. Schuwirth et al., 2019). The performance of the fitted lichen cover model was evaluated by both stand-level and regional-level simulations. At the stand level, the lichen cover model was applied in predicting the long-term development of lichen cover along with the development of individual stands, and at regional level, the effects of alternative forest management scenarios on the average lichen cover in three regions representing hemiboreal and middle and northern boreal forests in Finland (i.e. Southwest Finland, Northern Ostrobothnia and Lapland) were evaluated. In simulations, the evolution of stand characteristics was predicted

2.3. Validation of the lichen cover model In model evaluation, several elements are examined both 3

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

using the simulators (Motti and MELA) and then the values for these characteristics were used as the predictors of the lichen model. The predictions for lichen cover were computed based on the fixed part of the mixed-effects model by setting the random effects equal to zero. Such marginal predictions were expected to be biased. This bias was then corrected by applying an empirical ratio estimator suggested by Snowdon (1991), so that the marginal predictions were multiplied by the ratio of the mean of the lichen cover observed in the modelling data and the mean of the marginal predictions.

soils in three regions in Finland. The following scenarios based on the optimisation tasks described were used in regional-level calculations (http://mela2.metla.fi/mela/tupa/tupaindex-en.htm, accessed 4 Sep 2019): MAX_ECO: The maximum economic removal scenario defines roundwood and energy wood harvest levels that maximise the net present value without taking into account any sustainability constraints or the demand for wood products. This scenario will harvest all forest stands as soon as they are eligible for harvesting according to the silvicultural recommendations and which do not fulfil the economic prerequisite (5% discount rate) for further growing. This scenario defines the maximum short-term economical timber supply. MAX_SUS: The maximum sustainable removal is defined by maximising the net present value with a 4% discount rate subject to nondeclining periodic total roundwood and energy wood removals, saw log removals and net income. There are no sustainability constraints concerning tree species, harvest methods, age classes or the growth/drain -ratio in order to efficiently use the dynamics of forest structure. Energy wood removal can consist of stems, harvest residues, stumps and roots. BAU: The realised harvest removal scenario outlines the development of forest resources if the roundwood and energy wood harvesting levels realised in 2011–2013 are carried out in the future. In this scenario, the net present value is maximised with a 4% discount rate subject to realised harvest removals by timber assortments. Roundwood removals consist of industrial roundwood and roundwood for households harvested from the studied regions. Energy wood removals consist of forest chips for industrial heat and power used in the region and fuelwood for households. There are no restrictions concerning the removals e.g. by harvest methods, diameter classes or site fertility classes. The target constraints are allowed to vary ± 0.5%.

2.4.1. Stand-level simulations Stand-level simulations were conducted using the Motti simulator developed by Natural Resources Institute Finland (Salminen et al., 2005; Hynynen et al., 2014). Simulations were conducted for hemiboreal and middle and northern boreal forests in Finland, with the temperature sum set to 1300 dd (Turku, Southwest Finland), 1040 dd (Oulu, Northern Ostrobothnia) and 740 dd (Sodankylä, Lapland), respectively. Note that Sodankylä is located inside the reindeer management area (RMA). The development of Norway spruce stands was simulated on herb-rich (OMT, Oxalis-Myrtillus type) and mesic (MT, Myrtillus type) sites and that of Scots pine (Pinus sylvestris L.) stands on mesic (MT, Myrtillus type), sub-xeric (VT, Vaccinium type), xeric (CT, Calluna type) and barren (ClT, Cladonia type) sites. On sub-xeric, xeric and barren sites, pine stands were naturally regenerated, whereas artificial regeneration was used on herb-rich and mesic sites. On herbrich, mesic and sub-xeric sites, soil preparation (other soil preparation in Eq. (1)) was simulated when the initial stands were generated in the Motti simulator. On nutrient-poor sites (CT and ClT), no soil preparation was used. In the Motti simulations, the management of the initial stands followed the silvicultural recommendations applied in privately owned forests in Finland (Äijälä et al., 2019). Depending on the site fertility and geographic location, the stands were pre-commercially thinned and commercially thinned zero to three times before regeneration harvest at the age of 63–135 years. The percent cover of ground lichens was predicted as a function of the site and stand characteristics simulated by Motti. The bias-corrected, marginal predictions at the sample plot level were used as standlevel estimates. In long-term stand-level simulations, we evaluated the effects of site fertility, geographic location (described by temperature sum and RMA) and, particularly, of stand management on the predictions of lichen cover.

3. Results 3.1. Validation of the fitted lichen cover model The most important predictors of the percent cover of ground lichens were site fertility class, spruce proportion of total basal area, stand basal area, stand age, time since regeneration harvests and soil preparation in the NFI85 data (Table 2). The highest lichen cover was found on xeric and barren forests, whereas on fertile sites (herb-rich sites) the cover was scarce. Lichen cover was higher at lower temperature sums (i.e. the more northern locations), irrespective of the site fertility level. In Finnish forest simulators, the long-term average temperature sum is commonly used to describe the climatic conditions at the location of the stand where the temperature sum decreases gradually from southwest towards north and east. However, inside the reindeer management area, lichen cover was lower, especially on xeric and barren sites subjected to reindeer grazing. Increasing stand basal area and proportion of spruce of the total stand basal area decreased lichen abundance, whereas stand age had a positive effect on lichen cover. In forest regeneration, soil preparation methods such as ditching and ploughing negatively affected lichens, but after lighter soil preparation (disc trenching, patching or mounding), the cover of ground lichens was higher than on unprepared sites. In the modelling data (NFI85), the model for the percent cover of lichens gave positively biased predictions (0.98%-units) (Table 3). However, no systematic trends were observed in residuals. The standard deviation of the prediction errors (SD(err)) and the root mean square error (RMSE) were 6.82%-units and 6.89%-units, respectively. The empirical ratio estimator suggested by Snowdon (1991) was 1.41 (=3.36/2.38); this ratio was used to correct the bias of the marginal predictions calculated in stand-level and regional-level simulations. For the resurveyed sample plots (NFI95), model predictions were almost unbiased (Table 3). It should be noted that between the inventories NFI85 and NFI95, the percent cover of lichens decreased on average by 1.1%-units (SD = 5.6%, range = −49.4–28.6%) (cf.

2.4.2. Regional-level simulations The large-scale simulations were carried out using the Finnish forestry dynamics model MELA, which integrates stand-level simulation and forest-level optimisation (Packalen et al., 2017). The MELA stand simulator consists of empirical tree-level models for ingrowth, growth and mortality, similar to those applied in the Motti simulator. The lichen cover model fitted in this study was linked to the MELA program. Sample plot data of the 11th National Forest Inventory (NFI), gained in the years 2009–2013 and covering three regions of Finland, i.e. Southwest Finland, Northern Ostrobothnia and Lapland, representing hemiboreal, middle and northern boreal forests, were used in our analyses. The stand simulator automatically generates several alternative management schedules for the modelling units; here, the NFI sample plots representing the stands. The regional-level estimates of the average lichen cover were calculated by turning the predicted percent cover of lichens into area occupied by lichens on each NFI sample plot and comparing lichen occupied area to the total area represented by the NFI sample plots on each region. Using the MELA program and the NFI data, various harvest scenarios have been defined to calculate the production potential of Finnish forests (Packalen et al., 2017). The existing harvest scenarios were used to illustrate and evaluate the effects of alternative forest management scenarios on the average lichen cover in forests on mineral 4

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

Table 2 The multi-level Poisson model (Eq. (1)) estimated for the mean percent cover of ground lichens on four 2-m2 quadrats on each sample plot in the modelling data (NFI85, N = 1721). χ2 is the joined Wald χ2 test of the categorical fixed effects (type III test), the degrees of freedom in parentheses. Variable

Estimate

Std error

Intercept Harvests (ref. no harvests in 10 years) – Regeneration harvest in 0–5 years – Regeneration harvest in 6–10 years – Intermediate harvest in 0–5 years – Intermediate harvest in 6–10 years Site fertility (ref. OMaT or OMT, herbrich sites) – MT (mesic sites) – VT (sub-xeric sites) – CT or ClT (xeric or barren sites) Water management (ref. no management) – Drainage in 50 years – Paludified sites Stand age, years Stand basal area, m2ha−1 Proportion of spruce, % Temperature sum, dd Soil preparation (ref. no preparation in 25 years) – Ploughing in 25 years – DNMa in 25 years – Other soil preparation in 25 years Reindeer management area (RMA) RMA·Site fertility – RMA·MT – RMA·VT – RMA·CT or ClT Variance component var(u)

0.3087 −1.3174 −0.4584

0.7037 0.44 χ2 (4) =61.07 0.1817 −7.25 0.2094 −2.19

0.015 < 0.001 < 0.001 0.029

−0.0100 −0.0015

0.0929 0.0899

0.914 0.986

a

t-value

−0.11 −0.02

p-value

χ2 (3) =515.61

< 0.001

1.6566 3.2907 4.7726

0.3628 4.57 0.3551 9.27 0.3644 13.10 2 χ (2) =35.35

< 0.001 < 0.001 < 0.001 < 0.001

−0.2997 −1.1770 0.0037 −0.0280 −0.0147 −0.0018

0.2590 −1.16 0.2016 −5.84 0.0007 5.33 0.0049 −5.72 0.0018 −8.25 0.0005 −3.51 2 χ (3) =41.95

0.248 < 0.001 < 0.001 < 0.001 < 0.001 0.001 < 0.001

−1.1255 −1.1981 0.3470 −1.2431

0.2764 −4.07 0.2961 −4.05 0.1432 2.42 2.5717 −0.48 2 χ (3) =19.84 2.5684 0.64 2.5657 0.49 2.5718 0.20

< 0.001 < 0.001 0.016 0.629 < 0.001 0.519 0.622 0.843

1.6554 1.2639 0.5106 0.9522

DNM = ditch network maintenance

Table 3 Main statistics, bias and precision of predictions of the fixed part of the model for the percent cover of ground lichens (Eq. (1), Table 2), as well as Pearson’s correlation coefficients between observations and predictions. Only the sample plots available both in modelling (NFI85) and validation (NFI95) data sets were used in validation (N = 1692). The relative error statistics are in parentheses. Variable

Predicted cover (%) in NFI85

Mean 2.38 SD 4.59 Range 0.0–36.4 Pearson r 0.67 R2 42.1% Bias 0.98* (29%) SD(err) 6.82 (203%) RMSE 6.89 (205%)

Predicted cover (%) in NFI95

Predicted change in cover (%) from NFI85 to NFI95

2.31 4.46 0.0–31.0 0.59 33.9% −0.05 ns (2%) 5.61 (248%) 5.61 (248%)

−0.07 1.38 −9.5–20.9 0.15 1.3% −1.02* (93%) 5.60 (−511%) 5.69 (−519%)

Note: *significant at the 0.05 level;

ns

Fig. 1. Measured vs. predicted percent covers (above) and changes in percent cover (below) of lichens on the sample plots (N = 1692) measured in 1985 (NFI85) and remeasured in 1995 (NFI95). RMA = reindeer management area. Note different scale on axes. R2 is a goodness-of-fit measure for linear trendline, and dashed grey line is the diagonal where measured and predicted values are equal.

units), the bias of the same magnitude and direction was found in the predicted change in lichen cover between the inventories (Table 3). The correlation (0.15) between measured and predicted changes in lichen cover was low (Fig. 2), but there was no systematic bias in relation to predicted change in lichen cover or the predictors of the model.

non-significant at the 0.05 level.

Tonteri et al., 2016). The prediction bias in the modelling data was not found in the validation data, and thus, it seems that the model was not able to predict the average decrease in lichen cover from NFI85 (1985–1986) to NFI95 (1995) (Tables 1 and 3). In general, the model behaved rather similarly, showing positive Pearsońs correlation coefficients of 0.67 and 0.59 between measurements and predictions in the modelling and validation data sets, respectively (Fig. 1). Due to asymmetric distribution and zero or close-to-zero values of the response variable, the model was not able to predict either very low or high values for lichen cover (Fig. 2). Due to the observed decline in the average lichen cover (1.1%-

3.2. Stand-level simulations with lichen cover model The rotation-length development of lichen cover was predicted along with the development and management of stands in different site fertility classes at different locations (Fig. 3). The model for lichen cover logically predicted an increasing cover with decreasing site fertility. The superiority of poor site fertility classes (xeric and barren sites) over other site fertilities, in regard to the cover of ground lichens, was 5

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

Fig. 3. Simulated percent cover of lichens in stands in different fertility classes in Sodankylä (Lapland), Oulu (Northern Ostrobothnia) and Turku (Southwest Finland) using Eq. (1) (Table 2). The stand development was simulated using the Motti simulator, and thinnings and regeneration harvest were simulated according to the silvicultural recommendations. Site fertility (tree species): OMT = herb-rich (spruce), MT = mesic (pine), VT = sub-xeric (pine), CT = xeric (pine) and ClT = barren (pine). Grey and black arrows indicate precommercial and commercial thinnings, respectively.

Fig. 2. Distributions of measured and predicted percent covers and changes in percent cover of lichens on the sample plots (N = 1692) measured in 1985 (NFI85) and remeasured in 1995 (NFI95).

evident at each location studied. In fully stocked managed stands, predicted lichen cover was highest in stands aged 15–50 years, depending on site fertility and location. Especially at southern locations, cover decreased with increasing stand basal area and was only temporarily affected by thinning. On mesic sites, lichen cover in spruce stands was only a quarter of that in pine stands (results not shown). Compared to more fertile sites or more southern locations, on barren sites in the northern region the cover increased with stand age, mainly because forests grow extremely slowly

there and stand density (basal area) in this region does not affect lichen cover as much. On all sites, regeneration harvest decreased the cover to the same level as it was in the beginning of the rotation. 3.3. Effects of forest management scenarios on regional lichen cover At the beginning of the 30-year planning horizon, the average percent cover of ground lichens in forests on mineral soils was estimated to be 3.1, 5.0 and 18.1% in the Southwest Finland, Northern 6

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

Fig. 5. Simulated development of the volume of growing stock (above) and removals (below) in forests on mineral soils in the Lapland, Northern Ostrobothnia and Southwest Finland regions, according to alternative 30-year regional forest scenarios: Maximum economic removal (MAX_ECO), Maximum sustainable removal (MAX_SUS) and Realised harvest removals (BAU).

Finland (4.4%-unit). On xeric and barren sites in Southwest Finland, due to high removals of the MAX_ECO and MAX_SUS scenarios in the beginning of the 30-year period, lichen cover first decreased sharply, but finally was higher than that in the BAU scenario.

Fig. 4. Simulated development of the average lichen cover (%) in forests in the poorest site fertility classes and all fertility classes on mineral soils in the Lapland, Northern Ostrobothnia and Southwest Finland regions in Finland, according to alternative 30-year regional forest scenarios: Maximum economic removal (solid line), Maximum sustainable removal (dashed line) and Realised harvest removals (dotted line).

4. Discussion The model for percent cover of ground lichens was linked to the existing stand-level and regional-level forest simulators. The lichen model was based on data collected from permanent sample plots covering the entire country and quantitatively validated by using the data collected from the same sample plots 10 years later. The performance of the lichen cover model was evaluated – quantitatively and graphically – by predicting and illustrating the development of lichen cover along with the development of forest stands at both stand and regional levels. Thus, the study presents an approach to evaluate ecological models needed to quantify the effects of forest management scenarios on ground lichens. This approach represents a firm basis on which to implement other ecological models in decision support systems for forest management. There is no single approach to evaluate the performance of a mathematical model, but the most appropriate one should be selected depending on intended use. In this study, the evaluation of the lichen model followed an approach involving both quantitative and qualitative aspects as suggested by Vanclay and Skovsgaard (1997). The quality of the lichen model, i.e. the degree to which the model prediction corresponds to the observed lichen cover, was studied in the validation data set. Unfortunately, the data used in model validation

Ostrobothnia and Lapland regions, respectively (Fig. 4). In all scenarios and regions, the average lichen cover decreased by 0.01–0.5%-units during the 30-year simulation period, except in the MAX_ECO scenario in Lapland, in which after a temporary decrease, the average cover increased by 0.1%-units. On average, the development of lichen cover over 30 years varied only slightly among the alternative scenarios (Fig. 4). In all three regions studied, the average lichen cover decreased most in the BAU scenario, in which, due to lower removals, the growing stock volume over the 30-year period increased more than in other scenarios (Fig. 5). Although there were only small changes in the average lichen cover over 30 years, a clear decline in lichens was found on xeric and barren sites, favourable for ground lichens (Fig. 4). Infertile, unproductive sites (rocks, sands, fells, etc.) most favourable for ground lichens are mainly unmanaged, and therefore, the lichen cover on these sites (rocks & sands in Fig. 4) was almost unaffected in alternative scenarios. On xeric and barren sites in Lapland, the most prominent decline (4.1%-unit) in lichens was found in the MAX_ECO scenario and correspondingly in the BAU scenario in Northern Ostrobothnia (5.5%-unit) and in Southwest 7

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

were not completely independent from the modelling data set, because such data were not, and are not, commonly available. The validation data used here most probably provided a more optimistic result than a comparison with completely independent data. However, the validation data set consisting of remeasurements of the sample plots used in model fitting enabled us to validate also the predicted 10-year change in lichen cover. The strong effect of site conditions on lichen cover was expected: the cover on sites of lower fertility and in further northern locations was higher than that of more fertile sites and in the south. Similarly, Kumpula et al. (2014) and Akujärvi et al. (2014) showed a positive relationship between lichen cover and low site fertility in the north. In poor and dry site conditions, lichens can successfully compete with vascular plants (Cornelissen et al., 2001). In the reindeer management area of Finland, decline in lichen cover has been observed earlier by e.g. Kumpula et al. (2014) and Colpaert and Kumpula (2012). We found an acceptable relationship between the model predictions and observations, although the lichen model (marginal predictions) tended to underpredict the cover of lichens (bias 0.98%-units), especially on xeric and barren sites with high lichen cover. The most considerable disagreement between predicted and observed data was related to extreme values, i.e. zero or close-to-zero and very high lichen covers. In our model, site fertility classes were the most important predictors identifying lichen-rich forests. This is the most interesting finding from an ecosystem services point of view, as lichens serve as crucial reindeer pastures, and from the viewpoint of forest biodiversity. Besides the quasi-Poisson distribution used in this study, a binomial distribution with dispersion parameter or zero-inflated negative binomial distribution (using integer values for the response variable) could be applied to ecological modelling (Brooks et al., 2017). An alternative approach could be to fit separate models for zero and non-zero covers. Such hurdle model via truncated Poisson or negative binomial distribution is used, especially if the zeros and non-zeros are assumed to come from the separate data-generating processes (Bolker et al., 2012). We fitted alternative models to our data during the modelling process (using glmmADMB and glmmTMB model builders in R). Models with zero-inflation coefficients for handling excessive zeros did not give considerably better predictions or did not converge. Furthermore, the estimation of a hurdle model failed when we increased the number of predictors by adding variables found to be significant and ecologically meaningful in the model presented in Table 2. The model based on the cross-sectional data collected in 1985 (NFI85) showed an ability to describe the long-term variation in abundance of ground lichens along with the development and management of forest stands in different growing conditions. Model predictions were in accordance with the previous ecological knowledge of ground lichens. It must be noted that our model was fitted for ground lichens as a total group consisting of various species and morphological subgroups. Lichen species have different habitat requirements (Oksanen and Ahti, 1982) and may respond differently to harvesting. Due to the low number of species-level observations, we only present the model for the total cover. Lichen cover may also vary in a different way than does the corresponding lichen biomass (Mattila, 1981). This indicates that the effects of stand management on lichen biomass and the consequent effects on forest ecosystem functions need further research. The selected management schedule (e.g. thinning or regeneration harvest) for each stand will greatly affect the simulated development of the lichen cover (Fig. 3). For example, precommercial thinning, longer rotations and higher thinning intensities during rotation, particularly in stands growing on poor sites, would favour lichen growth (cf. Kivinen et al., 2010). Among the alternative regional forest management scenarios, the average lichen cover varied only slightly during the 30-year planning horizon. Compared to the overall declining trend observed in lichen cover, the effects of the alternative regional forest management scenarios were only minor. The average lichen cover was predicted to decrease in all scenarios and regions (except in the MAX_ECO scenario

in Lapland). This is in accordance with the result of Tonteri et al. (2016). In our study, lichen cover changed with forest development over the rotation period, which confirms earlier observations. In the beginning of the rotation, timber harvesting and soil preparation destroy ground lichens and increase soil nutrient levels due to the decomposition of harvest residues, thereby favouring herbs and graminoids, but suppressing lichens (Kellner and Mårshagen, 1991). However, lichen cover seemed to have an ability to recover even from serious disturbances within 15–20 years (Fig. 3). This was also observed in clear-cuts in west-central Alberta, Canada, where the total cover of lichens peaked 20–30 years after harvesting (Snyder and Woodard 1992). A closing canopy cover suppresses lichen growth by limiting the amount of light on the ground (Kivinen et al., 2010), and thus, thinnings seemed to favour lichen growth temporarily (Fig. 3). Besides shading, litterfall increases during forest succession, which promotes the growth of dwarf shrubs (Sulyma and Coxson, 2001; Stefanska-Krzaczek et al., 2018). Probably due to forest management and increasing growing stock volume, lichen cover, especially on nutrient-poor sites, was predicted to decrease in the forest management scenarios studied. A similar decline in lichen cover has been reported by Reinikainen et al. (2000) and Kumpula et al. (2014). Although the effects of the alternative regional forest management scenarios on the average lichen cover were minor, lichen cover decreased on xeric and barren sites. Growing forests more sparse and applying longer rotation and lighter site preparation might secure ground lichens on poor sites that are important for the quality of reindeer pastures in northern Finland and for preserving and increasing biodiversity also elsewhere. 5. Conclusions The lichen cover model fitted in this study was biologically realistic and acceptably predicted the 10-year change in lichen cover, despite of a bias not related to stand management. Evaluation of the long-term predictions suggested that the model can describe the development of lichen cover in accordance with the development and management of forests in different site fertility classes and in different locations in Finland. The model enabled us to quantify the effects of alternative forest management scenarios on the lichen cover, which is an important ecological indicator of e.g. reindeer pastures and forest biodiversity. Besides determining the predictive accuracy, the performance of the lichen cover model was quantitatively and graphically evaluated in long-term and large-scale simulations. This approach is recommended for evaluating ecological models fitted to be linked to forest simulators, where the models are used to describe the impacts of forest management strategies on ecosystem services. CRediT authorship contribution statement Jari Miina: Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing - original draft, Writing review & editing. Ville Hallikainen: Formal analysis, Investigation, Methodology, Software, Writing - review & editing. Kari Härkönen: Data curation, Formal analysis, Investigation, Methodology, Software, Writing - review & editing. Päivi Merilä: Validation, Writing - review & editing. Tuula Packalen: Conceptualization, Methodology, Writing review & editing. Pasi Rautio: Validation, Writing - review & editing. Maija Salemaa: Validation, Writing - review & editing. Tiina Tonteri: Data curation, Validation, Writing - review & editing. Anne Tolvanen: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing - review & editing. 8

Forest Ecology and Management 460 (2020) 117912

J. Miina, et al.

Conflict of interest

Hynynen, J., Salminen, H., Ahtikoski, A., Huuskonen, S., Ojansuu, R., Siipilehto, J., Lehtonen, M., Rummukainen, A., Kojola, S., Eerikäinen, K., 2014. Scenario Analysis for the Biomass Supply Potential and the Future Development of Finnish Forest Resources. Working Papers of the Finnish Forest Research Institute, pp. 106 p. Kauppi, P., Posch, M., Pirinen, P., 2014. Large impacts of climatic warming on growth of boreal forests since 1960. PLoS One 9, e111340. Kellner, O., Mårshagen, M., 1991. Effects of irrigation and fertilization on the ground vegetation in a 130-year-old stand of Scots pine. Can. J. For. Res. 21, 733–738. Kivinen, S., Moen, J., Berg, A., Eriksson, Å., 2010. Effects of modern forest management on winter grazing resources for reindeer in Sweden. Ambio 39, 269–278. Korosuo, A., Sandström, P., Öhman, K., Eriksson, L.O., 2014. Impacts of different forest management scenarios on forestry and reindeer husbandry. Scand. J. For. Res. 29, 234–251. Kumpula, J., Kurkilahti, M., Helle, T., Colpaert, A., 2014. Both reindeer management and several other land use factors explain the reduction in ground lichens (Cladonia spp.) in pastures grazed by semi-domesticated reindeer in Finland. Reg. Environ. Change 14, 541–559. Mattila, E., 1981. Survey of reindeer winter ranches as a part of the Finnish National Inventory in 1976–1978. Commun. Ins. For. Fenn. 99 (6), 78 p. McMullin, R.T., Wiersma, Y.F., 2019. Out with OLD growth, in with ecological continNEWity: new perspectives on forest conservation. Front. Ecol. Environ. 17, 176–181. McMullin, R.T., Thompson, I.D., Newmaster, S.G., 2013. Lichen conservation in heavily managed boreal forests. Conserv. Biol. 27, 1020–1030. Miina, J., Hotanen, J.-P., Salo, K., 2009. Modelling the abundance and temporal variation in the production of bilberry (Vaccinium myrtillus L.) in Finnish mineral soil forests. Silva Fenn. 43, 577–593. Odland, A., Sundstøl, S., Bjerketvedt, D., 2018. Alpine lichen-dominated heaths: ecology, effects of reindeer grazing, and climate change. A review. Oecol. Mont. 27, 30–50. Økland, T., Rydgren, K., Økland, R.H., Storaunet, K.O., Rolstad, J., 2003. Variation in environmental conditions, understorey species number, abundance and composition among natural and managed Picea abies forest stands. For. Ecol. Manage. 177, 17–37. Oksanen, J., Ahti, T., 1982. Lichen-rich pine forest vegetation in Finland. Ann. Bot. Fennici 19, 275–301. Packalen, T., Korhonen, K.T., Salminen, O., 2017. Finland. In: Barreiro, S., Schelhaas, M.J., McRoberts, R.E., Kändler, G. (Eds.), Forest-Inventory-Based Projection Systems for Wood and Biomass Availability. Managing Forest Ecosystems, pp. 139–148. Pohjanmies, T., Triviño, M., Le Tortorec, E., Salminen, H., Mönkkönen, M., 2017. Conflicting objectives in production forests pose a challenge for forest management. Ecosyst. Serv. 28, 298–310. Reinecke, J., Klemm, G., Heinken, T., 2014. Vegetation change and homogenization of species composition in temperate nutrient deficient Scots pine forests after 45 yr. J. Veg. Sci. 25, 113–121. Reinikainen, A., Mäkipää, R., Vanha-Majamaa, I., Hotanen, J.-.P. (Eds.), 2000. Kasvit muuttuvassa metsäluonnossa. [Changes in the frequency and abundance of forest and mire plants in Finland since 1950]. Tammi, Jyväskylä, pp. 384. Salminen, H., Lehtonen, M., Hynynen, J., 2005. Reusing legacy FORTRAN in the MOTTI growth and yield simulator. Comp. Electr. Agri. 49, 103–113. Sandström, P., Cory, N., Svensson, J., Hedenås, H., Jougda, L., Borchert, N., 2016. On the decline of ground lichen forests in the Swedish boreal landscape: implications for reindeer husbandry and sustainable forest management. Ambio 45, 415–429. Schuwirth, N., Borgwardt, F., Domisch, S., Friedrichs, M., Kattwinkel, M., Kneis, D., Kuemmerlen, M., Langhans, S.D., Martínez-López, J., Vermeiren, P., 2019. How to make ecological models useful for environmental management. Ecol. Model. 411, 108784. Snowdon, P., 1991. A ratio estimator for bias correction in logarithmic regressions. Can. J. For. Res. 21, 720–724. Snyder, J., Woodard, P.M., 1992. Lichen Regeneration Rates in Alberta Following Various Types of Logging and Wildfires Disturbances. Dept. of Renewable Resources, University of Alberta, Edmonton, Alberta, Canada, pp. 118. Stefanska-Krzaczek, E., Fałtynowicz, W., Szypuła, B., Kącki, Z., 2018. Diversity loss of lichen pine forests in Poland. Eur. For. Res. 137, 419–431. Sulyma, R., Coxson, D.S., 2001. Microsite displacement of terrestrial lichens by feather moss mats in late seral pine-lichen woodlands of north-central British Columbia. Bryologist 104, 505–516. Tonteri, T., Salemaa, M., Rautio, P., Hallikainen, V., Korpela, L., Merilä, P., 2016. Forest management regulates temporal change in the cover of boreal plant species. For. Ecol. Manage. 381, 115–124. Turtiainen, M., Miina, J., Salo, K., Hotanen, J.-P., 2013. Empirical prediction models for the coverage and yields of cowberry in Finland. Silva Fenn. 47, 22 article id 1005. Turtiainen, M., Miina, J., Salo, K., Hotanen, J.-P., 2016. Modelling the coverage and annual variation in bilberry yield in Finland. Silva Fenn. 50, 12 article id 1573. Vanclay, J.K., Skovsgaard, J.P., 1997. Evaluating forest growth models. Ecol. Model. 98, 1–12. Venables, W.N., Ripley, B.D., 2002. Modern Applied Statistics with S, 4th ed. Springer, New York, pp. 495. Yang, Y., Monserud, R.A., Huang, S., 2004. An evaluation of diagnostic tests and their roles in validating forest biometric models. Can. J. For. Res. 34, 619–629.

We have no conflicts of interest to disclose. Declaration of competing interest None. Acknowledgements This study was supported by the projects ‘Multiple Forest Use and Ecosystem Services’ and ‘Sustainable Use of Natural Resources’ funded by Luke, and ‘New Products from Forests’ (project A72143) funded by the European Regional Development Fund and executed by Luke, the University of Eastern Finland, Lapland University of Applied Sciences, the Finnish Forest Centre and Oulu University of Applied Sciences. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.foreco.2020.117912. References Äijälä, O., Koistinen, A., Sved, J., Vanhatalo, K., Väisänen, P. (Eds.), 2019. Metsänhoidon suositukset [Forest management recommendations]. Tapion julkaisuja, pp. 253 (In Finnish). Akujärvi, A., Hallikainen, V., Hyppönen, M., Mattila, E., Mikkola, K., Rautio, P., 2014. Effects of reindeer grazing and forestry on ground lichens in Finnish Lapland. Silva Fenn. 48, 18 article id 1153. Asplund, J., Wardle, D.A., 2017. How lichens impact on terrestrial community and ecosystem properties. Biol. Rev. 92, 1720–1738. Augusiak, J., Van den Brink, P.J., Grimm, V., 2014. Merging validation and evaluation of ecological models to ‘evaludation’: a review of terminology and a practical approach. Ecol. Model. 280, 117–128. Berg, A., Östlund, L., Moen, J., Olofsson, J., 2008. A century of logging and forestry in a reindeer herding area in northern Sweden. For. Ecol. Manage. 256, 1009–1020. Bolker, B., Skaug, H., Magnusson, A., Nielsen, A., 2012. Getting started with the glmmADMB package. http://glmmadmb.r-forge.r-project.org/glmmADMB.pdf. Brooks, M.E., Kristensen, K., van Benthem, K.J., Magnusson, A., Berg, C.W., Nielsen, A., Skaug, H.J., Mächler, M., Bolker, B.M., 2017. glmmTMB Balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9 (2), 378–400. Colpaert, A., Kumpula, J., 2012. Detecting changes in the state of reindeer pastures in northernmost Finland, 1995–2005. Polar Rec. 48, 74–82. Cornelissen, J.H.C., Callaghan, T.V., Alatalo, J.M., Michelsen, A., Graglia, E., Hartley, A.E., Hik, D.S., Hobbie, S.E., Press, M.C., Robinson, C.H., Henry, G.H.R., Shaver, G.R., Phoenix, G.K., Gwynn Jones, D., Johansson, S., Chapin III, F.S., Molau, U., Neill, C., Lee, J.A., Melillo, J.M., Sveinbjönrsson, B., Aerts, R., 2001. Global change and arctic ecosystems: Is lichen decline a function of increases in vascular plan biomass? J. Ecol. 89, 984–994. Dettki, H., Esseen, P.A., 2003. Modelling long-term effects of forest management on epiphytic lichens in northern Sweden. For. Ecol. Manage. 175, 223–238. Dirnböck, T., Grandin, U., Bernhardt-Römermann, M., Beudert, B., Canullo, R., Forsius, M., Grabner, M.-T., Holmberg, M., Kleemola, S., Lundin, L., Mirtl, M., Neumann, M., Pompei, E., Salemaa, M., Starlinger, F., Staszewski, T., Uzieblo, A.K., 2014. Forest floor vegetation response to nitrogen deposition in Europe. Glob. Change Biol. 20, 429–440. Eldegard, K., Scholten, J., Stokland, J.N., Granhus, A., Lie, M., 2019. The influence of stand density on bilberry (Vaccinium myrtillus L.) cover depends on stand age, solar irradiation, and tree species composition. For. Ecol. Manage. 432, 582–590. Gauch Jr., H.G., Hwang, J.T.G., Fick, G.W., 2003. Model evaluation by comparison of model-based predictions and measured values. Agron. J. 95, 1442–1446. Haara, A., Pykäläinen, J., Tolvanen, A., Kurttila, M., 2018. Use of interactive data visualization in multi-objective forest planning. J. Environ. Manage. 210, 71–86. Hart, S.A., Chen, H.Y.H., 2006. Understory vegetation dynamics of North American boreal forests. Crit. Rev. Plant Sci. 25, 381–397. Heggberget, T.M., Gaare, E., Ball, J.P., 2002. Reindeer (Rangifer tarandus) and climate change: Importance of winter forage. Rangifer 22, 13–31. Hughes, L., 2000. Biological consequences of global warming: is the signal already. Trends Ecol. Evol. 15, 56–61.

9