Processes > Mixing
Mineral Dissolution with mixing
achichat:
Hi,
I am looking to simulate dissolution of olivine in seawater, and I want to perform mixing with an amount of seawater after each step.
--- Code: --- TITLE Example 1. Olivine dissolution with kinetic rate law
SOLUTION 1-2
-water 1
units ppm
pH 8.1
pe 8.451
density 1.023
temp 15
Ca 400
Mg 1350
Na 10500
K 380
Si 4.28
Cl 19000
Alkalinity 142 as HCO3
S(6) 2700
Fe 0.02
Mn 0.002
EQUILIBRIUM_PHASES 1
CO2(g) -3.5 1000000000 # mass percent, and quantity in mol
O2(g) -0.7 1000000000
# Equilibrate solutions with atmospheric gases
USE solution 1
USE equilibrium_phases 1
RUN_CELLS
-cells 1
SAVE solution 1
SAVE equilibrium_phases 1
END
USE solution 2
USE equilibrium_phases 1
RUN_CELLS
-cells 1
SAVE solution 2
SAVE equilibrium_phases 1
END
PHASES 1
Olivine
Mg1.86Fe0.14SiO4 + 4 H+ = 2 H2O + 1 SiO2 + 1.86 Mg++ + 0.14 Fe++
-log_k 26.448
-delta_h 10.710 kJ
RATES 1 #Data from Hermanska
Olivine
-start
10 FM_Oli = 149.53/1000 # Formula mass (149.53 g.mol, or XXXX g/mol)
20 mass = m0*FM_Oli # Initial mass of olivine,
30 d0 = 250/1000000 # Diameter, m (cannot be below 10um!)
40 p_Oli = 3300 # Particle density, kg/m3
50 SA_i = (3*mass)/(p_Oli*(d0/2)) # Initial total geometric surface area olivine, m2/g
60 d_min = 10/1000000 # Particle diameter minimum
70 SA_max = 3*mass/(p_Oli*(d_min/2)) # Maximum possible surface area due to attrition, m2
80 SA = (SA_max-(SA_max-SA_i))*((m/m0)^(2/3)) # Total current olivine surface area, m2
90 Aa = 1.48*10^5 # Acidic pre-exponential factor, mol/(m2.s)
100 Ab = 2.20*10^2 # Basic pre-exponential factor, mol/(m2.s)
110 Ea = -70400 # Acidic activation energy, J/mol
120 Eb = -60900 # Basic activation energy, J/mol
130 na = 0.4 # Acidic eaction order, -,
140 nb = 0.2 # Basic eaction order, -,
150 R = 8.314 # Universal gas constant, J/(mol*K)
160 rate = SA*(Aa*((ACT("H+"))^na)*EXP(Ea/(R*TK))+Ab*((ACT("H+"))^nb)*EXP(Eb/(R*TK))) # Rate expression, mol/s
170 moles = rate*TIME # Current number of moles of olivine, mol
180 SAVE moles # Save moles for PHREEQC to track reaction
-end
KINETICS 1
Olivine
-m0 1
-tol 1e-4
-runge_kutta 6
END
--- End code ---
After initialising these settings, I am trying to run code which mixes the solution with olivine (solution 1) using raw seawater (solution 2), following by performing the reaction (through run_cells). I've made a program that generates the iterations (repeats of solution_mix and run_cells) based on the time_steps and mixing fractions I desire. Here is an example of the first few steps:
--- Code: ---
USE solution 1
USE solution 2
USE equilibrium_phases 1
USE kinetics 1
SOLUTION_MIX 1
1 1.0 # Current solution
2 0.0 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.25919999999999993
END
SOLUTION_MIX 1
1 0.7941550190597205 # Current solution
2 0.2058449809402795 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.011089857429480876
END
SOLUTION_MIX 1
1 0.9890317785822964 # Current solution
2 0.010968221417703565 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.011564336356202853
END
SOLUTION_MIX 1
1 0.9885678686559283 # Current solution
2 0.011432131344071719 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.012059115837131429
END
--- End code ---
This proceeds until timesteps as large as 1e7 seconds (logarithmically spaced from time 0 until 120 months), every second we are mixing a constant amount of seawater. However, the problem is that despite diluting in increasing quantities, the concentrations and pH of the solution do not return to that of seawater! I believe there is an issue with how the solutions are being mixed, I've tried changing the pH of solution 2 to a value of 5 and this didn't alter my results.
dlparkhurst:
I'm not quite sure what you are doing. Here is a simplified version of your scripts that, I think, captures the sequence of the calculations.
I'm a bit confused by the tiny time steps, but I assume you have a sequence of longer and longer calculations.
Note that EQUILIBRIUM_PHASES 1 and KINETICS 1 definitions will be used in every RUN_CELLS calculation.
Now, if my revision is equivalent to yours, please be more specific about the progression of further calculations, and what exactly (pH, element concentration, ...) is not matching your expectation.
--- Code: ---TITLE Example 1. Olivine dissolution with kinetic rate law
PHASES
Olivine
Mg1.86Fe0.14SiO4 + 4 H+ = 2 H2O + 1 SiO2 + 1.86 Mg++ + 0.14 Fe++
-log_k 26.448
-delta_h 10.710 kJ
RATES #Data from Hermanska
Olivine
-start
10 FM_Oli = 149.53/1000 # Formula mass (149.53 g.mol, or XXXX g/mol)
20 mass = m0*FM_Oli # Initial mass of olivine,
30 d0 = 250/1000000 # Diameter, m (cannot be below 10um!)
40 p_Oli = 3300 # Particle density, kg/m3
50 SA_i = (3*mass)/(p_Oli*(d0/2)) # Initial total geometric surface area olivine, m2/g
60 d_min = 10/1000000 # Particle diameter minimum
70 SA_max = 3*mass/(p_Oli*(d_min/2)) # Maximum possible surface area due to attrition, m2
80 SA = (SA_max-(SA_max-SA_i))*((m/m0)^(2/3)) # Total current olivine surface area, m2
90 Aa = 1.48*10^5 # Acidic pre-exponential factor, mol/(m2.s)
100 Ab = 2.20*10^2 # Basic pre-exponential factor, mol/(m2.s)
110 Ea = -70400 # Acidic activation energy, J/mol
120 Eb = -60900 # Basic activation energy, J/mol
130 na = 0.4 # Acidic eaction order, -,
140 nb = 0.2 # Basic eaction order, -,
150 R = 8.314 # Universal gas constant, J/(mol*K)
160 rate = SA*(Aa*((ACT("H+"))^na)*EXP(Ea/(R*TK))+Ab*((ACT("H+"))^nb)*EXP(Eb/(R*TK))) # Rate expression, mol/s
170 moles = rate*TIME # Current number of moles of olivine, mol
180 SAVE moles # Save moles for PHREEQC to track reaction
-end
END
SOLUTION 1-2
-water 1
units ppm
pH 8.1
pe 8.451
density 1.023
temp 15
Ca 400
Mg 1350
Na 10500
K 380
Si 4.28
Cl 19000
Alkalinity 142 as HCO3
S(6) 2700
Fe 0.02
Mn 0.002
END
EQUILIBRIUM_PHASES 1
CO2(g) -3.5 100 # mass percent, and quantity in mol
O2(g) -0.7 100
END
USE solution 1
USE equilibrium_phases 1
SAVE solution 1-2
END
KINETICS 1
Olivine
-m0 1
-tol 1e-4
-runge_kutta 6
END
SOLUTION_MIX 1
1 1.0 # Current solution
2 0.0 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.25919999999999993
END
SOLUTION_MIX 1
1 0.7941550190597205 # Current solution
2 0.2058449809402795 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.011089857429480876
END
SOLUTION_MIX 1
1 0.9890317785822964 # Current solution
2 0.010968221417703565 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.011564336356202853
END
SOLUTION_MIX 1
1 0.9885678686559283 # Current solution
2 0.011432131344071719 # Seawater addition
END
RUN_CELLS
-cells 1
-time_step 0.012059115837131429
END
--- End code ---
achichat:
Hi dlpkarhust,
So what I'm attempting to do two things here: (1) dissolve olivine in seawater and (2) dilute the reacted solution in more seawater
The time steps were made smaller at the start so that the rates wouldn't fluctuate so much. And yes we do have a much longer sequence of calculations (repeats of the SOLUTION_MIX and RUN_CELLS)
So we do want to use everything labelled 1 in RUN_CELLS, and would avoid using solution 2
Your revision creates the same output as my original calculations.
My problem at the moment is the reaction will slow or stop at some point (once all the olivine is consumed/reacted), and the concentration of solution 1 should return to its initial specification since it is constantly being refreshed with seawater (diluted using solution 2). What I've found instead is that the pH of the solutions does not return to its initial specification, and you would expect it to do so if it is repeatedly mixed with solution 2.
Would a longer sequence help understand this outcome? I can provide the python code used to generate it.
dlparkhurst:
So, the pH does not return to 8.237, which is the solution 1-2 pH after reaction with O2 and CO2?
Maybe post the output from the final reaction (including the SOLUTION_MIX definition), and I'll take a look.
achichat:
That's right, it doesn't return to the original pH (obtained after equilibrating with air), this is true for up to 120 months.
I've posted the final output here, notice that the final pH is 8.981 despite mixing with a very large amount of solution 2 (99% seawater!). I'm struggling to attach the full input file I've used here, I could send a dropbox link with this.
--- Code: ------------------------------------------
Reading input data for simulation 1001.
---------------------------------------
SOLUTION_MIX 1
1 8.171066525397563e-06 # Current solution
2 0.9999918289334746 # Seawater addition
END
------------------
End of simulation.
------------------
---------------------------------------
Reading input data for simulation 1002.
---------------------------------------
RUN_CELLS
-cells 1
-time_step 12761815.36247772
END
--------------------------
Beginning of run as cells.
--------------------------
-----------------------------------------
Beginning of batch-reaction calculations.
-----------------------------------------
Reaction step 1.
Using solution 1.
Using pure phase assemblage 1. Pure-phase assemblage after simulation 1000.
Using kinetics 1.
Kinetics 1.
Time step: 1.27618e+07 seconds
Rate name Delta Moles Total Moles Reactant Coefficient
Olivine -5.984e-03 9.840e+00 Olivine 1
-------------------------------Phase assemblage--------------------------------
Moles in assemblage
Phase SI log IAP log K(T, P) Initial Final Delta
CO2(g) -3.50 -11.26 -7.76 1.000e+09 1.000e+09 -1.411e-02
O2(g) -0.70 -3.51 -2.81 1.000e+09 1.000e+09 -2.086e-04
-----------------------------Solution composition------------------------------
Elements Molality Moles
C 1.638e-02 1.637e-02
Ca 1.034e-02 1.034e-02
Cl 5.551e-01 5.551e-01
Fe 8.383e-04 8.381e-04
K 1.007e-02 1.007e-02
Mg 6.867e-02 6.866e-02
Mn 3.771e-08 3.770e-08
Na 4.731e-01 4.730e-01
S 2.912e-02 2.912e-02
Si 6.059e-03 6.058e-03
----------------------------Description of solution----------------------------
pH = 8.981 Charge balance
pe = 12.376 Adjusted to redox equilibrium
Activity of water = 0.981
Ionic strength (mol/kgw) = 6.521e-01
Mass of water (kg) = 9.999e-01
Total alkalinity (eq/kg) = 2.551e-02
Total CO2 (mol/kg) = 1.638e-02
Temperature (oC) = 15.00
Electrical balance (eq) = 3.122e-03
Percent error, 100*(Cat-|An|)/(Cat+|An|) = 0.27
Iterations = 59
Total H = 1.110529e+02
Total O = 5.570187e+01
--- End code ---
Navigation
[0] Message Index
[#] Next page
Go to full version