kinetics for PFAs sorption

[oatteia]

Hello,
for PFAs, there is sorption both on the solid phase and on the air-water interface in the soil. This sorption depends on the water saturation, and can be kinetic. So in order to do that I made first two rates Pfas_aw and Pfas_s, using immobile dissolved species. It works. But it is quite slow, I thought that it was because the Pfas species were at disequilibrium when one rate is calculated and then tries to equilibrate ith the other rate.
Therefore I tried to use one rate (for dissolved Pfas) writitng the equations and solving them first (on paper), it gives a second order polynomial. When I use this (enclosed) it is strange, at the frist tep it gives the correct numbers, with Pfas_aw amount calculated as positive, but when it is used iwth the get, it seems that is becomes negative and then phreeqc spend time to try to recover.
I someone can help
attachements do not work so i copy the files below

intichem1.pqi
DATABASE D:\iPht3d\Exv2\PFAs\wallisPFAs1\phreeqc.dat
Solution 1
 units mol/kgw
Cl 1e-4
Na 1e-4
Pfas 1.6e-5
pH 7
pe 4
Pfas_aw 0.0
Pfas_s 0.0

Kinetics 1-70

Pfas
 -parms 10 0.059 0.12 90
-formula Pfas -1

Pfas_aw
-formula Pfas_aw 1

Pfas_s
-formula Pfas_s 1

-tol 1e-8

############## Phreeqc.dat
SOLUTION_MASTER_SPECIES
H        H+             -1.     H               1.008
H(1)     H+             -1.     H
H(0)    H2        0      H
E        e-             0.0     0.0             0.0
O        H2O            0.0     O               16.00
O(0)     O2             0.0     O
O(-2)    H2O            0.0     H2O
Na       Na+            0.0     Na              22.9898
Cl       Cl-            0.0     Cl              35.453
N     NO3-           0.0     NO3               14.
N(5)     NO3-           0.0     NO3      56.       
N(0)     N2             0.0     N2        28.     
S        SO4-2          0.0     SO4             32.064
S(+6)     SO4-2          0.0     SO4   96.
S(-2)    HS-            1.0     HS     34.
C        CO3-2          2.0     HCO3            12.0111
C(+4)    CO3-2          2.0     HCO3 
Tracer Tracer  0  Tracer  1
Pfas_s Pfas_s 0 Pfas_s 1
Pfas_aw Pfas_aw 0 Pfas_aw 1
Pfas Pfas 0 Pfas 1

SOLUTION_SPECIES
Tracer = Tracer ; log_k 0.0
Pfas_s = Pfas_s; log_k 0.0
Pfas_aw = Pfas_aw; log_k 0.0
Pfas = Pfas ; log_k 0.0

H+ = H+
 log_k  0.000
 -gamma 9.0000    0.0000
e- = e-
 log_k   0.000
H2O = H2O
 log_k  0.000
Na+ = Na+
 log_k  0.000
Cl- = Cl-
 log_k  0.000
NO3- = NO3-
 log_k  0.000
CO3-2 = CO3-2
 log_k  0.000
SO4-2 = SO4-2
 log_k  0.000
H2O = OH- + H+
 log_k   -14.000
2 H2O = O2 + 4 H+ + 4 e-
 log_k  -86.08
2 H+ + 2 e- = H2
 log_k     -6.
Na+ + H2O = NaOH + H+
 log_k    -14.180
CO3-2 + H+ = HCO3-
 log_k   10.329
CO3-2 + 2 H+ = CO2 + H2O
 log_k  16.681
SO4-2 + H+ = HSO4-
 log_k   1.988
HS- = S-2 + H+
 log_k  -12.918
SO4-2 + 9 H+ + 8 e- = HS- + 4 H2O
 log_k  33.65
HS- + H+ = H2S
 log_k  6.994
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
 log_k  207.080

RATES
Pfas
-start
1   # Pfas sorb,
2   k= parm(1)
3   a = parm(2)
4   b = parm(3)  # given in cm2/cm3
5   Koc = parm(4)
6   R = 8.314
7   T = 298.15
8   rhob = 1.6  # kg/dm3
9   rhos = 2.65  # solid kg/dm3
10   dgrain = 0.13416  # mm geom mean calc from paper 1/2 clay 1/2 sand
12   moles=0
14   C = Tot("Pfas")  # mol/L = mmol/cm3
16   if C < 1e-12 then goto 200
18     # print " C " C
20  toto=0 # to use phreeqc (0) and rm (1)
22  if toto<1 then poro = 0.4 else poro = CALLBACK(CELL_NO, 0, "PORO")
24  if toto<1 then sw=0.6 else sw = CALLBACK(CELL_NO, 0, "WSAT")
26  if toto<1 then foc=0.01 else sw = CALLBACK(CELL_NO, 0, "FOC")
28   #if sw=0 then goto 200
32   Aw=600*rhob/(dgrain*rhos)*(1-sw)  # in dm2/dm3
34   thetaw = sw*poro
36      #print "poro " poro thetaw foc 1/thetaw
40   Gam = 0.0757-0.0001515*(T-273.15) #here T in ?C ! (don't use Gamma, it is reserved for phreeqc)
42   sumC = TOT("Pfas_s")+TOT("Pfas_aw")+C
44      print "totals diss sorb aw " C TOT("Pfas_s") TOT("Pfas_aw")
46   B0 = thetaw/Aw*R*T/(Gam*b*10)
47      ### Kaw = Gam/(R*T)*b/(a+C*1000)*10 #fact 1000 : C-> umol/cm3, fact 10: cm3/cm2 -> dm3/dm2
48      ### Caw = TOT("Pfas_aw")*thetaw # Pfas_aw is immobile but expressed in the water phase: mol/L*thet-> (mol/dm3)
60   if foc>0 then B1 = thetaw/(rhob*Koc*foc) else B1 = 1
64   if foc>0 then e0 = -(1+1/B1)*1000 else e0 = -1000
65      print "gam a b Aw B1" Gam a b Aw B1
66   if foc>0 then e1 = sumC*1000-a*(1+1/B1)-1/B0 else e1 = sumC*1000-a-1/B0
68   delta = e1^2 - 4*e0*a*sumC
70   Ceq = (- e1 - sqrt(delta))/(2*e0) ###Ceq = a*B0/(1-B0*1000) #### Ceq = Caw/(Kaw*Aw)
71   Kaw = Gam/(R*T)*b/(a+Ceq*1000)*10
72      print " Kaw Caw C Ceq " Kaw Caw C Ceq
100  if C>=Ceq then rate = 1-exp(C-Ceq) else rate=exp(C-Ceq)-1  #rate * time
102  #moles = rate * time/24/3600
106     print "Pfas rate mol " rate moles # dont'know where it writes
110  if moles>C then moles=C*0.999
120  if (C + moles) < 0 then moles = -C*0.999
130  Cnow = C-moles
134  if foc>0 then moles_s = Cnow/B1-Tot("Pfas_s") else moles_s = 0
136  put(moles_s,1)  # for sorbed
144  moles_aw = Cnow*(Kaw*Aw)/thetaw-Tot("Pfas_aw")
147     print "cnow sorb aw" Cnow moles_aw moles_s
150  put(moles_aw,2)  # for aw
200  save moles
-end

################
Pfas_s
-start
10   rate = 1
20   moles = get(1)
22   print "in pfaa_s mole sorb " moles
30   save moles
-end

################
Pfas_aw
-start
20   moles = get(2)
22   print "in pfa_aw mole aw " moles
30   save moles
-end

[dlparkhurst]

Sorry, I can't debug something that is hooked to your transport model. If you can make a case that runs with PHREEQC, I might be able to help.

I don't think that you should have aqueous species for _aw and _s if they are sorbed. It seems like you would only need two reactions that represent the transfer to or from the sorbed pools.  The signs are tricky, but SAVE moles should be negative to increase a pool of sorbed Pfas; consequently, you should use Pfas +1 in  the formula so that the product of moles*coefficient is negative, which removes Pfas from solution.

You might also need to use -cvode, but I think the first problem is your KINETICS and RATES definitions.

[oatteia]

Thank you,
I think I will try that way. I use immobile species (like we often do for BTEX) as the sorption was on the air-water interface, which is not a solid. The other advantage of using a dissolved but immobile species is that my transport code can have access to it directly (and not through selected output which I found more complex...and may be slower no?)
olivier

[oatteia]

Me again, sorry,
There is something that I don't understand. As I made the analytical calculation, I have calculated the formula to get directly the correct moles for each species (solving my 2nd order polynomial), so ideally there should only be one kinetic step. But this doest not work:
in the first step (results given by the print in phreeqc)
the entering concentrations are :
totals dissolved   1.6000e-05     sorbed 0      aw  0
internal calculations (they are ok gam a b Aw B1   7.1913e-02   5.9000e-02   1.2000e-01   1.0801e+03   1.6667e-01 )
here C is initial, and Ceq is the conc to be reached:
Kaw 5.7381e-04            C_aw 0   C 1.6000e-05   Ceq 1.6697e-06
Moles transfer (from Pfas species)
   1.4330e-05 this is ok
Then I use the iniital C and the moles to get the resulting conc (Cnow) and calculate the sorbed and aw:
cnow sorb aw   1.6698e-06   4.3121e-06   1.0019e-05
this is correct, and thus the following mass transfer are also ok
in pfa_aw mole aw    4.3121e-06
in pfaa_s mole sorb    1.0019e-05

This is the next step, and my problem is here although the moles from above are transfered through the put and get, the above amounts are not really transfered, but it seems that only 20% of this amount is really transfered
totals dissolved 1.3134e-05 (instead of 1.67e-6) 2.0038e-06 (instead of 1e-5)   8.6243e-07 (instead of 4.31e-6)
I can understand that for convergence purpose the kinetics can be done by splitting the steps, but as when I solved my problem the equaiton is non linear, intermediate steps may not be correct.
so is there a way to tell phreeqc to stop at the first kinetic step?
sincerely
Olivier

[dlparkhurst]

Unless I can see a pure PHREEQC calculation, I don't think I can help. I simply don't understand your calculation well enough. Which of the pools of PFAs are transported? Tony has done a lot with SURFACES that can be mobile; perhaps you can treat some of your pools as surfaces, but the transport would be complicated. It may be that you are simply trying to do something for which PhreeqcRM is not suited.

The rate equations are integrated, either by Runge-Kutta or CVODE in multiple sub time steps. There is no way to do just one complete step. The sequence of times by which the rates are evaluated are not monotonic and the final evaluation may not be at the end of the time step. Rather, a calculation is made at the end of the sequence with the multiple evaluations to arrive at an overall integrated rate for the time step. So, I think you are telling me that using GET and PUT within RATES is not working as you hoped, but the integration methods have no options for alternative schemes.

If you make a PHREEQC calculation, I will look at it.

[oatteia]

Hi David,
thanks for your continuous help. I understand where the problem is : the rate and moles are calculated correclty for the calculation that I want to do, but as phreeqc take this moles amount and divide it for the purpose of reaching the correct rate. the problem for me is that when phreeqc divides, it divides all moles for all species by the same factor, which is not consistent with the double equilibrium that I need.
My problem is that I need this to be kinetic o I cannot really use surfacs may be kinetic minerals but in that case it won't change the problem I think.
so i need to think in a different way
by the way the files I sent are phreeqc ones
sincerely
Olivier

below is the solution that is correct and works (and need to be fastened as it is slowed because it produces negative values for aw sorbed species)

Pfas_s
-start
1   # Pfas sorbed,
3   k = parm(1) # kinetic transfer rate
4   Koc = parm(2)  # cm3/g, or dm3/kg
8   rhob=1.6 #g/cm3 or kg/dm3
10  moles=0
12  C = TOT("Pfas")
14  if C<1e-12 then goto 200
20  poro = 0.4
22  sw=0.6
24  foc= 0.01
32  thetaw = sw*poro
40  Cs = TOT("Pfas_s")*thetaw/rhob
48   Kd=Koc*foc
52   Ceq = Cs/Kd  # Kd=Cs/Cliq
60   rate = k*(C-Ceq)/24/3600
100  moles = rate * time
110  if moles>C then moles=C*0.999
200  save moles
-end

##############

Pfas_aw
-start
1   # Pfas sorb,
3   k= parm(1)
4   a = parm(2)
5   b = parm(3)  # given in cm2/cm3
6   R = 8.31
7   T = 298.15
8   rhob = 1.6  # kg/dm3
9   rhos = 2.65  # solid kg/dm3
10   dgrain = 0.13  # mm geom mean calc from paper 1/2 clay 1/2 sand
12   moles=0
14   C = Tot("Pfas")  # mol/L = mmol/cm3
16   if C < 1e-12 then goto 200
18   print " C " C
22  poro = 0.4
24  sw=0.6
33   Aw=600*rhob/(dgrain*rhos)*(1-sw)  # in dm2/dm3
34   thetaw = sw*poro
36   print "poro " 1/thetaw
40   Gam=0.0757-0.0001515*(T-273.15)
42   sumC = TOT("Pfas_aw")+C
46   B0 = thetaw/Aw*R*T/(Gam*b*10)
52   delta = (1000*sumC-a-1/B0)^2 + 4*1000*a*sumC
54   Ceq = (a+1/B0-1000*sumC - sqrt(delta))/(-2*1000)
55   print "Ceq " Gam Kaw Caw C Ceq
60   rate = k*(C-Ceq)/24/3600
100  moles = rate * time
102  print "Pfas rate " rate
110  if moles>TOT("Pfas") then moles=C*0.999
120  if (C + moles) < 0 then moles = -C*0.999
200  save moles
-end

[dlparkhurst]

I would need an entire functioning calculation in PHREEQC to help you, including SOLUTION, RATES, and KINETICS blocks and either one or two batch reactions that demonstrate the problem or an ADVECTION calculation (which would be simpler than TRANSPORT).

I'm not sure what you mean about "divides all moles for all species". I thought the -formula in the KINETICS definitions would only contain "Pfas", so there would not be multiple elements involved (KINETICS only deals with changes in element moles in solution), only changes in the Pfas concentration.