Disconnecting simulations for USER_GRAPH

[MichaelZ20]

Hi David!
I failed in disconnecting simulations of kinetic dissolution for getting a graph of two separate lines starting at zero.
Could you please help?
Thank you in advance.

Code: [Select]

DATABASE C:\phreeqc\database\core10.dat
RATES
Forsterite
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=8.38E+04
12 E1=67206
13 n1=0.470
20 rem neutral solution parameters
21 a2=1.58E+03
22 E2=79000
30 rem base solution parameters
31 a3=1.00E-07
32 E3=56637
33 n2=-0.600
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("forsterite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  # acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               # neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2  # base rate expression
90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
END

SOLUTION 1
units mol/L
Density 1 calculate
pH     7   charge
temp  25
pressure 1
Na 0.5
Cl 0.5
END

KINETICS 1
Forsterite
       -formula  Mg2SiO4  1.0
       -M0    0.2     
       -parms 3  0.1    # total surface area per kg of water (m2/kgw) and the scaling factor
-steps 12 days in 60 steps   
-cvode true
END

INCREMENTAL_REACTIONS true
END

REACTION_PRESSURE 1
40
END

USE SOLUTION 1
REACTION_TEMPERATURE 1
75
USE REACTION_PRESSURE 1
USE KINETICS 1

USER_GRAPH 1
    -headings  temperature  Fosterite-75oC
    -axis_titles            "Day" "mol" "mol"
    -chart_title            "Forsterite Dissolution"
    -initial_solutions      false
    -connect_simulations    false
  -start
10 GRAPH_X SIM_TIME
20 GRAPH_SY 0.2 - KIN("Forsterite")
  -end
  -active true
END

USER_GRAPH 1
-active false
-connect_simulations    false
END

USE SOLUTION 1
REACTION_TEMPERATURE 1
95
USE REACTION_PRESSURE 1
USE KINETICS 1

USER_GRAPH 1
    -active true
    -connect_simulations false
    -headings   temperature  Fosterite-95oC
END

[MichaelZ20]

Dear David, could you please help with my USER_GRAPH problem?

[dlparkhurst]

I don't see the problem. When I run the script I get a straight, red line with square symbols that runs from X=0.2 to 12, Y=0.58e-10 to 3.42e-8, and a separate (higher), straight green line with circle symbols that has the same X range and Y=3.673e-8 to 1.785e-7.

I am using Phreeqc Interactive 3.8.2.16969.

[MichaelZ20]

I want that both lines will start from the same initial mass

[dlparkhurst]

Like this?

Code: [Select]

RATES
Forsterite
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=8.38E+04
12 E1=67206
13 n1=0.470
20 rem neutral solution parameters
21 a2=1.58E+03
22 E2=79000
30 rem base solution parameters
31 a3=1.00E-07
32 E3=56637
33 n2=-0.600
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("forsterite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  # acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               # neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2  # base rate expression
90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
END

SOLUTION 1
units mol/L
Density 1 calculate
pH     7   charge
temp  25
pressure 1
Na 0.5
Cl 0.5
END

KINETICS 1
Forsterite
       -formula  Mg2SiO4  1.0
       -M0    0.2     
       -parms 3  0.1    # total surface area per kg of water (m2/kgw) and the scaling factor
-steps 0 60*0.2   
-cvode true
END

INCREMENTAL_REACTIONS true
END

REACTION_PRESSURE 1
40
END

USE SOLUTION 1
REACTION_TEMPERATURE 1
75
USE REACTION_PRESSURE 1
USE KINETICS 1

USER_GRAPH 1
    -headings  temperature  Fosterite-75oC
    -axis_titles            "Day" "mol" "mol"
    -chart_title            "Forsterite Dissolution"
    -initial_solutions      false
    -connect_simulations    false
  -start
10 GRAPH_X SIM_TIME
20 GRAPH_SY 0.2 - KIN("Forsterite")
  -end
  -active true
END

USER_GRAPH 1
-active false
-connect_simulations    false
END

USE SOLUTION 1
REACTION_TEMPERATURE 1
95
USE REACTION_PRESSURE 1
USE KINETICS 1

USER_GRAPH 1
    -active true
    -connect_simulations false
    -headings   temperature  Fosterite-95oC
END

[MichaelZ20]

Thank you, David, but 1-st line starts at zero, the 2-nd - from the value at what the 1-st simulation ends (3.452e-08).
For comparison, I would like that both simulation will start from the same m0.

[dlparkhurst]

Sorry, didn't check the details. How about this?

Code: [Select]

#DATABASE C:\phreeqc\database\core10.dat
RATES
Forsterite
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=8.38E+04
12 E1=67206
13 n1=0.470
20 rem neutral solution parameters
21 a2=1.58E+03
22 E2=79000
30 rem base solution parameters
31 a3=1.00E-07
32 E3=56637
33 n2=-0.600
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("forsterite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  # acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               # neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2  # base rate expression
90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
END

SOLUTION 1
units mol/L
Density 1 calculate
pH     7   charge
temp  25
pressure 1
Na 0.5
Cl 0.5
END

KINETICS 1
Forsterite
       -formula  Mg2SiO4  1.0
       -M0    0.2     
       -parms 3  0.1    # total surface area per kg of water (m2/kgw) and the scaling factor
-steps 0 2*0.2 #60*0.2   
-cvode true
END

INCREMENTAL_REACTIONS true
END

REACTION_PRESSURE 1
40
END

USE SOLUTION 1
REACTION_TEMPERATURE 1
75
USE REACTION_PRESSURE 1
USE KINETICS 1
USER_PRINT
10 PRINT "Total time: ", TOTAL_TIME
USER_GRAPH 1
    -headings  temperature  Fosterite-75oC
    -axis_titles            "Day" "mol" "mol"
    -chart_title            "Forsterite Dissolution"
    -initial_solutions      false
    -connect_simulations    false
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_SY 0.2 - KIN("Forsterite")
  -end
  -active true
END

USER_GRAPH 1
-active false
-connect_simulations    false
END

USE SOLUTION 1
REACTION_TEMPERATURE 1
95
USE REACTION_PRESSURE 1
USE KINETICS 1

USER_GRAPH 1
    -active true
    -connect_simulations false
    -headings   temperature  Fosterite-95oC
END

[MichaelZ20]

I have checked both with phreeqc3 and phreeqc interactive, but the two lines don't start at the same point. Is not this a bug?

[dlparkhurst]

No, it is not a bug, I am just slow this morning.

KINETICS is updated continuously as it is used, so, yes the second simulation starts with the KINETICS 1 at the end of the first simulation. To avoid this behavior, you need to reinitialize KINETICS 1. I have defined KINETICS 2, and then copied KINETICS 2 to KINETICS 1 before each calculation.

Code: [Select]

RATES
Forsterite
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=8.38E+04
12 E1=67206
13 n1=0.470
20 rem neutral solution parameters
21 a2=1.58E+03
22 E2=79000
30 rem base solution parameters
31 a3=1.00E-07
32 E3=56637
33 n2=-0.600
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("forsterite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  # acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               # neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2  # base rate expression
90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
END

SOLUTION 1
units mol/L
Density 1 calculate
pH     7   charge
temp  25
pressure 1
Na 0.5
Cl 0.5
END

KINETICS 2
Forsterite
       -formula  Mg2SiO4  1.0
       -M0    0.2     
       -parms 3  0.1    # total surface area per kg of water (m2/kgw) and the scaling factor
-steps 0 60*0.2 
-cvode true
END

INCREMENTAL_REACTIONS true
END
COPY kinetics 2 1
END

REACTION_PRESSURE 1
40
END
# initialize kinetics 1
COPY kinetics 2 1
END

USE SOLUTION 1
REACTION_TEMPERATURE 1
75
USE REACTION_PRESSURE 1
USE KINETICS 1
USER_GRAPH 1
    -headings  temperature  Fosterite-75oC
    -axis_titles            "Day" "mol" "mol"
    -chart_title            "Forsterite Dissolution"
    -initial_solutions      false
    -connect_simulations    false
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_SY 0.2 - KIN("Forsterite")
  -end
  -active true
END

USER_GRAPH 1
-active false
-connect_simulations    false
END
# re-initialize kinetics 1
COPY kinetics 2 1
END
USE SOLUTION 1
REACTION_TEMPERATURE 1
95
USE REACTION_PRESSURE 1
USE KINETICS 1

USER_GRAPH 1
    -active true
    -connect_simulations false
    -headings   temperature  Fosterite-95oC
END


[MichaelZ20]

David, thank you very much for your help!