Kinetics of anhydrite dissolution

[piet]

Hi everyone,

I am having trouble simulating the dissolution kinetics of anhydrite in water. The problem is relatively straightforward. I am interested in knowing how much sulphate can be released from volcanic ash upon exposure to water (solid:water ratio 1 g:100 ml) over a period of 3600 s (1 h). I know the ash material contains about 10 wt.% anhydrite, roughly corresponding to 0.01 mol anhydrite per kgw. I am using the RATES block for anhydrite as defined by Zhang et al. (2019). the script appears below.

i)   PHREEQC returns the following warning message:

WARNING: Negative moles in solution 1 for Ca, -1.399155e-02. Recovering...
WARNING: Negative moles in solution 1 for Ca, -7.645802e-03. Recovering...

ii)   The same total number of moles of anhydrite is dissolved at each time step, i.e. there does not seem to be a kinetically-controlled reaction going on. Anhydrite is predicted to remain unsaturated. Increasing the initial mole of anhydrite (m0) in the KINETICS block does not solve the problem.

Many thanks for your time.  Pierre


PHASES
Anhydrite
Ca1S1O4 = 1.0000 SO4-2 + 1.0000 Ca+2
   -analytic   9.835770e+003   2.965925e+000   -4.073817e+005   -3.829436e+003   1.835113e+007   -9.989717e-004
   -Vm   45.9400

RATES
Anhydrite
# from Palandri and Kharaka 2004
# experimental condition range T=25-50C, pH=~7

-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=0      
12 E1=0
13 n1=0
20 rem neutral solution parameters
21 a2=2.05E-01   
22 E2=14299
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("anhydrite")
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

END

SOLUTION 1
pH 7 charge
temp 20
C 1 CO2(g) -3.5
   
END

USE solution 1

INCREMENTAL_REACTIONS true

KINETICS 1
Anhydrite
   -formula   Ca1S1O4 1.0 
   -m0      1e-3    # initial moles of anhydrite per kg of water
   -parms   10   1 #total surface area per kg of water (m2/kgw) and the scaling factor
   -tol      1e-8
   -steps      900 1800 2700 3600 s
   

USER_GRAPH 1
    -headings               Time TOT(S) SI(Anhydrite)
    -axis_titles            "Time, seconds" "Molality" "Saturation index"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y TOT("S")
30 GRAPH_SY SI("Anhydrite")
  -end
    -active                 true
END

[dlparkhurst]

Two things. The given rate is very fast (equilibrium in a few seconds), and the solubility (~0.02 mol/kgw) is more than the specified amount of anhydrite available (0.001 mol/kgw).

Using an exponential distribution of time for the kinetics integration and a greater amount of anhydrite shows the progression to equilibrium. You can consider whether the very fast rate is appropriate or not.

Code: [Select]

PHASES
Anhydrite
Ca1S1O4 = 1.0000 SO4-2 + 1.0000 Ca+2
   -analytic   9.835770e+003   2.965925e+000   -4.073817e+005   -3.829436e+003   1.835113e+007   -9.989717e-004
   -Vm   45.9400

RATES
Anhydrite
# from Palandri and Kharaka 2004
# experimental condition range T=25-50C, pH=~7

-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=0     
12 E1=0
13 n1=0
20 rem neutral solution parameters
21 a2=2.05E-01   
22 E2=14299
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("anhydrite")
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

END

SOLUTION 1
pH 7 charge
temp 20
C 1 CO2(g) -3.5
   
END

USE solution 1

#INCREMENTAL_REACTIONS true

KINETICS 1
Anhydrite
   -formula   Ca1S1O4 1.0
   -m0      1e-1 #1e-3   # initial moles of anhydrite per kg of water
   -parms   10   1 #total surface area per kg of water (m2/kgw) and the scaling factor
   -tol      1e-8
   -steps    1e-4 1e-3 1e-2 1e-1 1 10 100 1000 # 900 1800 2700 3600 s
   

USER_GRAPH 1
    -headings               Time TOT(S) SI(Anhydrite)
    -axis_titles            "Time, seconds" "Molality" "Saturation index"
    -axis_scale x_axis      auto auto auto auto log
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y TOT("S")
30 GRAPH_SY SI("Anhydrite")
  -end
    -active                 true
END

[piet]

Thank you very much for your helpful and swift reply!