RATESCalcite -start1 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/K30 k1 = k1T25*EXP((Ea1/Rc)*((1/298.15)-(1/TK))) #new rate constants accounting for temperature31 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 250120 rate = (k1*aH + k2*aH2CO3 + k3)*(1-((aCa*aCO3)/Ksp)) #Specific rate (mol/m2/s)130 moles = rate * TIME * A#calculate new porosity140 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 volume148 new_por = por + dPor #New porosity [-]#149 CHANGE_POR(new_por,CELL_NO) #Changing porosity to new porosity149 CHANGE_POR(Vf,CELL_NO)250 SAVE moles-endUSER_PUNCH 3 -headings -porosities -start10 PUNCH GET_POR(CELL_NO) -endSELECTED_OUTPUT 3 -file porosities_temp.txt -high_precision true -reset false -user_punch true

-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

########################################################################### 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 CO2SOLUTION 0 Injection fluid in eq. with CO2 temp 50 pH 7 charge pe 4 redox pe units mol/kgw density 1 -water 1 # kgEQUILIBRIUM_PHASES 0 CO2(g) 2 10 # pCO2 = 10^2 atmSAVE SOLUTION 0ENDSOLUTION 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 # kgEQUILIBRIUM_PHASES 200 Calcite 0 # SI(Calcite) = 0SAVE 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 RATESCalcite -start1 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/K30 k1 = k1T25*EXP((Ea1/Rc)*((1/298.15)-(1/TK))) #new rate constants accounting for temperature31 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 250120 rate = (k1*aH + k2*aH2CO3 + k3)*(1-((aCa*aCO3)/Ksp)) #Specific rate (mol/m2/s)130 moles = rate * TIME * A#calculate new porosity: method 1139 dmol=moles140 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 volume148 new_por = por + dPor #New porosity [-]149 CHANGE_POR(new_por,cell_no) #Changing porosity to new porosity150 PUT(new_por,1)151 PUT(dPor,2)152 PUT(Vf,5)153 PUT(dmol,7)#calculate new porosity: method 2160 ini_por = PARM(CELL_NO+40) #Import initial porosity from the parameters in KINETICS161 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 time162 dmol2=M0-M #Define dmol as the total amount of dissolved calcite since t=0163 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=0166 new_por2 = ini_por + dPor2#167 CHANGE_POR(new_por2,cell_no) #Changing porosity to new porosity168 Vf2=new_por2*Vtot170 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.dmpTRANSPORT -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 -start10 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) -endSELECTED_OUTPUT 1 -file OUTPUT.txt -high_precision true -reset false -user_punch trueUSER_PUNCH 2 -headings k_Calcite dk_Calcite -start10 PUNCH M20 PUNCH M-M0 -endSELECTED_OUTPUT 2 -file delCalcite.txt -high_precision true -reset false -user_punch true PRINT -reset false END

RATES Halite10 k = 120 rate = k*(1 - SR("Halite"))30 moles = rate * time40 save moles50 PRINT "TOTAL_TIME: ", TOTAL_TIME, " TIME: ", TIME," KIN: ", M, " moles: ", molesSOLUTIONREACTIONNaCl 11e-9SAVE solution 1ENDUSE solution 1KINETICS 1-step 0.1Halite-M 100USER_PRINT10 PRINT "KIN: ", KIN("Halite")END

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