Economic impacts of ambient ozone pollution on wood production in Italy

Worldwide, tropospheric ozone (O3) is a potential threat to wood production, but our understanding of O3 economic impacts on forests is still limited. To overcome this issue, we developed an approach for integrating O3 risk modelling and economic estimates, by using the Italian forests as a case study. Results suggested a significant impact of O3 expressed in terms of stomatal flux with an hourly threshold of uptake (Y = 1 nmol O3 m−2 leaf area s−1 to represent the detoxification capacity of trees), i.e. POD1. In 2005, the annual POD1 averaged over Italy was 20.4 mmol m−2 and the consequent potential damage ranged from 790.90 M€ to 2.85 B€ of capital value (i.e. 255–869 € ha−1, on average) depending on the interest rate. The annual damage ranged from 31.6 to 57.1 M€ (i.e. 10–17 € ha−1 per year, on average). There was also a 1.1% reduction in the profitable forest areas, i.e. with a positive Forest Expectation Value (FEV), with significant declines of the annual national wood production of firewood (− 7.5%), timber pole (− 7.4%), roundwood (− 5.0%) and paper mill (− 4.8%). Results were significantly different in the different Italian regions. We recommend our combined approach for further studies under different economic and phytoclimatic conditions.

Tropospheric ozone (O 3 ) pollution affects large areas of the world 1,2 . Ozone is strongly phytotoxic and is considered as a serious issue for the health and productivity of forests 3 . Ozone risk assessment may use different O 3 metrics 4 or models 1 . The most common approach is the use of exposure-based metrics e.g. AOT40 i.e. the accumulation of hourly O 3 concentrations above 40 ppb for daylight hours during the growing season 4 . However, exposure-based metrics do not incorporate the effects of stomata, that are the only way of O 3 entry into the plant 5 . Therefore, a flux-based approach is recommended e.g. by the Convention on Long Range Transboundary Air Pollution of the United Nations 6 and by the National Emission Ceilings Directive of the European Union 7 , where the stomatal O 3 uptake is estimated through models integrating the effects of climatic factors and vegetation characteristics on stomata (e.g. the DO3SE model 8 ). Such flux metric is called Phytotoxic Ozone Dose, defined as the amount of O 3 absorbed into the leaves or needles through stomata over the growing season, and above a threshold Y of uptake (PODY).
Among the many ecosystem services provided by forests, only wood production losses have been estimated so far to assess the economic impact of O 3 pollution on forests [9][10][11] , because experimental dose-response relationships are available for estimating biomass losses [12][13][14] . Feng et al. 11 used the AOT40-response relationship and reported that current levels of O 3 across China may cause economic losses of forest production equivalent to 52.2 billion US$ in 2015. Karlsson et al. 9 similarly applied the AOT40-based approach to Swedish forests and simulated that the potential annual economic loss of forest production due to O 3 was 56.0 million € over the time period 1993-2003. Felzer et al. 10 predicted the O 3 effects on carbon sequestration in crops, pastures and forests at global scale by using the empirical function of AOT40 for gross primary production (GPP) and found that reduced CO 2 uptake due to O 3 exposure would increase the macroeconomic consumption cost of the greenhouse gas policy by 4.5 trillion US$ in 2100. Wood is central in a sustainable bioeconomy, e.g. forestry and extended wood-based value chains employed 4.5 million people in the European Union in 2018 15 .
Given the limited knowledge on the economic value of O 3 -induced wood production losses, in particular according to a flux-based approach, the aim of this study was to develop a combined modelling approach for realistic estimates of the economic impacts of POD1-based wood losses, by using Italy as a case study. Italy is a well-known hot-spot of O 3 pollution, given its central position in the Mediterranean basin in Southern

Results
POD1 values across Italy ranged between 0.3 and 100 mmol m −2 , with a spatial distribution showing lower values in the Alps mountains in the North and higher values over the peninsula and the islands where the climate is typically Mediterranean (Fig. 1A). The average value of POD1 was 20.4 mmol m −2 ranging from 8.0 mmol m −2 in Piedmont to 41.4 mmol m −2 in Sardinia (Fig. 1S).
Results suggested a significant impact of O 3 (POD1) on Italian wood losses, that reached almost half of the expected wood production in the no ozone scenario (WO) at individual grid points (Fig. 1B) as well as a significant influence of the interest rate (r) on economic trend. Indeed, the total capital value of Italian forests in the WO scenario (no ozone) ranged from 8.0 to 29.5 B€ and from 2578 to 8987 € ha −1 with 2-4% of interest rate ( Table 1). The POD1 scenario caused a total potential damage from 790.9 M€ to 2.85 B€ of capital value (255-869 € ha −1 , on average). The annual damage ranged from 31.6 to 57.1 M€ (10-17 € ha −1 year −1 , on average). The relative economic impact from WO to POD1 scenario was about 10% (total damage) and about 9% (average damage) for all interest rates. In addition, the total forest area with FEV > 0 decreased by around 1.1% from WO to POD1 (Table 2). The reduction ranged from 35,358 ha with r = 3% to 37,774 ha with r = 4% (1.1% and 1.2% of total forest area, respectively).
The O 3 impact on potential national wood production is reported in Table 3. According to the forest characteristics of Italy (widespread broadleaved stands, coppice management, private properties, etc.), the most represented wood assortment was firewood (from 5,449,258 to 5,603,228 m 3 year −1 in WO scenario), followed by timber pole (from 818,390 to 824,274 m 3 year −1 ), roundwood (from 587,011 to 672,446 m 3 year −1 ) and paper mill (from 288,511 to 365,174 m 3 year −1 ). Firewood and timber pole were strongly affected by moving from WO to POD1 scenario (on average the reduction was 7.5% and 7.4%, respectively). The decrease of roundwood and paper mill (on average 5.0% and 4.8%, respectively) was lower than that for timber pole and firewood.

Discussion
We merged an open-source add-on GIS software tool for estimating the economic value of forests 19 with a classic O 3 risk assessment approach 6,20 , and used forest inventory and WRF-CHIMERE outputs to spatially estimate the O 3 -induced wood losses in the year 2005. We used the accumulated stomatal O 3 uptake (POD1) as a metric of O 3 damage because it is considered as a better index than only O 3 concentrations in the air 21 . However, estimating POD1 at fine scale is challenging as it requires hourly inputs 22 . For the first time, an economic valuation of wood losses was based on POD1 at high spatial horizontal resolution (12 km 2 ).
Italy is known to be subject to elevated O 3 pollution 16 and is thus an ideal case study for applying this combined approach. In fact, the average POD1 value of Italy was 20.  28 for the protection of deciduous broadleaves (birch and beech) and conifers (Norway spruce), respectively, stressing that most of Italian forests are exposed to severe O 3 risks.
Modelling O 3 impact on Italian forests showed a marked decrease of their economic value. According to applied interest rate, annual damage ranged from 31.6 to 57.0 M€ per year with a loss of capital value of about 10%. Among the few previous papers on the economic impact of O 3 on forests, none used a methodology comparable with the present research [9][10][11]. However, Karlsson et al. 9 indicated similar results as the potential annual economic loss for Sweden due to negative impacts of O 3 on forest production was of the order of 56 M€ that is about 2.4 € ha −1 year −1 , while in Italy it ranged from 3.02 to 5.46 € ha −1 year −1 (r: 2-4%). The variability of potential economic impact among Italian administrative regions reflected the strong dissimilarity of geomorphological, logistic, vegetational as well as socio-economic conditions of Italian forests 29 . We innovatively defined the economic damage by O 3 also in terms of reduction of forest area with positive FEV. This reduction in the POD1 scenario was about 1.1-1.2% of the total forest area. This loss of stands with economic profitability may make active forest management no longer meaningful, which would result in an indirect negative effect due to the worsening of other ecosystem services 30 . As an example, hydrogeological problems or fire risk-which are relevant issues for the Italian territory-can increase due to a decline of silvicultural practices 31 . In general, trade-offs among provisioning and other forest benefits (regulating, supporting and cultural services 32 ) can happen and they may be measured not only in biophysical but also in economic terms 33 . A POD1-induced loss from about 535,000 to 511,000 m 3 year −1 of wood assortments was also estimated. Firewood was the most impacted product (from 421,408 to 407,915 m 3 year −1 of reduction) but the structure of Italian forest chain and the high www.nature.com/scientificreports/ added value for the other timber outputs (e.g., roundwood) suggested a potential negative cascade effect on the whole forest chain and on ancillary activities.

Conclusions
This work is one of the first fine-scale combined models to quantify the economic impacts of O 3 at national level. Results highlighted that the reduction of forest area with active management is limited to the most severely O 3 polluted areas, even though significant negative effects on timber production occur all across Italy. Consequences on other forest ecosystem services and socio-economic deterioration of the forest chain, such as occupational consequences and cascade effect on satellite activities, should be evaluated. This aspect indicates that the elevated economic impact (loss of capital and annual values of forest) here presented, is still an underestimation of the total losses. The open source software facilitates replicability as well as sensitivity analysis. Spatial analysis can take into account local peculiarities, thus helping to improve models and results, and finally resulting into guidelines to cope with the potential negative impacts of O 3 pollution on forests and the forestry sector. This GIS-based application is thus a valuable tool to quantify and localize potential negative impacts. Ozone damages can interfere with the vitality of species of plant communities, as well as that of the animals, fungi, bacteria and insects that live in close association with plants or in nearby soils 34 . Changes induced by O 3 impact on many ecological processes, affecting ecosystem services, flows, goods and values. Further activities in the definition of ecosystem-scale models suitable to extrapolate effects of O 3 on productivity of trees and entire ecosystems might be addressed to economically quantify also the loss of other ecosystem services such as biodiversity, resource allocation and/ or seed production.

Methodology
Study area and forest data. Data of biomass availability were obtained from the pan-European map of forest biomass increment 35 . The map was validated using information from the most recent and free georeferenced Italian National Forest Inventory, i.e. for the year 2005 36 . In Italy, the total forested area was 10,467,533 ha in 2005. There was a certain variability of forest cover (in percentage) as well as vegetation, geomorphologic, logistics and socio-economic characteristics among regions (NUTS-2 level). Private forests accounted for 63.5% of the total, and coppice was the prevalent forest management (54.0%). The most common forest typology was broadleaved species (83.6% of total area), mainly oaks (23.9% of broadleaved, mostly Quercus robur and Q. cerris), as well as Fagus sylvatica (11.8%) and Castanea sativa (9%). Among the conifers (16.4% of total area), Picea abies prevailed (34.2%), in particular in the Alpine forests 36 .

Modelling ozone pollution.
Hourly O 3 concentrations and hourly meteorological data needed to calculate POD1 (i.e. solar radiation, air temperature, relative humidity, soil water content, wind speed) were simulated over the domain at 12 × 12 km of horizontal resolution by the WRF-CHIMERE modelling system, as described in 26 . The O 3 concentrations at 20-25 m above ground level (top of the canopy) provided by CHIMERE were used to calculate PODY. For PODY, a threshold Y of 1 nmol m −2 s −1 per leaf area as recommended by 6 for forest protection was applied, and computed as in 6,37 : where [O 3 ] is hourly O 3 concentrations (ppb), dt is time step (1 h), SGS and EGS are the start and end date of the growing season computed as described in 38 , R b is the quasi-laminar resistance (s m −1 ), R c is the leaf surface resistance (s m −1 ), and g sto is the hourly value of stomatal conductance to O 3 (mmol O 3 m −2 PLA s −1 , where PLA is the Projected Leaf Area) computed as following: where g max is the maximum stomatal conductance to O 3 of a plant species expressed on a total leaf surface area (mmol O 3 m −2 PLA s −1 ). The maximum stomatal conductance (g max ) is experimentally obtained as average above the 90 th or 98 th percentile of g sto measurements under optimum environmental conditions for stomatal opening 6,14 . The functions f phen , f light , f temp , f VPD , and f SWC are the variation in g max with leaf age, photosynthetically flux density at the leaf surface (PPFD, μmol photons m −2 s −1 ), surface air temperature (T, °C), vapour pressure deficit (VPD, kPa), and volumetric soil water content (SWC, m 3 m −3 ), respectively. The function f min is the minimum stomatal conductance. These species-specific functions vary between 0 and 1, and are expressed as: www.nature.com/scientificreports/ where light a is an adimensional constant; PPFD is hourly photosynthetic photon flux density estimated through the solar radiation; T opt , T min , and T max represent the optimum, minimum, and maximum temperature for g sto , respectively; VPD min and VPD max are minimum and maximum VPD for g sto ; and WP and FC are the soil water content (SWC) at wilting point and field capacity, respectively 6 . WP and FC are constant and depend on the soil type obtained from a module included into the WRF-CHIMERE model. We assumed that f phen was 1 throughout the growing season (0 otherwise). Six types of parameterization were used based on the six dominant forest types identified by land cover and climate over Italy, i.e. Alpine (for which the Boreal parameterizations in 6 was used), Continental and Mediterranean with either Deciduous or Evergreen species (Table 5). After calculation of POD1, the dose-response functions specific per each forest type were applied in each 12 × 12 km grid point according to Table 6.
Modelling ozone-induced wood losses. The economic value of Italian forests was quantified through a spatial-based analysis centered on the r.green.biomassfor model 39 . The tool is available as open-source add-on in GRASS Geographic Information System (GIS) software (https ://grass .osgeo .org/grass 78/manua ls/addon s/r. green .bioma ssfor .html). The model allows a quantification of wood assortments as well as their economic value by a multistep procedure. First, the so-called "technical availability" of material is quantified considering logistic (distance from forest/main roads or landing site) and geomorphological conditions (slope and terrain roughness). The combination of the above variables allows defining forest production process in terms of organization of cutting, processing and extraction. The second step introduces the economic parameters to calculate revenues from sell of assortments and full costs. Potentially available material is finally quantified on forest surfaces with economic profitability of production process. In the original version of the model, profitability was expressed as positive stumpage value (difference between revenues and costs of final harvesting / thinning). Here-in order to evaluate a long-term impact-the positive forest expectation value (FEV) (capital value of bare land plus timber at year y) was considered as index of economic profitability. FEV corresponds to the present value of cashflows arising from both the land and the tree, in perpetuity 40 . Input data to run r.green.biomassfor in Italy were derived from 19 and are provided in Tables S1 and S2. The complete procedure was developed on raster basis with a resolution of 1 ha (squared pixel of 100 × 100 m). Wood losses due to O 3 were computed taking into account the reduction of forest increment. In the present work, two scenarios were developed to quantify FEV and-as  Table 5. Forest-type parameterization of DO3SE model according to 6 . g max , maximum stomatal conductance; f min minimum stomatal conductance; f light_a parameter determining the shape of the hyperbolic relationship of stomatal response to light (dimensionless); T max , T opt and T min are maximum, optimal and minimum temperature for calculating the function f temp that expresses the variation of g max with temperature; VPD min and VPD max are the vapor pressure deficit for attaining minimum and full stomatal aperture calculating the function f VPD that expresses the variation of g max with vapor pressure deficit. The parameters for the soil water content (f SWC ) and phenological functions (f phen ) are obtained by the WRF model and vary with the latitude.  Table 6. Forest-type dose-response functions to estimate total biomass losses (L) based on the cumulated stomatal ozone flux above a threshold of 1 mmol m −2 (POD1) according to 6 . www.nature.com/scientificreports/ a consequence-forest surface with positive economic value and wood production: (i) a hypothetical scenario without O 3 impact (WO) and (ii) a scenario based on POD1 limiting biomass production. The software r.green. biomassfor was then ran for both scenarios and results were reported at national and administrative (regions) level. FEV for each scenario s (FEVs) was calculated as in 41 based on the "future revenues" approach: where SV is stumpage value of final harvesting, T is stumpage value of intermediate thinning, t is rotation period, m is year of thinning, n is age of the forest, v and e are yearly income and cost, respectively, q = 1 + r with r interest rate, LEV is the land expectation value (bare soil) calculated as in 42 : FEVs were calculated at national and region level for both total (€) (Eq. 9) and average (€/ha) (Eq. 10) values. As POD1-based damage is accumulated over the growing season, total and average FEVs were annualized (aFEV) as expressed in Eq. (11) and Eq. (12), respectively.
where ω is ω-th administrative level, i is i-th pixel in the map, the expression i ∈ FEV > 0 ∧ ∀i ∈ ω represents forest surfaces with positive FEV for ω-th administrative level, and r is interest rate.
In each scenario s, the amount of timber potentially obtained was calculated on the basis of forest surface with positive FEV, forest type and partitioning (percentage) of biomass increment in each wood assortment. The analysis of economic losses was focused on forest surfaces where FEV > 0 in WO scenario because damage evaluation must compare ex-ante (WO scenario) and ex-post (POD1) situations. Therefore, in the hypothesis of timber products investigation, an economic damage exists if and only if forest area can be economically processed also in WO scenario. The main assortments considered as typical for Italy are roundwood, timber pole, paper mill and firewood 19 .
Long-term economic valuation (that is typical in forest sector) is very variable relative to the interest rate applied for capitalization. Therefore, a sensitivity analysis based on variation of the interest rate r was performed at national level. In other terms, economic metrics (outputs, see Table 1) were computed for three different interest rates (2%, 3% and 4%) that are appropriate in Italy on the basis of forest and forest owners' characteristics 41 . Sensitivity analysis facilitates depiction of potential range of results when input variables can cause high level of uncertainty in the outputs. Regional evaluation considered-for brevity-only r = 3% 43 for quantification of FEV as well as forest surface with positive FEV.
To link biomass from dose response functions to economic losses, the O 3 -induced biomass reductions (POD1 scenario) were attributed to the forest map by overlaying spatial locations. The forest map was created based on Corine Land Cover (CLC) polygons (Table S1). Each CLC forest polygon reports one forest category and was linked to one biogeographical region of Italy (i.e. alpine, continental or Mediterranean 44 ). The CLC polygons were associated, through specific alphanumeric rules, to the forest types applied to compute dose-response functions (alpine deciduous, alpine evergreen, continental deciduous, continental evergreen, Mediterranean deciduous, Mediterranean evergreen). In case of lack of spatial matching among CLC polygons and POD1 forest types, a geographical extension of POD1 information was performed by means of proximity analysis (i.e. Voronoi tassellation of POD1 centroids 45 ). The results were finally mapped by QGIS (https ://www.qgis.org/it/site/).