Combined effects of natural disturbances and management on forest carbon sequestration: the case of Vaia storm in Italy

Abstract

• Key message

Natural disturbances and management are key drivers for forest carbon balance. We modelled the impact of Vaia storm on forest sink at national scale in Italy. We demonstrate that after Vaia, carbon fluxes among pools and through harvested wood products from salvage logging limit the carbon losses. Our findings can improve the effectiveness of mitigation actions under disturbance scenarios.

• Context

Climate change increasingly modifies frequency and magnitude of extreme events, such as windstorms, with subsequent strong impacts not only on forest health and stability but also on the forest carbon balance.

• Aims

We aim to assess the combined impact of natural disturbances and forest management on the overall forest carbon accounting, including the mitigation potential from harvested wood products.

• Methods

We modelled the impact of Vaia storm on the evolution of forest carbon balance at national scale until 2030. We considered the effect of Vaia storm in combination with current management practices and salvage logging.

• Results

Our results suggest that the overall carbon sink decreased only by 4% due to Vaia, because of internal carbon transfers among forest pools (about 3.1 Mt C from living biomass to dead organic matter), and that the potential negative effects of salvage logging, removing about 1.2 Mt C from dead organic matter, can be counterbalanced by long-term carbon accumulation in harvested wood products.

• Conclusion

Based on our findings, there is an increasing need to robustly consider, through novel approaches (e.g. comprehensive and integrated modelling framework), the effects of natural disturbances in current accounting frameworks, with the final purpose to improve the effectiveness of mitigation strategies in the forestry sector.

1 Introduction

Natural disturbances, such as fires, windstorms, and insect outbreaks, increasingly affect forest health, growth and stability (Thom and Seidl 2016), both at national and local scale. Fires and windstorms are considered the most hazardous events in Mediterranean and Central and Northern European countries, respectively (Gardiner et al. 2010). Natural disturbances, and in particular windstorms, also determine year-by-year variations of about 10–15% of the forest C sink at EU level, and even more, if considered at national scale (Pilli et al. 2017). In particular, wind is projected to damage more than 40 · 10⁶ m³ growing stock year⁻¹ in 2030 in Europe—three times more than fires and insect outbreaks—of which about 3 · 10⁶ m³ growing stock year⁻¹ in the Alpine region (Seidl et al. 2014).

Climate change can modify intensity, frequency and geographical pattern of disturbances, as in the case of extra-tropical cyclones in Europe (Leckebusch et al. 2006). In autumn 2018, an extra-tropical windstorm, named Vaia, broke down at least 9.5 · 10⁶ m³ of merchantable wood in Northern Italy (Motta et al. 2018). Even if not comparable with the same type of events recorded in other European countries, Vaia is considered the major windstorm affecting the Italian territory since the last war (Motta et al. 2018). Indeed, Vaia damages corresponds to more than 70% of the total roundwood removed in Italy in the same year.

The available studies mainly focused on wildfires in Mediterranean countries (e.g. FAO and Plan Blue 2018) and windstorms in Central and Northern European countries (e.g. Gardiner et al. 2013; Forzieri et al. 2020). However, despite some studies addressed the impact of natural disturbances on forest carbon balance in Europe (e.g. Pilli et al. 2016), the combined impact of windstorms and fires on carbon balance at national scale is poorly investigated.

Moreover, there is poor understanding on how to face the immediate and counterfactual effects of natural disturbances and forest management on both ecological and social-economic systems (wood value-chain). This is the case of salvage logging (i.e. the harvesting of trees after natural disturbances; Lindenmayer and Ough 2006), whose interaction with forest ecosystem functions needs to be considered in enhanced policies and management strategies (e.g. Leverkus et al. 2018). Indeed, in the case of climate change mitigation, salvage logging may reduce the stand carbon stock in situ and improve the carbon storage in harvested wood products (e.g. Dobor et al. 2020). Along the forest-value chain, harvested wood products may potentially contribute to substitute non-wood materials for construction and energy uses (Jonsson et al. 2021). However, the mitigation potential of such substitution effect is context-dependent and characterized by large variability and uncertainty (Leskinen et al. 2018).

According to the above-mentioned issues, the present study aims to quantify the potential effect of Vaia storm on the Italian forest carbon sink, by comparing, through a modelling approach, a business-as-usual scenario, based on the historical (until 2018) and theoretical (from 2019 to 2030) evolution of the forest carbon sink, with a second scenario (named Vaia Storm, VS) including (i) the direct effect of the windstorm on the forest carbon pools (living biomass, dead organic matter (DOM), including dead wood and litter, and soil) and (ii) the indirect effects of salvage logging on the harvested wood products (HWP) pool and on the overall forest mitigation potential at national level. Our modelling framework focuses on the immediate contribution of forest management and forest carbon pools, including HWP, to mitigation after a large-scale disturbance, without specifically assessing the mitigation potential that may derive from the substitution effects of harvested wood products over the whole forest value-chain, which is further considered when discussing the overall results.

2 Material and methods

2.1 Study context

The study is downscaled from the national context to the area affected by Vaia on November 2018 in Northern Italy (Fig. 1). Vaia damaged a total forest area of more than 45 kha, distributed between five administrative regions: 7% in Friuli Venezia Giulia (FVG), 10% in Lombardia (LO), and 27% in Veneto (VE) regions, and 13% in Bolzano (BZ) and 43% in Trento (TN) autonomous provinces (Aichner et al. 2019; Chirici et al. 2019; PAT 2019). Vaia mostly affected pure Norway spruce stands, and mixed Norway spruce—Silver fir and Norway spruce—European beech stands. The calculated total damaged growing stock volume corresponded to more than 0.6% of the total volume of growing stock of all Italian forests (e.g. Motta et al. 2018; Chirici et al. 2019).

Fig. 1

Map of Italy, reporting the main regions affected by Vaia, as considered within the present study

2.2 Input data

Main input data are collated from the Italian National Forest and Carbon Inventory 2005 (INFC; Gasparini and Tabacchi 2011). For complete information on the dataset for the whole Italian territory, please refer to Pilli et al. 2013.

In the VS scenario, we used the damaged area as main model input and the amount of biomass affected by storm for model calibration (see Table 1), through combining different sources of information such as remote sensing techniques, field surveys, and expert judgments. From the above-reported information, we obtained, for each administrative unit and forest type (FT, defined by leading species), the following: (i) the damaged forest area (further distinguished between 4 classes), (ii) the amount of merchantable biomass directly damaged, and (iii) the amount of biomass removed by salvage logging in 2019. Table 1 reports the main information used as model input and for model calibration.

Table 1 Input data of area, damaged volume, and amount of salvage logging affected by Vaia as reported by different data sources ([1] PAT (2019); [2] Aichner et al. (2019); [3] Chirici et al. (2019); [4] Bruschini (2019); [5] GeoLab, Florence University, as reported by Cozzarini (2018), [6] by Forzieri et al. (2020), and based on our assumptions [7] within the present study (further information in Appendixes 1, 2, 3, and 4)

As additional inputs, we used the historical harvest intensity—the main driver affecting the forest C sink (Pilli et al. 2016)—and the area affected by wildfires—the main natural disturbance affecting the Italian forests. Both these parameters were defined according to a preliminary, critical assessment of the information provided by the literature (see Appendix 2).

2.3 Modelling framework

We applied the Carbon Budget Model (CBM) developed by the Canadian Forest Service (Kurz et al. 2009), which is fed by data provided by INFC (Gasparini and Tabacchi 2011). CBM was already applied and calibrated for the Italian forests (Pilli et al. 2013) and validated at regional level (for the Trento province, also part of the present study; Pilli et al. 2014). Since we use the same general assumptions applied in Pilli et al. (2013), we provide in Appendix 3 a short description of the model, and hereafter, we highlight only the modelling assumptions specifically tailored to the present study. Input raw data used to run CBM within this specific study are publicly available within a data repository (Pilli and Vizzarri, 2020).

The overall purpose of CBM is to simulate the main fluxes among each forest C pool, starting from the C inflow from the atmosphere (i.e. the Net Primary Production), accounting for all outflows due to natural processes and human activities (Fig. 2). For simulating the effect of ordinary silvicultural treatments, we defined a set of activities which can be carried out for each FT and administrative unit, and calibrated their intensity according to the total amount of harvest reported at country level (Pilli et al. 2013). Within the present study, natural disturbances include both fires and the Vaia storm. The first ones were defined according to the total forest area affected by fires at national level (as reported in Appendix 2, Fig. 6) and simulated by assuming that wildfires affect, on average, 50% of the living biomass with direct CO₂ emissions mainly concentrated on non-merchantable wood components and dead wood. The effect of Vaia was simulated through four different disturbance events, with an increasing share of living biomass moved to DOM pools (Table 1), further distributed between different administrative units and FTs according to the area reported in Table 1. Salvage logging was simulated through the area affected by a specific silvicultural treatment, carried out in 2019, 2020, and 2021 (Table 1), moving to the products pool 80% of the merchantable standing dead-biomass made available by the windstorm occurred in 2018.

Fig. 2

General model framework distinguished between the forest land system simulated through the CBM (grey background) and the products’ system simulated within the HWP module (light brown background). Boxes represent the main carbon pools, including the corresponding average (2000–2030) share of C estimated within the business-as-usual scenario. Arrows represent the main C inflows and outflows among pools due to natural processes (mortality and disturbance events) and human activities (ordinary silvicultural treatments and salvage logging). The average percentage outflow from each pool, and the Net Primary Production to Living Biomass, is reported as proportional to the corresponding C stock attributed to each pool under the business-as-usual scenario for the entire period

Part of the forest biomass that moves to the products pool—fed both by ordinary silvicultural treatments and salvage logging—will be used for energy production (i.e. fuelwood) and part as industrial roundwood. While the C stored within the fuelwood will be directly released to the atmosphere (assuming an instantaneous oxidation, i.e. attributing an half-life decay rate equal to 0, as reported in Fig. 2), a fraction of the C stored within industrial roundwood will be stocked within these products, for a certain number of years, depending on the expected life cycle of each product (defined according to the half-life values reported in Fig. 2, Appendix 4). To account the C stock changes in HWP, we applied the IPCC Tier 2 method (Appendix 4; IPCC 2014). Since 2019, under the business-as-usual scenario, input data required by the IPCC Tier 2 method were assumed as constant, and equal to the average values reported for the last period, consistently with the constant harvest rate applied as input to CBM.

Under the VS scenario an additional amount of harvest mainly used as industrial roundwood (SL_IRW) (at least for the fraction of merchantable biomass) will be available through salvage logging in the period 2019–2021. For this purpose, the biomass removed through salvage logging was further distinguished between branches and other wood components, directly used as fuelwood, broadleaves merchantable dead trees, and coniferous merchantable dead trees. Only this last component, largely predominant according to the preliminary information reported by literature, was transferred to the HWP pool. This amount could be fully accounted as part of the domestic production (increasing the Italian HWP C stock) or partially excluded from this pool, if a fraction of the SL_IRW is directly exported to other countries, as roundwood material.

Because a fraction of this material was exported to other countries (see for example, PAT 2019), we considered different possible scenarios, assuming that an increasing fraction of the SL_IRW (0 to 50%) is exported. This will also indirectly affect the amount of import and export, as well as the fraction of domestic production, as estimated with the IPCC Tier 2 method (further methodological assumptions are reported in Appendix 4).

3 Results

According to our analysis, the total volume damaged by Vaia is equal to about 11.0 · 10⁶ m³, including 9.2 · 10⁶ m³ of merchantable biomass and 1.8 · 10⁶ m³ of other wood components (including branches and tops). This outcome is consistent with the available preliminary estimates (between 10.5 and 14.7 · 10⁶ m³; Table 1), and with the available input data defined at local level (within a range equal to ± 5% at regional level). About 89% of this amount is concentrated within three administrative units: TN (46%), VE (27%), and BZ (16%).

The total amount of salvage logging estimated for 2019 is equal to 3.6 · 10⁶ m³, including about 0.5 · 10⁶ m³ of other wood components. Overall, this amounts to 33% of the total damaged volume estimated in the present study. About 36% of this amount is removed in TN, 29% in BZ, 22% in VE, and the remaining part within the other administrative units. These values are consistent with the values reported for TN (= 27% of the damaged biomass removed; PAT 2019) and BZ (= 62% of the damaged biomass removed; Aichner et al. 2019). According to our assumptions, an additional share of biomass is expected to be removed in 2020 (about 22% of the damaged trees) and in 2021 (about 14% of the damaged trees), while about 30% of the total biomass damaged is to be left on the site.

We estimate that the C sink for the historical period 2000–2017 varies between − 26.8 Mt CO₂ year⁻¹ in 2002 and − 22.4 Mt CO₂ year⁻¹ in 2017 (see Fig. 3a). The average C sink estimated by CBM within this period equals to − 24.7 Mt CO₂ year⁻¹ (± 1.3 Mt CO₂ year⁻¹), consistently with the value reported by ISPRA (2019), i.e. − 26.0 Mt CO₂ year⁻¹ (± 4.7 Mt CO₂ year⁻¹). Nevertheless, the overall trend is in line with the total amount of harvest, and major inter-annual variations are mainly due to the effect of the wildfires occurred within the same period (see for example 2007, 2012, and 2017 in Fig. 6, Appendix 2).

Fig. 3

Carbon sink: plot a reports the evolution of the total forest C sink in the historical period (TOTAL Hist., until 2017; solid red line) and in the future, for both business-as-usual and VS scenario (TOTAL business-as-usual and TOTAL VS, until 2030; dashed red lines), further distinguished between living biomass, dead organic matter (DOM, including dead wood and litter) and soil. Plot b reports the historical (Hist., until 2017) and future (until 2030, under the business-as-usual and VS scenarios) total forest C sink including the HWP mitigation potential (assuming that the total amount of salvage logging is accounted as domestic production). All values are referred to the category Forest Land remaining Forest Land (FL-FL, as reported in (ISPRA 2019), equal to about 7932 kha in 2017, and reported in Gg CO₂ year⁻¹, with negative values highlighting removals from the atmosphere

Under the business-as-usual scenario, the total C sink gradually decreases to − 16.2 Mt CO₂ year⁻¹ in 2030, due to the effect of forest aging, combined with the assumed constant amount of harvest (see Sect. 3.2). However, based on our estimates, in 2018, the living biomass C sink, equal to − 20.2 Mt CO₂ year⁻¹ within the business-as-usual scenario, has been halved by the Vaia storm. The simultaneous transfer of 3.1 Mt C (i.e. about 0.6% of the total living biomass C stock as reported in Fig. 2) to DOM increases the C sink attributed to this last pool, from about − 2.5 Mt CO₂ year⁻¹ to about − 11.0 Mt CO₂ year⁻¹. As a consequence, the overall C sink estimated at national level for 2018 only decreases by 4% between the VS and business-as-usual scenario.

In 2020, while the C sink attributed to living biomass is completely regained (− 1% within the VS scenario compared with the business-as-usual scenario), the DOM pool becomes a C source due to the effect of salvage logging (moving about 1.2% of the DOM C stock to HWP pool, as highlighted in Fig. 2), and subsequently the overall C sink decreases by − 14%, compared with the business-as-usual scenario. For the same reason, even during the following years, we estimate a reduction of the total C sink, in comparison with the business-as-usual scenario, corresponding to − 10% and − 7%, in 2021 and 2022, respectively. Since 2023, the overall effect of Vaia at national level is negligible. For both scenarios, the C stock within the soil pool is quite stable (see Fig. 3), at least within our time horizon, and only partially affected by Vaia in 2018.

While salvage logging decreases the DOM C sink, removing a fraction of the damaged biomass, the same activity has an opposite effect on the HWP pool. Assuming that the total amount of salvage logging is accounted as domestic production (0% of additional export of roundwood material to other countries), in 2020, the HWP mitigation potential increases from + 0.01 Mt CO₂ year⁻¹ (i.e. a C source under the business-as-usual scenario) to − 0.67 Mt CO₂ year⁻¹, under the VS scenario. This C sink further increases to − 1.5 Mt CO₂ year⁻¹ and − 2.3 Mt CO₂ year⁻¹ in 2021 and 2022, respectively. Because these figures compensate the decreasing C sink accounted under the forest pools, at least in 2022, the final C sink estimated at national level including the HWP mitigation potential is 12% and 8% higher than the values attributed to the VS scenario excluding HWP, and to the business-as-usual scenario including the HWP contribution (see Fig. 3b).

If an increased direct export of the roundwood material removed through salvage logging was considered, the HWP mitigation potential would gradually decrease, from − 2.3 Gg CO₂ year⁻¹ excluding any export in 2022 to − 1.7 Gg CO₂ year⁻¹ (i.e. − 29% compared with the previous scenario) assuming 50% of roundwood removed through salvage logging then exported to other countries (see Fig. 4). Nevertheless, even in the latter case, within the VS scenario, the overall C sink at national level including HWP, at least in 2022, is still 5% higher than the business-as-usual scenario.

Fig. 4

HWP mitigation potential: historical (until 2017) and future harvested wood product (HWP) mitigation potential estimated by our study under the business-as-usual scenario, assuming a constant harvest rate until 2030, and the VS scenario, assuming that an increasing fraction of roundwood is exported to other countries, from 0% (VS-Exp0, i.e. no export) to 50% (VS-Exp50)

4 Discussion

We demonstrate the cascading effect induced by forest windstorms on different C pools: carbon in living biomass moves to DOM—partially compensating the decrease in living biomass C sink—and through salvage logging to HWP—partially compensating a decrease in C sink in DOM (see Fig. 3). When fully accounted, the overall effect of this flux is a stabilization of the total C sink, which may slightly decrease (i.e. − 4% in 2018 when the event occurred) or increase (+ 8% in 2022, accounting the maximum compensation from the HWP pool), compared with the business-as-usual scenario. In this last case, where we estimate a reduction of the total C sink assuming a constant amount of harvest to 2030, our results are well in line with other studies, reporting a decreasing C sink due to the ongoing forest aging, as highlighted in other countries in Europe by Nabuurs et al. (2013) and Pilli et al. (2017). For the soil, we estimated a stable C stock. After the natural disturbances, C in living biomass instantaneously moves to dead wood and litter, and from these pools, through a slower process, to soil. Within the CBM model, this decay dynamic is simulated through temperature-dependent decay rates, which determine the fraction of organic matter that decomposes every year into more stable DOM pools (Kurz et al. 2009). Even though natural and human disturbances, including salvage logging and any other silvicultural treatments, can certainly affect the stabilization of the soil C pool at local level (Jandl et al. 2007), at national level, this pool is generally considered in equilibrium within the majority of the European countries (EEA 2019) because the inflow from living biomass and forest residues compensate the outflow to the atmosphere due to the annual decay rate (see Fig. 2).

Of course, our estimates strictly depend on the actual amount of biomass damaged by the Vaia storm and on the final fraction of salvage logging, including the amount that will be removed during the following years. As highlighted in Table 1, these values are quite uncertain at least in some regions (i.e. VE). Although salvage logging reduces C stocks in downed woody debris similarly to fire-impacted areas (e.g. Bradford et al. 2012), it has an important post-disturbance mitigation rule, because of contributing to a lateral C flux through the HWP (see e.g. in beetle-impacted stands; Lamers et al. 2014). At the same time, however, salvage logging may have several interacting effects with natural disturbances, since they have an influence on the risk for subsequent disturbances (e.g. insect outbreaks), the conservation of biodiversity (e.g. preservation of important habitats on-site), and tree regeneration (e.g. timing for forest canopy recovering) (Leverkus et al. 2018).

In such a framework, the use of remotely sensed images could help in assessing more precisely the extent of Vaia damages, and the spatial distribution of salvage cuttings, already put in place or carried out in the future. To this aim, both optical and radar imagery can be used. Optical Sentinel2 images were used by Chirici et al. (2019) to support a first fast assessment of damages, but in the period between the 29 October 2018 and 29 January 2019, more than 50% of the area was always covered by clouds during satellite passages, so no cloud free images were available. For this reason, SAR images (for example from Sentinel1) can be used alternatively for a rapid assessment of windthrow damages (Rüetschi et al. 2019).

Apart from the C sink in HWP, the merchantable biomass removed through salvage logging also provides an additional mitigation potential, due to the substitution benefits realized by using wood products instead of other GHG-intensive materials. This potential is generally accounted through a substitution factor, highly variable (Leskinen et al. 2018), which needs to be assessed case-by-case. Jonsson et al. (2021) demonstrate that HWP and material substitution cannot fully compensate, over the short run, a reduced net C sink in forests, mostly because of the half-life values assigned to HWP (see Fig. 2), generally much shorter than the rotation periods of managed forests in EU (e.g. 35 years vs. 80–120 years). It is also important to highlight that in general, the substitution effects of reduced GHG emissions in other sectors (e.g. energy, industry) may be more effective towards mitigation than the carbon accumulation in HWP pool (e.g. Lippke et al. 2011). An additional amount of biomass, including branches, tops, and other wood residues removed together with the merchantable components, is generally used for energy production and can provide a further mitigation potential when substituting fossil fuels. In order to fully account for all these mitigation potentials and the cascade effects within the entire production chain—i.e. the substitution effects—an integrated model assessment is generally required (Jonsson et al. 2018). A recent study from Leturcq (2020) argued that the carbon footprint of wood products is often underestimated because of a misinterpretation of the concept of “carbon neutrality” attributed to harvested wood products. In reality, a comprehensive modelling framework, such the one applied within our study, can fully account both for the direct (e.g. due to fires or forest residues) and indirect (e.g. consequent to a reduction of the biomass C stock) emissions due to harvest removals, without any “aprioristic” assumption of C neutrality.

The overall estimate of the historical and future C sink is also linked to the total amount of harvest at national level. We assessed this last figure according to the best available data reported in literature (see Appendix 2). However, because of the lack of reliable statistics, especially for the last years (Pilli et al. 2018), these data sources are highly uncertain and likely underestimated (Pra and Pettenella 2016). Even in this case, remote sensing information can complement data collected at national level, to analyze inter-annual variations, especially where detailed statistics are not available (Ceccherini et al. 2020).

The total amount of harvest also affects the relative share of salvage logging compared with ordinary silvicultural treatments and the possible additional effects on the HWP pool. For example, MZP (2019) reports that the salvage logging caused by bark beetle in Czech Republic was 90% of the total harvest in 2018, possibly increasing in 2019, with a reduction of the projected sink in living biomass in the medium term. Unfortunately, to our knowledge, there are not updated official statistics which would allow a comparison between planned and unplanned (salvage) fellings as driven by Vaia storm at national scale in Italy. Based on our findings, however, we highlighted that the relative impact of Vaia storm on the living biomass C stock is comparable with the average annual amount of C moved to HWP due to ordinary silvicultural treatments, as considered within the period 2000–2030, under the Businnes as usual scenario (see Fig. 3).

We showed that the relative contribution of the HWP to the total C sink can compensate the negative C balance attributed to DOM pool, due to salvage logging. The annual inflow to the HWP pool at first depends on the total amount of removals (salvage logging included), secondly on the relative share of industrial round wood (i.e. the fraction of removals used for material products), further distinguished by commodity (sawn wood, wood panels and paper with paperboard), and finally on the fraction of domestic production attributed to each commodity. Usually, depending on scale and time, natural disturbances may generate a distortion of the wood market because of an abrupt amount of wood to be commercialized, which may cause an immediate reduction of related prices (e.g. Holmes 1991). To our knowledge, there are not updated official statistics about the effect of Vaia storm on the wood market (balance between supply and demand, and import/export ratio) and associated prices due to increased salvage logging at national scale. However, regional statistics from TN show that wood prices dropped in 1 year of about 42% and 76% for standing forests and felled trees, respectively (Zanotelli et al. 2019). National and local authorities could implement strategies to reduce this negative effect, such as for example, stocking and gradually allocating wood products in the market over the years after the disturbance event (e.g. wood chips; Pieratti et al. 2019).

The uncertainty related to the adopted specific methodological assumptions, for example on the future amount of salvage logging and on the share of harvest attributed to different commodities, could influence our estimates. Nevertheless, the sensitivity analysis on different level of harvest provided by salvage logging and directly exported to other countries, highlights the need to carefully assess all these quantities, at least at national level. As argued by Johnston and Radeloff (2019) and Sato and Nojiri (2019), due to the assumptions behind the IPCC methods, the traded feedstock (i.e. the amount of salvage logging directly exported to other countries) will not be accounted within the HWP mitigation potential reported in Italy or in other countries.

The overall magnitude of Vaia is not comparable with the major windstorms occurred on other European countries, such as France and Germany, where Lothar and Martin directly damaged between 184 and 204 · 10⁶ m³ at the end of 1999, or Gudrun and Erwin, which damaged about 75 · 10⁶ m³ in northern Europe, in 2005 (Gardiner et al. 2010). In our case, we estimate that in 2018, Vaia reduced the merchantable net annual increment by 36% at national level and damaged about 11 · 10⁶ m³ of merchantable biomass and other wood components. Similar events affected other European countries. For example, in Austria, about 6, 9, and 10 · 10⁶ m³ were damaged by windstorms in 2002, 2007, and 2008, respectively (BNMT 2015). In Czech Republic in the period between 1998 and 2017, the average amount of salvage logging due to windstorm was equal to about 4 · 10⁶ m³ per year (Zahradník and Zahradníková 2019). Despite these countries already accounted for the effect of forest windstorms within their GHG inventories, in case of Italy, such as for other Mediterranean countries, to now the national GHG inventories have mainly focused on the impact of wildfires (ISPRA 2019). According to our estimates, the total CO₂ fires’ emissions affecting the Italian forests between 2000 and 2017 were equal on average to about 0.1 Mt CO₂ year⁻¹, but in 2007 and 2017, when major fires occurred, they amounted to 0.2 and 0.4 Mt CO₂ year⁻¹, respectively. These latter values refer to a burned area comparable to the area affected by Vaia (about 35,000 ha) and correspond to one third of the overall reduction of the national C sink attributed to Vaia in 2018, i.e. 1.0 Mt CO₂ year⁻¹. On the overall GHG balance at national scale, accounting for the impact associated with windstorms may counterbalance the direct emissions from wildfires, depending on the adopted modelling framework (e.g. fraction of biomass burned).

Although not considered in our modelling exercise, assessing the indirect consequences of probable insect outbreaks following wildfires or windstorms on the vitality of survived trees is extremely important. Bark beetle attacks have been always reported after major windstorms in Germany in 1999 (Hanewinkel et al. 2008), Sweden in 2005 and 2007 (Långström et al. 2009), Austria from 1998 to 2007 (Foglar-Deinhardstein et al. 2008), Lithuania between 1994 and 1997 (LRAM and APA 2017), and more recently, in Czech Republic (Cienciala et al. 2019). In these cases, after a windstorm occurred in the winter season, a bark beetle outbreak escaped from the fallen trees during the following summer season, and generally produced the major damages during the second and third seasons following the primary event. The magnitude of these attacks, however, is highly uncertain and usually driven by local climatic conditions during the next spring and summer seasons and by the amount of dead biomass left on site. Disturbance management approaches and high-intensity salvage logging may prevent future insect outbreaks and preserve the C sink in living biomass (e.g. Dobor et al. 2020).

5 Conclusion

We demonstrate that the impact of a natural disturbance on the carbon balance at national scale largely depends on the management following the event (e.g. salvage logging), the carbon fluxes among the pools, and between the pools and the atmosphere. If lateral carbon fluxes (i.e. HWP inflow) are considered, the potential for mitigation may be improved. Advanced modelling approaches are required to fully account for carbon gains and losses due to natural disturbances in both forest lands and wood products. Our analysis does not extend to the transformation and end-use phases in the forest-value chain. However, these aspects need to be incorporated in future research to further understand the implications of a large-scale disturbance not only from the biophysical perspective but also from the social-economic one.

The recent inclusion of forests into the EU climate targets highlights the need to provide policy makers with concrete studies, based on a system perspective on the forest mitigation potential (Grassi et al. 2018). Moreover, the EU Regulation on the LULUCF sector (EU 2018) offers the opportunity to EU Member States of excluding from their future accounts, from 2021 onward, the GHG emissions resulting from natural disturbances that exceed the average within the period 2001–2020. At the same time, forests health and vitality are at risk because of changing climate, along with the benefits they provide, depending on local conditions and management. This means that, within a more comprehensive system perspective, traditional management tools should be integrated with more flexible approaches, such as the example proposed within the present study, based on innovative and multifunctional forest management strategies, taking into account also the unpredictable but increasing effect of natural disturbances, as recently experienced in many European countries.

Appendices

Appendix 1. Input data

1.1 Assumptions on damaged area

Because of missing or incomplete on-ground information, we assumed that the damaged forest area distribution by forest type and damage intensity class for FVG, LO, BZ, and VE follows the same pattern as reported by PAT (2019) for TN (see Tables 2 and 3 below).

Table 2 Relative share of damaged area by forest type (source: PAT 2019)

Table 3 Relative share of damaged area by damage intensity class (source: PAT 2019)

Based on the above-reported parameters, for each administrative region, forest type, and damage intensity class, we estimated the damaged forest area as follow (Eq. 1):

$${\mathrm{DA}}_{i,j,k}={\mathrm{DA}}_{i}\bullet {({\mathrm{FT}}_{\mathrm{sh}})}_{j}\bullet {({\mathrm{DC}}_{\mathrm{sh}})}_{k}$$

(1)

where ${\mathrm{DA}}_{i, j, k}$ is the damaged forest area for the region $i$, the forest type $j$, and the damage intensity class $k$ (ha); ${\mathrm{DA}}_{i}$ is the damaged forest area in region $i$ (ha); ${\mathrm{FT}}_{\mathrm{sh}}$ is the relative share of damaged area by forest type $j$ (%) (see Table 2); ${\mathrm{DC}}_{\mathrm{sh}}$ is the relative share of damaged area by damage intensity class $k$ (%) (see Table 3). ${\mathrm{DA}}_{i}$: Aichner et al. (2019) for BZ; Chirici et al. (2019) for FVG, LO, and VE; and PAT (2019) for TN. ${\mathrm{FT}}_{\mathrm{sh}}$ and ${\mathrm{DC}}_{\mathrm{sh}}$: PAT (2019).

1.2 Assumptions on damaged biomass

Because of missing or incomplete on-ground information, we assumed that the damaged merchantable biomass distribution by damage intensity class for FVG, LO, BZ, and VE follows the same pattern as reported by PAT (2019) for TN (see Table 4). Based on the damage cover considered to define the damage intensity classes, we set a certain amount of damaged volume over the total volume of growing stock.

Table 4 Share of damaged volume and attributed damaged volume over the total volume of growing stock by damage intensity class, as reported for TN (PAT 2019)

To calculate the damaged merchantable biomass, at first, we calculated the expected damaged volume of growing stock per hectare for each region and damage intensity class as follow [Eq. 2]:

$${\mathrm{DV}}_{i,k}={{({\mathrm{DV}}_{\mathrm{REP}})}_{i}\bullet ({\mathrm{DV}}_{\mathrm{sh}})}_{k}$$

(2)

where ${DV}_{i,k}$ is the expected damaged volume of growing stock in the region $i$ for the damage intensity class $k$ (m³), ${\mathrm{DV}}_{\mathrm{REP}}$ is the damaged volume reported in literature for the region $i$ (m³) (see Table 1 for information sources), and ${\mathrm{DV}}_{\mathrm{sh}}$ is the relative share of damaged volume for a specific damage intensity class $k$ (%);

Then, we calculated the attributed damaged volume per hectare for each region and damage intensity class, as follow [Eq. 3]

$${a\mathrm{DV}}_{i,k}=\frac{{\mathrm{DV}}_{i,k}\bullet {\left({H}_{\mathrm{sh}}\right)}_{k}}{{(\sum_{j=1}^{n}{\mathrm{DA}}_{j})}_{i,k}}$$

(3)

where ${a\mathrm{DV}}_{i,k}$ is the attributed damaged volume per hectare in the region $i$ and for the damage intensity class $k$ (m³ ha⁻¹), ${\mathrm{DV}}_{i,k}$ is the expected damaged volume of growing stock in the region $i$ for the damage intensity class $k$ (m³), ${H}_{\mathrm{sh}}$ is the attributed damaged volume over the total volume of growing stock (%), and ${\mathrm{DA}}_{j}$ is the damaged area related to a specific forest type j (ha), obtained from ${\mathrm{DA}}_{i}\bullet {\left({\mathrm{FT}}_{\mathrm{sh}}\right)}_{j}.$ ${\mathrm{DA}}_{i}$ is the damaged forest area in region $i$ (ha). ${\mathrm{DA}}_{i}$: Aichner et al. (2019) for BZ; Chirici et al. (2019) for FVG, LO, and VE; and PAT (2019) for TN.

Finally, we assumed that the distribution of damaged biomass by damage intensity class for BZ, FVG, LO, and VE follows the distribution of damaged biomass by damage intensity class as provided by PAT (2019) for TN. See Table 5 below for details.

$${\mathrm{DV}}_{TN,k}\bullet \mathrm{WBD}\bullet 0.5$$

(4)

Table 5 Biomass (tons C ha⁻¹) by damage intensity class

where ${\mathrm{DV}}_{\mathrm{TN},k}$ is the damaged volume per hectare in TN administrative unit and for the damage intensity class $k$ (m³ ha⁻¹); WBD is the wood basic density (tons m⁻³), assumed to be 0.48 tons m⁻³ as weighted average of individual WBD values by forest type; and 0.5 is the IPCC standard biomass-carbon conversion factor. In the model calibration phase, we used the damaged biomass values from TN (${\mathrm{DB}}_{\mathrm{TN},k}$) to correct the thresholds for maximum allowed biomass to be disturbed for each administrative unit and damage intensity class, as follow [Eq. 5]:

$${\mathrm{DB}}_{i,k}={\mathrm{DB}}_{\mathrm{TN},k}\bullet \left(1+{CF}_{i}\right)$$

(5)

where ${\mathrm{DB}}_{i,k}$ is the maximum allowed biomass per hectare to be disturbed for administrative unit $i$ and damage intensity class $k$ (tons ha⁻¹), ${\mathrm{DB}}_{\mathrm{TN},k}$ is the damaged biomass per hectare in TN and for the damage intensity class $k$ (tons ha⁻¹), and ${\mathrm{CF}}_{i}$ is the correction factor in the administrative region $i$ regardless of damage intensity class $k$, with values ranging from − 0.4 to + 0.5.

1.3 CBM modelling assumptions

Tables 6 and 7 show forest types’ classes and modelling assumptions about carbon transfers among pools.

Table 6 Forest types (FT) and corresponding leading species considered within the present study, modified from Pilli et al. (2013)

Table 7 Main parameters defining the transition matrix applied by CBM model for simulating the effect of the windstorm, with an increasing intensity of the disturbance event, and salvage logging

Appendix 2. Historical harvest intensity and area affected by wildfires

To search for robust estimates of the overall amount of harvest at national level, we compared the data series reported by FAOSTAT (FAOSTAT 2019) with data available from the Italian National Inventory Report 2018 (ISPRA 2018a; i.e. no data in ISPRA 2019). We found the two data sources not fully consistent and comparable along the entire time series (Fig. 5). For this reason, we decided to use the historical harvest intensities for the period 2000–2016 (only available until 2016) as reported by the Italian National Inventory Report 2018 (ISPRA 2018a), also consistent with the Italian National Forestry Accounting Plan (NFAP, ISPRA 2018b). From 2017 onward, under the BaU scenario, we consider a constant harvest rate, equal to the harvest intensity as for 2016 (i.e. about 11 mil m³ year⁻¹), and under the VS scenario, an additional amount of harvest due to the salvage logging carried out from to 2019 to 2021, as reported by literature and estimated in the present study (see Fig. 5). The extent of burned area considered within the model run is reported in Fig. 6.

Fig. 5

total amount of harvest (m³) reported by (ISPRA 2018a, ISPRA 2018b) (solid red line) and FAOSTAT (solid black line, FAOSTAT, December 2019). Based on the amount of harvest reported in 2016 by the first data series used as input within the present study (solid red line, including primary forest), within the BaU scenario, we assumed a constant harvest rate from 2017 onward (dashed red line). Under the VS scenario, an additional amount of harvest is provided through salvage logging in 2019, 2020, and 2021 (triangles). In order to compare the relative impact of Vaia, against forest wildfires, which are the main natural disturbances affecting Italian forests, we also included in our assessment the total amount of burned area as reported by the Italian National Inventory Report (CRF Table 4(V), ISPRA 2019) in the historical period 2000–2017, then assumed as constant, and equal to the average historical burned area, until 2030, for both our scenarios (see Fig. 6)

Fig. 6

extent of burned area (ha) in the historical period (red bars) (CRF Table 4(V); (ISPRA 2019)), as reported for the category FL-FL) and projected (blue bars) at national scale

Appendix 3. Carbon budget model

The CBM is an inventory-based, yield-data driven model that simulates the stand- and landscape-level C dynamics of above- and below-ground biomass, dead organic matter (DOM), and mineral soil (Kurz et al. 2009). The spatial framework applied by the model conceptually follows Reporting Method 1 (IPCC 2014) in which the spatial units are defined by their geographic boundaries and all forest stands are geographically referenced to a spatial unit (SPU). Within a SPU, each forest stand is characterized by age, area, and 7 classifiers that provide administrative and ecological information, the link to the appropriate yield tables (YTs), and other parameters defining the forest composition (corresponding to 20 forest types—FTs—based on the leading species reported by INFC—see Table 6), the management type (MT, i.e. even-aged high forests, uneven-aged high forests, coppices) and the specific silvicultural system applied to each MT and FT (such as clear-cuts, thinnings, shelterwood systems, partial cuttings).

Species-specific, stand-level equations (Boudewyn et al. 2007) convert merchantable volume production into aboveground biomass, partitioned into merchantable stemwood, other components (tops, branches, sub-merchantable size trees) and foliage components (Kurz et al. 2009). At this purpose, each of the Italian FT has been associated to an appropriate Canadian species following the approach described in Pilli et al. (2013). The merchantable volume at the beginning of the model run is derived from the standing volumes per age class reported by the INFC for each FT and administrative unit. During the simulation, the stand-level volume accumulation is derived from the current annual increment reported by INFC, further distinguished by FT and administrative unit ( Pilli et al. 2013). Belowground biomass is calculated using the equations provided by Li et al. (2003), and the annual dead wood and foliage input is estimated as a pool-specific turnover rate (percentage) applied to the standing biomass stock.

The model uses an initialization process to estimate the size of all DOM pools at the start of the simulation and then, following IPCC guidance, links DOM dynamics to biomass dynamics. Inputs from biomass to DOM pools result from biomass litterfall and turnover as well as natural and human-caused disturbances. The DOM parameters have been calibrated within a previous study (see Pilli et al. 2013), then validated at regional level (Pilli et al. 2014). The model used as input the age class distribution reported by INFC for 2005, further distinguished between FTs and administrative units, according to the approach already described in Pilli et al. (2013). During the model run, anthropogenic and natural disturbances are defined through the amount (area or C target), type and intensity of each disturbance by year and SPU (Kull et al. 2019). Eligibility criteria, such as FT, age, or other classifier values can be used to define the eligible stands for each disturbance. Disturbance impacts are defined using a “disturbance matrix” that describes the proportion of C transferred between pools (i.e. from living biomass to DOM, due to the effect of windstorms), transferred to the product pool (within ordinary silvicultural treatments or through salvage logging after windstorms) or released to the atmosphere (with fires) for each disturbance type (Kurz et al. 2009, Kull et al. 2019). The Archive Index Data Base customized for EU countries and applied within this specific study is available at the following URL (Pilli et al. 2018): https://ec.europa.eu/jrc/en/scientific-tool/eu-archive-index-database-customised-carbon-budget-model-cbm-cfs3

Appendix 4. HWP modelling assumptions

The so-called IPCC Tier 2 method, considered within our study, applies first-order decay functions based on default half-lives numbers to the share of HWP originating from domestic products, generally distinguished between the following commodities: (i) sawn wood, (ii) wood-panels, and (iii) paper or paperboards (see Pilli et al. (2015) for further details).

Inflow–outflow methods estimate the changes in carbon stocks by counting the amount of wood products into and out of the stock. Changes in carbon stocks in year $i$ are estimated on the basis of information (i) on the inflow of wood products into the stock and (ii) of assumed lifetimes and (iii) decay factors of these products (Lifetime analysis), according to the following equations (see equations two-eight-five in IPCC (2014)):

$$C\left(i+1\right)={e}^{-k}C\left(i\right)+\left[\frac{1-{e}^{-k}}{k}\right]*\mathrm{Inflow}(i)$$

(6)

where $\mathrm{C}\left(\mathrm{i}\right)$ is the carbon stock in the particular HWP category at the beginning of year $i$, in Gg C (default conversion factors are in the table below); $k$ is the first-order decay constant for each HWP category, equal to $\mathrm{ln}\left(2\right)/\mathrm{HL}$, where $\mathrm{HL}$ is the default half-lives value of each HWP category in years (see table below); and $\mathrm{Inflow}(i)$ is the inflow to the particular HWP category during the year $i$, in Gg C year⁻¹.

HWP categories	Sawn wood		Wood based panels	Paper and paper boards
	Coniferous	Non-coniferous
Conversion factor per air dry density	0.225 tons C m−3	0.280 tons C m−3	0.269 tons C m−3	0.386 tons C tons−1
Default half-lives	35 years		25 years	2 years

The carbon stock change ($\Delta C\left(i\right)$ in Gg C year⁻¹) of the HWP category during the year $i$ is equal to

$$\Delta C=C\left(i+1\right)-C(i)$$

(7)

Equations (6) and (7) were applied separately for each semi-finished wood products category (sawn wood coniferous and non-coniferous, wood-based panels, and paper and paperboards).

All the input data used within this approach are based on production, import, and export quantities of each commodity (IPCC, 2014) The fraction of domestic production (f_DP) is computed considering, for each year (i), the ratio between the IRW production of each commodity originating from domestic forests (IRW_P) and the sum of production and import (IRW_IM), both curtail of the quota of export (IRW_EX), as reported below:

$${f}_{\mathrm{DP}}\left(i\right)=\frac{{\mathrm{IRW}}_{P}\left(i\right)-{\mathrm{IRW}}_{\mathrm{EX}}\left(i\right)}{{\mathrm{IRW}}_{P}\left(i\right)+{\mathrm{IRW}}_{\mathrm{IM}}\left(i\right)-{\mathrm{IRW}}_{\mathrm{EX}}\left(i\right)}$$

(8)

The denominator of Eq. (8) equals the consumption. For the historical period 2000–2018, all the parameters reported in Eq. (8) are based on the time series provided by FAOSTAT.

Since we assumed that only the coniferous merchantable dead trees will be transferred to the HWP pool through salvage logging, the following methodological assumptions only concern the coniferous component of this pool. For all the components not specified below, such as for the general computation of the overall pool, we applied the default IPCC approach (see Eq. (6) and (7)), as described in IPCC (2014), and in Pilli et al. (2015).

The parameters reported in Eq. (8) were modelled according to the following methodological assumptions (summarized in Table 8):

• For BaU scenario: all parameters were maintained constant (also for broadleaves) and equal to the last values reported by FAOSTAT in 2018 (last historical value).

• For all VS scenarios: for each year from 2020 to 2023 (assuming 1-year delay for transferring the roundwood material to the product pool, after removals) the coniferous IRW_P was increased proportionally to the relative variation of harvest due to salvage logging, compared with the BaU scenario. Broadleaves’s parameters were not modified.

• For VS-Exp₀: for conifers, IRW_EX was increased and IRW_IM was decreased, according to the same proportions applied to IRW_P. Broadleaves’s parameters were not modified.

• For VS-Exp₁₀ to VS-Exp₅₀: for conifers, IRW_EX was increased and IRW_IM was decreased, according to the same proportions applied to IRW_P, further reduced to account for an increasing share of roundwood directly export to other countries. We assumed that this share decreases each year, as reported in Table 8. Broadleaves’s parameters were not modified.

Table 8 For each modelling scenario and year (assuming 1-year delay for transferring roundwood to the product pool, after removals), the table reports the fraction of roundwood directly exported to other countries and the consequent relative variation on the amount of industrial roundwood production (IRWP), export (IRWEX), and import (IRWIM) applied to Eq. (8), in comparison with the value applied in BaU scenario

Fig. 7

from top to bottom, the evolution of industrial roundwood production (IRW_P), export (IRW_EXP), and import (IRW_IM) for coniferous species reported by FAOSTAT for the historical period (i.e. until 2018), and estimated in the present study. Dots (and gradations of grey) represent the variations in the assumptions for the Vaia storms effects (VS-Exp), ranging from 0% (light grey dots) to 50% of industrial roundwood directly exported (see also Table 8)

Fig. 8

The evolution of the fraction of domestic production as reported by FAOSTAT for the historical period (i.e. until 2018; bolded dashed line), and estimated in the present study, according to the BaU scenario (dashed line). Dots (and gradations of grey) represent the variations in the assumptions for the Vaia storms effects (VS-Exp), ranging from 0% (light grey dots) to 50% of industrial roundwood directly exported (see also Table 8)