Dissolution of phosphates in a granite

[RadiantRocks23]

Dear Mr Parkhurst, dear Phreeqc users,

I?m very new to PHREEQ. I tried to learn with the forum and with the documentation.

I am trying to model the dissolution of phosphates (such as fluorapatite) in a material with a granitic composition assuming equilibrium. My SOLUTION block is representative of the experimental conditions. I added the fluorapatite PHASE block from Zhang et al., 2019 Computer & Geoscience. In the EQUILIBRIUM_PHASES block, I assume that the material is made of 4 minerals. I calculated the number of moles based on the mineral proportions.

Code: [Select]

# Title Apatite Dissolution
SOLUTION 1
    temp      25
    pH        0.5
    pe        4
    redox     pe
    units     mol/l
    density   1
    N(5)      0.31
    -water    0.01 # kg
 
PHASES 1
Fluorapatite
Ca5P3O12F1 + 3.0000 H+ = 3.0000 HPO4-2 + 1.0000 F- + 5.0000 Ca+2
      -analytic   2.906406e+004     8.907085e+000     -1.192683e+006     -1.134608e+004    5.304647e+007     -3.060213e-003
      -Vm   157.5600
 
EQUILIBRIUM_PHASES 1
    Fluorapatite 0 0.0003
    K-mica    0 0.0007
    Quartz    0 0.003
    Albite 0 0.0007

In my experiments, fluorapatite has almost completely dissolved. However, in the ? Phase Assemblage ? output, albite has completely dissolved while most fluorapatite remains, which is not consistent with my experimental observations. I therefore supposed that I should include kinetics effects to my input and study the dissolution in function of time. In order to do so, I used the RATES block of Fluorapatite from Zhang et al., 2019 Computer & Geoscience.

Code: [Select]

# Title Apatite Dissolution
SOLUTION 1
    temp      25
    pH        0.5
    pe        4
    redox     pe
    units     mol/l
    density   1
    N(5)      0.31
    -water    0.01 # kg
 
PHASES 1
Fluorapatite
Ca5P3O12F1 + 3.0000 H+ = 3.0000 HPO4-2 + 1.0000 F- + 5.0000 Ca+2
      -analytic   2.906406e+004     8.907085e+000     -1.192683e+006     -1.134608e+004    5.304647e+007     -3.060213e-003
      -Vm   157.5600
 
RATES 1
Fluorapatite
# from Palandri and Kharaka 2004
# experimental condition range T=25C, pH=2.2-2.8
# uncertainty for Ea is large
# Hydroxyapatite activation energy is assumed
 
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=1.18E+40
12 E1=250000
13 n1=0.613
20 rem neutral solution parameters
21 a2=6.32E+35
22 E2=250000
30 rem base solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Fluorapatite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
 
EQUILIBRIUM_PHASES 1
    K-mica    0 0.0007
    Quartz    0 0.003
    Albite 0 0.0007
 
KINETICS 1
Fluorapatite
    -formula  Fluorapatite
    -m        0.04
    -parms    50 1
-steps       100000 in 100 steps # seconds
INCREMENTAL_REACTIONS True

Simulation with the latter input takes a while and then shows the following error: "Bad RK steps >500 cells in 1. Please decrease timestep.". However, I don't see why I would need to decrease the timestep. Are there any ways I could improve this code ?

Thanks a lot, I hope it will help future users as well.

[dlparkhurst]

You can solve your immediate problem by using -cvode in KINETICS. Cvode is an implicit ordinary differential equation integration scheme that can be more robust for some problems than the default Runge Kutta scheme.

However, albite still dissolves completely. You might want to start with only fluorapatite as an equilibrium reaction, and then as a kinetic reaction to see whether that reaction alone is sufficient to match your data. Next you could add additional reactions as either equilibrium or kinetics to refine your results.

If you have the evolution of the solution composition, you may want to consider INVERSE_MODELING. The exercise would be to determine a set or sets of mineral reactions that can account for the water chemistry evolution. It could be a faster way to quantify the reactions in your system, but it can be harder than it seems to get plausible reactions. Still, I think that the exercise is worthwhile because it makes you consider mineral and gas sources and sinks for the element concentrations that you are trying to match, and if successful, it matches your solution exactly.

[RadiantRocks23]

Dear David,

Thank you for your response. I paused my simulations because I needed more data to constrain my system. I now have the necessary results.

I tried running inverse geochemical modeling as you suggested but encountered some issues. I will post a new thread in the appropriate section to address this.

Best,