Speciation Error using Phreeplot

[FlouretAlex]

Hello, everyone.
I have an issue using the Phreeplot software.
By using the code below at the end of the run I have the following error: "Speciation error for data point 1 (fit data file line 2)."


I can't see what the problem with this one if someone can have an idea this will help me a lot.

Thanks,
Alex

PHREEPLOT SCRIPT.

# Illite fitting to soil-e isotherm using unit weighting
#   basic fit with isotherm plot

SPECIATION
  jobTitle                             "Test fitting: absolute errors (unit weighting)"
  calculationType                      fit
  calculationMethod                    1
  database                BRGM_database_phreeqc_ThermoddemV1.10_06Jun2017.dat
 
FIT
  dataFile                             soil_e_imput.dat                       # contains a list of observations and independent variables
  onepass                              FALSE                         # does a separate PHREEQC simulation for each observation -
                                                                     #   slow but easy to set up
  dependentVariableColumnObs           Cssorbed                      # name of column in fit data file containing the observations (dep variable)
  dependentVariableColumnCalc          sorbCs                        # name of column in selected output file containing the calcd values of the dep variable
  fitWeightingMethod                   0                             # 0 = unit weights
  fitFiniteDiffStepSize                1.0E-7                        # initial step size for each adjustable parameter
  mainLoopColumn                       sim                           # column header in fit data file for which PHREEQC simulation(s) to use for each observation
  numberOfFitParameters                1
  fitParameterNames                    M1
  fitLogParameters                     0                          # 0 = parameters on a linear scale, 1 = log10 parameters before fitting
  fitAdjustableParameters              1                          # 1 = adjustable
  fitParameterValues                   1.0E-06                   # initial values (starting point)

PLOT
  plotTitle                            "Cs sorption on Soil-E (illite only)"
  xtitle                               "Cs_ini (mol/L)"
  ytitle                               "Cs sorbed (mol)"
                                                                     # plot isotherm (x = Csconcn, y = sorbed)
  lines                                calculated                    # y = line from calculated values from out file with column heading = 'calculated'
  points                               observed                      # y = points from observed values from out file with column heading = 'observed'

  lineWidth                            0.4
  lineColor                            red
  labelSize                            0.0                           # no labels on plot
  legendTextSize                       0.0                           # no legend (key)
  pointSize                            4.0
  customXcolumn                        Csconcn                       # x = Csconcn column in out file

CHEMISTRY

PHASES
Fix_H+
  H+ = H+
  log_k 0.0

SURFACE_MASTER_SPECIES
   Illite_  Illite_OH+0.5
   
SURFACE_SPECIES
##### Illite_SOH

Illite_OH+0.5 = Illite_OH+0.5
Log_k 0.0

Illite_OH+0.5 =     Illite_O-0.5  +  H+
Log_k -3.46


Illite_OH+0.5 + Cs+ = Illite_OCs+0.5 + H+
Log_k 5.2

Illite_OH+0.5 + Na+   = Illite_ONa+0.5  +  H+
 Log_k -1.8

Illite_OH+0.5 + K+ = Illite_OK+0.5    +  H+
Log_k 0.6

Illite_OH+0.5 + Ca+2  = Illite_OCa+1.5     + H+
Log_k  -5

#   Illite_OH+0.5 + NH4+ = Illite_ONH4+0.5    +  H+
#   Log_k 1.5

SELECTED_OUTPUT
 -high_precision true
 -reset false
PRINT
    -reset false

USER_PUNCH
-headings  sorbCs pH molCs step_no                                  # these are the column headings used for tag names: NB this is the output pH not the input pH
10 sorbedCs=SURF("Cs","illite")
20 if sorbedCs>0 THEN punch sorbedCs, -la("H+"), tot("Cs"), step_no        # NB variable name (used internally) is distinct from column header

SOLUTION 1
  -pH <pHobs>
  -units mol/L
  Ca   0.033                                                          # 0.033 M CaCl2 background electrolyte
  Cl   0.033 charge
  Cs   <Csconcn>                                                     # <Csconcn> from soil_e_imput.dat

SURFACE
-sites_units density
  Illite_ <M1>    98   50                                                         # this is where the max number of sites (updated parameter) is substituted
  -equil 1
  -no_edl

EQUILIBRIUM_PHASES
  Fix_H+ -<pHobs>  NaOH 10                                              # pHobs from the fit data file
  -force_equality true
END

[John Mahoney]

Can you send the observation data file also.  I will try and run it and see what I get.

[FlouretAlex]

Hello,
 here the input datafile.
Thanks for your help.

Cssorbed   Csconcn   pHobs   wt   sim
0.000000218496341   0.000000010946363   7.5   0   1
0.000002232285254   0.000000111725621   7.5   0   1
0.000021709269448   0.000001087103094   7.5   0   1
0.000199251702351   0.000010581012339   7.5   0   1
0.001394322383157   0.000112011070891   7.5   0   1
0.008067309494316   0.000965002238133   7.5   0   1
0.037278023393657   0.009151374603158   7.5   0   1

[John Mahoney]

After several tries at figuring this out and a couple of lost transmittals (all for the good it appears, as I was wrong previously in my thinking). I think I found your problem.

Your speciation error is related to how you define your Illite surface.  Get rid of the underscore in the illite names. I used Illitex for the surface.  It worked eventually, some trial and error.   I actually went back to just PHREEQC and worked on your problem there I included a user_print to track the sorbed Cs and was seeing values of zero listed using the SURF("Cs","surface")  but the Surface summary block was correct, it was seeing sorption. 

Apparently your line   

10 sorbedCs=SURF("Cs","illite")   

 in user_punch was always producing a 0.0 value,  so line 20 could not handle that.
 I tried a few other things that can now be ignored.

HERE is one that sort of work,  well it avoids your Speciation error and plots a figure (Not great) 

SPECIATION
  jobTitle                             "Test fitting: absolute errors (unit weighting)"
  calculationType                      fit
  calculationMethod                    1
  database                PHREEQC_ThermoddemV1.10_06Jun2017.dat
    #debug 3
FIT
  dataFile                             soil_e_imput.dat                       # contains a list of observations and independent variables
  onepass                              FALSE                         # does a separate PHREEQC simulation for each observation -
                                                                     #   slow but easy to set up
  dependentVariableColumnObs           Cssorbed                      # name of column in fit data file containing the observations (dep variable)
  dependentVariableColumnCalc          1                        # name of column in selected output file containing the calcd values of the dep variable
  fitWeightingMethod                   0                             # 0 = unit weights
  fitFiniteDiffStepSize                1.0E-6                        # initial step size for each adjustable parameter
  mainLoopColumn                       sim                           # column header in fit data file for which PHREEQC simulation(s) to use for each observation
  numberOfFitParameters                1
  fitParameterNames                    "M1"
  fitLogParameters                     0                          # 0 = parameters on a linear scale, 1 = log10 parameters before fitting
  fitAdjustableParameters              1                          # 1 = adjustable
  fitParameterValues                   1.0E-04                   # initial values (starting point)

PLOT
  plotTitle                            "Cs sorption on Soil-E (illite only)"
  xtitle                               "Cs_ini (mol/L)"
  ytitle                               "Cs sorbed (mol)"
                                                                     # plot isotherm (x = Csconcn, y = sorbed)
  lines                                calculated                    # y = line from calculated values from out file with column heading = 'calculated'
  points                               observed                      # y = points from observed values from out file with column heading = 'observed'

  lineWidth                            0.4
  lineColor                            red
  labelSize                            0.0                           # no labels on plot
  legendTextSize                       0.0                           # no legend (key)
  pointSize                            4.0
  customXcolumn                        Csconcn                       # x = Csconcn column in out file
png true

CHEMISTRY

PHASES
Fix_H+
  H+ = H+
  log_k 0.0

SURFACE_MASTER_SPECIES
   Illitex  IllitexOH+0.5
   
SURFACE_SPECIES
##### Illite_SOH

IllitexOH+0.5 = IllitexOH+0.5
Log_k 0.0

IllitexOH+0.5 =     IllitexO-0.5  +  H+
Log_k -3.46


IllitexOH+0.5 + Cs+ = IllitexOCs+0.5 + H+
Log_k 5.2

IllitexOH+0.5 + Na+   = IllitexONa+0.5  +  H+
 Log_k -1.8

IllitexOH+0.5 + K+ = IllitexOK+0.5    +  H+
Log_k 0.6

IllitexOH+0.5 + Ca+2  = IllitexOCa+1.5     + H+
Log_k  -5

#   IllitexOH+0.5 + NH4+ = IllitexONH4+0.5    +  H+
#   Log_k 1.5

SELECTED_OUTPUT
-file select.sel
 -high_precision true
 -reset false
 
PRINT
    -selected_output true

USER_PUNCH
-headings  sorbCs pH molCs step_no                                  # these are the column headings used for tag names: NB this is the output pH not the input pH
-start
10 sorbedCs = SURF("Cs","Illitex")
20 if sorbedCs > 0.0  THEN  punch sorbedCs, -la("H+"), tot("Cs"), step_no         # NB variable name (used internally) is distinct from column header
-end
END

SOLUTION 1
  pH <pHobs>
  -units mol/L
  Na 0.001
  Ca   0.033                                                          # 0.033 M CaCl2 background electrolyte
  Cl   0.066 charge
  Cs   <Csconcn>                                                     # <Csconcn> from soil_e_imput.dat

SURFACE
-sites_units density
  Illitex <M1>    98   50                                                         # this is where the max number of sites (updated parameter) is substituted
  -equil 1
  -no_edl

EQUILIBRIUM_PHASES
  Fix_H+ -<pHobs>  NaOH 10                                              # pHobs from the fit data file
  -force_equality true
END
 

I used a slightly different database but there  tends not to be too much for Cs complexation so it worked well enough to solve your problem.
 
This works but the fit is poor  - that is up to you to figure that out.  Run it and see the figure it makes (PNG).  also look at the Standard deviations for your values.   

You might want to fit the log K for the speciation reaction, as well as site density function.   You might want to think about the site density definition and where you start from

I tried a couple of modifications playing with the Surface I got some pretty messed up results, so I will not include them.

[FlouretAlex]

Hello John,

Thank you for your help!
This is great I didn't think about this problem has it was working fine with PHREEQC.
I will look into this with a more complex model now that was a test with a simplified model.
Thank you!
For the problem in the basic statement where line 20produce this type of error "WARNING: Zero divide in BASIC line". I think is due to the initialization step before the sorption calculation that produces this.
You can observe this by plotting the simulation with USER_GRAPH data block with "-initial_solutions TRUE"

Alex