Question about Kinetic 1D modeling

[Jeonghwan Hwang]

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

[dlparkhurst]

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
...