Understanding Beaker Reaction Results

[swhit]

I have built a model of an experiment I am going to run and hopefully will be able to calibrate a model to later.

Physical Experiment:

Start with acid mine drainage  in beaker, well agitated with atmospheric gas bubbled in (solution details in code)
Add increments of CaCO3 powder (e.g. 5g 4 times)
Log pH, ORP, Specific Conductace, Temperature, and Dissolved Oxygen as reaction progresses

I am confused about what happens each time I add Calcite to the beaker.  What I expect is the following:

- add solid calcite
- calcite starts to dissolve.  as calcite goes into solution:
- pH goes up
- Ferrihydrite ppts.
- Gypsum ppts. by common ion effect, lowering the sulfate concentration

What I see in the results:
- sharp increases in calcite and then slow declines as it dissolves (as expected)
- sharp decrease in Ferrihydrite and Gypsum solid phases upon addition of calcite (not expected, expected an increase here)
- slower return to previous amount of solid phase

I am wondering if the code is not representing my physical experiment like I think it is.  Any advice/thoughts/interpretations appreciated.

Plot of results attached.


Code:

CALCULATE_VALUES #calc variables inside a basic program to be used elsewhere
molCaCO3(mg)
-start
  10 mg_added = mg
  20 mol_added = mg_added/GFW("CaCO3")
  30 SAVE mol_added
-end

SC_T25
-start
  10 sc25 = SC/(1 + 0.021 * (TC - 25))
  20 SAVE sc25
-end

#roughly follows experimental steps for first beaker reaction test
#use solution from pit 4 add p200 alluvium (calcite)
#in future, update database with more of the elements or use another database (conductivity requires phreeqc.dat)
SOLUTION 1# Analytical Results for Pond 3, from Table 20 of EE/CA
pH 2.5 # Default: 7.
pe 4.0 #11.59 # Default: 4., 11.59 from tailings beaker test
redox pe
-units mg/L # Default: mmol/kgw
Al 100 mg/L
#As 13.0 mg/L
Ba 1.7 mg/L
#Be 0.018 mg/L
Cd 1.7 mg/L
Ca 500.0 mg/L
#Cr 3.5 mg/L
#Co 0.13 mg/L
Cu 8.9 mg/L
Fe 3100 mg/L
Mg 97 mg/L
#Hg .00047 mg/L
#Ni .3 mg/L
#V  .49 mg/L
Zn 610 mg/L

Na 0.0 mg/L# Chemical symbol from the 1st column in SOLUTION_MASTER_SPECIES, concentration, concentration is adapted to charge balance
Cl 5.8
F 12.0 mg/L
S(6) 15000 charge  # concentration is adapted to charge balance
C(4) 0.6 as HCO3 # Concentration is in mg of HCO3 = 0.6 / 0.61 = 0.98 mmol/L

Si 1.0 mg/L

EQUILIBRIUM_PHASES 1
#atmospheric gases
CO2(g) -3.5
O2(g) -0.67

#soil phases (update to match XRD)
Dolomite 0 0
Quartz 0 0
Chalcedony 0 0

#secondary phases
Fe(OH)3(a) 0 0
Gypsum 0 0

SAVE SOLUTION 1

END

SELECTED_OUTPUT # add this block where writing to output file should begin
-file D:\GoogleDrive\GradSchool\Research\Caselton\PHREEQC\BatchRxnSulfideTailingsAlluvium.csv # file name (plot in jupyter notebook)
-selected_out true
-user_punch true
-high_precision false
-reset false
-pH true
-pe true
-water true
-ionic_strength true
-time true
-solution true
-totals Ca Fe S(6)
-equilibrium_phases Fe(OH)3(a) Gypsum
-kinetic_reactants Calcite
-saturation_indices Fe(OH)3(a) Calcite Gypsum
-calculate_values SC_T25

END

USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1


KINETICS 1
Calcite
   -m       .05
   -m0      .05
   -parms   125000.0      30
   -tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS false

SAVE SOLUTION 2
END

USE SOLUTION 2
USE EQUILIBRIUM_PHASES 1
KINETICS 2
Calcite
   -m       .05
   -m0      .05
   -parms   125000.0      30
   -tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS false

SAVE SOLUTION 3
END

USE SOLUTION 3
USE EQUILIBRIUM_PHASES 1
KINETICS 3
Calcite
   -m       .05
   -m0      .05
   -parms   125000.0      30
   -tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS false

SAVE SOLUTION 4
END

USE SOLUTION 4
USE EQUILIBRIUM_PHASES 1
KINETICS 4
Calcite
   -m       .05
   -m0      .05
   -parms   125000.0      30
   -tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS false

SAVE SOLUTION 5
END

RATES
Calcite
-start
1   rem M = current number of moles of calcite
2   rem M0 = number of moles of calcite initially present
3   rem PARM(1) = A/V, cm^2/L
4   rem PARM(2) = exponent for M/M0
10  si_cc = SI("Calcite")
20  if (M <= 0 and si_cc < 0) then goto 200
30    k1 = 10^(0.198 - 444.0 / TK )
40    k2 = 10^(2.84 - 2177.0 / TK)
50    if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
60    if TC > 25 then k3  = 10^(-1.1 - 1737.0 / TK )
70    t = 1
80    if M0 > 0 then t = M/M0
90    if t = 0 then t = 1
100   area = PARM(1) * (t)^PARM(2)
110   rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120   rem 1e-3 converts mmol to mol
130   rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
140   moles = rate * TIME
200 SAVE moles
-end

[dlparkhurst]

I don't really want to sort through all of your code. Try something like this. RUNCELLS "USEs" and "SAVEs" solutions, equilibrium_phases, and other reactants (KINETICS is different in that the number of moles, M, is continuously updated in all simulations.)

Your number for parameter 2 in KINETICS should be around 1.

Your RATES block needs to be defined before KINETICS. Here I use the default Calcite rate from phreeqc.dat.

Code: [Select]

SOLUTION 1
    temp      25
    pH        2.5
    pe        4
    redox     pe
    units     mg/l
    density   1
    C(4)      0.6 as HCO3
    Ca        500
    Cl        5.8
    Fe        3100
    Mg        97
    Na        0
    S(6)      15000 charge
    Si        1
    Zn        610
    -water    1 # kg

EQUILIBRIUM_PHASES 1
#atmospheric gases
CO2(g) -3.5
O2(g) -0.67

#soil phases (update to match XRD)
Dolomite 0 0
Quartz 0 0
Chalcedony 0 0

#secondary phases
Fe(OH)3(a) 0 0
Gypsum 0 0

SAVE SOLUTION 1

END
USER_GRAPH 1
10 GRAPH_X (total_time)
20 GRAPH_Y SI("Calcite")
30 GRAPH_SY -LA("H+")
END

KINETICS 1
Calcite
   -m       0.05
   -m0      0.05
   -parms   125000.0    0.6 #  30
   -tol     1.e-8
-steps 10000 in 10 steps #seconds
END
RUN_CELLS
-cell 1
-initial_time 0
END

KINETICS 1
Calcite
   -m       0.05
   -m0      0.05
   -parms   125000.0    0.6 #  30
   -tol     1.e-8
-steps 10000 in 20 steps #seconds
RUN_CELLS
-cell 1
-initial_time 10000
END


[swhit]

Thanks for the suggestions David. I like the more succinct code you shared, and it seems to produce a geochemical behavior that is more in line with what I expected.  I also figured out why mine wasn't working as expected -- I wasn't saving and updating EQUILIBRIUM_PHASES as I was with SOLUTION, so with each SAVE/USE to start a new addition, I was starting with the original EQUILIBIRIUM_PHASES.

A couple questions/observations:

- I think you meant to type "-start_time" instead of "-initial_time" in the RUN_CELLS code block
- Even when -initial_time is changed to -start_time each RUN_CELLS call still starts at time zero when the results are output to a .csv via SELECTED_OUTPUT (though the results seem to plot correctly on USER_GRAPH). I am addressing this with post processing and the "sim" number but it would good to know if there is a way to change this
- Does the "END" at the end of the RUN_CELLS block do anything, since the RUN_CELLS ends itself?  I took it out and it seems to not make a difference

- I think that in the real experiment, the dissolution rate and total amount of dissoltion per addition may have been reduced by precipitation of Ferrihydrite on the surface of the calcite particles.  I think I have a decent handle on how to integrate the kinetics for Ferrihydrite precipitation (and Fe2+ oxidation), but am wondering if you have any suggestions for how to incorporate precipiation of the Ferrihydrite coating on the calcite grains.

Thanks!

Example of original code, modified to work correctly for anyone to look at later (David's way is much cleaner but this works too):

Code: [Select]

SOLUTION 1# Analytical Results for Pond 3, from Table 20 of EE/CA
pH 2.5 # Default: 7.
pe 4.0 #11.59 # Default: 4., 11.59 from tailings beaker test
redox pe
-units mg/L # Default: mmol/kgw
Al 100 mg/L
#As 13.0 mg/L
Ba 1.7 mg/L
#Be 0.018 mg/L
Cd 1.7 mg/L
Ca 500.0 mg/L
#Cr 3.5 mg/L
#Co 0.13 mg/L
Cu 8.9 mg/L
Fe 3100 mg/L
# Pb non-detect, confirmed in EECA text, but look at new data
Mg 97 mg/L
Mn 1200 mg/L   #should this be ug/L ??? listed as mg/L in EECA
#Hg .00047 mg/L
#Ni .3 mg/L
#V  .49 mg/L
Zn 610 mg/L

Na 0.0 mg/L# Chemical symbol from the 1st column in SOLUTION_MASTER_SPECIES, concentration, concentration is adapted to charge balance
Cl 5.8
F 12.0 mg/L
S(6) 15000 charge  # concentration is adapted to charge balance
C(4) 0.6 as HCO3 # Concentration is in mg of HCO3 = 0.6 / 0.61 = 0.98 mmol/L

Si 1.0 mg/L

EQUILIBRIUM_PHASES 1 # cell number or a range of cells. Default: 1.
#atmospheric gases
CO2(g) -3.5
O2(g) -0.67

#soil phases (update to match XRD)
#Calcite 0 0 # name (must be defined in PHASES), Saturation Index (Default: 0), initial amount (moles, default: 10 moles)
Dolomite 0 0 #diss # Dolomite can dissolve only
Quartz 0 0 #precipitate # Quartz can precipitate only if uncommented
Chalcedony 0 0

#secondary phases
Fe(OH)3(a) 0 0
Gypsum 0 0
#Goethite  0 0

SAVE SOLUTION 1

END

SELECTED_OUTPUT # add this block where writing to output file should begin
-file D:\GoogleDrive\GradSchool\Research\Caselton\PHREEQC\BatchRxnSulfideTailingsAlluvium.csv # file name (plot in jupyter notebook)
-selected_out true
-user_punch true
-high_precision false
-reset false
-pH true
-pe true
-water true
-ionic_strength true
-time true
-solution true
-totals Ca Fe S(6)
-equilibrium_phases Fe(OH)3(a) Gypsum
-kinetic_reactants Calcite
-saturation_indices Fe(OH)3(a) Calcite Gypsum
-calculate_values SC_T25

RATES
Calcite
-start
1   rem M = current number of moles of calcite
2   rem M0 = number of moles of calcite initially present
3   rem PARM(1) = A/V, cm^2/L
4   rem PARM(2) = exponent for M/M0
10  si_cc = SI("Calcite")
20  if (M <= 0 and si_cc < 0) then goto 200
30    k1 = 10^(0.198 - 444.0 / TK )
40    k2 = 10^(2.84 - 2177.0 / TK)
50    if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
60    if TC > 25 then k3  = 10^(-1.1 - 1737.0 / TK )
70    t = 1
80    if M0 > 0 then t = M/M0
90    if t = 0 then t = 1
100   area = PARM(1) * (t)^PARM(2)
110   rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120   rem 1e-3 converts mmol to mol
130   rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
140   moles = rate * TIME
200 SAVE moles
-end

USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1


KINETICS 1
Calcite
-m       .02
-m0      .02
-parms   125000.0      0.6
-tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS true

SAVE SOLUTION 2
SAVE EQUILIBRIUM_PHASES 2
END

USE SOLUTION 2
USE EQUILIBRIUM_PHASES 2
KINETICS 1
Calcite
-m       .02
-m0      .02
-parms   125000.0      0.6
-tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS true

SAVE SOLUTION 3
SAVE EQUILIBRIUM_PHASES 3
END

USE SOLUTION 3
USE EQUILIBRIUM_PHASES 3
KINETICS 1
Calcite
-m      .02
-m0      .02
-parms   125000.0      0.6
-tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS true

SAVE SOLUTION 4
SAVE EQUILIBRIUM_PHASES 4
END

USE SOLUTION 4
USE EQUILIBRIUM_PHASES 4
KINETICS 1
Calcite
-m       .02
-m0      .02
-parms   125000.0      0.6
-tol     1.e-8
-steps 10000 in 100 steps #seconds
INCREMENTAL_REACTIONS true

SAVE SOLUTION 5
SAVE EQUILIBRIUM_PHASES 5
END


-

[dlparkhurst]

(1) Yes, should be -start_time.

(2) The Basic function TOTAL_TIME should give the time accounting for the start_time. You can use it in USER_PUNCH.

SELECTED_OUTPUT 1
...
USER_PUNCH 1
-heading total_time_sec
10 PUNCH TOTAL_TIME

(3) END is implied at the end of the input file. It would be required if you continue with more keywords and calculations.