PAGE 1
LETTER: CALCITE DISSOLUTION UNDER TURBULENT FL OW CONDITIONS: A REMAINING CONUNDRUM P ROBLEM NASTANKA FASET IN RAZTAPLJANJA KALCITA V TURBULENTNEM TOKU Matthew D. C OVINGTON 1 UDK 551.435.8:549.742.111 1 University of Arkansas Department of Geosciences Fayetteville, Arkansas, USA, email: mcoving@uark.edu Received/Prejeto: 26.11.2013 COBISS: 1.01 ACTA CARSOLOGICA 43/1, 195, POSTOJNA 2014 I NTRODUCTION A signicant 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). e relationships derived from this work have been applied extensively to the develop ment of mechanistic models of speleogenesis (e.g. Drey brodt 1996; Dreybrodt et al. 2005; Birk e t 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 ow conditions. Recent e orts have begun to develop mechanistic models for processes governing later stages of cave evolution, considering factors such as turbulent ow structures (Hammer et al. 2011) and open channels (Perne 2012). However, such studies remain limited, in part due to signicant quantitative uncertainties in a va riety of processes that become important beyond the in cipient speleogenesis stage (Covington et al. 2013). While speleogenetic models have not typically been run much beyond the transition from laminar to turbulent ow conditions, experimental and theoretical studies of carbonate dissolution have constrained dis solution rates and ratecontrolling mechanisms under turbulent ow (e.g. Rickard & Sjberg 1983; Dreybrodt & Buhmann 1991; Liu & Dreybrodt 1997). However, a direct application of these results and comparison to eld observations leads to an apparent conundrum. Ac cording to the theory, solutional forms such as scallops and utes should not exist in limestone; however, they clearly do exist at a wide variety of sites and scales. is 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 ow conditions are typically calculated by dividing the uid into two re gions, a turbulent core that is wellmixed, and a di usion boundary layer (DBL) that lies between the turbulent core and the dissolving wall (Dreybrodt & Buhmann 1991). Within the DBL, ow is presumed to be laminar and species cross the layer via the process of molecular di usion. ere are two endmember regimes in which overall dissolution rates are limited by either di usion rates (the transportlimited regime),
PAGE 2
ACTA CARSOLOGICA 43/1 â€“ 2014 196 MATTHEW D. C OVINGTON or surface reaction rates (the reactionlimited regime), where F D and F s are the transport, and reactionlimit ed uxes, respectively, D is the di usion constant, C eq is the equilibrium concentration of calcium, C b is the concentration of calcium in the bulk solution, is the DBL thickness, and k 1 is the kinetic rate constant for cal cite dissolution in the linear kinetic regime. e mixed kinetic regime occurs when the di usionlimited and surfacereactionlimited rates are similar (Rickard & Sjberg 1983). In this case, both processes are important. For simplicity, I only consider linear dissolution kinetics here, as is common below about 80% saturation, though our conclusions are not sensitive to this choice. In the case of a suciently thick DBL, or a suciently thin tur bulent core, conversion of CO 2 can also be ratelimiting (Dreybrodt & Buhmann 1991). ough nonlinear ki netics and CO 2 conversion limitation are not explicitly considered here, if anything, they would exacerbate the presented conundrum by providing further means to re duce 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 diusion rates for a wide variety of head gradients and hydraulic diameters using the DarcyWeisbach equa tion and ColebrookWhite relation. For these calcula tions, the ducial parameter values given in Dreybrodt et al. (2005) were used, with k 1 = 4 10 mol cm s , D = 10 cm 2 s , and C eq = 2 10 mol cm . For ease of comparison, Covington et al. (2012) converted the â€œsur faceâ€ rate law to the form F s = s ( C eq C ), (3) where s = 2 10 cm s was determined by dividing k1 by C eq . Here I have slightly modied the notation of Covington et al. (2012), replacing with s to clarify that the constant applies to the chemical reaction rate at the surface. In many previous manuscripts (e.g. Buhmann & Dreybrodt 1985) without a subscript represents the combined e ects of di usion and reaction. Covington et al. (2012) compared the rates predict ed by Equations 1 and 3 by nondimensionalizing the two equations by dividing by the term s ( C eq C ), such that the dimensionless dissolution rates are independent of chemistry and become F s = 1, (4) and where the stars denote dimensionless rates. en was calculated for a wide variety of head gradients and hy draulic diameters to examine the relative magnitude of the reactionand transportlimited equations. e 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 calcu lations a fractional roughness of 0.05 was assumed. e 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 ow F ig. 1: Di usional (solid) and surface reaction (dashed) dissolu tion rates as a function of hydraulic diameter, D H , depicted for di erent hydraulic head gradients, h. Sharp jumps in di usion rates occur at the onset of turbulent ow. e laminar/turbulent transition is indicated with dotted lines. F or the parameter space shown, the theory suggests that di usional rates only control dis solution for lowgradient conduits just below the laminar/turbu lent transition. Elsewhere, di usional rates are much higher than surface reaction rates. F igure reproduced with permission from Covington et al. (2012). F ig. 2: Scallops formed on a limestone surface within a cave, in dicating contrasts in dissolution rates as a result of turbulent ow structures. M ajor divisions on the ruler are in cm (P hoto: M atija P erne).
PAGE 3
ACTA CARSOLOGICA 43/1 â€“ 2014 197 LETTER CALCITE DISSOLUTION UNDER TURBULENT FL OW CONDITIONS: A REMAINING CONUNDRUM 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 di usive rates are predicted to be for the majority of the turbulent region of the parameter space (Fig. 1). For the convenience of future studies, a tting function for PWP rates as a function of tempera ture and P CO 2 is provided in the Appendix. is relation is not used explicitly here, but provides a quick means of estimating PWP rates in a given setting. S URFACE RATES FROM THE PWP E QUATION e above analysis makes use of ducial kinetic constants from Dreybrodt et al. (2005). However, these constants are motivated by typical values from experiments, where the e ects of transport and reaction are both present. erefore, a more accurate approach is to use the full PlummerWigleyParkhurst (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 re sult from a temperature and partial pressure of CO 2 that approximately reproduce the ducial value of C eq = 2 10 mol cm used above. is equilibrium concentration is roughly obtained (for an open system) with T = 10 C and P CO 2 = 0.01 atm, values which are also quite typical in the cave environment. e PWP rates are depicted in Fig. 3, alongside the simplied lin ear relationship with two di erent values of s . e 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 ducial value of s used to create Fig. 1 above, with the bestt value for the linear region of the PWP equation (Fig. 3), suggests that the value of s 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 shied up by a similar factor of 3 to 4. While this does imply that the surface rates are higher F ig. 3: e dissolution rate as a function of Calcium concentra tion for PCO 2 = 0.01 atm and T = 10 C. is example is con structed to have a similar value of C eq and temperature to the ducial case from Dreybrodt et al. (2005). C OMPARING SURFACEREACTION AND DIFFUSION RATES e equations given above for di usionlimited (Equa tion 1) or surfacelimited rates (Equation 3) are only val id under conditions where dissolution is truly limited by the respective process. In steady state, di usion and sur regime show transportlimited rates. In fact, di usion rates under turbulent ow are typically several orders of magnitude larger than the surface reaction rate. is conclusion leads to the conundrum. e model suggests that dissolution rates should almost always be controlled by the surface reaction rate in the turbulent ow regime, and, therefore, that dissolution rates should be independ ent of any spatial variations in DBL thickness. On the con trary, caves and channels formed in limestone frequently contain scallops (Fig. 2) and utes that apparently form as the result of systematic contrasts in DBL thickness as a result of turbulent ow structures (Curl 1966; Good child & Ford 1971; Blumberg & Curl 1974; Curl 1974). e remainder of this note will explore this conundrum in more detail, and discuss possible resolutions.
PAGE 4
ACTA CARSOLOGICA 43/1 â€“ 2014 198 face reaction rates must be equal, and the rate of di usion is given by where C s is the concentration at the surface, and C b is the concentration in the bulk ow. For the di usionlimited case C s C eq , and we recover Equation 1 above. Similarly, surface reaction rates are dependent on the concentra tion at the surface and are given by F s = s ( C eq C s ). (7) Again, if dissolution is reactionlimited C s C b and we recover Equation 3. Setting di usion and reaction rates equal, solving for C s , and plugging this back into Equa tion 7, one can see that the dissolution rate accounting for both reaction and di usion is given by where d = D / . e lines depicting the di usionlimited rates in Fig. 1 are equivalent to d / s , so the analysis and conclusions of Covington et al. (2012) can also be cast in terms of the relative magnitudes of d and s . ere fore, if dissolution rates are wellrepresented by the linear relation in Equation 7, then the relative importance of di usion and reaction in determining the overall rates can be quantied by comparing d and s . A critical di usion boundary layer thickness, at which di usion and reaction are equally important can be calculated by setting d = s , leading to crit = D / s . (9) When >> crit then rates are di usionlimited, and when << crit rates are reactionlimited. If ~ crit then disso lution will occur via mixed kinetics. At = crit di usion surpresses dissolution rates to half of what they would be from surface reaction alone. e 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 di usionlimited equation can be applied, while in the limit of thin DBL the PWP surface rate can be applied. An estimate of crit can be obtained by extrapolating the di usionlimited equation until it intercepts the PWP surface rate (Fig. 4). Quanti tatively, this results in the relation 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 crit then dissolution rates will be limited by either di usion or surface reac tion, respectively. Using this relation, and the full PWP equation, one can calculate critical DBL thicknesses for a range of temperatures, P CO 2 values, and calcium con centrations. is calculation (Fig. 5) shows that the criti cal DBL thickness is on the order of magnitude of 1 mm for essentially the entire range of temperature, P CO 2 , and dissolved load found in natural karst systems, except in highly undersaturated conditions where C 0 and crit << 1 mm. F ig. 4: e critical D BL thickness, below which dissolution rates are surface reactionlimited, can be estimated by extrapolating the di usionlimited rate until it intersects the surfacelimited rate determined by the PWP Equation. is leads to the rela tion in Equation 10. ough this approximation scheme can be employed more generally, the lines plotted in this gure represent exactly the case where the rates can be approximated using D and s . F ig. 5: e critical D BL thickness, below which dissolution rates are surface reactionlimited, shown for a wide range of temper ature, P CO 2 values, and dissolved loads. Lines terminate on the righthand side at saturation. MATTHEW D. C OVINGTON
PAGE 5
ACTA CARSOLOGICA 43/1 â€“ 2014 199 C ONCLUDING THOUGHTS AND POTENTIAL RESOLUTIONS A more careful analysis reinforces the conundrum. e theory suggests that limestone dissolution rates under turbulent ow conditions should be typically limited by surface reaction rates. However, eld observations of scallops and other solutional forms clearly suggest other wise. A recent model developed using PWP dissolution equations and computational uid dynamics also found that limestone dissolution utes were unstable under the calculated dissolution rates, though utes in gypsum, which has much higher surface rates, proved to be some what more stable (Hammer et al. 2011). Supercially, this result is in agreement with the analysis presented here, and may provide further evidence of difculties with the current theory. e 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 tur bulent dissolution model is the semiempirical equa tion used to calculate Sherwood Number (Eqn. 13), and consequently DBL thickness (Eqn. 12). is 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 utes in gypsum suggest a typical DBL thickness of = 0.0089 L Sc /3 , where L is the ute length and Sc is the Schmidt Number, with Sc 1000 for water at the relevant temperatures. is results in a DBL thickness that is roughly 0.1% of the ute 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. is again is in conict with eld observations, where scallops and utes oen form on length scales of centimeters (Fig. 2). e approximation crit ~ 1 mm can also be repro duced using a rough estimate of PWP rates. In most natural settings, the forward reaction rate is dominat ed by reaction III ( 3 in the PWP equation) (Dreybrodt 1988). is rate is simply a function of temperature and is on the order of 10 mol cm s (Dreybrodt 1988, p. 127). Critical DBL thickness can then be esti mated as which is quite similar to the result of 1 mm obtained from the larger parameter study (Fig. 5). Within speleogenesis models DBL thickness is usu ally estimated from a relationship for pipe ow that em ploys the Sherwood Number (Sh), = D H /Sh, (12) where Sh is given by an empirical relationship such as (Dreybrodt et al. 2005). Here, Re is the Reynolds number, Sc is the Schmidt number, and f is the Cole brookWhite friction factor. Using these relations, one can calculate typical values of DBL thickness, ( ), for a selection of head gradients and hydraulic diameters (Fig. 6). Again, for this calculation a relative pipe rough ness of 0.05 is assumed, though the result is not particu F ig. 6: D BL thickness ( ) under turbulent ow conditions as esti mated from Sherwood number for a wide range of hydraulic pa rameters. Each line represents a choice of head gradient. e lines terminate on the le at Re = 4000, where the ColebrookWhite formulation breaks down. T ypical values of D BL thickness as esti mated from Sherwood Number are much less than the 1 mm mag nitude of the critical value determined from the P W P equations. LETTER CALCITE DISSOLUTION UNDER TURBULENT FL OW CONDITIONS: A REMAINING CONUNDRUM 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 mag nitude estimate of a critical DBL thickness of 1 mm. is is true for all turbulent ow cases except those at very low head gradients (i.e. 10 ).
PAGE 6
ACTA CARSOLOGICA 43/1 â€“ 2014 200 A CKNOWLEDGMENTS: I would like to acknowledge Franci Gabrovek, Matija Perne, Joe Myre, Wolfgang Dreybrodt, and Douchko Romanov for stimulating discussions and useful com ments. is material is based upon work supported by the National Science Foundation under Grant No. EAR 1226903. Any opinions, ndings, and conclusions or rec ommendations expressed in this material are those of the author and do not necessarily reect the views of the Na tional Science Foundation. Another possibility is that the rates estimated by the PWP equation are signicantly 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). ese 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 appar ent increase in surface reaction rates that would reduce the critical DBL thickness. However, the observed dis crepancy of a factor of two is not sucient to explain away the conundrum by itself. If natural surfaces are still signicantly rougher at scales below the DBL thickness than the disks used in experiments, this could potential ly resolve the problem. Other factors, such as microbial lms, could also inuence surface reaction rates, and perhaps in some cases increase them. One might ask whether a very small contrast in dis solution rates would be sucient to form utes and scal lops. However, ute experiments show that dissolution rates vary by roughly a factor of two over the length of a ute (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 sup pressed by roughly a factor of two. A thinner DBL would not allow sucient contrast in rates. Forms somewhat similar to scallops (mechanical erosion utes) also are found in bedrock channels in relatively insoluble rocks, though the variety and prev alence of such forms is greatest in highly soluble rocks (Richardson & Carling 2005). One possible explanation is that socalled solutional forms, such as scallops and utes, 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 suciently small (Fig. 5). However, the authorâ€™s experi ence would suggest that scallops are also present in lo cations without such highly aggressive water. Little sys tematic attempt has been made to study the locations in which scallops form, and whether forms di er accord ing to hydrological or lithological settings. Such studies might also provide clues as to the correct resolution of the current conundrum. Birk, S., Liedl, R., Sauter, M. & G. Teutsch, 2005: Simu lation of the development of gypsum maze caves.Environmental Geology, 48, 3, 296. Blumberg, P. & R.L. Curl, 1974: Experimental and theo retical studies of dissolution roughness.Journal of Fluid Mechanics, 65, 4, 735. Buhmann, D. & W. Dreybrodt, 1985: e kinetics of cal cite dissolution and precipitation in geologically relevant situations of karst areas: 1. Open system.Chemical Geology, 48, 1, 189. MATTHEW D. C OVINGTON REFERENCES 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. Covington, M.D., Prelovek, M. & F. Gabrovek, 2013: Inuence of CO 2 dynamics on the longitudinal variation of incision rates in soluble bedrock chan nels: Feedback mechanisms.Geomorphology, 186, 85.
PAGE 7
ACTA CARSOLOGICA 43/1 â€“ 2014 201 Liu, Z. & W. Dreybrodt, 1997: Dissolution kinetics of cal cium carbonate minerals in H 2 0â€“CO 2 solutions in turbulent ow: the role of the di usion boundary layer and the slow reaction H 2 O + CO 2 H + + HCO 3 .Geochimica Cosmochimica Acta, 61, 14, 2879. Morse, J. & R. Arvidson, 2002: e dissolution kinetics of major sedimentary carbonate minerals.EarthScience Reviews, 58, 51. Perne, M., 2012: M odelling speleogenesis in transition from pressurised to free surface ow.Ph.D. thesis, University of Nova Gorica. Plummer, L., Wigley, T. & D.L. Parkhurst, 1978: e Ki netics of Calcite Dissolution in CO 2 Water Systems at 5 to 60 C and 0.0 to 1.0 ATM CO 2 .American Journal of Science, 278, 179. Rehrl, C., Birk, S. & A.B. Klimchouk, 2008: Conduit evo lution in deepseated settings: Conceptual and nu merical models based on eld observations.Water Resources Research, 44, 11, 1. Richardson, K. & P. Carling, 2005: A typology of sculpted forms in open bedrock channels. Special P ublication 392.Geological Society of America, pp. 108, Boul der, Colorado. Rickard, D. & E. Sjberg, 1983: Mixed kinetic control of calcite dissolution rates.American Journal of Sci ence, 283, 815. Szymczak, P. & A.J. Ladd, 2011: e initial stages of cave formation: beyond the onedimensional paradigm.Earth and Planetary Science Letters, 301, 3, 424â€“ 432. LETTER CALCITE DISSOLUTION UNDER TURBULENT FL OW CONDITIONS: A REMAINING CONUNDRUM Curl, R.L., 1966: Scallops and Flutes.Transactions of the Cave Research Group of Great Britain, 7, 2, 121. Curl, R.L., 1974: Deducing ow velocity in cave conduits from scallops.National Speleological Society Bul letin, 36, 2, 1. Dreybrodt, W., 1988: P rocesses in Karst Systems: P hys ics, Chemistry, and G eology. Springer, pp. 288, New Y ork, USA. Dreybrodt, W., 1996: Principles of early development of karst conduits under natural and manmade condi tions revealed by mathematical analysis of numeri cal 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 Ge ology, 90, 1, 107. Dreybrodt, W., Gabrovek, F. & D. Romanov, 2005: P ro cesses of Speleogenesis: A M odeling Approach.ZRC Publishing, pp. 376, Ljubljana, Slovenia. Goodchild, M. & D. Ford, 1971: Analysis of scallop pat terns by simulation under controlled conditions.Journal of Geology, 79, 1, 52. Hammer, ., Lauritzen, S.E. & B. Jamtveit, 2011: Stability of dissolution utes under turbulent ow.Journal of Cave and Karst Studies, 73, 3, 181. Kaufmann, G., 2009: Modelling karst geomorphology on di erent time scales.Geomorphology, 106, 1, 62.
PAGE 8
ACTA CARSOLOGICA 43/1 â€“ 2014 202 e PWP equation is roughly a linear function of ( C eq C ) between 10% and 90% of saturation for a wide range of tem peratures and P CO 2 values. To facilitate future work with PWP rates, I present a tting function for the rate constant s , where the best t 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 di erent values of P CO2 sampled evenly in log space in the range 3 P CO2 0.1 and for 11 di erent values of temperature in the range 0 C T 24 C. For each choice of T and P CO2 a best t value of s for the relation F pwp = s ( C eq C ) was calculated using least squares. is resulted in 1100 values of s for a range of T and P CO2 conditions. Typical residuals between the linear relation and the full PWP equation are about 10%. e maximum residuals for any of the ts are around 70% and occur near saturation for cases with high P CO2 . e parameters in the best t relation above (Eqn. 14) were determined via least squares tting of these 1100 val ues of s . Eqn. 14 reproduces the best t values of s within 10%. erefore, 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 tting 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. A PPENDI X: A N APPRO X IMATION FOR THE PWP EQUATION T able 1: B est t parameters for the equation for the linear approximation to the PWP dissolution rate equation. Parameter Best t value A .30 B 1 1.40 10 B 2 0.150 C 1 .83 10 C 2 8.76 10 MATTHEW D. C OVINGTON
