Granite dissolution and Clay precipitation

[N_Harpers]

Hello everybody,

I’m new to PHREEQC and am trying to model experiments that we will conduct in the lab which are related to a deep geothermal context in granites. It will be a flow-through setup so I tried to create a 1D TRANSPORT model in PHREEQC and I want to use variable fluid compositions (i.e. deionized water for a leaching experiment, and waters close to meteoric and deep brine compositions).

My main interests are:
1)   The development of the fluid composition
2)   The dissolution of the minerals making up the granite
3)   The precipitation of different clays and other minerals related to granite alteration

All of these at variable pressures, temperatures and fluid velocities over a time period of about a month. I’m taking into account the reaction kinetics following data from Palandri and Kharaka (2004), Marty et al. (2015) and Zhang et al. (2019) and use ThermoddemV1 as database. Furthermore, I also use a simplified formulation of the nucleation code by Fritz et al (2009) to also take into account nucleation and growth of the platy clays (without keeping track of nuclei populations).

The simulations at the moment work for runs at 25 °C but lead to problems at higher temperatures (I tried 50, 75 and 90 °C) using Runge-Kutta. When increasing The bad step max (to make runs possible at al) the calculations at some point reach some kind of endless loop: Stuck in one time step for one cell, the bad steps don’t count higher, while the OK-steps count upwards unendingly (can’t say if really infinitely, but the highest count I had there was something over 300000, after I had left the computer for a while).  I assume because of that that I should look again into the way I defined the kinetics, but also wanted to see if someone might have encountered something similar. Is there a strategy to maybe identify the problem? Using CVODE as solver led to the simulations at higher T to run to the end, but the results look extremely chaotic which makes me doubt their credibility.
I attached the input-deck.

Example for clay dissolution and precipitation in RATES Block:

Code: [Select]

Kaolinite
# kinetic data extracted from 06has/vil 03tou/nea 03koh/duf 05koh/bos
# kinetic data extracted from Kaolinite precipitation rate: 08yan/ste 93nag/las 97dev/sch
# surface area data extracted from 15mar/cla
# Confidence level: 4
-start
1 SRmin = SR("Kaolinite") # mineral saturation ratio
10 moles = 0 # initial moles of mineral
20 If (m = 0) and (SRmin <= 1) Then GoTo 400 # if neither dissolution nor precipitation possible then end
30 If (SRmin = 1) Then GoTo 400 # if system in equilibrium then end
40 S = 11.8 # average BET; suggested value in m2/g
50 Mm = 258.16 # molar mass in g/mol
60 Mv = 9.876E-05 # molar volume in m3/mol
70 If (SRmin > 1) Then GoTo 160 # if supersaturated go to precipitation
########## start dissolution bloc ##########
80 knu = 1.1E-14 * exp((-38000 / 8.314) * ((1 / TK) - (1 / 298.15)))
90 k1 = 7.5E-12 * exp((-43000 / 8.314) * ((1 / TK) - (1 / 298.15))) * (ACT("H+") ^ 0.51)
100 k2 = 2.5E-11 * exp((-46000 / 8.314) * ((1 / TK) - (1 / 298.15))) * (ACT("OH-") ^ 0.58)
110 k = knu + k1 + k2
120 Rsrf = S * m * Mm * k * ((1 - SRmin))
130 moles = Rsrf * Time
#### Do not dissolve more than there is solid present ####
140 IF (moles > M) THEN moles = M
150 GoTo 400
########## end dissolution bloc ##########
########## start precipitation bloc ##########
160 surf_energy = 0.1 # interfacial energy for lateral surface in J m-2 (after Fritz et al 2009; lateral s chosen to be 2 times basal))
165 sheet_thick = 7e-10 # thickness of one mineral layer in m
170 molvol = 1.64E-28 # molecular volume in m3
180 cntc_angl = 90 # mineral-substrate contact angle in degrees (assume hemi-spherical)
190 u = ( 2 * sqrt(3) * sheet_thick * surf_energy^2 * molvol ) / ( (1.38E-23 * TK)^2 ) #
200 J0 = 1E20  # nucleation rate in nuclei m-2 sec-1 (Based on Fritz et al 2009)
210 SRcrit = exp( u / LOG( J0 ) ) # critical saturation threshold
220 IF ( SRmin < SRcrit ) Then GoTo 290 # if saturation index below threshold no nucleation occurs
230 rate_nuc = J0 * exp( -u / LOG(SRmin) ) # nucleation rate in nuclie sec-1 multiplied by time to account for nucleation over time
240 IF ( rate_nuc < 1 ) THEN GoTo 290 # condition that rate needs to be bigger than 1 nucl per sec (Arbitrary)
250 nj = u / (LOG(SRmin)^2) # number of growth units in critical nuclie radius
260 vol = nj * molvol # critical nuclie volume
270 moles_nuc = ( -rate_nuc * Time * vol ) / Mv # moles of nuclie formed
280 GoTo 300
290 moles_nuc = 0
#### kinetic data for growth ####
300 IF ( m = 0 ) THEN GoTo 350
310 knu = 5.5E-13 * exp((-66000 / 8.314) * ((1 / TK) - (1 / 298.15)))
320 kpre = (-1) * knu
330 Rsrf = S * m * Mm * kpre * (ABS((SRmin ^ 0.06) - 1) ^ 1.68)
350 moles = moles_nuc + Rsrf * Time
#### Do not precipitate more than the elements in solution ####
360 maxMol = TOT("Al")/2
370 IF (maxMol > TOT("Si")/2) THEN maxMol = TOT("Si")/2
380 IF (maxMol > TOT("H")/6) THEN maxMol = TOT("H")/6
390 IF (maxMol < -moles) THEN moles = -maxMol
########## end precipitation bloc ##########
400 Save moles
-end


[dlparkhurst]

You need to include your RATES file (MINERAL KINETICS DATABASE_NH_MOD_NUCL3) (and database if possible) for me to be able to run your script.