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: [Select]

        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
       

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: [Select]

       
                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

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: [Select]

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


[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: [Select]

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


[dlparkhurst]

You are still dissolving olivine to the tune of 6 millimoles in this calculation. There are still 9.84 moles of olivine available to react. You won't return to initial conditions until the olivine is depleted.

[achichat]

Hi dlparkhurst,

Thank you for pointing this out, there is in fact an issue with the specified conditions, I have reduced the quantity of solid and the particle size to verify the behaviour of the code. We did eventually manage to get the result we expected (peak in pH followed by return to pH = 8.237), however, this required us to also change how we specified the time_steps of RUN_CELLS. Initially I presumed that it was meant to show the increment in time for each calculation but it only worked when I specified the point in time the calculation would be made for (and the latter made it work properly).

Any chance you'd be able to clarify this? I'm still doubting whether my calculations are integrating from 0 to the specified time_step or if it is doing it in an incremental fashion (meaning over the time increment of the next RUN_CELLS). Note that I have removed any specification for INCREMENTAL_REACTIONS in my input file.

Code: [Select]


        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 1 # mass percent, and quantity in mol
        O2(g)      -0.7 1
       

        PHASES
        Olivine
        Mg1.86Fe0.14SiO4 + 4 H+ = 2 H2O + 1 SiO2 + 1.86 Mg++ + 0.14 Fe++
        -log_k 26.448
        -delta_h 10.71 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 = 20/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
       
        USE solution 1
        USE equilibrium_phases 1
        SAVE solution 1-2
        END
       
        KINETICS 1
        Olivine
        -m0       0.1
        -tol 1e-4
        -runge_kutta    6     
       
        INCREMENTAL_REACTIONS False
       
        SELECTED_OUTPUT
            -reset false
            -selected_out true
            -high_precision true
            -time true
            -pH true
            -alkalinity true
            -totals C(4) Ca Mg Si Fe
            -molalities C_di HCO3- CO3-2
            -saturation_indices Calcite Aragonite
            -kinetic_reactants Olivine
            -equilibrium_phases Olivine CO2(g)

  USER_GRAPH 1

-headings conversion pH
-axis_scale x_axis auto
-axis_scale y_axis auto
-axis_scale sy_axis auto
-axis_titles "Time" "Olivine conversion, %" "pH"
-start
10 PLOT_XY time, ((m0-m)/m0)*100, y-axis = 1
20 PLOT_XY time, -la("H+"), y-axis = 2
-end
           
        END
           
                         
                USE solution 1
                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.9974147010947624  # Current solution
                    2 0.0025852989052376233    # Seawater addition
                    END
                   
                    RUN_CELLS
                    -cells 1
                    -time_step 0.2702898574294808
                    END
                   
                    SOLUTION_MIX 1
                    1 0.9998891137228353  # Current solution
                    2 0.00011088627716476233    # Seawater addition
                    END
                   
                    RUN_CELLS
                    -cells 1
                    -time_step 0.28185419378568366
                    END
                   
                    SOLUTION_MIX 1
                    1 0.9998843700082791  # Current solution
                    2 0.00011562999172085728    # Seawater addition
                    END
                   
                    RUN_CELLS
                    -cells 1
                    -time_step 0.2939133096228151
                    END


Final output (After multiple iterations) yields the desired pH. In the previous examples I provided, the time_steps were the increments, notice how the first time_step is larger than the second time_step.

[dlparkhurst]

You need to use the Basic function TOTAL_TIME, not TIME. Maybe not my best moment, but TIME is an internally derived sub time step, whereas TOTAL_TIME is the elapsed time from the beginning of the simulations.