Dreybrodt.indd 129 � AB STRA CT 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, ex- tending 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 corre- sponding 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. Keywords: speleogenesis, porous limestone, mixing corrosion, caves, Florida Evolution of caves in porous limestone by mixing corrosion: A model approach � Wolfgang Dreybrodt1, Douchko Romanov2 and Georg Kaufmann2 1Karst Processes Research Group, Institute of Experimental Physics, University of Bremen, 28359 Bremen, Germany; (dreybrodt@ifp.uni-bremen.de) 2Institute of Geological Sciences, FU Berlin, Malteserstr. 74-100, Building D, 12249 Berlin, Germany; (dromanov@zedat.fu-berlin.de, gkaufma@zedat.fu-berlin.de) doi: 104154/gc.2010.09 Geologia Croatica 63/2 129–135 6 Figs. Zagreb 2010 Geologia CroaticaGeologia Croatica 1. 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. Conse- quently, 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 Geologia Croatica 63/2Geologia Croatica 130 the rim of the island. These isolated caves at the freshwater lens, originate as chambers of moderate size without en- trances 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), and FLOREA et al. (2007) have described caves of moderate passage length with blind ending pas- sages. 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 inte- grated 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 ho- rizons of cave development are related to prior phreatic wa- ter tables (FLOREA, 2006; FLOREA et al., 2007) of an un- confi 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 CO2 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 vy, be- tween 10–7 to 10–6 m/s, (HANSHAW et al., 1965; DAVIS, 1996) the mixed solution again becomes saturated after sev- eral tens of centimetres. Therefore initial changes of poros- ity 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 af- fected by molecular diffusion and transversal dispersion across the dividing streamline (DS). The diffusion-mixing zone between surface water and the phreatic water zone ex- tends along the z-axis where the surface water enters. Mix- ing 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 horizon- tal, parallel to the water table in contrast to the fl ow ascend- ing 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. 2. THEORETICAL BACKGROUND Figure 1 represents the geometry of the two domains of phre- atic and surface water. It assumes a local input of water from a lake or a river, much larger than the input by meteoric pre- cipitation, which therefore is neglected. Such geological set- tings are reported in lakes in the Central Lake District, Flor- ida, 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 Darcy- fl 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 en- counters 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 CO2 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. 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 l- tration, 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 direc- tions. The rectangle MD represents the modelling domain, shown in Fig.2. a) b) Dreybrodt et al.: Evolution of caves in porous limestone by mixing corrosion: A model approach Geologia Croatica 131 Reactive transport in such a geochemical situation is de- scribed by the dispersion-advection-reaction equations of the species involved, which are Ca2+ and CO2. These equations read (DREYBRODT & ROMANOV, 2007; PHILLIPS, 1991) (3) where Φ is the porosity of the rock, D is the constant of dis- persion, ci the concentration of species i, q→ the Darcy – fl ow rate and Qi the source term of dissolved material. In the limit of fast dissolution ensuring equilibrium of Ca2+ with respect to calcite one gets (4) ceq 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 be- tween 1 and 0. Therefore it can be described by a dispersion- advection equation without a source term (5) From equations 4 and 5 Qceq and from this the resulting change in porosity dΦ/dt in 1/s [volume (m3) of dissolved rock per (m3/s) of the aquifer] can then be calculated by a simple relation (PHILLIPS, 1991; ROMANOV & DREY- BRODT, 2006). (6) M is the molar weight of CaCO3 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 ceq mix/ dm2 = d 2Dceq/ dm2. This purely geochemical entity is the second de- rivative of Δceq(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Δceq/dm2 as a func- tion of the local coordinates. This method has been success- fully applied to the evolution of fl ank – margin caves in the saltwater – freshwater diffusion-mixing zone of coastal car- bonate aquifers (ROMANOV & DREYBRODT, 2006). encounter of the two water bodies. The y-axis represents the initial part of the dividing streamline. Although the veloci- ties could vary locally, in a fi rst approach we assume a con- stant fl ow velocity vy everywhere in the domain. The corre- sponding boundary condition is a constant flow rate dis tributed evenly to all cells of the input boundary. The up- per 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 bound- ary 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 ceq phr for x < 0.5m and ceq 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 vy·t = y after time t, the sharp boundary has been widened to a mixed fringe with width Wx = 2(D·t)1/2 in the x-direction, where D=Dm+vy·Dm is, the trans- verse dispersion coeffi cient resulting from molecular diffusion with coeffi cient Dm and transversal mechanical dispersion by fl ow with velocity vy through the pores with diameter d. This mixing produces undersaturation with respect to calcite and consequently dissolution of carbonate rock. The concentration cmix of the fresh mixed solution depends on the degree of mixing m, which can be written as: cmix(m)=m·ceq sur+(1–m)·ceq phr (1) with m= Vsur/ (Vsur+ Vphr), where Vphr and Vsur 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 Δceq. The equilibrium con- centration of such mixed water can be written as ceq mix(m)=cmix(m)+Dceq(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. Fi gu re 2: Modelling domain derived from Fig. 1. A constant fl ow with velocity vy is injected into the left hand side boundary. For 0 < x < 0.5m the infl owing water stems from the phreatic aqui- fer, 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 out- fl ow. The upper and lower boundaries are regarded as imper- meable (see text). a) b) Geologia Croatica 63/2Geologia Croatica 132 In the next step we have to fi nd d 2Δceq/dm for realistic geochemical scenarios. Ground waters in carbonate rocks are calcium-bicarbonate dominated. Therefore in a good ap- proximation, for simplicity we assume both waters are pure calcium-bicarbonate waters in the system CaCO3-CO2 –H2O, saturated with respect to calcite and in equilibrium with some low PCO2 – value characteristic for the phreatic water, and with a PCO2 – value signifi cantly higher for the water enter- ing 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 PCO2 = 10–3 atm for the phreatic water and 10–2 atm for the water from the surface. This is a case dominated by an inor- ganic chemistry. Fig. 3a shows the results of PHREEQC at a temperature of 25 °C. It depicts Δceq(m) as a function of m. Δceq(m) can be fi tted by a quadratic expression. Therefore the second derivative becomes a constant d2Δceq/dm2 = –0.3 mol/m3. In the second scenario b) we assume a signifi cantly higher PCO2-value in the surface water, which may have been caused by organic CO2-producing processes in the vadose zone. Therefore the PCO2 of the surface water is raised to PCO2 = 0.1 atm. Δceq(m) also shown in Figure 3a, is no longer sym- metric. In the region between m = 0.4 to m = 0.6 values are about –3 mol/m3 one order of magnitude higher than in sce- nario a). It is just this region, where Ñ → (m)2 exhibits its max- imum. 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 PCO2 = 3 10–4 atm of the at- mosphere and the surface water now has a PCO2 of 0.05 atm. The second derivative has values of –2.17 mol/m3 at m = 0.5. Also this scenario exhibits enhanced porosity change com- pared to scenario a). Fig. 3b shows the second derivative for the three scenarios. To calculate the distribution m(x,y) we use the model- ling domain in Fig. 2. But, we now replace the two concen- trations 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 (7) To obtain the local distribution of m we create a numer- ical 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). 3. RESULTS We use the parameters Φ0 = 0.3, vy = 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 diffu- sion-mixing zone widens with distance from the entrance. Fi gu re 3: a) Renewed solution capacity Δceq for surface and phreatic wa- ter as a function of the mixing ratio m for three geochemical scenarios. a) (inorganic) PCO2 of phreatic water 0.001 atm, PCO2 of surface water water 0.01 atm, b) (decay of organic matter). PCO2 of phreatic water 0.001 atm, PCO2 of surface water water 0.1 atm, c) (decay of organic matter). PCO2 of phreatic water 0.0003 atm, PCO2 of surface water 0.1 atm. b) second de- rivatives d 2Δceq/dm2 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. vy = 10–6 m/s. a) b) Dreybrodt et al.: Evolution of caves in porous limestone by mixing corrosion: A model approach Geologia Croatica 133 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 (9) Using the positive feedback of increasing porosity ex- pressed 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 V0 (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·τ. 3.1. Geochemical scenario A: The inorganic case In the following we will fi rst discuss the results for the geo- chemical scenario a), where the second derivative is con- stant, because in this case the distribution of porosity change is proportional to (Ñ → (m(x,y))2. Therefore all fi gures show- ing 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 poros- ity changes increase inside the region bordered by the iso- line. 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 dis- tance 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 sev- eral 10–7 m/s. Therefore in Fig. 6 we show the results for a fl ow velocity of 5·10–7m/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- Fi gu re 5: a) Spatial distribution of initial porosity change dΦ/dt/0. The iso- lines represent porosity change in year–1 for the geochemical scenario a (inorganic). vy = 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 iso- lines represent porosity change in year –1 a) the geochemical scenario a (inorganic). vy =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. a) a) b) b) Geologia Croatica 63/2Geologia Croatica 134 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. 3.2. 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 CO2 concentration. We have taken this into account by increasing the PCO2 of the surface waters in the geochemical scenarios b) and c). As can be seen in Fig. 3b the second de- rivatives 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 sce- narios. Due to the asymmetry of the second derivatives the maximum of porosity evolution shifts upwards into the sur- face water region. These results indicate that caves with di- mensions of about 1m width and 0.2 m in height can be cre- ated along the border where surface water with higher PCO2 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. Fu- ture modelling should take into account changes of hydrau- lic conductivity due to the change in porosity. The corre- sponding 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 PCO2 – levels by decomposition of organic matter in the water seep- ing 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 maxi- mum times of constant sea level stand in the Quaternary. For the inorganic scenario, signifi cant changes of poros- ity 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 sce- narios, 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. 4. CONCLUSION When two carbonate waters saturated with respect to calcite, but with different concentrations of calcium meet in an aqui- fer, they are mixed, and the spatial distribution of the mixing ratio m(x,y) describes the degree of mixing. Depending on m an amount Dceq(m) is dissolved in the diffusion-mixing zone by mixing corrosion and the porosity of the rock in- creases. 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 pas- sages 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 PCO2 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 iso- lated 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 under- standing of speleogenesis of isolated caves in porous rock, which are so much different from the caves arising by dis- solution of limestone with negligible porosity by aggressive solutions fl owing along fractures from an input to an out- put. REFERENCES CHIANG, W.H. & KINZELBACH, W. (2001): 3D-Groundwater Mod- eling with PMWIN.− Springer-Verlag, Berlin, Heidelberg, New York, 397 p. DAVIS, H. (1996): Hydrogeologic investigation and simulation of ground- water fl ow in the Upper Floridan aquifer of north-central Florida and southwestern Georgia and delineation of contributing areas for se- lected City of Tallahassee, Florida, water-supply wells: U.S. Geolog- ical Survey Water-Resources Investigations Report, 95−4296. DREYBRODT, W. (1988): Processes in Karst Systems.− Physics, Chem- istry, and Geology, Springer, Berlin, 288 p. DREYBRODT, W. & GABROVŠEK, F. (2003): Basic processes and mechanisms governing the evolution of karst. / Speleogenesis and Dreybrodt et al.: Evolution of caves in porous limestone by mixing corrosion: A model approach Geologia Croatica 135 Evolution of Karst Aquifers 1 (1), re-published from: GABROVŠEK, F. (ed.): Evolution of karst: from prekarst to cessation. Postojna– Ljubljana: Zalozba ZRC, 115−154. DREYBRODT, W. & ROMANOV, D. (2007): Time scales in the evo- lution of solution porosity in porous coastal carbonate aquifers by mixing corrosion in the saltwater – freshwater transition zone.− Acta Carsologica, 36/1, 25−34. DREYBRODT, W., GABROVŠEK, F. & ROMANOV, D. (2005): Proc- esses of speleogenesis: A modeling approach.– Karst Research Institute, ZRC SAZU, Postojna, Slovenia. FLOREA, L.J. (2006): Architecture of air-fi lled caves within the karst of the Brooksville Ridge, west-central Florida.− Journal of Cave and Karst Studies, 68/2, 64–75. FLOREA, L.J., VACHER, H.L., DONAHUE, B. & NAAR, D. (2007): Quaternary cave levels in peninsular Florida.− Quaternary Science Reviews, 26, 1344−1361. FLOREA, L.J. & VACHER, L.H. (2004): Cave levels, marine terraces, paleoshorelines, and the water table in peninsular Florida.− Ar- chives of climate change in karst. Karst Waters Institute Special Publication, 10. FORD, D.C. & WILLIAMS, P.W. (2007): Karst Hydrogeology and Ge- omorphology.− John Wiley & Sons Ltd., 562 p. FRUMKIN, A. & FISCHHENDLER, I. (2005): Morphometry and dis- tribution of isolated caves as a guide for phreatic and confi ned pale- ohydrological conditions.− Geomorphology, 67, 457−471. HANSHAW, B.B., BACK, B. & RUBIN, M. (1965): Radiocarbon De- terminations for Estimating Groundwater Flow Velocities in Central Florida.−. Science, 23, 494–495, doi: 10.1126/sci ence.148.3669. 494. KAUFMANN, G. (2005): Structure and evolution of karst aquifers: a fi nite-element numerical modeling approach.− In: DREYBRODT, W., GABROVŠEK, F. & ROMANOV, D. (eds.): Processes of Spe- leogenesis: A modeling approach. Karst Research Institute, ZRC SAZU, Postojna, Slovenia, 323−371. KLIMCHOUK, A. (2007): Hypogene Speleogenesis: Hydrological and Morphogenetic Perspective, National Cave and Karst Research In- stitute. Special Paper No.1. McDONALD, M.C.& HARBAUGH, A.W. (1988): MODFLOW, A modular three – dimensional fi nite difference ground – water fl ow model.− U. S. Geological Survey, Open – fi le report, 83−875. MOTZ, L.H. (1998): Vertical leakage and vertically averaged vertical conductance for karst lakes in Florida.− Water Resources Research, 34, 159−167. MOTZ, L.H., SOUSA, G.D.& ANNABLE, M.D. (2001): Water budget and vertical conductance for Lowry (Sand Hill) Lake in north cen- tral Florida, USA.− Journal of Hydrology, 250, 134−148. MYLROIE, J.E. & CAREW, J.L. (1990): The fl ank margin model for dissolution cave development in carbonate platforms. Earth Surface Processes and Landforms, 15, 413–424. MYLROIE, J.E. & CAREW, J.L. (2000): Speleogenesis in coastal and oceanic settings.− In: KLIMCHOUK, A., FORD, D.C., PALMER, A. & DREYBRODT, W. (eds.): Speleogenesis: Evolution of Karst Aquifers. National Speleological Society, Huntsville, 226–233. MYLROIE, J.E. & CAREW, J.L. (2003): Karst development on carbon- ate islands. Speleogenesis and Evolution of Karst Aquife & HAR- RIS, P.A. (eds.), 1995. Unconformities in Carbonate Strata – Their Recognition and the Signifi cance of Associated Porosity, AAPG Memoir, 63, 55–76. MYLROIE, J.E. & MYLROIE, J.R. (2007): Development of the Car- bonate Karst Model.– Journal of Cave and Karst Studies, 69, 59−75. PALMER, A.N. (1991): Origin and morphology of limestone caves.− Geological Society of America BULLETIN, 103/1, 1−21. PALMER, A.N. (2007): Cave Geology.− Cave Books, 454 p. PARKHURST, D.L. & APPELO, C.A.J. (1999): User's guide to PHREE- QC (version 2). A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations: U.S. Geological Survey Water-Resources Investigations Report 99−4259, 312 p. PHILLIPS, O.M. (1991): Flow and Reactions in Permeable Rocks.− Cambridge University Press, New York, 285 p. ROMANOV, D. & DREYBRODT, W. (2006): Evolution of porosity in the saltwater-freshwater diffusion-mixing zone of coastal carbonate aquifers: An alternative modeling approach.− Journal of Hydrology, 329, 661−673. WERTH, C.J., CIRPKA, O.A. & GRATHWOHL, P. (2005): Enhanced mixing and reaction through fl ow focusing in heterogeneous po- rous media Water Resour. Res., 42, 12, W12414 10.1029/2005WR 004511 ZHENG, C. & WANG, P.P. (1999): MT3DMS, A modular three-dimen- sional multi-species transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and user’s guide, U.S. Army Engineer Re- search and Development Center Contract Report SERDP-99-1, Vicksburg, MS, 202 p. Manuscript received November 24, 2009 Revised manuscript accepted March 19, 2010 Available online May 31, 2010 Geologia Croatica 63/2Geologia Croatica 136