A Combined Measurement and Modelling Approach to Assess the Sustainability of Whole-Tree Harvesting—A Swedish Case Study

: The demand of renewable energy has increased the interest in whole-tree harvesting. The sustainability of whole-tree harvesting after clear-cutting, from an acidiﬁcation point of view, depends on two factors: the present acidiﬁcation status and the further loss of buffering capacity at harvesting. The aims of this study were to investigate the relationship between these two factors at 26 sites along an acidiﬁcation gradient in Sweden, to divide the sites into risk classes, and to examine the geographical distribution of them in order to provide policy-relevant insights. The present status was represented by the acid neutralizing capacity (ANC) in soil solution, and the loss of buffering capacity was represented by the estimated exceedance of critical biomass harvesting (CBH). The sites were divided into three risk classes combining ANC and exceedance of CBH. ANC and exceedance of CBH were negatively correlated, and most sites had either ANC < 0 and exceedance (high risk) or ANC > 0 and no exceedance (low risk). There was a geographical pattern, with the high risk class concentrated to southern Sweden, which was mainly explained by higher historical sulfur deposition and site productivity in the south. The risk classes can be used in the formulation of policies on whole-tree harvesting and wood ash recycling. whole-tree harvesting sustainability, based on present acidiﬁcation status neutralizing capacity (ANC) in soil solution) and exceedance of CBH.


Introduction
Emissions of sulfur (S) and nitrogen (N) from the burning of fossil fuels have led to acidic deposition, peaking in the second half of the 20th century [1,2] and resulting in acidified soils and waters in the forest ecosystems of Europe [3]. These emissions have been reduced considerably, but both measurements [4][5][6] and modelling [7] indicate that the recovery of forest soils, which is required for surface waters to fully recover, is slow compared to the emission reductions. Recovery may slow down more, or even be reversed, if branches, tops, and stumps are harvested, since harvesting is an acidifying process [8,9] which is reinforced if other parts than the stems are harvested [10][11][12][13]. Removal of branches and tops means a significantly increased removal of base cations, due to the higher base cation concentrations in those parts [14,15].
The demand for renewable fuel is expected to increase to replace fossil fuel and thus mitigate climate change [16]. In Sweden, the harvesting of branches and tops were notified on 35-41% of the notified final fellings during the five years between 2015 and 2019, whereas stump harvesting was very uncommon [17]. Estimations of the future biomass potential in Sweden have shown that the harvesting of branches and tops has the potential to increase [18]. It is therefore important to be able to assess the sustainability of whole-tree harvesting in regions with different acidification history and site properties.
The sustainability of whole-tree harvesting at a site, from an acidification point of view, depends on two factors: (1) how the site has been affected by previous acidification and (2) how whole-tree harvesting will affect the buffering capacity of the soil. The effects of anthropogenic acidification on soils are often described by soil solution chemistry below the root zone [4][5][6]19]. Acid neutralizing capacity (ANC) is an indicator often used to provide a quantitative estimate of the acidification status at a site as well as the quality of the water leaving the root zone [20]. However, it does not give any indications about the effect of whole-tree harvesting. The effect of whole-tree harvesting on the buffering capacity can be estimated using acidity budget calculations, e.g., through the concept of critical biomass harvesting (CBH) [13]. Such estimates can give indications of how the buffering capacity is affected by whole-tree harvesting, but since they only give information about the direction and a rate of the development, and not about the starting point, they do not give a complete assessment of the risks related to whole-tree harvesting. By combining ANC data based on soil solution chemistry measurements with acidity budget calculations, both previous acidification and the effect of whole-tree harvesting on the buffering capacity of the soil can be taken into account.
In this study we perform a refined mapping of the risks of acidification of forest soils caused by whole-tree harvesting, in this article defined as the harvesting of branches and tops but not stumps after clear-cutting, across a strong deposition gradient in Sweden. We achieve this by combining data on the present acidification status, represented by ANC in soil solution, with estimates of the potential effect of whole-tree harvesting on the buffering capacity, represented by exceedance of CBH, at 26 well investigated spruce sites. The aims were to (1) investigate the relationship between ANC in soil solution and exceedance of CBH, (2) identify risk classes for the sustainability of whole-tree harvesting, taking both factors into account, (3) investigate the geographical distribution of the risk classes and (4) discuss the results in relation to policies about whole-tree harvesting and wood ash recycling.

Study Sites and Measurements
The study was performed at 26 spruce sites from the Swedish Throughfall Monitoring Network (SWETHRO) [21], most of them active during all five years in the period 2014-2018 (Table 1). The 26 sites are geographically distributed across all three regions of Sweden-South, Central and North ( Figure 1a)-to cover the steeply decreasing gradient of atmospheric S deposition from southwest to the north (Figure 1b). The site density is highest in the south of Sweden where deposition is highest, which is reflected in the number of sites in the three regions (Table 1).
The SWETHRO sites are 30 by 30 m squared plots within managed forest stands, where atmospheric deposition and soil solution chemistry are measured continuously. Throughfall deposition, i.e., precipitation that passed through the canopies, is collected monthly using open buckets in the winter, and polyethylene bottles with funnels threaded into the lid in the summer. Ten throughfall collectors are placed on each site, in an "L" shape along two of the borders of the monitoring area. The water samples are merged to one composite sample at the end of each month for chemical analysis at the accredited laboratory at IVL Swedish Environmental Research Institute and analyzed as described in previous studies [5]. The amount of throughfall precipitation per hectare is estimated based on the volume of water and the known diameter of the open bucket/funnel, and the throughfall deposition is derived by multiplying the amount of throughfall precipitation with measured concentrations. Site index is directly linked to site productivity (m 3 ha −1 y −1 ) and is a measure of the optimal growth of a stand. The letter "G" means that the tree species is Norway spruce and the number is the projected tree height when the stand is 100 years old. 2 Field assessments in 2020. 3 From the ICP level II database (International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests, http://www.icp-forests.org (accessed on 19 February 2021)). 4 Assessment based on satellite data and information from a nearby wood-land key habitat. 5 From forest owner/forest manager/Forest Agency. 6  For some of the sites, there are also wet deposition measurements on an open field close by, which are used to separate between wet and dry deposition. Wet deposition is measured using a 1.5 m high collector made of polyvinylchloride. A polyethylene plastic bag is placed inside the cylinder and changed at every sampling occasion. In addition, a robust plastic ring with a smooth rim is attached to the top of the cylinder. A net made of polyolefin is placed at the opening of the cylinder during spring, summer and autumn, to prevent debris to enter the plastic bag. During 2016 and 2017, the wet deposition samplers were equipped with "birdrings", positioned horizontally around the openings of the sampler to prevent birds from sitting on the rim of the sampler.
To be able to quantify dry deposition for the chemical species that interact with the tree crowns, e.g., N, Ca, Mg, and K, string samplers are placed at the open field of some of the sites. The string samplers are placed under a transparent roof made of polycarbonate. The samplers consist of Teflon strings. They are sprayed with deionized water once a month, all year round, and the samples are sent to the laboratory for chemical analysis. During wintertime the string samplers sometimes have to be brought indoors before being sprayed with deionized water [23]. At each of the 26 sites, soil solution samples from a depth of 50 cm in the mineral soil are collected using 3-6 separate suction lysimeters with ceramic cups (P 80) ( Table 1). Samples are taken three times a year: before (February-May), during (May-September), and after (September-December) the vegetation period. At the time of the sampling, negative pressure is initiated so that the lysimeter can suck in water, and after two days the water is collected. For each occasion, the water samples from the lysimeters are merged into one composite sample for analysis. If there is no collectible water, the site is still not revisited until the next scheduled sampling occasion. Drier sites therefore have fewer measurements of soil water concentrations.
Starting in 2010, a soil sampling campaign was performed at the SWETHRO sites. Soil samples from 4-5 distinct soil horizons, down to 50 cm in the mineral soil, were collected and sent for chemical analysis of, e.g., total elemental content and grain size distribution. Cylinders (three per soil sample) were used to estimate soil bulk density. Stoniness and moisture were assessed in field at all sites.

Estimation of ANC in Soil Solution
Acid neutralizing capacity (ANC in meq L −1 ) in soil solution was calculated for each site and each sampling occasion as: where all element concentrations are given in meq L −1 . The median ANC for the available samples in 2014-2018 (Table 1) was calculated and used for analyses in this study.

Estimation of Exceedance of Critical Biomass Harvesting (CBH) 2.3.1. Concept
The concept critical biomass harvesting (CBH [13]) has been used in Swedish policies since 2018, when a new indicator "Acidification from forestry" was introduced in the Swedish Environmental Objective Framework. CBH builds on the same principle as critical loads (CL), which was a successful concept in reducing acidic deposition [19].
The basis for both CL and CBH calculations is the SMB formula (Equation (2), [24]) S dep + N dep + Cl dep + BC harv +Alk leach = BC dep + BC weath + N imm + N harv + N de (2) where dep = deposition, BC = base cations (Ca, Mg, Na and K), harv = net losses at harvesting, Alk= Alkalinity, leach = leaching, weath = weathering, imm = immobilization and de = denitrification (meq m −2 y −1 ). In CL calculations, the critical load of acidity, i.e., the maximum amount of acid deposition (S + N) that can be allowed without negative ecosystem effects, is estimated through Equation (3), where all factors except S dep and N dep have been moved to the right side of the equation and the Alk leach is the critical alkalinity leaching (Alk leach(crit) ).
To calculate CL at a site, all factors on the right side of the equation except Alk leach(crit) are quantified based on site data. A chemical criterion and a critical limit are defined, linking acid deposition to ecosystem effects. They are used to calculate Alk leach(crit) , which is a measure of the acidification of runoff water. Exceedance is then calculated according to Equation (4), where the actual deposition at the site is compared with the CL.
In a national study in Sweden, CBH was defined as the maximum biomass extraction that does not lead to an ANC value less than zero in the water leaving the root zone, and therefore does not lead to export of acidity [13]. Biomass harvesting contributes to the acidity flows by removing BC from the forest ecosystem. In the CBH calculations, the critical removal of BC through harvesting was quantified and compared with the actual biomass harvesting, based on Equation (2) and the same logics as in the CL calculations in Equations (3) and (4), but with a few simplifications. The chemical criterion used was ANC, and since the critical limit was set to 0, the Alk leach(crit) factor disappeared. Furthermore, the N cycling was greatly simplified based on three assumptions:

•
Only N that leaches is affecting acidity: NO 3 -N leaching is acidifying, one equivalent, and NH 4 -N leaching counteracts acidification, one equivalent, based on theories in [25]. • N leaching is not affected by whole-tree harvesting. • N accumulating in soil organic matter will not contribute to acidification in the future.
By implementing those assumptions in Equation (1), CBH could be calculated using Equation (5), and exceedance could be calculated using Equation (6).
In this study, the same formulas and assumptions as in the national CBH calculations were applied [13], and calculations were made for the upper 50 cm in the mineral soil, which is assumed to be the root zone for spruce. The calculations assume a steady state, meaning that all fluxes are given as annual averages over an entire forest rotation, varying according to site productivity. Therefore, the acidification history of a site is not taken into account, and the results only show a direction and a rate of the change.

The Weathering Model PROFILE
Whereas deposition, leaching, and removal of BC by harvesting can be based on site measurements, weathering rates were modelled using the PROFILE model [26]. PROFILE is a steady state soil geochemical model, in which weathering rates are calculated using transition state theory. Geochemical properties of the soil, e.g., mineral composition, mineral surface area, and soil moisture, are required as input, along with atmospheric deposition, climate, and biological parameters. An overview of the input data required for PROFILE and for calculations of CBH and its exceedance is given in Table 2, and a detailed description of the data and how they were processed is given in the following paragraphs.  3 Weath, Exc Estimated based on site index from sites (Table 1) Net harvest loss of N 3 Weath Estimated based on site index from sites (Table 1) 1 Swedish Meteorological and Hydrological Institute. 2 The A2M program [27]. 3 Net loss at harvesting of biomass (stems, branches, and tops).

Input Data
Deposition was required both for PROFILE modelling and for calculating CBH and its exceedance. The average deposition of S, Cl, N and BC for the years 2014-2018 was based on SWETHRO measurements (Appendix A). For S and Cl, where the canopy exchange was assumed to be very small, throughfall deposition was used. For N and the base cations Ca, Mg, and K, the canopy exchange is substantial, and thus the throughfall deposition could not be used as it is. For N, wet deposition was measured on all 26 sites. Dry deposition was estimated for ten SWETHRO sites in Sweden, based on measurements and calculations included in the surrogate surface method [23]. These results were generalized to all 26 sites, based on a method where the share of dry to total deposition of inorganic N (NO 3 -N + NH 4 -N) could be correlated with geographical position [23]. Wet deposition of Ca, Mg, and K for 2014-2018 was available for 14 of the 26 sites (Table 1). For the sites missing measurements, wet deposition was derived from interpolation using the inverse distance method. For this interpolation, the 14 sites, and additionally 12 sites from the SWETHRO network not included in this study but for which wet deposition data were available, were used. Dry deposition of Ca, Mg, and K was estimated based on the relationship between wet and dry deposition in an extensive modelling study based on data from 1998 [28,29], assuming that the relationship between wet and dry deposition has been constant over time. Na was assumed to come only from sea salt, and was calculated based on Cl deposition and the sea salt composition.
Climate data was required for the weathering modelling. Temperature, precipitation, and runoff for the years 2014-2018 were obtained from Swedish Meteorological and Hydrological Institute (SMHI) (Appendix A). Temperature and precipitation were extracted from the database PTHBV (data delivered in August 2020). Runoff was derived from SMHI:s water web, with data for 40,000 subcatchments in Sweden.
Leaching of N in the form of NH 4 -N and NO 3 -N was required for the CBH calculations, and was based on measurements of concentrations in soil solution at 50 cm depth and runoff (Table 3; Appendix A). Median concentrations from the years 2014-2018 were multiplied with runoff from the sites. For the results to be representative for a forest rotation, NO 3 -N after clearcutting was added. For that, an empirical function from was used, where NO 3 -N leaching in the clearcut phase was related to site productivity [30]. Table 3. Soil solution chemistry, medians for the period 2014-2018 (meq L −1 ). The sites are arranged from north to south. The last column "n" is the number of samples. Net harvest losses of N and BC by whole-tree harvesting were required inputs for weathering modelling with the PROFILE model, and net harvest losses of BC were also important in the calculations of exceedance. The net harvest losses were calculated based on site productivity from the sites, derived from site index ( Table 1). The site productivity is the growth of a stand under optimal conditions, and to imitate real conditions the site productivity was reduced by 20%, in accordance with earlier mass balance studies [13,31]. All stems and 60% of the branches, accompanied by 75% of their needles, were assumed to be harvested, in accordance with a scenario from the Swedish Forest Agency [32]. The calculations were performed in the same way as in a national CBH study [13], and the methodology along with densities and nutrient concentrations in different tree parts are given there.

Site
Data on soil properties were important inputs to the weathering model PROFILE. Soil data in this study came from soil samplings at the sites, performed between 2010 and 2016. One of the required inputs was mineralogy for all soil layers, which was not measured, but can be calculated from measured total elemental content, quantified by plasma emission spectrometry analysis (ICP-AES). The total elemental content was used in the A2M program [27], model version 1.41, to calculate all possible mineral modes given a set of allowed minerals. However, for the organic upper layer the total elemental content from the second layer was used since the organic material in the organic layer substantially affects the total elemental content. For each layer at each site, the arithmetic means from A2M were then selected to represent the mineralogy at each layer. In absence of any other information, this was regarded as the best solution as it represents the barycentre of the solution polyhedron spanned by the extreme modes [27]. Mineralogies used for the sites are presented in Appendix B.
Soil bulk density was derived by weighing the dried samples collected with cylinders with known volume. Specific surface area was estimated based on grain size distribution and an empirically derived formula [33]. Stoniness (volume fraction of stones and boulders) and moisture were estimated based on visual assessments on the sites. Moisture was translated from a moisture class to m 3 water per m 3 soil [22]. Soil input data are described in Appendices C and D.

Risk Classification
As a basis for the risk classification, factor two, exceedance of CBH was plotted against factor one, ANC in soil solution. The sites were then grouped into three classes, based on their position in the graph (Table 4). Sites where both factors indicated risks related to acidification, i.e., sites where ANC was negative and whole-tree harvesting led to an exceedance of CBH, were placed in Class 1. Sites for which one of the factors indicated a risk but not the other, i.e., either negative ANC and no exceedance, or positive ANC and exceedance, were placed in Class 2. Finally, sites where none of the factors indicated risks related to acidification, i.e., sites for which ANC was positive and whole-tree harvesting did not lead to an exceedance of CBH, were placed in Class 3. Table 4. Risk classes for whole-tree harvesting sustainability, based on present acidification status (acid neutralizing capacity (ANC) in soil solution) and exceedance of CBH.

Class
Outcome for the Two Risk Factors Risk Related to Acidification from Whole-Tree Harvesting 1 ANC < 0 and Exc CBH > 0 High 2 ANC < 0 and Exc CBH < 0 or ANC > 0 and Exc CBH > 0 Medium 3 ANC > 0 and Exc CBH < 0 Low

Present Soil Status: ANC in Soil Water (Factor 1)
ANC in soil solution spanned from −0.29 meq L −1 to 0.39 meq L −1 , along a geographical gradient with generally higher ANC towards the north, although several deviations occurred at the regional scale ( Figure 2, Table 3). In the northern region, ANC was clearly positive at all five sites. In the central region, ANC was positive at two sites and negative at one. In the southern region only four of the eighteen sites showed a positive ANC.

Effect of Whole-Tree Harvesting: Exceedance of CBH (Factor 2)
The CBH calculations showed an exceedance of critical harvesting for 21 of the 26 sites ( Figure 3, Table 5). Just as for ANC in soil solution, there was a clear geographical gradient, with a generally decreasing exceedance towards the north. In the northern region, there was no exceedance for two of the five sites, and a low exceedance (up to 25 meq m −2 y −1 ) for the others. In the central and southern regions, the exceedance was above 25 meq m −2 y −1 for most of the sites, and at three of the sites in southernmost Sweden the exceedance was over 100 meq m −2 y −1 . In the central and southern regions, three sites stood out with negative exceedance.

Risk Classification Based on the Two Factors
A comparison between ANC in soil solution and the exceedance of CBH for the 26 sites showed that exceedance was negatively and significantly (p = 0.00025) correlated with ANC in soil solution (Figure 4a). A slightly weaker correlation was also found between ANC and CBH (p = 0.0046). The risk classes (Table 4) were unevenly distributed in the country (Figure 4b, Table 6). Risk Class 1 was highly dominating in the southern region, whereas the sites in the northern region belonged to either Class 2 (three sites) or Class 3 (two sites). In the central region, with only three sites, all three classes were represented.

Geographical Distribution of Risk Classes
The geographical pattern, with Risk Class 1 highly dominating in the southern region and Risk Class 2 and 3 dominating in the other regions, could be explained by the S deposition gradient with decreasing deposition towards the north (Figure 1b), and the site productivity gradient going in the same direction (represented by site index in Table 1). The historical S deposition, which has been higher than today but always followed a similar gradient [2], could be regarded as the main explanatory factor for the gradient for the first factor in the risk classification: ANC in soil solution. S deposition causes BC to leach from the soils [3], and the effect could be seen in the relationship between S and the base cations Ca, Mg, and K in soil solution, with generally more S compared to BC towards the south (Table 3). For the second factor, exceedance of CBH, site productivity was the most important reason for the gradient. Higher site productivity in the south leads to more removal of base cations at whole-tree harvesting (Table 5), which increases the risk of exceedance of CBH. This is consistent with the conclusions from a study at nine forest sites with different site productivity in Sweden and Scotland where fast growing forests led to more acidification than slow growing forests [8].
The sites in Risk Class 1 had all been exposed to high S deposition and all had high site productivity. The two sites with the lowest ANC, M16 and L18, situated in the southernmost part of Sweden, were the only sites with elevated NO 3 -N concentrations in soil solution. This had a high impact both on ANC in soil solution and exceedance of CBH, and highlights the considerable potential of acidification from N deposition in areas where the N retention capacity is being exceeded [25,34].
Sites in Risk Class 3 appeared in the whole country, except for in the far south. Three of the sites, BD06, S05, and G22, had high weathering rates in comparison with most of the other sites (Table 5), which can partly explain both the high ANC and the negative exceedance. The site productivities for BD06 and S05 were low, which was also the case for the fourth site, AC04, and low site productivities were an important contributing factor for the negative exceedance at those sites.
Additionally, the sites in Risk Class 2 were spread over the country. Several of the sites had values close to 0 for ANC and/or for exceedance of CBH, and for those sites, uncertainties have a great impact in the risk classification. Seven of the eight sites in Class 2 had positive ANC, but the CBH calculations showed exceedance at whole-tree harvesting. Three of these sites, in the northern parts of the southern region, E21, R09, and P95, were situated outside the area with the highest historical S deposition (Figure 1b), but still in the part of southern Sweden where site productivities are generally high, which could explain the combination of positive ANC and exceedance of CBH. In northern Sweden, the exceedances were generally close to 0, and small differences in weathering rates, deposition or site productivity determined if the estimated exceedance was positive or negative. The site AC34 had an exceedance just above 0, but very high ANC. The high ANC indicates that there are sources of base cations not accounted for in the CBH calculations, e.g., that weathering rates are underestimated.
One of the eight sites in Class 2, F18, had a slightly negative ANC, but CBH was not exceeded at whole-tree harvesting. Negative exceedance at a site in this part of Sweden, where most sites show exceedance, could be explained by a combination of higher estimated weathering rates and a somewhat lower site productivity than at most surrounding sites.
The low density of sites in the central and northern regions in this study limited the possibilities to perform geographical generalizations. However, earlier studies of ANC in soil solution [21,35] and exceedance of CBH [13] at higher resolution can contribute to the interpretation of the results in this study. The ANC gradient over Sweden has been presented for 2006-2008 [21] and 2017-2019 [35], where ANC on all SWETHRO sites active A high resolution national mapping of exceedance of CBH [13] showed that CBH was generally exceeded in spruce forests in Sweden at whole-tree harvesting, except for the inner part of northern Sweden, but it also showed that there were several exceptions for some areas in central Sweden and also a few areas in southern Sweden where there was no exceedance. This is consistent with the results from this study.

Policy Implications
At sites in Risk Class 1, from which the water leaving the root zone has no buffering capacity (ANC < 0) and where whole-tree harvesting leads to an exceedance of CBH, whole-tree harvesting is not sustainable from an acidification point of view. The sites in Risk Class 1 were situated in the southern part of Sweden, where recovery can be expected to be slow or non-existing even when whole-tree harvesting is not applied [5]. Whole-tree harvesting will lead to further loss of buffering capacity, which may inhibit recovery from acidification and thereby exacerbate the acidification status of the forest site. Wood ash recycling means that nutrients are returned to the soil, but there is a risk that the lost nutrients are only partly compensated for and/or that there is a time lag before the effect of wood ash recycling appears. Based on this, our assessment is that whole-tree harvesting at sites belonging to Risk Class 1 is not compatible with the Swedish environmental objective about acidification, even if the removal of base cations is compensated for by wood ash recycling.
For Risk Class 2, ANC and/or exceedance of CBH was often close to 0. Due to uncertainties in measurements and calculations, the class affiliation in itself was more uncertain than for Class 1. Sites that have a positive ANC, but an exceedance of CBH at whole-tree harvesting, have a better starting point than those in Risk Class 1, but the exceedance indicates that whole-tree harvesting is not sustainable in the long term; i.e., it will lead to loss of buffering capacity and a decrease in ANC. If ANC is decreasing at sites within this ANC interval, pH can be expected to be substantially affected. This was highlighted in a study on the effects of whole-tree harvesting and wood ash recycling on ANC in surface waters [36]. By compensating for the losses through wood ash recycling, the loss of buffering capacity caused by whole-tree harvesting can be inhibited. Sites in Risk Class 2 that have a negative ANC but no exceedance of CBH at whole-tree harvesting have the potential for recovery, but whole-tree harvesting will slow it down. We suggest that whole-tree harvesting at sites belonging to Risk Class 2 should be accompanied with wood-ash recycling.
In Risk Class 3, ANC was positive and the CBH was not exceeded at whole-tree harvesting. Our assessment is that the risk of negative effects of whole-tree harvesting on the acidification status is small, and that wood-ash recycling is not necessary at those sites.
In Risk Classes 1 and 2, where runoff water is already acidified or where acidification is in progress, there is a risk that nutritional imbalances will arise, which has been indicated in European studies [34,37]. This can, in turn, have a negative effect on tree vitality and tree growth [38,39]. Thus, although the risk assessment in this study focuses on the quality of the runoff water leaving the root zone, with potential effects on surface water, it is also highly relevant for soil acidification and tree nutrition. This study focuses on acidification, but there are other environmental aspects that can place restrictions on sustainable biomass harvesting, e.g., biodiversity [40]. These aspects should be taken into account in an overall risk assessment for sustainable biomass harvesting.

Conclusions
We investigated two factors for the assessment of the sustainability of whole-tree harvesting after clear-cutting: present acidification status represented by ANC in soil solution and exceedance of CBH. The two factors were negatively correlated, and the three risk classes, based on the two factors, showed a clear geographical pattern, although some deviations occurred. The gradients of historical S deposition and the site productivity were identified as the main causes of the geographical pattern. The high potential impact of acidification from N was highlighted at the two most extreme sites in the risk classification, where the N retention capacity was exceeded and nitrification caused acidification. The impact of soil weathering could be seen at single sites with high weathering rates, which were assigned in the low risk class, although historical S deposition and site productivity were high. The results strengthen the Swedish recommendations that wood ash recycling should be applied after whole-tree harvesting, except in the areas with low site productivity corresponding mainly to the inner part of northern Sweden. However, by including the present acidification status, we have also identified sites in southern Sweden where we assess that whole-tree harvesting is not compatible with the Swedish acidification objective.

Informed Consent Statement: Not applicable.
Data Availability Statement: Input data for weathering modelling and CBH calculations are available in Appendix A (climate and deposition), Appendix B (mineralogy), and Appendices C and D and Table 1 (soil properties). Site productivity is available in Table 1, and nutrient concentrations and tree densities used in the calculations of removal of nutrients are published in papers referred to in the methods section. ANC in soil solution used in the risk classification is presented in Table 3. Deposition and soil solution chemistry presented in this paper is average or median values. The original data can be downloaded from the IVL website, www.ivl.se (accessed on 19 February 2021).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix C
For all sites used in this study except BD06, clay and silt were not separated in the grain size analysis. Other sites in the SWETHRO database and in other databases had separate analyses of clay and silt, and the average fraction from these soil samples was used for the sites with missing data (Table A4).  [21]. 3 [42].
Appendix D Table A5. Layer thickness and density.