A deductive model of karst evolution based on hydrological probability

zyxwv EARTH SURFACE PROCESSES AND LANDFORMS, VOL. 13, 271-288 (1988) A DEDUCTIVE MODEL OF KARST EVOLUTION BASED O N HYDROLOGICAL PROBABILITY C. C. SMART zyxwvutsr Department of Geography, University of Western Ontario. London, Ontario, N 6 A 5C2, Canada Received 2 December 1986 Revised 30 March 1987 zyxwvu zyxwvutsrq ABSTRACT Combination of a conduit flow law with the exceedance probability of stream discharge allows the estimation of the exceedance probability of water level in a simple conduit aquifer. The probability of water levels higher than available relief is interpreted as the probability of surface as opposed to underground runoff. Very high probability of surface runoff implies a fluvial environment, whereas very low probabilities define a mature karst or ‘holokarst’. Intermediate probabilities identify ‘fluviokarst’. Overflow probability depends o n available relief, mean discharge, and especially conduit radius. Growth rate of the underground conduit depends on saturation deficit which thus controls the rate of evolution of the fluvial landscape, through fluviokarst t o holokarst. However, variations of discharge and sedimentation through time can cause dramatic reversion of karst drainage into less mature states. Landscapes experiencing such periodic rejuvenation will have a confused morphology. A functional definition of landscape may be more objective and pertinent than arbitrary interpretation of form. KEY WORDS Karst Landscape evolution Modelling INTRODUCTION Deductive models of landscape evolution developed by Davis (1909) have long fallen into disfavour, largely due to an unacceptable predominance of intuition over factual data concerning landscape processes. Both theoretical and empirical understanding of geomorphic processes has improved over the past 30 years and renewed attempts at evolutionary modelling are now possible (e.g. Kirkby, 1978). Karst is particularly amenable to modelling because the dominant solutional processes are reasonably well understood (e.g. White, 1984).The major gaps in understanding are related to hydrology and geological control which are responsible for much of the local idiosyncracy of karst (e.g. Jennings, 1985). There are two major types of karst landscape: holokarst and fluviokarst. Holokarst is characterized by pure depression forms, principally dolines and the much larger poljes. There are very few stream channels, and little if any surface water. Fluviokarst areas exhiba’intermittent surface streams and segments of perennial rivers, and contain elements of both karst and fluvial forms such as blind and dry valleys. Most authors consider there to be an evolutionary progression from an initially fluvial landscape, through fluviokarst to holokarst (e.g. Cvijic, 1918; Roglic, 1964),although Grund (1914) has suggested that holokarst develops directly from an exposed limestone surface. Miller (1981) has presented strong evidence for the former case through the successful reconstruction of palaeovalley networks through areas of tropical holokarst. There is unresolved debate concerning the influence of climate on karst. Distinct morphoclimatic styles are no longer widely recognized (e.g. Ford, 1980), although solutional processes may be closely controlled by zyxwvutsrqp zyxwvu 0 197-9337/88/03027 1 - 18$09.00 0 1988 by John Wiley & Sons, Ltd. 272 zyxwvutsrqpon zyxwvutsr C. C . SMART climate (e.g. Drake, 1980). Climatic change is considered to influence both runoff and sedimentation in karst (Jennings, 1985). Sawicki (1909) pointed out that periodic sedimentation can cause a reversion of holokarst to fluviokarst through the plugging of cave systems. This suggests that a karst ‘cycle of erosion’ is inherently reversible, unlike the classical Davisian Model. Rather than basing an evaluation of these notions on geomorphological intuition, or the interpretation of form, it is intended here to develop a simple hydrological model of a karst aquifer, and to apply this model to landscape evolution. First, evolution under steady state conditions will be considered, then climatic change will be investigated in terms of runoff, solution potential, and sedimentation. THE AQUIFER MODEL The evolution of karst reflects the pattern of solutional erosion which is directed by the hydrological characteristics of the karst aquifer. Water in karst aquifers occurs both in discrete solutional conduits and in a more widely distributed diffuse form. Most karst aquifers consist of complex admixtures of these two components, but for the sake of simplicity, the pure conduit aquifer will be considered in this case. Such aquifers generally develop in the most massive carbonates, and characteristically receive recharge from discrete sinking streams. Groundwater is conveyed through conduits to a resurgence, and there is little hydrological interaction with the host carbonate. Most aquifer models are developed to provide accurate hydrological forecasts, and so are evaluated in terms of the real variance explained by the model. The most successful models are generally stochastic, but they are seldom physically realistic, and may not be interpreted physically. The model presented here is not a forecasting model, but is suitable for straightforward numerical experiments designed to explore aquifer behaviour. A single, circular section, water-filled, horizontal conduit is taken to run between an inflow shaft and a resurgence (Figure 1). Water entering the shaft causes a rise in head, driving water through the conduit to the resurgence. Flow through the system is governed by the Darcy-Weisbach equation: zyxwvutsr zyxwvut zyxwvut where Q is system discharge. r is the radius of the circular section conduit. n is a constant, 3.14159. y is the acceleration due to gravity. f is the Darcy-Weisbach friction factor. L is the length of the conduit. Qin zyxwvutsrqp , Figure 1. Cross-section through the conduit model. Flow through the system is controlled by h which varies instantaneously so that Qm Qmt zyxwvutsrqpon zyxwvuts zyxwvutsr 273 A DEDUCTIVE MODEL OF KARST EVOLUTION For present purposes,fand L are taken to be constant at 0.33 and 1000 m respectively. A natural conduit may be truncated by stream incision, or extended by underground capture. Such complications are largely ignored in this analysis. Some evolution in the friction factor would be expected as it is dependent on relative roughness of the conduit and to some extent on flow regime (e.g. Lauritzen et al., 1985). Extremely high values for the friction factor determined for karst conduits probably result from unrecognized constrictions (e.g. Atkinson et al., 1983). Length, roughness, and conduit radius are all form variables, but the overwhelming significance of the latter allows length and roughness to be considered constant as a first approximation. To allow a straightforward adoption to the probability domain, it is assumed that inflow to the conduit is immediately matched by outflow. This implies that changes in head occur instantaneously, and that the volume (and length) of the system is constant regardless of discharge. Smart (in press) has explored the constraints of this assumption which is not critical to the present model. Inspection of equation 1 shows that in the short term most variation in discharge will be attained through changes in head. A simplified operational equation then becomes: zyxwv zyxwvutsrqpo zyxwvutsrq zyxwv zyxw 0.5 where c' = r2'57c($) . Over long periods, r becomes a variable rather than a parameter, giving a long-term operational equation of the form: Q = c r 2 5 h0-5 (3) where c = c ' r - 2 5 . The exceedance probability of discharge expresses the likelihood of a reference discharge ( Q )being equalled or exceeded by actual discharge (q), providing that q is an independent random variable. This provides a suitable means of describing natural flow variability. The exponential distribution provides a convenient expression for the exceedance probability of flow: P(q 3 Q ) = ~ X (-1.Q) P (4) where the parameter 1. is equal to 1/Q and Q is the mean discharge. This is a reasonable approximation for many smaller rivers, or rivers with a sustained low flow season and a near zero minimum discharge (Smart, in press). Equation ( 2 ) may be substituted for Q on the right hand side of (4),providing the exceedance probability distribution of head within the aquifer. Thus: p ( h 2 H ) = exp ( - Ac'h05) (5) This relationship applies over the short term, and over the long term as the conduit enlarges, it may be replaced by: p ( h 3!* H ) = exp ( - Acr23h05) (6) I Figure 2 shows the exceedance probability distribution of head for several values of r for an aquifer with a mean flow discharge of 1 m3 s- I , a conduit length of 1000 m and a friction factor of 0.33. The conduit radius has a profound effect on the head distribution. When conduits are small, enormous values of head are necessary to drive commonly occurring flood discharges through the system. For example, for a conduit radius of 0.5 m, and a mean discharge of 1 m3 s - ',equation (6)implies that a discharge of 5.9 m3 s - will on average be equalled or exceeded on one day per year ( p = 0.0027). The head associated with this discharge is about 950 m, giving a hydraulic gradient of 1.0. Such extreme conditions are unlikely to be realized in nature, primarily because of insufficient relief. Other constraints on the actual head developed may be the absorption of the flood impulse by groundwater storage (both diffuse and conduit types), and the activation of intermittent conduits which changes the geometry of the system (Smart, in press). Here, however, only the inadequate relief effect will be considered. ' 274 C. C. SMART 0 500 zyxwvuts zyxwvutsr zyxwv zyxwvut Head (m) Figure 2. The exceedance probabilityof hea# on a conduit aquifer for conduit ra&i:0.251.5 rn. The probabilityof head being over 100 rn is given for a range of conduits in Table I zy zyxwvut zyxwvutsrq THE PROBABILITY OF SURFACE OVERFLOW The maximum head which can be developed on a karst conduit is given by the level which water can reach before spilling off over the surface. The probability of this head being exceeded is therefore given by p(h 3 H,) = exp (- kr2'5H,05) (7) where H a is the maximum possible head, or available relief. p(h 3 H a ) is the probability of surface overflow, p ( 0 F ) . A very high probability of surface overflow implies that surface rivers will be virtually perennial, a very low probability implies that surface flows will be very rare or nonexistent. This implies that karst may be described in terms of its hydrological function rather than its appearance, with the probability of surface overflow providing a simple index of behaviour. A high probability of surface flow indicates the predominance of fluvial processes, while very low probability of surface flow typifies a holokarst. Intermediate probabilities represent fluviokarst. In this formulation, there is clearly no discrete boundary delimiting the karst types, but for convenience in discussion the fluvial-fluviokarst transition will be arbitrarily taken at p = 0.1, (or 36.5 days per year), while the transition to holokarst will be taken at p = lo-' (or 365 days per 1,OOO years). For example, equation (7) has been evaluated for 1 km long conduits of various radii, an available relief of 100 m and a mean discharge of 1 m3 s- (Figure 2). The probability of overflow, and the inferred landscape are summarized in Table I. The table shows that by the arbitrary criteria assigned above, it is only those areas Table I. Probability of overflow and inferred hydrological state for modelaquifers with 100 mqygiilablereliefand 1 m3s - ' mean discharge Radius (m) 0.01 0.10 0.25 0.50 0.75 1 .oo 1.50 2430 3.00 Probability of Overflow Landscape 1.00 x 100 9.66 x lo-' 7.13 x lo-' 1.47 x lo-' 5.11 x 10-3 1.97 x 10-5 1.08 x 243 x lo-*' 4.58 x Fluvial Fluvial F1uvia1 Fluvial Fluviokarst Fluviokarst Holokarst Holokarst Holokarst zyxwvutsrqpon zyxwvutsr A DEDUCTIVE MODEL OF KARST EVOLUTION 275 drained by conduits somewhat larger than 0.5 m radius which become karstic. Large conduits are necessary to provide sufficient drainage to maintain a holokarst, although this also depends on other factors as will be seen shortly. Karst does not develop spontaneously according to this model. Rather, it requires the evolution of effective underground drainage to produce true karst. While groundwater conduits are initiated, the landscape will be moulded largely by fluvial processes, providing the initial permeability of the host carbonate is low. The latter assumption is probably justified for the massive carbonates likely to support pure conduit aquifers. This result supports the fluvial-fluviokarst-holokarst evolutionary sequence outlined above. It follows that if conduit size controls the hydrological function of the landscape, then conduit growth rates will control the time taken to pass through the sequence. The next section focusses on the assignment of a time base to landscape evolution. THE RATE OF LANDSCAPE EVOLUTION The rate of solutional conduit growth in karst aquifers is a complex function of hydrology and hydrochemistry, and has been approached through considerations of hydraulics, chemical thermodynamics and kinetics, and also direct measurement (e.g. Palmer, 1981, 1984; White, 1977, 1984). It is assumed that a simple C02-H20-CaC03 system applies here. Other processes of conduit growth such as corrasion or the action of organic and sulphuric acids are not considered. The most sustained yet least understood phase in the evolution of a karst aquifer is the initiation of flow routes from tight primary fractures. White and Longyear (1962), Atkinson (1968), and White (1977) have discussed this problem and consider that a conduit radius of about 5 mm marks the point at which conventional development begins. For present purposes it is assumed that the conduit has already reached this critical size and has fully developed turbulent flow. The zero origin on time axes is therefore rather arbitrary, but marks the initiation of the karst aquifer from primary fractures. Palmer (198 1, 1984) has applied present understanding of the solution process to the problem of conduit growth. As a first approximation, solution rate can be considered approximately independent of conduit radius, discharge, and hydraulic gradient, providing that: 1. waters are at more than 50 per cent of saturation throughout the aquifer,and 2. the ratio 2Q/rL is greater than m/s. The water passing through the conduit is unlikely to be less than 50 per cent saturated unless the sinking stream flows from a non-carbonate, which will be assumed here. As a worst case for the second condition, ifa conduit radius of 3 m (the largest considered here, and a rough upper limit to a structurally stable form) is considered to be 1000 m in length, then the assumption remains reasonable down to discharges of 0.015 m3/s. If the discharge is exponentially distributed with a mean of 1 m3/s then 99 per cent of flows will be greater than this limit. Accepting these assumptions implies that saturation deficit is the critical control on conduit growth rate and this permits adoption o f a simple solution rate equation basedentirely on the saturation deficit (Palmer, 1981). zyxwvuts zyxwvut zyxwv 6r = K (AC)2 here 6r is the growth rate of the conduit radius in m a - l K is a constant (3.51 x lop8 m7 g-’ a) derived from the reaction constant (0.0948 m4g-’ a) divided by rock density (2.7 x lo6 gm-3). At greater than 90 per cent saturation K decreases dramatically (White, 1977) providing a further constraint on the present analysis. AC is the saturation deficit, the difference between actual solution concentration of C a C 0 3 and saturation concentration in g m - 3 (or mgl-’). The conduit radius (r)at time t may be simply defined as an integral function of the saturation deficit over time: 276 zyxwvutsrqp zyxwvutsr zyxwvut zyxwv C . C . SMART Figure 3 shows how the saturation deficit controls the rate of conduit growth, using the number of years required for a conduit to grow by one metre. The range runs over 5 orders of magnitude, encompassing periods of both relative environmental stability and dramatic geomorphic change. It is unlikely that the higher saturation deficits plotted in Figure 3 would be found in sinking streams and here an upper limit of 40 mgl- is employed. Figure 4a shows the growth curves derived from equation (9) for conduits with constant saturation deficits of 10,20 and 40 mgl- '. The first and second reach 0.7 m and 2.8 m radius respectively after 200,000 years, while the last reaches 3 m radius and likely structural instability in 53,000 years. These rates are compatible with ages inferred from direct measurement (e.g. Coward, 1971)and indirect dating control (e.g. Gascoyne, 1981) and geomorphological inference (e.g. Tratman, 1969). The following discussion will use this model, first to consider the influence of available relief, saturation deficit, and mean discharge over time. Then cyclic variations of discharge, saturation deficit, and sedimentation will be considered, and finally a more complex model incorporating all of these effects will be investigated. STEADY STATE EVOLUTION The probability of surface overflow, p ( O F ) , provides a convenient index of the functional order of a karst landscape, and it is primarily dependent on the available relief, the mean discharge of the sinking stream, the conduit radius, and other form variables. The conduit radius increases through time as a function of the saturation deficit. The governing equation under steady conditions is obtained by modifying equation (9) for steady state conditions and substituting into equation (7). zyx zyxwvu p ( O F ) = exp ( - I C ( K A C ~ ~ ) * ' ~ H : ~ ) (10) The controls on landscape evolution will be considered in terms of this equation. Available relief Assuming a saturation deficit of 20 mg 1 and a mean discharge of 1 m3 s - ' various values of H , may be substituted into equation (10). Figure 4b shows the evolution of the probability of surface overflow with available relief of 10, 100 and 1000 m. Table 11 gives the time required to reach fluviokarst and holokarst status. This is also indicated by the horizontal dotted lines in the probability plots. Despite the 2 orders of magnitude contrast in available relief, the low relief area takes only 2.5 times longer to reach holokarst maturity. Furthermore, the time differences might be masked in nature by enhanced mass movement in high relief areas which may constrain conduit development. Sustained surface overflow will cause stream incision, reducing H , over time and slightly increasing the probability of surface overflow. If surface erosion rates are based on equation (9) the effect is trivial except for small H,. Faster surface erosion rates might be expected, however, due to the higher discharges anticipated in ~ ~ ~ " " ' " " " 100" " ' ' ' ~ 200 Saturation Deficit, mg/l CaCO3 Figure 3. The minimum time required to enlarge a conduit by 1 m as a function of saturation deficit zyxwvutsrqpo zyxwvu zyxwvut zyxwvutsrq 211 A DEDUCTIVE MODEL OF KARST EVOLUTION A ln 3- 0 200 100 Time (ka) Time (ka) 0 zyxwvuts 200 100 Time (ka) 0 U=lOO -6 0 100 Time (ka) zyxwvu 200 Figure 4. (a) The growth rate of conduits for saturation deficits of 1420, and 40 mg 1- ';(b) The probability of overflow against time for available relief of 10, 100 and lo00 m. Saturation deficit is 20 mgl-' and mean discharge is 1 m3 s-'. The transitions between fluvial-fluviokarst and fluviokarst to holokarst are marked by horizontal dotted lines. See also Table II; (c) The probability of overflow against time for saturation deficit of 10, 20, and 40 mgl-'. Available relief is 100 m and mean discharge is 1 m3s - I . See also Table II; (d) The probability of overflow against time for mean discharge of 0.01,0.1, 1.0, 10, and 100 m3 s - '. Available relief is 100 m and the saturation deficit is 20 mgl-'. See also Table I1 zyxwvutsr zyxwvuts zyxwvutsr 278 C. C . SMART Table 11. Parameters of models used in Figure 4 and times of transitions between fluvial, fluviokarst, and holokarst. A indicates the respective amplitude of cycles of mean discharge and saturation deficit Parameter Fig. 4b Line 6) (ii) (iii) 4c 4d A 0 0 0 (ii) (iii) 0 0 0 (i) (ii) (iii) (iv) (v) 0 0 0 0 0 (1) Time of transition (ka) Fluvial-Fluviokarst-Holokarst Q A AC Hu 1 1 1 0 0 0 20 20 20 10 100 1000 61 39 25 106 67 43 1 1 1 0 0 0 10 20 40 100 100 100 154 39 10 267 67 17 001 01 1.0 10 100 0 0 0 0 0 20 20 20 20 20 100 100 100 7 16 39 97 243 12 30 73 184 46 1 zyxw zyxwvuts 100 100 overflows (Smart, in press) and correspondingly greater sediment transport. Incision of surface channels may also shorten the conduit at its outlet at a rate roughly inversely proportional to the surface gradient ( H J L ) , thus decreasing the probability of overflow. These two effects completely cancel one another out, providing they follow the rules outlined. Under more complex conditions either may predominate. Because of the low sensitivity of overflow to available relief and incision, H , , will be held constant at 100 m in all subsequent analyses. Saturation deficit Different values of the saturation deficit (AC) have been substituted into (lo),and Figure 4c shows dramatic differences in the rate of karst evolution for values of 10,20, and 40 mgl-'. Mean discharge is 1 m3 s - and available relief is 100 m. The low saturation deficit model only evolves into fluviokarst after 154 ka (Table II), whereas the stream with large solution potential becomes fluviokarstic in 10 ka, and holokarstic in 17 ka. These results may be compared to karst features in western Ireland. The outlets of extremely long caves draining limestone catchments have little if any solution potential (e.g. the Fergus River Cave, Tratman, 1969),and they show very little development and flood frequently despite their great age (P. Mills, personal communication). In contrast, acidic streams draining directly from shale catchments have developed significant caves in 10 ka of post-glacial time (e.g. Cullaun I, Tratman, 1969). However, these respective systems are also distal and proximal extremities and reflect downstream changes in solution potential not considered in the present model. The saturation deficit is clearly of great significance in determining the rate of karst evolution. ' Mean discharge Large karst features are often associated with large streams, yet large streams will be associated with a high probability of surface overflow. Mean discharges covering four orders ofmagnitude have been substituted into 1, in (10) and plotted in Figure 4d. The times of evolution are givon in Table 11. The discharge exceedance probability distribution is assumed exponential for all streams, and the saturation deficit is constant 20 mg 1It is clear that large streams take far longer to evolve into karst than small trickles. This finding demands some exploration. The imprint of a particular stream on a landscape reflects the magnitude of the stream and the duration of flow, i.e. the overall amount of geomorphic work accomplished. The largest rivers will remain on the surface, developing significant valleys, and possibly cutting through the available relief to frustrate development of '. zy zyxwvutsrqpon zyxwvuts A DEDUCTIVE MODEL OF KARST E V O L U T I O N 279 underground conduits. At the extreme they may produce the deeply incised gorges which are a common feature of karst areas (e.g. Sweeting, 1973, p. 104). Intermediate size streams develop valley forms of significant size before karst capture, whereas the smallest streams sink underground early in their development and will not leave significant fluvial forms. Classical karst depressions are essentially surface catchments with corresponding surface stream systems (Williams, 1972). The spatial character of a karst landscape is controlled by the size of the formative river system which depends on the extent of the sustainable surface catchment. The extent of a surface catchment in karst is determined by the spacing of penetrable fissures. The most dense frequency of sinkpoints will have correspondingly minute surface catchments and will fail to develop any surface fluvial forms. Widely distributed solutional lowering will occur and inherited landscape elements may be preserved. Examples of such landscapes might be those developed on the highly permeable Cretaceous chalk (e.g. The Chalk Downs of northwestern Europe), or on Pleistocene carbonates in the Caribbean. It is doubtful if these areas would be described as true ‘karst’. More massive limestone will support larger depressions, and larger surface streams. These in turn will overflow more frequently, leaving a stronger fluvial imprint on the resulting landscape. Non-limestone (or allogenic) catchments have no ‘penetrable fissures’,and are usually associated with the larger sinking rivers, blind valleys, and gorges, i.e. the most fluvial karst forms. The karst landscape can therefore be viewed as one permitting discrete underground drainage or surface catchments. The character of the landscape reflects the size ofthe catchment elements as well as the stage in overall evolution of underground drainage. The simple steady state model has provided a systematic. approach to understanding karst landscape evolution. The size and hydrochemistry of karst streams appear to be of greater significance than available relief in the evolution of karst landscapes. However, the mean discharge of a stream, its hydrochemistry and sediment load are variable over time, especially in response to climatic change. This provides the focus for the next section. zyx zyx zyxwv zyxwvu UNSTEADY CONDITIONS AND CLIMATIC CHANGE Karst landscapes are considered to evolve over periods of from lo4 to lo5 years (e.g. Atkinson et al., 1978). Throughout the Tertiary such timescales have been associated with profound climatic variations. However, the influence of climatic change on karst has received less attention than regional climatic differences. Rather than considering the dramatic impact of glaciers on karst, attention will focus on modelling more moderate, but less well understood effects caused by changes in runoff, temperature, and sediment loads. The general equation (10) developed for steady state conditions needs a little modification to give the probability of overflow at time t . The mean discharge and thus (1,) will become a function of time. The size of the conduit reflects the cumulative growth under variable saturation deficit and must be modelled by equation (9). The only other variable considered here is sediment accumulation S [ t ] . This is modelled as a simple reduction in conduit radius, regardless of effects on the cross-sectional form or the friction factor. The general equation has the form: p ( ~ ~ ) [=t ]exp (-- i . [ t ] ~ r [ t ] ~ . ~ ~ , O . ~ ) (1 1) where r[t] = (k j:I:AC[Tl2 1 d T - S[t] S[t] is the reduction in radius due to sedimentation at time t . Functions defining [t], C [ t ] , and S[r] will be described below as each factor is considered in turn. Runoff variations The mean discharge of rivers is widely accepted as variable over hundreds of years. Rather than attempting a reconstruction of global palaeohydrology, a simple sinusoidal hydrograph may be employed (Figure 5a) providing a moderately varying flow far below the extremes claimed by Dury (e.g. 1964): 280 zyxwvu zyxwv zyxwvut zyxwvutsrqp zyx zyxwvut c. c. SMART a, 0 0 100 200 Time (ka) Time (ka) "................... ..._ x=10 zyxwv .. . ... . .. C, 100 200 Time (ka) Figure 5 (a) Model variation of mean discharge over 200 ka; (b) The growth rate of condults under saturation deficits of 10, 15, and 40 mgl-'; (c) The probability of overflow against tlme for variable discharge and saturatlon deficits of 10, 15, and 40 mg I - ' . Available relief is 100 m See also Table 111 where Q is overall mean discharge (lm3s - ' ) . A is amplitude (0.5 m3 s - '). t is time in years. and P is period of variation (50 ka). The steady state model may be modulated by substitution of equation (12) into (1 1).This function is plotted in Figure 5b for steady saturation deficits of 10,15, and 40Bmg 1- and S [ t ] = 0. Table Ill summarizes the times of hydrological transition. High saturation deficits allow karst to develop so rapidly that climatic changes in runoff are of negligible importance. The conduit with intermediate saturation deficit exhibits some fluctuations in overflow probability within the 200 ka window, briefly adopting a fluviokarst habit after 25 ka of holokarst evolution. The stream with low saturation deficit varies between fluvial and fluviokarst in Figure 5b. Table 111 shows repeated modulation in the likelihood of surface runoff for this stream for 300 ka because the rate ofgrowth of the conduit is slow in comparison to the periodicity of runoff. This variability indicates periodic manifestation of less mature behaviour in fluviokarst and holokarst. Relatively minor flow variation has a profound effect on the probability of surface overflow. In effect, the karst has been rejuvenated, reverting in function, but not necessarily in form. Such landscapes will often be significantly decoupled from their hydrological function. Furthermore, the convenient concept of irreversible evolution has been refuted with only a very minor modulation of climatic conditions. I ' zyxwvuts zyxwvu zyx zyxwvuts zyxwv zyxwvuts 281 A DEDUCTIVE MODEL OF KARST EVOLUTION Table 111. Parameters of models used in Figure 5 and times of transitions between fluvial, fluviokarst, and holokarst. Underlined figures indicate reverse evolution Parameter (m3s - ') Fig. 5 Line (i) Time of transition AC Q fmg 1 ') ~ A Q A 0.5 1.o 0 rc 10 (ka) Fluvial-Fluviokarst-Holokarst 122 133 161 224 229 266 - 287 310 5 (ii) 0.5 1.o 0 15 65 117 142 152 - 5 (iii) 0.5 1.o 0 40 11 17 zyxw Variations in saturation deficit Despite the lack of headway made in understanding the influence of climatic change on karst evolution, karst processes are known to be profoundly influenced by climate (Smith and Atkinson, 1976).Temperature is especially important because of its influence on biological production of C 0 2 (Drake, 1980). It is difficult to synchronize saturation deficit fluctuations with the discharge cycles, because the influence of runoff on saturation deficit is not well understood in a long-term or spatial sense, although short-term dischargehardness effects are known (e.g. drought conditions cause a fall in biological productivity, and floods can either flush out saturated pre-flood water, or dilute pre-storm flows). These observations are not helpful in determining long-term variations of saturation deficit in a karst stream, however. Here it is assumed that the previously modelled pluvial episodes are associated with declining temperatures and increased flushing, causing reduction of soil C 0 2 .The outcome is that the saturation deficit is inversely correlated with discharge and has the same periodicity. This relationship is quite possibly spurious, nevertheless, three levels of fluctuation are modelled (Figure 6a) according to zyxw zyxwvuts ( 6"') A C [ t ] = z + A Cos n+- where c is the long-term average saturatien deficit, and takes values of 10, 15, and 30 mgl-'. A is the amplitude of the 50 ka cycles, and has respective values of 5, 10, and 20 mgl- '. The conduit growth curves are given by substitution of equation (13) into equation (9)and they are shown in Figure 6b. The slope of the line is strongly influenced by saturation deficit. Substitution of equation (13) into equation (1 1) and holding mean discharge constant at 1 m3 s-' gives the sequence of overflow probability shown in Figure 6c and the threshold times listed in Table IV. The major difference between this model and that with a constant saturation deficit (Figure 4c, Table 11) is the periodic retardation and acceleration of development caused by the saturation deficit cycles. The effect is most significant in delaying evolution of the stream with the highest saturation deficit. Otherwise, apart from slight modulation, there is little contrast in evolution. There are two reasons for the lack of apparent impact. First, variations in the saturation deficit cannot increase the probability of surface overflow, unlike the discharge cycles outlined above, or sedimentation cycles below. Second, there is an association between 282 zyxwvutsrqpo zyxwvuts zyx C . C . SMART a 100 -0 200 Time (ka) b /------ zyxwvu ............................ 100 Time (ka) 0 100 Time (ka) .................... ................ 200 zyxwvu 200 zyxwvutsr zyxwvu zyxwvutsrqpon Figure 6. (a) Model cycles of saturation deficit against time; (b) The growth rate of conduits under variable regimes of saturation deficit; (c) The probability of overflow against time for variable saturation deficits. Mean discharge is 1 m3 s- ' and available relief is 100 m. See also Table IV Table IV. Parameters of models used in Figure 6 and times of transitions between fluvial, fluviokarst, and holokarst Parameter Q (m3s-') Fig. 6 6 6 Line (i) (ii) (iii) A 0 0 0 Q 1 .o 1.o 1.o AC (mg I A@-'. ~ '1Ac Time of transition (ka) Fluvial-Fluviokarst-Holokarst 5 10 10 13 1 266 15 65 115 20 30 20 26 amplitude and average saturation deficits. High saturation deficits might cause strong periodicity but the associated systems evolve so rapidly that they are not vulnerable. Lower saturation deficits sustain evolution, but are insufficiently varied to generate more than minor oscillation about the long term mean evolution path. In the absence of a reversion mechanism, karst evolution will be progressive, sequential, and irreversible. zy zyxwvutsrqpon zyxwvutsrq 283 A DEDUCTIVE MODEL OF KARST EVOLUTION Periodic sedimentation One ofthe most widely recognized geomorphic effects of climatic change is a variation in sediment yield (e.g. Schumm, 1969). Excess sediment has a repressive effect on karst processes. This may be partly due to an absorption of solution potential by carbonate sediments, but is largely due to the plugging of sinkholes and underground conduits. Only the last effect will be considered here. The principal cause of sedimentation is karst climatic disturbance, but volcanic ash falls (e.g. Smith, 1975)and human disturbance (e.g. Sweeting, 1973, p. 297; Drew, 1983) may also generate sporadic episodes of infilling. Sawicki (1909) suggested that periodic sedimentation occurred naturally throughout karst evolution. Sediment accumulation in caves results in a reduction in passage radius and often an increase in the friction factor. Here, sedimentation is modelled through a simple progressive reduction in passage radius, no friction effects are considered. Following sedimentation, a passage may be abandoned, or eventually flushed out to continue its evolution, as is modelled here. Although the ages of sedimentation events in caves have been established, there is no generalized chronology, some events are catastrophic, others gradual. In the present model, events are arbitrarily triggered at 35 ka and at 100 ka intervals thereafter. Maximum infilling occurs 5 ka after initiation with a 1 m reduction in passage radius. Sediment is removed exponentially over the subsequent 20 ka. The sedimentation events are plotted in Figure 7a and each has the form: zyxwvuts S [ t ] = 0.3075t" exp (- 0408 t') (14) Time (ka) zyxwvu zyxwvut /.----- ........?.?.??.?. ......... b, 100 Time (ka) 0 ---. -6 0 ........ I $1' I I . ................................ 200 ..... .. , I 100 Time (ka) 200 Figure 7. (a) Model periodic sedimentation events. Sedimentation is modelled as a reduction in conduit radius; (b) The growth rate of conduits with periodic sedimentation. Saturation deficit is constant at 10,15, and 35 mg I I , mean discharge is 1 m3 s - and available relief is 100 m;(c) The probability of overflow against time with periodic sedimentation. See also Table V ~ 284 zyxwvutsrqp zyxwvutsr zyxw C . C . SMART where t' is the time from initiation. Solutional enlargement is continued throughout the infilling episodes, giving the predicted evolution of passage radius over time shown in Figure 7b under three regimes of saturation deficit (10,15, and 35 mgl-'). Clearly, the impact on the large conduit is far less than on the smaller conduits which are closed completely. This is somewhat misleading, because the degree of infilling in nature may be largely determined by the local conduit form and catchment morphology so that great size may not necessarily prevent complete occlusion. Equation 14 may be substituted into equation (11) to give the overflow probability plotted in Figure 7c and the landscape transitions listed in Table V. Mean discharge is constant at 1 m3 s - '. It is clear that sediment infilling can have a dramatic, but short term impact under the model conditions. Smaller conduits are closed by sedimentation and even holokarst systems may revert to fluvial behaviour. Although the larger conduits are not closed, they show a significant reversion in function, until they are sufficiently large to have immunity from the infilling episodes. In contrast, the slowest growing conduit shows four reversions in 540 ka (Table V). Within the range considered, saturation deficit cycles, discharge periodicity, and sedimentation events are of increasing significance in disrupting the steady rate of karst development. Only the latter two are capable of causing rejuvenation, while saturation deficit determines how rapidly the system becomes independent of discharge and sedimentation cycles. Unfortunately, the sensitivity of the system to sedimentation means that errors in identifying its character may have a profound effect on system evolution. Despite these concerns, the three 'climatic' variables can be integrated into a single model. AN INTEGRATED MODEL Under natural circumstances the fluctuations in discharge, saturation deficit, and sediment loading will interact in a more or less complex manner. Assuming a common dependence on climatic change, then the phase zyxwvu zyxwvu zyxwv zyx Table V. Parameters of model in Figure 7 and times of transitions between fluvial, fluviokarst, and holokarst. Underlined figures indicate reverse evolution. S indicates periodic sedimentation Parameter Q A Q AC (mg 1 - I L A AC 0 1.0 0 (m3S C ' ) Fig. Line 10 Time of transition (ka) Fluvial-Fluviokarst-Holokarst S 158 237 248 292 338 344 - 336 349 437 - 445 / 539 542 I .- 7 (ii) 0 1.o 0 15 S 69 130 137 145 136 - - 150 2 39 242 - 7 (iii) 0 1.o 0 35 s 13 24 38 43 zyxwvutsrqpon zyxwvut 285 A DEDUCTIVE MODEL OF KARST EVOLUTION relations between these variables can be assigned. Combining Figures 5a, 6a and 7a, the overall time series can be pictured as in Figure 8a. The sediment events are associated with incipient pluvials (e.g. Schumm, 1969),and the saturation deficit is out of phase with discharge by pi radians. For clarity, only one out of the three levels of saturation deficit used (as portrayed in Figure 6a) is illustrated in Figure 8a. Figure 8 has been extended to 500 ka to show that the conduit growth rate is repeatedly modulated by fluctuations in saturation deficit and punctuated by episodes of sedimentation. Substitution of equations (12), (13) and (14) into equation (11) produces the pattern of karst evolution plotted in Figure 8c and listed in Table VI. All three model scenarios produce reversions in evolution, many of them forcing holokarst right back to the fluvial condition. The most powerful reversion process is sedimentation, and this effect is amplified by in-phase discharge variations. In isolation, the saturation deficit is unable to cause reversion, but it does control the conduit growth rate and so becomes the principal control on the long-term sensitivity to reversion. zyxwvu DISCUSSION The insight provided by this model must not be stretched beyond the limits of the simplifying assumptions. These assumptions are of several distinct orders. First, selected variables have been artificially constrained for computational simplicity. For example, the friction factor will vary with discharge and passage radius over time, and the conduit will vary in length and volume with discharge. Such assumptions will distort. but not materially affect the implications of the model. 0 100 zyxwvutsrqp zyxwvu 200 300 400 500 300 400 500 300 400 500 Time (ka) -0 100 200 Time (ka) 0 100 200 Time (ka) Figure 8. (a) A 500 ka synthesis of periodic controlling variables; discharge, saturation deficit, and sedimentation.Only one saturation deficit curve is drawn to preserve clarity; (b) The growth rate of conduits under variable saturation deficit, discharge, and sedimentation regimes; (c) The probability of overflow against time for conditions in (a) and (b). See also Table VI 286 zyxwvutsrqpo zyxwvutsr C. C. S M A R T zyxwvuts zyxwvu zyxwvuts Table VI. Parameters of models used in Figure 8 and times of transitions between fluvial, fluviokarst, and holokarst. Underlined figures indicate reverse evolution. S indicates periodic sedimentation Para meter Q AC ( m g l -') A Q A Ac Time of transition (ka) Fluvial-Fluviokarst-Holokarst 0.5 1 .o 5 10 S (m's '1 Fig. 8 Line (i) 119 136 159 - 237 241 220 235 264 296 304 - 339 343 337 350 438 445 - 8 (ii) 0.5 1.0 10 15 S 64 78 81 115 143 137 149 19 23 139 8 (iii) 0.5 1.0 20 35 s - 38 44 - Other assumptions are less explicit. For example, few limestones will maintain conduits of 3 m radius without collapse causing an increase in friction factor, trapping of transported sediment, and a decrease in conduit radius. However, we can now imagine that such an event will be not unlike the sedimentation episodes pictured in Figure 7. Conduit breakdown could be readily incorporated in a rational way with a random algorithm generating collapse according to a radius dependent probability. The associated impact of local constrictions on the conduit has been described by Atkinson et al. (1983). Similar improvements incorporating more complex conduit systems and a diffuse flow component may simply require more processing providing a rational structure can be determined. More difficult are improvements in the overall geomorphic context of the model. In terms of form, much karst evolution is predicated by external controls of hydraulic gradients through incision and capture. Simple modelling of such attributes is usually futile as model outcome is often explicit in model design. More distressing difficulties concern our understanding of processes. How does conduit erosion vary with discharge and downstream distance? What are the relative rates of surface and underground channel erosion? How does solutional potential vary with changing climate? These factors could be readily incorporated into a more realistic model if our understanding were better. zyxwvutsrqpon zy A DEDUCTIVE MODEL OF KARST EVOLUTION 287 The present model describes hydrological behaviour not landscape evolution. Although there is a close association between karst hydrology and geomorphology, equilibrium requires considerable time. Instability of the kind portrayed by the present model suggests that landscapes will frequently be out of equilibrium with their hydrology. There has been no attempt to incorporate landscape evolution here, because, once again, we have no suitable data on the nature of landscape response to the various changes in hydrology suggested by the model. CONCLUSIONS The application of exceedance probability concepts to a physical model of the karst aquifer has provided a hybrid model applicable to the modelling of karst evolution. Karst may be conveniently defined on the basis of hydrological function rather than apparent form. The model first implies that karst underlain by a conduit aquifer will evolve from an initially fluvial condition, through fluviokarst to holokarst. Stream size has a profound influence on hydrological function and this depends upon the size of the sustainable surface catchment. The rate of evolution is primarily dependent on the saturation deficit of the inflowing water. Periodic increases in discharge or sediment deposition in conduits can cause a dramatic reversion in the hydrological function of the karst. This suggests that karst may be characterized by more or less frequent rejuvenation, rather than by a progressive, sequential, and irreversible evolution. Thus karst landscapes may frequently be out of equilibrium with hydrological function which rather complicates simple classification on the basis of form alone. The model has been developed in a simple manner to illustrate a few significant concepts. More realistic models will require information which may not be available. Furthermore, they will also be more complex, and significant conclusions may be more difficult to reach. ACKNOWLEDGEMENTS The Natural Sciences and Engineering Research Council of Canada provided both direct support for work summarized in this paper, and have also funded field studies which have given meaning to the modelling. T. C. Atkinson, A. N. Palmer, and D. C. Ford provided valuable comments on early drafts of this paper. zyxw zyx zyxwvutsr REFERENCES Atkinson, T. C. 1968. ‘The earliest stages of underground drainage in limestone-a speculative discussion’, Proceedings of the British Speleological Association, Papers presented to the annual conference, No. 6, 53-70. Atkinson, T. C., Harmon, R. S., Smart, P. L., and Waltham, A. C. 1978. ‘Palaeoclimatic and geomorphic implications of ThjU dates on speleothems from Britain’, Nature, 272, 24-28. Atkinson, T. C., Smart, P. L., and Wigley, T. M. L. 1983. ‘Climate and natural Radon levels in Castleguard Cave, Columbia Icefields, Alberta, Canada’, Arctic and Alpine Research, 15, 487-502. Coward, J. 1971. ‘Direct measure of erosion of a stream bed of a West Virginia Cave (Abstract)’, Caves and Karst, 13, 39. Cvijic, J. 1918. ‘L‘hydrographie souterraine et I’evolution morphologique du karst’, Revue de Ceographie Alpine, 6, 375-426. Davis, W. M. 1909. Geographical Essays, reprinted by Dover, New York 1954, 777 pp. Drake, J. J. 1980. ‘The effect of soil activity on the cbGmistry of carbonate groundwaters’, Water Resources Research, 16, 381-386. Drew, D. P. 1983. ‘Accelerated soil erosion in a karst area: the Burren, western Ireland’, Journal OfHydrology, 61, 113-126. Dury, G. H. 1964. ‘Principles of underfit streams’, US.Geological Survey Professional Paper, 452A. Ford, D. C. 1980. ‘Threshold and limit effects in karst geomorphology, in Coates, D. R. and Vitek, J. D. (Eds), Thresholdr in Geomorphology, George Allen and Unwin, London, 345-362. Gascoyne, M. 1981.‘Rates ofcave passage entrenchment and valley lowering determined from speleothem age measurements’, in Beck B. F. (Ed.), Proceedings of the VIII Int. Cong. Speleology, Bowling Green, Kentucky, 99-100. Gale, S. J. 1983. ‘The hydraulics of conduit flow in carbonate aquifers’, Journal of Hydrology, 70, 309-327. Grund, A. 1914. ‘Der Geograpische Zyklus im Karst’, Zeitschrift Ges Erdkunde Berlin, 52,621 4 4 0 . Translated in Sweeting, M. M. (Ed.) 1982. Karst Geomorphology, Benchmark Papers in Earth Sciences, 54-61. Jennings, J. N. 1985. Karst Geomorphology, Basil Blackwell, Oxford, U.K., 293 pp. Kirkby, M. J. 1978. ‘Soil development models as a component of slope models’, Earth Surface Processes, 2, 203-230. Lauritzen, S. E., Abbott, J., Arnesen, R., Crossley, G., Grepperud, D., Ive, A,, and Johnson, S. 1985. ‘Morphology and hydraulics of an active phreatic conduit’, Caue Science, 12, 139-146. Miller, T. 198I. Hydrochemistry. hydrology, and morphology of the Caves Branch Karst, Belize, Unpublished PhD. Thesis, McMaster University, 280 pp. 288 zyxwvutsrqp zyxwv C. C. SMART zyxwvut zyxwvutsrqp zy Palmer, A. N. 1981. ‘Hydrochemical factors in the origin of limestone caves’, in Beck, B. F. (Ed.), Proceedings ofthe VIII h i . Cong. Speleology, Bowling Green, Kentucky, 120-122. Palmer, A. N. 1984. ‘Geomorphic interpretation of karst features’, in La Fleur, R. G. (Ed.), Groundwater as a Geomorphic Agent, Allen and Unwin, Boston 1984, 173-209. Roglic, J. 1964. ‘Karst Valleys in the Dinaric Karst’, Erdkunde, 18, 113-1 16. Sawicki. L. R. 1909. ‘Ein beitrag zum geographischen Zyklus im Karst’, Geographische Zeitschrijt, 15, 185-204, 259-281. Schumm, S. A. 1969. ‘Geomorphic implications ofclimatic changes’, in Chorley, R. C. (Ed.), Water, Earth and Man, Methuen, London, 525-534. Smart, C. C. in press. ‘The exceedance probability distributions of steady conduit flow in karst aquifers’, .Hydrological Processes. Smith, D. I. 1975. ‘The residual hypothesis for the formation of Jamaican bauxite-a consideration of the rate of limestone erosion’, Journal of the Geological Society of Jamaica, 11, 3-12. Smith, D. I. and Atkinson, T. C. 1976. ‘Process, landforms and climate in limestone regions’, in Derbyshire, E. (Ed.), Geomorphology and Climate, Arnold, London, 367-409. Sweeting, M. M . 1973. Karst Landjorms, Columbia University Press, New York, 362 pp. Tratman, E. K. 1969. The Caves ofNorthwest Clare, Ireland, David and Charles, Newton Abbot, U.K., 256 pp. White, W. B. 1977. ‘The role of solution kinetics in the development of karst aquifers’, in Tolson, J. S. and Doyle, F. L. (Eds), Karst ffydrogeology.Memoir 12 Int. Assoc. Hydrogeol., 503 -5 17. White, W. B. 1984. ‘Rate processes: chemical kinetics and karst landform development’, In La Fleur, R. G . (Ed.), Groundwater as a Geomorphic Agent, Allen and Unwin, Boston, 227-248. White, W. B. and Longyear, J. 1962. ‘Some limitations on speleogenetic speculation imposed by the hydraulics of groundwater flow in limestone’, Nittany Grotto Newsletter (National Speleological Society), 10, 155-1 67. Williams, P. W. 1972. ‘Morphometric analysis of polygonal karst in New Guinea’, Bulletin ofthe Geological Society ofAmerica, 83, 761 796. -