Computation time for advection with kinetics

[Flopi]

Dear forum,
I started recently to look into 1D models and I have a big problem of calculation time. If I run the following code on my (very powerful) machine, the first shift is not finished after several hours.
So basically, my questions is: did I do something wrong that slows the calculation and is there a way to smartly accelerate it ? My biggest source of confusion is that I'm not sure how the number of cells plays a role in all of this: from my understanding, if I increase the number of cells, the duration of each shift is shortened which leads to shorter calculation time per cell, but hardly shorter overall because of the increase of the number of cells. And the other hand, decreasing the number of cells increases the calculation time per cell...
I'm also using the graphical interface Phreeqc v3, not sure that impacts anything though.

Code: [Select]

PHASES
Olivine-Fo90
Mg1.796Fe0.2Ni0.004SiO4 + 4 H+ = 1.0H4SiO4 + 1.796Mg+2 + 0.2 Fe+2 + 0.004Ni+2
log_k 27.646
-analytic -7.3366336e2 -0.1113871 4.90932603e4 2.63819198e2 -2.0395157e6


SOLUTION 0
temp      200
pressure 200
pe        0
redox     pe
units mmol/kgw
    pH        7 charge
    Cl 1200
Na 1200
-density 1
   -water    1 # kg

SAVE SOLUTION 0
END

USE SOLUTION 0
COPY SOLUTION 0 1-10000
SAVE SOLUTION 1-10000
END

RATES
Olivine-Fo90 #From Summary Kinetics Rates -- Rimstidt 2012
-start
10 si_Fo = SI("Olivine-Fo90")
20 if si_FO > 0 then goto 200 #On evite la précipitation
30 F = 1
40 SSA = 2*59400
50 SA = F * SSA #Masse fraction * specific surface area (2 m2/g * 100 kg)
60 if (M<=0 and si_FO < 0) then goto 200
70 rf = 0.0217 * exp(-60900/(8.314*TK))*ACT("H+")^0.22 #Mécanisme basique mol/m2/s
80 rate = rf * (1 - 10^(si_FO)) * SA
90 moles = rate * TIME
#100 PRINT "moles Oliv" moles
200 SAVE moles
-end

Lizardite #Daval
-start
80 F = 1
90 SSA = 2*59400
100 SA = F*SSA
110 if ACT("H+") > 10^(-6.7) then rf = 10^(-2.27)*exp(-42000/8.314/TK)*ACT("H+")^0.53 else rf = 10^(-12.08)*exp(-42000/8.314/TK)
120 rate = rf * (1 - SR("Lizardite")) *SA
200 SAVE rate * time
-end
END


ADVECTION
    -cells 10000
    -shifts 10
    -time_step 10000 # seconds
    -print_cells 1 1000 2000 3000 4000 5000 6000
                 7000 8000 9000
    -punch_cells 1 1000 2000 3000 4000 5000 6000
                 7000 8000 9000
    -warnings false

SELECTED_OUTPUT
-file Advection1D-2
-distance true
-time true
-temperature true
-ionic_strength true
-water true
-totals H(0) Ni Fe(2) Fe(3) C(4)
-equilibrium_phases Magnetite NiFe2O4 Brucite
-Saturation_indices Olivine-Fo90 Lizardite
-kinetic_reactants Olivine-Fo90 Lizardite

EQUILIBRIUM_PHASES 1-10000
Brucite 0 0
Magnetite 0 0
NiFe2O4 0 0

KINETICS 1-10000 #quantities are for 29.7 kg or rock (10% porosity, 3.3kg/L)
Lizardite
-m 1e-6
Olivine-Fo90
-m 194


Thanks a lot for the help,

EDIT: just ran the same simulation at 25°C, and it ran for 760s. It seems the higher temperature slows down the calculations. Would that make sense and is there a way to accelerate it ?

[dlparkhurst]

I admire your patience, but I would never use 10000 cells. Let's start with a batch reaction, where the solution, equilibrium phases, and kinetics react in a beaker for a short period of time. I am using wateq4f.dat with some PHASES definitions that I found. If you use a sit or pitzer database, calculations will be slower.

I have used -cvode in KINETICS, which works best when reactions have widely differing rates. The results show that olivine comes to equilibrium in less than 10 minutes.

Code: [Select]

PHASES
Olivine-Fo90
Mg1.796Fe0.2Ni0.004SiO4 + 4 H+ = 1.0H4SiO4 + 1.796Mg+2 + 0.2 Fe+2 + 0.004Ni+2
log_k 27.646
-analytic -7.3366336e2 -0.1113871 4.90932603e4 2.63819198e2 -2.0395157e6
Lizardite
Mg3Si2O5(OH)4     = 3.000Mg+2     - 6.000H+     + 2.000H4SiO4     + 1.000H2O
     log_k    33.100     
     delta_h -247.218    #kJ/mol       
     # Enthalpy of formation:           -4362        #kJ/mol        #04EVA
     -analytic -1.02107E+1 0E+0 1.29131E+4 0E+0 0E+0
NiFe2O4
NiFe2O4 + 8 H+ = Ni+2 + 2 Fe+3 + 4 H2O
log_k 9.7876
-delta_H -215.338 kJ/mol
# deltafH -1081.15 kJ/mol
-analytic -1.4322e2 -2.9429e-2 1.4518e4 4.5698e1 2.4658e2
# Range 0-200
-Vm 44.89 # Webmineral.com
# Extrapol Constant H approx
# Ref RHF79
RATES
Olivine-Fo90 #From Summary Kinetics Rates -- Rimstidt 2012
-start
10 si_Fo = SI("Olivine-Fo90")
20 if si_FO > 0 then goto 200 #On evite la précipitation
30 F = 1
40 SSA = 2*59400
50 SA = F * SSA #Masse fraction * specific surface area (2 m2/g * 100 kg)
60 if (M<=0 and si_FO < 0) then goto 200
70 rf = 0.0217 * exp(-60900/(8.314*TK))*ACT("H+")^0.22 #Mécanisme basique mol/m2/s
80 rate = rf * (1 - 10^(si_FO)) * SA
90 moles = rate * TIME
#100 PRINT "moles Oliv" moles
200 SAVE moles
-end

Lizardite #Daval
-start
80 F = 1
90 SSA = 2*59400
100 SA = F*SSA
110 if ACT("H+") > 10^(-6.7) then rf = 10^(-2.27)*exp(-42000/8.314/TK)*ACT("H+")^0.53 else rf = 10^(-12.08)*exp(-42000/8.314/TK)
120 rate = rf * (1 - SR("Lizardite")) *SA
200 SAVE rate * time
-end
END
SOLUTION 1
temp      200
pressure 200
pe        0
redox     pe
units mmol/kgw
    pH        7 charge
    Cl 1200
Na 1200
-density 1
   -water    1 # kg
END
EQUILIBRIUM_PHASES 1
Brucite 0 0
Magnetite 0 0
NiFe2O4 0 0
KINETICS 1 #quantities are for 29.7 kg or rock (10% porosity, 3.3kg/L)
-cvode
Lizardite
-m 1e-6
Olivine-Fo90
-m 194
-step 1 1e1 1e2 1e3 1e4 1e5 1e6
END
INCREMENTAL_REACTIONS true
USE solution 1
USE equilibrium_phases 1
USE kinetics 1
USER_GRAPH 1
    -headings               time Olivine_si Lizardite_si Brucite Magnitite NiFe2O4
    -axis_titles            "Time, minutes" "SI" "Moles"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X (TOTAL_TIME/60)
20 GRAPH_Y SI("Olivine-Fo90"), SI("Lizardite")
30 GRAPH_SY EQUI("Brucite"), EQUI("Magnetite"),EQUI("NiFe2O4")
  -end
    -active                 true
END
If the rate of olivine is that fast, it is numerically much simpler if it is moved to EQUILIBRIUM_PHASES from  KINETICS. By removing one of the reactions from KINETICS, it is possible to use the default Runge-Kutta method, which should be faster than -cvode. This second script runs the batch reaction 100,000 years using olivine in equilibrium phases. It takes about 10,000 years for lizardite to precipitate to equilibrium, as brucite, magnetite, and NiFe2O4 are precipitated.

Code: [Select]

PHASES
Olivine-Fo90
Mg1.796Fe0.2Ni0.004SiO4 + 4 H+ = 1.0H4SiO4 + 1.796Mg+2 + 0.2 Fe+2 + 0.004Ni+2
log_k 27.646
-analytic -7.3366336e2 -0.1113871 4.90932603e4 2.63819198e2 -2.0395157e6
Lizardite
Mg3Si2O5(OH)4     = 3.000Mg+2     - 6.000H+     + 2.000H4SiO4     + 1.000H2O
     log_k    33.100     
     delta_h -247.218    #kJ/mol       
     # Enthalpy of formation:           -4362        #kJ/mol        #04EVA
     -analytic -1.02107E+1 0E+0 1.29131E+4 0E+0 0E+0
NiFe2O4
NiFe2O4 + 8 H+ = Ni+2 + 2 Fe+3 + 4 H2O
log_k 9.7876
-delta_H -215.338 kJ/mol
# deltafH -1081.15 kJ/mol
-analytic -1.4322e2 -2.9429e-2 1.4518e4 4.5698e1 2.4658e2
# Range 0-200
-Vm 44.89 # Webmineral.com
# Extrapol Constant H approx
# Ref RHF79
RATES
Olivine-Fo90 #From Summary Kinetics Rates -- Rimstidt 2012
-start
10 si_Fo = SI("Olivine-Fo90")
20 if si_FO > 0 then goto 200 #On evite la précipitation
30 F = 1
40 SSA = 2*59400
50 SA = F * SSA #Masse fraction * specific surface area (2 m2/g * 100 kg)
60 if (M<=0 and si_FO < 0) then goto 200
70 rf = 0.0217 * exp(-60900/(8.314*TK))*ACT("H+")^0.22 #Mécanisme basique mol/m2/s
80 rate = rf * (1 - 10^(si_FO)) * SA
90 moles = rate * TIME
#100 PRINT "moles Oliv" moles
200 SAVE moles
-end

Lizardite #Daval
-start
80 F = 1
90 SSA = 2*59400
100 SA = F*SSA
110 if ACT("H+") > 10^(-6.7) then rf = 10^(-2.27)*exp(-42000/8.314/TK)*ACT("H+")^0.53 else rf = 10^(-12.08)*exp(-42000/8.314/TK)
120 rate = rf * (1 - SR("Lizardite")) *SA
200 SAVE rate * time
-end
END
SOLUTION 1
temp      200
pressure 200
pe        0
redox     pe
units mmol/kgw
    pH        7 charge
    Cl 1200
Na 1200
-density 1
   -water    1 # kg
END
EQUILIBRIUM_PHASES 1
Brucite 0 0
Magnetite 0 0
NiFe2O4 0 0
      Olivine-Fo90 0 194
KINETICS 1 #quantities are for 29.7 kg or rock (10% porosity, 3.3kg/L)
#-cvode
Lizardite
-m 1e-6
# Olivine-Fo90
# -m 194
-step 3.15e5 3.15e6 3.15e7 3.15e8 3.15e9 3.15e10 3.15e11 3.15e12
END
INCREMENTAL_REACTIONS true
USE solution 1
USE equilibrium_phases 1
USE kinetics 1
USER_GRAPH 1
    -headings               time Olivine_si Lizardite_si Brucite Magnitite NiFe2O4
    -axis_titles            "Time, years" "SI" "Moles"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X (TOTAL_TIME/3.15e7)
20 GRAPH_Y SI("Olivine-Fo90"), SI("Lizardite")
30 GRAPH_SY EQUI("Brucite"), EQUI("Magnetite"),EQUI("NiFe2O4")
  -end
    -active                 true
END

It is a good idea not to rely only on graphs. You should look at the output as well. You will find at the end of the run, the system is extremely reducing with a large quantity of H2(a) and a small mass of water (0.1 kg). So, you may want to reassess whether the overall reaction is reasonable. A simple batch reaction showing results as a function of Lizardite precipitation may help in your evaluation.

Code: [Select]

PHASES
Olivine-Fo90
Mg1.796Fe0.2Ni0.004SiO4 + 4 H+ = 1.0H4SiO4 + 1.796Mg+2 + 0.2 Fe+2 + 0.004Ni+2
log_k 27.646
-analytic -7.3366336e2 -0.1113871 4.90932603e4 2.63819198e2 -2.0395157e6
Lizardite
Mg3Si2O5(OH)4     = 3.000Mg+2     - 6.000H+     + 2.000H4SiO4     + 1.000H2O
     log_k    33.100     
     delta_h -247.218    #kJ/mol       
     # Enthalpy of formation:           -4362        #kJ/mol        #04EVA
     -analytic -1.02107E+1 0E+0 1.29131E+4 0E+0 0E+0
NiFe2O4
NiFe2O4 + 8 H+ = Ni+2 + 2 Fe+3 + 4 H2O
log_k 9.7876
-delta_H -215.338 kJ/mol
# deltafH -1081.15 kJ/mol
-analytic -1.4322e2 -2.9429e-2 1.4518e4 4.5698e1 2.4658e2
# Range 0-200
-Vm 44.89 # Webmineral.com
# Extrapol Constant H approx
# Ref RHF79

SOLUTION 1
temp      200
pressure 200
pe        0
redox     pe
units mmol/kgw
    pH        7 charge
    Cl 1200
Na 1200
-density 1
   -water    1 # kg

EQUILIBRIUM_PHASES 1
Brucite 0 0
Magnetite 0 0
NiFe2O4 0 0
      Olivine-Fo90 0 194

REACTION 1
Lizardite 1
-10 in 10
USER_GRAPH 1
    -headings               rxn Water H2(aq) Magnetite Brucite
    -axis_titles            "Lizardite precipitated, moles" "Water, kg" "Moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X -RXN
20 GRAPH_Y TOT("water")
30 GRAPH_SY MOL("H2")*TOT("water"), EQUI("Magnetite"), EQUI("Brucite")
  -end
    -active                 true
END

But, if we forge ahead, we can then go back to an ADVECTION calculation. ADVECTION has no space dimension, so, if we use 20 cells, it may refer to 10 mm or 1 km, although the time step puts some reality on what distance is reasonable. Here is a simulation with 100-y time steps with 20 cells. So water at the exit of the column is about 2,000 years old. I think this simulation is a fairly complete description of the system. Olivine is ultimately dissolved completely and the water eventually equilibrates with lizardite. Although minerals continually dissolve and precipitate, the water composition reaches more-or-less a steady state between the point where olivine is present and the point where a cell is in equilibrium with lizardite.

Code: [Select]

PHASES
Olivine-Fo90
Mg1.796Fe0.2Ni0.004SiO4 + 4 H+ = 1.0H4SiO4 + 1.796Mg+2 + 0.2 Fe+2 + 0.004Ni+2
log_k 27.646
-analytic -7.3366336e2 -0.1113871 4.90932603e4 2.63819198e2 -2.0395157e6
Lizardite
Mg3Si2O5(OH)4     = 3.000Mg+2     - 6.000H+     + 2.000H4SiO4     + 1.000H2O
     log_k    33.100     
     delta_h -247.218    #kJ/mol       
     # Enthalpy of formation:           -4362        #kJ/mol        #04EVA
     -analytic -1.02107E+1 0E+0 1.29131E+4 0E+0 0E+0
NiFe2O4
NiFe2O4 + 8 H+ = Ni+2 + 2 Fe+3 + 4 H2O
log_k 9.7876
-delta_H -215.338 kJ/mol
# deltafH -1081.15 kJ/mol
-analytic -1.4322e2 -2.9429e-2 1.4518e4 4.5698e1 2.4658e2
# Range 0-200
-Vm 44.89 # Webmineral.com
# Extrapol Constant H approx
# Ref RHF79
RATES
Olivine-Fo90 #From Summary Kinetics Rates -- Rimstidt 2012
-start
10 si_Fo = SI("Olivine-Fo90")
20 if si_FO > 0 then goto 200 #On evite la précipitation
30 F = 1
40 SSA = 2*59400
50 SA = F * SSA #Masse fraction * specific surface area (2 m2/g * 100 kg)
60 if (M<=0 and si_FO < 0) then goto 200
70 rf = 0.0217 * exp(-60900/(8.314*TK))*ACT("H+")^0.22 #Mécanisme basique mol/m2/s
80 rate = rf * (1 - 10^(si_FO)) * SA
90 moles = rate * TIME
#100 PRINT "moles Oliv" moles
200 SAVE moles
-end

Lizardite #Daval
-start
80 F = 1
90 SSA = 2*59400
100 SA = F*SSA
110 if ACT("H+") > 10^(-6.7) then rf = 10^(-2.27)*exp(-42000/8.314/TK)*ACT("H+")^0.53 else rf = 10^(-12.08)*exp(-42000/8.314/TK)
120 rate = rf * (1 - SR("Lizardite")) *SA
200 SAVE rate * time
-end
END
SOLUTION 0-100
temp      200
pressure 200
pe        0
redox     pe
units mmol/kgw
    pH        7 charge
    Cl 1200
Na 1200
-density 1
   -water    1 # kg
END

KINETICS 1-100 #quantities are for 29.7 kg or rock (10% porosity, 3.3kg/L)
#-cvode
Lizardite
-m 1e-6
# Olivine-Fo90
# -m 194
-step 1e5
END
EQUILIBRIUM_PHASES 1-100
Brucite 0 0
Magnetite 0 0
NiFe2O4 0 0
      Olivine-Fo90 0 194
END
ADVECTION
-cells 10
-time_step 3.15e8
-shifts 15
-punch_frequency 5

USER_GRAPH 1
    -headings               time Olivine_si Lizardite_si Brucite Magnitite NiFe2O4
    -axis_titles            "Cell number" "SI" "Moles"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 GRAPH_X CELL_NO
20 GRAPH_Y SI("Olivine-Fo90"), SI("Lizardite")
30 GRAPH_SY EQUI("Brucite"), EQUI("Magnetite"),EQUI("NiFe2O4")
  -end
    -active                 true
END

[Flopi]

I am always amazed how fast and accurately you answer our random questions.
Thank you so very much !