UIS Commission on Karst Hydrogeology and Speleogenesis
Georg Kaufmann
Republished from: Gabrovšek, F. (Ed.). 2002. Evolution of karst: from prekarst to cessation. Postojna-Ljubljana, Zalozba ZRC, 243-258.
  PDF: /pdf/
Georg Kaufmann

Institute of Geophysics, University of Göttingen, Herzberger Landstrasse 180,
37075 Göttingen, Germany


We present results of a numerical study of karst denudation on limestone plateaux. The landscape evolution model used incorporates not only long-range fluvial processes and short-range hill-slope processes, but also large-scale chemical dissolution of limestone surfaces. The relative efficiencies of fluvial and chemical processes are of equal importance to the landscape evolution of a plateau dropping to sea level along an escarpment. While fluvial processes have an impact confined to river channels, the karst denudation process is more uniform, removing material also from the plateau surface. The combined effect of both processes results in a landscape evolution almost twice as effective as the purely erosional evolution of an insoluble landscape.

Keywords: Landscape evolution, dissolution kinetics, denudation, karst.


Landscape evolution is governed by a balance of forces; on the one hand, vertical tectonic movements resulting from the interaction between lithospheric plates, and on the other, erosion and deposition controlled by a range of processes whose relative importance depend on local climatic conditions, vegetation, and rock type. During the last decade, numerical models simulating landscape evolution from a large-scale, long-term perspective have become increasingly sophisticated. Most surface process models incorporate the effects of short-range processes such as local hillslope diffusion and long-range processes such as fluvial transport (Willgoose et al. 1991, Beaumont et al. 1992, Chase 1992, Howard et al. 1994, Tucker et al. 1994, Kooi et al. 1994, Braun et al. 1997), and they have been applied to the evolution of rifted margins, regions of continental convergence and mountain building. While these models describe landscape evolution satisfactorily in temperate climates and insoluble geological settings, other surface processes of regional importance have so far been largely neglected. For example, glacial erosion has been regarded as an important process in mid- to high-latitude active orogenic regions, which have experienced repeated large-scale glaciations during the last two million years. Recently, Braun et al. (1999) have incorporated a first-order parameterization of ice-bedrock interaction into a large-scale fluvial erosion model. They showed how the interplay between the two processes leads to enhanced rates of surface erosion in mountaineous areas affected by periodic climatic fluctuations.

In the present paper, another process governing landscape evolution in soluble rocks is considered. As surface runoff enriched in carbon dioxide becomes weakly agressive and is able to remove calcite by chemical dissolution, typical karst landscapes evolve, forming steep-walled valleys, enclosed depressions, and, finally, subsurface drainage through caves (e.g. Jennings 1985). While the evolution of subsurface drainage in a karst landscape had been successfully modeled in the past (e.g. Groves et al. 1994, Howard et al. 1995, Clemens et al. 1997, Siemers et al. 1998, Kaufmann et al. 1999, Kaufmann et al. 2000), no attempt has been made so far to parameterize chemical dissolution on the surface of a karst landscape, which is termed karst denudation (e.g. Jennings 1985). These two processes, subsurface drainage through evolving cave systems and surficial karst denudation, are coupled in that the surface discharge on a mature karst landscape quickly disappears underground and drainage is dominated by subsurface flow. As a consequence, valleys become dry, streams disappear at sinks and closed depressions, and surface processes driven by surface flow become less efficient. Herein, attention is focussed on the karst denudation process, while the effect of subsurface drainage and chemical enlargement of fractures within the rock is disregarded. Such a simplified model can thus be regarded as a first step in quantifying the evolution of a soluble landscape.

A karst landscape is often characterized by steep and prominent gorges, which separate limestone plateaux. Compared to insoluble landscapes, gorges occur more frequently and are more rugged. Examples include the Vicos Gorge in Greece, the Verdon Gorge in southern France, the Strickland Gorge in New Guinea, which are very deep and steep and all formed in relatively young mountain ranges, but also the Geikie, Windjana, and Galeru Gorges in Western Australia and the Tarn, Lot, and Dourbie Gorges in the Grandes Causses of southern France, which are prominent features in lower relief landscapes. In all cases, the plateaux dissected by these steep-walled gorges are denudated more evenly.

Measurements of long-term denudation rates on soluble landscapes are often based on the height of pedestals above the surrounding limestone surface. These pedestals are remnants of an older surface, protected by erratic boulders from the last glaciation. It is commonly assumed that the height of the pedestals is a measure of surface lowering since the retreat of the ice at the end of the last glaciation. Jennings (1985) has reviewed literature data on pedestals, and reports denudation rates between 15 and 40 mm/kyr. Short-term denudation rates can be measured by micro-erosion meters (Trudgill et al. 1981). Forti (1984) reports denudation rates from the Triestine karst in Italy, ranging from 20 mm/kyr in regions of relatively low (1442 mm/yr) precipitation, to 30 mm/kyr near Mount Canin characterized by a much higher (2800 mm/yr) precipitation rate. Additional short-term denudation rates reported by Jennings (1985) are around 5-17 mm/kyr on bare rock surfaces.

In the present paper, a simple parameterization of solutional processes is developed and incorporated into an existing surface process model (CASCADE). A series of numerical experiments is then performed in which an originally steep escarpment flanking a flat, high-elevation plateau is allowed to evolve through time. In these experiments, the parameters controlling the efficiency of each process (fluvial incision, hillslope diffusion, and chemical dissolution) are varied to investigate their relative influence on the form of the landscape and the rate at which it evolves.

Table 1 Reference model parameters





Net precipitation

[mm yr-1]


Diffusion constant

[m2 yr-1]


Erosion constant

[m yr-1]


Alluvium length scale



Bedrock length scale



Calcite density

[kg m-3]





CO2 pressure



Dissolution constant

[m yr-1]


Geomorphic model

Large-scale landscape evolution on tectonic timescales is controlled by a number of processes. Short-range hillslope processes, such as weathering, slope wash, mass wasting, and soil creep, redistribute mass over short distances, while long-range flow processes, such as fluvial erosion and sedimentation, and karst denudation, control landscape evolution over long distances. Other processes such as glacial erosion and landsliding are not discussed in the present formulation. The mathematical approach for both short- and long-range processes used in model simulations is outlined below.

Hillslope processes

Hillslope processes are modelled by means of linear downslope diffusion (e.g. Kooi et al. 1994, Braun et al. 1997),


which assumes that the rate of change of topography at nodeis proportional to the second spatial derivative of, is assumed constant and is termed diffusivity. The topography is the sum of bedrock height and sediment cover. Following Beaumont et al. (1992), the diffusivity can be expressed as , where is the thickness of a regolith layer and is the transport velocity of the eroded material. This implies that is a function of both the lithology (small for resistant bedrock, larger for detachable sediments) and climate (increased weathering increases the thickness of the regolith layer). Thus, short-range processes are important on steep landscapes with strong curvature; they reduce the high-frequency content of the topography spectrum. As a consequence, topography is smoothed, slopes decline, drainage divides are eroded and move towards the area of lower slope gradients (Kooi et al. 1994).

Fluvial processes

Erosion and deposition of sediments and bedrock is modeled as a channel-flow process. The discharge (in m3 s-1) at node i is defined as the sum of upstream discharge and local discharge ,

, (2)

where is the net precipitation (precipitation minus evapotranspiration , in m s-1), Ai is the surface area around node i, and is the amount of surface runoff already collected in the upstream catchment of node i. The maximum carrying capacity of the channel between node i and lower neighbour j is given as (e.g. Beaumont et al. 1992, Kooi et al. 1994)

, (3)

where is a dimensionless river erosion constant. Note that the maximum carrying capacity (in m3 s-1) depends both on local slope and discharge , but a river must not necessarily be at maximum carrying capacity. Sediment flux (including suspended sediment and bedload) is calculated as the cumulative sum of sediment transported from the upstream nodes. is compared to the maximum carrying capacity, . In cases where , the river deposits sediments; in cases where , the river incises. The fluvial erosion rate (in m s-1) is given by


In (4), is a characteristic length scale for erosion, which is defined differently for erosion in alluvials() or in bedrock (), is the channel width, with (in ), and as the channel length. For , material is easily eroded, as it is highly detachable. In this situation, the downstream graded river section of the river rapidly propagates upstream towards the headwater of the river, and the ungraded river section rapidly steepens. River profiles are characterized by a short steep upper section and a long low-gradient downstream section. In contrast, characterizes more resistant bedrock and longer erosion timescales. In these situations, grading is inefficient and river profiles are characterized by a more uniform gradient.

In summary, long-range fluvial transport processes result in landscapes in which drainage divides are not eroded and remain stationary, while river profiles steepen in the upstream part of the ungraded section, and slope gradients decline in the downstream graded section.

Table 2 Karst denudation parameters





Absolute temperature


Ion activity



Debye-Hückel coefficient


Debye-Hückel coefficient


Activity coefficient


Activity coefficient


Mass-balance coefficient

[mol l-1]

Mass-balance coefficient

[mol l-1]

Mass-balance coefficient

[mol2 l-2]

Mass-balance coefficient

[mol l-1 atm-1]

Chemical processes

On soluble landscapes such as limestone plateaux, the landscape is also altered by chemical dissolution, as surface water enriched in carbon dioxide is weakly aggressive and is able to dissolve the limestone surface. The denudation rate (in m s-1) can be defined as (e.g. Dreybrodt1988)


where is a dimensionless dissolution constant given by


Here, is the density of calcite (in kg m-3), and [Ca2+]eq the equilibrium concentration of calcium (in mol l-1), given as (Dreybrodt, 1988, pp. 27)


The coefficients , and are experimentally derived mass balance coefficients (e.g. Plummer et al. 1982), which depend on temperature (Table 2). The dimensionless ion activities of calcium and bicarbonate, and , can be derived from the extended Debye-Hückel equation (e.g. Truesdell et al. 1974, Robinson et al. 1955). The carbon dioxide pressure (in atm) may vary over several orders of magnitude in nature, from 10-3.5 atm over bare rock, 10-2.5 atm in temperate climate soils, to 10-1.5 atm in tropical soils (e.g. White, 1984). Hence, depends on climate (temperature), and strongly on lithology (CO2-pressure in atmosphere and soil).

Note that (5) depends on discharge only. The factor 40 is necessary to transform [Ca2+]eq from mol l-1 to kg m-3. In Fig. 1, denudation rates as a function of discharge are shown for three different temperatures and carbon dioxide pressures. It is observed that the carbon dioxide pressure is the dominant factor controlling denudation, while changes in temperature are less important.

Fig 1. Denudation rate as a function of discharge for different partial carbon dioxide pressures (in atm) and water temperatures (in oC).


All numerical experiments start with identical initial topographies. The model domain is a rectangular area of 100-km side length, discretized into an irregular grid of 2500 nodal points, with an average nodal spacing of 2 km. An initial topography is chosen resembling a high plateau with an average height of 1000 m, and dropping to sea level (0 m) along its southern boundary over a distance of 10 km. This drop in topography simulates an escarpment created by a recent continental rifting event. The plateau itself rises slightly towards the northern boundary to ensure that initial drainage is towards the south and all surface runoff drains through the escarpment. Both northern and southern boundaries have fixed elevation, but sediment is allowed to exit across these boundaries.

Reference experiments

Two reference runs based on fluvial erosion and diffusion alone where initiated to establish reference landforms to which soluble landscapes can be compared. Both the modelled topography after 1 Myr of evolution and a characteristic example of a river profile through time were examined in order to better define differences in morphology between the evolving landscapes. Model parameters are listed in Table 1 for reference run 1 shown in Fig. 2. To characterize each experiment, following Kooi et al. (1994), the ratio of diffusive to fluvial efficiency is given as


For run 1, R1 = 20 m; this low value of R1 indicates that long-range fluvial processes dominate the landscape evolution, which can be seen in the modeled topography after 1 Myr (Fig. 2a), dominated by eight steep rivers draining the plateau. All rivers follow fairly linear courses, and the valleys are characterized by steep sidewalls and are relatively narrow. The initial escarpment is left unchanged except where it has been incised by the rivers. A typical river profile is shown in Fig. 2b and is characterized by a convex shape near the escarpment edge and a relatively linear, less steep section near the base. This indicates that diffusive processes control the evolution in the ungraded river section. Fluvial erosion is responsible for the efficient removal of sediments, and the resulting graded river section is less steep and more linear. The escarpment retreats fast along the main rivers, but the headwaters of the rivers on the plateau remain at a constant height.

Results of the reference run 2 are shown in Fig. 3. Here, the diffusion constant is increased by a factor of 100 to m2 yr-1 to simulate a wetter climate. With a diffusive to fluvial efficiency ratio of R1=2000 m, diffusion is much more efficient in this run. As a consequence, the simulated topography is much smoother after 1 Myr. The short-range diffusional processes have filtered out the high-frequency components of the topography, especially across the escarpment. Rivers flow in wide valleys with gentle sidewalls. The ungraded river section extends much further downstream, as can be seen in the typical river profiles redrawn (Fig. 3b). The convex shape of the ungraded river section extents halfway down the escarpment, and the shorter graded river section is characterized by a steep gradient. Escarpment retreat is less effective, when compared to reference run 1.

Notably, the river erosion on the plateau itself is not very effective for both reference runs, as the erosion rate (4) depends on the product of channel slope and discharge. Hence, the lack of relief on the top of the plateau inhibits erosion.

Fig 2. Modeled topography after 1 Myr. The evolution is based of fluvial erosion and diffusion, parameters are given in table 1. The thick line on top of the topography represents the course of the river, whose profile is shown in the top panel as a function of time, for subsequent time steps 0.2, 0.4, 0.6, 0.8, and 1 Ma.

Fig 3. Same as Fig. 2, but based on fluvial erosion and increased hillslope diffusion.

Karst denudation experiments

Next, the additional effect on landscape evolution introduced by chemical dissolution of limestone is examined. According to (5), removal of material by chemical dissolution only depends on discharge, not on slope gradient (as for fluvial erosion), nor on curvature (as for hillslope processes). Therefore, karst denudation is also an efficient process on the gently sloping plateau area, where surface runoff increases in downstream direction to the south, and along the base of the escarpment. However, drainage divides cannot be eroded away by fluvial processes, or any process that is dependent on discharge. As a consequence, a different landform evolution is expected on a soluble landscape. A second ratio is introduced, relating diffusive efficiency to the sum of fluvial and chemical efficiency,


Results from the first of the two model runs are shown in Fig. 4. A low value of carbon dioxide pressure of PCO2=10-3.5 atm is assumed that is characteristic of barren karst surfaces. All other parameter values are as defined in Tables 1 and 2. For the given parameters, the dissolution constant is comparable to the fluvial constant ; therefore both processes are similarly effective. Hence, the efficiency ratio is reduced to R2=12 m (compared to R1=20 m for the reference run 1). Long-range processes are therefore more effective. This can be seen after 1 Myr of landscape evolution, by which time the plateau is deeply dissected by eight major rivers (Fig. 4b). Most of the rivers have wide flat-bottomed valleys in their downstream graded section, whose sidewalls are generally steeper than in the reference run 1. In parts, the old plateau surface is reduced to narrow necks between deep river gorges. Major differences in river profile geometry can be observed (Fig. 4a).

In contrast to the reference runs on insoluble landscapes, the ungraded upstream river section now is significantly eroded. Valley incision starts close to the headwaters along the northern end of the plateau. The ungraded river section is characterized by a concave gradient similar to gradients resulting from diffusion processes. Gradients in the downstream graded river sections are low and convex to linear in shape, and thus dominated by fluvial processes. With time, the graded river section is cutting back into the plateau and gradients decline, while the ungraded river section is steepening. However, in contrast to reference run 1, incision is significant also on the plateau itself; thus, the total vertical elevation drop across the escarpment is reduced with time.

In the next run (Fig. 5), the carbon dioxide pressure is increased to PCO2=10-2.5 atm, simulating a temperate climate soil cover. Surface runoff now percolates slowly through the soil and is enriched in carbon dioxide, thus becoming more aggressive. Consequently, chemical dissolution is enhanced, as can be seen in the larger value of . The efficiency ratio is further reduced to R2=8 m; thus, diffusion is even less efficient. Surface morphology after 1 Myr is similar to the previous run, with eight major rivers draining the plateau to the south (Fig. 5a). However, river valleys are much more incised, resulting in deep narrow gorges reaching far upstream. The effect of an increased carbon dioxide pressure can be seen more clearly in the chosen set of river profiles (Fig. 5a). While the general characteristics -- concave steepening gradients in the ungraded section and flat convex to linear gradients in the graded section – are similar to the previous run, valley incision has progressed much more. In fact, after 1 Myr the graded low-gradient river section occupies more than 80 percent of the river, leaving only a short but very steep ungraded river section close to the main drainage divide in the north.

Fig 4. Same as Fig. 2, but based on karst denudation, fluvial erosion and diffusion. The partial carbon dioxide pressure is 10-3.5 atm

Fig 5. Same as Fig. 2, but based on karst denudation, fluvial erosion and diffusion. The partial carbon dioxide pressure is 10-2.5 atm.

Hypsometric curves

The effect of denudation on landscape evolution can also be seen in a statistical quantity, the hypsometric curve (Fig. 6). This quantity is derived from the topographical elevations of the model domain by summing up the elevations for a particular elevation range and relating it to the total elevation. As it can be seen in Fig. 6a, the initial topography for all models is characterized by a peak in the hypsometric curve around the 1000 m elevation, which represents the initial plateau topography. As the escarpment is fairly steep, the initial hypsometric curve for elevations between sea level and 1000 m is almost uniform and orders of magnitude lower than the peak at 1000 m altitude. At sea level, a second peak in the hypsometric curve represents the outwash plain of the plateau.

For the reference models shown in Figs. 2 and 3, the final hypsometric curves are characterized by a uniform increase of elevation between sea level and 1000 m elevation, which represents the incised river sections. The plateau surface itself has been reduced in size, as the peak is slightly smaller, but the elevation has been kept, as no erosion is possible on the plateau. No significant differences between the low and high diffusion landscapes is visible.

However, for the two models with karst denudation shown in Figs. 4 and 5, the final hypsometric curves are different: The peaks for the plateaux has moved downward and decreased significantly, as the denudation is able to lower the plateaux height itself. For larger effects of karst denudation (higher partial carbon dioxide pressure), the peak has moved further downwards and the hypsometric curve has become more uniform.

Fig 6. Hypsometric curves for the initial (solid line) and the final (grey areas) models.


The numerical models presented herein are a first attempt to incorporate chemical processes into landscape models. Thus, only a qualitative comparison with topographical data can be made without statistical testing.

The modeled topography of the runs including karst denudation is similar to a low-relief plateau dissected by narrow river gorges such as the Grandes Causses in southern France, where the deep, steeply incised gorges of the Tarn, Lot, and Dourbie contrast with the relatively flat segments of the plateau of the Causses in between (Fig. 7). The evolving river valleys can have a significant sinuosity, and often only small remnants of the original plateau surface separate two gorges as narrow necks. As the karst denudation process is as effective as fluvial river erosion in eroding a landscape, evolution on a limestone plateau will downgrade the landscape surface twice as fast as on a similar sandstone plateau. The major difference in the karst denudation and the fluvial erosion processes are visible on the plateau itself: while fluvial erosion is less effective on the flat plateau, as fluvial erosion rate depends on the product of channel slope and river discharge, karst denudation can effectively remove material on the plateau itself, as denudation rate depends only on river discharge. Hence, karst denudation is an efficient process for surface lowering on plateau landscapes.

Fig 7. Topography of the Grandes Causses in southern France, with the deep, narrow gorges of the rivers Lot, Tarn and Dourbie, and the Causse de Sauveterre (CS), Causse Mejean (CM), and Causse Noir (CN).


The DEM data have been kindly provided by Francis Lucazeau from Universite de Montpellier II. The figures in this paper are drawn using the GMT graphics package (Wessel et al. 1991, Wessel et al. 1998).


  • Beaumont C., Fullsack P. and Hamilton J. 1992. Erosional control of active compressional orogens. In: McClay K. R. (Ed.), Thrust Tectonics. New York: Chapman and Hall, 1-18.
  • Braun J. and Sambridge M. 1997. Modelling landscape evolution on geological time scales: a new method based on irregular spatial discretisation. Basin Res. 9, 27-52.
  • Braun J., Zwartz D. and Tomkin J. H. 1999. A new surface processes model combining glacial and fluvial erosion. Ann. Glaciol. 28, 282-290.
  • Chase C. G. 1992. Fluvial landsculpturing and the fractal dimension of topography. Geomorphology 5, 39-57.
  • Clemens T., Hückinghaus D., Sauter M., Liedl R. and Teutsch G. 1997. Modelling the genesis of karst aquifer systems using a coupled reactive network model. In: Hard Rock Hydrosciences, Proceedings of Rabat Symp. S2., IAHS Publ., vol. 241, Colorado, 3-10.
  • Dreybrodt W. 1988. Processes in Karst Systems. Berlin: Springer.
  • Forti F. 1984. Messungen des Karstabtrages in der Region Friaul-J ulisch-Venetien (Italien). Die Höhle 35, 135-139.
  • Groves C. G. and Howard A. D. 1994. Early development of karst systems 1. Preferential flow path enlargement under laminar flow. Water Resour. Res. 30 (10), 2837-2846.
  • Howard A. D. and Groves C. G. 1995. Early development of karst systems 2. Turbulent flow. Water Resour. Res. 31 (1), 19-26.
  • Howard A. D., Dietrich W. E. and Seidl M. A. 1994. Modeling fluvial erosion on regional continental scales. J. Geophys. Res. 99 (B7), 13971-13986.
  • Jennings J. N. 1985. Karst geomorphology. Oxford: Blackwell.
  • Kaufmann G. and Braun J. 1999. Karst aquifer evolution in fractured rocks. Water Resources Res. 35 (11), 3223-3238.
  • Kaufmann G. and Braun J. 2000. Karst aquifer evolution in fractured, porous rocks. Water Resources Res. 36 (6), 1381-1392.
  • Kooi H. and Beaumont C. 1994. Escarpment evolution on high-elevation rifted margins: Insights derived from a surface processes model that combines diffusion, advection, and reaction. J. Geophys. Res. 99 (B6), 12191-12209.
  • Plummer L. N. and Busenberg E. 1982. The solubilities of calcite, aragonite and vaterite in CO2-H2O solutions between 0 and 90oC, and an evaluation of the aqueous model of the system CaCO3-CO2-H2O. Geochim. Cosmochim. Acta 46, 1011-1040.
  • Robinson R. A. and Stokes R. H. 1955. Electrolyte solutions. London: Butterworths Sci. Publ.
  • Siemers J. and Dreybrodt W. 1998. Early development of karst aquifers on percolation networks of fractures in limestone. Water Resour. Res. 34 (3), 409-419.
  • Trudgill S. T., High C. J. and Hana F. K. 1981. Improvements in the micro-erosion meter. Brit. Cave Res. Group Tech. Bull. 29, 3-17.
  • Truesdell A. H. and Jones B. F. 1974. WATEQ, a computer program for calculating chemical equilibria of natural waters. U.S. Geol. Surv. J. Res. 2 (2), 233-248.
  • Tucker G. E. and Slingerland R. L. 1994. Erosional dynamics, flexural isostasy, and long-lived escarpments: A numerical modeling study. J. Geophys. Res. 99, 12229-12243.
  • Wessel P. and Smith W. H. F. 1991. Free software helps map and display data. EOS 72, 441-446.
  • Wessel P. and Smith W. H. F. 1998. New, improved version of generic mapping tools released. EOS 79, 579.
  • White W. B. 1984. Rate processes: chemical kinetics and karst landform development. In: La Fleur R. G. (Ed.), Groundwater as a geomorphic agent. London, Boston, Sydney: Allen and Unwin, 227-248.
  • Willgoose G., Bras R. L. and Rodriguez-Iturbe I. 1991. A physically based coupled channel network growth and hillslope evolution model. 2. Applications. Water Resour. Res. 27 (7), 1671-1684.