Click here to donate to keep PhreeqcUsers open
Welcome,
Guest
. Please
login
or
register
.
Did you miss your
activation email
?
1 Hour
1 Day
1 Week
1 Month
Forever
Login with username, password and session length
Forum Home
Login
Register
PhreeqcUsers Discussion Forum
»
Conceptual Models
»
Kinetics and rate controlling factors
»
Question about Kinetic 1D modeling
« previous
next »
Print
Pages: [
1
]
Go Down
Author
Topic: Question about Kinetic 1D modeling (Read 658 times)
Jeonghwan Hwang
Top Contributor
Posts: 77
Question about Kinetic 1D modeling
«
on:
June 27, 2024, 09:00:59 AM »
I have a question about Kinetic 1D modeling.
Both Kinetics and Transport include time information through the -steps and -time_step parameters, respectively.
I am curious about how the calculations are performed in this case.
From my understanding, it seems that the calculation first proceeds with the Kinetics for 100 years, followed by a shift where the fluid flow for 100 years is calculated sequentially.
Is this correct? Or are all the calculations performed at once?
Additionally, if there are any issues with my model (in any part), I would appreciate your comments.
Thank you for reading
Sincerely,
Jeonghwan Hwang
My code is below
================================================================
PHASES
#LLNL based ToughReact
Na-Montmorillonite
Na0.33Mg0.33Al1.67Si4O10(OH)2 = + 0.33 Mg++ + 0.33 Na+ + 1.67 AlO2- + 0.66 H2O + 4 SiO2 + 0.68 H+
log_k -34.622 # at 25 cel
-delta_H 133.98 kJ/mol # Calculated enthalpy of reaction
Illite
K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2 = 0.6 K+ + 0.25 Mg++ + 2.3 AlO2- + 0.4 H2O + 3.5 SiO2 + 1.2 H+
log_k -42.6981 # at 25 cel
-delta_H 165.83 kJ/mol # Calculated enthalpy of reaction
END
RATES
Na-Montmorillonite # from ToughReact
-start
1 rem Toughreact code
2 rem Miner-2 montmorillonite 1 (kinetic) 3 (dissolution and precipitation) 1 (endmember of the solid solutions) 0 (mineral formed with wet block)
3 rem Miner-2.1 1.6596E-13 (coefficient K at 25 cel, mol/m2/sec) 2 (rate constant dependent on pH, additional equation) 1 (exponents n) 1 (exponents e) 35 (activation energy, kJ/mol) 0 0 0 (should be st to zero)
4 rem Miner-2.1.2 2 (number of additional mechanisms of rate constant)
5 rem Miner-2.1.2.1 1.0471E-11 (ki at 25 cel) 23.6 (Activation energy) 1 (number of species involved reaction) 'h+' (name of species) 0.34 (power term)
6 rem Miner-2.1.2.1 3.0200e-17 (ki at 25 cel) 58.9 (Activation energy) 1 (number of species involved reaction) 'h+' (name of species) -0.40 (power term)
7 rem Miner-2.2 1.6596E-13 2 1 1 35 0 0 0 (Same with Miner-2.1) 1E-6 (initial volume fraction, Vmineral/Vsolid, default 1E-15) 0 (precipitation law index)
8 rem Miner-3 0 (log q/k gap) 0 (temperature in cel in which to reducing gap) 0
#Miner-2.1
10 k1 = 1.6596E-13 #mol/m2/sec
11 Ea1 = 35 #kJ/mol
12 n1 = 1
13 e1 = 1
#Miner-2.1.2.1-1
20 k2 = 1.0471E-11
21 Ea2 = 23.6
22 n2 = 0.34
#Miner-2.1.2.1-2
30 k3 = 3.02e-17
31 Ea3 = 58.9
32 n3 = -0.40
#additional information
40 R = 8.31451 #J/K/mol
#Surface Area has a problem
41 rem Ar (surface area) = (Vfrac*Am(surface area of minerals) + Aprc(Precursor surface area))/Cw(wetted-surface conversion factor)
42 Vfrac = 1E-6 #Volume fraction
43 rem no sufficient information in Toughreact code. Thus, SSA should be considered as an adjustable parameter, See line 63
#Kinetic equation
50 Rate1 = k1*EXP(-Ea1/R/1000*(1/TK-1/298.15))
51 Rate2 = K2*EXP(-Ea2/R/1000*(1/TK-1/298.15))*ACT("H+")^n2
53 Rate3 = K3*EXP(-Ea3/R/1000*(1/TK-1/298.15))*ACT("H+")^n3
54 Rate = Rate1+Rate2+Rate3
#Kinetic Calculation
60 SR_Mineral = SR("Na-Montmorillonite") # SR>1 precipitation, SR<1 dissolution, SR=1 equilibrium
61 if (M<0) then goto 200 #stop calculation when current Montmorillonite mol (M) is negative
62 if (M=0 and SR_Mineral<1) then goto 200 #montmorillonite has to dissolved, but no current montmorillonite, so stop calculation
#63 if (M0<=0) then SSA=Parm(1) else SSA=Parm(1)*(M/M0)^0.67 #Calculate SSA, adjustable parameter
63 SSA=Parm(1) #Parm = 0.1516 m2/g - Tough
64 if (SSA<=0) then SSA = 1
65 Kinetic = Rate*SSA*(1-(SR_Mineral^n1))^e1
100 moles= Kinetic*Time
200 save moles
-end
Illite # from ToughReact
-start
1 rem Toughreact code # Illite Kinetic parameter is same as Na-Montmorillonite
2 rem Miner-2 Illite 1 (kinetic) 3 (dissolution and precipitation) 1 (endmember of the solid solutions) 0 (mineral formed with wet block)
3 rem Miner-2.1 1.6596E-13 (coefficient K at 25 cel, mol/m2/sec) 2 (rate constant dependent on pH, additional equation) 1 (exponents n) 1 (exponents e) 35 (activation energy, kJ/mol) 0 0 0 (should be st to zero)
4 rem Miner-2.1.2 2 (number of additional mechanisms of rate constant)
5 rem Miner-2.1.2.1 1.0471E-11 (ki at 25 cel) 23.6 (Activation energy) 1 (number of species involved reaction) 'h+' (name of species) 0.34 (power term)
6 rem Miner-2.1.2.1 3.0200e-17 (ki at 25 cel) 58.9 (Activation energy) 1 (number of species involved reaction) 'h+' (name of species) -0.40 (power term)
7 rem Miner-2.2 1.6596E-13 2 1 1 35 0 0 0 (Same with Miner-2.1) 1E-6 (initial volume fraction, Vmineral/Vsolid, default 1E-15) 0 (precipitation law index)
8 rem Miner-3 0 (log q/k gap) 0 (temperature in cel in which to reducing gap) 0
#Miner-2.1
10 k1 = 1.6596E-13 #mol/m2/sec
11 Ea1 = 35 #kJ/mol
12 n1 = 1
13 e1 = 1
#Miner-2.1.2.1-1
20 k2 = 1.0471E-11
21 Ea2 = 23.6
22 n2 = 0.34
#Miner-2.1.2.1-2
30 k3 = 3.02e-17
31 Ea3 = 58.9
32 n3 = -0.40
#additional information
40 R = 8.31451 #J/K/mol
#Surface Area has a problem
41 rem Ar (surface area) = (Vfrac*Am(surface area of minerals) + Aprc(Precursor surface area))/Cw(wetted-surface conversion factor)
42 Vfrac = 1E-6 #Volume fraction
43 rem no sufficient information in Toughreact code. Thus, SSA should be considered as an adjustable parameter, See line 63
#Kinetic equation
50 Rate1 = k1*EXP(-Ea1/R/1000*(1/TK-1/298.15))
51 Rate2 = K2*EXP(-Ea2/R/1000*(1/TK-1/298.15))*ACT("H+")^n2
53 Rate3 = K3*EXP(-Ea3/R/1000*(1/TK-1/298.15))*ACT("H+")^n3
54 Rate = Rate1+Rate2+Rate3
#Kinetic Calculation
60 SR_Mineral = SR("Illite") # SR>1 precipitation, SR<1 dissolution, SR=1 equilibrium
61 if (M<0) then goto 200 #stop calculation when current Montmorillonite mol (M) is negative
62 if (M=0 and SR_Mineral<1) then goto 200 #montmorillonite has to dissolved, but no current montmorillonite, so stop calculation
#63 if (M0<=0) then SSA=Parm(1) else SSA=Parm(1)*(M/M0)^0.67 #Calculate SSA, adjustable parameter
63 SSA=Parm(1) #Parm = 0.1516 m2/g - Tough
64 if (SSA<=0) then SSA = 1
65 Kinetic = Rate*SSA*(1-(SR_Mineral^n1))^e1
100 moles= Kinetic*Time
200 save moles
-end
END
Selected_output
-reset false
-file 3. Illitization_1D_Kin.txt
-pH true
-Solution true
-time true
User_Punch
headings Temp Naconc Kconc Mgconc Siconc Alconc SRMont SRIlli Mont Illi
-start
10 punch TC
20 punch TOT("Na") TOT("K") TOT("Mg") TOT("Si") TOT("Al")
30 punch SR("Na-Montmorillonite")
60 punch SR("Illite")
90 punch KIN("Na-Montmorillonite")
100 punch KIN("Illite")
-end
end
Solution 0
-units mol/kgw
pH 8
temp 130
K 0.1
Ca 1E-8
Mg 1E-8
Na 1E-8
Cl 1e-8 .charge
Si 1E-8 as SiO2
C 1E-8 as HCO3
S 1E-8 as SO4
Al 1E-8 as AlO2
Fe 1E-8
-water 1
End
SOLUTION 1-5
-units mol/kgw
temp 25
pH 8 .charge
K 1E-8
Ca 1E-8
Mg 1E-8
Na 1E-8
Cl 1e-8
Si 1E-8 as SiO2
C 1E-8 as HCO3
S 1E-8 as SO4
Al 1E-8 as AlO2
Fe 1E-8
-water 1
Kinetics 1-5
Na-Montmorillonite
-tol 1e-8
-m0 1E-2
-m 1E-2
-parms 151.6
Illite
-tol 1e-8
-m0 0
-m 0
-parms 151.6
-bad_step_max 1000
-steps 100 years in 10000
End
Transport
-cells 5
-shifts 5
-time_step 100 yr 10000.0
-flow_direction forward
-boundary_conditions flux flux
-diffusion_coefficient 0.0e-9
-punch_cells 2
-punch_frequency 1
PRINT; -reset false; -status false
Logged
dlparkhurst
Global Moderator
Posts: 3805
Re: Question about Kinetic 1D modeling
«
Reply #1 on:
June 27, 2024, 03:22:05 PM »
Please use the "#" button to insert code into your post.
The "END" statements are key to understanding the sequence of calculations. In the following excerpt from your script, there will be
(1) An initial solution calculation with solution 0.
(2) An initial solution calculation with solution 1 (which will be copied to solutions 2-5.
() Kinetics 1-5 are defined.
(3) A reaction calculation between solution 0 and kinetics 1 will occur for the time specified in kinetics 100 years in 10,000 steps. (0.01 year steps). If you are not sure about time steping, 10 steps or less would suffice. Note that kinetics 1 will change due to the reaction calculation, but kinetics 2-5 remain unaltered.
(4) A transport calculation will occur using the time step and shifts defined in TRANSPORT, 5 shifts of 100 years each. The 10000.0 following the time step defines the number of substeps to use in a multicomponent diffusion calculation; it is not used in your calculation.
The kinetics time step will be used for a reaction calculation involving kinetics; that is a calculation not involving TRANSPORT. The transport time step will be used in the transport calculation.
Code:
[Select]
SOLUTION 0
...
End
SOLUTION 1-5
...
Kinetics 1-5
...
End
Transport
...
Logged
Print
Pages: [
1
]
Go Up
« previous
next »
PhreeqcUsers Discussion Forum
»
Conceptual Models
»
Kinetics and rate controlling factors
»
Question about Kinetic 1D modeling