rate for surface complexation

[hamiderfani]

Hello everybody
hope you are doing very well and thanks for your reply in advance.
I have a question regarding the surface complexation in Phreeqc.
I have a SCM and it works very well. I wanted to know if there is any way to assign a rate equation (something like a langmuir type rate equation) to reactions in SCM?
In reality the response is never instantaneous on the surface so in this way I hope to be able to impose some delay to the effect on surface properties.
Any other way to impose a delay on the surface properties is also welcome.

Thanks

[dlparkhurst]

Yes, you can use KINETICS to simulate kinetic sorption, although it is not very easy.

Attached is a file that kinetically sorbs a constituent "A" using a Langmuir isotherm. The idea is to calculate what surface sorption would be given the dissolved concentration and use the difference between that and what is actually sorbed to formulate a rate equation.

Also attached is an example of kinetic sorption for a surface complexation example. It is not my example, and it has been a long time since I looked at it, so you are on your own to decipher it. Hopefully, it works.

[md.muniruzzaman]

Yes, you can use KINETICS to simulate kinetic sorption, although it is not very easy.

Attached is a file that kinetically sorbs a constituent "A" using a Langmuir isotherm. The idea is to calculate what surface sorption would be given the dissolved concentration and use the difference between that and what is actually sorbed to formulate a rate equation.

Also attached is an example of kinetic sorption for a surface complexation example. It is not my example, and it has been a long time since I looked at it, so you are on your own to decipher it. Hopefully, it works.

Hi David, I think the attachments are still not functional! Would you happen to have these phreeqc scripts for kinetic sorption / surface complexation?

[dlparkhurst]

I think this file may be similar to the KinSurf.pqi script.

Code: [Select]

DATABASE ../database/wateq4f.dat
SELECTED_OUTPUT 101
    -file                 kinsurf_101.sel
USER_PUNCH 101
    -headings Mu SC
    -start
10 PUNCH STR_F$(MU, 20, 12)
20 PUNCH STR_F$(SC, 20, 10)
    -end
SURFACE_SPECIES
#   All surface data from
#   Dzombak and Morel, 1990
#
#   weak binding site--Hfo_w

        Hfo_wOH = Hfo_wOH
        log_k  0.0

        Hfo_wOH  + H+ = Hfo_wOH2+
        log_k  7.29    # = pKa1,int

        Hfo_wOH = Hfo_wO- + H+
        log_k  -8.93   # = -pKa2,int

###############################################
#             ANIONS                          #
###############################################
#
#   Anions from table 10.6
#
#   Phosphate
        Hfo_wOH + PO4-3 + 3H+ = Hfo_wH2PO4 + H2O
        log_k   31.29

        Hfo_wOH + PO4-3 + 2H+ = Hfo_wHPO4- + H2O
#        log_k   25.39
        log_k   -25.39

        Hfo_wOH + PO4-3 + H+ = Hfo_wPO4-2 + H2O
#        log_k   17.72
        log_k   -17.72
#   Arsenate
        Hfo_wOH + AsO4-3 + 3H+ = Hfo_wH2AsO4 + H2O
        log_k   29.31

        Hfo_wOH + AsO4-3 + 2H+ = Hfo_wHAsO4- + H2O
#        log_k   23.51
        log_k   -23.51

        Hfo_wOH + AsO4-3 = Hfo_wOHAsO4-3
#        log_k   10.58
        log_k   -10.58

SOLUTION
pH  4
pe 10
As(5) 1
P  1
Na 10
Cl 10   charge
END
PHASES
fix_pH
H+ = H+
log_k 0
USE solution 1
EQUILIBRIUM_PHASES
# fix_ph  -7  HCl
# O2(g)   -0.7
SURFACE
-no_edl
Hfo_wOH   .001
#*SURFACE_SPECIES
#*#   All surface data from
#*#   Dzombak and Morel, 1990
#*#
#*#   weak binding site--Hfo_w
#*
#*        Hfo_wOH = Hfo_wOH
#*        log_k  0.0
#*
#*        Hfo_wOH  + H+ = Hfo_wOH2+
#*        log_k  7.29    # = pKa1,int
#*
#*        Hfo_wOH = Hfo_wO- + H+
#*        log_k  -8.93   # = -pKa2,int
#*
#*###############################################
#*#             ANIONS                          #
#*###############################################
#*#
#*#   Anions from table 10.6
#*#
#*#   Phosphate
#*        Hfo_wOH + PO4-3 + 3H+ = Hfo_wH2PO4 + H2O
#*        log_k   31.29
#*
#*        Hfo_wOH + PO4-3 + 2H+ = Hfo_wHPO4- + H2O
#*        log_k   25.39
#*
#*        Hfo_wOH + PO4-3 + H+ = Hfo_wPO4-2 + H2O
#*        log_k   17.72
#*#   Arsenate
#*        Hfo_wOH + AsO4-3 + 3H+ = Hfo_wH2AsO4 + H2O
#*        log_k   29.31
#*
#*        Hfo_wOH + AsO4-3 + 2H+ = Hfo_wHAsO4- + H2O
#*        log_k   23.51
#*
#*        Hfo_wOH + AsO4-3 = Hfo_wOHAsO4-3
#*        log_k   10.58
RATES


HFO_OH2+       # species 2
-start
10  k = 1/(24*3600)
20 s$ = "HFO_OH2+"
30 sites = 1
30 sites = KIN("HFO_OH") + KIN("HFO_OH2+") + KIN("HFO_O-") + KIN("HFO_H2PO4") + KIN("HFO_H2AsO4")
40 species = KIN(s$)   # moles
50 log_master_frac = LOG10(KIN("HFO_OH")/sites)
60 if (species > 0) then log_species_frac = LOG10(species/sites) else log_species_frac = -100
70 q = log_species_frac - LA("H+") - log_master_frac
80 satrat = 10^(q - 7.29)
85 put (q - 7.29, 2, 1)
90 rate = k*(1-satrat)
95 if (species <= 0) then rate = k
100 moles = rate*TIME
110 SAVE -moles
120 PUT(moles, 2)
-end

HFO_O-         # species 3
-start
10  k = 1/(24*3600)
20 s$ = "HFO_O-"
30 sites = 1
30 sites = KIN("HFO_OH") + KIN("HFO_OH2+") + KIN("HFO_O-") + KIN("HFO_H2PO4") + KIN("HFO_H2AsO4")
40 species = KIN(s$)   # moles
50 log_master_frac = LOG10(KIN("HFO_OH")/sites)
60 if (species > 0) then log_species_frac = LOG10(species/sites) else log_species_frac = -100
70 q = log_species_frac + LA("H+") - log_master_frac
80 satrat = 10^(q - -8.93)
85 put (q + 8.93, 3, 1)
90 rate = k*(1-satrat)
95 if (species <= 0) then rate = k
100 moles = rate*TIME
110 SAVE -moles
120 PUT(moles, 3)
-end

HFO_H2PO4      # species 4
-start
10  k = 1/(24*3600)
20 s$ = "HFO_H2PO4"
30 sites = 1
30 sites = KIN("HFO_OH") + KIN("HFO_OH2+") + KIN("HFO_O-") + KIN("HFO_H2PO4") + KIN("HFO_H2AsO4")
40 species = KIN(s$)   # moles
50 log_master_frac = LOG10(KIN("HFO_OH")/sites)
60 if (species > 0) then log_species_frac = LOG10(species/sites) else log_species_frac = -100
70 q = log_species_frac + LA("H2O") -3*LA("H+") - LA("PO4-3") - log_master_frac
80 satrat = 10^(q - 31.29)
85 put (q - 31.29, 4, 1)
90 rate = k*(1-satrat)
95 if (species <= 0) then rate = k
100 moles = rate*TIME
110 SAVE -moles
120 PUT(moles, 4)
-end

HFO_H2AsO4     # species 5
-start
10  k = 1/(24*3600)
20 s$ = "HFO_H2AsO4"
30 sites = 1
30 sites = KIN("HFO_OH") + KIN("HFO_OH2+") + KIN("HFO_O-") + KIN("HFO_H2PO4") + KIN("HFO_H2AsO4")
40 species = KIN(s$)   # moles
50 log_master_frac = LOG10(KIN("HFO_OH")/sites)
60 if (species > 0) then log_species_frac = LOG10(species/sites) else log_species_frac = -100
70 q = log_species_frac + LA("H2O") -3*LA("H+") - LA("AsO4-3") - log_master_frac
80 satrat = 10^(q - 29.31)
85 put (q - 29.31, 5, 1)
90 rate = k*(1-satrat)
95 if (species <= 0) then rate = k
100 moles = rate*TIME
110 SAVE -moles
120 PUT(moles, 5)
-end

HFO_OH   # species 1
-start
#
# HFO_OH is difference after all other surface kinetics has been calculates
#
10 sites = 1
10 sites = KIN("HFO_OH") + KIN("HFO_OH2+") + KIN("HFO_O-") + KIN("HFO_H2PO4") + KIN("HFO_H2AsO4")
20 frac = KIN("HFO_OH")/sites
30 moles = (get(2) + get(3) + get(4) + get(5))
#30 moles  = 0
40 save moles
-end

SOLUTION
pH  4
-pe 10
As(5) 1
P  1
Na 10
Cl 10 charge
END
PHASES
fix_pH
H+ = H+
log_k 0
USE solution 1
EQUILIBRIUM_PHASES
# fix_ph  -7  HCl
# O2(g)   -0.7
KINETICS
-cvode
-step 10000 in 1 steps
     HFO_OH2+
  -formula   OH2  1 
          -M 1e-8

     HFO_O-
#   -formula   O 1
  -formula   Na2O 1
          -M 1e-8

     HFO_H2PO4
#          -formula   H2PO4  1
          -formula   NaH2PO4  1
          -M 1e-8

     HFO_H2AsO4
#   -formula   H2AsO4  1
  -formula   NaH2AsO4  1
          -M 1e-8

     HFO_OH     
#   -formula   OH 1
  -formula   NaOH 1
          -M .001

INCREMENTAL_REACTIONS
PRINT
-warnings 1
USER_PRINT
10 sites = KIN("HFO_OH") + KIN("HFO_OH2+") + KIN("HFO_O-") + KIN("HFO_H2PO4") + KIN("HFO_H2AsO4")
20 print "Sum of species = sites:    ", sites
30 charge = KIN("HFO_OH2+") - KIN("HFO_O-")
40 print "Charge on surface (eq):    ", charge
50 print "log Q/K HFO_OH2+:          ", get(2,1)
60 print "log Q/K HFO_O-:            ", get(3,1)
70 print "log Q/K HFO_H2PO4:         ", get(4,1)
80 print "log Q/K HFO_H2AsO4:        ", get(5,1)
END


[md.muniruzzaman]

Thank you very much!

[khwee]

Hi David - Could you also post the code from KinLang.pqi? The attachment did not work for me as well. Thanks!

[dlparkhurst]

Isotherms are not my favorite, and I do not have the original file. Here is a script that I begrudgingly recreated for kinetic Langmuir sorption that may or may not be similar to the original post.

The RATES definition calculates the equilibrium distribution between solution and sorbed "A", and then kinetically approaches the equilibrium distribution.

Code: [Select]

RATES
A_langmuir
#Q0 = parm(1)  mg/g sediment
#KL = parm(2)  L/mg
#g  = parm(3)  g of sediment
#Ce            mg/L  adsorbate
#qe            mg/g  sorbed
# qe = Q0 * KL * Ce / (1 + KL*Ce)
10 PUT(parm(1), 1)
20 Q0 = parm(1)
30 PUT(parm(2), 2)
40 KL = parm(2)
50 PUT(parm(3), 3)
60 g = parm(3)
70 Ce = TOT("A") * GFW("A") * 1000 # mg/L
80 A_sys = TOT("A") + M
85 PRINT "x", A_sys, TOT("A"), M
90 FOR i = 1 TO 20
100 qe =  Q0 * KL * Ce / (1 + KL*Ce)
110 A_calc = (qe * g + Ce) / (GFW("A")*1000)
120 Ce_new = Ce * A_sys / A_calc
130 REM PRINT Ce_new, Ce
140 IF [ABS(Ce_new - Ce) / Ce < 1e-8] THEN GOTO 200
150 Ce = Ce_new
160 NEXT i
170 PRINT "A_langmuir rate calculation failed.----------------------"
200 k = 1/3600
210 rate = k * [M - qe * g / (GFW("A")*1000)]
220 moles = rate * TIME
230 SAVE moles


-end
SOLUTION_MASTER_SPECIES
   A  A 0 10 10
SOLUTION_SPECIES
   A = A
   log_k 0
SOLUTION 1
      pH  7
END
REACTION 1
      A  1
      1e-6 #in 20
END
USE solution 1
USE reaction 1
KINETICS
A_langmuir
-cvode
-m 0
-formula A 1
-parm 1 1 1
-time 36000 in 10
USER_PRINT
10 Q0 = GET(1)
20 KL = GET(2)
30 g = GET(3)
40 Ce = TOT("A") * GFW("A") * 1000 # mg / L
50 qe =  Q0 * KL * Ce / (1 + KL*Ce)
60 qe_check = KIN("A_langmuir") * GFW("A") * 1000 / g
70 PRINT qe, qe_check
USER_GRAPH 1
    -headings               hr A_equilibrium A_kinetic A_sorb_equilibrium A_sorb_kinetic
    -axis_titles            "Time, hours" "A Dissolved, mg/L" "A Sorbed, mg"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 20 Q0 = GET(1)
 40 KL = GET(2)
 60 g = GET(3)
 70 Ce = TOT("A") * GFW("A") * 1000
 80 A_sys = TOT("A") + KIN("A_langmuir")
 85 PRINT A_sys, TOT("A"), KIN("A_langmuir")
 90 FOR i = 1 TO 20
100 qe =  Q0 * KL * Ce / (1 + KL*Ce)
110 A_calc = (qe * g + Ce) / (GFW("A")*1000)
120 Ce_new = Ce * A_sys / A_calc
140 IF [ABS(Ce_new - Ce) / Ce < 1e-8] THEN GOTO 200
150 Ce = Ce_new
160 NEXT i
200 A_diss_eq_mg_L = Ce
220 A_sorb_eq_mg = qe * g
300 GRAPH_X TOTAL_TIME / 3600
310 GRAPH_Y A_diss_eq_mg_L, TOT("A") * GFW("A") * 1000
320 GRAPH_SY A_sorb_eq_mg, KIN("A_langmuir") * GFW("A") * 1000
  -end
    -active                 true
END
 

[khwee]

Thanks so much David! This sample script will help a lot :)

[dlparkhurst]

I'm not sure how, but some non-printing characters were added at the beginning of lines 20 through 90 of USER_GRAPH. I also forgot to comment some debugging prints in the Basic.

Code: [Select]

RATES
A_langmuir
#Q0 = parm(1)  mg/g sediment
#KL = parm(2)  L/mg
#g  = parm(3)  g of sediment
#Ce            mg/L  adsorbate
#qe            mg/g  sorbed
# qe = Q0 * KL * Ce / (1 + KL*Ce)
10 PUT(parm(1), 1)
20 Q0 = parm(1)
30 PUT(parm(2), 2)
40 KL = parm(2)
50 PUT(parm(3), 3)
60 g = parm(3)
70 Ce = TOT("A") * GFW("A") * 1000 # mg/L
80 A_sys = TOT("A") + M
#85 PRINT "x", A_sys, TOT("A"), M
90 FOR i = 1 TO 20
100 qe =  Q0 * KL * Ce / (1 + KL*Ce)
110 A_calc = (qe * g + Ce) / (GFW("A")*1000)
120 Ce_new = Ce * A_sys / A_calc
130 REM PRINT Ce_new, Ce
140 IF [ABS(Ce_new - Ce) / Ce < 1e-8] THEN GOTO 200
150 Ce = Ce_new
160 NEXT i
170 PRINT "A_langmuir rate calculation failed.----------------------"
200 k = 1/3600
210 rate = k * [M - qe * g / (GFW("A")*1000)]
220 moles = rate * TIME
230 SAVE moles


-end
SOLUTION_MASTER_SPECIES
   A  A 0 10 10
SOLUTION_SPECIES
   A = A
   log_k 0
SOLUTION 1
      pH  7
END
REACTION 1
      A  1
      1e-6 #in 20
END
USE solution 1
USE reaction 1
KINETICS
A_langmuir
-cvode
-m 0
-formula A 1
-parm 1 1 1
-time 36000 in 10
USER_PRINT
10 Q0 = GET(1)
20 KL = GET(2)
30 g = GET(3)
40 Ce = TOT("A") * GFW("A") * 1000 # mg / L
50 qe =  Q0 * KL * Ce / (1 + KL*Ce)
60 qe_check = KIN("A_langmuir") * GFW("A") * 1000 / g
70 PRINT qe, qe_check
USER_GRAPH 1
    -headings               hr A_equilibrium A_kinetic A_sorb_equilibrium A_sorb_kinetic
    -axis_titles            "Time, hours" "A Dissolved, mg/L" "A Sorbed, mg"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
20 Q0 = GET(1)
40 KL = GET(2)
60 g = GET(3)
70 Ce = TOT("A") * GFW("A") * 1000
80 A_sys = TOT("A") + KIN("A_langmuir")
#85 PRINT A_sys, TOT("A"), KIN("A_langmuir")
90 FOR i = 1 TO 20
100 qe =  Q0 * KL * Ce / (1 + KL*Ce)
110 A_calc = (qe * g + Ce) / (GFW("A")*1000)
120 Ce_new = Ce * A_sys / A_calc
140 IF [ABS(Ce_new - Ce) / Ce < 1e-8] THEN GOTO 200
150 Ce = Ce_new
160 NEXT i
200 A_diss_eq_mg_L = Ce
220 A_sorb_eq_mg = qe * g
300 GRAPH_X TOTAL_TIME / 3600
310 GRAPH_Y A_diss_eq_mg_L, TOT("A") * GFW("A") * 1000
320 GRAPH_SY A_sorb_eq_mg, KIN("A_langmuir") * GFW("A") * 1000
  -end
    -active                 true
END