Dissolution of minerals with rough surfaces

Accepted Manuscript Dissolution of minerals with rough surfaces Thiago A. de Assis, Fábio D.A. Aarão Reis PII: DOI: Reference:

S0016-7037(18)30108-X https://doi.org/10.1016/j.gca.2018.02.026 GCA 10669

To appear in:

Geochimica et Cosmochimica Acta

Received Date: Accepted Date:

14 October 2017 16 February 2018

Please cite this article as: de Assis, T.A., Aarão Reis, F.D.A., Dissolution of minerals with rough surfaces, Geochimica et Cosmochimica Acta (2018), doi: https://doi.org/10.1016/j.gca.2018.02.026

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

2

Dissolution of minerals with rough surfaces Thiago A. de Assis1 and F´abio D. A. Aar˜ao Reis2 1 Instituto de F´ısica, Universidade Federal da Bahia, Rua Bar˜ao de Jeremoabo s/n, 40170-115, Salvador, BA, Brazil [email protected] 2 Instituto de F´ısica, Universidade Federal Fluminense, Avenida Litorˆanea s/n, 24210-340 Niter´oi, RJ, Brazil [email protected]ff.br

3

February 5, 2018

4

Abstract

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

We study dissolution of minerals with initial rough surfaces using kinetic Monte Carlo simulations and a scaling approach. We consider a simple cubic lattice structure, a thermally activated rate of detachment of a molecule (site), and rough surface configurations produced by fractional Brownian motion algorithm. First we revisit the problem of dissolution of initial flat surfaces, in which the dissolution rate rF reaches an approximately constant value at short times and is controlled by detachment of step edge sites. For initial rough surfaces, the dissolution rate r at short times is much larger than rF ; after dissolution of some hundreds of molecular layers, r decreases by some orders of magnitude across several time decades. Meanwhile, the surface evolves through configurations of decreasing energy, beginning with dissolution of isolated sites, then formation of terraces with disordered boundaries, their growth, and final smoothing. A crossover time to a smooth configuration is defined when r = 1.5rF ; the surface retreat at the crossover is approximately 3 times the initial roughness and is temperature-independent, while the crossover time is proportional to the initial roughness and is controlled by step-edge site detachment. The initial dissolution process is described by the so-called

1

rough rates, which are measured for fixed ratios between the surface retreat and the initial roughness. The temperature dependence of the rough rates indicates control by kink site detachment; in general, it suggests that rough rates are controlled by the weakest microscopic bonds during the nucleation and formation of the lowest energy configurations of the crystalline surface. Our results are related to recent laboratory studies which show enhanced dissolution in polished calcite surfaces. In the application to calcite dissolution in alkaline environment, the minimal values of recently measured dissolution rate spectra give rF ∼ 10−9 mol/(m2 s), and the calculated rate laws of our model give rough rates in the range 10−6 -10−5 mol/(m2 s). This estimate is consistent with the range of calcite dissolution rates obtained in a recent work after treatment of literature data, which suggests the universal control of kink site dissolution in short term laboratory works. The weak effects of lattice size on our results also suggest that smoothing of mineral grain surfaces across geological times may be a microscopic explanation for the difference of chemical weathering rate of silicate minerals in laboratory and in the environment.

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

42

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

1

Introduction

Mineral dissolution is important for several environmental and industrial processes, which explains the large number of efforts to determine the relation between the dissolution rate and the physical and chemical conditions involved in that process. However, the connection of laboratory studies, field observations, and modeling approaches is frequently difficult due to the discrepancy between estimates of the dissolution rates (Fischer et al., 2014). With the advance in experimental methods in the last decades, it was shown that the dissolution rate may have significant variations in different facets of the same solid; etch pits, grain boundaries, and other defects may also have significant contributions to the rate variability (Liang and Baer, 1997; Jordan and Rammensee, 1998; Gautier et al., 2001; Luttge et al., 2003; Ruiz-Agudo and Putnis, 2012; Levenson and Emmanuel, 2013; Smith et al., 2013; Godinho et al., 2014; Fischer et al., 2014; Emmanuel and Levenson, 2014; Laanait et al., 2015; Pollet-Villard et al., 2016a,b; Feng et al., 2017; Zareeipolgardani et al., 2017; Fischer and Luttge, 2017; Michaelis et al., 2017; Saldi et al., 2017). In many cases, etch pits were shown to dissolve much faster than the surrounding planar surfaces, such as in calcite (Liang and 2

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

Baer, 1997; Fischer et al., 2012; Smith et al., 2013), dolomite (Luttge et al., 2003; Saldi et al., 2017), and orthoclase (001) (Pollet-Villard et al., 2016b). However, this is not a universal feature; for instance, etch pits dissolve with approximately the same rate of other parts of the surface of quartz (Gautier et al., 2001) and in the long term calcite dissolution (Smith et al., 2013). A recent advance is the measurement of local dissolution rates in different regions of a single sample, which provides dissolution rate spectra (Fischer and Luttge, 2017; Michaelis et al., 2017). This approach confirms the observation of previous works that dissolution of low energy surfaces (i. e. in which all molecules are strongly bonded) is slow when compared with surfaces with rough steps, kink sites, and other high energy configurations. These results help to understand whether a given morphological feature will have lower or higher contribution to an average rate. There are also efforts to determine common features in a large number of measurements of dissolution rates. For instance, some methods were proposed to adjust published values of mineral dissolution rates considering different hydrodynamic and chemical conditions (Colombani, 2008, 2016); for calcite, final values within a single order of magnitude were recently obtained (Colombani, 2016). Pignatelli et al. (2016) proposed to relate the average dissolution rate with the the number of topological constraints between the atoms of the bulk solid and showed that this relation could be extended to some glasses; this approach suggests that surface reorganization controls the dissolution rate. Kinetic Monte Carlo (KMC) simulations are very useful for modeling mineral dissolution and may contribute to this debate (Lasaga and Blum, 1986; McCoy and LaFemina, 1997; Lasaga and Luttge, 2001; Zhang and Luttge, 2009; Kurganskaya and Luttge, 2013; Silveira and Reis, 2013; Bandstra and Brantley, 2015; Kurganskaya and Luttge, 2016; Fischer and Luttge, 2017). In KMC, the attachment and detachment rates of single atoms or molecules depend on the interaction energies with their neighborhoods and on temperature (Lasaga and Blum, 1986). The energetic description of a given surface may be obtained from theoretical approaches, such as density functional theory or molecular dynamics simulation, and from experimental methods, such as atomic force microscopy or X-ray reflectivity. For instance, Liang and Baer (1997) observed different reactivity of acute and obtuse steps of calcite surfaces, which motivated the KMC work of McCoy and LaFemina (1997); subsequently, a combination of various ab initio methods provided improved results (de Leeuw et al., 1999; Lardge et al., 2010; Wolthers et al., 2012) 3

98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135

which were recently used in KMC simulations (Kurganskaya and Luttge, 2016). With KMC methods, simulations of systems with sizes . 1µm are usual; these sizes are much larger than those of molecular dynamics simulations, which typically consider ∼ 103 molecules. Despite this advantage, most crystalline grains have sizes 100µm or larger, thus the extrapolation of KMC simulation results may be necessary for the application to macroscopic samples. Other applications of KMC methods show its relevance for relating the energetics of dissolution and the surface morphology, such as modeling dissolution stepwaves in etch pits and the consequent nonlinear variation of dissolution rate near equilibrium (Lasaga and Luttge, 2001), the description of morphological evolution of dissolving feldspar nanoparticles with various shapes (Zhang and Luttge, 2009), and modeling the rate spectrum of calcite (Fischer and Luttge, 2017). KMC simulation may also predict universal features of some classes of systems; for instance, the detachment of non-dissolved clusters observed in a model of dissolution of an inhomogeneous solid (Silveira and Reis, 2013) helped to explain micron-scale grain detachment in limestone weathering (Emmanuel and Levenson, 2014). In the present work, KMC simulations are used to model the dissolution of crystalline solids with initial rough surfaces in conditions of high undersaturation, i. e. in fully irreversible dissolution. A simple cubic lattice (Kossel crystal) is considered, the activation energy is modeled with a bond counting procedure, the detachment rate has an Arrhenius form, and initial surface configurations are produced by a fractional Brownian motion algorithm. The first aim of this study is to answer two questions: (i) What is the difference between the dissolution rate of a rough surface (e. g. a polished surface) and of a flat surface? (ii) How does the dissolution rate evolve in time with these different initial configurations? From the answers to these questions, connections with other problems involving mineral dissolution rates are discussed. In the case of an initial flat surface, we will show that the dissolution rate rapidly reaches a constant value rF . However, for an initial rough surface and small detachment rates (typical of room temperature), the initial rate is much larger than rF and then decreases by some orders of magnitude across several time decades; this is related to the evolution of the surface from higher to lower energy configurations. That time evolution is then treated by a scaling approach, which shows that crossover times and crossover surface retreats are proportional to the initial roughness. For describing the initial dissolution process, we define rough dissolution rates as averages for constant ratios between surface retreat and initial roughness. The rough rates vary with 4

153

temperature but weakly depend on initial morphological features; moreover, they are much larger than rF , which agrees with laboratory results on flat and polished surfaces. Since the dependence of our results with the size of the simulation box is weak and mineral grains frequently have rough boundaries, we suggest that our microscopic interpretation for the decay of the dissolution rate across several time decades may explain the differences in the rates measured in the laboratory and in the environment for several minerals. Calcite dissolution in alkaline environment is considered for application of the model throughout this work, and the results suggest that the narrow range of laboratory rates obtained by Colombani (2016) is an evidence of a universal control by kink site dissolution. The rest of this paper is organized as follows. In Sec. 2, we review basic features of rough surfaces, present our model, and the simulation methods. In Sec. 3, we analyze the dissolution of initial flat surfaces. In Sec. 4, we study the evolution of the dissolution of initial rough surfaces. Sec. 5 presents the scaling approach to rough surface dissolution. Sec. 6 discusses the connection of our results with recent works. Sec. 7 summarizes our results and conclusions.

154

2

155

2.1

136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

156 157 158 159 160 161 162 163 164

165 166 167

Model and basic concepts Description of rough surfaces

The solids have simple cubic structure (Kossel crystals) with sites of edge a. The order of magnitude of the size a is that of the size of one molecule, which is ∼ 1nm for a large variety of minerals. All lengths defined in this work (height, roughness, correlation length, etc) are dimensionless quantities defined in units of a. The definition of dimensionless time t is presented in Sec. 2.3, where the dissolution model is defined. The global roughness of a surface at time t is defined as the root-meansquare (rms) fluctuation of the surface height z measured relative to the reference plane xy: [⟨ ⟩]1/2 w (t) ≡ (z − z)2 , (1) where the overbars indicate a spatial average (along the x and y directions) and the angular brackets indicate an average over different configurations at a given time. Hereafter we refer to w simply as surface roughness; in the 5

168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183

184 185

186

187 188 189 190 191 192 193 194 195 196 197 198 199 200 201

context of kinetic roughening, it is usually called interface width (Barab´asi and Stanley, 1995; Krug, 1997). The autocorrelation function characterizes the lateral correlations of surface fluctuations in a distance s at time t (Zhao et al., 2001): Γ (s, t) ≡ ⟨[˜ z (⃗s0 + ⃗s, t) z˜ (⃗s0 , t)]⟩, where |⃗s| ≡ s, ⃗s has different orientations along the surface, z˜ ≡ z − z, and the averages are taken over different initial positions ⃗s0 in the xy plane and different configurations. For a fixed time, the correlation length ξ (t) is defined as the smallest distance s in which Γ (s, t) = 0. This definition is suitable for surfaces with large hills and valleys, which is the case of our samples. The height fluctuations in lengthscales “s” smaller than the lateral size of the sample are characterized by the local roughness wl . We let a square box of side “s” to glide along the directions of the reference plane and, at each position, we compute the rms height fluctuation; the average value is defined as the local roughness wl (s) for that box size (Zhao et al., 2001). For s ≪ ξ, the local roughness scales as wl ∼ s H , (2) where H is called Hurst exponent or roughness exponent. For large s, we have wl → w.

2.2

Construction of the initial surfaces

Our samples have (dimensionless) length l in the x and y directions. The samples with initial flat surfaces have solid sites at all points with z ≤ 0 and solution at z > 0, and periodic boundary conditions in x and y directions. In the following, positions and lengths are also defined in units of the site edge a. The samples with initial rough surfaces are constructed with chosen values of the Hurst exponent H and with global roughness wI . We use the fractional Brownian motion (fBm) algorithm and the midpoint displacement method to generate those surfaces (Russ, 1994), considering free boundary conditions. The algorithm starts with an square in the plane, with vertices placed at (x1 , y1 ),(x2 , y2 ),(x3 , y3 ), and (x4 , y4 ), and corresponding heights z1 ,z2 ,z3 , and z4 , selected by a Gaussian random variable with mean 0 and variance σ 2 . The value of σ is adjusted to provide the global roughness wI . In the first step of the iterative process (i = 0), the height zC0 of the point at the square centroid (xC0 , yC0 ) is vertically displaced from the average of {z1 , z2 , z3 , z4 } by a 6

202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218

{ } Gaussian random offset d0 , with mean 0 and variance (∆0 )2 = 1/[22H ] σ 2 . In the second step (i = 1), we have four centroids of four squares in the plane with sides’ length corresponding to half of the square of previous step. Each centroid is then vertically displaced from the average coordinate z of the vertices of the square a} Gaussian random offset d1 , with mean 0 { by 2 4H and variance (∆1 ) = 1/[2 ] σ 2 . In the i-th step, the number of centroids will be 22i and the same procedure is repeated for their heights z, which are displaced by a }Gaussian random variable with mean 0 and vari{ 2 2(i+1)H ance (∆i ) = 1/[2 ] σ 2 . It is clear that, for larger i, the displacements are comparatively smaller. In this work, i = 12 has been used. Figures 1(a) and 1(b) show surfaces generated by this fBm algorithm with H = 0.3 and H = 0.7, respectively, and with the same global roughness wI = 100. Fluctuations of long wavelength give the largest contributions to wI . Surfaces with H = 0.3 have local asperity higher than that of surfaces with H = 0.7; this means that the short wavelength oscillations in surfaces with H = 0.3 have larger amplitudes. The correlation lengths of these surfaces are not very different: ξI ≈ 203 (H = 0.3) and ξI ≈ 268 (H = 0.7).

Figure 1: Sections of lateral size 200 × 200 of solids whose surfaces were generated by fBm algorithm with (a) H=0.3 and (b) H=0.7. These sections were extracted from lattices with lateral size l = 512 and global roughness wI = 100. Note the differences in short wavelength fluctuations in (a) and (b). 219 220 221

The values of H were set by the algorithm of construction of the surfaces. However, calculation of accurate values of H from Eq. (2) is difficult due to the finite lattice size; see e. g. Chame and Reis (2004). The same difficulty 7

222

is observed if H is calculated from the height-height correlation function.

223

2.3

224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239

Dissolution model

We model dissolution in contact with a highly unsaturated solution, which is an irreversible process in which atoms or molecules detach from the solid surface and are transported by a solution, but the reverse process (viz. precipitation of atoms/molecules from the solution at the solid surface) does not occur. The solid surface is defined as the set of l2 top sites of all columns (x, y). The dissolution process is represented by a sequential detachment (removal) of surface sites. The rate K of detachment of a given site has the thermally activated form D = ν exp [Eb / (kB T )], where ν is a characteristic frequency (physical dimension T−1 ), Eb is the interaction energy between that site and the rest of the solid, T is the temperature, and kB is the Boltzmann constant. We consider a simple bond-counting rule to represent the dependence of Eb on the local surface morphology: Eb = −nE, where n is the number of nearest neighbor solid sites (coordination number) and E > 0 is the absolute value of the binding energy with each neighbor. A detachment probability per nearest neighbor is defined as ϵ = exp [−E/ (kB T )]

240 241

and is used hereafter as the model parameter (instead of E or T ). The detachment rate D of a site with n NN can then be written as D = νϵn .

242 243 244 245 246 247 248 249 250 251 252

(3)

(4)

Figure 2 shows a surface with a terrace where sites with all values of n are present. Some terms frequently used to refer to those sites are: kink and 2-bonded chain sites (n = 3; only kink shown in Fig. 2); step-edge sites (n = 4); terrace sites (n = 5). We refer to sites with n = 1 as isolated sites and to sites with n = 2 as 1-bonded sites. The solid surface does not have overhangs because only the top site of a column (x, y) can be dissolved at each step. Consequently, aggregates with unreacted sites cannot detach from the rest of the solid. A related model by Silveira and Reis (2013) included the possibility of detachment of non-dissolved aggregates, but showed that this phenomenon is negligible for ϵ . 0.1 (which is the case in this work). Also note that our model does not 8

n=1 n=5

n=3 n=2 n=4

Figure 2: Part of a surface with two terraces and a variety of sites with the indicated coordination numbers. Exposed faces of those sites are in bright colors and covered faces have edges indicated with dashed lines.

254

consider surface relaxation during the dissolution; the only changes in the surface configuration are those due to site detachment.

255

2.4

253

256 257 258 259

Simulation method

The simplest algorithm to simulate the model is the random select and test method, in which random surface sites are sequentially chosen and removed with probability P = ϵn . From Eq. (4), the attempt to detach l2 sites (i. e. one attempt per column) takes place in the time interval τ1 ≡

260 261 262 263 264 265 266 267 268 269 270 271

1 . ν

(5)

However, this algorithm is very time consuming for small ϵ because most attempts are not executed. For this reason, here we adopt an algorithm that keeps track of the probabilities of removal of all surface sites and uses this set to choose a site to be removed with probability 1. For details, the reader may consult Meakin and Rosso (2008), who call it the divide and conquer algorithm. A summary of its basic steps is: 1) Five lists of sites, labeled according to the coordination number n = {1, . . . , 5}, are created; each list contains the coordinates (x, y) of all sites with the corresponding n. All sites in each list have the same dissolution rate. 2) The sum of the dissolution rates of all surface sites, Σ, is calculated (directly from the number of elements in each list and the corresponding rates). 3) One of the lists is randomly chosen; the probability of 9

285

choosing list n is proportional to the number of elements of this list multiplied by ϵn . 4) One of the sites of the chosen list is now selected at random (equal probability to all sites in the list) and this site is removed. 5) The time is incremented by 1/Σ. 6) The column (x, y) of the removed site and the NN columns have their positions in the lists updated. Most simulations consider lateral size l = 512; some samples with l = 256 and l = 1024 were also used to check for possible finite-size effects. All results in initial flat surfaces were obtained in l = 1024. The range 10−4 ≤ ϵ ≤ 10−1 was considered. At room temperature (T = 300K), it corresponds to binding energy E from 0.06eV (6kJ/mol) to 0.24eV (23kJ/mol). Thus, the set of model parameters includes ϵ, the lateral size l, and the initial surface parameters wI and H. For each set, dissolution was simulated in 50-100 different initial configurations with a maximal number of dissolved layers Nmax = 103 .

286

2.5

272 273 274 275 276 277 278 279 280 281 282 283 284

287 288

Calculation of primary quantities

The simplest quantity measured during the dissolution process is the number of detached sites NS . The number of dissolved molecular layers is defined as NS . (6) l2 Thus, detachment of one layer corresponds to the dissolution of an average of one molecule per column. The average length of surface retreat is N a. The number of dissolved moles is n = NS /NA , where NA = 6.02 × 1023 is the Avogadro number. The dissolution process may be interpreted as an attack of the solution to randomly chosen surface sites, which are detached with the probability P defined in Sec. 2.4. The time τ1 [Eq. (5)] is the time interval for attempting to dissolve one molecular layer. If ϵ = 1, all detachment attempts are executed and we obtain ∆N = 1 in τ1 ; for ϵ < 1, ∆N < 1 in this time due to the rejected attempts. The total dissolution time is denoted as τ . In our simulations, the dimensionless time t is measured in units of τ1 : τ t≡ = ντ. (7) τ1 N≡

289 290 291 292 293 294 295 296 297 298 299 300

301 302

Thus, the time t is numerically equal to the number of molecular layers attacked by the solution (but not necessarily dissolved). 10

303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321

322 323 324 325 326 327 328

The value usually adopted for ν is the intermolecular vibrational frequency of bulk water, 1012 s−1 , which leads to τ1 = 10−12 s. This frequency is expected to be a reasonable approximation for the rate of attacks of the water molecules to a silica surface (Pelmenschikov et al., 2001) and is also adopted in simulations of calcite dissolution (Fischer et al., 2014; Kurganskaya and Luttge, 2016). We will consider this value in the quantitative applications of our model. The algorithm used here implies a uniform variation of NS and N after each simulation step. However, the time interval between two consecutive site detachments is 1/Σ (Sec. 2.4), which is not uniform during the simulation; this inteval depends on ϵ and on the numbers of elements in the simulation lists. For this reason, for a given parameter set, we calculate the average time of dissolution for each integer value of N , where this average is taken over the different initial configurations. The average dissolution time is measured in units of τ1 and is also denoted as t. For this reason, in the presentation of our results, the intervals of N are uniform, but the intervals of t are not uniform. The instantaneous dissolution rate is defined as the number of dissolved moles per unit area per unit time (physical dimension L−2 T−1 ): ∆n ∆N R≡ = , (8) 2 NA a2 ∆τ (la) ∆τ where ∆NS is the number of dissolved layers in the time ∆τ . This rate is normalized by the geometric (flat) area of the sample. An alternative definition considers the normalization by the exposed surface area AE , which includes top and lateral faces of the cubic sites. This corresponds to the area measured by adsorption of small molecules, such as BET adsorption. Here, the dimensionless exposed area AE is measured in units of a2 (area of one face of the cubic site). The corresponding dissolution rate is ∆n l2 = R. (9) AE a2 ∆τ AE In this definition, one may consider a constant value of AE measured before dissolution (as usual in experimental works) or consider the instantaneous value of AE , which changes in the course of the dissolution. In our simulations, we measure dimensionless dissolution rates r and k, which are related to R and K as ∆N r≡ = NA a2 νR, (10) ∆t K≡

329 330 331 332 333

11

l2 ∆N = NA l2 a2 AE νK. (11) AE ∆t The rate r may be interpreted as a dimensionless velocity of surface retreat. Other microscopic and macroscopic quantities were measured during the dissolution: the surface roughness w, the lateral correlation length ξ, and the fractions of surface sites with each coordination number n, which are denoted ∑ as fn (n = 1, . . . , 5); normalization gives 5n=1 fn = 1. k≡

334 335 336 337 338

339

3

340

3.1

341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365

Dissolution of initial flat surfaces Long time dissolution

Figs. 3(a)-(b) show the number of dissolved layers as a function of time t for two values of ϵ. The data fit straight lines with good accuracy, which means that the dissolution rate attains an approximately constant value at short times; this is confirmed by estimates of the slopes in various time intervals. Thus, the dissolution rate apparently reaches a steady state after dissolution of a few layers. The average rate rF is calculated considering the retreat of N = 103 layers, i. e. the full time ranges shown in Figs. 3(a)-(b). Since the surface has very low roughness during the dissolution, kF ≈ rF . Using ν = 1012 s−1 , the maximal dissolution time in Fig. 3(a) (ϵ = 10−2 ) is approximately one tenth of second, but in Fig. 3(b) (ϵ = 10−4 ) it is almost one year. This shows the large sensitivity of the dissolution rate with the physical parameters that control ϵ, namely the temperature and the bond energy. The apparent steady state of dissolution can be explained by kinetic roughening concepts. We first recall that Silveira and Reis (2013) simulated this model for ϵ ∼ 0.1 and showed that the roughness has Kardar-ParisiZhang (KPZ) scaling (Kardar et al., 1986). Universality concepts indicate that KPZ scaling is expected for all values of ϵ (Barab´asi and Stanley, 1995; Krug, 1997). As shown by Krug (1990), the velocity of surface retreat in a KPZ process varies as r ≈ r∞ (1 + bN −0.77±0.02 ), where b is a constant of order 1 and r∞ is the asymptotic value of that velocity. For N ∼ 102 -103 , the second term in this expression is small, which explains the small variations of r at long times. This justifies the use of the average rate rF as representative of r∞ ; however, the power law variation of r rules out the formal use of the 12

1000

-2

=10

=10

-4

N

800 600 400 200 0

(b)

(a)

0

4

8 t

12

0

5 10 15 20 25

-10

t

10

-18

10

-2

10

-6

rF

10

-10

10

-14

10

(c) -18

10

-4

10

-3

10

-2

10

-1

10

Figure 3: (a), (b) Number of dissolved molecular layers (dimensionless surface retreat) as a function of time for the values of ϵ indicated and initial flat surfaces. (c) Average dissolution rate rF as a function of ϵ; the dashed line is a least squares fit of the data with slope 4.153 ± 0.002. 366 367 368

369 370 371 372

term “steady state” for this process. The estimates of rF are shown in Fig. 3(c) as a function of ϵ. The linear fit gives rF ≈ r0 ϵα = r0 exp [−αE/ (kB T )], (12) where r0 = 1.86 (2) and α = 4.153 (2). The estimate of exponent α is consistent with (and quantitatively improves) the value 4.0 obtained in Wehrli (1989) for a kinetic Ising model which is equivalent to our model in the limit of small ϵ.

13

399

These results partly answer question (ii) proposed in Sec. 1: the dissolution rate of an initially rough surface rapidly attains an approximately constant value (informally termed steady state), with temperature dependence given by Eq. (12). The physical interpretation of Eq. (12) is that the thermally activated detachment at molecular scale [Eqs. (3) and (4)] implies the thermally activated (act) dissolution at macroscopic scale, with an activation energy EF = αE. The value α ≈ 4 indicates that step-edge dissolution is the rate limiting step of the macroscopic process. This occurs because large terraces with long steps appear in the solid surface, with a very small number of kink sites or other low energy sites. This interpretation is also consistent with that of Lasaga and Luttge (2004) for a similar model, which showed that the activation energy is related to the average number of bonds broken during dissolution. At room temperature and for the range 10−4 ≤ ϵ ≤ 10−1 simulated here, we (act) have 24kJ/mol≤ EF ≤ 95kJ/mol. We apply these results to the alkaline dissolution of calcite, whose rate is weakly dependent on pH (Arvidson et al., 2003). In a recent work, Fischer and Luttge (2017) showed dissolution rate spectra of calcite in which the minimum measured value was ≈ 3 × 10−9 mol/(m2 s) (the spectra spanned more than two orders of magnitude of the rate). This minimum is expected to correspond to dissolution of the local surface configurations with the lowest energy (in the range where the equipment can measure). Here we assume that this minimum is an approximation of the flat surface dissolution rate RF . Using Eq. (10) with a = 0.35nm, we obtain rF ≈ 2 × 10−16 , and substitution in Eq. (12) gives ϵ ≈ 1.5 × 10−4 and E/kB T ≈ 8.8. At room temperature, this energy is near the lowest bound obtained by Raiteri et al. (2010) using force field methods.

400

3.2

373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398

401 402 403 404 405 406 407 408

Short time dissolution

Figures 4(a) and 4(b) show the evolution of the rate r and of the surface roughness w, respectively, for ϵ = 10−3 at very short times. In this case, the values of these quantities were calculated in very short time intervals, which correspond to variations ∆N = 0.002 in the number of dissolved layers. The results are shown for a single realization and for an average among 20 realizations. Each oscillation in Figs. 4(a),(b) corresponds to dissolution of one layer. The largest values of r and w are obtained when there is approximately half 14

409 410 411

coverage of the top layer; the smallest values are obtained when the surface has a very large terrace with a few fluctuations. Similar oscillations in r were recently shown in Fischer et al. (2014).

r/10

-12

1 sample

2

20 samples

(a)

1

0

0.6

(b)

0.3

0.0 0.0

1.2

2.4

t

10

3.6

-13

Figure 4: Short time evolution of (a) the instantaneous dissolution rate and (b) the surface roughness in one sample (red) and in an average over 20 samples (blue); simulations are for ϵ = 10−3 .

419

An important feature in Figs. 4(a),(b) is that the oscillations for one and for 20 realizations are in phase when the number of dissolved layers is 10 or smaller. Thus, a single realization contains a large number of microscopic environments, which are the same environments represented in many realizations. As time increases, deviations appear, which means that the number of such environments increased and cannot be represented by a single configuration; in this case, averages over several configurations are necessary to obtain reliable estimates for a given set of parameters.

420

4

421

4.1

412 413 414 415 416 417 418

422 423 424

Dissolution of rough surfaces Evolution of the dissolution rate

In Figs. 5(a)-(d), we show the number of dissolved layers as a function of time for two values of ϵ and several initial configurations. The results for flat surfaces with the same value of ϵ are also shown. 15

3

10

H=0.3

H=0.3 I

2

10

I

= 50

=10

N

I

= 200

I

-2

= 50 = 200

=10

-4

1

10

0

10

(a)

2

10 3

10

6

10

10

t

(b)

3

10

2

10

N

I

=10

11

15

19

10 10 10 10 10 t H=0.7

H=0.7 I

7

= 50

I

= 200

I

-2

=10

= 50 = 200 -4

1

10

0

10

(c)

2

10

6

10

10

t

(d)

3

10

7

11

15

19

10 10 10 10 10 t

Figure 5: Number of dissolved layers as a function of time for different values of Hurst exponent, initial roughness, and detachment probability. All data were measured in l = 512. Dashed lines are results for flat surfaces with the same value of ϵ of each plot and l = 1024. 425 426 427 428 429 430 431 432 433 434 435

When the number of dissolved layers N is not very large (typically N . 100), Figures 5(a)-(d) show that the dissolution of a flat surface occurs in a much longer time than in the rough surfaces; this means that rough surface dissolution is much faster than flat surface dissolution. This effect is enhanced as ϵ decreases. For instance, consider a retreat of N ∼ 100 layers: compared to the flat surface, the dissolution times are smaller by a factor up to ∼ 102 for ϵ = 10−2 and up to ∼ 104 for ϵ = 10−4 . Figures 5(a)-(d) also show that the larger the initial roughness wI , the faster the initial dissolution. However, for the largest values of N , dissolution times of initial flat and initial rough surfaces are not very different. This suggests significant smoothing of the initial rough configurations. 16

436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473

Figure 6(a) shows the time evolution of the (instantaneous) dissolution rate for wI = 100 and two values of H, with ϵ = 10−4 ; the corresponding value for initial flat surfaces, rF , is also indicated. Figures 6(b)-(f) show the evolution of several quantities during the dissolution: fractions of sites with each value of coordination number n, exposed surface area, roughness, and lateral correlation length. These data were obtained in lattices with l = 512. Results in l = 256 for the same parameters are presented in the electronic supplement and have small quantitative differences. The dissolution of the surfaces with H = 0.3 begins much earlier than the dissolution of surfaces with H = 0.7. This is explained by the high density of small peaks (needle-like structures) in the surfaces with H = 0.3 [Fig. 1(a)], which facilitates detachment. However, after the dissolution of a few tens of layers, the values of r for different H differ by less than one order of magnitude; this typically corresponds to t > t1 in Fig. 6(a). The most striking result in Fig. 6(a) is the decay of r in time: it decreases four orders of magnitude when the number of dissolved layers increases from ∼ 10 to ∼ 103 , for both H, and this decay takes place in six orders of magnitude of dissolution time. For ν = 1012 s−1 , the time range is from ≈ t1 = 2.4s to 1 year. For calcite (a ≈ 0.35nm, ν ≈ 1012 s−1 ), the dissolution rate decreases from ∼ 10−5 mol/(m2 s) to ∼ 10−9 mol/(m2 s). Note that these values of dissolution rates and these ranges of dissolution times are strongly dependent on ϵ and wI . The variations in the fractions fn , shown in Fig. 6(b) (H = 0.3) and Fig. 6(c) (H = 0.7), can be related to the variations in r. As time increases, the surface is dominated by sites with increasing coordination number. Indeed, sites with small coordination are rapidly dissolved, then leaving sites with larger coordination exposed at the surface. For small ϵ, the ratio between the probability of dissolution of a site with coordination n and another site with coordination n + 1 is very large. As time increases, this explains the increase in the time interval necessary for detachment of a single layer. The ratio between exposed area and geometrical area, AE /l2 , is shown in Fig. 6(d). For H = 0.3, it decreases by a factor ∼ 20 after dissolution of 103 layers. For H = 0.7, a decrease by a factor ∼ 3 is observed. These factors are very small compared to the variations of the dissolution rate, which reach factors ∼ 104 in the same conditions [Fig. 6(a)]. They also imply that the dissolution rate k [Eq. (11)], which is normalized by the exposed area AE , varies of several orders of magnitude in the course of the dissolution. This is confirmed by simulation results presented in the electronic supplement. 17

-12

2

t

(d)

3

t

10

4

t

-15

5

10

1.00

n

f

n

(b)

0.75

(e)

1

3

35

4

0.25

5

0.00 1.00

f

n

0.75

(c)

(f)

n

1 2

0.50

4

7

10

10

13

16

10

10

19

10

t

7

10

10

10

13

10

16

10

19

10

t

Figure 6: Time evolution of quantities measured in dissolution of initial rough surfaces with wI = 100, ϵ = 10−4 , two values of H [H = 0.3 (green squares); H = 0.7 (red triangles)], and l = 512. (a) Dissolution rate; the rate rF in the flat surface in shown as a blue line. (b), (c) Fraction fn for each coordination number in surfaces with H = 0.3 and H = 0.7, respectively. (d) Ratio between exposed surface area and geometrical area l2 . (e) Global roughness. (f) Correlation length.

18

150

0

5

0.00

10

0 225

75

3

0.25

105 70

2

0.50

2

(a)

E

r

t

1

/ l

-9

10

25 20 15 10 5 0

A

t

474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495

496

497 498 499 500 501 502 503 504 505 506 507 508 509

Roughness and correlation length are the most frequently used quantities to characterize the changes in the global morphology of a surface (Krug, 1997). However, Figs. 6(a) and (e) show that w is reduced by a factor smaller than 10 (relatively to wI ) while r has decreased more than 3 orders of magnitude. The variation of the correlation length [Fig. 6(f)] is much slower; it decreases by a factor ≈ 2 in the same conditions. We conclude that variations of macroscopic quantities (e. g. roughness, exposed area, and correlation length) cannot be used to explain the variations of dissolution rate. Our results also show that the instantaneous dissolution rate has a large variation in any time range that exceeds one decade. Thus, dissolution of rough surfaces cannot be described by an average rate equivalent to rF . Moreover, r ≫ rF in a long time range. These results do not change if we consider the rate k normalized by the exposed area AE because this area has much slower variation [Fig. 6(d)]; see electronic supplement. Qualitatively, questions (i) and (ii) proposed in Sec. 1 were answered here. The dissolution rates obtained in different lattices sizes are similar to the ones presented here for l = 512. This can be confirmed by comparison with the plots in the electronic supplement. Due to the free boundary conditions of the rough surfaces, the fractions of sites with n = 3 and n = 4 may have significant differences in l = 256 and l = 512 in certain time ranges, but r differs by a factor ≈ 3 or smaller.

4.2

Surface morphology during the dissolution

Figures 7(a)-(f) show snapshots of a dissolving surface with wI = 100, H = 0.3, and ϵ = 10−4 . The time t0 corresponds to the initial (not dissolved) surface and the times {t1 , . . . , t5 } are indicated in Fig. 6(a). The first plateau of r appears at t1 ∼ 1012 , when isolated sites and 1bonded sites are completely removed and the fraction of kink sites has a peak [Figs. 6(a),(b)]. The surface at t1 is visually not very different from its initial configuration [Figs. 7(a),(b)]. At t2 ∼ 1014 , close inspection of Fig. 7(c) reveals the increase in the number of flat regions with disordered boundaries at hills and valleys. At this time, the fraction of kink sites has decreased and the fractions of step-edge sites and terrace sites have increased; this is the beginning of a regime of terrace growth. At time t3 , approximately 100 layers have been dissolved, i. e. the surface retreat is numerically equal to the initial roughness. The rate r is decreasing 19

Figure 7: Snapshots of a dissolving surface with wI = 100, H = 0.3, l = 512 and ϵ = 10−4 , at the times (a) t0 = 0 (initial condition) , (b) t1 = 2.39 × 1012 , (c) t2 = 8.42 × 1013 , (d) t3 = 3.78 × 1014 , (e) t4 = 1.05 × 1016 , and t5 = 1.20 × 1018 , as indicated in Fig. 6(a). The color codes refer to the local surface height; observe the different code in (f). 510 511 512

[Fig. 6(a)], the fraction of kink sites is small, f4 fluctuates and f5 continues to increase [Fig. 6(b)]. Comparison of Figs. 7(c) and 7(d) shows that this is a process of terrace broadening at the lowest parts of the surface; the 20

530

highest parts have smaller terraces. We still observe large height differences between the hills and valleys in Fig. 7(d); indeed, the roughness at t3 is large (w ≈ 70%wI ) and the correlation length has a small change in comparison with the initial value [Figs. 6(e),(f)] . This means that only short wavelength fluctuations were suppressed during the dissolution until ∼ t3 ; long wavelength fluctuations remain and give the main contributions to w and ξ. At time t4 ∼ 1016 , another plateau of the dissolution rate is attained [Fig. 6(a)]. The surface shown in Fig. 7(e) is smooth due to terrace coarsening at the lowest parts. Only step-edge and terrace sites have significant fractions [Fig. 6(b)], and the roughness is much smaller than the initial value [Fig. 6(e)]. Some peaks remain, which are reminiscent of the initial large fluctuations. The correlation length did not have a significant decrease because it is an estimate of the distance between the largest surface fluctuations. Finally, at t5 ∼ 1018 , Fig. 7(f) shows a surface with large terraces separated by long step edges, similar to the ones obtained in dissolution of an initially flat surface. The fraction of step-edge sites is much smaller than that of terrace sites [Fig. 6(b)] and the dissolution rate has reached a value close to rF [Fig. 6(a)].

531

5

513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529

532

533

534 535 536 537 538 539 540 541 542 543 544

5.1

Scaling approach to rough surface dissolution Crossover to low roughness configurations

The results of Sec. 4 show the crossover from a regime with large roughness and large rate to a regime with relatively smooth surfaces and small rate. This is a type of smoothing process, which is expected to obey scaling laws predicted by kinetic roughening concepts, as formerly done in the context of film deposition in rough substrates (Majaniemi et al., 1996; Nguyen et al., 2010; de Assis and Reis, 2015). Our first step to describe the crossover is to analyze the variation of r/rF with the number of dissolved layers N , as shown in Fig. 8 for initial surfaces with H = 0.3, several values of wI , and two values of ϵ. For a given set of parameters, we define the crossover time tc and the crossover surface retreat Nc as the values of these quantities when rc = 1.5rF , 21

(13)

545 546 547 548 549

i. e. when the rate is only 50% larger than that in a flat surface. When this condition is attained, the surface has large terraces and a small global roughness, similarly to that shown in Fig. 7(e). The crossover points are indicated in Fig. 8 for each wI : Nc is of the same order as wI and is weakly dependent on ϵ. Results for H = 0.7 also show that ϵ weakly affects Nc . 2

10

25

(a)

50

1

100

10

r/r

F

200

0

10

8

10

(b)

5

10

2

10

-1

10

0

200

400

600

N

800

1000

Figure 8: Ratio r/rF as a function of N for initial surfaces with H = 0.3, several values of wI , (a) ϵ = 10−2 and (b) ϵ = 10−4 . In all cases, the lattice size is l = 512. The crossover points are indicated by arrows for each data set. 550 551 552 553 554 555

556 557 558 559 560 561 562

Figures 9(a) and 9(b) show Nc as a function of the initial roughness wI for H = 0.3 and H = 0.7, respectively, with several values of ϵ. Most results are for l = 512, but some results for l = 256 are also shown and confirm the small effect of the lattice size. For small ϵ, the data collapse in a single line shows that Nc is proportional to wI (deviations are observed only for ϵ = 10−1 ). Linear fits give Nc ≈ AwI , (14) with A ≈ 3.6 for H = 0.3 and A ≈ 2.8 for H = 0.7. Eq. (14) does not depend on temperature nor on the binding energy. Moreover, for different roughness patterns (i. e. different H), the relation is valid with approximately the same amplitude A ∼ 3. In Fig. 10, we present a scheme that illustrates the consequences of this result; in summary, it suggests that the ration N/wI controls the morphological evolution of the dissolving surface. After a surface retreat N ≈ wI , terraces are formed, 22

-4

=10

-4

=10

1000 (l=256)

-3

800

-2

600

Nc

=10 =10

-1

=10

(a)

1000

200 0

800

-4

=10

-4

=10

(l=256)

-3

600

=10

400

=10

-2

-1

=10

200 0

400

(b)

0

40

80

120

160

200

I

Figure 9: Dimensionless retreat at the crossover, Nc , as a function of the initial roughness wI , for several values of ϵ and Hurst exponents (a) H = 0.3 and (b) H = 0.7. The lattice size is l = 512 except where indicated. Dashed lines are linear fits of the data for ϵ ≤ 10−2 . 563 564 565 566 567 568 569 570 571 572 573 574 575 576

577 578 579 580

but are not very wide; the roughness is still large, in agreement with the discussion in Sec. 4.2. However, when the surface retreat is N ≈ 3wI , the process is in the crossover regime and the surface is almost flat. For this reason, after dissolution of a thickness ∼ 3 times larger than the initial roughness, the dissolution rate becomes very close to that of the surface with the lowest energy. Note that this result does not depend on the time necessary to reach the configurations with N ≈ wI or with N ≈ 3wI . Now we analyze the crossover time tc . For fixed ϵ, tc is expected to be proportional to Nc and, consequently, to be proportional to wI . Thus, in Figs. 11(a) and 11(b) we plot tc /wI as a function of ϵ for H = 0.3 and H = 0.7, respectively, with several values of wI . The plots include some data from a smaller lattice (l = 256), which shows that finite-size effects are small. The linear fits in Figs. 11(a),(b) give a universal relation for the crossover time as tc ∼ wI ϵ−3.9 . (15) The exponent in Eq. (15) is close to −α, which is the exponent of rF in Eq. (12). This expresses the fact that the crossover regime is also controlled by dissolution of step-edge sites (n = 4), which is the longest regime necessary to reach smooth surface configurations. 23

WI (a)

WI

L c =3WI

(b)

2WI

Figure 10: Schematic view of three stages of dissolution: an initial rough surface, a configuration when the surface retreat is equal to the initial roughness, and a configuration at the crossover time. 17

10 (a) =200

I

I

9

t /

=100

10

c

I

=50

=200 (l=256)

I

5

10

I

=25

=100 (l=256)

I

17

13

10

I

1

10

10

13

(b)

10

=200

I

9

10

=100

I

5

10

I

=100 (l=256)

=25

1

10

=200 (l=256)

=50

I

I

I

-4

10

-3

-2

10

10

-1

10

Figure 11: Scaled crossover time as a function of the detachment probability for several values of wI and Hurst exponents (a) H = 0.3 and (b) H = 0.7. The lattice size is l = 512 except where indicated. Dashed lines are least squares fits of all plotted data in each panel. 581

582 583 584 585 586

5.2

The rough dissolution rates

Sec. 5.1 showed universal relations for Nc and tc , which are measured in conditions where the surface retreat is proportional to the initial roughness wI . Here we extend this analysis to regimes much earlier than the crossover times, in which the surface retreat N is smaller than or near wI . As shown in Fig. 10, this is a regime in which the surface roughness does not have a

24

587 588 589 590

significant decrease in comparison with wI . We begin with the definition of two average dissolution rates. The first one is the average rate of dissolution of a thickness wI . Letting tW be the average time of this process, i. e. N (tW ) = wI , this rate is wI . tW

rW ≡ 591 592

The second rate is the average for dissolution of a thickness wI /3 and tW 3 is the corresponding dissolution time, i. e. N (tW 3 ) = wI /3: rW 3 ≡

593 594 595 596 597 598 599 600

602 603 604 605 606 607 608 609 610 611 612 613 614 615

wI . 3tW 3

(17)

Hereafter we refer to rW and rW 3 as rough rates. Figures 12(a)-(d) show the rough rates as a function of ϵ for initial surfaces with the two values of H and several values of wI . We observe an excellent data collapse for different values of wI ; some data for l = 256 are also shown to confirm small finite-size effects. The rates for H = 0.3 are larger than those for H = 0.7 due to the larger number of low coordination sites in the initial samples. The linear fits in Figs. 12(a),(b) give the scaling of the rough rates as rW ≈ r1 ϵβ

601

(16)

,

rW 3 ≈ r3 ϵγ .

(18)

For H = 0.3, we obtained β = γ = 3.0, r1 = 0.8, and r3 = 1.1; for H = 0.7, we obtained β = 3.4, γ = 3.1, r1 = 1.3, and r3 = 0.6. Note that the prefactors r1 and r3 are not very different. The only significant differences in the rough rates at a given temperature are due to the small differences in the exponents β and γ. In the electronic supplement, we show the rough rates kW and kW 3 , which are defined for the same time intervals of rW and rW 3 but are normalized by the initial exposed surface area AE , analogously to Eq. (11). These rates also show a good data collapse for different wI and have the same power-law dependence on ϵ of Eq. (18); the corresponding values of γ, β, r1 , and r3 are close to the above estimates for rW and rW 3 . Eq. (18) shows that the thermally activated detachment at molecular scale also implies thermal activation of the rough rates, with activation en(act) (act) ergies EW = βE and EW 3 = γE. The values of β are slightly different for the different initial patterns (i. e. different H). For H = 0.3, it indicates 25

(a)

-3

10

H=0.3

-7

r

W

10 =50

=200

I

I

=100

=25

I

-3

10

I

=200 (l=256)

I

-11

10

=100 (l=256)

I

-15

10 (b)

H=0.7

-7

10

-11

10

=50

=200

I

10

I

=100

-15

=25

I

-4

10

(c)

-3

I

=200 (l=256)

I

=100 (l=256)

I

-2

10

-1

10

10

-3

H=0.3

10

-7

r

W3

10 =50

=200

I

I

=100

=25

I

=200 (l=256)

I

=100 (l=256)

-11

10

I

-15

10

-3

10

I

(d)

H=0.7

-7

10

-11

10 10

=50

=200

I

I

=100

=25

I

-15 -4

10

-3

I

-2

10

10

=200 (l=256)

I

=100 (l=256)

I

-1

10

Figure 12: Rough rates as a function of ϵ: (a), (b) rW for H = 0.3 and H = 0.7, respectively; (c), (d) rW 3 for H = 0.3 and H = 0.7, respectively. The lattice size is l = 512 except where indicated. The dashed lines are least squares fits of each data set. 616 617 618 619 620 621 622 623 624

control by kink site dissolution, which is characteristic of the process in which terraces are nucleated and begin to grow. For H = 0.7, it indicates a transition from the kink site to the step-edge site control; this is characteristic of the growth and coarsening of terraces, in which the fraction of kink sites is decreasing. On the other hand, the γ ≈ 3 suggests that the control by kink site dissolution is universal in the regime where N ≈ wI /3. The variation of the rough rates with ϵ [Eq. 18] is much slower than the variation of rF [Eq. 12] because β < α and γ < α. The ratio rF /rW 3 is typically of order ϵ; at room temperature, this means that rough rates are

26

650

typically much larger than dissolution rates of cleaved surfaces of the same material, in which the number of defects is small. These results quantitatively answer question (i) proposed in Sec. 1. In the present model, the kink sites are the weakest bonds present during the nucleation and initial growth of terraces with disordered boundaries, which are the surface configurations with the lowest energy. This interpretation may be extended to dissolution of minerals with other lattice structures. In the nucleation and initial growth of the most stable surface configurations of a given lattice, there are some molecules with high energy, typically at the boundaries of those configurations. The activation energy of a rough rate is expected to be related to the detachment of those molecules. This energy can then be related to the microscopic bond energy by the factors β and γ, which also depend on the lattice. As a quantitative application, consider that a rough calcite surface is modeled by a fBm surface with H = 0.7, which has small local asperity. Our application of the flat surface dissolution model in Sec. 3.1 gave an estimate ϵ ≈ 1.5 × 10−4 ; substitution in Eq. (18) yields rW ≈ 1.3 × 10−13 and rW 3 ≈ 8.4 × 10−13 . Using the estimates of ν and a, the dissolution rates are RW ≈ 2 × 10−6 mol/(m2 s) and RW 3 ≈ 1 × 10−5 mol/(m2 s). These values are normalized by the flat area of the sample. Instead, normalization by the initial exposed surface area gives KW ≈ 0.7 × 10−6 mol/(m2 s) and KW 3 ≈ 2 × 10−6 mol/(m2 s). The values of RW , RW 3 , KW , and KW 3 are very close to the laboratory estimates for calcite dissolution in alkaline solution (Arvidson et al., 2003; Colombani, 2016). The consequences of this agreement are discussed in Sec. 6.3.

651

6

652

6.1

625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649

653 654 655 656 657 658

Discussion Dissolution of minerals with rough surfaces

Fischer et al. (2012) studied dissolution of cleaved calcite surfaces with initial roughness ∼ 1nm and of polished surfaces with global roughness in the range 120nm-180nm. At pH 4.7, they showed that the dissolution rate of the polished samples was larger than the rate of the cleaved surface by a factor ∼ 7, even with the presence of some etch pits in the latter, which tend to fasten dissolution. Our results are qualitatively consistent with these observations. 27

682

The largest surface retreats of the polished samples were 3 − 4 times larger than their initial roughness, but the change in the global roughness was not significant during this process [compare e. g. the roughness in the largest observation windows shown in Figs. 3B and 3C of Fischer et al. (2012)]. Thus, the scheme of Fig. 10 is not applicable, i. e. the roughness is still large after a retreat ∆N ∼ 3wI ; for this reason, a quantitative comparison with the rates of our model (e. g. rough rates) is not possible. Levenson and Emmanuel (2013) studied dissolution of a limestone (98% calcite, 2% quartz) with micron-sized pores and showed that dissolution in those pores was slower than in the polished surface surrounding them. They proposed that the faster dissolution in the polished region was due to a larger density of high curvature features. Our results support this interpretation. In particular, note that Figs. 7(a)-(e) clearly show that the dissolution of the central valley is much slower than the dissolution at the borders because the valley rapidly evolves to an approximately flat configuration; this evolution may be compared to Figs. 3, 4, and 6 of Levenson and Emmanuel (2013). However, note that there are cases in which rougher configurations are more stable than flat configurations. For instance, fluorite dissolution shows increase of roughness during dissolution (Godinho et al., 2014), and the lattice model of silicon etching of Mirabella et al. (2014) has bond counting rules that privilege formation of rough patterns. In these cases, dissolution of rough surfaces is certainly slower than dissolution of flat surfaces, but the association of high (low) rate with large (small) surface energy is still applicable.

683

6.2

659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681

684 685 686 687 688 689 690 691 692 693 694

Time decrease of dissolution rates

A compilation of laboratory and field data by White and Brantley (2003) showed that dissolution rates of silicate minerals (e. g. plagioclase, Kfeldspar, biotite, and hornblende) measured in the laboratory are larger than the rates measured in the field by 3-4 orders of magnitude. For instance, the laboratory rate of plagioclase dissolution is 10−12 -10−11 moles/(m2 s) and the field rate is 10−16 -10−15 moles/(m2 s). However, the laboratory works typically last 0.1-1yr and the environmental dissolution lasts 105 -107 yr, i. e. they are separated by 6-7 time decades. Based on these data, a fit as R ∼ t−0.61 was proposed. White and Brantley (2003) also performed long term (6yr) column experiments with previously weathered (kaolinized) granite and obtained dissolution rates consistent with this fit if the time t was taken as the age of the 28

695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732

samples in the regolith. They suggested some factors to explain the small dissolution rate in the environment: changes in mineral surface area, progressive depletion of energetically reactive surfaces, accumulation of leached layers or secondary precipitates, and the conditions close to equilibrium for the weathering reactions. Recent results of long term (13.8 years) flow-through column experiments with granitoids also help to fill the gap between the laboratory and environmental time scales (White et al., 2017). In these experiments, plagioclase dissolution reached an apparent steady state in most samples in the last decade of the work, with a rate ∼ 10−13 moles/(m2 s). This value is one order of magnitude lower than those reported in short duration experiments, but is still 2-3 orders of magnitude higher than the field values. Less than 1% of the total mineral mass was lost in this experiment. Simulations of a reactive transport model (CrunchFlow) by White et al. (2017) show that laboratory and environmental results are consistent only if the exposure of mineral surfaces is significantly limited in the environment compared to column experiments. Our model shows that the apparent laboratory-field discrepancy can be explained by other factor proposed by White and Brantley (2003): the progressive depletion of energetically reactive surfaces, as will be justified below. First, note that unweathered mineral grains with disordered shapes are shown in several works. Some examples are: in Fig. 7 of Hausrath et al. (2011), backscattered electron (BE) images show plagioclase, pyroxene, and iron minerals as disordered domains separated by very rough boundaries; the scanning electron microscopy (SEM) images of Fig. 2 of Bazilevskaya et al. (2013) show plagioclase grains with rough boundaries in granite; in Figs. 2 and 3 of da Silva et al. (2016), macro-photographs (with sizes near 1mm) of granites show grains of quartz, muscovite, biotite, and other minerals whose boundaries are also disordered. The main difference from our samples is the size, which is . 1µm in the simulations, while real crystalline grains usually have sizes 100µm-1mm. However, the rough rates calculated here, which are averages over long periods of time, have approximately the same value in two different lateral sizes, l = 256 and l = 512; moreover, at long times, the rate r measured in those sizes converges to the flat rate rF , which was measured in lattices with l = 1024 and with periodic boundaries. Thus, there is no significant finite-size effect in our simulation data; this suggests that the main features and the physical interpretation may be extended to larger sizes. We believe that the long-term decrease of the dissolution rate discussed 29

754

in White and Brantley (2003) and White et al. (2017) is, at least partially, a consequence of long-term smoothing of mineral grain surfaces which initially had rough surfaces. Grains with disordered shapes initially have a high density of microscopic configurations with high energy, similarly to the kink sites in the Kossel crystal. This is the condition of pristine samples used in the column experiments. However, in geological time scales, the grain surfaces evolve to configurations of lower energy and the dissolution rates decreases. The dissolution of low energy surfaces is the slow step that controls the average dissolution rate in the environment. In the context of our model, the field rate is closer to rF , while the laboratory rate is closer to a rough rate, which differ by a factor ∼ ϵ ≪ 1 at room temperature. For instance, our data in Fig. 6(a) for H = 0.7 (initial surface with small local asperity) shows dissolution of ∼ 103 layers taking place in 7 time decades, with the dissolution rate decreasing by approximately 4 orders of magnitude; these factors are surprisingly close to those shown in White and Brantley (2003) for plagioclase dissolution. This interpretation does not exclude other mechanisms that delay the environmental dissolution. However, several of these mechanisms were already considered in reactive transport models such as the model of White et al. (2017). It would be interesting if these models could also account for the time-dependent dissolution rates associated with the evolution of mineral surfaces from high to low energy configurations.

755

6.3

733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753

756 757 758 759 760 761 762 763 764 765 766 767 768

Calcite dissolution rates

In Secs. 3.1 and 5.2, we obtained calcite dissolution rates with the following procedure: using the minimum rate RF ∼ 10−9 mol/(m2 s) of the dissolution rate spectra in Fischer and Luttge (2017), the fits of rF (flat initial surfaces, l = 1024), and reasonable values of ν and a, we estimated the parameter ϵ; using the fits of rough rates (obtained in l = 256 and 512) and those values of ϵ, ν, and a, we estimated the rough rates RW ≈ 2 × 10−6 mol/(m2 s) and RW 3 ≈ 1×10−5 mol/(m2 s). The corresponding rates normalized by the initial exposed surface area, KW and KW 3 , are in the range 0.7-2 × 10−6 mol/(m2 s). This summary highlights the combination of results in various lattice sizes. Colombani (2016) recently performed an analysis of data for calcite dissolution in alkaline environment and showed that suitable rescaling of published values provides estimates between 0.5×10−6 mol/(m2 s) and 6×10−6 mol/(m2 s) in the limit of highly unsaturated solution. This is a narrow range if com30

786

pared with the spread of estimates found in the literature. Colombani (2016) also proposes that the minimal value corresponds to flat, low energy topography, and that the maximal value corresponds to polished (more energetic) surfaces. The rough rates calculated here are of the same order of magnitude of the rates obtained by Colombani (2016). However, our estimate of RF is 3-4 orders smaller than the rough rates and at least 2 orders of magnitude smaller than the minimal value of Colombani (2016). Thus, we understand that the minimal rate of Colombani (2016) is not representative of low energy surface dissolution. Instead, we interpret the range obtained in that work as representative of a variety of high energy configurations of the initial calcite surfaces; this is a reasonable hypothesis for a variety of grain shapes and sizes. Both the lowest and the highest values in that range are consistent with the rough rates obtained here in a regime dominated by kink site dissolution. Within this interpretation, the data treatment of Colombani (2016) shows the universality of the kink site detachment control in laboratory dissolution (in relatively short time scales); such universality was hidden in the previous spread of literature values.

787

7

769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785

788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803

Conclusion

We studied dissolution of a mineral with a simple cubic lattice structure and initial rough surfaces using a model in which the probability of detachment of a site has Arrhenius form and the energy is a sum over the nearest neighbors. First we analyzed dissolution of initial flat surfaces, in which the dissolution rate rF reaches an approximately constant value. The temperature dependence of rF shows that the process is dominated by dissolution of step edge sites. For initial rough surfaces, during dissolution of 103 molecular layers, no apparent steady state is observed; instead, the dissolution rate varies by several orders of magnitude while the time also varies by several orders. The evolution of the dissolving surfaces show configurations of decreasing energy, beginning with dissolution of isolated sites, then passing to nucleation of terraces at the lowest regions, to the growth of terraces and finally to flattening of the whole sample; in this final regime, the dissolution rate is close to rF . This evolution is qualitatively similar to that of limestone surfaces with pores of micrometer size, in which the pore bottom is flatter and dissolves slower (Levenson and Emmanuel, 2013). The results are also consistent with recent 31

830

observations of faster dissolution in polished calcite surfaces when compared to that of cleaved surfaces (Fischer et al., 2012). Moreover, the time evolution of dissolution rates observed here suggests a microscopic explanation to the decrease of dissolution rates across several time scales (from laboratory scale to geological scale) reported in White and Brantley (2003) and White et al. (2017); this explanation assumes that dissolving mineral grains initially have rough boundaries, which are flattened in geological times. Similarly to other smoothing processes, the crossover from rough to flat surface dissolution observed in our model is characterized by a crossover time tc . We showed that tc and the crossover surface retreat Nc are proportional to the initial roughness. The variation of tc with the temperature is consistent with dominant step-edge dissolution (as is rF ), but Nc is weakly dependent on temperature. These relations are weakly dependent on the lattice size, which suggests their use to describe dissolution of macroscopic samples, at least qualitatively. In order to describe short term dissolution, we defined rough dissolution rates as averages for surface retreats proportional to wI . They were calculated for proportionality factors 1 and 1/3 and showed temperature dependence controlled by kink site detachment, which is characteristic of the nucleation and formation of terraces. In systems where dissolving crystalline grains have irregular shapes but narrow size distributions, these results suggest a narrow range of the order of magnitude of the dissolution rates in highly unsaturated solution. This is a possible explanation for the narrow range of calcite dissolution rates obtained by Colombani (2016), which reflects the universal control by kink site dissolution in laboratory works. We expect that the combination of numerical simulations and scaling concepts presented here can also be useful for future studies of mineral dissolution, possibly with models for specific materials.

831

Acknowledgment

804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829

832 833 834

This work was supported by CNPq, FAPERJ, and CAPES (Brazilian agencies). FDAAR thanks Sue Brantley, Diogo Bolster, and Jean Colombani for helpful discussions.

32

835

Appendix

837

Supplementary data associated with this article can be found in a separate document.

838

References

836

839 840 841

842 843 844

845 846

847 848 849 850

851 852

853 854

855 856

857 858 859 860

861 862

R. S. Arvidson, I. E. Ertan, J. E. Amonette, and A. Luttge. Variation in calcite dissolution rates: A fundamental problem? Geochim. Cosmochim. Acta, 67:1623–1634, 2003. J. Z. Bandstra and S. L. Brantley. Understanding the mechanisms of solidwater reactions through analysis of surface topography. Phys. Rev. E, 92: 062114, 2015. A.-L. Barab´asi and H. E. Stanley. Fractal Concepts in Surface Growth. Cambridge University Press, New York, USA, 1995. E. Bazilevskaya, M. Lebedeva, M. Pavich, G. Rother, D. Y. Parkinson, D. Cole, and S. L. Brantley. Where fast weathering creates thin regolith and slow weathering creates thick regolith. Earth Surf. Process. Landforms, 38:847–858, 2013. A. Chame and F. D. A. A. Reis. Scaling of local interface width of statistical growth models. Surf. Sci., 553:145–154, 2004. J. Colombani. Measurement of the pure dissolution rate constant of a mineral in water. Geochim. Cosmochim. Acta, 72:5634–5640, 2008. J. Colombani. The alkaline dissolution rate of calcite. J. Phys. Chem. C Lett., 7:2376–2380, 2016. Y. J. A. B. da Silva, C. W. A. do Nascimento, C. M. Biondi, P. van Straaten, V. S. de Souza Jr, and T. O. Ferreira. Weathering rates and carbon storage along a climosequence of soils developed from contrasting granites in northeast brazil. Geoderma, 284:1–12, 2016. T. A. de Assis and F. D. A. A. Reis. Smoothening in thin-film deposition on rough substrates. Phys. Rev. E, 92:052405, 2015.

33

863 864 865

866 867

868 869 870

871 872

873 874

875 876

877 878 879

880 881 882

883 884 885

886 887 888 889

890 891

892 893

N. H. de Leeuw, S. C. Parker, and J. H. Harding. Molecular dynamics simulation of crystal dissolution from calcite steps. Phys. Rev. B, 60: 13792–13799, 1999. S. Emmanuel and Y. Levenson. Limestone weathering rates accelerated by micron-scale grain detachment. Geology, 42:751–754, 2014. P. Feng, A. S. Brand, L. Chen, and J. W. Bullard. In situ nanoscale observations of gypsum dissolution by digital holographic microscopy. Chem. Geol., 460:25–36, 2017. C. Fischer and A. Luttge. Beyond the conventional understanding of waterrock reactivity. Earth. Planet. Sci. Lett., 457:100–105, 2017. C. Fischer, R. S. Arvidson, and A. Luttge. How predictable are dissolution rates of crystalline material. Geochim. Cosmochim. Acta, 98:177–185, 2012. C. Fischer, I. Kurganskaya, T. Schafer, and A. Luttge. Variability of crystal surface reactivity: What do we know? Appl. Geochem., 43:132–157, 2014. J. M. Gautier, E. H. Oelkers, and J. Schott. Are quartz dissolution rates proportional to BET surface areas? Geochim. Cosmochim. Acta, 65:1059– 1070, 2001. J. R. A. Godinho, S. Piazolo, and T. Balic-Zunic. Importance of surface structure on dissolution of fluorite: Implications for surface dynamics and dissolution rates. Geochim. Cosmochim. Acta, 126:398–410, 2014. E. M. Hausrath, A. K. Navarre-Sitchler, P. B. Sak, J. Z. Williams, and S. L. Brantley. Soil profiles as indicators of mineral weathering rates and organic interactions for a pennsylvania diabase. Chem. Geol., 290:89–100, 2011. G. Jordan and W. Rammensee. Dissolution rates of calcite (1014) obtained by scanning force microscopy: Microtopography-based dissolution kinetics on surfaces with anisotropic step velocities. Geochim. Cosmochim. Acta, 62:941–947, 1998. M. Kardar, G. Parisi, and Y.-C. Zhang. Dynamic scaling of growing interfaces. Phys. Rev. Lett., 56:889–892, 1986. J. Krug. Universal finite-size effects in the rate of growth processes. J. Phys. A: Math. Gen., 23:L987–L994, 1990. 34

894 895

896 897 898

899 900

901 902 903

904 905 906

907 908

909 910

911 912

913 914 915

916 917

918 919 920

921 922 923

J. Krug. Origins of scale invariance in growth processes. Adv. Phys., 46: 139–282, 1997. I. Kurganskaya and A. Luttge. A comprehensive stochastic model of phyllosilicate dissolution: Structure and kinematics of etch pits formed on muscovite basal face. Geochim. Cosmochim. Acta, 120:545–560, 2013. I. Kurganskaya and A. Luttge. Kinetic Monte Carlo approach to study carbonate dissolution. J. Phys. Chem. C, 120:6482–6492, 2016. N. Laanait, E. B. R. Callagon, Z. Zhang, N. C. Sturchio, S. S. Lee, and P. Fenter. X-ray driven reaction front dynamics at calcite-water interfaces. Science, 349:1330–1334, 2015. J. F. Lardge, D. M. Duffy, M. J. Gillan, and M. Watkins. jab initio( simu-) lations of the interaction between water and defects on the calcite 1014 surface. J. Phys. Chem. C, 114:2664–2668, 2010. A. C. Lasaga and A. E. Blum. Surface chemistry, etch pits and mineral-water reactions. Geochim. Cosmochim. Acta, 50:2363–2379, 1986. A. C. Lasaga and A. Luttge. Variation of crystal dissolution rate based on a dissolution stepwave model. Science, 291:2400–2404, 2001. A. C. Lasaga and A. Luttge. Mineralogical approaches to fundamental crystal dissolution kinetics. Am. Mineral., 89:527–540, 2004. Y. Levenson and S. Emmanuel. Pore-scale heterogeneous reaction rates on a dissolving limestone surface. Geochim. Cosmochim. Acta, 119:188–197, 2013. Y. Liang and D. R. Baer. Anisotropic dissolution at the CaCO3 (1014) water interface. Surf. Sci., 373:275–287, 1997. A. Luttge, U. Winkler, and A. C. Lasaga. Interferometric study of the dolomite dissolution: A new conceptual model for mineral dissolution. Geochim. Cosmochim. Acta, 67:1099–1116, 2003. S. Majaniemi, T. Ala-Nissila, and J. Krug. Kinetic roughening of surfaces: Derivation, solution, and application of linear growth equations. Phys. Rev. B, 53:8071–8082, 1996. 35

924 925 926

927 928 929

930 931

932 933 934

935 936 937 938

939 940 941

942 943

944 945 946 947

948 949 950

951 952 953 954

J. M. McCoy and J. P. LaFemina. Kinetic monte carlo investigation of pit formation at the CaCO3 (1014) surface-water interface. Surf. Sci., 373: 288–299, 1997. P. Meakin and K. M. Rosso. Simple kinetic Monte Carlo models for dissolution pitting induced by crystal defects. J. Chem. Phys., 129:204106, 2008. M. Michaelis, C. Fischer, L. C. Ciacchi, and A. Luttge. Variability of zinc oxide dissolution rates. Environ. Sci. Technol., 51:4297–4305, 2017. D. A. Mirabella, G. P. Su´arez, M. P. Su´arez, and C. M. Aldao. Silicon wet etching: Hillock formation mechanisms and dynamic scaling properties. Physica A, 395:105–111, 2014. T. T. T. Nguyen, D. Bonamy, L. P. Van, J. Cousty, and L. Barbier. Scaling and universality in the kinetic smoothening of interfaces: Application to the analysis of the relaxation of rough vicinal steps of an oxide surface. EPL, 89:60005, 2010. A. Pelmenschikov, J. Leszczynski, and L. G. M. Pettersson. Mechanism of dissolution of neutral silica surfaces: Including effect of self-healing. J. Phys. Chem. A, 105:9528–9532, 2001. I. Pignatelli, A. Kumar, M. Bauchy, and G. Sant. Topological control on silicates’ dissolution kinetics. Langmuir, 32:4434–4439, 2016. M. Pollet-Villard, D. Daval, P. Ackerer, G. D. Saldi, B. Wild, K. G. Knauss, and B. Fritz. Does crystallographic anisotropy prevent the conventional treatment of aqueous mineral reactivity? a case study based on K-feldspar dissolution kinetics. Geochim. Cosmochim. Acta, 190:294–308, 2016a. M. Pollet-Villard, D. Daval, B. Fritz, K. G. Knauss, G. Schafer, and P. Ackerer. Influence of etch pit development on the surface area and dissolution kinetics of the orthoclase (001) surface. Chem. Geol., 447:79–92, 2016b. P. Raiteri, J. D. Gale, D. Quigley, and P. M. Rodger. Derivation of an accurate force-field for simulating the growth of calcium carbonate from aqueous solution: A new model for the calcite-water interface. J. Phys. Chem. C, 114:5997–6010, 2010.

36

957

E. Ruiz-Agudo and C. V. Putnis. Direct observations of mineral-fluid reactions using atomic force microscopy: the specific example of calcite. Mineral. Mag., 76:227–253, 2012.

958

J. H. Russ. Fractal Surfaces. Springer US, New York, USA, 1994.

955 956

959 960 961

962 963 964

965 966 967

968 969

970 971 972

973 974 975 976

977 978 979 980

981 982

983 984 985 986

G. D. Saldi, M. Voltolini, and K. G. Knauss. Effects of surface orientation, fluid chemistry and mechanical polishing on the variability of dolomite dissolution rates. Geochim. Cosmochim. Acta, 206:94–111, 2017. F. A. Silveira and F. D. A. A. Reis. Detachment of non-dissolved clusters and surface roughening in solid dissolution. Electrochim. Acta, 111:1–8, 2013. M. E. Smith, K. G. Knauss, and S. R. Higgins. Effects of crystal orientation on the dissolution of calcite by chemical and microscopic analysis. Chem. Geol., 360-361:10–21, 2013. B. Wehrli. Monte Carlo simulations of surface morphologies during mineral dissolution. J. Coll. Int. Sci., 132:230–242, 1989. A. F. White and S. L. Brantley. The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field? Chem. Geol., 202:479–506, 2003. A. F. White, M. S. Schulz, C. R. Lawrence, D. V. Vivit, and D. A. Stonestrom. Long-term flow-through column experiments and their relevance to natural granitoid weathering rates. Geochim. Cosmochim. Acta, 202: 190–214, 2017. M. Wolthers, D. Di Tommaso, Z. Du, and N. H. Leeuw. Calcite surface structure and reactivity: molecular dynamics simulations and macroscopic surface modelling of the calcitewater interface. Phys. Chem. Chem. Phys., 14:15145–15157, 2012. B. Zareeipolgardani, A. Piednoir, and J. Colombani. Gypsum dissolution rate from atomic step kinetics. J. Phys. Chem. C, 121:9325–9330, 2017. L. Zhang and A. Luttge. Morphological evolution of dissolving feldspar particles with anisotropic surface kinetics and implications for dissolution rate normalization and grain size dependence: A kinetic modeling study. Geochim. Cosmochim. Acta, 73:6757–6770, 2009. 37

987 988 989

Y. Zhao, G.-C. Wang, and T.-M. Lu. Characterization of Amorphous and Crystalline Rough Surface: Principles and Applications. Academic Press, San Diego, CA, USA, 2001.

38