Phreeqc + Pest data fitting problem

[waynehwm]

Hi all,

I am trying to use pest to optimized my reaction contacts for my surface complexation model. My phreeqc + pest works great for optimizing one initial arsenic concentration at different pH conditions, but my script to fit more initial arsenic concentrations curves at the same time is giving me error as below (but it runs fine in phreeqc):

Error condition prevents continued PEST execution:-
Unexpected end to model output file case.prn. Instruction line follows -
"l1 [c31]1:25"

I've been trying to fix this for a while and not luck yet. .......Please let me know if you can find the source of this error, thank you all so much for your help!



My phreeqc script is:

Code: [Select]

ptf @ # PEST fills parameters in between @ @

 
Title sorption of U(VI) on STx-1b
 
SURFACE_MASTER_SPECIES
Hfo_s  Hfo_sOH
Hfo_w  Hfo_wOH
 
SURFACE_SPECIES
Hfo_wOH  + H+ = Hfo_wOH2+
log_k  7.29    # = pKa1,int

Hfo_wOH = Hfo_wO- + H+
log_k  -8.93   # = -pKa2,int
 
Hfo_wOH + AsO4-3 + 3H+ = Hfo_wH2AsO4 + H2O
log_k   @k1 @

Hfo_wOH + AsO4-3 + 2H+ = Hfo_wHAsO4- + H2O
log_k   @k2 @

Hfo_wOH + AsO4-3 = Hfo_wOHAsO4-3
log_k   @k3 @
 
 
SURFACE 1
Hfo_wOH 7e-5 2.1 5

 
Phases; fix_pH; H+ = H+; log_k 0
fix_pe; e- = e- ; log_k 0
Solution 1; units mmol/l; temp 25; pH 2.5; As(+5) 5.78E-02; Na 50; N(+5) 50
END
 
SELECTED_OUTPUT; -file case.prn; -reset false; -high_p true
USER_PUNCH
 -start
 10 Qadsorb = mol("Hfo_wH2AsO4") + mol("Hfo_wHAsO4-") + mol("Hfo_wOHAsO4-3")
 20 Punch Qadsorb/(TOT("As(+5)")+TOT("As(+3)")+Qadsorb)*100
 -end
 
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -2.50 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -3.33 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -4.71 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -4.74 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.25 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.57 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.90 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -6.15 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -7.35 NaOH 10; end
use solution 1; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -8.00 NaOH 10; end
 
Solution 2; units mmol/l; temp 25; pH 2.43; As(+5) 1.18E-02; Na 50; N(+5) 50
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -2.43 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -3.48 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -3.72 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -4.29 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -4.70 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.67 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -6.70 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -7.04 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -7.14 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -8.08 NaOH 10; end
use solution 2; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -8.22 NaOH 10; end
 
Solution 3; units mmol/l; temp 25; pH 2.53; As(+5) 1.44E-03; Na 50; N(+5) 50
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -2.53 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -2.97 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -4.57 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -4.70 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.08 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.25 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.48 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -5.90 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -6.00 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -6.30 NaOH 10; end
use solution 3; use surface 1; EQUILIBRIUM_PHASES 1; CO2(g) -3.5; fix_pH -7.10 NaOH 10; end
My pest pcf file script is:

Code: [Select]

pcf
* control data
restart estimation
3 32 1 0 3
# no of parameters,observations, group parameters, prior infos , groups
1 1 single point 1 0 0
# no of template files, no of instruction files, ….
10.0 2.0 0.3 0.03 10
3.0 3.0 0.001
0.1
30 0.001 3 3 0.01 3
1 1 1

* parameter groups
k1 relative 1e-2 0.000001 switch 2 parabolic

* parameter data
k1 none relative 29.31 3.0 50 k1 1.0 0.0 1
#param_name,type, initial, min, max, ..….
k2 none relative 23.51 3.0 50 k1 1.0 0.0 1
k3 none relative 10.58 1.0 30 k1 1.0 0.0 1

* observation groups
group_1
group_2
group_3

* observation data
c1 99.98 1 group_1
# obs_name, value, weight, group no
c2 99.84 1 group_1
c3 79.81 1 group_1
c4 78.62 1 group_1
c5 57.02 1 group_1
c6 46.28 1 group_1
c7 38.7 1 group_1
c8 37.56 1 group_1
c9 32.98 1 group_1
c10 32.54 1 group_1
c11 97.36 1 group_2
c12 99.52 1 group_2
c13 99.26 1 group_2
c14 99.46 1 group_2
c15 82.89 1 group_2
c16 72.64 1 group_2
c17 60.11 1 group_2
c18 57.49 1 group_2
c19 54.52 1 group_2
c20 38.37 1 group_2
c21 35.32 1 group_2
c22 94.65 1 group_3
c23 98.68 1 group_3
c24 98.40 1 group_3
c25 98.40 1 group_3
c26 97.92 1 group_3
c27 96.32 1 group_3
c28 90.42 1 group_3
c29 87.36 1 group_3
c30 85.90 1 group_3
c31 63.44 1 group_3
c32 52.12 1 group_3

* model command line
phreeqc.bat case.phrq case.out C:/phreeqc/database/wateq4f.dat scr.out
* model input/output
case.tpl case.phrq
case.ins case.prn
* prior information
And my pest pif script is:

Code: [Select]

pif @ # read line 1 or 2 (l1 or l2); data to be read between position 1 and 25 i
l2 [c1]1:25
l1 [c2]1:25
l1 [c3]1:25
l1 [c4]1:25
l1 [c5]1:25
l1 [c6]1:25
l1 [c7]1:25
l1 [c8]1:25
l1 [c9]1:25
l1 [c10]1:25
l2 [c11]1:25
l1 [c12]1:25
l1 [c13]1:25
l1 [c14]1:25
l1 [c15]1:25
l1 [c16]1:25
l1 [c17]1:25
l1 [c18]1:25
l1 [c19]1:25
l1 [c20]1:25
l1 [c21]1:25
l2 [c22]1:25
l1 [c23]1:25
l1 [c24]1:25
l1 [c25]1:25
l1 [c26]1:25
l1 [c27]1:25
l1 [c28]1:25
l1 [c29]1:25
l1 [c30]1:25
l1 [c31]1:25
l1 [c32]1:25
Experimental Data:

Code: [Select]

99.98
99.84
79.81
78.62
57.02
46.28
38.70
37.56
32.98
32.54

97.36
99.52
99.26
99.46
82.89
72.64
60.11
57.49
54.52
38.37
35.27

94.65
98.68
98.40
98.40
97.92
96.32
90.42
87.36
85.90
63.44
52.12

Model outpu:

Code: [Select]

 
 8.351437813887e+01
  8.532573646886e+01
  7.496505695011e+01
  7.434860128424e+01
  6.172737656355e+01
  5.247644165240e+01
  4.286862944892e+01
  3.611249822419e+01
  1.589915648597e+01
  9.882746002669e+00
  2.537659649104e-92
  9.887817208241e+01
  9.886365219428e+01
  9.863620542418e+01
  9.761401090096e+01
  9.623998565600e+01
  8.851483430185e+01
  6.936248827317e+01
  6.020162677238e+01
  5.724112119731e+01
  2.628890062291e+01
  2.181124686832e+01
  2.079472718407e-91
  9.933216210958e+01
  9.938473142527e+01
  9.944843698765e+01
  9.950172736885e+01
  9.960341269357e+01
  9.962407409635e+01
  9.962991928194e+01
  9.956771468829e+01
  9.953411084012e+01
  9.935764608979e+01
  9.667402849929e+01

[dlparkhurst]

Looks like you only have 30 simulations (lines with "USE solution ..."), and you are trying to read results for 32.

[waynehwm]

Looks like you only have 30 simulations (lines with "USE solution ..."), and you are trying to read results for 32.

Thank you dlparkhurst!
Yes, indeed I made that mistake. I've correct it along with a little more modifications on the script. Now the models runs, but the output has two additional data point "2.54E-92", "2.08E-91" in between data groups. Do you know how to remove these data points?

Many thanks!

[dlparkhurst]

I need to see the script to know what you did, but there will be a line for every SOLUTION definition, and for every "USE solution ..." line from the point that SELECTED_OUTPUT and USER_PUNCH are defined.

You are probably performing a SOLUTION calculation, (or possibly another calculation like surface initialization, -equil). You may be able to move the calculations you don't want to show up in the output before you define SELECTED_OUTPUT, or you can disable selected output while the additional calculations are performed (PRINT; -selected_output false; END and PRINT; -selected_output true; END). Liberal use of END will help separate calculations for which SELECTED_OUTPUT is turned on and off.

(Alternatively, you can adjust your PEST instruction file to skip the lines you don't want.)

[waynehwm]

I need to see the script to know what you did, but there will be a line for every SOLUTION definition, and for every "USE solution ..." line from the point that SELECTED_OUTPUT and USER_PUNCH are defined.

You are probably performing a SOLUTION calculation, (or possibly another calculation like surface initialization, -equil). You may be able to move the calculations you don't want to show up in the output before you define SELECTED_OUTPUT, or you can disable selected output while the additional calculations are performed (PRINT; -selected_output false; END and PRINT; -selected_output true; END). Liberal use of END will help separate calculations for which SELECTED_OUTPUT is turned on and off.

(Alternatively, you can adjust your PEST instruction file to skip the lines you don't want.)

Thank you dlparkhurst, I think I know what you meant, so there is a line for my "Solution 2 XXXX" and "Solution 3 XXX" statement, and we can used selected_output to by pass these lines?

I have updated the previous post with my latest script and result. Can you show me how to use selected_output to by pass these lines?

Thank you so much!

[waynehwm]

I need to see the script to know what you did, but there will be a line for every SOLUTION definition, and for every "USE solution ..." line from the point that SELECTED_OUTPUT and USER_PUNCH are defined.

You are probably performing a SOLUTION calculation, (or possibly another calculation like surface initialization, -equil). You may be able to move the calculations you don't want to show up in the output before you define SELECTED_OUTPUT, or you can disable selected output while the additional calculations are performed (PRINT; -selected_output false; END and PRINT; -selected_output true; END). Liberal use of END will help separate calculations for which SELECTED_OUTPUT is turned on and off.

(Alternatively, you can adjust your PEST instruction file to skip the lines you don't want.)

Thank you dlparkhurst, I think I know what you meant, so there is a line for my "Solution 2 XXXX" and "Solution 3 XXX" statement, and we can used selected_output to by pass these lines?

I have updated the previous post with my latest script and result. Can you show me how to use selected_output to by pass these lines?

Thank you so much!

I think I got it now. I moved the "Solution XXXXX" lines to the Phases section, and the extra data points are removed. Thank you for your help dlparkhurst!!