LETTER: CALCITE DISSOLUTION UNDER TURBULENT FLOW CONDITIONS: A REMAINING CONUNDRUM PROBLEM NASTANKA FASET IN RAZTAPLJANJA KALCITA V TURBULENTNEM TOKU Matthew D. COVINGTON1 UDK 551.435.8:549.742.111 INTRODUCTION A significant body of experimental and theoretical work has examined the dissolution rates of calcite, and other carbonate minerals, under varying chemical and hydro-dynamic conditions (see Morse & Arvidson (2002) for a comprehensive review). The relationships derived from this work have been applied extensively to the development of mechanistic models of speleogenesis (e.g. Drey-brodt 1996; Dreybrodt et al. 2005; Birk et al. 2005; Rehrl et al. 2008; Kaufmann 2009; Szymczak & Ladd 2011). However, the primary focus of these models has been on the early stages of cave formation, with less attention toward the later stages of cave evolution and turbulent flow conditions. Recent efforts have begun to develop mechanistic models for processes governing later stages of cave evolution, considering factors such as turbulent flow structures (Hammer et al. 2011) and open channels (Perne 2012). However, such studies remain limited, in part due to significant quantitative uncertainties in a va- riety of processes that become important beyond the incipient speleogenesis stage (Covington et al. 2013). While speleogenetic models have not typically been run much beyond the transition from laminar to turbulent flow conditions, experimental and theoretical studies of carbonate dissolution have constrained dissolution rates and rate-controlling mechanisms under turbulent flow (e.g. Rickard & Sjöberg 1983; Dreybrodt & Buhmann 1991; Liu & Dreybrodt 1997). However, a direct application of these results and comparison to field observations leads to an apparent conundrum. According to the theory, solutional forms such as scallops and flutes should not exist in limestone; however, they clearly do exist at a wide variety of sites and scales. This conundrum suggests that there may be problems with the theory, problems with our understanding of scallop formation, or both. THE CONUNDRUM Dissolution rates under turbulent flow conditions are typically calculated by dividing the fluid into two regions, a turbulent core that is well-mixed, and a diffusion boundary layer (DBL) that lies between the turbulent core and the dissolving wall (Dreybrodt & Buhmann 1991). Within the DBL, flow is presumed to be laminar and species cross the layer via the process of molecular diffusion. There are two end-member regimes in which overall dissolution rates are limited by either diffusion rates (the transport-limited regime), 1 University of Arkansas Department of Geosciences Fayetteville, Arkansas, USA, e-mail: mcoving@uark.edu Received/Prejeto: 26.11.2013 Fü (1) or surface reaction rates (the reaction-limited regime), F, = /c,(l-Cb/Ce,), (2) where FD and Fs are the transport-, and reaction-limited fluxes, respectively, D is the diffusion constant, Ceq is the equilibrium concentration of calcium, Cb is the concentration of calcium in the bulk solution, e is the DBL thickness, and fcj is the kinetic rate constant for cal-cite dissolution in the linear kinetic regime. The mixed kinetic regime occurs when the diffusion-limited and surface-reaction-limited rates are similar (Rickard & Sjöberg 1983). In this case, both processes are important. For simplicity, I only consider linear dissolution kinetics 10^ in 9 'ai C s S S CO 1 DC 10 10" Q 10" ---SurtaceHeaction Hate 10"^ 10-" 10-' 10"^ 10" Hydraulic Diameter (m) 10' Fig. 1: Diffusional (solid) and surface reaction (dashed) dissolution rates as a function of hydraulic diameter, DH, depicted for different hydraulic head gradients, Vh. Sharp jumps in diffusion rates occur at the onset of turbulent flow. The laminar/turbulent transition is indicated with dotted lines. For the parameter space shown, the theory suggests that diffusional rates only control dissolution for low-gradient conduits just below the laminar/turbulent transition. Elsewhere, diffusional rates are much higher than surface reaction rates. Figure reproduced with permission from Covington et al. (2012). here, as is common below about 80% saturation, though our conclusions are not sensitive to this choice. In the case of a sufficiently thick DBL, or a sufficiently thin turbulent core, conversion of CO2 can also be rate-limiting (Dreybrodt & Buhmann 1991). Though non-linear kinetics and CO2 conversion limitation are not explicitly considered here, if anything, they would exacerbate the presented conundrum by providing further means to reduce surface reaction rates. In order to determine the typical conditions under which rates were limited by either reaction or transport, Covington et al. (2012) calculated surface reaction and diffusion rates for a wide variety of head gradients and hydraulic diameters using the Darcy-Weisbach equation and Colebrook-White relation. For these calculations, the fiducial parameter values given in Dreybrodt et al. (2005) were used, with k1 = 4 x 10-11 mol cm-2 s-1, D = 10-5 cm2 s-1, and Ceq = 2 x 10-6 mol cm-3. For ease of comparison, Covington et al. (2012) converted the "surface" rate law to the form Fs = as(Ceq-C), (3) where as = 2 x 10-5 cm s-1 was determined by dividing k1 by Ceq. Here I have slightly modified the notation of Covington et al. (2012), replacing a with as to clarify that the constant applies to the chemical reaction rate at the surface. In many previous manuscripts (e.g. Buhmann & Dreybrodt 1985) a without a subscript represents the combined effects of diffusion and reaction. Fig. 2: Scallops formed on a limestone surface within a cave, indicating contrasts in dissolution rates as a result of turbulent flow structures. Major divisions on the ruler are in cm (Photo: Matija perne). Covington et al. (2012) compared the rates predicted by Equations 1 and 3 by nondimensionalizing the two equations by dividing by the term as (Ceq - C), such that the dimensionless dissolution rates are independent of chemistry and become F? = 1, (4) and P*. D aß (5) where the stars denote dimensionless rates. Then e was calculated for a wide variety of head gradients and hydraulic diameters to examine the relative magnitude of the reaction-and transport-limited equations. The DBL thickness was calculated using equations 2.13 and 2.14 from Dreybrodt et al. (2005), which are equivalent to Equations 9 and 10 of this manuscript. For these calculations a fractional roughness of 0.05 was assumed. The results of this calculation (Fig. 1) show that, according to the equations used above, the surface reaction rate is almost always limiting, and no cases in the turbulent flow regime show transport-limited rates. In fact, diffusion rates under turbulent flow are typically several orders of magnitude larger than the surface reaction rate. This conclusion leads to the conundrum. The model suggests that dissolution rates should almost always be controlled by the surface reaction rate in the turbulent flow regime, and, therefore, that dissolution rates should be independent of any spatial variations in DBL thickness. On the con- trary, caves and channels formed in limestone frequently contain scallops (Fig. 2) and flutes that apparently form as the result of systematic contrasts in DBL thickness as a result of turbulent flow structures (Curl 1966; Good-child & Ford 1971; Blumberg & Curl 1974; Curl 1974). The remainder of this note will explore this conundrum in more detail, and discuss possible resolutions. SURFACE RATES FROM THE PWP EQUATION The above analysis makes use of fiducial kinetic constants from Dreybrodt et al. (2005). However, these constants are motivated by typical values from experiments, where the effects of transport and reaction are both present. Therefore, a more accurate approach is to use the full Plummer-Wigley-Parkhurst (PWP) equation (Eqn. 6, Plummer et al. 1978) to calculate true theoretical surface rates. If the PWP rates were high enough, then they could explain away the conundrum. For comparison, I calculate PWP rates that result from a temperature and partial pressure of CO2 that approximately reproduce the fiducial value of Ceq = 2 X 10-6 mol cm-3 used above. This equilibrium concentration is roughly obtained (for an open system) with t = 10 °C and Pco2 = 0.01 atm, values which are also quite typical in the cave environment. The PWP rates are depicted in Fig. 3, alongside the simplified linear relationship with two different values of as. The form of the PWP equation over the relevant range is approximately linear for this example, except for the highly undersaturated end, suggesting that use of a linear function is reasonable in this case, except at nearly zero dissolved Ca. However, comparing the fiducial value of a^ used to create Fig. 1 above, with the best-fit value for the linear region of the PWP equation (Fig. 3), suggests that the value of as from Dreybrodt et al. (2005) is too small by about a factor of 3 to 4 in order to approximate the PWP rate. Consequently, the surface reaction line in Fig. 1 should be shifted up by a similar factor of 3 to 4. While this does imply that the surface rates are higher than would be predicted by the equation and parameter values given in Tab. 2.1 of Dreybrodt et al. (2005), the more accurate surface rates obtained from PWP are still much lower than the diffusive rates are predicted to be for the majority of the turbulent region of the parameter Fig. 3: The dissolution rate as a function of Calcium concentration for pCO2 = 0.01 atm and T = 10 °C. This example is constructed to have a similar value of Ceq and temperature to the fiducial case from Dreybrodt et al. (2005). space (Fig. 1). For the convenience of future studies, a fitting function for PWP rates as a function of temperature and PCO2 is provided in the Appendix. This relation is not used explicitly here, but provides a quick means of estimating PWP rates in a given setting. COMPARING SURFACE-REACTION AND DIFFUSION RATES The equations given above for diffusion-limited (Equation 1) or surface-limited rates (Equation 3) are only val- id under conditions where dissolution is truly limited by the respective process. In steady state, diffusion and sur- face reaction rates must be equal, and the rate of diffusion is given by D(Cs-Q) Pß--e-' (6) where Cs is the concentration at the surface, and Cb is the concentration in the bulk flow. For the diffusion-limited case Cs ^ Ceq, and we recover Equation 1 above. Similarly, surface reaction rates are dependent on the concentration at the surface and are given by Fs = as(Ceq-Cs). (7) Again, if dissolution is reaction-limited Cs ^ Cb and we recover Equation 3. Setting diffusion and reaction rates ecrit = D/as. (9) When e >> ecrit then rates are diffusion-limited, and when e << ecrit rates are reaction-limited. If e ~ ecrit then dissolution will occur via mixed kinetics. At e = ecrit diffusion surpresses dissolution rates to half of what they would be from surface reaction alone. The above analysis relies on a surface reaction rate of the linear form given by Equation 7. However, it can be generalized using the full PWP dissolution rate. In the limit of a thick DBL the diffusion-limited equation can be applied, while in the limit of thin DBL the PWP surface rate can be applied. An estimate of ecrit can be 10' a> [E (D I 10° 3 CO 10" ffrit approximated by extrapolation Surface-limited rale ju uy B) \ Diffusion-limited Rate : \Mixed Kinetic Rate Diffusion Surface Controlled i Controlled Tf^ 10" 10^ Fig. 4: The critical DBL thickness, below which dissolution rates are surface reaction-limited, can be estimated by extrapolating the diffusion-limited rate until it intersects the surface-limited rate determined by the pwp Equation. This leads to the relation in Equation 10. Though this approximation scheme can be employed more generally, the lines plotted in this figure represent exactly the case where the rates can be approximated using aD and as. equal, solving for Cs, and plugging this back into Equation 7, one can see that the dissolution rate accounting for both reaction and diffusion is given by (8) where ad = D/e. The lines depicting the diffusion-limited rates in Fig. 1 are equivalent to ad/as, so the analysis and conclusions of Covington et al. (2012) can also be cast in terms of the relative magnitudes of ad and as. Therefore, if dissolution rates are well-represented by the linear relation in Equation 7, then the relative importance of diffusion and reaction in determining the overall rates can be quantified by comparing ad and as. A critical diffusion boundary layer thickness, at which diffusion and reaction are equally important can be calculated by setting ad = as, leading to 0,25 0.20 0,15 S Ü.1U 0,05 — P002-0.0003. T=5.0 -- P002-0.0003.T=15-0 ... P002=0.0003. T=a5-0 K — P002-0,003, T=5,Ü I ^^ -- P002=0,003,T=15.0 PC02=0,003, T=25.0 — PC02=0,03, T=5.0 -- PC02=0,03, T=15,0 . . PC02=0,03, T=25,0 - PC02=0,1, T=5,0 -- PC02=0.1,T=15.0 -.-. PC02^.1, T=25.0 O.ÜÜÜ 0.002 0.004 0.005 (molL-^) 0.008 0.010 Fig. 5: The critical DBL thickness, below which dissolution rates are surface reaction-limited, shown for a wide range of temperature, Poo., values, and dissolved loads. Lines terminate on the right-hand side at saturation. obtained by extrapolating the diffusion-limited equation until it intercepts the PWP surface rate (Fig. 4). Quantitatively, this results in the relation ^crit- D(C„ - C) p (10) When the DBL thickness is near this critical value, dissolution will occur via mixed kinetics, but when the DBL is much larger or smaller than ecrit then dissolution rates will be limited by either diffusion or surface reaction, respectively. Using this relation, and the full PWP equation, one can calculate critical DBL thicknesses for a range of temperatures, PCO2 values, and calcium concentrations. This calculation (Fig. 5) shows that the critical DBL thickness is on the order of magnitude of 1 mm for essentially the entire range of temperature, PCO2, and dissolved load found in natural karst systems, except in highly undersaturated conditions where C ^ 0 and ecrit << 1 mm. tte approximation ecrit ~ 1 mm can also be reproduced using a rough estimate of PWP rates. In most natural settings, the forward reaction rate is dominated by reaction III (k3 in the PWP equation) (Dreybrodt 1988). ttis rate is simply a function of temperature and is on the order of 10 10 mol cm (Dreybrodt 1988, p. 127). Critical DBL thickness can then be estimated as £crit~ DC«, (10-5cm2s-')(2xl0-®molcm-') = 0.2 cm (11) lO-x'molcm-^ which is quite similar to the result of 1 mm obtained from the larger parameter study (Fig. 5). Within speleogenesis models DBL thickness is usually estimated from a relationship for pipe flow that employs the Sherwood Number (Sh), e = ÖH/Sh, (12) where Sh is given by an empirical relationship such as Sh- y/8)(Re - 1000)Sc - 1) (Dreybrodt et al. 2005). Here, Re is the Reynolds number, Sc is the Schmidt number, and f is the Cole-brook-White friction factor. Using these relations, one can calculate typical values of DBL thickness, (e), for a selection of head gradients and hydraulic diameters (Fig. 6). Again, for this calculation a relative pipe roughness of 0.05 is assumed, though the result is not particu- Fig. 6: DBL thickness (e) under turbulent flow conditions as estimated from Sherwood number for a wide range of hydraulic parameters. Each line represents a choice of head gradient. The lines terminate on the left at Re = 4000, where the Colebrook-White formulation breaks down. typical values of DBL thickness as estimated from Sherwood Number are much less than the 1 mm magnitude of the critical value determined from the pwp equations. larly sensitive to this choice. Typical values of the DBL thickness under turbulent conditions, as calculated from Sherwood Number, are much less than the order of magnitude estimate of a critical DBL thickness of 1 mm. This is true for all turbulent flow cases except those at very low head gradients (i.e. 10-5). CONCLUDING THOUGHTS AND POTENTIAL RESOLUTIONS A more careful analysis reinforces the conundrum. The theory suggests that limestone dissolution rates under turbulent flow conditions should be typically limited by surface reaction rates. However, field observations of scallops and other solutional forms clearly suggest otherwise. A recent model developed using PWP dissolution equations and computational fluid dynamics also found that limestone dissolution flutes were unstable under the calculated dissolution rates, though flutes in gypsum, which has much higher surface rates, proved to be somewhat more stable (Hammer et al. 2011). Superficially, this result is in agreement with the analysis presented here, and may provide further evidence of difficulties with the current theory. The central purpose of this letter is to clearly state the problem, rather than to solve it, in hopes of stimulating some discussion on the subject. However, I here discuss a few potential resolutions. Perhaps the most suspicious component of the turbulent dissolution model is the semi-empirical equa- tion used to calculate Sherwood Number (Eqn. 13), and consequently DBL thickness (Eqn. 12). This equation is based on experiments in smooth pipes, and therefore ought to be used with caution when calculating DBL thickness for the rough and irregular surfaces found on natural bedrock channel walls. However, the experiments of Blumberg & Curl (1974) provide a means of checking this as a potential resolution. Mass transfer data from their experiments with flutes in gypsum suggest a typical DBL thickness of e = 0.0089 L Sc-1/3, where L is the flute length and Sc is the Schmidt Number, with Sc a 1000 for water at the relevant temperatures. This results in a DBL thickness that is roughly s 0.1% of the flute length. In combination with the estimate of critical DBL thickness of 1 mm for limestone, this result suggests that such so-lutional forms should develop, but only on length scales of a meter or greater. This again is in conflict with field observations, where scallops and flutes often form on length scales of centimeters (Fig. 2). Another possibility is that the rates estimated by the PWP equation are significantly lower than rates on natural surfaces. In fact, several studies of dissolution rates measured on rotating disks produce rates that are approximately twice those given by the PWP equation (Dreybrodt & Buhmann 1991; Liu & Dreybrodt 1997). These studies have suggested that increased surface roughness, and therefore surface area, on the disks may result in the higher apparent dissolution rates. Similarly, roughness on natural limestone surfaces on a smaller scale than the DBL thickness could result in an apparent increase in surface reaction rates that would reduce the critical DBL thickness. However, the observed discrepancy of a factor of two is not sufficient to explain away the conundrum by itself. If natural surfaces are still significantly rougher at scales below the DBL thickness than the disks used in experiments, this could potentially resolve the problem. Other factors, such as microbial films, could also influence surface reaction rates, and perhaps in some cases increase them. One might ask whether a very small contrast in dissolution rates would be sufficient to form flutes and scallops. However, flute experiments show that dissolution rates vary by roughly a factor of two over the length of a flute (Blumberg & Curl 1974), a number that is thought to be relatively constant across the parameter space (Curl 1966). A constrast in dissolution rates of a factor of two requires that the thickest portion of the DBL within the scallop is greater than or equal to the critical thickness, as this is the DBL thickness when surface rates are suppressed by roughly a factor of two. A thinner DBL would not allow sufficient contrast in rates. Forms somewhat similar to scallops (mechanical erosion flutes) also are found in bedrock channels in relatively insoluble rocks, though the variety and prevalence of such forms is greatest in highly soluble rocks (Richardson & Carling 2005). One possible explanation is that so-called solutional forms, such as scallops and flutes, actually form by a combination of solutional and mechanical processes. For example, chemical processes could loosen individual grains that are later plucked from the surface by mechanical processes. Finally, it could be that scallops form only in highly aggressive waters, where the critical DBL thickness is sufficiently small (Fig. 5). However, the author's experience would suggest that scallops are also present in locations without such highly aggressive water. Little systematic attempt has been made to study the locations in which scallops form, and whether forms differ according to hydrological or lithological settings. Such studies might also provide clues as to the correct resolution of the current conundrum. ACKNOWLEDGMENTS: I would like to acknowledge Franci Gabrovšek, Matija Perne, Joe Myre, Wolfgang Dreybrodt, and Douchko Romanov for stimulating discussions and useful comments. This material is based upon work supported by the National Science Foundation under Grant No. EAR 1226903. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation. REFERENCES Birk, S., Liedl, R., Sauter, M. & G. Teutsch, 2005: Simulation of the development of gypsum maze caves.-Environmental Geology, 48, 3, 296-306. Blumberg, P. & R.L. Curl, 1974: Experimental and theoretical studies of dissolution roughness.- Journal of Fluid Mechanics, 65, 4, 735-751. Buhmann, D. & W. Dreybrodt, 1985: The kinetics of cal-cite dissolution and precipitation in geologically relevant situations of karst areas: 1. Open system.-Chemical Geology, 48, 1-4, 189-211. Covington, M.D., Luhmann, A.J., Wicks, C.M. & M.O. Saar, 2012: Process length scales and longitudinal damping in karst conduits.- Journal of Geophysical Research, 117, F1, 1-19. Covington, M.D., Prelovšek, M. & F. Gabrovšek, 2013: Influence of CO2 dynamics on the longitudinal variation of incision rates in soluble bedrock channels: Feedback mechanisms.- Geomorphology, 186, 85-95. Curl, R.L., 1966: Scallops and Flutes.- Transactions of the Cave Research Group of Great Britain, 7, 2, 121-160. Curl, R.L., 1974: Deducing flow velocity in cave conduits from scallops.- National Speleological Society Bulletin, 36, 2, 1-5. Dreybrodt, W., 1988: Processes in Karst Systems: Physics, Chemistry, and Geology.- Springer, pp. 288, New York, USA. Dreybrodt, W., 1996: Principles of early development of karst conduits under natural and man-made conditions revealed by mathematical analysis of numerical models.- Water Resources Research, 32, 9, 2923. Dreybrodt, W. & D. Buhmann, 1991: A mass transfer model for dissolution and precipitation of calcite from solutions in turbulent motion.- Chemical Geology, 90, 1-2, 107-122. Dreybrodt, W., Gabrovšek, F. & D. Romanov, 2005: processes of Speleogenesis: A Modeling Approach.- ZRC Publishing, pp. 376, Ljubljana, Slovenia. Goodchild, M. & D. Ford, 1971: Analysis of scallop patterns by simulation under controlled conditions.-Journal of Geology, 79, 1, 52-62. Hammer, 0., Lauritzen, S.E. & B. Jamtveit, 2011: Stability of dissolution flutes under turbulent flow.- Journal of Cave and Karst Studies, 73, 3, 181-186. Kaufmann, G., 2009: Modelling karst geomorphology on different time scales.- Geomorphology, 106, 1-2, 62-77. Liu, Z. & W. Dreybrodt, 1997: Dissolution kinetics of calcium carbonate minerals in H20-CO2 solutions in turbulent flow: the role of the diffusion boundary layer and the slow reaction H2O + CO2 ^ H+ + HCO-.- Geochimica Cosmochimica Acta, 61, 14, 2879-2889. Morse, J. & R. Arvidson, 2002: The dissolution kinetics of major sedimentary carbonate minerals.- Earth-Science Reviews, 58, 51-84. Perne, M., 2012: modelling speleogenesis in transition from pressurised to free surface flow.- Ph.D. thesis, University of Nova Gorica. Plummer, L., Wigley, T. & D.L. Parkhurst, 1978: The Kinetics of Calcite Dissolution in CO2- Water Systems at 5° to 60 °C and 0.0 to 1.0 ATM CO2.- American Journal of Science, 278, 179-216. Rehrl, C., Birk, S. & A.B. Klimchouk, 2008: Conduit evolution in deep-seated settings: Conceptual and numerical models based on field observations.- Water Resources Research, 44, 11, 1-13. Richardson, K. & P. Carling, 2005: A typology of sculpted forms in open bedrock channels. Special publication 392.- Geological Society of America, pp. 108, Boulder, Colorado. Rickard, D. & E. Sjöberg, 1983: Mixed kinetic control of calcite dissolution rates.- American Journal of Science, 283, 815-830. Szymczak, P. & A.J. Ladd, 2011: The initial stages of cave formation: beyond the one-dimensional paradigm.-Earth and Planetary Science Letters, 301, 3-4, 424432. APPENDIX: AN APPROXIMATION FOR THE PWP EQUATION The PWP equation is roughly a linear function of (Ceq - C) between 10% and 90% of saturation for a wide range of temperatures and PCo2 values. To facilitate future work with PWP rates, I present a fitting function for the rate constant as, log(aJ ^A + BJ+B2log(Pco2)+ + C2[log(i'co2)]' (14) where the best fit values of the parameters are given in Tab. 1. To determine this relationship, I calculated PWP rates for 100 calcium concentration values from 10% to 90% saturation for each of 100 different values of PCO2 sampled evenly in log space in the range 3x10-4 < PCO2 < 0.1 and for 11 different values of temperature in the range 0 °C < T < 24 °C. For each choice of t and PCO2 a best fit value of a^ for the relation Fpwp = as(Ceq - C) was calculated using least squares. This resulted in 1100 values of a^ for a range of T and PCO2 conditions. Typical residuals between the linear relation and the full PWP equation are about 10%. The maximum residuals for any of the fits are around 70% and occur near saturation for cases with high pCO2 . Parameter Best fit value A -4.30 B1 1.40 x 10-2 B2 0.150 C1 -2.83 x 10-5 C2 8.76 x 10-2 Table 1: Best fit parameters for the equation for the linear approximation to the PWP dissolution rate equation. The parameters in the best fit relation above (Eqn. 14) were determined via least squares fitting of these 1100 values of as. Eqn. 14 reproduces the best fit values of a^ within 10%. Therefore, this relation is not intended for precision work, but can allow quick estimation of PWP rates according to a simpler relation. In addition to the fitting relation, the Python code used to calculate PWP rates is available online at http://www.speleophysics.com. All calculations shown above employed the full PWP equation rather than the approximation given by Eqn. 14.