Weathering rates in Swedish forest soils

Soil and water acidification was internationally recognised as a severe environmental problem in the late 1960s. The interest in establishing “critical loads” led to a peak in weathering research in the 1980s and 1990s, since base cation weathering is the long-term counterbalance to acidification pressure. Assessments of weathering rates and associated uncertainties have recently become an area of renewed research interest, this time due to demand for forest residues to provide renewable bioenergy. Increased demand for forest fuels increases the risk of depleting the soils of base cations produced in situ by weathering. This is the background to the research programme Quantifying Weathering Rates for Sustainable Forestry (QWARTS), which ran from 2012 to 2019. The programme involved research groups working at different scales, from laboratory experiments to modelling. The aims of this study were to (1) investigate the variation in published weathering rates of base cations from different approaches in Sweden, with consideration of the key uncertainties for each method; (2) assess the robustness of the results in relation to sustainable forestry; and (3) discuss the results in relation to new insights from the QWARTS programme and propose ways to further reduce uncertainties. In the study we found that the variation in estimated weathering rates at single-site level was large, but still most sites could be placed reliably in broader classes of weathering rates. At the regional level, the results from the different approaches were in general agreement. Comparisons with base cation losses after stem-only and whole-tree harvesting showed sites where whole-tree harvesting was clearly not sustainable and other sites where variation in weathering rates from different approaches obscured the overall balance. Clear imbalances appeared mainly after whole-tree harvesting in spruce forests in southern and central Sweden. Based on the research findings in the QWARTS programme, it was concluded that the PROFILE/ForSAFE family of models provides the most important fundamental understanding of the contribution of weathering to long-term availability of base cations to support forest growth. However, these approaches should be continually assessed against other approaches. Uncertainties in the model approaches can be further reduced, mainly by finding ways to reduce uncertainties in input data on soil texture and associated hydrological parameters but also by developing the models, e.g. to better represent biological feedbacks under the influence of climate change.

1 Introduction

Acidification of soils and water, caused by long-range transport of acidic compounds, was recognised as an environmental problem in Europe in the late 1960s (Odén, 1968). In subsequent decades, extensive research examined processes that acidify and counteract acidification (Reuss and Johnson, 1986). Two key research programmes were the Surface Water Acidification Programme (1985–1990, Mason, 1990) funded by the UK (GBP 5 million) and the National Acid Precipitation Assessment Program (1980–1990, Irving, 1991) funded by the US government (USD 17 million). At the end of the 1980s, the critical load concept was developed as an effect-based approach for emission reductions (Nilsson and Grennfelt, 1988) and served as a link between science and policy within the framework of the UNECE–CLRTAP (the United Nations Economic Commission for Europe – Convention on Long-Range Transport of Air Pollutants) (Lidskog and Sundqvist, 2002). A critical load is defined as “a quantitative estimate of an exposure to one or more pollutants below which significant harmful effects on specified elements do not occur according to present knowledge” (UNECE, 1994). To calculate critical loads of acidity and their exceedance, mass balance calculations of acidity are used together with a critical limit for a chemical criterion, defining the maximum acidity of soil/runoff water that can be allowed without a risk of negative effects on a chosen biological indicator (Sverdrup and de Vries, 1994).

Figure 1 Base cation (BC) mass balance for the rooting zone of a well-drained forest soil. ΔSoil BC is the sum of the net change in BC in soil solution, the net change in soil exchangeable BC and the net change in soil organic material BC.

Estimates of base cation weathering (Ca, Mg, K, Na) play a key role in all kinds of mass balance calculations related to acidity, e.g. the mentioned critical load calculations, as weathering is an important long-term natural source of base cations and as such a sink of acidity. The net accumulation or depletion of soil base cations in the soil is the result of the mass balance between inputs (atmospheric deposition and weathering) and outputs (losses through leaching and harvesting) of base cations (Fig. 1). Atmospheric deposition depends on external factors and can vary in time (Hedin et al., 1994). The leaching term is directly dependent on the mass balance as it mirrors the aqueous pool of base cations in the soil. Harvest losses are predefined following forest management. Finally, weathering is the long-term source of base cations, depending largely on soil mineral content and soil texture. In areas that have been covered by ice, like Scandinavia, the soil properties at a specific site depend on the parent material from which the glacial till originates. When trees and forest soils are in focus, the weathering term refers to the weathering products in the rooting zone, i.e. the base cations available for trees.

Due to its central role in mass balance and acidity calculations, weathering was studied extensively during the 1980s and 1990s, to enable accurate weathering quantifications. Weathering rates are, however, difficult to quantify through direct measurements in the field due to the complexity of base cation dynamics in soil. There are several different pools of base cations in soil and also several different flows, e.g. decomposition, uptake, ion exchange and weathering, and it is difficult to distinguish between these different sources and sinks (Rosenstock et al., 2019) and to define the pools accurately (van der Heijden et al., 2018). Therefore, a number of indirect methods to quantify weathering rates have been developed: process-based modelling (Sverdrup and Warfvinge, 1993); soil measurements where the depletion of weathering products in different soil layers is determined in order to assess average weathering rates since soil formation (Olsson et al., 1993); and budget calculations where the flows in the mass balance, except weathering, are measured (Lundström, 1990; Jacks and Åberg, 1987; Wickman and Jacks, 1991; Sverdrup et al., 1998).

The political and scientific agreement on the critical load concept as a basis for managing acidic deposition was a major factor in subsequent policy success on limiting acidifying emissions (Lidskog and Sundqvist, 2002). Major acidic deposition reductions occurred, and some recovery of soils and surface waters has been noted by the monitoring operations put in place by the UNECE–CLRTAP (Graf Pannatier et al., 2011; Pihl Karlsson et al., 2011; Akselsson et al., 2013). The uncertainties in weathering rates were, during the times of high deposition, less important. Therefore, the interest waned in further weathering research that might revise these weathering estimates.

As the severity of climate change became fully recognised, policies for mitigation of climate change led to increased demand for renewable fuels, thereby increasing the pressure on forests. Whole-tree harvesting, here defined as harvesting of stems and branches, was seen as an important source of renewable fuel. Since 2000 in Sweden, the proportion of clearcuts involving whole-tree harvesting has increased from around 15 % to 25 %–35 %, according to statistics from the Swedish Forest Agency, except for 2014–2016, when the proportion was temporarily between 15 % and 25 % due to lower energy prices. Demand is likely to increase in the future (Börjesson et al., 2017). Harvesting of branches and the nutrient-rich needles also means a substantially increased removal of base cations compared to conventional stem harvesting, which may counteract the recovery from acidification. Thus, whereas the weathering rates in the past mainly were compared with deposition, which at that time (1970–1980) was much greater than estimated weathering and critical loads of acidic deposition, the interesting comparison in today's policy context is between base cation weathering rates and base cation losses through harvesting, to assess sustainable forest management. Accordingly, mass balance calculations, with weathering and deposition as inputs and harvest losses and losses through leaching as outputs (Akselsson et al., 2007) (Fig. 1), as well as simplified calculations, where weathering rates are compared with base cation losses through harvesting (Olsson et al., 1993), have been used during the last decades for forest sustainability assessments. The conclusions from these studies were that base cations may be depleted in spruce forests in certain regions if whole-tree harvesting is applied.

The prospect of increasing demand for forest biomass, and particularly concerns about the effects of whole-tree harvesting on nutrient sustainability, renewed the interest in weathering in Scandinavia for new forest policy issues. The accuracy of the weathering calculations and the conclusions about the long-term sustainability of forests were questioned by Klaminder et al. (2011), who compared estimations of weathering rates for Ca and K from different approaches at a site in northern Sweden. These estimates differed widely, and the study concluded that nutrient budgets, based on calculations including weathering rates, are too uncertain to be useful in shaping forest policies regarding harvest practices. Futter et al. (2012) suggested that some of this variation is due to differences in boundary conditions, for example the depth to which weathering had been calculated. They examined weathering estimates from 82 sites, with up to eight different weathering estimates per site, and found considerable variability in weathering rates estimated with different methods, often with results differing on the same site by several hundred percent. They identified uncertainties in input data as the largest contributor to the variability, but differences in the soil depth for which weathering was calculated also contributed, in the same way as in Klaminder et al. (2011). Futter et al. (2012) concluded that the uncertainties are large and that at least three independent methods should be used when making management decisions.

In 2012, a SEK 25 million research programme, Quantifying Weathering Rates for Sustainable Forestry (QWARTS), was started in Sweden. The programme, which focused on weathering rates for the base cations Ca, Mg, Na and K in Sweden, included approaches covering the whole spectra from laboratory-scale experiments, through plot- and catchment-scale experiments in the field, to extensive weathering modelling. Different approaches were compared at different scales and were in some cases refined (Stendahl et al., 2013; Casetou-Gustafson et al., 2018, 2019a, b; Akselsson et al., 2016; Belyazid et al., 2019; Kronnäs et al., 2019; Erlandsson et al., 2016). The input data were also examined for uncertainties relating to the generalisations made when estimating normative mineralogy based on total chemical analyses (Casetou-Gustafson et al., 2018, 2019a). Two other potential sources of uncertainties that have been explored, but that are still widely discussed, were revisited: (1) the role of biological weathering that might generate weathering not included in the current generation of biogeochemical models (Banfield et al., 1999; Finlay et al., 2009, 2019) and (2) model simplifications related to base cation exchange and aluminium complexation (Gustafsson et al., 2018; van der Heijden et al., 2018). Furthermore, the weathering kinetics used in models were revisited (Sverdrup et al., 2019), and weathering rates representing not only the rooting zone but the full catchment scale were studied to assess the export of weathering products to surface waters (Ameli et al., 2017; Erlandsson Lampa et al., 2019).

The aims of this study were to (1) investigate the variation in published weathering rates from different approaches in Sweden, with consideration of the key uncertainties for each method; (2) assess the robustness of the results in relation to sustainable forestry; and (3) discuss the results in relation to new insights from the QWARTS programme and propose ways to further reduce uncertainties. While weathering is important for understanding the acidification of both soils and surface waters, this paper focuses on soils, and specifically the rooting zone (approximately 50 cm).

Table 1 Sites included in the study, where at least two well-documented approaches for estimating weathering rates have been applied to the same depth.

^a Including the O layer.
^b Not including the O layer.
^c For MAGIC the weathering rate was calculated to 0.6 m, including the O layer.

Figure 2 Well-investigated sites where weathering rates have been calculated for the same soil depth with at least two different approaches. See also Table 1.

2 Methods

Weathering rates from different approaches were compared on a site level on well-investigated sites, and on a regional level on a larger number of sites but with more generalised input data. For the comparison of weathering rates on single sites, weathering estimates of base cations (Ca, Mg, K and Na) from Swedish forest sites, where at least two well-described approaches had been applied to the same soil depth, were compiled from literature and compared (Table 1; Fig. 2). The 23 sites found were located on till soils, with a mineralogical composition characterised by granitic and gneissic bedrock, i.e. with mainly quartz, orthoclase and plagioclase, and small amounts of mafic minerals such as amphibole and epidote (Table 2). A total of 13 of the sites were taken from Stendahl et al. (2013). Of the originally 16 sites in that study, 3 were excluded, since site conditions were not appropriate for using the depletion method (two of the sites) or PROFILE (one of the sites) (Stendahl et al., 2013). For each of the 23 sites found, medians of the different approaches were calculated along with maximum deviation from the median (Table 3). The different approaches are described in Sect. 2.1–2.4.

Table 2 Mineral content in soil at 50 cm depth (weight %) used as input for modelling weathering rates for the different sites: Qz (quartz), Or (orthoclase), Pl (plagioclase), Am (amphibole), Ep (epidote). Bi (biotite), Ap (Apatite), Mu (Muscovite), Ch (chlorite), Il (illite), Ve (vermiculite), and Hy (hydrobiotite). Minerals occurring in very small amounts, with minor effect of weathering rates, are not included in the table. Input data were not found for Gårdsjön C. The mineralogy has in some cases been slightly simplified to make it fit in one table. For detailed mineralogy, see original references (Table 1).

^a Mineralogy has been calculated with multiple models, but only the X-ray powder diffraction (XRPD) results are presented here.
^b Muscovite and illite could not be separated with the XRPD method.

Table 3 Weathering rates ($\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and statistics at the sites described in Table 1. Values in parentheses after the intervals represent the middle of the intervals.

^a Maximum difference from median (+ indicates that the maximum difference is higher than the median, – indicates that it is lower than the median, and ± indicates that it is as big a difference to the maximum and minimum values).

For the regional-level comparison, three published approaches for calculating weathering rates of the nutrient base cations (Ca, Mg, K) on a regional level in Sweden were revisited, harmonised and compared: PROFILE, ForSAFE and the depletion method combined with the total analysis regression approach (Sect. 2.1–2.3). Na was not included, since there were no estimates for Na from the latter approach. Weathering rates were compared in 346 sites in the SAFE database (Alveteg, 2004), which were the sites of the in total 640 sites in the database for which all methods could be applied successfully and where data on stones and boulders were available (Stendahl et al., 2009). For a regional comparison, the sites were divided into seven climate regions, simplified from 19 weather forecast regions used by the Swedish Meteorological and Hydrological Institute (SMHI) (Fig. 4). One of the regions, northwestern Sweden, only contained one site and was therefore excluded from the analysis. The regional approaches are further described in Sect. 2.1 and 2.3.

To put the weathering estimates in a sustainability perspective, simplified base cation mass balance calculations were performed for the single sites, where weathering rates, the most important natural long-term source of base cations, were compared with harvest losses of base cations, which is the one of the outputs that can be anthropogenically controlled (Fig. 1; Sect. 2.5). Finally, the results from the studies in QWARTS, on biological weathering and on the representation of base cation exchange and aluminium complexation in the models, were synthesised; main uncertainties were highlighted; and ways to reduce them were proposed.

2.1 Modelling based on weathering kinetics

Due to the difficulties in measuring field weathering rates, weathering kinetics have been frequently studied in laboratory environments (Brantley et al., 2008). Mechanistic modelling of weathering rates, based on laboratory-determined weathering kinetics, is one of the most widely used approaches for estimating field weathering rates (Warfvinge and Sverdrup, 1992; Godderis et al., 2006; Maher et al., 2009). The PROFILE model (Warfvinge and Sverdrup, 1992) is a steady-state soil chemistry model where weathering is derived from the breakdown of minerals, based on process-oriented descriptions of chemical weathering and solution equilibrium reactions. Weathering rates are calculated for different layers, with different soil properties, using transition state theory and the geochemical properties of the soil system, such as soil wetness, temperature, mineral surface area and mineral composition, and organic acid concentrations. Deposition of sulfur, nitrogen and base cations, as well as net losses of base cations and nitrogen through harvesting, is used as input for modelling of pH and base cation concentrations in soil water, which is required for the weathering modelling. Weathering rates are calculated for each mineral separately, using rate coefficients from laboratory studies for four reactions: with H⁺, water, CO₂ and dissolved organic carbon (DOC) (Sverdrup and Warfvinge, 1993). PROFILE has been used widely for estimating weathering in Europe (Akselsson et al., 2004, 2016; Stendahl et al., 2013; Holmqvist et al., 2003; Koptsik et al., 1999; Langan et al., 1995), the US (Phelan et al., 2014) and Asia (Fumoto et al., 2001). At all 23 sites in the single-site comparison in this paper, weathering estimates from PROFILE were available (Table 3).

The weathering submodel in PROFILE was later built in to the dynamic version SAFE (Alveteg et al., 1995), which was mainly used for acidification assessments, but has also been used for studying the dynamics of weathering rates (Warfvinge et al., 1995). Later, the SAFE model was coupled with the tree growth model PnET (Aber and Federer, 1992), the decomposition model DECOMP (Walse et al., 1998; Wallman et al., 2004) and the hydrological model PULSE (Lindström and Gardelin, 1992), resulting in the forest ecosystem model ForSAFE (Wallman et al., 2005; Belyazid et al., 2006). The ForSAFE model simulates the integrated biogeochemical processes of a forest ecosystem. It covers the processes of photosynthesis, allocation and growth, water and nutrient uptake, litterfall, organic matter decomposition and mineralisation, ion exchange, chemical speciation of and reactions between different elements, and hydrological transport. All process rates are internally regulated by microenvironmental conditions such as acidity, water availability, temperature and element concentrations. The model requires inputs of external drivers in the form of climate, atmospheric deposition and forest management, as well as inputs on the properties of the forest ecosystem, such as soil texture, mineralogy and tree species. ForSAFE is used for studying the effects of climate change, atmospheric deposition and forest management on tree growth, soil chemistry and runoff water quality. Although the weathering module is the same in PROFILE and ForSAFE, some differences can be expected since ForSAFE includes dynamics, which means that weathering is affected by other processes over time and that soil moisture, which is an input in PROFILE, is dynamically modelled in ForSAFE. In the single-site comparisons in this paper, weathering estimates from ForSAFE were available at two sites (Table 3).

For the regional comparison of weathering rates of nutrient base cations (Ca, Mg, K), regional runs from both PROFILE and ForSAFE where included. Regional PROFILE weathering estimations for Sweden were taken from Akselsson et al. (2016), where weathering rates for the upper 50 cm of the soil (including the organic layer) have been modelled based on data from 17 333 Swedish National Forest Inventory (NFI) sites (Fridman et al., 2014). Within QWARTS, weathering rates to 50 cm depth (including the organic layer) were modelled also with the ForSAFE model (Belyazid et al., 2019) on sites in the SAFE database. The SAFE database is a subset, consisting of 640 sites, of the NFI sites used for PROFILE modelling.

2.2 The depletion method

Another widely used approach for estimating weathering rates is the depletion method (e.g. Olsson et al., 1993; Starr et al., 1998; Stendahl et al., 2013). The method estimates historical weathering, i.e. the average weathering rate since the last deglaciation, of mobile (weatherable) elements, based on element concentrations in weathered soil horizons as compared to unweathered parent material. The method accounts for the general losses of soil material in a horizon by including an immobile (recalcitrant) element in the estimation. Concentrations of mobile elements will decrease as a result of weathering, while the immobile element will be enriched towards the soil surface. The concept has a long history (Marshall and Haseman, 1942), while the theoretical framework was later formalised by Brimhall and Dietrich (1987) and Brimhall et al. (1991). The most commonly used immobile element is zirconium, which is found in the resistant mineral zircon (ZrSiO₄) with negligible weathering (Hodson, 2002). The assumptions for the depletion method are as follows: (1) there is no weathering of the immobile element; (2) the soil pedon consists of homogeneous soil, where the deep soil constitutes the parent material; and (3) no weathering occurs beyond a certain depth. The average annual weathering rate is calculated from the soil age, i.e. the time since deglaciation or since the land rose from the sea due to glacio-isostatic uplift. The average rate may deviate from current levels depending on the variation in weathering rates over time (Taylor and Blum, 1995). Weathering rates from the depletion method were available for 18 of the 23 sites in the single-site comparison (Table 3).The depletion method has been used, in combination with the total analysis regression approach, for regional applications (see Sect. 2.3).

2.3 Total analysis regression approach

The total analysis regression approach is a derivative from the depletion method, which requires much less soil data than the depletion method. It is based on the fact that weathering rates of different elements have been found to correlate with total content of the elements in soil and temperature (Olsson and Melkerud, 1990). Based on weathering estimates from the depletion method, linear regressions containing total chemical contents for base cations in the C horizon (either separately or lumped together), and temperatures or temperature sums (i.e. daily mean temperature above a threshold value, summarised for the growing season), have been produced for a number of sites, and the regressions have then been applied to other sites (Sverdrup et al., 1998; Maxe, 1995; Olsson et al., 1993). In the single-site comparison, estimates based on the total analysis regression approach were available for three sites (Table 3).

On a regional level in Sweden, weathering rates for Ca, Mg and K have been calculated based on the depletion method in combination with the total analysis regression approach, as a basis for assessments of nutrient sustainability after whole-tree harvesting (Olsson et al., 1993). In the study from 1993, regressions between weathering rates calculated with the depletion method and different site factors were analysed on 11 sites. The strongest relationships were found between weathering rates of an element and the product of the concentration of the element in the C horizon and the temperature sum. In the present study, the regression functions in Olsson et al. (1993) were used to calculate weathering rates of Ca, Mg and K at the sites from the SAFE database that were modelled with both PROFILE and ForSAFE (see above), for a proper comparison. The temperature sum was calculated based on latitude and altitude according to Morén and Perttu (1994). Some of the calculations gave negative results for one of the base cations. This can be explained by the fact that the regressions are used for a new dataset, covering a broader range of temperature sums than the original dataset, which illustrates a limitation of the method. The estimations apply down to the weathering depth, which means different soil depths on different sites but normally between 40 and 70 cm in the mineral soil (Mats Olsson, personal communication, 2002).

2.4 Budget approach

Weathering rates can also be estimated through the budget approach (Paces, 1986; Lundström, 1990; Sverdrup and Warfvinge, 1991; Sverdrup et al., 1998). In the budget approach, sources and sinks of base cations are considered, and weathering is calculated as the difference between sinks and sources. However, one difficulty is to distinguish between weathering, changes in the exchangeable pool and net mineralisation, so a steady state is often assumed, and the weathering rates are calculated as leaching + net uptake − deposition. The budget approach has been applied at five of the sites in the site-level comparison (Table 3). In Svartberget in northern Sweden the simplified approach assuming steady state was used (Lundström, 1990). On two sites in Gårdsjön in southwestern Sweden, net mineralisation was estimated and considered in the weathering estimations, but changes in the exchangeable pool were disregarded (Sverdrup et al., 1998). Casetou-Gustafson et al. (2019b) estimated weathering rates based on measurements of atmospheric deposition, leaching, accumulation in biomass and changes in the soil exchangeable pool for control plots of long-term fertilisation experiments in young Norway spruce forests in Asa and Flakaliden, in southern and northern Sweden, respectively. When using weathering estimates from the budget approach, the assumptions used must be carefully evaluated (Sverdrup and Warfvinge, 1991; Rosenstock et al., 2019).

The budget approach can also be applied by using the MAGIC model (Cosby et al., 2001), which was developed to predict effects of acidic deposition on surface water acidification. Weathering is not mechanistically modelled as in the models described above; instead, weathering rates are calculated internally using mass balances (Maxe, 1995; Köhler et al., 2011). MAGIC handles input fluxes – atmospheric deposition and base cation weathering – and output fluxes – net uptake in biomass and runoff losses. These fluxes govern processes in the soil, e.g. cation exchange, with the pool of exchangeable base cations in the soils at the centre. When the fluxes change over time, it affects the chemical equilibria between soil and soil solution, which has an impact on surface water chemistry. Observed values of surface water and soil chemistry are used to calibrate the model. Weathering rates from MAGIC were available at two of the sites in the site-level comparison (Table 3).

Another way of applying the budget approach is to use the strontium (Sr) isotope ratio. In weathering estimations based on Sr isotope ratios, the difference in the ratio ⁸⁷Sr∕⁸⁶Sr in bedrock and in atmospheric deposition is used (Wickman and Jacks, 1991). Since soil water is a mixture of what comes from deposition and what comes from weathering, the weathering rate of Sr can be estimated. Ca and Sr follow each other closely in forests (Wickman and Jacks, 1991), so the weathering of Ca is assumed to be linearly correlated to the weathering of Sr in the calculations. This method was used at three sites in the site-level comparison (Table 3). The base cation∕Ca fraction was assumed to be constant in all three studies.

2.5 Assessment of forest sustainability

Simplified budget calculations based only on weathering rates and harvest losses (Olsson et al., 1993; Klaminder et al., 2011; Stendahl et al., 2013) were performed for a selection of the sites in Table 3. The criteria that had to be fulfilled for inclusion of a site were (1) availability of weathering rate assessments for the root zone which in Swedish forest soils is often defined as 0.5 m (Rosengren and Stjernquist, 2004) but in this study included depths down to 0.7 m and (2) access to data on site quality. Calculations were made for sites in spruce and pine forest.

The calculations of harvest losses were based on the site quality of the forest (average growth rate per year during a forest rotation at optimal conditions), reduced by 20 % to mimic actual conditions, and generalised densities and nutrient base cation concentrations in different tree parts. Two types of harvesting were considered, conventional stem-only harvesting and whole-tree harvesting, where, in addition to stem, tops and branches are removed for biofuel. It was assumed that, in whole-tree harvesting, all stems and 60 % of the branches were harvested and that 75 % of the needles were removed with the harvested branches. The methodology along with densities and base cation concentrations used is more thoroughly described in Akselsson et al. (2007).

Figure 3 Base cation weathering rates (sum of Ca, Mg, Na and K) for sites where different methods have been applied for the same depth on the same site (Fig. 2). The soil depth is around 0.50 m (with or without organic layer; see Table 1), except for a few cases where greater depths are given. The four spans to the right are intervals that were commonly used in the critical load work of CCE (Coordination Centre for Effects) within the UNECE Convention on Long-range Transboundary Air Pollution (LRTAP convention) (de Vries, 1994; Umweltsbundesamt, 1996). The intervals correspond to weathering rates for different parent material classes – acidic, intermediate and mafic – and different texture classes – coarse, medium (including the mix between medium and course material) and fine (including the mix between fine and medium material).

3 Results and discussion

3.1 Weathering rate comparisons at the site level

The single-site comparison enabled us to quantify the span of base cation weathering rates produced by the different approaches (Fig. 3, Table 3). The median weathering rates for the 23 sites spanned between 22 and 61 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 3). For two of the six sites where at least three approaches had been applied, the weathering rate spans were narrow: Gårdsjön B (49–62 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ and a maximum deviation from the median of +15 %) and Gårdsjön C (36–40 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ and a maximum deviation from the median of ±5 %). The span was somewhat wider in two of the other well-investigated sites: Stubbetorp (35–67 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ and a maximum deviation from the median of +56 %) and Flakaliden (34–61 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ and a maximum deviation from the median of +42 %). The weathering rate spans in Svartberget B and Asa, with four and three weathering estimates respectively, were remarkably wide: 31–85 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ and a maximum deviation from the median of +121 % in Svartberget B and 11–131 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ and a maximum deviation from the median of 254 % in Asa. The wide spans in Asa and Svartberget B were mainly due to substantially higher weathering rates according to the budget approach as compared with the other methods. Also, in Flakaliden, the span was expanded by the budget approach.

In Svartberget B, the weathering estimates from the budget approach can be expected to be overestimated for several reasons (Lundström, 1990). The method does not distinguish between weathering, base cation exchange and base cation release through decomposition. The measurements were carried out in the 1980s when the acidification process was taking place, leading to base cation release from the exchangeable pool, although this effect was much more pronounced in southern Sweden. Finally, dry deposition was not included in the calculations due to lack of data. In Flakaliden and Asa, base cation accumulation in biomass was the major sink in the mass balances. The high weathering rates produced by the mass balances, especially in Asa, were largely explained by the measured depletion of base cations in the soil being much lower than the accumulation in the young Norway spruce stands (Casetou-Gustafson et al., 2019b). The very high estimated weathering rates, especially in Asa, indicate that flows are described inadequately. Uptake occurring below the defined rooting zone could be one contributing factor (Casetou-Gustafson et al., 2019b). In Asa, very low weathering rates from the depletion method contributed to the wide span. The low weathering rates originated from a fairly flat Zr depth gradient in the soil, which indicates that the soil had probably been disturbed, so the necessary assumptions for the depletion method were not satisfied (Casetou-Gustafson et al., 2019b).

In Stubbetorp, the PROFILE estimates were substantially higher (67 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) than the estimates from the total analysis regression approach and MAGIC (35–43 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Maxe (1995) noted that the PROFILE weathering rate was higher than expected, given the soil properties in Stubbetorp, and argued that it may be due to unreasonably high specific surface area of the soil as input to PROFILE on the site. Specific surface area has been determined by BET analysis and could, according to Maxe (1995), be overestimated due to a large occurrence of Al and Fe precipitates.

Whereas weathering rates in the different sites in Gårdsjön (A, B and C) and Svartberget (A and B) were generally on the same level, except for the budget approach in Svartberget B as discussed above, the two sites in Risfallet (A and B) gave quite different weathering rates: 29–68 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in Risfallet A (0.5 m) and 25–29 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in Risfallet B (1 m). It could be expected that the weathering rates were higher in the deeper soil profile, but instead it was the other way around. The PROFILE-modelled weathering rate in Risfallet A is one of the highest of all sites reported in Stendahl et al. (2013), which can be explained by a relatively high clay content (7 %) and high soil bulk density. In contrast, Risfallet B (1 m) has a very low specific surface area (Sverdrup and Warfvinge, 1993), which can explain the low weathering rate produced by the model.

In five of the 17 cases where only two methods per site were applied, the maximum difference between the calculated median and the estimated weathering rates was less than ±10 % (Table 3). In the other end, five sites showed a corresponding difference between ±40 % and ±64 %. The results show that the width of the span varies substantially between sites, and to explain these differences the sites and the methods need to be studied in detail.

The weathering rates calculated with the depletion method for the 13 sites from Stendahl et al. (2013) were generally lower than the PROFILE-modelled weathering rates (Table 3). The two sites with largest difference, Skånes Värsjö and Kloten, with a maximum difference between the median and the estimated weathering rate of around ±60 %, distinguished themselves with very low rates estimated with the depletion method: 8 and 11 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. This is contrary to what was expected, since the weatherability of the soil is believed to decrease over time due to the depletion of more easily weathered minerals and formation of resistant coatings on the mineral surfaces (Taylor and Blum, 1995). The reasons for this discrepancy are not fully known, but one reason could be that the original till partly comprises already weathered till from previous glaciations (Stendahl et al., 2013). Moreover, it is likely that declining weatherability over time is less pronounced in these young glacial till profiles, where the easily weathered minerals remain in the profile. Furthermore, some drivers of weathering, such as forest growth, are more prominent today, which may overshadow long-term decline in soil weatherability. For the total sum of base cations there was a tendency towards higher modelled rates when the rates from the depletion method were higher, but the relationship was weak (r²=0.20) on the 13 sites (Stendahl et al., 2013). However, for Ca and Mg there was a much stronger relationship: r²=0.46 for Ca and r²=0.64 for Mg. Thus, based on current knowledge and models, it seems possible to identify a narrow weathering rate interval for Ca and Mg weathering rates, but it seems difficult for K and Na, as discussed in Stendahl et al. (2013).

At Västra Torup and Hissmossa, weathering rates produced by PROFILE and ForSAFE were compared in detail for the first time. The results at Västra Torup showed that the maximum difference between the modelled weathering rates and the median was ±4 %. At Hissmossa, the corresponding difference was ±9 %. Much of the difference could be attributed to the difference between soil moisture input data in PROFILE and modelled soil moisture in ForSAFE. The sandy soil in Hissmossa gave substantially lower modelled soil moisture than the moisture estimates based on field observations used as inputs to PROFILE, resulting in lower weathering rates (Kronnäs et al., 2019).

The intervals of all sites were compared with four reference weathering intervals, based on weathering rate approximations frequently used in the critical load work (Fig. 3; de Vries, 1994; Umweltsbundesamt, 1996). The majority of weathering rate estimates was within or close to the interval outlined for acidic/intermediate parent material with coarse texture.

The weathering rates for the sites with 0.5 m soil depth (including Gårdsjön G1 where the soil depth is 0.47 m) can roughly be divided into four different groups depending on the intervals, except for four sites with contradictory results (Table 4).

Table 4 Classification of the sites with soil depth of approximately 0.5 m, in four classes, based on intervals used in the critical load work of CCE (Coordination Centre for Effects) within the UNECE Convention on Long-range Transboundary Air Pollution (LRTAP convention) (de Vries, 1994; Umweltsbundesamt, 1996): very low, low, and intermediate weathering rates and a group with non-conclusive results, i.e. that did not fit into any of the other groups. Sites were placed in one of the three weathering groups if the median fell within the main interval given and if the maximum and minimum values fell within the extended interval (±5). A and C refer to different profiles, with different soil depths.

^a Corresponding to the lowest span in Fig. 3, “acidic/intermediate parent material, coarse-textured”.
^b Corresponding to the lower part of the two spans “acidic/intermediate parent material, medium-textured” and “basic parent material, all grain sizes” in Fig. 3.
^c Corresponding to the upper part of the span “acidic/intermediate parent material, medium-textured” and the lower-intermediate part of the span “basic parent material, all grain sizes” in Fig. 3.

The assessment of whether the weathering rate intervals are accurate enough depends on their intended use. The weathering rates in relation to forest sustainability assessments are analysed in Sect. 3.3.

Figure 4 Weathering rates on a regional scale on 346 National Forest Inventory (NFI) sites, calculated with the depletion method/total analysis regression approach (a), modelled with PROFILE (b) and with ForSAFE (c), in seven climate regions in Sweden, delimited by black lines.

Figure 5 Box plots for the three methods and for the climate regions: (2) inner part of northern Sweden, (3) coastal part of northern Sweden, (4) western part of central Sweden, (5) eastern part of central Sweden, (6) southwestern Sweden, and (7) southeastern Sweden. The Northwestern mountain region (1) was excluded since it only contained one site.

3.2 Weathering rate comparisons at a regional scale

The weathering rates for the nutrient base cations Ca, Mg and K varied widely within the regions for all methods, but there were no large systematic differences between the medians or ranges for the different methods (Figs. 4–5). However, PROFILE gave generally somewhat lower weathering rates than ForSAFE and the depletion method/total analysis regression approach, with overall medians of 14.3 $\mathrm{mekv}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (PROFILE) and 17.8 $\mathrm{mekv}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (ForSAFE and the depletion method/total analysis regression approach).

For ForSAFE, this difference was most distinct in the northern regions. One explanation regarding the difference between PROFILE and ForSAFE could be differences in the method for estimating mineralogy from total chemistry, where possible minerals have to be set by the modeller. Since qualitative data on mineral contents in soils in most cases are not available at the modelling sites, the mineralogy has to be estimated based on a number of assumptions. In the ForSAFE database, limestone seems to have been set as a possible mineral more often than in the PROFILE database, due to differences in the assumptions made. Since even a small amount of limestone has a large effect on weathering rates, the modelled Ca weathering rates were substantially higher at some sites in the ForSAFE results.

The difference between PROFILE and the depletion method/total analysis regression approach was contradictory to the site-level comparisons, where PROFILE generally gave substantially higher weathering rates than the depletion method (Table 3). This can partly be explained by the site-level comparisons being made for the same depths (50 cm), whereas the regional calculations with the depletion method/total analysis regression approach gave the weathering to the maximum weathering depth, which is often more than the 50 cm (including a 10 cm organic layer) used for PROFILE. Methodological differences between the old and the new calculations for the depletion method/total analysis regression approach can also be part of the explanation, e.g. concerning how the weathering depth has been defined based on curves of the elemental variation with depth, which is partly a subjective operation. Finally, the fact that the regional calculations combine two methods (depletion method and total analysis regression approach), whereas the site-level estimates are based only on the depletion method, may contribute to the differences.

In most cases, the width of the weathering rate intervals showed no major differences between the regions. The difference between the 25th and the 75th percentile was generally 10 to 20 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. 5). An exception was region 5, the eastern part of central Sweden, where the corresponding interval was up to 100 $\mathrm{mekv}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. In region 5, lime-rich soils are common, which can explain this pattern. Other minor differences were that the weathering rates in the southern regions were generally on a somewhat higher level than the others (medians: 14–24 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), whereas the weathering rates in the western part of central Sweden were towards the lower end (medians: 9–14 $\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$).

Figure 6 Weathering rates of base cations calculated with different methods on spruce sites, compared with harvest losses of base cations at whole-tree harvesting (100 % of the stems and 60 % of the branches harvested, and 75 % of the needles on the branches removed). The horizontal dashed lines in the harvesting bars show the levels for stem-only harvesting. The sites are ordered from north to south.

Figure 7 Weathering rates of base cations calculated with different methods on pine sites, compared with harvest losses of base cations at whole-tree harvesting (100 % of the stems and 60 % of the branches harvested, and 75 % of the needles on the branches removed). The horizontal dashed lines in the harvesting bars show the levels for stem-only harvesting. The sites are ordered from north to south.

3.3 Weathering in a sustainability perspective

Although the difference in weathering rates between methods is large on several sites, a number of general conclusions could be drawn from the comparison with harvest losses. For the stem-only harvesting scenario, the harvest losses were generally at the same level or lower than PROFILE weathering rates (Figs. 6–7). The harvest losses at stem-only harvesting were lower or at the same level as weathering rates according to the depletion method in the northern sites but in most cases higher in the southern sites.

Table 5 Weathering rates estimated with different approaches and harvest losses at conventional harvesting (CH, only stems) and whole-tree harvesting (WTH, stems and branches) ($\mathrm{meq}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$).

In five of the seven spruce forests in southern and central Sweden (the sites to the right in Fig. 6), the harvest losses in the whole-tree harvesting scenario were higher than the weathering rates (5 %–740 %), regardless of the method used to calculate weathering rates (Table 5, Fig. 6). The exceptions were Asa, where the weathering rates from the budget approach gave a many times higher weathering rate than the other methods, and Gårdsjön, where the weathering rates from most of the methods were of a similar size to the harvest losses. Despite the variation between methods, the results clearly indicate that whole-tree harvesting is not sustainable in the long term in spruce forests in southern and central Sweden, since the weathering rates generally are substantially lower than the base cation losses at whole-tree harvesting.

On the four spruce sites in northern Sweden (to the left in Fig. 6), PROFILE gave 20 %–130 % higher weathering rates than harvest losses after whole-tree harvesting, whereas the depletion method gave 20 %–50 % lower weathering rates than harvest losses for three of the four sites (Table 5, Fig. 6). The budget approach in Flakaliden gave 220 % higher weathering rates than the harvest losses after whole-tree harvesting. Despite the difference between the methods, the results clearly indicate that the effects of whole-tree harvesting in spruce forests in northern Sweden are substantially smaller than for spruce forests in southern and central Sweden.

All pine sites where comparisons could be made in this study were situated in northern or central Sweden (Fig. 2; Table 1). For all three sites, the PROFILE weathering rates were substantially higher than the harvest losses, both for the stem-only (150 %–240 %) and whole-tree harvesting scenario (90 %–170 %) (Table 5, Fig. 7). Weathering rates calculated with the depletion method were of a similar size to the harvest losses at whole-tree harvesting, except for Kloten where the weathering rate was 50 % lower. Thus, the conclusions about pine forests are similar to those for spruce forests in northern Sweden, i.e. that long-term losses are less of a concern, although the variation in weathering rates makes it difficult to say whether the weathering rates are higher or lower than the harvest losses in those forests.

In the above assessment, the extent to which whole-tree harvesting itself affected the weathering rates was not explicitly considered. As an increased forest harvest intensity leads to slightly more acidic conditions, it could be hypothesised that increased intensity leads to an increased proton-promoted dissolution of minerals, thereby providing a feedback mechanism in which increased weathering could partially alleviate the effect on soil acidity and base cation status. However, according to recent HD-MINTEQ modelling in which PROFILE was used to simulate weathering, the weathering rate was largely unaffected by soil solution pH and by the harvesting method used (McGivney et al., 2019). This was explained as being the net result of the opposing effects of pH and dissolved Al on the weathering rate. While a decreased pH itself leads to an increased weathering rate, it also leads to increased levels of dissolved Al, which is a potent weathering “brake”, offsetting the pH effect. Another source of uncertainties is the potential effect of whole-tree harvesting on mycorrhiza activity, which is further discussed in Sect. 3.4 and in Finlay et al. (2019).

In the assessments of base cation sustainability, it is not only important to focus on uncertainties in the actual soil weathering rates. Other important topics are how much of the weathered material the tree roots can reach; the size of the base cation deposition; the uncertainties in the assessment of base cations through harvesting; and how the base cation losses are distributed between soil, biomass and runoff water. The use of a constant and static rooting depth introduces uncertainties in the sustainability assessments, since root depth varies both spatially and temporally, depending on variations in site conditions (Hodge, 2013; Rosengren and Stjernquist, 2004). The base cation deposition in Sweden is assessed to be of a similar size to the base cation weathering (Akselsson et al., 2007), but a national survey of total base cation deposition, including dry deposition, is not available, so uncertainties of base cation deposition are large. The uncertainties in the assessments of base cation losses at harvesting can be divided in uncertainties in the amount of biomass extracted and the concentration of BC in biomass (Akselsson, 2005). A sensitivity analysis for Ca showed that the lack of site-specific nutrient concentration data was the main source of uncertainties in calculations of harvest losses of Ca, whereas the estimations of biomass available for extraction, and the amount of branches left on the ground, contributed less to the uncertainties (Zetterberg et al., 2014). Finally, the effect of base cation losses in soil is reduced by the fact that the rates of tree growth and leaching decline after whole-tree harvesting, mitigating some of the impacts of harvest on soil base cation status (Zetterberg et al., 2013; Egnell, 2016).

3.4 Potential for biological weathering

Biological weathering often takes place in conjunction with physical and chemical processes, but there is still disagreement over the extent of its quantitative contribution to overall weathering (Finlay and Clemmensen, 2017; Leake and Read, 2017; Smits and Wallander, 2017). Insufficient representation of biological weathering in weathering models such as PROFILE, and its effects of weathering rate uncertainties, has been frequently discussed (Finlay et al., 2009). Below, a description of how biological weathering is presently represented in the PROFILE/ForSAFE models is given, followed by a discussion about potential shortcomings in the light of the latest research about biological weathering. A more thorough description of the state of knowledge and a more comprehensive discussion can be found in the article by Finlay et al. (2019).

Although the four weathering pathways, upon which PROFILE is built (the reaction with H⁺, CO₂ and DOC) are chemical (the dismantling of mineral matrices by charged or dissolving particles to produce free elements), their drivers are strongly dependent on biological activity in the soil (Sverdrup and Warfvinge, 1993). Soil solution H⁺ is determined by the charge balance resulting from uptake, ion exchange, mineralisation of organic matter (solid and dissolved) and hydrological transport, all of which are affected by biological activity. Water availability is directly controlled by water uptake. The partial pressure of dissolved CO₂ stems from root and root symbiont respiration, decomposition and hydrological transport. Finally, DOC is directly and indirectly produced by plants. The transition state theory, governing the weathering kinetics in PROFILE/ForSAFE, dictates that the net weathering rates should decline towards zero near equilibrium. This is represented in the model by retardation factors that increase in strength with the concentrations of the weathering products (Erlandsson et al., 2016). These concentrations are in turn dependent on biological activity, such as uptake reducing nutrient base cation concentrations or the mobilisation of aluminium through biological acidification.

Although the weathering process is strongly affected by biological processes in the current generation of the PROFILE/ForSAFE family of models, the models still fail to capture the biological feedback mechanisms in their entirety. The possible roles of fungi, especially ectomycorrhizal fungi, in biological weathering in boreal forests were summarised by Finlay et al. (2009) and have been the subject of many subsequent studies (Finlay et al., 2019). These fungi can acidify their surrounding environment and release organic acids and siderophores, which may enhance weathering. They can also exert biomechanical forcing and alter interlayer spacing associated with depletion of potassium from biotite (Bonneville et al., 2009). Furthermore, recent work with atomic force microscopy has demonstrated nanoscale alteration of surface topography of minerals and attachment and deposition of organic biolayers by fungal hyphae (McMaster, 2012; Gazzè et al., 2012, 2013; Saccone et al., 2012). Many fungal hyphae produce extracellular polysaccharides (EPS) at their hyphal tips, providing an interface that ensures intimate contact between the hyphae and mineral substrates. The contact area between hyphae and mineral surfaces is increased by EPS haloes (Gazzè et al., 2013), and many fungal exudation products such as organic acids and siderophores may be released into polysaccharide matrices (Flemming et al., 2016) in close proximity to mineral surfaces. Here, they are effectively isolated from the bulk soil solution and may be protected from microbial decomposition by antibiotic compounds also produced by the fungi. This is in contrast to the assumption in the models that soil solution is homogeneous at any given depth, not discerning bulk solution from the said EPS haloes. This may increase the effective concentrations of organic weathering agents at sites of active weathering and structure the bacterial communities associated with particular mycorrhizal fungi (Marupakula et al., 2016).

The potential mechanisms for biological enhancement of mineral weathering and the current debate about the importance of these processes for overall weathering are discussed in detail by Finlay et al. (2019). The biological activity of symbiotic ectomycorrhizal fungi and the evolution of their interactions with their tree hosts have led to systems that are highly adapted to efficient recycling of plant nutrients from organic matter, as well as the release of base cations from mineral substrates through weathering. Ectomycorrhizal mobilisation of N and P through decomposition of organic residues is dependent on carbon supplied from tree hosts. Mycorrhizal weathering of minerals is also dependent on carbon supply from trees, and ongoing experiments (Finlay et al., 2019) suggest that depletion of organic substrates (containing N) will restrict tree growth and therefore also reduce the carbon supply to ectomycorrhizal fungi colonising mineral substrates, with concomitant, negative effects on base cation release from biological weathering. Existing models are therefore probably sufficient to give guidelines about sustainable forestry, including the prediction that, under intensive forestry with removal of organic residues, base cation supply will not be sustainable in the long term. However, the biological feedbacks during transition from one state to another may not be fully covered by the models. For instance, a forest exposed to N deposition may pass from N limitation to limitation by another nutrient, which may have consequences for belowground carbon allocation. Ectomycorrhizal fungi are dependent on carbon supplied from their host plants and can be expected to exert a stronger effect on mineral weathering if they have access to more carbon, which may influence mobilisation of nutrients that are limiting. A possible way to include the biological effects on mineral weathering would be to better describe belowground carbon allocation in the models. Enhanced weathering rates of apatite have been seen when host trees suffer from P shortage, which is known to enhance belowground carbon allocation (Smits et al., 2012). More elasticity in carbon allocation in the models is needed to capture these empirical observations. Furthermore, carbon allocation will also regulate exudation from roots and associated mycorrhizal fungi, which is another process involved in mineral weathering that is not covered in the models.

3.5 Implications of improved model descriptions of base cation exchange and aluminium complexation

Aluminium (Al) and base cation concentrations are the primary weathering brakes in unsaturated soil (Warfvinge and Sverdrup, 1992). Higher concentrations of these elements have a negative effect on the dissolution rates of the minerals containing the elements (Sverdrup et al., 2019). It is therefore imperative to correctly simulate the concentrations of Al and base cations in the soil solution.

Different soil chemical models simulate the dynamics of inorganic Al and base cations in different ways. These can be classified into two categories: (1) simpler ion-exchange equations (e.g. Gaines–Thomas or Gapon) that conceptualise sorption and desorption of Al³⁺, H⁺ and base cations as a series of ion-exchange reactions; and (2) more advanced organic complexation models such as WHAM, NICA-Donnan or SHM (Tipping, 2002; Kinniburgh et al., 1999; Gustafsson, 2001) that treat organic matter as the main cation sorbent, where proton dissociation over a wide pH range drives complexation and exchange of Al and base cations.

In general, the use of organic complexation models to simulate base cation and Al dynamics is strongly supported by empirical evidence (e.g. Tipping, 2002), but for a long time, the simpler ion-exchange equations have been more widely used in popular biogeochemical models such as MAGIC, PROFILE and ForSAFE. However, in 1996, the CHUM model was introduced, which incorporates a version of WHAM (Tipping, 1996), and today SMARTml and HD-MINTEQ provide additional examples of (bio)geochemical codes that employ organic complexation models (Bonten et al., 2011; Löfgren et al., 2017).

Gustafsson et al. (2018) investigated the implications of using the two model approaches on the dynamics of Al, base cations and acidity. Overall, the two model approaches provided the same type of response to changes in input chemistry, implying that, in many cases, there may be a rather limited benefit from using organic complexation models when calculating weathering rates. However, although these results suggest that the current model setup in for example ForSAFE may be sufficient in many cases, certain differences remain between the two categories of models. The Gaines–Thomas and Gapon exchange equations produce a relatively stronger buffering of soil solution pH over a relatively narrow pH range. Together with the general oversimplification of the cation binding process this also causes the ion-exchange equations to overestimate the historical levels of exchangeable base cations (Gustafsson et al., 2018). Consequently, it may be necessary to include organic complexation under such conditions as prolonged or substantial changes in acidic input, such as in the case of sea spray events. Not explicitly simulating organic complexation may require additional coefficients that account for temporal changes in cation selectivity to correctly predict pH, base cations and Al, thereby entailing more uncertainty. Excluding organic complexation can bring into question the ability of biogeochemical models to predict the effect of large changes in acidic input on weathering rates (Gustafsson et al., 2018).

3.6 Prospects for method development

Since quantifying uncertainties for weathering rates is difficult, the use of multiple methods is often proposed as a way of increasing the robustness of weathering rate estimates. However, the number of available methods is low, and all are burdened with different types of limitations and uncertainties. In the review of weathering studies in this paper, only six locations could be found where at least three methods had been implemented and where the criterion of the same depth was fulfilled. Thus, the recommendation in Futter et al. (2012), that at least three independent methods should be used to quantify weathering rates on a site for sustainability assessments, is unrealistic. Nevertheless, comparisons between weathering rates from different approaches for the same sites, and continuous development of the different approaches, will contribute to more robust sustainability assessments. In the next sections, the main uncertainties and development potential related to process descriptions and input data for PROFILE/ForSAFE are discussed, followed by uncertainties and potential development areas for the depletion method/total analysis regression approach and the budget approach.

3.6.1 PROFILE/ForSAFE – process descriptions

By far the most widely used, and most evaluated, method for estimating weathering rates for soils in Sweden is the PROFILE model. The successful testing of weathering rate modelling with ForSAFE (Kronnäs et al., 2019) in QWARTS will be the starting point for more studies on weathering dynamics, using ForSAFE. To better represent ectomycorrhizal fungi processes in the PROFILE/ForSAFE models, and thus reduce uncertainties in modelled weathering rates, three main improvements need to be made: (1) the EPS microenvironments, described in Sect. 3.4, need to be determined in the field and considered in models; (2) methods to distinguish between roots and mycorrhizal hyphae need to be developed, to be able to better represent the process of nutrient uptake and translocation towards the plant root; and (3) more elasticity in carbon allocation in the models is needed to be able to better describe the carbon availability for fungi and to represent the regulation of exudation from roots and associated hyphae.

The uncertainties in the simplified description of base cation exchange and aluminium complexation were generally small, according to studies in QWARTS (Gustafsson et al., 2018). However, we note that a modification would be desirable concerning (1) long-term simulations over hundreds of years when large changes occur in the chemical drivers and (2) sites experiencing frequent or strong sea salt episodes causing large changes in the chemical composition of the influent water.

Based on the assumption that weathering retardation is mainly caused by elevated concentrations of base cations and aluminium, PROFILE/ForSAFE produces reasonable weathering rates in the unsaturated rooting zone (Sverdrup and Warfvinge, 1993; Erlandsson et al., 2016). However, moving into the saturated zone, the strength of the usual weathering brakes fails to slow down the mineral dissolution, which leads to grossly overestimated rates of weathering (Stendahl et al., 2013; Erlandsson Lampa et al. (2019). In this environment, soil solution silicate concentrations play a central role in hindering mineral dissolution (Sverdrup et al., 2019). For this reason, the kinetics of silicate release from mineral dissolution has been added to the traditional elements, as well as the dynamics of silicate concentrations in the soil solution. Erlandsson Lampa et al. (2019) tested a prototype of this addition, and the results proved promising in keeping weathering rates within observation levels in the saturated zone, but this is yet to be implanted and tested in PROFILE/ForSAFE.

Lateral flow has recently been included in ForSAFE, and a new version, ForSAFE-2D, has been developed (Zanchi et al., 2016). The model has been evaluated on the basis of hydrological flows and chloride concentrations and transport, with good results. Evaluating the modelled base cation concentrations in surface water highlighted the need for adjusting the weathering brakes (see discussion above about silicate brakes) and also a need to revisit the decomposition process descriptions, thereby validating them for the saturated zone. Further development of ForSAFE-2D has the potential to provide a mechanistic tool for assessing weathering rates also for surface water applications. The importance of correctly defining the flow pathways and residence times for the delivery of weathering products to the surface waters, and the potential value of concentration–discharge relationships for calibrating biogeochemical models, was explored by Ameli et al. (2017).

3.6.2 PROFILE/ForSAFE – input data

Although continuously improved process descriptions are desirable to get more robust weathering estimates, improvements related to input data are more urgent. Mineralogy, specific surface area and soil moisture are of key importance in weathering modelling but are often burdened with high uncertainties. To reduce input data uncertainties, a focus should be placed on those three parameters.

Mineralogy inputs to PROFILE/ForSAFE are often estimated from total chemistry with the A2M model (“Analysis to Mineralogy”, Posch and Kurz, 2007), since direct mineralogy measurements are not available on most sites. To accurately estimate a probable mineralogy, not only are good soil chemistry measurements required, but also information about which minerals can be expected in the soil. In Sweden, four different geographical mineralogy regions have been used since the 1990s to assign qualitative mineralogy to a site (Warfvinge and Sverdrup, 1995). Casetou-Gustafson et al. (2019a) compared weathering rates calculated based on three sets of mineralogies: one based on direct measurements of quantitative mineralogy, one based on normative modelling with A2M using direct measurements of qualitative mineralogy, and one based on normative modelling with A2M using data from the regions mentioned above. It could not be concluded that the A2M runs based on direct measurements of qualitative mineralogy gave better results. Although these results strengthen the credibility for the normative mineralogy regions, Casetou-Gustafson et al. (2019a) recommend continued work to reduce uncertainties related to mineralogy, mainly by revisiting and, if appropriate, updating mineral rate coefficients. More comparisons of weathering rates from normative mineralogies based on generalised and site-specific quantitative mineralogy are needed to adequately assess whether the regional divisions need to be revised and refined in order to further reduce the uncertainties in the mineralogy estimates.

A2M gives as output a multidimensional space of solutions, all of which have the same probability. Often, the centre point of the space is used for weathering calculations. However, the span can be quite broad, which leads to uncertainties in the calculated weathering rates (Casetou-Gustafson et al., 2019a). Future research focusing on constraints that could help to narrow the space of possible solutions that A2M creates, e.g. based on the grain size distribution, could reduce those uncertainties.

Minerals are assumed to be evenly distributed among grain sizes in PROFILE and ForSAFE. The effect of this assumption has not been fully analysed. The most obvious example showing that minerals are not evenly distributed among grain sizes is clay minerals, which are found in the clay fraction. The extremely high surface area of clays leads to very high base cation weathering rates when the clay fraction is high, although the content of base cations is low. Due to this, Phelan et al. (2014) introduced a correction factor. A thorough analysis of all grain size fractions can help to further refine these methods.

The surface area of soils is often calculated with regressions based on old BET measurements (Warfvinge and Sverdrup, 1995). The regressions reveal that the uncertainties are large. Revisions of the regressions, based on a larger data material, could reduce the uncertainties.

The soil moisture is one of the most important factors that introduces large uncertainties in the results, both in PROFILE where it is an input (Rapp and Bishop, 2003) and in ForSAFE where it is modelled based on hydrological parameters (Kronnäs et al., 2019). Improved input data quality for soil moisture would substantially reduce uncertainties in PROFILE, and, even more importantly, soil moisture modelled by ForSAFE needs to be evaluated, and the sensitivity to soil input data needs to be examined.

3.6.3 The depletion method and the total analysis regression approach

Next to the PROFILE model, the depletion method is the most used method in Sweden, often in combination with the total analysis regression approach. To further evaluate the accuracy of results from the depletion method, as a proxy for the weathering rates of today, the reliability of the assumptions needs to be further evaluated, and the relationship between the average weathering rate since the last glaciation and today's weathering rate needs to be assessed. The latter can be done by performing ForSAFE modelling on a site where the depletion method has been applied. A similar exercise has been done with the SAFE model (Warfvinge et al., 1995), but the inclusion of tree growth and decomposition in ForSAFE can be expected to improve the results. Furthermore, standardised methods for setting the weathering depth based on the elemental content curve and for the analysis of fulfilment of the requirements in the soil profile would enable objective and comparable estimates, including requirements that must be fulfilled for soil profiles to be regarded as undisturbed. The total analysis regression approach will give more robust results if more depletion method estimates are available for the regressions.

3.6.4 The budget approach

Different applications of the budget approach handle the distinction between sources of base cations in the soil in different ways, affecting the uncertainties in the estimated weathering rates. The uncertainties in base cation deposition add to the overall uncertainties. For a fair comparison between weathering rates from the budget approach and from other methods, ways to distinguish between different sources and sinks need to be further developed.

An advantage of the budget approach using the Sr isotope ratio is that it can distinguish between weathering and release from the exchangeable pool. As for all budget approaches, deposition and leaching measurements are required as inputs. The few comparisons made in this study show promising results, and we therefore encourage estimates on more sites to enable evaluation of the budget approach based on the Sr isotope ratio.

In the MAGIC model, the release from the exchangeable pool is thoroughly modelled, but for other sources and sinks of base cations, the same problems apply as for other budget approaches, e.g. uncertainties in base cation deposition. These uncertainties exacerbate the uncertainties in weathering rates that derive from the mass balances in MAGIC. Nevertheless, MAGIC theoretically provides a good basis for conducting independent weathering rate assessments. On sites with relatively small input data uncertainties, our recommendation is to carry out such comparisons.

4 Conclusions

Uncertainties in weathering rates have often been presented as an obstacle in the assessment of sustainable forestry. The comparison between approaches in this paper, on a regional level as well as on a site level, suggests that both weathering rate gradients and approximate weathering rate levels can be captured with available methods. Although the variation in weathering estimates was large on single sites, most of the sites could be grouped into broader classes representing very low, low and intermediate weathering rates, which can be used for general, but not specific, weathering rate assessments at the site level. The more and better input data that are available, and the more methods that are applied and compared for a single site, the more robust the overall assessments that can be done at the site level are, provided that the conceptual differences, boundary conditions and assumptions between methods are kept in mind.

Based on the results from this study, we argue that modelled weathering rates can be used for sustainability assessments, as long as the uncertainties, i.e. the intervals on single sites presented in this paper, are recognised. The ability to draw conclusions about sustainable forestry at the site level depends not only on uncertainties in weathering rates, but also on other site properties, related to forest properties and other base cation flows, such as base cation deposition, and the associated uncertainties. Irrespective of the uncertainties related to the sustainability assessments, a robust conclusion was that weathering rates in spruce forests in southern and central Sweden generally were substantially lower than the harvest losses at whole-tree harvesting, indicating that whole-tree harvesting without nutrient compensation is not sustainable in these areas. There is less risk of a negative effect for spruce forests in northern Sweden, as well as pine forests in central and northern Sweden.

The research performed in the five years of the QWARTS programme supports the continued use of the PROFILE/ForSAFE models. ForSAFE is the only method that gives time-resolved results, i.e. the only method that can be used to study dynamic effects of changing climate and changing management methods. Although there is still scope for improving process understanding and incorporation of that understanding into PROFILE and ForSAFE, e.g. regarding biological weathering and weathering brakes, the most important way to reduce uncertainties in modelled weathering rates is to reduce input data uncertainties, mainly regarding soil texture and associated hydrological parameters. However, it is also important to continue to compare with results from the depletion method and the budget approach.