Climate change due to anthropogenic emissions of greenhouse gases is predicted to increase the average surface temperature. The most evident soil changes in the Alps will occur in proglacial areas where already existing young soils (with an age in most cases of up to 150 years) will continuously develop and new soils will form due to glacier retreat. Based on existing soil chronosequence data and statistical analyses in the proglacial area Morteratsch (Switzerland), the present-day state of the soils as well as their future development in the next 100 years in the existing and new proglacial area has been modeled taking the retreat of the glacier into consideration. The present-day as well as the future soil distribution was modeled using a probabilistic approach. Several soil characteristics have been modeled such as the pH value, the skeleton content, and the soil depth relevant to plant growth. To model soil properties in a future proglacial area (that is now covered by ice), the glacier-bed morphology had to be modeled. The calculations were performed using the cubic Non-Uniform Rational B-Spline (NURBS) curve to parameterize the course of a branch in flow direction. With the help of the ice cap and relief factor the thickness of the glacier was modeled. Climate change was introduced numerically by changing the mass balance of the glacier. For the area of interest a temperature increase of 1.6°C by the year 2050 and 3°C by the year 2100 can be assumed (according to the scenario A1B of IPCC). In the upper part of the proglacial area mostly Skeletic/Lithic Leptosols and Humi-Skeletic Leptosols will be found. In flat parts close to the main river, additional Fluvisols will develop. A considerable part of the upper proglacial area does not have any soil cover. Lithic/Skeletic to Humi-Skeletic Leptosols are modeled on the young lateral moraines. Chronosequences were vital to make any (4D) predictions of soil evolution in the proglacial area. The statistically and probabilistically based model also had, however, its weaknesses. The problems are related to the sediment properties in the glacier bed and the stability of new moraines.
Easily recognizable traces of dramatic climatic variations make high mountain areas unique geotopes and “storytellers” about past as well as potential future climate change effects on landscape dynamics and living conditions in regions of rugged topography. Long-term observations of glaciers have provided convincing evidence of rapid global climatic change; the worldwide retreat of mountain glaciers during the 20th century was striking (Haeberli et al., 1999; Meier and Bahr, 1996). The apparent homogeneity of the signal at the secular time scale, however, contrasts with the large variability at local/regional scales and over time scales of years to decades (Letréguilly and Reynaud, 1990). Based on an increased theoretical and empirical knowledge of processes and feedback effects, it has been possible also to predict possible climate changes on a regional scale (IPCC, 2001; OcCC, 2002). Direct and dramatic ecological responses to this impending warming are expected (Peters and Lovejoy, 1992), in the form of feedbacks that could modify transfer rates of energy, water, and trace greenhouse gases at the earth's surface (Rosenberg et al., 1983). Landscapes may respond very noticeably and differentially to climate change as they integrate all ecological and historical factors (Theurillat et al., 1998). A key or “interface” function must thereby be attributed to soils. Soil-landscape patterns result from the integration of short- and long-term pedogeomorphic processes (Friedrich, 1996; Klingl, 1996). Despite numerous studies on the effects of climate warming on single processes, little is known about the reaction of a whole soil ecosystem (Rustad et al., 2001).
Many of the soil properties change continuously with time. The soil, however, can only be measured at a finite number of places and times, and any statement concerning the soil at other places or times involves prediction. Variation in soil is so complex that no description of it can be complete, and so prediction is inevitably uncertain (Heuvelink and Webster, 2001). Only appropriate observation will provide the basic knowledge necessary for assessing the development in reality, and statistically calibrated numerical spatial/gray-box models (rather than sophisticated deterministic models as applied in detailed scientific process studies) must help in the prediction of consequences and possible future scenarios. The results of Mendonça Santos et al. (2000) and Herrmann et al. (2001) demonstrate the utility of GIS technology for the facilitation of data-set management, for spatialization, analysis, and visualization.
Anthropogenic emission of greenhouse gases is predicted to increase the earth's average surface temperature during the next 50–100 years. For the area of interest (proglacial area of Morteratsch, Switzerland), a temperature increase of +1.6°C by the year 2050 and +3°C by the year 2100 can be assumed (according to the scenario A1B of IPCC, 2000). Therefore, additional areas will become ice-free and subject to weathering and soil formation. The most evident soil changes in the Alps will occur in proglacial areas where already-existing young soils will continuously develop and new soils will form due to the glacier retreat. The rate of the reactions that are of fundamental interest in the understanding of the soil system and its interaction with the environmental surrounding conditions has been deduced by a chronosequence in this proglacial area (Egli et al., 2006 [this issue]).
The main aim of this investigation was to predict future soil development for the next 100 years in the existing and new proglacial area associated with the retreat of the glacier Morteratsch and using the previous findings from chronosequences and statistical analyses (Egli et al., 2006 [this issue]).
The studied soils lie within the proglacial area (1900–2150 m a.s.l., with a present-day area of 1.8 km2) of Morteratsch in the Upper Engadine (Switzerland). The investigation area is described in more detail in Egli et al. (2006 [this issue]).
Modeling of Soil Properties
Spatial modeling was performed with ArcGIS 8.3 (ESRI) with modules programmed in Visual Basic for Applications (VBA). Input data sets were the digital soil map, the glacial states, the digital terrain model (DTM; raster of 20 m) within the proglacial area, the digital elevation after the subtraction of the glacier surfaces (Biegger, 2004; see below), and simulated extents of the Morteratsch glacier for the years 2050 and 2100 (Biegger, 2004). The calculations were done raster-based (GRID; 20 m resolution).
Soil types were evaluated using the statistical trends of the individual soils as a function of time, landscape form, exposure, and slope according to Egli et al. (2006 [this issue]) and to Tables 1 and 2.
Regressions between topographical features and time for the soil type Skeletic/Lithic Leptosol.
Regressions between topographical features and time for the soil type Humi-Skeletic Leptosol.
The regressions are a mathematical expression of the probability of a certain soil type as a function of time and a given topographic feature. Equations 1–3 show how the probabilities of a certain soil type on a specific grid with certain topographic characteristics were calculated. The probabilities of a soil type were calculated as a function of the slope WS, the exposure WE, and shape of the landscape WL. These probabilities are, furthermore, related to the time of the soil formation. The semi-empirical factors a, b, and c are derived from the regression equations.
The investigated factors are independent. To calculate the probability of a soil type (WT) with a certain age at a specific grid, the probabilities according to Equations 1–3 have to be multiplied.
The minimum probability for the occurrence of a soil type has been determined in an iterative process (learning process; cf. also Behrens et al., 2005), i.e. an optimization of the modeled soil types with the soil map. The output of the model finally produces values ranging from 0 to 1. The closer the values are to 1, the more likely is the occurrence of the considered soil type. The classification returning the highest prediction accuracy is used for the final prediction.
The soil properties are strongly dependent on the soil type, shape of the landscape, topography, and age. A more detailed description of the modeled soil properties is given in Table 3. The modeling of the soil properties, consequently, was carried out heuristically. The modeled soil types served as a base that was related to the soil age and topography to derive the individual soil properties. The calculation was according to the entity-relationship principle (Klingl, 1996). The implementation was done by a query of the databank using Boolean functions such as “and,” “or,” “not” as well as logical operators such as “lower than,” “greater than,” “equal,” “not equal” to get the corresponding properties (Klingl, 1996). The pH values in the topsoil (surface soil layer) and subsoil, the soil depth relevant for plant growth, and the skeleton content in the top- and subsoil could be modeled with this procedure (see Table 4, which lists the model structure for the calculation of the specific soil properties).
Modeled soil properties.
Modeling of soil properties (entity relationships containing Boolean and logical functions) derived from soil types, soil age, and topography.
Soil cartography and classification (also of soil properties) was made according to the FAL system (Brunner et al., 1997). Soil types are translated into the WRB (FAO, 1998) system. The dataflow is given in Figure 1. Some specific properties were modeled according to the FAO-UNESCO system (1990) as the WRB system does not fully meet the requirements for a soil classification in Alpine areas; e.g. the soil type “Ranker” indicates an intermediate stage between Humi-Skeletic Leptosols and Dystric Cambisols (endoskeletic).
Modeling the Surface Area Underneath the Glacier
The present topographical situation will change with a further retreat of the Morteratsch glacier. To model soil properties in a future proglacial area (that is presently covered by ice), it is essential to have or infer information about the surface underneath the glacier because one important model input is the topography (see above). To calculate this topography the glacier was regarded as the sum of several surfaces, each of them representing an individual branch. A cubic Non-Uniform Rational B-Spline (NURBS) (Piegl and Tiller, 1997) curve was applied to parameterize the course of a branch in flow direction (Biegger, 2004). A NURBS surface S(u,v) of degree p in the u direction and of degree q in the v direction is a piecewise rational functionHaeberli et al. (1999). The glacier thickness along the flow-line Cf can be expressed as a combination of the ice cap factor kx and the relief factor rx. The maximum thickness dmax of the profile would result at the equilibrium line (ELA = equilibrium line altitude of glaciers). The ice cap factor can be written as (Biegger, 2004) Biegger, 2004). To geometrically model the branch surface, the longitudinal profile had to be extended with nx cross sections. Although a glacier may consist of more than one branch, each branch had to be treated separately. Each branch runs from the source point Ps to the terminal point Pt. In order to compute the glacier bed, the surface S(u,v) was discretized by projecting each grid point of the original DEM onto the branch surface and subtracting the thickness from the DEM at each grid point. For the deformation of a glacier due to climate change, each branch was considered individually and its reaction was parameterized as proposed by Haeberli and Hoelze (1995). Climate change was introduced numerically by changing the mass balance of the glacier. Given a step-wise temperature change ΔTAir assuming a certain climate scenario, the change Δb of the mass balance could be calculated by Kuhn, 1990). Experimental studies have shown a good accordance with historical measurements applying a value of 160 m/°C for the Morteratsch glacier. The experimentally evaluated ratio of Δb/Δh was 0.68 (Biegger, 2004).
The spatial distribution of the simulated glacier thickness is visualized in Figure 2. It can be seen that there are two regions of increased ice thickness. Both regions are located in a plane just below steep slopes where the thickness is less than 100 m. The applied spatial sampling was 20 m according to the horizontal resolution of the DTM. The surface area with and without the Morteratsch glacier (glacier bed) is plotted in Figure 3. The future evolution of the Morteratsch glacier distinctly depends on the climate scenario. The sensitivity study by Gyalistras (2000) for estimating the future trends of air temperature and precipitation in the Alps showed that air temperature will possibly rise more than the global average. A probable climate scenario could be the one according to the scenario A1B of IPCC (2000) with a temperature increase of +1.6°C by the year 2050 and +3°C by the year 2100. The impact of the applied climate scenario on the Morteratsch glacier is visualized in Figure 4. The most obvious changes can be recognized at the tongue, whereas the accumulation area remains more or less unaffected. In the year 2100 the branches formerly contributing to the Morteratsch glacier will become isolated. Within the next 100 years a much greater area than the present-day proglacial area will become ice-free where soils will be able to form. According to the surface area calculations, new proglacial lakes are not supposed to be formed in that time span.
The basis for soil modeling included the time since deglaciation and topographic features. The slope, shape, and exposure of the landscape were calculated using the DTM25 (with 20 m resolution) and the modeled surfaces according to Biegger (2004). Additionally, the water net and the flow direction (modeled surfaces) were taken into account. The soil model calculates in the present proglacial area larger areas having Skeletic/Lithic Leptosols after about 30 to 50 years. Already after 90–100 years of soil evolution a transition (during a time span of 20–40 years) of Skeletic/Lithic Leptosols into Humi-Skeletic Leptosols is modeled which agrees well with the soil map. Along the main river system, the distribution of the Fluvisols generally agrees well with reality. A lower agreement was, however, found along small river branches. A comparison between the modeled area with soil cover and mapped soil area gave an agreement of about 74% (see also Fig. 5). Problems arose especially where the sediment bed of the glacier was extremely thin or even absent.
Several soil characteristics have been modeled and compared with the soil map. Among the modeled characteristics are the pH value (CaCl2) in the topsoil (uppermost, surface layer) and in the subsoil, skeleton content in the top- and subsoil, and the soil depth relevant for plant growth (according to the FAL system: soil depth = soil volume – skeleton volume – groundwater volume; result is related to depth instead of volume). In the correctly modeled soil area, the agreement of the modeled characteristics with the mapped ones was in all cases very high (Table 5); only very small differences to the soil mapping could be discerned.
Comparison between modeled and mapped soil properties (values refer to the correct modeled soil area).
The stream network of the future proglacial area had to be calculated using the modeled surface area under the glacier. Together with the existing stream network the future situation was derived. The stream network was a basis for the calculation of the distribution of Fluvisols.
During the first 25 years of soil formation only patches having a significant soil cover could be found or are modeled. Thus, a certain area in front of the glacier tongue remains without any soil cover. The area with soil cover then steadily increases with Skeletic/Lithic Leptosols that appear first and transform later into Humi-Skeletic Leptosols (Fig. 6). The older Humi-Skeletic Leptosols show signs of an initial B horizon formation (here called Ranker; see Fig. 6). This soil type will increasingly appear in the year 2050. Ranker will first be formed in rather flat areas with a soil age of >150 years. The eastern part of the investigation sites should develop slightly more quickly due to having more north-facing slopes. After 200 years (maximum age in the year 2050) of soil development, nearly 96% of the glacier bed will be covered by a (thin) soil layer. A major problem in soil modeling is the consideration of the lateral moraines. They are usually very unstable and disturb soil formation. The periphery of the lateral moraines has to be considered as an area more for potential than for effective soil formation. Ranker will be the dominant soil type in the front of the proglacial area at the end of the 21st century (Fig. 6). In the upper part of the proglacial area mostly Skeletic/Lithic Leptosols and Humi-Skeletic Leptosols will be found. In flat parts close to the main river, additional Fluvisols should develop. A considerable part of the upper proglacial area does not have any soil cover. Lithic/Skeletic to Humi-Skeletic Leptosols are modeled on the young lateral moraines. Any soil development there will depend on the stability of these moraines ( = potential areas of soil formation).
There was no model for the glacier tongue states available for the year 2150. Assuming that the glaciers will not have readvanced to the original position of the year 1850 ( = lowest part of the proglacial area), which is highly improbable, soil evolution can be calculated at least for the lower part of the proglacial area. If we assume that the glaciers will retreat also in the 22nd century, then a soil development can also be modeled for the upper part of the proglacial area (Fig. 6). After 300 years of soil evolution (lower part of the proglacial area, year 2150), Dystric Cambisols (endoskeletic) will probably appear. Because this soil type has not been found in the present-day proglacial area (after 150 years of soil evolution), its first appearance was derived from chronosequences in the central Alps (Egli et al., 2001). This makes the prediction for the year 2150 even more speculative and was, therefore, only done for the soil types.
Discussion and Conclusion
Soil evolution will relatively rapidly follow the glacier retreat. However, it takes about 25 years until significant signs of soil formation will be observable. Within the next 100–150 years the soils in the present-day proglacial area will thicken; the skeleton content will be reduced by being subjected to physical and chemical weathering (and also because of the accumulation of organic matter; Fig. 7), will continuously acidify (Fig. 8), and will progressively melanize (darken) with time and distance from the ice. Soil types will grade finally into Rankers (Humi-Skeletic Leptosols with traces of a B horizon) and most probably not before the year 2100 into Dystric Cambisols (endoskeletic). Soil evolution, although quick in the Alps, seems to be slower when compared to similar sites in Norway where Regosols were observed already after 47 years of deglaciation which then graded into Brunisols (Cambisols) in less than 120 years (Mellor, 1985, 1986). Frost disturbance as well as slope instabilities can, however, impede or delay such a progression of soil types (Haugland, 2004). This can lead to highly patterned features of the spatial soil distribution associated with an abrupt threshold between genesis and stabilization. As a consequence, modeling of soils on a small scale in existing and future proglacial areas can hardly be very precise.
Modeling alpine soil properties can be done using several methods. Huber (1994) made a process-oriented modeling of soils with GIS in the Bavarian Alps using topographic indications, comparisons with existing soil maps, and the formulation of geo-ecological interrelationships based on the principles of Leser and Klink (1988). A further possibility is the derivation of soil relevant data using remote-sensing techniques (e.g. Mikhaylov, 1990; Gauthier and Tabbagh, 1994). If enough point data are available, then geostatistical principles can also be used for predicting soil properties (e.g. Ahn et al., 1999; Grunwald et al., 2000; Zhu et al., 2001; Lark, 2003). The applicability to larger regions of detailed soil surveys in small reference areas was successfully assessed by Lagacherie and Voltz (2000). The mapping method consisted of a soil classification in a reference area. The probabilities of occurrence of soil classes at a site were used to predict soil properties at unvisited sites using digital data sets such as the DTM, geology, and hydrology. Gessler et al. (2000) were able to account for between 52 and 88% of soil property variance using easily computed terrain variables such as slope and flow accumulation. Empirical models developed from DTM were also successful in predicting horizon depths of the topsoil (Thompson et al., 2001). Promising results in mapping soil properties are also obtained from fuzzy functions (Hannemann, 2005) or neuronal networks (Behrens et al., 2005).
One drawback of all these mapping techniques is the absence of the time scale that was important for the case study presented. Chronosequences and therefore the relationship of soil properties to time as a function of topographic properties are vital in making any (4D) predictions of soil evolution in proglacial areas. The statistically and probabilistically based model, however, also has its weaknesses. The modeling of the present-day situation of soil distribution gave “only” an agreement of about 73% with the soil map, although very significant correlations of the individual soil types with topographic features and time could be found. One major problem occurs because the sediment bed of the glacier was partially extremely thin or even absent. In such a case, soil modeling failed. Another problem is related to the lateral moraines that may be unstable and consequently hinder soil development. Erosion or small debris flows from such moraines occur randomly and were not really predictable. Finally, the interaction of geomorphodynamics and soil formation makes soil modeling difficult in high Alpine environments as described in Egli et al. (2006 [this issue]).
Nevertheless, the modeling of soils will, for instance, allow a more precise prediction of vegetation changes and thus enable a more comprehensive prediction of future landscape evolution in a rapidly changing high Alpine environment.
This research was supported by grants of the National Research Programme 48 “Landscapes and Habitats of the Alps” (Swiss National Foundation), project number 4048-064352. We are, furthermore, indebted to Aldo Mirabella for his helpful comments on an earlier version of the manuscript.