ACTA CARSOLOGICA	34/1	1	9-23	LJUBLJANA 2005
COBISS: 1.01
THEORETICAL INVESTIGATION OF THE DURATION OF KARSTIC DENUDATION ON BARE, SLOPING LIMESTONE
SURFACE
HITROST KRAŠKE DENUDACIJE NA GOLEM APNENČASTEM POVRŠJU - TEORETIČNI PRISTOP
GÄBOR SZUNYOGH1
1 Dr. Szunyogh Gabor, Berzsenyi Daniel Foiskola, Fizika Tanszek, Szombathely, Karolyi Gäspär ter 4, H 9700, Email:szgabor@nadir.hmst.hu, Tel.: 00 36 94 504 352
Abstract	UDC: 551.311:551.44
Gabor Szunyogh: A Theoretical Approach to Establish the Duration of Denudation on Limestone Surface without Soil Cover
This paper deals with the question of how the duration of karstic denudation depends on the dip angle, the annual amount of precipitation, the rain intensity, and the prevailing wind direction and speed in case of an initially plane, sloping limestone surface without soil cover. The answer is given by the solution of a differential equation system describing the lowering speed of the rock surface. It turns out that the rate of the denudation does not increase in proportion to the intensity of precipitation and that it can never exceed a maximal value. Furthermore, long, soft rains result in higher annual denudation than short, abundant downpours. With increasing wind-speed the corrosion rate also increases, but above a certain wind speed the dissolution does not become faster. This paper presents numerical examples with diagrams about how these factors affect the expected duration of denudation. Key words: karstic denudation, lowering speed, rain intensity, slope angle, theoretical model, wind.
Izvleček	UDK: 551.311:551.44
Gabor Szunyogh: Hitrost kraške denudacije na golem apnenčastem površju - teoretični pristop
V članku obravnavam hitrost kraške denudacije na primeru nagnjene apnenčaste plošče, neposredno izpostavljene dežju. Zanima me, kako je hitrost denudacije odvisna od parametrov kot so naklon površine, količina in intenziteta padavin ter smer in hitrost vetra. Odgovor poiščem v rešitvi diferencialne enačbe, ki opisuje hitrost zniževanja površine. Izkaže se, da se hitrost denudacije ne veča sorazmerno z intenziteto padavin in da nikoli ne preseže določene mejne vrednosti. Dolgotrajno deževje majhne intenzitete povzroča višjo denudacijsko stopnjo, kot niz kratkotrajnih padavin. Hitrost denudacije narašča s hitrostjo vetra, vendar le do določene hitrosti. Ključne besede: kraška denudacija, hitrost zniževanja površja, intenziteta padavin, naklonski kot površja, teoretični modeli, veter.
INTRODUCTION
In case of field experiments a frequently asked question is how much time is necessary for the karstic denudation of a sloping limestone surface to reach the actual shape. The answer is usually uncertain since many factors (e.g. the aspect and the dip angle of the slope, the amount and the intensity of rainfall, the direction and the force of the wind, etc.) influence the denudation time and large differences are possible. It is well-known that several research teams perform regular „in situ" measurements in order to get the rate of denudation (High and Hanna 1970; White 2000). The different ways how the bare limestone surface is corroded and its morphology are the subjects of long and widespread research (Ford and Williams 1989; Veress 2000, 2004; Veress & Zentai et al. 2003; Telbisz 1999a, 1999b; Klimchouk & al. 2000, etc). Although most of the ideas take qualitatively into account that there is a strong correlation between the karst corrosion rate and the flow properties of the water resulting from the slope conditions of the area, quantitative equations have as yet not been formulated. Based on the limestone solution kinetics equations of Dreybrodt (1989), a mathematical model of this type of denudation has been developed by the author (Szunyogh 2000).
THE EQUATIONS OF LANDFORM EVOLUTION WITH CONSIDERATION TO WIND SPEED AND DAILY HOURS OF RAINFALL
The model is based upon the general equation system describing the denudation of bare (without soil cover) limestone surfaces. The aim of this system is to determine the function
z = f (x,y,t)	(1)
that gives the rock surface at time t, where z is the height above sea level, x and y are horizontal coordinates of the surface points (Fig. 1). The limestone is dissolved by carbonic acid in the rain water, consequently the shape of its surface keeps changing. The lowering speed of the rock surface can be interpreted by the following derivative:
dz
W = dt (2)
(The minus sign is given to express the fact that the surface is lowered due to karstic denudation, i.e. z(x,y,t) decreases as time passes.)
In order to determine the sinking speed or the shape of the surface, the following laws are built in the model:
1. The law of mass conservation is applied for the water. Let us consider theoretically a small volume of water that is in connection with the limestone surface (Fig. 2). The sides are vertical and parallel with the x and y axes, and the side lengths are dx and dy, respectively. The height is equal to the depth of water moving on the rock surface, and the top coincides with the water surface. According to the law of mass conservation the amount of water flowing into this small volume through the sides and the top during a time unit should be equal with the amount of water flowing out of it. It is mathematically expressed by the following formula:
d( mvx) + 3(mvy) = 1 q^^^n^	(3)
dx	dy	cos a
where is the density of water, m(x,y,t) is the depth of the thin water layer moving on the rock
surface, v (x,y,t) and v (x,y,t) are the x and y directional components of the velocity vector of the flowing
x	y
water, qran is the mass-flow density of rainfall referring to a unit area, n is the normal vector of the limestone surface (directed outward from the rock mass), is the slope angle of the limestone surface. (a has no sign: 0<a<90°). Physically 0° has no meaning, since the water does not flow upslope). The functions m, vx and vy are at the moment unknown, qrain can be expressed by multiplication of the velocity vector of raindrops (v ,) and the total mass of raindrops in a unit volume (p ,):
qrain Prain Vrain
(4)
The horizontal component of the velocity vector of raindrops is determined by the wind speed, vw, the vertical component by the falling speed of raindrops (relative to the air), v^^^^ (Fig. 3). öw indicates the angle of the wind direction relative to the x axis. (The wind direction is the direction from where the wind blows. 0<öw<360°.) Using these notations
Vrain = -Vw C0S ^w^ - Vw Sin ^wi - Vdrk
therefore
Qrain = -PrainVw C0S Swi - Pra^Vw Sin Sw j - PrainVdrk
(5)
(6)
Figure 1: The structure of the model.
In order to determine prai.n we can use the definition of the annual precipitation, Qa, that is the volume of rain water in a horizontal unit area over a year. Moreover, this water amount is directly proportional to the vertical component of the rain mass-flow density, so
pwater Qa = prainVdrtdNa ,	(7)
where td marks the daily sum of rainfall hours and Na is the number of days during a year. Expressing p , from Eq. (7) we get the following formula:
p _ P water Q
r rain
VdrtdNa
(8)
The values of n and cos a on the right of Eq. (3) can be obtained by partial differentiation of the function describing the shape of the rock surface:
n _
d ^ . d z . ,
1- — J + k
a x d y
	fd z ^	+	fd z ]	2 +1
	d x v J		y J	
(9)
Figure 2: Illustration of how the law of mass conservation is used.
and cos a =
	fd z 1	2 +	fd z 1	2 +1
	d x V y		d y V y	
(10)
Substituting the formulae (6), (8), (9) and (10) into (3) we get the following differential equation for the law of mass conservation
d( mvx) + d( mvy) = _Qa
dx
dy
VdrtdNa
d z	cd z . _
^ ^ vw cos 8w ^ ^ vw sin 8w + vc
a x	d y
. (11)
Eq. (11) gives a relationship between the velocity and depth of the water moving on the rock surface and the direction and force of the wind.
2. The Navier-Stokes equation. According to hydraulic calculations, water moves downwards on a rock surface in a thin layer as laminar flow with friction, therefore the velocity profile is parabolic as given by the Navier-Stokes equation (Fig. 4). The mean velocity along the profile is
„_P watergh2,,-
3n
-sin a,
(12)
where g is the gravitational acceleration, n is the dynamic viscosity factor of water and h is the thickness of the water layer. There exists a simple relationship between h and m:
h = mcos a.	(13)
Figure 3: The movement of raindrops.
Figure 4: The velocity profile of water in a thin water layer.
The water flows down along the slope lines of the limestone surface, consequently the direction of its velocity vector is determined by the gradient of the surface. Calculating the appropriate partial derivates of z(x,y) the following formulae are drawn:
-P watergh2sir 3n
(dz >	+	(dz 1
l^x J		dy V /
2 dz dx
(14)
and v„ = —

3ti
sin oc cos a
dx
\2

\2

(15)
(The minus sign refers to the fact that if the limestone surface is higher at larger x and y coordinates, i.e. the partial derivates of z are positive, then the water flows into the opposite direction towards smaller x and y coordinates.)
3. Relationship between the amount of dissolved CaCO3 and the lowering speed of the rock surface. Since precipitation water is aggressive due to the CO2 content of the air, it is able to dissolve the limestone (Fig 5). Therefore the rock surface lowers at a speed of
w = qr°ck
ProckCOS a
where prock is the density of rock, marks the mass-flow density of the limestone.
(16)
Figure 5: The relationship between the amount of dissolved limestone, the lowering of the surface and the concentration change of the water.
x
4. According to the law of reaction kinetics the larger is the difference between the equilibrium concentration, c and the actual calcium-carbonate concentration, c, of the water, the faster is the
' eq
solution process is. The mass-flow density of the limestone ^,moving" into the solution is proportional to this difference (Gabrovšek 2000):
k(Ceq -c), if c < Cs,
kn (Ceq - C)n,lf C > Cs,	(17)
qrock =
where k and k are solution velocity constants, n»4, c »0.9c . The values of c and k depend
n	s	eq	eq	^
on the temperature of the water and on the CO2-content of the air (Dreybrodt 1988). Since in the present model, c is always smaller than cs, hereafter only the upper part of formula (17) is taken into account.
5. The law of mass conservation applied for the CaCO. The law of mass conservation has its distinct validity also for CaCO3 (cf. Fig. 5). Applying it for the abovementioned small imaginary volume, the difference between the amount of (dissolved) CaCO3 input and the amount of the output through the sides of this volume is equal to the amount of CaCO3 in the solution dissolved from the limestone surface during a time unit. Mathematically:
d(mvxc) + d(mvyc)= qrock .	(18)
dx	dy cos a
The number of unknown variables (z, a, m, h, vx, vy, w, qrock and c) of equations (2), (10), (11), (13), (14), (15), (16) (17) and (18) is 9, therefore this system can - theoretically - be worked out. The independent variables of the equations are the x and y space coordinates and the t time. In order of solve this equation system one may start from a known initial limestone surface and so determine the future shape of the area at an arbitrary time, and also determine how much time is necessary for a given amount of denudation.
THE LAWS OF DENUDATION APPLIED FOR AN INITIALLY PLANE, SLOPING ROCK SURFACE
In the present paper we study the case of a limestone surface that is initially a plane surface sloping to the x direction with an angle of ao; therefore the function describing the shape of this surface may be considered to be independent from y (Fig. 6). Consequently, in the above equations the partial derivatives with respect to y become 0, and the partial derivate of the z function (describing the shape of the surface) with respect to x is proportional to the slope gradient:
^^ = - tg a.	(19)
dx
(The minus sign indicates the fact that with increasing x the limestone surface will decrease.)
Taking into consideration all of these, the equations (10)-(18) will result (after mathematical transformations) in the followings:
Gabor Szunyogh: A Theoretical Approach to Establish the Duration of Denudation on Limestone Surface without Soil Cover P water gm2
v=
3n
-sin a cos a,
dx (mvx )^:-QV (tg a cos ^wvw + vdr) , d( mvxc )= k
vdrtdNa
dx cos a
( ceq - c ) ,
w=
k ceq - c
Prock cos a
Initial and boundary conditions are given as follows: z (x,t )|t=0 = zo - tg ao • x .
vx (x,t) = 0 ,h (x,t) = 0.
^ x=0 , and ^ ' x=0
(20)
(21) (22)
(23)
(24)
(25)
The values of zo, ao, öw, vw, vdr, k, ceq, prock, n and g are given. The solution for w(x, t) and z(x, t) are to be determined.
Figure 6: Illustration of the intially plane slope model.
Substituting vx from (20) into (21):
d
dx
p watergm sin a cos3 a
3n
Qa
VdrtdNa
(tg a COS 8wVw + Vdr )
(26)
Let us localize (26) to the starting time of the denudation when t=0. According to the present model, the shape of the limestone surface is plane with a constant slope angle (ao), which is independent from x. Using this fact and taking into account the boundary condition (25), equation (26) can be integrated:
m (x,t)
1 I 3n sinao cos6wvw + cos aovdr Q^
cos P waterg
vdrtdNa
(27)
From this equation it can be stated that at the beginning of the denudation process, the depth of the flowing water layer on the rock surface is proportional to the cubic root of the distance from the top of the slope. Writing m into Eq. (20) the following relationship is drawn for the velocity of the flowing water:
vx (x,t)
p g
= wate^ sin ao cos ao
^ o o
3n
3 n sin aotg 8wvw + cos aovdr Qa
P waterg
vdrtdNa
(28)
Substituting m and vx into Eq. (22) and after lengthy but fundamental transformations the following differential equation can be drawn for the function c(x):
Q.
/ ■ s \ Qax dc , (sin ao cos8wvw + cosaovw)-^--_ kceq -c
k + -
- (sin ao cos 8wvw + cos ao vw)
.(29)
Eq. (29) can be separated according to the variables:
(sin ao cos 8wvw + cos aoVdr )-
Qa -dc
k^ - c
VdrtdNa
dx
k + —— (sin ao cos 8wvw + cos aovdr)
VdrtaNa
(30)
which can be integrated into the following formula:
k
c (x,t)
t_0	Q
k +-^ (sin ao cos 8wVw + cos aoVdr)
c+
VdrtdNa
_ kVdrstdN, +Q.(sinao cos8w -Vw +cosao -Vdr)
+ A . x	Qa (sinaocos 8w-Vw + cosao-Vdr )
(31)
where A is an integration constant. Since the second term of the right side of the equation converges to ^ as x^0, the physical interpretation of (31) requires that A=0. Therefore
t_0
sin ao cos ao
3
sin a cos a
oo
c ( x,t)
kvdrtdNa
k vdrtdNa + Qa (sin ao cos 8wvw + cos aovw ) eq '
(32)
The lowering speed of the limestone surface can be drawn if c is substituted into (32):
w (x,t)
k
Qa (sin ao cos 8wvw + cos ao vdr )
Prock cos ao k vdrtdNa + Qa (sin ao cos 8wvw + cos ao vdr )
. (33)
As the right side of this expression is independent from x, it means that at the beginning all points of the limestone surface have the same lowering speed, i.e. there is a parallel slope retreat. But this fact implies a very important consequence: if the rock surface was plane at the start it remains plane later as well. However, if the limestone surface preserves its plane shape with an angle ao then the formulae (27), (32) and (33) are valid for not only time t=0 but for any time later on, consequently w is constant in both space and time (Fig. 7). Similarly, according to (32) the concentration of water flowing on the surface is also independent from x. It means that there is a dynamic equilibrium between the 'diluting effect' of fresh rain water and the 'concentration-increasing' effect of the dissolution process.
The expression (33) does not give any direct information about the long-term mean lowering speed of the limestone surface because it is valid only for the duration of rainfall. The total annual lowering of the surface can be calculated by multiplying the surface lowering of a given time unit as expressed by (33) and the total annual duration of rainfalls:
w
= dztN = dttdNa
(34)
Figure 7: The shape of the intially plane limestone surface as a function of time.
t=0
c
t=0
a
therefore wa =
Qa (sin ao c0S SwVw + c0S ao Vdr ) tdNa
Prock coS aokvdrtdNa + Qa (Sin ao coS Sw ' Vw + coS ao ' Vw )
. (35)
For calculations it is recommended to build in some conversion constants by substituting parameters as used in everyday practice. Let Sh be the hour-second, Mmm the millimeter-meter conversion value. Adding it to the above formula we obtain:
= k ceq Qa (sin aocos SwVw + cos ao Vir ) tdShNaMmm
Prock coS ao k VdrtdShNa + Qa (Sin ao coS Sw Vw + coS ao Vdr ) '
(36)
The shape of the surface can be drawn by the integration of (2), which is easy in the present case since the function w at the left side is constant, as it was stated above:
z (x,t ) = Zo - wa • t - tg ao
(37)
The Eq. (36) gives the vertical velocity of surface lowering. However, the following question my arise: how thick is the layer that is peeled from the layer in a given time. Its annual rate can be calculated as follows:
ba = cos ao • wa
(38)
where ba is the annual thickness of the peeled rock layer.
Figure 8: The lowering speed of the surface.
a)	as a function of the annual precipitation (Td = 10 hours/day, vw = 0 m/s)
b)	as a function of the daily rainfall hours (Q = 2000 mm/year, v = 0 m/s)
k
c
x
Finally, the At time necessary for the Az denudation of a karstifying slope is given:
= Prock cos ao k vdrtdNa + Qa (sin ao cos ^w vw + cos ao vdr )
k ceq Qa (sin aocos 8wvw + cos aovdr ) tdShNaMmm
Az
(39)
QUALITATIVE AND QUANTITATIVE ANALYSIS OF THE DERIVED
RELATIONSHIPS
The parameter values in the aforementioned formulae are as follows: the dynamic viscosity of water: n =1,910-3 kg/ms; the density of limestone prock =2300 kg/m3; the gravitational acceleration g=9,81 m/s2. In case of open air carbon-dioxide content and T=5° C the velocity constant of the dissolution k=3,2.10-7 m/s, the equilibrium concentration of calcium-carbonate ceq=0,0546 kg/m3 (Dreybrodt & Eisenlohr 2000). The falling speed of raindrops vdr=3-8 m/s (Budo 1972). The numerical values of the conversion parameters: N =365 days/year, S=3600 s/hour, M =1000 mm/m.
A	^	h	mm
In Fig. 8a the lowering speed is presented as a function of annual precipitation at different slope angles (assuming t^^ =10 hours/day rainfall time). It can be observed that in case of ga=2-6000 mm/ year precipitation the lowering speed is 50-150 ^m/year. It is also clear that steeper slopes have a faster lowering speed. It is noteworthy that the increase of wa with growing Qa is not linear.
In Fig. 8b w^^ is represented as a function of daily rainfall hours at constant annual precipitation (Q=2000 mm/year). It is remarkable that a longer duration of rainfall results in faster denudation. It suggests that if the precipitation is mainly in forms of short, intense showers the denudation is less effective, whereas long, soft rains result in larger denudation. This can be explained by the fact that
Figure 9: The lowering speed of the surface as a function of a) the wind speed (dw= 0) and b) the wind direction (vw = 2 m/s). (Qa = 2000 mm/year, td = hours/day).
although intense rains dissolve more CaCO3 during the duration of rainfall, there is no corrosion during long inter-rain periods.
In Fig. 9a the denudation rate is shown as a function of wind-speed above the area (at constant values of Q and t ). It can be seen that wind increases the rate of denudation. This increase in cor-
^a	d'
rosion is, of course, larger as the slope becomes steeper because in this case the surface exposed to the wind is greater. But it is also observable that the increase of with growing wind-speed is not linear but converges asymptotically. Consequently, above a threshold-limit further increase in wind speed will practically not result in a higher lowering speed.
Fig. 9b illustrates the influence of the direction of the wind (at constant Qa, td and vw). It is remarkable how sensitive the denudation rate is to the aspect of the slope. Highest denudation values are expected at winds blowing against the slope. As the wind direction turns away, w decreases gradually and at öw>90° (i.e. the wind blows from „upstream", down the slope) the rate of denudation quickly drops and it may even become 0 depending on the slope angle. Fig. 9b clarifies how the extreme karren forms of stormy areas (Veress & al. 2003) can develop where the orientation of the karren is definitely in connection with the prevailing wind direction.
According to Fig. 10a the time necessary to cause corrosion decreases with growing annual precipitation, but it does not diminish below a certain limit. The time necessary for corrosion is shorter in case of long, soft rains (as evidenced by Fig. 10b) than in case of short, intense showers (assuming equal annual precipitation).
The diagrams demonstrate that the presented theoretical model is numerically in agreement with field experiment data, since the average denudation rate values are in the range of 50-500 ^m/a according to earlier field research (High & Hanna 1970, White 2000).
Figure 10: The time necessary for the corrosion of a 1 mm thick rock layer as a function of a) the annual precipitation (Td = 10 hours/day, vw = 0 m/s), b) the daily rainfall hours (Qa = 2000 mm/year, vw = 0 m/s).	w	a
CONCLUSIONS
According to the derived theoretical model, denudation of an initially plane limestone surface with a constant slope angle due to precipitation occurs in a way that the angle of the slope remains the same for longer periods, that is the limestone surface retreats parallel to itself. The rate of denudation is not proportional to rain intensity, i.e. if annual precipitation occurs mainly in the form of long, soft fall of rain the annual lowering speed of the limestone surface is higher than in the case of short and intense showers. The equations take also into account the direction and force of the wind during the rainfall. Increased wind-speed results in a higher corrosion rate on steeper rock slopes exposed to the wind. Corrosion rate, however, cannot surpass a certain limit in wind of any strength. The concentration of the thin water layer flowing down the rock surface is constant along the slope (i.e. it is in dynamic equilibrium). Inserting real values of precipitation in the equations, the denudation rate of the rock surface is in the range of 50-300 ^m/a.
REFERENCES
Budo, A. 1971: K^serleti fizika I. (Experimental Physics) - Tankönyvkiado, Budapest, 517.p. Dreybrodt, W. 1988: Processes in Karst Systems. Physics, Chemistry and Geology. - Springer,
Berlin, New York, Heidelberg, 294. p. Dreybrodt, W., Eisenlohr, L. 2000: Limestone dissolution rates in karst environments. - In: A. Klimchouk, D.C. Ford, A.N. Palmer & W. Dreybrodt (edit.), Speleogenesis: Evolution of Karst Aquifers. NSS, Huntsville. Ford, D.C., Williams, P. 1989: Karst geomorphology and hydrology. - Unvin Hyman, London, England, 601.p.
Gabrovšek, F. 2000: Evolution of early karst aquifers: from simple principles to complex models.
-	Založba ZRC, Ljubljana, 150.p.
High, C., Hanna, F. K. 1970: A method for the direct measurement of erosion on rock surfaces.
-	British Geomorphological Research Group Technical Bulletin, no. 5, 24 p. Klimchouk, A., D.C. Ford, A.N. Palmer & W. Dreybrodt (edit.), Speleogenesis^ Evolution of Karst
Aquifers. NSS, Huntsville. 527.p. Plummer, T., Wigley, T., Parkhurst, D. 1978: The dissolution of calcite on CO2-water systems at 5°
to 60°C and 0.0 to 1.0 atm CO2, - American journal of science, 278. p. 179-795 Szunyogh, G. 2000: The Theoretical-Physical Study of the Process of Karren Development - Karszt-
fejlodes IV. Szombathely. p. 125-150. Telbisz, T. 1999a: The qualitative and quantitative effect of precipitation on karst corrosion. In: Zambo
L., Telbisz, T. (eds.): Investigations of karst corrosion. - ELTE TTK, Budapest, p. 11-13. Telbisz, T. 1999b: Szäm^t6gepes szimulacio a felsz^nalaktanban. (Computer simulation in geomorphology) - Földr. Közl., 123. k. p. 151-162. Veress, M. 2000: The main types of karren development of limestone surface without soil covering
-	Karsztfejlodes IV. Szombathely. p. 7-30
Veress M.-T6th, G.-Zentai Z.-Czöpek I. 2003: Vitesse de recul d'un escarpement lapiaze (Ile Diego
de Almagro, Patagonie, Chili) - Karstologia No 41-1/2003, p. 23-26. Veress, M., T6th, G. 2004: Types of meandering karren - Z. Geomorph. N. F. 48. 1. p. 53-77 White, W. B. 2000: Dissolution of limestone from field observations. In: A. Klimchouk, D.C. Ford, A.N. Palmer & W. Dreybrodt (edit.), Speleogenesis: Evolution of Karst Aquifers. NSS, Huntsville.