Critical Loads of Sulphur and Nitrogen for Terrestrial Ecosystems in Europe and Northern Asia Using Different Soil Chemical Criteria

Abstract

A critical load data base was developed for Europe and Northern Asia using the latest data bases on soils, vegetation, climate and forest growth. Critical loads for acidity and nutrient nitrogen for terrestrial ecosystems were computed with the Simple Mass Balance model. The resulting critical loads are in accordance with critical loads from previous global empirical studies, but have a much higher spatial resolution. Critical loads of acidity are sensitive to both the chemical criterion and the critical limit chosen. Therefore a sensitivity analysis of critical loads was performed by employing different chemical criteria. A critical limit based on an acid neutralizing capacity (ANC) of zero resulted in critical loads that protect ecosystems against toxic concentrations of aluminium and unfavourable Al/Bc ratios, suggesting that ANC could be an alternative to the commonly used Al/Bc ratio. Critical loads of nutrient nitrogen are sensitive to the specified critical nitrate concentration, especially in areas with a high precipitation surplus. If limits of 3–6 mg N l⁻¹ are used for Western Europe instead of the widely used 0.2 mg N l⁻¹, critical loads double on average. In low precipitation areas, the increase is less than 50%. The strong dependence on precipitation surplus is a consequence of the simple modelling approach. Future models should explore other nitrogen parameters (such as nitrogen availability) instead of leaching as the factor influencing vegetation changes in terrestrial ecosystems.

1 Introduction

During the past two decades, critical loads of sulphur (S) and nitrogen (N) have been used as an indicator of ecosystem sensitivity to acidification and eutrophication under the Convention of Long-range Transboundary Air Pollution (LRTAP) within the United Nations Economic Commission for Europe (UNECE). Starting with the Sulphur Protocol of 1994, critical loads, as part of integrated assessment modelling, were used in the negotiations on emission reductions (Bull et al. 2001; Gregor et al. 2001). More recently critical loads have also been used in the revision of the European Commission’s National Emission Ceilings (NEC) directive (Amann et al. 2006). Critical load maps used for policy support in Europe consist of data submitted by National Focal Centres as well as of critical loads based on models applied to a general, European-wide background database (Hettelingh et al. 1995a; Hettelingh et al. 2001, 2007).

Several studies have determined critical loads for nitrogen and acidity for terrestrial ecosystems at the European (De Vries et al. 1994; Kuylenstierna et al. 1998), SE-Asian (Hettelingh et al. 1995b), Northern Asian (Bashkin et al. 1995; Semenov et al. 2001) and global scale (Kuylenstierna et al. 2001; Bouwman et al. 2002). The global scale studies assigned low resolution critical loads based on literature data on ecosystem sensitivity, whereas the study by De Vries et al. (1994) used the simple mass balance model. In recent years, high-resolution data bases have become available for soils (JRC 2006), land cover (Bartholome et al. 2002), climate (New et al. 1999, 2000) and forest growth (Schelhaas et al. 1999). These data bases provide much more spatial detail for Europe and Northern Asia than those used in previous studies. A harmonized land cover map has been recently completed by combining Corine land cover (version 12/2000 extended coverage) with additional sources from the Stockholm Environment Institute (SEI; Cinderby et al. 2007) to provide a Pan-European coverage on a 100 × 100 m grid. Furthermore, EUNIS ecosystem codes (Davies et al. 2004) have been assigned to all land cover classes (Slootweg et al. 2005).

In addition to the need to use more recent data, it has been suggested that chemical criteria other than the widely used Al/Bc ratio could be used for the computation of critical loads of acidity for soils (Holmberg et al. 2001; UBA 2004). Several national studies have shown that the choice of the chemical criterion can have a strong influence on critical load values (e.g. Hall et al. 2001; Aherne et al. 2001, for the UK and Ireland), and on critical load exceedances (e.g. Holmberg et al. 2001 for Finland). Similarly a recent study has shown that the widely used nutrient nitrogen limit of 0.2 mg N l⁻¹ related to vegetation changes in forests is applicable for Scandinavia but not for Western Europe (De Vries et al. 2007). These authors suggested new limits for several forest and vegetation types. Finally, in recent years the (mostly) Asian part of the UNECE region—termed EECCA (Eastern Europe, Caucasian and Central Asian) countries—as well as Cyprus and Turkey have become more involved in the work under the LRTAP Convention, necessitating the extension of the European background data base to these countries.

In this study we describe the latest data bases and methodologies and use them to compute critical loads of S and N as well as of nutrient N for terrestrial ecosystems in Europe and Northern Asia, comprising the successor states of the former Soviet Union. In addition, uncertainties in the regional patterns of critical loads were explored using a range of chemical criteria. Furthermore, we show how various chemical criteria are interconnected introducing the concept of equivalent criteria.

2 Methods

2.1 The Critical Load Model

Critical loads were computed with the so-called Simple Mass Balance (SMB) model, which links deposition to a chemical variable in the soil or soil solution, which can be associated with ecosystem effects; and the violation of specific values of such a variable (the ‘critical limit’) can be linked to ecosystem damage. In this way deposition(s) are linked to a ‘harmful effects’: a low critical load implies a high sensitivity of the ecosystem to deposition, and vice versa. The critical load of S and N acidity is not a single value, but a trapezoidal function in the N–S-deposition plane characterised by three quantities, CL_max(S), CL_min(N) and CL_max(N). For the eutrophying effect of N, the critical load is given as a single number, CL_nut(N). These quantities are derived in Appendix A and are illustrated in Fig. 1. In this paper we specifically look at the distribution of CL_max(S)—also called the critical load of acidity—and CL_nut(N), the main quantities characterising ecosystem sensitivity to S and N deposition.

Fig. 1

Critical load function of sulphur and nitrogen, defined by the three quantities CL_max(S), CL_min(N) and CL_max(N) of the acidity critical loads and the critical load of nutrient N, CL_nut(N). The grey area below the CL function denotes deposition pairs (N _dep , S _dep ) resulting in an ANC leaching smaller than ANC _le,crit and a nitrate leaching less than Q∙[N] _acc , i.e. non-exceedance of critical loads

2.2 Chemical Criteria

Critical loads link deposition to ecosystem effects via soil chemical criteria (critical limits). These limits are based on dose–response relationships between chemical characteristics and ecosystem functioning. A critical load equals the deposition that results in a steady state in an ecosystem compartment (e.g. soil, groundwater, plant) that does not exceed the selected critical limit, thus preventing ‘significant harmful effects on specified sensitive elements of the environment’ (Nilsson and Grennfelt 1988). Consequently, the selection of the chemical criterion and its critical limit is a crucial step in deriving a critical load, and has to be guided by the negative effect(s) one wants to avoid.

To date mostly soil chemical criteria (e.g. nitrate and aluminium concentrations or aluminium to base cation ratios) have been used to derive critical loads with simple steady-state models. The relationship between these critical limits and the ‘harmful effects’ is one of the largest sources of uncertainty. For surface waters, there is a clear relationship between damage (fish dieback) and critical loads exceedance (Jensen and Snekvik 1972; Dickson 1978; Henriksen et al. 1989), but for terrestrial ecosystems the correlation is less convincing, although damage to tree crowns has been recorded in association with exceedance of critical loads (e.g. Nelleman and Frogner 1994). One reason may be that effects are mainly invisible, such as effects of aluminium (Al) on fine root growth. Several authors have doubted the validity of the widely-used critical Al/Bc ratio as an indicator for harmful effects on forests, as no field evidence of such a relationship on mature trees has been demonstrated (Løkke et al. 1996; De Wit et al. 2001; Göransson and Eldhuset 2001; Nyberg et al. 2001; Nygaard and de Wit 2004). Some authors have therefore suggested to use other criteria such as a critical value for base saturation (Bsat) or acid neutralising capacity (ANC) leaching (Holmberg et al. 2001) or to preserve stable pools of aluminium hydroxides (De Vries 1993). There is more empirical evidence for relationships between atmospheric N deposition and effects on plant species diversity, and these form the basis for (empirical) critical loads of nutrient nitrogen (e.g. Bobbink et al. 1988; Nordin et al. 1998; Clark and Tilman 2008). However, use of a single critical N concentration for Europe to derive critical N loads with the SMB model is not a valid option: De Vries et al. (2007) demonstrated that it is possible to define several associated ‘critical’ N concentration ranges for various ecosystems in Europe to protect for vegetation change, in contrast to the single values previously used. They also showed that effects are better related to other predictors such as the N mineralization flux (De Vries et al. 2007). In this study several criteria were used to investigate the effect of the choice of the criterion on the (patterns in) critical loads (Table 1).

Table 1 Chemical criteria for computing critical loads

For a given site a fixed relationship exists between each of the acidification criteria, referred to hereafter as ‘equivalent criteria’ (see Appendix and Fig. 2). The widely used criterion of Al/Bc = 1 leads to (strongly) negative ANC values in the soil solution, except for soils with a very low base cation (Bc) concentration (Fig. 2). The decrease in ANC equivalent to Al/Bc = 1 for increasing Bc concentrations is due to the fact that with increasing Bc concentration the Al concentration has to increase to keep the Al/Bc ratio constant, and the increasing Al concentration leads to a decreasing ANC. A critical ANC of zero results in an equivalent base saturation between about 5% for soils poor in base cations (such as podzols) to about 30% for richer soils. Aiming at a base saturation for all soils of, e.g., 15% requires an (unrealistically) high ANC especially in soils poor in base cations.

Fig. 2

Functional relationship between ANC and molar Al/Bc ratio (left) and base saturation (right). The graphs on the left are for [Bc] = 1 (top curve), 5, 20 and 50 meq m⁻³; the ones to the right for K _Gap √[Bc] = 1 (top curve), 2, 4 and 8 meq m⁻³. All graphs are made with $K_{HCO_3 } p_{\operatorname{CO} _2 } = 200$ (mmol m⁻³)², m∙ DOC = 20 mmol m⁻³ and pK ₁  = 4.5 (for notation see Appendix). The parameters used to compute these curves cover the majority of soils in Europe and Northern Asia

2.3 Geographical Data Bases

The required input data for critical load calculations consist of spatial information describing climatic variables, base cation deposition and weathering, nutrient uptake and N transformations and were derived combining maps of soils, land cover and forest growth regions. To cover the entire geographical area of interest, several thematic maps had to be combined:

• (a) Land cover: The harmonised land cover map produced by the CCE and SEI was used for Europe (Slootweg et al. 2005, 2007). For the EECCA countries we used the Global Land Cover 2000 project map at 1 km resolution (Bartholome et al. 2002). Only forests (Eunis code ‘G’) and (semi-)natural vegetation (codes ‘D’, ‘E’ and ‘F’) were considered in this study.

• (b) Soils: The European Soil Database v2 polygon map (JRC 2006) at a scale 1:1 M was used for Europe including the entire Russian territory, Belarus and Ukraine. For the other CIS (Commonwealth of Independent States) countries, Turkey and Cyprus the less detailed FAO 1:5 M soil map (FAO 1988) was selected.

• (c) Forest growth: Average forest growth was derived from an updated data base of the European Forest Institute (EFI), which contains growth data for a variety of species and age classes in about 250 regions in Europe (Schelhaas et al. 1999). This map was combined with a map of 74 administrative regions in Russia for which forest stock data are provided by Alexeyev et al. (2004). For other CIS States, Cyprus and Turkey a map was used that delineates the forest growth regions in these countries.

Overlaying these maps and merging polygons with common soil, vegetation and region characteristics within blocks of 10 × 10 km (a subdivision of the EMEP 50 × 50 km grid) resulted in about 3.8 million computational units with a total area of 16.6 M km². In the standard model runs, we used only computational units larger than 1 km² reducing their number to 1.3 M but occupying 96% of the total area. One simulation was performed with all units included to determine the effect of leaving out the very small units on the distribution of critical loads.

The soil maps are composed of so-called soil associations, each polygon on the map representing one soil association. Every association, in turn, consists of several soil typological units (soil types) that each occupy a known percentage of the soil association, but with unknown location within the association. The soil units on the maps are classified into more than 200 soil types (European Soil Bureau Network 2004), with associated attribute data such as soil texture, parent material class and drainage class. Six texture classes are defined, based on clay and sand content (FAO-UNESCO 2003). The drainage classes, which are used to estimate the denitrification fraction, are derived from the dominant annual soil water regime (FAO-UNESCO 2003; European Soil Bureau Network 2004).

Combining soil- and land cover maps shows that podzols are most frequent (about 18%), especially in the north-western part of the area, followed by gleysols (17%), cambisols (about 13%) and regosols and podzoluvisols (Table 2). Whereas in Europe podzols are by far the most important soils (De Vries et al. 1994; Posch and Reinds 2005), large natural areas in the EECCA countries are occupied by wet soils such as gleysols. Soils in natural areas mainly occur on coarse (texture class 1, 41%) and medium soil textures (class 2, 48%). Natural vegetation on fine textures (classes 3–5) is rare and occupies only about 10% of the area. About 4% of the vegetation is located on peat soils.

Table 2 Area and share of the 10 most common vegetation-soil combinations

There are inaccuracies in these estimates, because the soil map consists of soil associations. The map overlay thus gives an area for each association, not for each soil type. Vegetation has been assigned evenly to all soil types within the association, which in reality will not always be the case. Nevertheless, a previous study showed that an uneven allocation of forests, in which they are assigned to poor and steep soils first, yields an almost identical distribution of forest-soil combinations in Europe as an even distribution (De Vries et al. 1994).

2.4 Meteorology and Hydrology

The annual water flux through the soil at the bottom of the rooting zone is required to compute the concentration and leaching of compounds. The bottom of the root zone was set at 50 cm, except for lithosols which have a soil depth of 10 cm only. The leaching rate was estimated from meteorological data and soil properties. Long-term (1961–1990) average monthly temperature, precipitation and cloudiness were derived from a high resolution European data base (New et al. 1999) that contains monthly values for the years 1901–2001 for land-based grid-cells of 10′ × 10′ (approximately 15 × 18 km in central Europe). For sites east of 32°, a coarser 0.5° × 0.5° global database from the same authors was used (New et al. 2000).

Evapotranspiration was calculated with a sub-model used in the IMAGE global change model (Leemans and Van den Born 1994) following the approach by Prentice et al. (1993). Potential evapotranspiration was computed from temperature, sunshine and latitude. The effect of snow cover on evapotranspiration was included by simulating accumulation and melting of a snow layer at each site using the temperature and precipitation data. Actual evapotranspiration was then computed using a reduction function for potential evapotranspiration based on the available water content in the soil described by Federer (1982). Soil water content is in turn estimated using a simple bucket-like model that uses water holding capacity and precipitation data. The available water content (AWC) and the water content at wilting point were estimated as a function of soil type and texture class according to Batjes (1996). Batjes (1996) provides texture class dependent AWC values for FAO soil types based on an extensive literature review and developed a transfer function to compute wilting point from soil texture and soil organic carbon content. A complete description of the hydrological model (without the snow module) can be found in Reinds et al. (2001).

The critical load models implicitly assume free draining soils. Large areas of Northern Russia, however, have shallow permafrost where the critical load model in its present form cannot be applied. Therefore, soils in areas with an average monthly temperature below zero for at least eight months of the year have been excluded from the simulations. This area corresponds well with the areas with shallow permafrost reported in FAO-UNESCO (2003).

2.5 Base Cation Deposition and Weathering

Base cation deposition for Europe was taken from simulations with an atmospheric dispersion model for base cations (Van Loon et al. 2005). For northern Asia calcium (Ca) deposition was taken from a global map computed with a model of Tegen and Fung (1995) using estimates of soil Ca content (Bouwman et al. 2002). Comparing the spatial patterns of Ca deposition of the European map and the global map, it is clear that the global Ca map underestimates the deposition in Europe by at least a factor of two. This underestimation was also recognized by Lee et al. (1999) who attribute the difference mainly to the fact that they did not include all important, local, sources of Ca in their modelling. Using both maps would thus lead to a non-smooth transition in BC deposition in eastern Europe. Since the European map does include all sources (natural and anthropogenic) it was taken as the reference and the Ca from the global were multiplied by two to generate a consistent deposition pattern over the entire modelling area in the combined map.

Magnesium (Mg) and potassium (K) deposition are also needed in the EECCA countries; relationships between Ca deposition and Mg and K deposition were derived on the basis of measurements at 95 EMEP/CCC monitoring stations in Europe (Hjellbrekke et al. 1997). Because of the different origin of base cations, the spatial patterns in these ratios are far from constant. In southern areas Bc input is dominated by Ca from Saharan dust, whereas in Northern Europe Mg and K become more important with Mg dominating the Bc input in coastal regions. Mg deposition was modelled as a function of Ca deposition, Ca _dep (eq ha⁻¹ per year), and the distance to coast:

$$Mg_{{dep}} = a\, \times Ca_{{dep}} \, + \,b{\left( {fDistToCoast} \right)}$$

Regression (r ² = 0.375) yields a = 0.4748 and b = 239.6, 171.6, 104.9, –18.2 and –39.4 for distance-to-coast classes <10, 10–20, 20–50, 50–100 and >100 km, respectively. Mg deposition thus increases with increasing Ca deposition and strongly decreases with distance from the coast. K deposition was estimated as a function of Ca deposition, Ca _dep (in eq. ha⁻¹ year⁻¹), and latitude, Ylat (in degrees):

$$K_{{dep}} = a\, + \,b \cdot Ca_{{dep}} \, + c \times Ylat$$

Regression yields a = −95.1, b = 0.2419 and c = 1.731 with r ² = 0.552; K deposition increases with increasing Ca deposition and with latitude.

Weathering of base cations was computed as a function of parent material class and texture class and corrected for temperature, as described in (UBA 2004). For Europe and all of Russia parent material was obtained from the 1:1 M soil map and for the rest of the territory as a function of soil type as described by De Vries (1991). The texture class attribute is common to both soil maps. From the total BC weathering, the weathering rates of Ca, Mg, K and Na were estimated as a function of clay and silt content for texture classes 2 to 5 (following Van der Salm 1999) and as fixed fractions of total weathering for texture class 1 (De Vries 1994).

2.6 Nutrient Uptake, Nitrogen Immobilization and Denitrification

The net growth uptake of Bc and N by forests was computed by multiplying the estimated annual average growth of stems and branches with the element contents of base cations and N in these compartments based on an extensive literature review by Jacobsen et al. (2002). The average nutrient contents of spruce and pine were assigned to conifers and the average of oak and beech to deciduous forests. An average of these two values was used for mixed forests. Wood densities of 450 and 650 kg.m^-3 as well as branch-to-stem ratios of 0.15 and 0.20 for coniferous and deciduous trees, respectively, have been used (Kimmins et al. 1985). For mixed forests the averages of these values were applied. For other (semi-)natural vegetations, net uptake was set to zero assuming that no net removal of N and Bc occurs.

Forest growth for Europe was derived from the EFI database (Schelhaas et al. 1999) that provides measured growth data for about 250 regions in Europe for various species and age classes. Growth was assessed by computing the area-weighted average growth over all age classes for each combination of region and tree species group. Forest growth for Russia was estimated from data by Alexeyev et al. (2004) who compiled statistical data on growing stock and areas of stocked land from available data sources, tabulated for 74 administrative regions within Russia. Alexeyev et al. (2004) provide areas per region of conifers forest, deciduous hardwood and deciduous softwood forests for the age classes young, middle-aged, maturing and mature/over mature forests as well as the standing biomass per region for these species and age classes. Net growth was estimated by computing the standing volume per hectare (using total volumes and stocked areas) per age class, assuming ages of 30, 60, 90 and 140 years and fitting a logistic growth curve to the volume-age data. Finally the average growth was obtained from this growth curve. For the other CIS states, growth rates were obtained from Prins and Korotkov (1994), who provide the growing stock per hectare. Assuming an average stand age of 60 years gives an approximation of average forest growth in these regions.

For Turkey, growth rates were kindly supplied by the Turkish ICP Forest National Focal Centre as growth rates for thirty species and two forest states (degraded and non-degraded). Furthermore, for a few tree species, growth rates were supplied for coppice and high forest separately. These data were combined with a map showing the distribution of species over Turkey to arrive at growth rates per region per species group (conifers, broadleaves). For Cyprus a crude approximation of an average growth rate of 0.8 m³.ha^-1 was made based on the average standing biomass of 43 m³.ha^-1 given by FAO (2000) and assuming an average stand age of 60 years.

The denitrification fraction, f _de , was computed as a function of the soils’ drainage status (Reinds et al. 2001) and varies between 0.1 for well drained soils to 0.8 for peaty soils. The long-term net N immobilization was set at 1 kg N ha⁻¹a⁻¹, which is at the upper end of the estimated annual accumulation rates for the build-up of stable C–N compounds in soils (UBA 2004).

2.7 Al-H Relationship and Organic Acids

The Al concentration is computed from a gibbsite equilibrium and the equilibrium constant K _gibb is estimated as a function of soil texture class based on simultaneous measurements of Al concentration and pH at about 150 European forest monitoring plots (De Vries et al. 2003). The dependence of K _gibb on the (soil) temperature T (°C) is modelled according to the Van’t Hoff equation:

$$K_{{gibb}} {\left( T \right)}\, = \,K_{{gibb}} {\left( {T_{0} } \right)}\, \cdot \exp {\left[ {\frac{{\Delta H}}{R} \cdot \,{\left( {\frac{1}{{273 + T_{0} }} - \frac{1}{{273 + T}}} \right)}} \right]}$$

where ΔH is the reaction enthalpy (=−95,490 J mol⁻¹), R the gas constant (=8.314 J mol⁻¹ K⁻¹) and T ₀ (=10°C) is a reference temperature. The same database was used to derive a relationship between DOC concentration and soil pH and texture. In turn, estimates of soil pH were obtained from an extensive soil data base with about 6,000 soil profiles in Europe (Van Mechelen et al. 1997) and a data base for the Russian territory (Stolbovoi and Savin 2002). The same datasets were used to estimate soil organic carbon contents needed for estimating soil water holding capacity (see section 2.4).

3 Results

The input data have been derived for each of the 1.3 M receptors. In the following sections results are presented as maps showing the median values within each 50 × 50 km EMEP grid cell for input data. For critical loads, however, fifth percentiles are shown to indicate the most sensitive ecosystems. The area in Northern Russia with no data shown is the area excluded from modelling because of shallow permafrost.

3.1 Input Data

Leaching fluxes vary from less than 100 mm per year in arid regions such as central Spain, central Turkey and large parts of the southern CIS states to >300 mm per year in areas with high precipitation such as along the west coast of Europe and along many mountain ranges (Fig. 3a).The uncertainty in the leaching flux is linked to the reliability of the climate data: values in western Europe are more certain than those in the EECCA area, as the density of rainfall stations used to estimate the grid rainfall is much higher in Europe than in the EECCA area (New et al. 1999). Median base cation deposition shows a strong north-south gradient with values >600 eq ha⁻¹ per year in southern Europe and the southern parts of the CIS states, caused by high dust input from nearby desert areas, and very low values of <200 eq ha⁻¹ per year in the northern part of the modelled area (Fig. 3b). The map also shows that a reasonably consistent spatial pattern is achieved even though two data sources were used.

Fig. 3

a Grid-median leaching flux from root zone (mm per year); b base cation deposition (eq ha⁻¹ per year)

Very low weathering rates (<150 eq ha⁻¹ per year) are found in most of Scandinavia where poor soils prevail and temperatures are low (Fig. 4a). The same holds for large parts of Northern Russia. Very low weathering rates also occur in central and western Spain in areas dominated by acid dystric regosols developed on granites. High weathering rates (>1,000 eq ha⁻¹ per year) are confined to regions with soils developed on volcanic materials and especially in areas dominated by calcareous soils that occupy parts of Spain, France, Hungary, large parts of Turkey and most of the areas with forests and/or natural vegetation of e.g. Kazakhstan, Uzbekistan and Turkmenistan. The high weathering rates in eastern Russia occur in areas dominated by rendzic leptosols, soils with an organic rich topsoil overlying parent material with at least 40% calcium carbonate equivalent (FAO 1988). This is also reflected in the forest growth rates in this region in eastern Russia which are somewhat higher than in the surrounding areas (see Fig. 4b), showing the influence of site quality on growth.

Fig. 4

a Grid-median base cation weathering (eq ha⁻¹ per year); b grid-median forest growth rate (m³ ha⁻¹ per year)

Net uptake of nitrogen and base cations in forest ecosystems is determined by nutrient contents and growth rate. Median growth rates of forests show the well-known pattern in Europe were high growth rates are found in central Europe were climate, site quality and intensive management allow highly productive forests (Fig. 4b). Low growth rates in Europe are confined to arid regions such as central Spain, parts of France and Turkey and to areas with low temperature and poor site quality (shallow, poor soils) such as northern Scandinavia. In the EECCA territory growth rates are generally low (1–3 m³ ha⁻¹ per year), although relative high forest productivity can be found in the areas west of the Ural mountains and in Georgia. Very low growth rates (<0.5 m³ ha⁻¹ per year) occur in arid regions such as Tajikistan, Turkmenistan and Uzbekistan. Growth rates for Russia are somewhat uncertain as they were derived indirectly from growing stock data. For the other CIS states and Cyprus, where only one growth rate per country could be assigned, the growth data do not represent the spatial variability in growth rates. Although the European data are based on a large data set, it is clear that some border effects occur, probably due to the fact that for some countries (e.g. Ukraine and Romania) the area-representation of the supplied data is relatively poor.

3.2 Critical Loads of Acidity (Sulphur)

The critical load of acidity, CL_max(S), based on Al/Bc = 1, is computed as the sum of base cation input through weathering and deposition minus the removal of base cations by uptake minus ANC leaching (see Appendix). Highest critical loads are thus found in areas with high base cation weathering and/or base cation deposition such as along the Mediterranean coast, parts of eastern Europe and the southern parts of the CIS states (Fig. 5a). For calcareous soils CL_max(S) has been set to 10,000 eq ha⁻¹ per year, representing the very high weathering rates in such soils. Furthermore, high critical loads are found in areas with high acidity leaching (areas with a high precipitation surplus) such as along the coast of north-western Europe. Such high acidity leaching due to high base cation (especially Mg) input from sea salts is considered inappropriate by some authors, who have used alternative criteria, e.g. pH in Ireland (Aherne et al. 2001) and Ca/Al ratio in the UK (Hall et al. 2001). Lowest critical loads are found in areas with low weathering rates associated with coarse soils on acid parent material such as central Spain and/or low temperatures (Scandinavia and northern Russia).

Fig. 5

a Fifth percentile critical load CL_max(S) (eq.ha⁻¹ per year); b cumulative frequency distribution of CL_max(S) for three vegetation classes (eq.ha⁻¹ per year)

The cumulative frequency distributions of CL_max(S) for the three vegetation groups: forests, grasslands and heath-lands/tundra show that the critical load distribution for forests and heath land/tundra are very similar (Fig. 5b). For the heath-lands/tundra a larger fraction is located on calcareous soils, illustrated by the larger fraction of CL_max(S) of 10,000 eq ha⁻¹ per year. Grasslands generally have higher critical loads because they occur on richer soils than forests and heath-lands/tundra: average computed weathering rates for grassland soils are about twice as high as those for forest and heath land soils. About 35% of the grasslands are located on calcareous soils.

3.3 Critical Loads of N

The minimum critical load of N (see Fig. 1) consists of the long-term immobilization and net uptake. Because we assumed no removal of growing material from natural vegetations such as grasslands and heath lands, the minimum critical N load consists of the fixed N immobilisation of 1 kg N per year (=71 eq ha⁻¹ per year) only. For forests, net N uptake is accounted for and CL_min(N) ranges between 150 eq ha⁻¹ per year in low productive areas to about 600 eq ha⁻¹ per year in regions with high forest growth.

The critical load of nutrient N, CL_nut(N), is computed from CL_min(N) by adding a critical N leaching and denitrification. Average critical nitrogen concentrations were set to 0.2–0.4 mg N l⁻¹ for forests, 3 mg N l⁻¹ for natural grassland (UBA 2004) and 4 mg N l⁻¹ for heath land (De Vries et al. 2007). Figure 6a show the spatial patterns in CL_nut(N). Highest critical loads are confined to regions dominated by grass and heath lands with a high precipitation surplus, leading to high N leaching rates. Such regions are Ireland, the western parts of the UK and Norway, northern Spain and the region along the Adriatic coast. Low critical loads are found in arid regions were N leaching is low such as in the southern part of the CIS states, most of Turkey and parts of Central Europe. As expected, there is a marked difference in distribution of CL_nut(N) between various ecosystems (Fig. 6b). Because the critical N concentration for forests is much lower than for the other ecosystems, also the critical leaching rate is much lower. In low precipitation areas, critical loads for all ecosystems are about equal as the sum of the low leaching rate in forest and the N uptake in forests is about equal to the higher leaching rate for grasslands and heath lands. In high precipitation areas, leaching becomes the dominant process and the CL_nut(N) for grasslands and heath land is much higher than that for forest.

Fig. 6

a Fifth percentile critical load CL_nut(N) (eq.ha⁻¹ per year); b cumulative frequency distribution of CL_nut(N) for three vegetation classes (eq.ha⁻¹ per year)

3.4 Sensitivity of Critical Loads

3.4.1 Sensitivity to the Selection of Receptors

To reduce computing time and to eliminate spurious polygons from the map overlay that may be due to map inaccuracies, only receptors of at least 1 km² were included in the modelling. To determine the effect of this cut-off, one simulation was made with all 3.8 million receptors. Results show that only in about 4–5% of the EMEP grid cells the fifth percentiles of CL_max(S) and CL_nut(N) differ by more than 5% from the fifth percentile values based on all 3.2 million data points and that in more than 90% of the cells the difference is less than 1%, thus justifying the limitation to receptors larger than 1 km².

3.4.2 Sensitivity to Criteria

Critical Load of Acidity: CL_max(S)

To test the sensitivity of the model to the chemical criterion used, a number of simulations were made with the criteria listed in Table 1. Results are shown in Table 3 that lists the 5th percentiles and median values of CL_max(S) per region for non-calcareous soils. For most regions the Al/Bc and combined Al&Al/Bc criteria yield about the same critical loads. Exceptions are the fifth percentile critical loads in Russia and Scandinavia which have a low base cation input, and therefore an Al/Bc based critical load can be very low because of low base cation concentrations, even though the associated Al concentration is (far) below the critical value of 0.2 eq m⁻³. In such cases the critical load is determined by the Al concentration alone, resulting in a higher critical load than the Al/Bc-based critical load.

Table 3 5th percentile and median values for CL_max(S) for non-calcareous soils in six regions for various criteria (Central Europe consists of Austria, Belarus, Bulgaria, Czech Republic, Hungary, Moldavia, Poland, Rumania, Slovakia, Switzerland and Ukraine; the CIS except Russia, Ukraine and Belarus comprises Armenia, Azerbaijan, Georgia, Kazakhstan, Kyrgyzstan, Tajikistan, Turkmenistan and Uzbekistan; Scandinavia is defined as Norway, Sweden, Finland and the Baltic states; western Europe consists of Belgium, Germany, Denmark, France, Ireland, Luxemburg, the Netherlands and the United Kingdom; Southern Europe comprises the remaining European countries)

Using a critical Al concentration alone instead of Al/Bc to compute critical loads gives mostly lower critical loads. The only exception is Scandinavia were the fifth percentile is higher than the Al/Bc-based CL because of the very poor soils in these regions: the low Bc input causes very low Al/Bc-based critical loads. Using Al/Bc = 1 to compute the critical load for these sensitive ecosystems will lead in steady state to Al concentrations (far) above 0.2 eq m⁻³ for more than half of the ecosystems, base saturations below 15% and to negative ANC concentrations for almost all grid cells (Fig. 7): an Al/Bc ratio of 1 is mostly equivalent to negative ANC concentrations and low base saturation values (see also Fig. 2).

Fig. 7

Cumulative frequency distributions of grid-average Al/Bc = 1 equivalent values of [Al], Bsat and ANC for sensitive receptors in an EMEP grid cell defined as having a CL_max(S) within 50% of the fifth percentile CL_max(S) of that cell (based on Al/Bc = 1)

Critical loads based on a base saturation of 15% as critical limit are even lower than those computed from a critical Al concentration (Table 3). Especially in Scandinavia the base saturation criterion leads to very low critical loads due to the very low base cation input: in these poor, sandy soils the equilibrium base saturation with no acid input at all would be less than 15%. The average base saturation for the sensitive receptors per EMEP cell equivalent to an ANC concentration of zero is shown in Fig. 8. For Scandinavia and Northern Russia the equivalent base saturation at ANC_le = 0 would be about 10% or lower, due to natural soil acidification processes: aiming at 15% base saturation in such poor acid soils leads to negative critical loads.

Fig. 8

Average base saturation per grid cell equivalent to ANC = 0 for sensitive ecosystems (-)

Critical loads for ANC_le = 0 show strong similarities with those based on the [Al] and Bsat criteria, but as with the Bsat-based CLs, ANC-based CLs can be very low for the most sensitive ecosystems in Russia, Scandinavia and Western Europe. The highest critical loads are computed when aiming at a stable Al_ox pool. This shows that all other criteria protect the soils from loosing their Al buffering capacity.

Critical Load for Eutrophication: CL_nut(N)

If the denitrification rate is low, the critical load of nutrient N mainly consists of N removal through leaching, net growth and immobilization. In the past CL_nut(N) was mostly computed with a critical limit of 0.2–0.4 mg N l⁻¹ which was considered to be representative for forests. This limit was recommended as the concentration in soil solution above which vegetation changes in the under story of forests could occur (UBA 2004). A recent study on critical limits for critical loads of nitrogen (De Vries et al. 2007) revealed that such low limits are mostly related to vegetation changes in Scandinavia, but effects elsewhere in Western Europe probably occur at higher values. Based on literature data and on model results from The Netherlands, the authors therefore suggest that higher limits of about 3–6 mg N l⁻¹ for Western Europe may be used. The choice of the critical concentration strongly determines the relative importance of N removal and N leaching in the critical load. The ratio (N _u + N _i)/N _le as a function of forest growth rate and critical N concentration, [N] _acc , is shown in Fig. 9. In this example Q = 300 mm per year, f _de  = 0, N _i = 1 kg N ha⁻¹ and the N content in stem wood is set at 1.15 g kg⁻¹ representing conifers forest. At [N] _acc  = 0.0143 eq m⁻³ (=0.2 mg N l⁻¹) the ratio between removal and leaching steeply increases from three at a growth rate of 2 m³ ha⁻¹ per year to 12 at a growth rate of 10 m³ ha⁻¹ per year. This illustrates that in areas with high growth rates such as central Europe, N leaching is unimportant compared to N uptake at such low critical N concentrations. For [N]_acc = 0.3 eq m⁻³ (= 4.2 mg N l⁻¹) leaching and removal are almost equal at a growth rate of 10 m³ ha⁻¹ per year, but leaching dominates at lower growth rates. At a growth rate of 2 m³ ha⁻¹ per year, N leaching is five times higher than N removal.

Fig. 9

Ratio (N _u + N _i )/N _le as a function of forest growth rate and critical N concentration

Comparing the cumulative frequency distributions of CL_nut(N) for forests in western Europe, using the critical limits of 0.2 mg N l⁻¹ for conifers forest and 0.3 mg N l⁻¹ for deciduous forests, with CL_nut(N) using the limit of 3 mg N l⁻¹ (specifically suggested for non-Nordic forests) shows an obvious increase in CL_nut(N) with higher concentrations (Fig. 10): the median value over Europe is 930 eq ha⁻¹ per year which is about twice as high as the median value computed with the lower critical concentrations. The ratio between the fifth percentile CL_nut(N) for each EMEP grid cell computed with a critical concentration of 3 mg N l⁻¹ and the 5th percentile CL_nut(N) computed with 0.2–0.3 mg N l⁻¹ is shown in Fig. 10b. Large differences are found in regions with a high precipitation surplus such as Ireland, UK and the northern part of the Alps. In areas with a low leaching rate such as the eastern part of Germany, and southern and central France, the increase in the fifth percentile critical load is less than 50%. This shows that increasing the critical N concentration by a factor of 10 does not always mean a commensurate increase in critical N load.

Fig. 10

a Cumulative frequency distributions of CL_nut(N) for forests in western Europe; b ratio between the fifth percentile CL_nut(N) with different values for [N]_acc

4 Discussion and Conclusions

Combining the latest data bases on soil, land cover, climate and forest growth provided a detailed map with almost 4 million receptors for Europe and Northern Asia suitable for spatially highly disaggregated critical load calculations. The patterns in critical loads show a similarity with those shown in Kuylenstierna et al. (2001), but since their empirical critical loads for acidity are based on soil sensitivity alone (determined by weathering rate), they do not show the influence of the removal of nitrogen or precipitation excess on critical load patterns.

Critical loads are sensitive to the criterion chosen. This study shows that for the most sensitive ecosystems critical loads based on an Al/Bc = 1 or [Al] = 0.2 eq m⁻³ are comparable but that critical loads based on an ANC = 0 are substantially lower. These conclusions can also be drawn by looking at the concept of equivalent criteria: Al/Bc = 1 leads to positive ANC in soils with a very low base cation concentration only. Very low critical loads are also computed when using a critical base saturation of 15%. The same result was found by Holmberg et al. (2001) who computed much higher exceedances of Bsat-based critical loads than of Al/Bc-based critical loads.

Computing Bsat from Al/Bc shows that in the most sensitive ecosystems Al/Bc = 1 will be equivalent to base saturations (far) below 15%. To obtain 15% base saturation in poor sandy soils requires an ANC far above zero which is not very realistic for soils that undergo natural acidification due to leaching of bicarbonates and organic acids (e.g. De Vries and Breeuwsma (1986).

Using ANC = 0 yields critical loads lower than with Al/Bc = 1 but higher than with Bsat = 15%. This seems to contradict the results of Holmberg et al. (2001) who found lowest critical loads when using ANC = 0. However, they included neither bicarbonate leaching nor organic acids in the computation of ANC, leading to an underestimation of actual ANC and very low critical loads. Bicarbonate leaching is important at higher pH values and organic acid leaching is important in soils rich in organic matter (Holmberg et al. 2001); both terms have been added to the ANC calculation (UBA 2004).

This study shows that as an alternative to Al/Bc = 1, using ANC = 0 with bicarbonate and organic acid leaching included is a more realistic option than using a Bsat criterion, as the latter criterion forces naturally acid soils to unrealistically high base saturations and ANC concentrations in the soil solution, thus leading to extremely low critical loads. An advantage of using ANC instead of Al/Bc is that it provides a more general protection against acidification, not based on a weak relationship with forest health alone. However, studies providing limits for ANC leaching related to concrete effects on terrestrial ecosystems are lacking.

Critical loads for nutrient N computed for Western Europe using limits of 3–6 mg N l⁻¹ for the N concentration in soil solution, as suggested in De Vries et al. (2007) for forests in this region, are much higher than the critical loads based on the generally applied limits of 0.2–0.4 mg N l⁻¹, that may be reasonable values for Northern Europe. This sensitivity to the N criterion was already reported by De Vries et al. (1994). Over the whole of western Europe, the median critical load for nutrient N with the revised N limits is about 1,000 eq ha⁻¹ per year (=14 kg N per year). This is in better accordance with the empirical critical loads of 10–15 kg provided by Bobbink et al. (2002) related to vegetation changes in forests than the 500 eq ha⁻¹ per year obtained with the limit of 0.2 mg N l⁻¹. The strong dependence of the critical load on leaching rate does however show the limitations of the SMB model for the assessment of nutrient N critical loads, as it assumes a relationship between nitrate concentration and plant species sensitivity, whereas N availability may be more important (De Vries et al. 2007). Moving to more elaborate models that compute critical N loads based on limits for e.g. pH and N availability related to occurrence of plant species (see e.g. Van Dobben et al. 2006) may improve estimates of nutrient N critical loads.

In this study we have only investigated the uncertainty in critical loads due to the use of different soil chemical criteria. It is clear that the overall uncertainty of the critical loads is also determined by model uncertainty and parameter uncertainty (see, e.g., Zak and Beven 1999; Skeffington et al. 2007). Parameter uncertainty can be substantial and strongly contribute to uncertainty in critical loads as shown by, for example, De Vries et al. (1994). A study for a forested catchment in Belgium showed that using various local and general estimates for nutrient removal by the forest lead to a wide range in critical N loads (Bosman et al. 2001). Others have shown that different methods of estimating weathering rates can yield strongly varying results (Hodson and Langan 1999). Model structure is another important source of uncertainty. The SMB model by definition is a very simple model of reality. The way denitrification is modelled and the fact that we assume a homogenous soil layer contributes to uncertainty. Earlier studies showed that there is a strong effect of the depth at which the chemical criterion should be met especially when using the Al/Bc criterion (De Vries et al. 1994).

Appendix

1.1 Simple Mass Balance (SMB) models

In this Appendix we summarise the critical load models used in this paper. The basic reference is the Mapping Manual (UBA 2004), where also original references are given, but see also Posch and De Vries (1999) and De Vries and Posch (2003).

1.2 Critical Load of Nutrient Nitrogen

Critical loads of N for terrestrial ecosystems can be derived from the N mass balance. Neglecting adsorption, volatilisation and fixation this reads (in eq ha⁻¹ per year):

$$N_{{dep}} = N_{i} + N_{u} + N_{{de}} + N_{{le}} $$

(A1)

The subscript dep refers to total deposition, i to net immobilisation, u to net growth uptake, de to denitrification, and le to leaching. Denitrification is modelled as a fraction of the net input of N (De Vries et al. 1994):

$$N_{de} = f_{de} \cdot \left( {N_{dep} - N_i - N_u } \right)$$

(A2)

where f _de (0 ≤ f _de  < 1) is the denitrification fraction. From Eq. A1 a critical load is obtained by defining a critical (acceptable) limit to the leaching of N being the product of a critical concentration [N] _acc and the water flux Q (m per year). Combining Eqs. A1 and A2 we thus obtain the critical load of nutrient N, CL _nut (N):

$$CL_{{nut}} {\left( N \right)} = N_{i} + N_{u} + Q \cdot {{\left[ N \right]}_{{acc}} } \mathord{\left/ {\vphantom {{{\left[ N \right]}_{{acc}} } {{\left( {1 - f_{{de}} } \right)}}}} \right. \kern-\nulldelimiterspace} {{\left( {1 - f_{{de}} } \right)}}$$

(A3)

In this equation N _i stands for a long-term acceptable N immobilisation and N _u for the long-term average net growth uptake.

1.3 Critical Loads of Acidity

Critical loads of acidity, induced by deposition of N and S, can be derived from the steady-state charge balance for the ions in the soil leachate (in eq ha⁻¹ per year) leaving the root zone (modelled as a single homogeneous layer):

$$H_{{le}} + Al_{{le}} + BC_{{le}} + NH_{{4,le}} = SO_{{4,le}} + NO_{{3,le}} + Cl_{{le}} + HCO_{{3,le}} + Org_{{le}} $$

(A4)

where Org _le is the leaching flux of the sum of organic anions. Neglecting OH⁻ and ${\text{CO}}_3^{2 - }$ (a reasonable assumption even for calcareous soils), the (leaching of) alkalinity or ANC (Acid Neutralizing Capacity) can be defined as:

$$ANC_{{le}} = HCO_{{3,le}} + Org_{{le}} - H_{{le}} - Al_{{le}} $$

(A5)

A steady-state situation with respect to acidification implies a constant pool of exchangeable base cations. Consequently the following mass balance holds for base cations:

$$BC_{{le}} = BC_{{dep}} + BC_{w} - Bc_{u} $$

(A6)

Note, that BC _dep and BC _w include all four base cations (BC = Ca + Mg + K + Na), whereas sodium is not taken up by vegetation (Bc = BC − Na). Since sulphur is regarded a tracer one has:

$$SO_{{4,le}} = S_{{dep}} $$

(A7)

Combining Eqs. A8–A7 and the N balance derived in Eq. A1 (NO _3,le = N _le ) yields for the charge balance (Eq. A4):

$$S_{{dep}} + N_{{dep}} = BC_{{dep}} - Cl_{{dep}} + BC_{w} - Bc_{u} + N_{i} + N_{u} + N_{{de}} - ANC_{{le}} $$

(A8)

It is assumed that there are no sources or sinks of chloride in the soil compartment, and therefore leaching equals deposition. Knowledge of the deposition terms, weathering and net uptake of base cations as well as nitrogen uptake, immobilisation and denitrification allows to calculate the ANC leaching, and thus to assess the acidification status of the soil. Conversely, critical loads of S and N can be computed by defining a critical (or acceptable) ANC leaching, ANC _le,crit , which is set to avoid “harmful effects” on a “sensitive element of the environment” (e.g. damage to fine roots). Using also the equation for the deposition-dependent denitrification (Eq. A2), one obtains for the critical loads of sulphur, CL(S), and acidifying nitrogen, CL(N):

$$CL{\left( S \right)} + {\left( {1 - f_{{de}} } \right)} \cdot CL{\left( N \right)} = BC_{{dep}} - Cl_{{dep}} + BC_{w} - Bc_{u} + {\left( {1 - f_{{de}} } \right)} \cdot {\left( {N_{i} + N_{u} } \right)} - ANC_{{le,crit}} $$

(A9)

Note, that these critical loads of S and N are not unique; every pair of deposition (N _dep ,S _dep ) which fulfils Eq. A9 are critical loads of acidity. However, when comparing S and N deposition to critical loads one has to bear in mind that the N sinks cannot compensate incoming S acidity, i.e. the maximum critical load of sulphur is given by:

$$CL_{{max}} {\left( S \right)} = BC_{{dep}} - Cl_{{dep}} + BC_{w} - Bc_{u} - ANC_{{le,crit}} $$

(A10)

Furthermore, if

$$N_{{dep}} \leqslant N_{i} + N_{u} = CL_{{\min }} {\left( N \right)}$$

(A11)

all deposited N is consumed by uptake and immobilisation, and S can be considered alone. The maximum amount of allowable N deposition (in case of zero S deposition) is given by:

$$CL_{{max}} {\left( N \right)} = CL_{{min}} {\left( N \right)} + {CL_{{max}} {\left( S \right)}} \mathord{\left/ {\vphantom {{CL_{{max}} {\left( S \right)}} {{\left( {1 - f_{{de}} } \right)}}}} \right. \kern-\nulldelimiterspace} {{\left( {1 - f_{{de}} } \right)}}$$

(A12)

1.4 Derivation of the Critical ANC Leaching

To compute acidity critical loads, the critical leaching of ANC has to be specified. With the aid of equilibrium reactions the concentration of Al, bicarbonate and organic acids can be expressed as function of [H] and thus ANC becomes a function of [H] alone. The concentration of Al is modelled by a gibbsite equilibrium:

$${\left[ {Al} \right]} = K_{{gibb}} \cdot {\left[ H \right]}^{3} $$

(A13)

where K _gibb is the gibbsite equilibrium constant. The concentration of HCO₃ is derived from the dissociation of CO₂ according to:

$${\left[ {HCO_{3} } \right]} \cdot {\left[ H \right]} = K_{{HCO_{3} }} \cdot p_{{CO_{2} }} $$

(A14)

where $K_{{HCO_{3} }}$ is the temperature-dependent dissociation constant and $p_{{CO_{2} }}$ the partial pressure of CO₂. The dissociation of organic acids is modelled assuming that only monoprotic organic anions are produced:

$${\left[ {Org} \right]} = \frac{{m \cdot DOC \cdot K_{1} }}{{K_{1} + {\left[ H \right]}}}$$

(A15)

where DOC is the concentration of dissolved organic carbon (molC m⁻³), m is the charge density and K ₁ the dissociation constant. A charge density of m = 0.023 mol molC⁻¹ and pK ₁ = 4.5 was used throughout. The leaching of ANC (Eq. A5) is computed as Q∙[ANC].

Eqs. A13–A15 allow to compute the critical ANC leaching for any critical limit defined in terms of [H] (or [Al]). If a critical molar Al/Bc ratio is defined, the corresponding critical Al concentration is obtained as:

$$[\operatorname{Al} ]_{crit} = 1.5 \cdot (\operatorname{Al} /Bc)_{crit} \cdot [\operatorname{Bc} ]$$

(A16)

with [Bc] = Bc _le /Q; and the factor 1.5 arises from the conversion of moles to equivalents (assuming K as divalent). In addition to soil solution variables, also base saturation has been suggested as a criterion. To link base saturation to soil solution chemistry, we use the Gapon model with simulating the exchange between Al, protons and base cations (Bc), as used in the dynamic soil model SAFE (Warfvinge et al. 1993). Then the following relationship between [H] and base saturation E _Bc can be derived:

$$\left[ H \right] = K_{Gap} \cdot \sqrt {\left[ {Bc} \right]} \cdot \left( {\frac{1}{{E_{Bc} }} - 1} \right)$$

(A17)

with

$$K_{{Gap}} = \frac{1}{{k_{{HBc}} + k_{{AlBc}} \cdot K^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}_{{gibb}} }}$$

(A18)

where k _HBc and k _AlBc are the two selectivity coefficients describing cation exchange.