Converging error about cation exchange and surface complexation coupling model

[Jeonghwan Hwang]

Hello, this is Jeonghwan Hwang.
I studied to simulate uranium sorption and transport using cation exchange and surface complexation bonding models.

While running the sorption model, I found an error for an unknown reason.
When I looked at the error, I confirmed that the capacity of the cation exchange changed very seriously.
Then, this caused the severe increasing for Na concentration than Cl concentration, and the model did not converge.

But I can't quite find the solution of this error.
I hope to take an advise for solving this problem.

I attached my example code below;
=============================================================================

DATABASE c:\phreeqc\database\Sit.dat


SURFACE_MASTER_SPECIES # b_Bentonite, b_illite = Missana (2021)
        Bentonite_s     Bentonite_sOH #s
        Bentonitea_w     Bentonitea_wOH #w1
   Bentoniteb_w   Bentoniteb_wOH #w2
   
   Illite_s     Illite_sOH #s
        Illitea_w     Illitea_wOH #w1
   Illiteb_w   Illiteb_wOH #w2

SURFACE_SPECIES
#Montmorillonite SWY-1
        Bentonite_sOH = Bentonite_sOH
        log_k  0.0
   
        Bentonitea_wOH = Bentonitea_wOH
        log_k  0.0   

   Bentoniteb_wOH = Bentoniteb_wOH
        log_k  0.0   
#bentonite properties with tirtation curve   
        Bentonite_sOH  + H+ = Bentonite_sOH2+
        log_k  4.5
 
        Bentonite_sOH = Bentonite_sO- + H+
        log_k  -7.9
   
   Bentonitea_wOH  + H+ = Bentonitea_wOH2+
        log_k  4.5
 
        Bentonitea_wOH = Bentonitea_wO- + H+
        log_k  -7.9
   
   Bentoniteb_wOH  + H+ = Bentoniteb_wOH2+
        log_k  6.0
 
        Bentoniteb_wOH = Bentoniteb_wO- + H+
        log_k  -10.5   
#Strong site      
        Bentonite_sOH + UO2+2 = Bentonite_sOUO2+ + H+
        log_k  3.1
   
   Bentonite_sOH + UO2+2 + H2O = Bentonite_sOUO2OH + 2H+
        log_k  -4.6

   Bentonite_sOH + UO2+2 + 2H2O = Bentonite_sOUO2(OH)2- + 3H+
        log_k  -12.6

   Bentonite_sOH + UO2+2 + 3H2O = Bentonite_sOUO2(OH)3-2 + 4H+
        log_k  -20.9
#Weak site
        Bentonitea_wOH  + UO2+2 = Bentonitea_wOUO2+ + H+
        log_k  0.5
 
        Bentonitea_wOH  + UO2+2 + H2O = Bentonitea_wOUO2OH + 2H+
        log_k  -5.7


#Illite   
        Illite_sOH = Illite_sOH
        log_k  0.0
   
        Illitea_wOH = Illitea_wOH
        log_k  0.0   

   Illiteb_wOH = Illiteb_wOH
        log_k  0.0   
   

#Illite properties with tirtation curve   
        Illite_sOH  + H+ = Illite_sOH2+
        log_k  4.0 #Bradbury 2009a, Geochi et Cosmo
 
        Illite_sOH = Illite_sO- + H+
        log_k  -6.2 #Bradbury 2009a, Geochi et Cosmo
   
   Illitea_wOH  + H+ = Illitea_wOH2+
        log_k  4.0 #Bradbury 2009a, Geochi et Cosmo
 
        Illitea_wOH = Illitea_wO- + H+
        log_k  -6.2 #Bradbury 2009a, Geochi et Cosmo
   
   Illiteb_wOH  + H+ = Illiteb_wOH2+
        log_k  8.5 #Bradbury 2009a, Geochi et Cosmo
 
        Illiteb_wOH = Illiteb_wO- + H+
        log_k  -10.5 #Bradbury 2009a, Geochi et Cosmo
   
   
#Strong site      
        Illite_sOH + UO2+2 = Illite_sOUO2+ + H+
        log_k  2.0     #Bradbury and Baeyens, 2009b
   
   Illite_sOH + UO2+2 + H2O = Illite_sOUO2OH + 2H+
        log_k  -3.5    #Bradbury and Baeyens, 2009b 3.5

   Illite_sOH + UO2+2 + 2H2O = Illite_sOUO2(OH)2- + 3H+
        log_k  -10.6   #Bradbury and Baeyens, 2009b 10.6

   Illite_sOH + UO2+2 + 3H2O = Illite_sOUO2(OH)3-2 + 4H+
        log_k  -19.0   #Bradbury and Baeyens, 2009b  19

#Weak site
        Illitea_wOH  + UO2+2 = Illitea_wOUO2+ + H+
        log_k  0.1   #Bradbury 2005, Geochi et Cosmo
 
        Illitea_wOH  + UO2+2 + H2O = Illitea_wOUO2OH + 2H+
        log_k  -5.3  #Bradbury 2005, Geochi et Cosmo
   
   
   
EXCHANGE_MASTER_SPECIES
   Bentonite_ex    Bentonite_ex-   
   Illite_ex    Illite_ex-
      
EXCHANGE_SPECIES

#Montmorillonite   
   Bentonite_ex- = Bentonite_ex-
   log_K 0
   
   Na+ + Bentonite_ex- = NaBentonite_ex
   log_K 0
   
   2 NaBentonite_ex + UO2+2 = UO2(Bentonite_ex)2 + 2Na+
   log_K 0.45
   
#Illite   
   Illite_ex- = Illite_ex-
   log_K 0 #Bradbury 2009a, Geochi et Cosmo
   
   Na+ + Illite_ex- = NaIllite_ex
   log_K 0 #Bradbury 2009a, Geochi et Cosmo

       2 NaIllite_ex + UO2+2 = UO2(Illite_ex)2 + 2Na+
   log_K 0.65   #Bradbury 2009b, Geochi et Cosmo   


       Surface 1
          sites_units absolute
          Bentonite_sOH 1.0873e-07
          Bentonitea_wOH 0.0013726
          Bentoniteb_wOH 0.0013726
          Illite_sOH 2.6246e-05
          Illitea_wOH 0.15348
          Illiteb_wOH 0.15348
          no_edl
          end
          Exchange 1
          NaBentonite_ex 0.034842
          UO2(Bentonite_ex)2 1.794e-05
          NaIllite_ex 0.89208
          UO2(Illite_ex)2 0.00072753
          -exchange_gammas false
          end
          SOLUTION 1
          -units  mmol/kgw
         -water 0.4
          temp 25
           pH  7
          U(+6) 0.0063776
          Na 0.015109
          Cl 0.017122. charge
         EQUILIBRIUM_PHASES 1
        save solution 1
          end

     USE solution 1
         USE surface 1
         USE Exchange 1
         EQUILIBRIUM_PHASES 1
         end


=============================================================================

Thank you.

Sincerely,

Jeonghwan Hwang

[dlparkhurst]

There seems to be a bug, I think in the SAVE method, but I'll have to look at the code.

If you remove EQUILIBRIUM_PHASES 1, as in the script below, your script will run. Note that EQUILIBRIUM_PHASES 1 has no minerals, so no reaction actually occurs and the SOLUTION is the same anyway. Oddly, if a mineral is included in EQUILIBRIUM_PHASES, the script runs. Also, if you set -exchange_gammas to true, the script also runs.

My suggestion is to remove the do-nothing reaction with EQUILIBRIUM_PHASES (as below). At least in this case, the script will run. If you have another convergence problem for a non-trivial reaction, then I would like to see the input file.

Code: [Select]

SURFACE_MASTER_SPECIES # b_Bentonite, b_illite = Missana (2021)
   Bentonite_s     Bentonite_sOH #s
   Bentonitea_w    Bentonitea_wOH #w1
   Bentoniteb_w    Bentoniteb_wOH #w2
   Illite_s        Illite_sOH #s
   Illitea_w       Illitea_wOH #w1
   Illiteb_w       Illiteb_wOH #w2

SURFACE_SPECIES
#Montmorillonite SWY-1
   Bentonite_sOH = Bentonite_sOH
        log_k  0.0
   
   Bentonitea_wOH = Bentonitea_wOH
        log_k  0.0   

   Bentoniteb_wOH = Bentoniteb_wOH
        log_k  0.0   
#bentonite properties with tirtation curve   
   Bentonite_sOH  + H+ = Bentonite_sOH2+
        log_k  4.5
 
   Bentonite_sOH = Bentonite_sO- + H+
        log_k  -7.9
   
   Bentonitea_wOH  + H+ = Bentonitea_wOH2+
        log_k  4.5
 
   Bentonitea_wOH = Bentonitea_wO- + H+
        log_k  -7.9
   
   Bentoniteb_wOH  + H+ = Bentoniteb_wOH2+
        log_k  6.0
 
   Bentoniteb_wOH = Bentoniteb_wO- + H+
        log_k  -10.5   
#Strong site     
   Bentonite_sOH + UO2+2 = Bentonite_sOUO2+ + H+
        log_k  3.1
   
   Bentonite_sOH + UO2+2 + H2O = Bentonite_sOUO2OH + 2H+
        log_k  -4.6

   Bentonite_sOH + UO2+2 + 2H2O = Bentonite_sOUO2(OH)2- + 3H+
        log_k  -12.6

   Bentonite_sOH + UO2+2 + 3H2O = Bentonite_sOUO2(OH)3-2 + 4H+
        log_k  -20.9
#Weak site
   Bentonitea_wOH  + UO2+2 = Bentonitea_wOUO2+ + H+
        log_k  0.5
 
   Bentonitea_wOH  + UO2+2 + H2O = Bentonitea_wOUO2OH + 2H+
        log_k  -5.7


#Illite   
   Illite_sOH = Illite_sOH
        log_k  0.0
   
   Illitea_wOH = Illitea_wOH
        log_k  0.0   

   Illiteb_wOH = Illiteb_wOH
        log_k  0.0   

#Illite properties with tirtation curve   
   Illite_sOH  + H+ = Illite_sOH2+
        log_k  4.0 #Bradbury 2009a, Geochi et Cosmo
 
   Illite_sOH = Illite_sO- + H+
        log_k  -6.2 #Bradbury 2009a, Geochi et Cosmo
   
   Illitea_wOH  + H+ = Illitea_wOH2+
        log_k  4.0 #Bradbury 2009a, Geochi et Cosmo
 
   Illitea_wOH = Illitea_wO- + H+
        log_k  -6.2 #Bradbury 2009a, Geochi et Cosmo
   
   Illiteb_wOH  + H+ = Illiteb_wOH2+
        log_k  8.5 #Bradbury 2009a, Geochi et Cosmo
 
   Illiteb_wOH = Illiteb_wO- + H+
        log_k  -10.5 #Bradbury 2009a, Geochi et Cosmo
   
#Strong site     
   Illite_sOH + UO2+2 = Illite_sOUO2+ + H+
        log_k  2.0     #Bradbury and Baeyens, 2009b
   
   Illite_sOH + UO2+2 + H2O = Illite_sOUO2OH + 2H+
        log_k  -3.5    #Bradbury and Baeyens, 2009b 3.5

   Illite_sOH + UO2+2 + 2H2O = Illite_sOUO2(OH)2- + 3H+
        log_k  -10.6   #Bradbury and Baeyens, 2009b 10.6

   Illite_sOH + UO2+2 + 3H2O = Illite_sOUO2(OH)3-2 + 4H+
        log_k  -19.0   #Bradbury and Baeyens, 2009b  19

#Weak site
   Illitea_wOH  + UO2+2 = Illitea_wOUO2+ + H+
        log_k  0.1   #Bradbury 2005, Geochi et Cosmo
 
   Illitea_wOH  + UO2+2 + H2O = Illitea_wOUO2OH + 2H+
        log_k  -5.3  #Bradbury 2005, Geochi et Cosmo
   
EXCHANGE_MASTER_SPECIES
   Bentonite_ex    Bentonite_ex-   
   Illite_ex       Illite_ex-
     
EXCHANGE_SPECIES
#Montmorillonite   
   Bentonite_ex- = Bentonite_ex-
   log_K 0
   
   Na+ + Bentonite_ex- = NaBentonite_ex
   log_K 0
   
   2 NaBentonite_ex + UO2+2 = UO2(Bentonite_ex)2 + 2Na+
   log_K 0.45
   
#Illite   
   Illite_ex- = Illite_ex-
   log_K 0 #Bradbury 2009a, Geochi et Cosmo
   
   Na+ + Illite_ex- = NaIllite_ex
   log_K 0 #Bradbury 2009a, Geochi et Cosmo

       2 NaIllite_ex + UO2+2 = UO2(Illite_ex)2 + 2Na+
   log_K 0.65   #Bradbury 2009b, Geochi et Cosmo   
END
SURFACE 1
    Bentonite_sOH    1.0873e-07
    Bentonitea_wOH    0.0013726
    Bentoniteb_wOH    0.0013726
    Illite_sOH    2.6246e-05
    Illitea_wOH    0.15348
    Illiteb_wOH    0.15348
    -no_edl
END
EXCHANGE 1
NaBentonite_ex 0.034842
UO2(Bentonite_ex)2 1.794e-05
NaIllite_ex 0.89208
UO2(Illite_ex)2 0.00072753
-exchange_gammas false
END
SOLUTION 1
    temp      25
    pH        7
    Cl        0.017122 charge
    Na        0.015109
    U(6)      0.0063776
    -water    0.4 # kg
#EQUILIBRIUM_PHASES 1
#SAVE solution 1
END

USE exchange 1
USE solution 1
USE surface 1
EQUILIBRIUM_PHASES 1
END

[pawanagra0@gmail.com]

Hello,

I am doing a similar job. I am also trying to fit the experimental distribution coefficient by using the available TSMs. This difference is that I am doing the batch sorption experiment with biotite at 25 with 0.001M ionic strength where pH is varying from 5-9.

Since I have run this experiment for 3 months. The ICP-MS analysis suggests that there is a dissolution of biotite.

Can you please suggest by giving a rough example that how can I include the dissolution of biotite in the surface complexation model?

Have a good weekend

Thank you

[dlparkhurst]

I don't have enough information to work with.

If you are modeling the surface of biotite, the change in surface sorption sites may not be significant and you could ignore it. If not, then you need to link the surface to either a KINETIC reaction of biotite or an EQUILIBRIUM_PHASE of biotite. See the manual for the SURFACE data block.

If you just want to dissolve some biotite, you can do it with REACTION, KINETICS, or EQUILIBRIUM_PHASES.

[pawanagra0@gmail.com]

Hello Dr.

Thank you so much for your reply.

I could not find out the dissolution reaction for the biotite. Therefore I am using annite mineral for modeling the biotite dissolution result.

Below is the data that I got from ICPMS measurement of mineral dissolution.

exp Data ICPMS      
pH     Si(mol(0.005 l)    Fe(mol(0.005 l)
5,67      1,3E-04           7,4E-06
6,08      1,0E-04           3,3E-06
7,10      1,0E-04           1,5E-06
8,05      1,5E-04           1,7E-06
8,58      1,9E-04           2,2E-06

Hers is my code which is not giving me results anywhere close to what I got from experiment. Can you please look at it and tell me how can I get the correct results?

SELECTED_OUTPUT 1
    -file                 C:\My_Folder\fre\selected_output_1.txt
    -reset                false
    -pH                   true
    -molalities           Mg+2  SiO2
USER_PUNCH
    -start
10 PUNCH MOL("Mg+2")
20 PUNCH MOL("SiO2")
    -end
END
SOLUTION_MASTER_SPECIES
    Cl            Cl-              0     Cl              35.4527
    Cl(-1)        Cl-              0     Cl             
    Cl(1)         ClO-             0     Cl             
    Cl(3)         ClO2-            0     Cl             
    Cl(5)         ClO3-            0     Cl             
    Cl(7)         ClO4-            0     Cl
#    Ti         TiO2       0.0    Ti      47.88         
SOLUTION_SPECIES
Cl- + 0.5O2 = ClO-
    log_k     -15.1014
    delta_h   -66.0361 kJ
Cl- + O2 = ClO2-
    log_k     -23.108
    delta_h   -112.688 kJ
Cl- + 1.5O2 = ClO3-
    log_k     -17.2608
    delta_h   -81.3077 kJ
Cl- + 2O2 = ClO4-
    log_k     100
    delta_h   -62.0194 kJ
ClO- + H+ = HClO
    log_k     7.5692
ClO2- + H+ = HClO2
    log_k     3.1698
#TiO2 =  TiO2
#log_k     0   

END

SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1
    Cl(7)     0.001
    Na        0.001
    -water    0.005 # kg

PHASES
Fix_H+
    H+ = H+
    log_k     0

USE solution 1
EQUILIBRIUM_PHASES 1
Fix_H+ -5.67  HCl  1.0
annite
end
USE solution 1
EQUILIBRIUM_PHASES 2
Fix_H+   -6.08  HCl  1.0
annite 
end
USE solution 1
EQUILIBRIUM_PHASES 3
Fix_H+   -7.25  NaOH  1.0
annite   
end

USE solution 1
EQUILIBRIUM_PHASES 3
Fix_H+   -8.05  NaOH  1.0
annite   
end
USE solution 1
EQUILIBRIUM_PHASES 4
Fix_H+   -8.58  NaOH  1.0
annite   
end

I will be very thankful to you!

Have a good day

Pawan
 

[dlparkhurst]

I'm not sure what the units are for you analytical data, but there is less Fe than Si. If you simply dissolve Annite, with no secondary reactions, then you will get equimolar Fe and Si.

One choice you have to make is the with perchlorate. If you want perchlorate to oxidize the iron, your definitions are correct. If you want perchlorate to be an inert electrolyte, then for simplicity, use Cl(-1) instead.

With the calculation as you have defined it, most of the iron is oxidized to Fe(3) and is supersaturated with iron oxyhydroxide phases. If you include, say, goethite in your equilibrium phases, then most of the iron will precipitate, which is more consistent with your experiment. Is your solution turning red from formation of iron oxyhydroxides?

Your reaction solutions are also supersaturated with Kaolinite and probably other clays, so secondary aluminosilicate precipitation is also possible.

Finally, because you don't know the phase or its log K, you can adjust the log K for Annite as part of the fitting process.

[pawanagra0@gmail.com]

Pardon me, but the units are in "moles". Yes, there is less iron than Si.

Actually, I did my batch sorption experiment with an inert 0.001M sodium perchlorate solution inside the glove box at the given pH for two months.

when I am defining C(-1). I still can see the Fe+3 formation.
 
No, there was no change in the color of the solution at the end of my experiments. which means that even if there is an oxidation of Fe then it has to be very low.

I choose the mineral "annite" because biotite's chemical compositional range lies between annite and phlogopite (Deer et al. (1992)).

Yes, I can optimize the constant for annite. Do you think that constant can be a problem coz it is given in the llnl.dat

Do you think there has to be something else that needs to be added to the program to get theoretical data similar to the experiential?

[dlparkhurst]

Where is the iron? There are lots of possibilities. Do you know the stoichiometry of Fe:Si your mineral? Even if you do, there may be selective dissolution such that iron is retained in the solid or iron precipitates as a secondary mineral. The iron could be sorbed on the surface. There is no guarantee that Si and Al are not forming secondary minerals too.

You don't know whether the experiment reaches equilibrium with the mineral or some kind of steady state, so combined with the uncertainty in stoichiometry, it is difficult to estimate a log K. 

I am glad you are doing the experiments, but I think you have a lot of work to determine what is happening in your system.