Evolution of caves in porous limestone by mixing corrosion : A model approach

When water from the surface, e.g. from a lake, fl ows down through porous carbonate rocks, through a region with high hydraulic conductivity and encounters the water table of a phreatic aquifer, both waters mix by diffusion along their boundary. In a carbonate aquifer, where both surface and phreatic waters are saturated with respect to calcite, mixing corrosion causes renewed dissolution capacity Δceq of the carbonate rock in the diffusion-mixing zone, extending from the boundary separating the phreatic water from the surface water encountering it. A numerical model is presented from which the initial change of porosity in such a diffusion-mixing zone is obtained. The initial change of porosity can be calculated from the local distribution of the mixing ratio, and the second derivative of Δceq with respect to m. m(x,y)is the spatial distribution of the mixing ratio m= Vsurf /(Vsuf+ Vprh), where the V’s assign the corresponding volumes of surface and phreatic water. The second derivative has been calculated for three geochemical scenarios with differing CO2-concentrations of surface and phreatic water by the use of PHREEQC-2. The spatial distribution m(x,y) is obtained by using MODFLOW and MT3DMS in a modelling domain with constant hydraulic conductivity for various fl ow velocities of the phreatic aquifer. The time scale of cave evolution is estimated from the results. Passages of dimensions of about one metre in width and several 10 cm in height, extending in length along the boundary line, where surface and phreatic waters meet, can be created in time scales of 10 000 years. These caves are horizontal with blind ending passages and closely resemble the isolated caves observed in Central West Florida.


INTRODUCTION
Caves in young porous rocks are remarkably different from those in ancient fractured rocks.In the latter, speleogenesis is governed by fl ow along fractures, such as bedding planes and joints, and the matrix conductivity is negligible.Consequently, conduits develop along the fl ow directions in these fractures and their walls are solid rock without any porosity.The speleogenesis of these caves was the focus of research for two decades (DREYBRODT et al., 2005;DREYBRODT & GABROVSEK, 2003;PALMER, 1991;PALMER, 2007;FORD & WILLIAMS, 2007;KAUFMANN, 2005).
Caves in porous rocks, however, where the fl ow is not guided along fi ssures but where water is transported by Darcy fl ow through the matrix of the rock are different.They are known from young carbonate islands.Here mixing of saltwater with freshwater creates fl ank margin caves along the rim of the island.These isolated caves at the freshwater lens, originate as chambers of moderate size without entrances and exits, MYLROIE & CAREW, (1990, 2000, 2003), ROMANOV & DREYBRODT (2006), DREYBRODT & ROMANOV (2007).
In the porous rocks of Central West Florida, isolated caves are also common.FLOREA (2006), FLOREA &VACHER (2004), andFLOREA et al. (2007) have described caves of moderate passage length with blind ending passages.Most of these are aligned along NW-SE and NE-SW regional discontinuities as revealed from photo lineament features on the surface.The cave passages are not longer than several hundred metres and do not represent an integrated system of passages between inputs and outputs of the aquifer.The walls of these caves consist of porous rock with many small-scale solution features.The principal horizons of cave development appear between 3 m and 5 m; 12 m to 15 m; and 20 m to 22 m above present sea level.These horizons of cave development are related to prior phreatic water tables (FLOREA, 2006;FLOREA et al., 2007) of an unconfi ned aquifer.The caves originate from the following speleogenetic processes.When water from the surface above, e.g. a lake, reaches such a water table, it is already saturated with respect to calcite.It is mixed with the phreatic water, also in equilibrium with respect to calcite.When the CO 2 concentrations of both waters differ, mixing corrosion boosts the dissolution capacity (DREYBRODT, 1988).Due to the low hydraulic gradients and the slow fl ow velocity v y, between 10 -7 to 10 -6 m/s, (HANSHAW et al., 1965;DAVIS, 1996) the mixed solution again becomes saturated after several tens of centimetres.Therefore initial changes of porosity are restricted to a region close to the boundary line where surface water and phreatic water meet and voids directed vertical to the fl ow of the phreatic water originate.
Since both bodies of water exhibit Darcy fl ow, mixing of the phreatic water with the surface water can only be affected by molecular diffusion and transversal dispersion across the dividing streamline (DS).The diffusion-mixing zone between surface water and the phreatic water zone extends along the z-axis where the surface water enters.Mixing corrosion is active in this region and causes increasing porosity along the z-axis where the surface water enters, and down the y-direction along the dividing streamline.This way an isolated cave passage could originate, similar to those passages observed in porous rocks in Florida.The conduits are orientated normal to the fl ow direction of water, which has created them.This is similar to the concept of transverse speleogenesis as fi rst suggested by KLIMCHOUK (2007), although fl ow is through the matrix of the rock and horizontal, parallel to the water table in contrast to the fl ow ascending through fractures from below in Klimchouk's concept.Here, we will investigate how and over what time scale, caves can originate in such a hydrological setup.

THEORETICAL BACKGROUND
Figure 1 represents the geometry of the two domains of phreatic and surface water.It assumes a local input of water from a lake or a river, much larger than the input by meteoric precipitation, which therefore is neglected.Such geological settings are reported in lakes in the Central Lake District, Florida, USA, where huge amounts (up to 4 m/year) of water leak down vertically through the bottom of the lake to the water table of the Upper Floridan Aquifer (MOTZ, 1998;MOTZ et al. 2001).The domain of surface water in Fig. 1 is idealized and can have other shapes.In any case some kind of a new water table will originate above the former water table of the phreatic aquifer.In this mound of water, fl ow is phreatic and governed by Darcy's law.Because in Darcyfl ow, fl ow lines can never cross, they must become parallel at the boundary between the water of the underlying aquifer and the surface water or in other words along the dividing streamline.
This can be described in the modelling domain shown in Fig. 2. The left hand side input boundary idealizes the fi rst Fi gu re 1: a) Geometry of the hydrological concept: Water from a lake at the surface fl ows along some discontinuity with high hydraulic conductivity, down to the water table of an unconfi ned aquifer.The surface water encounters this water table along the z-axis and fl oats with some depth on the phreatic water along the streamline (DS), dividing the two domains of water.Both waters are saturated with respect to calcite, but have diff ering CO 2 concentrations.A diff usion-mixing zone is created along the dividing streamline (DS) of both water bodies.Therefore mixing corrosion becomes active creating porosity along the z-axis close to the dividing streamlines.Reactive transport in such a geochemical situation is described by the dispersion-advection-reaction equations of the species involved, which are Ca 2+ and CO 2 .These equations read (DREYBRODT & ROMANOV, 2007;PHILLIPS, 1991) (3 where Φ is the porosity of the rock, D is the constant of dispersion, c i the concentration of species i, q → the Darcy -fl ow rate and Q i the source term of dissolved material.In the limit of fast dissolution ensuring equilibrium of Ca 2+ with respect to calcite one gets (4) c eq mix (m,(x,y)), the equilibrium concentration of calcium with respect to calcite is a function of the mixing ratio m.
The local distribution of m can be represented by the distribution of a non-reactive tracer with a concentration between 1 and 0. Therefore it can be described by a dispersionadvection equation without a source term (6) M is the molar weight of CaCO 3 and r its density.m(x,y) is the local distribution of mixing in the aquifer.Therefore it is only necessary to solve the diffusion -advection equation ( 6) for mixing in the aquifer.From Eq.2 one has d 2 c eq mix / dm 2 = d 2 Dc eq / dm 2 .This purely geochemical entity is the second derivative of Δc eq (m) with respect to m, which can be obtained by use of PHREEQC as shown in Fig. 3b.To obtain dΦ/dt one has to know the mixing ratio m(x,y) and d 2 Δc eq /dm 2 as a function of the local coordinates.This method has been successfully applied to the evolution of fl ank -margin caves in the saltwater -freshwater diffusion-mixing zone of coastal carbonate aquifers (ROMANOV & DREYBRODT, 2006).
encounter of the two water bodies.The y-axis represents the initial part of the dividing streamline.Although the velocities could vary locally, in a fi rst approach we assume a constant fl ow velocity v y everywhere in the domain.The corresponding boundary condition is a constant flow rate dis tributed evenly to all cells of the input boundary.The upper and lower boundaries are impermeable to the transport of dissolved ions.This is true in nature only for the upper boundary representing the new water table.The lower boundary in nature is permeable, because the aquifer reaches much deeper.This, however, is only relevant when the width of the diffusion-mixing zone becomes close to the width of the modelling domain.Therefore care has to be taken that this does not happen in the model calculations.Choosing the width of 1m warrants this.We use cell sizes of 1cm by 1cm.A sharp boundary exists at x = 0.5m, y = 0 m, with calcium concentrations c eq phr for x < 0.5m and c eq sur for x > 0.5m.This sharp boundary cannot be maintained downstream or y > 0m, because there is transversal diffusion and dispersion causing mixing of the two fl uids when they fl ow down the aquifer.When a parcel of water moves downstream along the y-axis and arrives at the position v y •t = y after time t, the sharp boundary has been widened to a mixed fringe with width W x = 2(D•t) 1/2 in the x-direction, where D=D m +v y ·D m is, the transverse dispersion coeffi cient resulting from molecular diffusion with coeffi cient D m and transversal mechanical dispersion by fl ow with velocity v y through the pores with diameter d.
This mixing produces undersaturation with respect to calcite and consequently dissolution of carbonate rock.The concentration c mix of the fresh mixed solution depends on the degree of mixing m, which can be written as: with m= V sur / (V sur + V phr ), where V phr and V sur are the volumes of phreatic and surface water respectively.Due to the renewed aggressivity of the mixed water, its equilibrium concentration with respect to calcite increases by Δc eq .The equilibrium concentration of such mixed water can be written as c eq mix (m)=c mix (m)+Dc eq (m) (2) It can be calculated by use of PHREEQC (PARKHURST & APPELLO, 1999) as a function of m.Since the groundwater velocities are so slow, the mixed solution has suffi cient time to come to the new equilibrium, and we can assume that this is true everywhere in the diffusion-mixing zone.In the next step we have to fi nd d 2 Δc eq /dm for realistic geochemical scenarios.Ground waters in carbonate rocks are calcium-bicarbonate dominated.Therefore in a good approximation, for simplicity we assume both waters are pure calcium-bicarbonate waters in the system CaCO 3 -CO 2 -H 2 O, saturated with respect to calcite and in equilibrium with some low P CO2 -value characteristic for the phreatic water, and with a P CO2 -value signifi cantly higher for the water entering from the surface.With these conditions, both waters are completely defi ned and the chemistry of solutions mixed from these waters can be calculated by geochemical codes such as PHREEQC.We have investigated three different geochemical scenarios.In the fi rst scenario a) we assume P CO2 = 10 -3 atm for the phreatic water and 10 -2 atm for the water from the surface.This is a case dominated by an inorganic chemistry.Fig. 3a shows the results of PHREEQC at a temperature of 25 °C.It depicts Δc eq (m) as a function of m.Δc eq (m) can be fi tted by a quadratic expression.Therefore the second derivative becomes a constant d 2 Δc eq /dm 2 = -0.3mol/m 3 .
In the second scenario b) we assume a signifi cantly higher P CO2 -value in the surface water, which may have been caused by organic CO 2 -producing processes in the vadose zone.Therefore the P CO2 of the surface water is raised to P CO2 = 0.1 atm.Δc eq (m) also shown in Figure 3a, is no longer symmetric.In the region between m = 0.4 to m = 0.6 values are about -3 mol/m 3 one order of magnitude higher than in scenario a).It is just this region, where Ñ → (m) 2 exhibits its maximum.As a consequence porosity change is enhanced in scenario b) by one order of magnitude compared to scenario a).This will be discussed later in detail.
In scenario c) we assume that the water of the phreatic zone is in equilibrium with the P CO2 = 3 10 -4 atm of the atmosphere and the surface water now has a P CO2 of 0.05 atm.
The second derivative has values of -2.17 mol/m 3 at m = 0.5.Also this scenario exhibits enhanced porosity change compared to scenario a).Fig. 3b shows the second derivative for the three scenarios.
To calculate the distribution m(x,y) we use the modelling domain in Fig. 2. But, we now replace the two concentrations by their mixing ratios, m = 1 for the surface water and m = 0 for the phreatic water.The equation governing mixing m is Eq. 5, which we rewrite as To obtain the local distribution of m we create a numerical model using PMWIN (CHIANG & KINZELBACH, 2001).This software provides a front-end to MODFLOW (McDONALD & HARBAUGH, 1988) and MT3DMS (ZHENG & WANG, 1999).We use these models to calcula te the distribution of m in our modelling domain (see Fig. 2).

RESULTS
We use the parameters Φ 0 = 0.3, v y = 10 -6 m/s and d = 100µm.Furthermore we assume a constant hydraulic conductivity in the entire modelling domain.Fig. 4 depicts the distri bution of mixing m throughout the domain.As expected, the diffusion-mixing zone widens with distance from the entrance.By calculating numerically the partial derivatives ¶m/ ¶x and ¶m/ ¶y and applying Eq. 7 one obtains the initial change in porosity dΦ/dt /0 year -1 .The increasing porosity causes a local positive feedback.From Eq. 7 we have (8) and from this Using the positive feedback of increasing porosity expressed in Eqs. 8 and 9, one can estimate an upper limit of the time T for the change of porosity from its initial value Φ 0 = 0.3 to a value V 0 (T)=1 as T = -Φ 0 •τ•ln (Φ 0 ).Therefore with Φ 0 = 0.3 the maximum time T needed to excavate the region bordered by the isoline 1/τ = dΦ/dt /0 is about T = 0.36•τ.

Geochemical scenario A: The inorganic case
In the following we will fi rst discuss the results for the geochemical scenario a), where the second derivative is constant, because in this case the distribution of porosity change is proportional to (Ñ → (m(x,y)) 2 .Therefore all fi gures showing dΦ/dt /0 also refl ect the local distribution of (Ñ → (m(x,y)) 2 .Fig. 5a illustrates the spatial distribution of porosity change by a gray scale.Furthermore it depicts isolines of dΦ/dt /0 .
From Fig. 5 one visualizes that over a maximum time span of 120 000 years, a cave is created, extending to about 3 m in width and 0.2 in height.This is seen from the isoline with 3·10 -6 year -1 .To excavate a void of about 1m width and 0.1 in height (isoline 1·10 -5 year -1 ), only 38 000 years are needed.These numbers are the upper limits because porosity changes increase inside the region bordered by the isoline.Fig. 5b illustrates the distribution of initial change of porosity along transects in the x-direction with y = 0.1, 0.5, 1 and 2.5 m.Porosity change decreases with increasing distance from the input in the y-direction, such that evolution of porosity is restricted to a region extending a few metres away from the input in the fl ow direction.In the transverse direction to fl ow along x this region is limited to several tenths of a metre.
Flow velocities in the Floridan aquifer are as low as several 10 -7 m/s.Therefore in Fig. 6 we show the results for a fl ow velocity of 5·10 -7 m/s, everything else being unchanged in comparison to Fig. 5. Due to the low fl ow velocities, the region where the diffusion-mixing zone exhibits high gradi- ents and consequently high changes of porosity, becomes smaller as can be seen by the isoline of 3·10 -6 year -1 .The cavity created in 120 000 years extends only to 1m from the input in the y -direction and has a height of about 0.1m.The distribution of porosity change in the x -direction along the transect in Fig. 6b closely resembles that of Fig. 5b, but the porosity change is lower in the centre of the diffusion-mixing zone.

Geochemical scenarios b) and c):
The organic case The times of cave evolution in the inorganic scenario are too long to explain the existence of isolated caves in Florida.
These have originated during Quaternary sea level changes, where times of constant sea level elevations lasted no longer than 10 000 years.MYLROYE & MYLROYE (2007) have suggested that biogenic decay of organic matter caught on the top of the phreatic mound of the surface water increases its CO 2 concentration.We have taken this into account by increasing the P CO2 of the surface waters in the geochemical scenarios b) and c).As can be seen in Fig. 3b the second derivatives are higher by almost a factor of ten for m ≈ 0.5.Therefore one expects the evolution of caves in time scales of 10 000 years, supporting Mylroye's thesis.
We have modelled these organic cases.Transects of the evolution of porosity in comparison to the inorganic scenario are illustrated in Figs.5b and 6b.They show that porosity changes are indeed signifi cantly higher for both organic scenarios.Due to the asymmetry of the second derivatives the maximum of porosity evolution shifts upwards into the surface water region.These results indicate that caves with dimensions of about 1m width and 0.2 m in height can be created along the border where surface water with higher P CO2 encounters the water table of the phreatic base fl ow.This way we obtain blind ending passages of several 10 m to100 m length.These results are in concordance with the height to width ratios observed by FLOREA (2006) in the caves of West Central Florida.Most of them (47%) have ratios of about 0.3, similar to our values of about 0.1.The time needed to create such caves is in the order of 10 000 years.Note that our modelling assumes a constant sea level during the entire time of evolution to which the palaeo water table of the phreatic aquifer was tied.It should be noted here that at present in this work we have only modelled initial porosity changes and have just considered the feedback by Eqs.8-11 to estimate an upper limit of the time for cave evolution.Future modelling should take into account changes of hydraulic conductivity due to the change in porosity.The corresponding changes in the fl ow pattern by focusing of fl ow lines and the corresponding increase in mixing gradients might enhance porosity change in time.
Our modelling results show that due to the rise of P CO2levels by decomposition of organic matter in the water seeping down to the phreatic aquifer signifi cant voids and caves similar to those observed in Florida can originate in a time scale of 10 000 years as an upper limit.These are the maximum times of constant sea level stand in the Quaternary.
For the inorganic scenario, signifi cant changes of porosity from 0.3 to 0.5 arise within time spans of about 10 000 years in regions with dΦ/dt /0 = 10 -5 year -1 .In the organic scenarios, even if the times of constant boundary conditions are too short for the evolution of caves, changes in porosity by about 50 % need time spans of about 1000 years.

CONCLUSION
When two carbonate waters saturated with respect to calcite, but with different concentrations of calcium meet in an aquifer, they are mixed, and the spatial distribution of the mixing ratio m(x,y) describes the degree of mixing.Depending on m an amount Dc eq (m) is dissolved in the diffusion-mixing zone by mixing corrosion and the porosity of the rock increases.The change dΦ/dt in porosity depends critically on the spatial distribution of m(x,y) and on the chemistry of the mixed solution.When the two water bodies meet along a one dimensional border, isolated caves with blind ending passages can develop in the diffusion-mixing zone extending along such boundary lines.As long as these are located in the phreatic zone, fl ow through them is directed transverse to these passages and horizontally along the water table.
We have modelled cave evolution for a hydrogeological scenario where surface water with elevated P CO2 of 0.1 atm, originating from biogenic decay of organic matter, enters the water table of an unconfi ned aquifer below.The results of these modelling attempts closely resemble caves in West -Central Florida (FLOREA, 2006;FLOREA et al., 2007).
This concept of horizontal transverse speleogenesis could be used in the future to explain the evolution of isolated caves, in other karst areas, such as those described at Mt Scopus Formation (Israel) by FRUMKIN & FISCHHEN-DLER (2005).Our concept should not be mixed with the concept of hypogenic ascending speleogenesis by KLIM-CHOUK (2007).In this case, aggressive water ascends from adjacent formations below the cave level through fractures, and fl ow crosses the cave passages vertically.The modelling approach presented here is a fi rst step into gaining an understanding of speleogenesis of isolated caves in porous rock, which are so much different from the caves arising by dissolution of limestone with negligible porosity by aggressive solutions fl owing along fractures from an input to an output.
The small rectangle close to the origin of the coordinate system represents the region where the modelling domain is located.b) Close up view of a transect showing the rectangle in the x-y plane of Fig.1a.Due to the infi ltration, a mound of surface water builds up.Its upper boundary creates the new water table.Its lower border is the dividing streamline (DS) separating phreatic water from surface water.Black triangles indicate the fl ow directions.The rectangle MD represents the modelling domain, shown in Fig.2.a) b) From equations 4 and 5 Q ceq and from this the resulting change in porosity dΦ/dt in 1/s [volume (m 3 ) of dissolved rock per (m 3 /s) of the aquifer] can then be calculated by a simple relation (PHILLIPS, 1991; ROMANOV & DREY-BRODT, 2006).

Fi gu re 2 :
Modelling domain derived from Fig. 1.A constant fl ow with velocity v y is injected into the left hand side boundary.For 0 < x < 0.5m the infl owing water stems from the phreatic aquifer, corresponding to a mixing ratio m = 0, whereas the region from 0.5m < x < 1m carries water from the surface, with mixing ratio m = 1.The right hand side boundary is open for the outfl ow.The upper and lower boundaries are regarded as impermeable (see text).a) b)

Fi
gu re 3: a) Renewed solution capacity Δc eq for surface and phreatic water as a function of the mixing ratio m for three geochemical scenarios.a) (inorganic) P CO2 of phreatic water 0.001 atm, P CO2 of surface water water 0.01 atm, b) (decay of organic matter).P CO2 of phreatic water 0.001 atm, P CO2 of surface water water 0.1 atm, c) (decay of organic matter).P CO2 of phreatic water 0.0003 atm, P CO2 of surface water 0.1 atm.b) second derivatives d 2 Δc eq /dm 2 for the three scenarios a), b) and c).Fi gu re 4: Spatial distribution of the mixing ratio m for a homogeneous aquifer in the modelling domain of Fig. 2. The curves give isolines of m. v y = 10 -6 m/s.

Fi
gu re 5: a) Spatial distribution of initial porosity change dΦ/dt /0 .The isolines represent porosity change in year -1 for the geochemical scenario a (inorganic).v y = 10 -6 m/s; b) Vertical distribution of initial porosity change for the three geochemical scenarios a, b, and c along transects shown in panel a as vertical lines.Fi gu re 6: a) Spatial distribution of initial porosity change dΦ/dt /0 .The isolines represent porosity change in year -1 a) the geochemical scenario a (inorganic).v y =5·10 -7 m/s; b) Vertical distribution of initial porosity change for the three geochemical scenarios a, b, and c along transects shown in panel a as vertical lines.