COBISS: 1.01 TIME SCALES IN THE EVOLUTION OF SOLUTION POROSITy IN POROUS COASTAL CARBONATE AqUIFERS By MIxING CORRO-SION IN THE SALTwATER-FRESHwATER TRANSITION ZONE. ČASOVNO MERILO RAZVOJA POROZNOSTI ZARADI KOROZIJE MEŠANICE V MEJNEM OBMOČJU SLADKOVODNIH LEČ V MEDZRNSKO POROZNEM KARBONATNEM OBALNEM VODONOSNIKU wolfgang DREyBRODT1 & Douchko ROMANOV2 Abstract UDC 556.3:552.54:539.217 Wolfgang Dreybrodt and Douchko Romanov: Time scales in the evolution of solution porosity in porous coastal carbonate aquifers by mixing corrosion in the saltwater-freshwater transition zone. Dissolution of calcium carbonate in the saltwater-freshwater mixing zone of coastal carbonate aquifers up to now has been treated by coupling geochemical equilibrium codes to a reactive-transport model. Te result is a complex nonlinear coupled set of diferential transport-advection equations, which need high computational eforts. However, if dissolution rates of calcite are sufciently fast, such that one can assume the solution to be in equilibrium with respect to calcite a highly simplifed modelling approach can be used. To calculate initial changes of porosity in the rock matrix one only needs to solve the advection-transport equation for salinity s in the freshwater lens and its transition zone below the island. Current codes on density driven fow such as SEAwAT can be used. To obtain the dissolution capacity of the mixed saltwater-freshwater solutions the calcium equilib-rium concentration ceq(s) is obtained as a function of salinity by PHREEqC-2. Initial porosity changes can then be calculated by a simple analytical expression of the gradient of the spatial distribution s(x, y) of salinity, the distribution of fow fuxes q(x,y) and the second derivative of the calcium equilibrium concentra-tion ceq(s) with respect to salinity s. Tis modelling approach is employed to porosity evolution in homogeneous and heterogeneous carbonate islands and coastal aquifers. Te geometrical patterns of porosity changes and the reasons of their origin will be discussed in detail. Te results reveal initial changes of porosity in the order of several 10-6 per year. Tis places the time scale of cavern evolution to orders from several tens of thousands to a hundred thousand years. Keywords: Calcite dissolution, mixing corrosion, saltwater-freshwater, mixing zone, coastal aquifer, evolution of porosity. Izvleček UDK 556.3:552.54:539.217 Wolfgang Dreybrodt and Douchko Romanov: Časovno merilo razvoja poroznosti zaradi korozije mešanice v mejnem območju sladkovodnih leč v medzrnsko poroznem karbonatnem obalnem vodonosniku Dosedanji modeli raztapljanja kalcijevega karbonata v območju mešanja sladke in slane vode temeljijo na združitvi geokemičnih ravnotežnih in reakcijsko transportnih modelov. Dobljeni sistem nelinearnih enačb zahteva veliko računske moči. Če je hitrost raztapljanja dovolj visoka in lahko predpostavimo, da je raztopina ves čas v ravnotežju glede na kalcit, rešimo problem z poenostavljenim modelskim pristopom. Začetno spreminjanje poroznosti v kamninski matriki določa advekcijsko tranportna enačbo, ki opisuje slanost v sladkovodni leči in prehodnem območju pod njo. Pri reševanju porabimo dostopne programske kode. Tokove nastale zaradi razlik v gostoti modeliramo s programom SEAwAT, topnost kalcita v mešanici sladke in slane vode v odvisnosti od slanosti pa izračunamo s programom PHREEqC-2. Začetno spreminjanje poroznosti lahko nato izračunamo z enostavnim analitičnim izrazom gradienta prostorske razporeditve slanosti s(x,y), razporeditve gostot toka q(x,y) in drugega odvoda ravnotežne koncentracije kalcija po slanosti. Tak modelski pristop uporabimo pri računanju razvoja poroznosti v homogenih in heterogenih karbonatnih otokih in obalnih vodonosnikih. Podrobno so prikazani vzroki in geometrijski vzorci spreminjanja poroznosti. Rezultati kažejo, da je začetna hitrost spremembe poroznosti reda velikosti 10-6 na leto. To postavi časovno merilo razvoja jam v območje nekaj deset tisoč do sto tisoč let. Ključne besede: Raztapljanje kalcita, korozija mešanice, območje mešanja sladke in slane vode, obalni vodonosnik, razvoj poroznosti. 1 Universitaet Bremen, FB1, Karst Processes Research Group, Bremen, Germany, e-mail: dreybrodt@ifp.uni-bremen.de 2 Freie Universitaet Berlin, Fachbereich Geowissenschafen, Berlin, Germany, e-mail: dromanov@zedat.fu-berlin.de Received/Prejeto: 21.12.2006 TIME in KARST, POSTOJNA 2007, 25–34 wOLFGANG DREyBRODT & DOUCHKO ROMANOV INTRODUCTION Carbonate islands consisting of porous rocks show typical karst features characterized by large dissolution chambers close to the coast, which have been created by mixing corrosion in the fresh-saltwater transition zone (Mylroie and Carew, 2000). Figure1 represents the basic concept. Due to meteoric precipitation a freshwater lens Fig. 1: Conceptual representation of a carbonate island from mylroie and Carew (2000). the chlorine concentration s, termed as chlorinity further on, of the mixture undersaturation or supersaturation may result. Figure 3 gives an example. It depicts the diference Fig. 3: Ace9(s) = ceg(s)-cmi[(j) as a function of chlorine concentration. Te curve extends from pure freshwater (right) to pure seawater (lef). builds up, foating on the denser saltwater (Vacher, 1988) Te transition from freshwater to seawater is not sharp. De-pending on many factors, such as tidal pumping, periodic-ity of annual recharge, and the heterogeneity of the rock’s properties in the aquifer it exhibits a transition zone. Tis zone can range from a few meters to half the depth of the lens. In this zone mixing between saltwater and freshwater activates mixing corrosion, which creates large chambers. Tese are called fank-margin caves. Figure 2 shows such a cave with its typical solutional features on its ceiling. Fig. 2: Flank-margin cave. when seawater mixes with a solution of H2O-CO2-CaCO3 saturated with respect to CaCO3 the mixture is no longer in equilibrium with respect to calcite. Depending on hceq(s) = ceq(s)-cmix(s) of the calcium concentration c ix (s) of the mixture and that of its corresponding equilibrium concentration c (s) as a function of s. Tis is the amount of calcium, which can be dissolved or precipitated, when the mixed solution is in contact with carbonate rock. Te HO-CaCO-CO solution used to calculate this data is in 2 32 equilibrium with a partial pressure of CO2 of 0.01 atm at a temperature of 20°C. Te seawater also is at 20°C. Te data in Fig. 3 were obtained by use of the code PHREEqC-2 (Parkhurst and Apello, 1999 ). From Figure 3 it is evident that mixtures with low content of seawater, chlorinity s ? 0.3 mol/^, can dissolve calcite, whereas mixtures with higher chlorinity may precipitate calcite. Renewed aggressivity due to mixing therefore occurs only at the freshwater side of the mixing zone where chlorinity is low. If one assumes that dissolution of calcite proceeds sufciently fast the solution there will be saturated with respect to calcite. Dissolution of minerals under such conditions is termed a gradient reaction (Phillips, 1991). Here we use this as a novel instrument to explain the evolution of porosity in carbonate islands. Dissolution rates of limestone are sufciently fast, such that afer mixing between saltwater and freshwater we assume saturation with respect to calcite in the entire lens. Afer attaining equilibrium the local distribution of calcium concentration c (s(x,z)) becomes stationary and exhibits gradients. Necessarily advection and difusion must transport the dissolved limestone to the outfow of the aquifer. 26 TIME in KARST – 2007 TIME SCALES IN THE EVOLUTION OF SOLUTION POROSITy IN POROUS COASTAL CARBONATE AqUIFERS ... DISSOLUTION IN THE MIxING ZONE Te advection term: In Figure 4 we consider a volume element dxdydz at position (x,y,z), into which fow, with components qx and qz, enters perpendicular to dydz or dxdy. Te fux q is defned by the volume of fuid per time unit entering through a unit of surface area and is given in [cm3/(cm2s) = cms-1]. afresh + (c, sea fresh J (5) Fig. 4: mass balance for the advection term. Te component qx transports solution from the neighbouring elementary cell at position (x-dx, y,z) via the area dydz into the cell dxdydz. Tis solution has al-ready attained equilibrium ceq(s(x-dx,y,z)) at position x-dx. when it enters into the volume element dxdydz it must dissolve or precipitate limestone to adjust its calcium concentration to equilibrium ceq(s(x,y,z)) at position x. On the other hand solution from the element dxdydz fows out into the neighbouring cell with fux qx(x,y,z). Mass conservation requires that the amount of limestone dis-solved per time unit in the element dxdydz must be equal to the diference of mass transported into the cell and that transported out of it. From this one fnds Ac«? is the increase of equilibrium concentration as given in Figure 3, s is chlorinity of seawater. sea Te difusion term: Our mass balance so far, however, is incomplete because gradients of c cause transport by difusion. Te rate qD of mass transport by difusion is given by Qd=~0'DW (cmlx+Aceq) (6) where D = qd/? + D is the m coefcient of dispersion. ? is the porosity of the rock and d its grain size. (Phillips, 1991). Dm is the constant of molecular diffusion (10-5cm s-1). Te total rate: Te total dissolution rate q is then given by qD+q d. ßw -qgradic^)- 0Z>V2(cmjJ)+ qgradfaj- 0Z)V (Ac^) (7) Due to the linearity of cmix with salinity s (eqn. 5) one fnds grad(cmix) proportional to grad(s). Te distribution of salinity is governed by the ad-vection-difusion equation ds I dt = qgrads - ®D(V)2s = 0 a (1) (qx(x, z) ¦ Ceq(x, z) - qx(x -dx,z)- Ceq(x - dx, z))dydz dxdydz An analogue equation exists for qz the amount of limestone dissolved by the fux component q entering via the surfaces (dx,dy). (8) because the distribution s is stationary. From the linearity of s with cmix we have dcmlIldt=qgradcmix-0D(y)2cm„= 0 (9) (q2 (x, z) • Ceq(x, z) - qz (x, z - dz) ¦ Ceq(x, z - dz))dxdy dxdydz Terefore Qadv = Q* + Qz = z)- dM (2) (3) Te total dissolution rate qtot is given by the master equation ßto(=?grarf(ACJ-0i5V2(ACe?) (10) Because the fux q follows the Darcy law of incompressible fuids, div(q) =0. Since ?ceq(s(x,z)) is a function of local distribution s(x,z) by diferentiating and using the chain rule, one fnds using equation 8 ßo*=? •L*•«<%„*+Ac^) (4) -0(qd/0+Dm)-(Vs(x,z))2 ds2 (11) whereby we have replaced c^ = cmtt + Aceq ¦ cmfe is the calcium concentration resulting from the mixing of sea- Tis master equation relates the amount of dissolved water and freshwater and is a linear function of Cl-con- material per unit volume of the rock matrix [mol cm-3 s-1] centration s. TIME in KARST - 2007 27 , , wOLFGANG DREyBRODT & DOUCHKO ROMANOV to the gradient of salinity s, to the second derivative of ceq(s) = cmjx(s)+Aceq(s), and the fux q. d2Aceq/ds2 can be obtained by diferentiating twice the data set of Figure 3. Tis data set was obtained by using the program PHRE-EqC2 and calculating about 50 closely spaced points to avoid numerical errors, when diferentiating twice. Te result is shown in Figure 5. Fig. 5: Second derivative 3\ ds2 Te function |As(jc)| and the fow distribution \q(x)\ can be obtained by the numerical hydrologic model SEA-wAT, as will be shown in the next sections. To calculate the initial change of porosity it is sufcient to obtain the fux and salinity distribution of an island without considering calcite dissolution, because the time to establish a stationary state of the lens is in the order of 100 years. It is a good approximation to assume that during this time the change of porosity is insignifcant. Equation. 11 can be written in terms of the change of porosity as dt p Qlo,=-0D(fs(x,z)j d2Ac^ M ds2 p (1/s) (12) M = 100 g/mol is the molecular weight at CaCO3, r=2.7 g/cm3 is the density of compact CaCO3. q the mass of CaCO3 dissolved per time from a unit volume of the rock matrix is given in mol s-1 cm-3. dipl dt is the amount of volume dissolved per time from a unit volume of the rock matrix (cm3s-1/cm3). By use of equation 12 it is now possible to construct a conceptual frame for the evolution of porosity. Tests of this approach on simple benchmark models have shown its reliability and have found agreement to experimental data (Romanov and Dreybrodt, 2006). INITIAL CHANGES OF POROSITy IN A HOMOGENEOUS ISLAND. To obtain the initial distributions of fux q and chlorinity s in the lens of a carbonate island we have used SEAwAT by USGS (Guo and Langevin, 2002). Te modeling domain is shown in Figure 6. Te island is a strip of 1 km Fig. 6: modeling domain of a carbonate island. width. Porosity 0 and the hydraulic conductivity K are uniform (0=O.3O,K = l0m/day). Te transversal disper-sivity is aT = d = 0.01 cm, the longitudinal dispersivity is 28 TIME in KARST - 2007 al = 0.1 cm. Infltration is 3 . 10-3m/day =1.11m/year. Tis way the maximal depth of the lens is about 50 m below sea level. Te lower border of the domain reaches down to 70 m. At that boundary an impermeable layer imposes no-fow conditions. Te grid size in the domain is 1 m x 1 m in the part below sea level. In the part above sea level (2 m) the grid size is 0.2 m by 1 m. In its initial state when the island emerges out from the sea the entire aquifer is flled with seawater. when the island receives recharge from meteoric freshwater the lens builds up. A stable stationary lens is obtained afer about 30 years. Fig.7 shows the results of the model run. Figure 7a shows the freshwater lens (white), the transition zone and its distribution of Cl-concentration by a color code. From this distribution of chlorinity one can extract the scalar value Vs(jc) and d2ceq(s(x))/d s2. Figure 7b shows the chlorinity in units normalized to its maximum value along several horizontal sections as depicted in Figure 7a. Te lowest section at -68 m is entirely in saltwater with maximum Cl-concentration. Te section at -55 m extends through the almost horizontal base of . TIME SCALES IN THE EVOLUTION OF SOLUTION POROSITy IN POROUS COASTAL CARBONATE AqUIFERS ... Fig. 7: homogeneous island. a) Local distribution of chlorinity s(x). Te white region designates the freshwater lens. b) chlorinity along horizontal sections as indicated in a). the lens and shows a wide zone where the concentration raises to that of seawater. Te upper sections cut through the mixing zone and there the rise in concentration from freshwater to seawater becomes steeper. Te square of the gradient |vs| is shown by Figures 8a,b also normalized to its maximum value in Figure 8b. Figure 8a illustrates its local distribution, which exhibits large values only in the region of the transition zone. Te horizontal distribution along horizontal sections is depicted in Figure 8b. Te second derivative d2c (s(x))/ d s2 obtained from the Cl-concentration in Figure 7a is given in Figure 9a. Its distribution is limited to that part of the transition zone with 0 < s < 0.03 mole/^. See Figure 5. Tis corresponds to a narrow fringe at the freshwater side of the transition zone with seawater content from zero up to about 4%. In any case creation of porosity is possible only in this restricted region. Figure 9b for completeness depicts some distributions of d2c (s(x))/d s2 along horizontal sections. To calculate the initial rate of change in porosity (conf. eqn 12) the Darcy fuxes q must be known. Tey are also obtained from the model run and shown in Fig 10. Te fux is low in the center of the island q^l m/year), but increases by orders of magnitude when the fuid Fig. 8: homogeneous island. a) Local distribution of the square of gradients |Vs(x)| , b) square of gradients along horizontal sections as indicated in a). moves coastward, where it becomes about 0.2 m/day at the outfow. Te dispersion coefcient D = qd/0 + D (conf. eqn. 11) depends on the fux q, but also on the coef-cient of molecular difusion D =105 cm2/s. For low fux m q<10 4 cms1 and particle diameters d?10-2 cm dispersion is dominated by molecular difusion. In the following scenarios we have used d=10-2 cm, a realistic value in porous limestone. Terefore in the range of fux, which can be read from Figure 10b the dispersion coefcient in the center of the island is D=10-5 cm2s-1. It increases by about 60% of this value at the coast. From the data given in Figures 7a, 8a, and 9a the initial porosity is obtained by use of eqn. 12. Figure 11 illustrates these results. Changes in porosity are restricted to a small fringe in the transition zone and are fairly even along it. Tey are in the order of 10-6 year1. Tis is sufcient to create substantial porosity within 100,000 years. At the outfow fank margin caves can develop in 10,000 years. One has to keep in mind, however, that the approximation as a homogeneous island is a high idealization. Any disturbances, which increase the width of the transition zone, will reduce the gradients of chlorinity and therefore on more realistic settings the initial porosity changes accordingly. TIME in KARST - 2007 29 wOLFGANG DREyBRODT & DOUCHKO ROMANOV Fig. 9: homogeneous island. Local distribution of the second derivative —j, b) second derivative along horizontal sections as indicated in a). As we have stated already, the second derivative is restricted to narrow regions in the freshwater side of the transition zone. It exhibits signifcant values only at locations where the water contains between zero and 4% saltwater (see Figure 5). On the other hand the gradient in salinity is maximal at mixtures of about 50% seawater, because it arises from a difusive process. In the region of maximal gradients, however, the second derivative is small. Vice versa in the region of high values of the second derivative, the gradients of salinity are low. Tis is illustrated in Figure 12. Tis fgure is an overlay of the horizontal distributions (grads)2 in Figure 8b (red curves), the second derivative in Figure 9b (green curves), and the initial porosity change in Figure 11b (black curves). All curves are normalized to their individual maximum values. Terefore their values are not comparable in this fgure. what can be compared, are the locations. Evidently the curves for gradients and second derivative are well separated. Te curves of porosity change are proportional to the product of the square of the gradient and the second derivative. Porosity change displays high values Fig. 10: homogeneous island. a) Local distribution of fux b) fux along horizontal sections as indicated in a). in between their maxima but close to the region of high values of the second derivative. Figure 13 further illustrates this qualitatively. Te red region depicts the locations of the modeling domain where (grads)2 exhibits values val a 10"2 val^ valma is the maximal value. Te green region shows these locations for the second derivative and fnally the black region shows the locations of signifcant changes of porosity. Tese fndings agree with those of Sanford and Konikow (1989) who also found that changes in porosity are restricted to regions where waters contain between 0.5% and 3% of seawater. It should be noted here that any mechanism, which changes the sigmoid shape of the salinity distribution to a linear profle would enhance evolution of porosity dramatically. In this case salinity gradients become constant in the entire mixing zone and their value is at least one order of magnitude higher at the maximal value of the second derivative. One could speculate that tidal pumping and fuctuations of the water table due to seasonal changes of infltration could cause such linear mixing zones. Present observations in boreholes give some evidence for such transition zones. 30 TIME in KARST – 2007 TIME SCALES IN THE EVOLUTION OF SOLUTION POROSITy IN POROUS COASTAL CARBONATE AqUIFERS ... Fig. 11: homogeneous island. a) Local distribution of initial change of porosity dldt. Fig. 20: heterogeneous island. Regions of high values of (grads)2 (red), -----r^ (green), and change of porosity (black) in the modeling domain. the freshwater lens. Tis can be also visualized from the second derivatives as shown in Figure 18. Figure 19 illustrates the initial change of porosity, which exhibits high values of 3 . 10-6 year1 (red) at only a few locations close to the freshwater side of the transition TIME SCALES IN THE EVOLUTION OF SOLUTION POROSITy IN POROUS COASTAL CARBONATE AqUIFERS ... zone. At some favorable locations (red and yellow) caves ure 21 depicts (grads)2 (red), d2ce /ds2 (green), and d(j>/dt may evolve there in several 10,000 to 100,000 years. (black) along selected horizontal sections. Tis is further illustrated by Figure 20, which shows In both fgures we fnd that the regions of (grads)2 the regions of high values for (grads)2 (red), d2c /ds2 (red), d2ceq /ds2 (green) are well separated and porosity (green), and dldt (black) in the modeling domain. Fig- develops in between. Due to the heterogeneity, however, the patterns become complex. Fig. 21: a) heterogeneous island. distributions of (grads)2 (red), ----^ (green), and porosity change (black) along selected horizontal sections. Number on the sets of curves give the depth of the section. INITIAL CHANGES OF POROSITy IN SALTwATER TONGUES. when impermeable strata underlay an island sufciently close to its surface the freshwater lens cannot extend be-low this layer and a saltwater tongue intrudes from the coastland inward until it reaches the impermeable layer. From thereon the freshwater lens is truncated by this layer. In this situation mixing of waters is restricted to the transition zone of the tongue and one expects high dissolution rates in this region. Fig. 22 shows the local distribution of chlorinity and the initial change of porosity using the statistical distribution in Figure 14 for the upper permeable part. Te mixing zone exhibits a structure similar to that of the heterogeneous island at the corresponding loca-tions. Porosity changes at the outfow are low, but we fnd values up to 10-6 1/year land inward at various lo- Fig. 22: Coastal aquifer with heterogeneous conductivity down to 29 m as used in Fig. 14. Te strata below 29 m are impermeable (grey). Local distribution of Cl-concentration s(x) and initial porosity change d