changing porosity during dissolution/precipitation reactions

[nnegash]

Hi all!

So i am trying to model reactive transport of low pH solution through a carbonate, but would like to change porosity as the reaction progresses. To do so I use GET_POR(CELL_NO) in my RATES datablock to first determine the appropriate  volumes, and then use CHANGE_POR(new_porosity,CELL_NO) to update the porosity. However, I'm having trouble with the GET_POR since when i multiply it with another variable something goes wrong; some values go to infinity when it shouldnt. I'll include my RATES datablock to make it a bit clearer:

Code: [Select]

RATES
Calcite
 -start
1   A = PARM(CELL_NO+40) #Specific reactive surface area [m^2/L]

10  k1T25 = 0.89 #rate constant mol/m2/s (Chou et al., 1989)
11  k2T25 = 5.01*10^(-4) #(Chou et al., 1989)
12  k3T25 = 6.6*10^(-7) #(Chou et al., 1989)

20 Ea1 = 8.37 #activation energy kJ/mol (Plummer et al., 1978)
21 Ea2 = 41.88 #(Plummer et al., 1978)
22 Ea3 = 33.08 #(Plummer et al., 1978)
23 Rc = 8.314e-3 #ideal gas constant kJ/mol/K

30 k1 = k1T25*EXP((Ea1/Rc)*((1/298.15)-(1/TK))) #new rate constants accounting for temperature
31 k2 = k2T25*EXP((Ea2/Rc)*((1/298.15)-(1/TK)))
32 k3 = k3T25*EXP((Ea3/Rc)*((1/298.15)-(1/TK)))

40 Ksp = 10^LK_PHASE("Calcite") #equilibrium constant calcite [-] (PHREEQC database)

50 aH = ACT("H+")
51 aH2CO3 = ACT("CO2")
52 aCO3 = ACT("CO3-2")
53 aCa = ACT("Ca+2")

100 si_cc=LOG10(aCa*aCO3/Ksp)
101 IF (si_cc = 0) THEN GOTO 250
#101 IF (si_cc > -1e-14 and si_cc < 1e-14) THEN GOTO 250

120 rate = (k1*aH + k2*aH2CO3 + k3)*(1-((aCa*aCO3)/Ksp)) #Specific rate (mol/m2/s)
130 moles = rate * TIME * A

#calculate new porosity
140 Vtot = PARM(81) #Total volume of a cell (as per PoreFlow, determined externally); [mm3]
141 por = GET_POR(CELL_NO) #current porosity [-]
142 Vf = por*Vtot #Total fluid volume (Vf); [mm3]
143 delcal = moles * (Vf / 1000000) #Amount of dissolved calcite adjusted for the volume of a cell in PoreFlow [moles]
144 molmass = 100.0869 #molar mass of calcite; [g/mol]
145 dens = 2.711 #density of calcite; [g/cm3]
146 dV = (delcal * (molmass/dens))*1000 #Volume of dissolved calcite; times by 1000 to get cm3 in mm3 [mm3]
147 dPor = dV/Vtot #Change in porosity; change in volume over total volume
148 new_por = por + dPor #New porosity [-]
#149 CHANGE_POR(new_por,CELL_NO) #Changing porosity to new porosity
149 CHANGE_POR(Vf,CELL_NO)

250 SAVE moles
-end

USER_PUNCH 3
-headings -porosities
 -start
10 PUNCH GET_POR(CELL_NO)
 -end

SELECTED_OUTPUT 3
-file porosities_temp.txt
-high_precision true
-reset false
-user_punch true
 
In line 149 of RATES I output the Vf just to demonstrate what the output looks like from SELECTED_OUTPUT 3 and where the error is:

Code: [Select]

         -porosities
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.328158710000e-01
  1.#INF00000000e+00
 1.758791866573e+155
  1.967841821970e+99
  6.769666555052e+01
  6.769666555052e+01
  6.769666555052e+01
  6.769666555052e+01
  6.769666555052e+01
  6.769666555052e+01
  6.769666555052e+01

So as you can see, the 10 cells for timestep = 0 are outputted correctly, with the initial porosity being the same for all cells. But when outputting Vf after the timestep = 1 the values are weird. I've also checked by multiplying "por" in line 142 by 2 and it also gives really strange results, but when multiplied by 1 then there is no problem oddly enough. Is there anyway around this? I fear it may be out of my control, but maybe it can be fixed (if it can) in the next update? Anyway, thanks in advance!

Naod

[dlparkhurst]

Can you provide a complete input file that demonstrates the problem, including the TRANSPORT and KINETICS? Timing might be good as we are working on updates now.

[nnegash]

Hi David,

Sorry it took me a while to get back to you. I'm making a model that combines PHREEQC and a pore-network model in batch, and so it took me some time to rewrite my PHREEQC input file so that it can run by itself. Anywho, i managed to do it, but i have come to realise that what thought before was incorrect; there doesn't seem to be a problem with the GET_POR function. I was being a bit stupid because i was adding Vf to CHANGE_POR just to visualise it in the output, but i didnt realise i was multiplying it by Vtot every time (with Vtot=~7 at the time) and runnning the model for 500 shifts, so it was no surprise the porosity went to infinity. So that was a little dumb of me.

Nevertheless, the reason why i thought there was a problem in the first place hasnt changed. As i said before my aim is to determine the porosity change within phreeqc based solely on the amount of dissolved calcite. I have attached an input file and another file (called sample_PHREEQC_CELL_INFO.TXT) which is included in the input file with INCLUDE$ demonstrating this problem. In the input file in the RATES datablock i show two methods for calculating the associated change in volume for a given change in moles of calcite, scaled to the volume of a cell in the pore-network model. I've added some comments in the input file to elucidate what i mean. As far as i understand, these two methods should give (close to) the same change in porosity, but they do not. Furthermore, when i calculate the volume and porosity change from output file delCalcite.txt it comes closer to the second method. So initially i was using method 1 and was surprised by the porosities i was getting in phreeqc compared with the porosities that were outputted from the pore-network model, which prompted me to write the first message. But now I opted for method 2 since it appears to give better results. So i've managed to get a way around it (i think), but i'm still very curious to hear what you have to say!

Input file:

Code: [Select]

##########################################################################
# 1D kinetic dissolution model of calcite due to low pH solution injection
# Naod Negash and Flor Wassing
# 15-04-2019
##########################################################################

#Create injection solution with pure water in equilibrium with CO2
SOLUTION 0 Injection fluid in eq. with CO2
    temp      50
    pH        7 charge
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
EQUILIBRIUM_PHASES 0
    CO2(g)    2 10 # pCO2 = 10^2 atm
SAVE SOLUTION 0
END

SOLUTION 200 Background solution in eq. with calcite initially filling column
    temp      50
    pH        7 charge
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
EQUILIBRIUM_PHASES 200
Calcite 0 # SI(Calcite) = 0
SAVE SOLUTION 1-10
END
 
KINETICS 1
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 2
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 3
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 4
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 5
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 6
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 7
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 8
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 9
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS 10
Calcite
-tol   1e-8
-bad_step_max 100000
-m0    100
-parms
 
KINETICS_MODIFY 1
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 2
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 3
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 4
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 5
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 6
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 7
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 8
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 9
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT
 
KINETICS_MODIFY 10
INCLUDE$ sample_PHREEQC_CELL_INFO.TXT


RATES
Calcite
 -start
1  A = PARM(CELL_NO+20) #total reactive surface area [m^2]

10 k1T25 = 0.89 #rate constant mol/m2/s (Chou et al., 1989)
11 k2T25 = 5.01*10^(-4) #(Chou et al., 1989)
12 k3T25 = 6.6*10^(-7) #(Chou et al., 1989)

20 Ea1 = 8.37 #activation energy kJ/mol (Plummer et al., 1978)
21 Ea2 = 41.88 #(Plummer et al., 1978)
22 Ea3 = 33.08 #(Plummer et al., 1978)
23 Rc = 8.314e-3 #ideal gas constant kJ/mol/K

30 k1 = k1T25*EXP((Ea1/Rc)*((1/298.15)-(1/TK))) #new rate constants accounting for temperature
31 k2 = k2T25*EXP((Ea2/Rc)*((1/298.15)-(1/TK)))
32 k3 = k3T25*EXP((Ea3/Rc)*((1/298.15)-(1/TK)))

40 Ksp = 10^LK_PHASE("Calcite") #equilibrium constant calcite [-] (PHREEQC database)

50 aH = ACT("H+")
51 aH2CO3 = ACT("CO2")
52 aCO3 = ACT("CO3-2")
53 aCa = ACT("Ca+2")

100 si_cc=LOG10(aCa*aCO3/Ksp)
101 IF (si_cc = 0) THEN GOTO 250

120 rate = (k1*aH + k2*aH2CO3 + k3)*(1-((aCa*aCO3)/Ksp)) #Specific rate (mol/m2/s)
130 moles = rate * TIME * A

#calculate new porosity: method 1
139 dmol=moles
140 Vtot = PARM(51) #Total volume of a cell (as per PoreFlow, determined externally); [mm3]
141 por = GET_POR(CELL_NO) #Define porosity [-]
142 Vf = por*Vtot #Total fluid volume (Vf); [mm3]
143 delcal = dmol * (Vf / 1000000) #Amount of dissolved calcite adjusted for the volume of a cell in PoreFlow [moles]
144 molmass = 100.0869 #molar mass of calcite; [g/mol]
145 dens = 2.711 #density of calcite; [g/cm3]
146 dV = (delcal * (molmass/dens))*1000 #Volume of dissolved calcite; times by 1000 to get cm3 in mm3 [mm3]
147 dPor = dV/Vtot #Change in porosity; change in volume over total volume
148 new_por = por + dPor #New porosity [-]
149 CHANGE_POR(new_por,cell_no) #Changing porosity to new porosity

150 PUT(new_por,1)
151 PUT(dPor,2)
152 PUT(Vf,5)
153 PUT(dmol,7)

#calculate new porosity: method 2
160 ini_por = PARM(CELL_NO+40) #Import initial porosity from the parameters in KINETICS
161 IF (TIME=0) THEN Vf2 = ini_por*Vtot ELSE Vf2=GET(6) #Definte Vf depenging on the time: initially = ini_por*vtot but gets updated each time
162 dmol2=M0-M #Define dmol as the total amount of dissolved calcite since t=0
163 delcal2 = dmol2 * (Vf2 / 1000000) #Amount of dissolved calcite since t=0 adjusted for the volume of a cell in PoreFlow [moles]
164 dV2 = (delcal2 * (molmass/dens))*1000 #Volume of dissolved calcite since t=0; times by 1000 to get cm3 in mm3 [mm3]
165 dPor2 = dV2/Vtot #Change in porosity since t=0; change in volume over total volume
#New porosity; dPor2 is added to ini_por because now we calculated the associated porosity (i.e. volume) change from total amount of dissolved Calcite since t=0
166 new_por2 = ini_por + dPor2
#167 CHANGE_POR(new_por2,cell_no) #Changing porosity to new porosity
168 Vf2=new_por2*Vtot

170 PUT(new_por2,3)
171 PUT(dPor2,4)
172 PUT(Vf2,6)
173 PUT(dmol2,8)

250 SAVE moles
 -end

#Reading TRANSPORT datablock from trans.dmp
TRANSPORT
    -cells                 10
    -shifts                5000
    -time_step             0.001 # seconds
    -lengths               0.0032 #meters
    -dispersivities        0.00032
    -correct_disp          true
    -punch_cells        1-10
    -punch_frequency    5000
    -multi_d               true 1e-09 0.13 0 1
 -porosities
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00
  0.132815871E+00

USER_PUNCH 1
-headings porosities1 dPor1 porosities2 dPor2 Vf1 Vf2 moles M0-M
 -start
10 PUNCH GET(1)
20 PUNCH GET(2)
30 PUNCH GET(3)
40 PUNCH GET(4)
50 PUNCH GET(5)
60 PUNCH GET(6)
70 PUNCH GET(7)
80 PUNCH GET(8)
 -end

SELECTED_OUTPUT 1
-file OUTPUT.txt
-high_precision true
-reset false
-user_punch true

USER_PUNCH 2
-headings k_Calcite dk_Calcite
 -start
10 PUNCH M
20 PUNCH M-M0
 -end

SELECTED_OUTPUT 2
-file delCalcite.txt
-high_precision true
-reset false
-user_punch true

PRINT
-reset false

END


Thanks again!
Naod

[dlparkhurst]

I was afraid it would be a bit of work for you to make an input file.

First, GET(6) may be an error. You are using storage site 6 for all cells, so in the first encounter of GET(6) for cell 2 it is probably picking up the value from cell 1. In general, you are using one set of storage locations for all of the cells. If you store with the cell number you could avoid the possibility of cross contamination. PUT(vf2, 6, cell_no), for example.

Second, you have if TIME > 0, but TIME is an internal time substep. (You write rate*TIME*A for example.) It is adjusted during the reaction integration. With Runge-Kutta, it takes on values so that the RATES function is evaluated at different intervals between the starting time of the time step and ending time of the time step from which the function is integrated and an error estimate of the integration accuracy is made.  In addition, if errors are too great, the integration is restarted using smaller sub time steps. You probably intended the function TOTAL_TIME which is the current time of the integration starting from the beginning time (usually 0) of the integration.

I think you need to be careful how you adjust the porosity. Look careful at the example below. The moles of reaction evaluated at each sub interval is for the total time step, so if you add it at each sub time step, you will be adding approximately 6 or 8 times too much. The functions EQUI and KIN give the amount of onsider not using an incremental change with each sub time step.

But, to dig a little deeper, the porosity does not affect the chemistry calculation. So it does not really matter what the porosity is while the rates are being integrated. I would just set the porosity based on the amount of calcite at any given time, which can be determined by the function EQUI("Calcite") or the change from the initial amount of calcite at the start of the time step EQUI_DELTA("Calcite"). At the end of the time step, the porosity should be correct, and it will then be used for the next transport step.

Code: [Select]

RATES
Halite
10 k = 1
20 rate = k*(1 - SR("Halite"))
30 moles = rate * time
40 save moles
50 PRINT "TOTAL_TIME: ", TOTAL_TIME, "  TIME: ", TIME,"  KIN: ", M, "  moles: ", moles

SOLUTION
REACTION
NaCl 1
1e-9
SAVE solution 1
END
USE solution 1
KINETICS 1
-step  0.1
Halite
-M 100
USER_PRINT
10 PRINT "KIN: ", KIN("Halite")
END


Code: [Select]

TOTAL_TIME:             0   TIME:    1.0000e-01   KIN:           100   moles:    1.0000e-01
TOTAL_TIME:    2.0000e-02   TIME:    1.0000e-01   KIN:    9.9980e+01   moles:    9.9999e-02
TOTAL_TIME:    3.0000e-02   TIME:    1.0000e-01   KIN:    9.9970e+01   moles:    9.9998e-02
TOTAL_TIME:    6.0000e-02   TIME:    1.0000e-01   KIN:    9.9940e+01   moles:    9.9994e-02
TOTAL_TIME:    1.0000e-01   TIME:    1.0000e-01   KIN:    9.9900e+01   moles:    9.9984e-02
TOTAL_TIME:    8.7500e-02   TIME:    1.0000e-01   KIN:    9.9913e+01   moles:    9.9987e-02
TOTAL_TIME:    8.7500e-02   TIME:    1.0000e-01   KIN:    9.9900e+01   moles:    9.9984e-02


[nnegash]

Ahhhh yes, that's where it went wrong. I didn't realise the moles of reaction calculated at each subinterval was for the entire timestep, and that indeed explains why the porosity i was calculating was too high; i just ran it with dividing "moles" by 7 in my code (line 139) and both methods now give pretty much the same result. I'm still a bit confused though by the total time outputted from the script you shared. I see how the TIME is the same for each subinterval, explaining why porosity increased too much, but then i don't understand what the total_time is and how it goes up to 0.1 (the final time?) and then retreats to 0.0875.. I guess i'll have to read the manual more thoroughly!

Also, you are absolutely right, it's way easier to just change the porosity based on the amount of moles of kinetics reactant at any given time, i don't know why i didn't think of that. I just implemented it and works out perfectly, so much easier. Also, thank you for your suggestions for GET(6,cell_no) and IF TOTAL_TIME =0! As i said i'm still a little confused with" total_time" and "time" at the subintervals, i thought TIME was the amount of time of the subinterval itself, but i'll read the manual more carefully! Anway, thank you so much for your help! I feel like i understand the workings of PHREEQC a bit better now

[dlparkhurst]

It's been a long while since I wrote the kinetics parts, and I admit that I may have explained it incorrectly at various times over the years.

The Runge-Kutta method evaluates the rate functions, I think, at 6 different times over the interval, which are the times shown by TOTAL_TIME. These times are not monotonic, as is also evident. Rates based previous evaluations are used to estimate concentrations at any TOTAL_TIME. After evaluating at all the TOTAL_TIMES the method then uses some of the information to get the overall rate that is fifth order in the differential equation sense, and others to get a fourth order overall rate. The difference between the two estimates is used to estimate the error. If the error is too large, then the time interval is subdivided into more time segments, and the process is restarted. So the total number of evaluations may be 6 or 7 or many times that if the error estimate requires many sub intervals. You can understand a simpler integration scheme of predictor/corrector that estimates rates, steps forward, estimates rates again, and then uses an average of the two rates for the time step.

[nnegash]

Hmm i think i understand, so at every (sub)interval (TOTAL_TIME) the rate is evaluated for the full step interval (TIME), which in turn is used to determine the rate at the intervals, and then uses information of all evaluations to determine the overall rate over the full timestep, or something like that.

[dlparkhurst]

I think that is close. You can look up Runge Kutta integration; it is a standard method for explicit numerical integration. There are different orders of accuracy, PHREEQC uses 4/5 order accuracy.

PHREEQC also has an implicit method, -cvode, which evalutes rate derivatives numerically as part of the integration.

[nnegash]

Will do! Again, many thanks!