Adding solid phases during kinetic reactions

[Rensburg]

Dear PHREEQC researchers,
   I am currently learning to use PHREEQC and have encountered some confusion, so I'm here seeking help.

  ①To simulate a scenario where a solution undergoes KINETICS reactions for a period under specific solid and gas phase conditions, and then the solid and gas phase environment is modified (e.g., by adding Al(OH)3(a) to the solution) before continuing with further KINETICS reactions, I modified ex9 as an example as following. I would like to ask if my approach is correct, or if there is a more suitable method. If the reaction rate involves KIN("") or M, should I accurately obtain 'm' and define it correctly in KINETICS 2?

Code: [Select]

TITLE Example 9.--Kinetically controlled oxidation of ferrous
                  iron. Decoupled valence states of iron.
SOLUTION_MASTER_SPECIES
Fe_di              Fe_di+2    0.0     Fe_di              55.847
Fe_tri             Fe_tri+3   0.0     Fe_tri             55.847
SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
        log_k   0.0
Fe_tri+3 = Fe_tri+3
        log_k   0.0
#
# Fe+2 species
#
Fe_di+2 + H2O = Fe_diOH+ + H+
        log_k   -9.5
        delta_h 13.20   kcal
#
#... and also other Fe+2 species
#
Fe_di+2 + Cl- = Fe_diCl+
        log_k   0.14
Fe_di+2 + CO3-2 = Fe_diCO3
        log_k   4.38
Fe_di+2 + HCO3- = Fe_diHCO3+
        log_k   2.0
Fe_di+2 + SO4-2 = Fe_diSO4
        log_k   2.25
        delta_h 3.230   kcal
Fe_di+2 + HSO4- = Fe_diHSO4+
        log_k   1.08
Fe_di+2 + 2HS- = Fe_di(HS)2
        log_k   8.95
Fe_di+2 + 3HS- = Fe_di(HS)3-
        log_k   10.987
Fe_di+2 + HPO4-2 = Fe_diHPO4
        log_k   3.6
Fe_di+2 + H2PO4- = Fe_diH2PO4+
        log_k   2.7
Fe_di+2 + F- = Fe_diF+
        log_k   1.0
#
# Fe+3 species
#
Fe_tri+3 + H2O = Fe_triOH+2 + H+
        log_k   -2.19
        delta_h 10.4    kcal
#
#... and also other Fe+3 species
#
Fe_tri+3 + 2 H2O = Fe_tri(OH)2+ + 2 H+
        log_k   -5.67
        delta_h 17.1    kcal
Fe_tri+3 + 3 H2O = Fe_tri(OH)3 + 3 H+
        log_k   -12.56
        delta_h 24.8    kcal
Fe_tri+3 + 4 H2O = Fe_tri(OH)4- + 4 H+
        log_k   -21.6
        delta_h 31.9    kcal
2 Fe_tri+3 + 2 H2O = Fe_tri2(OH)2+4 + 2 H+
        log_k   -2.95
        delta_h 13.5    kcal
3 Fe_tri+3 + 4 H2O = Fe_tri3(OH)4+5 + 4 H+
        log_k   -6.3
        delta_h 14.3    kcal
Fe_tri+3 + Cl- = Fe_triCl+2
        log_k   1.48
        delta_h 5.6     kcal
Fe_tri+3 + 2 Cl- = Fe_triCl2+
        log_k   2.13
Fe_tri+3 + 3 Cl- = Fe_triCl3
        log_k   1.13
Fe_tri+3 + SO4-2 = Fe_triSO4+
        log_k   4.04
        delta_h 3.91    kcal
Fe_tri+3 + HSO4- = Fe_triHSO4+2
        log_k   2.48
Fe_tri+3 + 2 SO4-2 = Fe_tri(SO4)2-
        log_k   5.38
        delta_h 4.60    kcal
Fe_tri+3 + HPO4-2 = Fe_triHPO4+
        log_k   5.43
        delta_h 5.76    kcal
Fe_tri+3 + H2PO4- = Fe_triH2PO4+2
        log_k   5.43
Fe_tri+3 + F- = Fe_triF+2
        log_k   6.2
        delta_h 2.7     kcal
Fe_tri+3 + 2 F- = Fe_triF2+
        log_k   10.8
        delta_h 4.8     kcal
Fe_tri+3 + 3 F- = Fe_triF3
        log_k   14.0
        delta_h 5.4     kcal
PHASES
Goethite
        Fe_triOOH + 3 H+ = Fe_tri+3 + 2 H2O
        log_k   -1.0
END
SOLUTION 1
        pH  7.0
        pe 10.0
        O2(g) -0.67
        Fe_di  0.1
        Na  10.
        Cl  10.  charge
EQUILIBRIUM_PHASES 1
        O2(g)           -0.67
RATES
Fe_di_ox
  -start
  10  Fe_di = TOT("Fe_di")
  20  if (Fe_di <= 0) then goto 200
  30  p_o2 = SR("O2(g)")
  40  moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2 * p_o2) * Fe_di * TIME
  200 SAVE moles
  -end
KINETICS 1
Fe_di_ox
        -formula  Fe_di  -1.0  Fe_tri  1.0
        -steps 100 400 3100 10800 21600 5.04e4
        -step_divide 1e-4
INCREMENTAL_REACTIONS true
SAVE SOLUTION 2
SELECTED_OUTPUT
        -file ex9_1.sel
        -reset false
USER_PUNCH
        -headings  Fe(2)  Fe(3)  pH  si_goethite
  10 PUNCH  TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6, -LA("H+"), SI("Goethite")
END

USE SOLUTION 2
EQUILIBRIUM_PHASES 2
       Al(OH)3(a)    0  10
        O2(g)           -0.67
KINETICS 2
Fe_di_ox
        -formula  Fe_di  -1.0  Fe_tri  1.0
        -steps 8.64e4 1.728e5 1.728e5 1.728e5 1.728e5
        -step_divide 1e-4
INCREMENTAL_REACTIONS true
SELECTED_OUTPUT
        -file ex9_2.sel
        -reset false
USER_PUNCH
        -headings Fe(2)  Fe(3)  pH  si_goethite si_Calcite
  10 PUNCH TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6, -LA("H+"),SI("Goethite"),SI("Calcite")
END
  ②If I want to implement this operation of adding a solid phase during kinetic reactions in phreeqcrm, what should I do? Is it necessary to first run the first part of the KINETICS (ex9_1.pqi), use GetConcentrations() to obtain the results, and then redefine the ex9_2.pqi file to rerun the code, as shown in the input below? Is there a better way to achieve this? Because this method is not very intelligent, and I don't know how to obtain pe and pH to accurately define the new solution.

Code: [Select]

#ex9_1.pqi,and phreeqcrm's step is 100 400 3100 10800 21600 5.04e4
SOLUTION_MASTER_SPECIES
SOLUTION_SPECIES
PHASES
END
SOLUTION 1
        pH  7.0
        pe 10.0
        O2(g) -0.67
        Fe_di  0.1
        Na  10.
        Cl  10.  charge
EQUILIBRIUM_PHASES 1
        O2(g)           -0.67
RATES
Fe_di_ox
  -start
  10  Fe_di = TOT("Fe_di")
  20  if (Fe_di <= 0) then goto 200
  30  p_o2 = SR("O2(g)")
  40  moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2 * p_o2) * Fe_di * TIME
  200 SAVE moles
  -end
KINETICS 1
Fe_di_ox
        -formula  Fe_di  -1.0  Fe_tri  1.0
        -step_divide 1e-4
SELECTED_OUTPUT
        -file ex9_1.sel
        -reset false
USER_PUNCH
        -headings  Fe(2)  Fe(3)  pH  si_goethite
  10 PUNCH  TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6, -LA("H+"), SI("Goethite")
END


Code: [Select]

#ex9_2.pqi
SOLUTION_MASTER_SPECIES
SOLUTION_SPECIES
PHASES
END
SOLUTION 1 #new value from getconcentrations()
        pH  7.0
        pe 10.0
        O  new_value
        Fe_di  new_value
        …
EQUILIBRIUM_PHASES 1
        Al(OH)3(a)    0   10
        O2(g)           -0.67
RATES

KINETICS 1
Fe_di_ox
        -formula  Fe_di  -1.0  Fe_tri  1.0
        -step_divide 1e-4
SELECTED_OUTPUT
        -file ex9_2.sel
        -reset false
USER_PUNCH
        -headings  Fe(2)  Fe(3)  pH  si_goethite
  10 PUNCH  TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6, -LA("H+"), SI("Goethite")
END
  ③I am using phreeqc.dat for both the software and phreeqcrm. Why does phreeqcrm display the warning: “WARNING: Could not find element in database, O2(g). Concentration is set to zero.”

[dlparkhurst]

I think it is best if you use RUN_CELLS to do this type of calculation. RUN_CELLS will react all reactants with the same user number. In the example below, all the reactants numbered 1 will be reacted, including solution 1, kinetics 1, and equilibrium phases 1.

The KINETICS definition has multiple time steps, so the reaction will be integrated over those time steps, processing output (SELECTED_OUTPUT and(or) USER_GRAPH) at each time step. After the first RUN_CELLS, solution 1, equilibrium phases 1, and kinetics 1 are in the state after the last reaction step.

You can then redefine equilibrium phases 1 to include Al(OH)3 before another set of reactions with RUN_CELLS.

There is not an easy way to define multiple, unequal time steps with RUN_CELLS. You do not want to use the same unequal time steps of the original KINETICS definition, so you can define RUN_CELLS for each additional time increment, or I have finessed 5 equal time steps by defining REACTION_TEMPERATURE to have five time steps and the RUN_CELLS; -time_step is set to the total time of the additional simulations. This will result in five additional calculations, each with a time step of 2e5 seconds (and temperature of 25 C).

This may not be exactly what you want, but it is close, and you can modify as needed.

You must define elements or redox states in SOLUTION. O2(g) is not an element, but O(0) is the redox state of dissolved oxygen. You can then specify the concentration of O(0) to be that which is in equilibrium with O2(g) at a partial pressure of 10^-0.67 as is done in the script.

Code: [Select]

TITLE Example 9.--Kinetically controlled oxidation of ferrous
                  iron. Decoupled valence states of iron.
SOLUTION_MASTER_SPECIES
Fe_di              Fe_di+2    0.0     Fe_di              55.847
Fe_tri             Fe_tri+3   0.0     Fe_tri             55.847
SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
        log_k   0.0
Fe_tri+3 = Fe_tri+3
        log_k   0.0
#
# Fe+2 species
#
Fe_di+2 + H2O = Fe_diOH+ + H+
        log_k   -9.5
        delta_h 13.20   kcal
#
#... and also other Fe+2 species
#
Fe_di+2 + Cl- = Fe_diCl+
        log_k   0.14
Fe_di+2 + CO3-2 = Fe_diCO3
        log_k   4.38
Fe_di+2 + HCO3- = Fe_diHCO3+
        log_k   2.0
Fe_di+2 + SO4-2 = Fe_diSO4
        log_k   2.25
        delta_h 3.230   kcal
Fe_di+2 + HSO4- = Fe_diHSO4+
        log_k   1.08
Fe_di+2 + 2HS- = Fe_di(HS)2
        log_k   8.95
Fe_di+2 + 3HS- = Fe_di(HS)3-
        log_k   10.987
Fe_di+2 + HPO4-2 = Fe_diHPO4
        log_k   3.6
Fe_di+2 + H2PO4- = Fe_diH2PO4+
        log_k   2.7
Fe_di+2 + F- = Fe_diF+
        log_k   1.0
#
# Fe+3 species
#
Fe_tri+3 + H2O = Fe_triOH+2 + H+
        log_k   -2.19
        delta_h 10.4    kcal
#
#... and also other Fe+3 species
#
Fe_tri+3 + 2 H2O = Fe_tri(OH)2+ + 2 H+
        log_k   -5.67
        delta_h 17.1    kcal
Fe_tri+3 + 3 H2O = Fe_tri(OH)3 + 3 H+
        log_k   -12.56
        delta_h 24.8    kcal
Fe_tri+3 + 4 H2O = Fe_tri(OH)4- + 4 H+
        log_k   -21.6
        delta_h 31.9    kcal
2 Fe_tri+3 + 2 H2O = Fe_tri2(OH)2+4 + 2 H+
        log_k   -2.95
        delta_h 13.5    kcal
3 Fe_tri+3 + 4 H2O = Fe_tri3(OH)4+5 + 4 H+
        log_k   -6.3
        delta_h 14.3    kcal
Fe_tri+3 + Cl- = Fe_triCl+2
        log_k   1.48
        delta_h 5.6     kcal
Fe_tri+3 + 2 Cl- = Fe_triCl2+
        log_k   2.13
Fe_tri+3 + 3 Cl- = Fe_triCl3
        log_k   1.13
Fe_tri+3 + SO4-2 = Fe_triSO4+
        log_k   4.04
        delta_h 3.91    kcal
Fe_tri+3 + HSO4- = Fe_triHSO4+2
        log_k   2.48
Fe_tri+3 + 2 SO4-2 = Fe_tri(SO4)2-
        log_k   5.38
        delta_h 4.60    kcal
Fe_tri+3 + HPO4-2 = Fe_triHPO4+
        log_k   5.43
        delta_h 5.76    kcal
Fe_tri+3 + H2PO4- = Fe_triH2PO4+2
        log_k   5.43
Fe_tri+3 + F- = Fe_triF+2
        log_k   6.2
        delta_h 2.7     kcal
Fe_tri+3 + 2 F- = Fe_triF2+
        log_k   10.8
        delta_h 4.8     kcal
Fe_tri+3 + 3 F- = Fe_triF3
        log_k   14.0
        delta_h 5.4     kcal
PHASES
Goethite
        Fe_triOOH + 3 H+ = Fe_tri+3 + 2 H2O
        log_k   -1.0
RATES
Fe_di_ox
  -start
  10  Fe_di = TOT("Fe_di")
  20  if (Fe_di <= 0) then goto 200
  30  p_o2 = SR("O2(g)")
  40  moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2 * p_o2) * Fe_di * TIME
  200 SAVE moles
  -end
END
SOLUTION 1
        pH  7.0
        pe 10.0
        O(0) 1 O2(g) -0.67
        Fe_di  0.1
        Na  10.
        Cl  10.  charge
END
EQUILIBRIUM_PHASES 1
        O2(g)           -0.67 10
END
KINETICS 1
Fe_di_ox
        -formula  Fe_di  -1.0  Fe_tri  1.0
        -steps 100 400 3100 10800 21600 5.04e4
        -step_divide 1e-4
END
USER_GRAPH 1
    -headings               time Fe_di Fe_tri
    -axis_titles            "Days" "Molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 86400
20 GRAPH_Y TOT("Fe_di"), TOT("Fe_tri")
  -end
    -active                 true
END
RUN_CELLS
-cell 1
END
EQUILIBRIUM_PHASES 1
       Al(OH)3(a)    0  10
        O2(g)           -0.67
END
RUN_CELLS
-cell 1
-time_step 1e6
REACTION_TEMPERATURE 1
25 25 in 5 steps
END


[Rensburg]

  Thank you for your help. I have just looked into the use of RUN_CELLS, and it nicely preserves the previous state. However, this method seems to only work in the software interface, not in phreeqcrm. But you have inspired me, and I have come up with a new idea, though I'm not sure if it's reasonable as following. 

Code: [Select]

RunFile (ex9_1.pqi)
RUN_CELLS()
input = "EQUILIBRIUM_PHASES 1 
       Al(OH)3(a)    0  10 
        O2(g)           -0.67"
RunString(input)
RUN_CELLS()
  Please don't mind that the code I wrote for running Phreeqcrm is quite rudimentary.

[dlparkhurst]

I did not know you were using PhreeqcRM. PhreeqcRM effectively uses RUN_CELLS for each cell for each time step, so you generally would not use RUN_CELLS within PhreeqcRM.

So, are you trying to define initial conditions with Phreeqc (with SOLUTION, EQUILIBRIUM_PHASES, and KINETICS) that are then transferred to PhreeqcRM, or are you trying to adjust the EQUILIBRIUM_PHASES definition in PhreeqcRM?

[Rensburg]

In the future, I will only use PhreeqcRM and make adjustments directly in PhreeqcRM.

[dlparkhurst]

In that case, you can run the new definition of EQUILIBRIUM_PHASES for a cell with the RUN_STRING (or RUN_FILE) method of PhreeqcRM. In this case, you cannot save the old moles of other equilibrium phases.

Another option is to use the EQUILIBRIUM_PHASES_MODIFY keyword described in the appendix of the PHREEQC manual to add a new phase to a cell. The amounts of the other phases would then be retained. You would construct this keyword data block and, again, run it with the RUN_STRING (or RUN_FILE) method of PhreeqcRM.

Code: [Select]

SOLUTION
END
EQUILIBRIUM_PHASES 1
O2(g) -0.7 0.001
END
RUN_CELLS
-cell 1
END
EQUILIBRIUM_PHASES_MODIFY 1
-component Gibbsite
# -moles 0
#     -si    0
END
RUN_CELLS
-cell 1
END

[Rensburg]

  I tried writing both the PHREEQC software version and the PhreeqcRM version, and both ran , yielding nearly identical results. However, I'm not entirely sure if I understood your instructions correctly, so I'd like to discuss it with you to confirm whether my approach is correct.

  Below is the general idea of my modified PHREEQC software version

Code: [Select]

TITLE Example 9.--Kinetically controlled oxidation of ferrous
                  iron. Decoupled valence states of iron.
SOLUTION_MASTER_SPECIES
Fe_di              Fe_di+2    0.0     Fe_di              55.847
Fe_tri             Fe_tri+3   0.0     Fe_tri             55.847
SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
        log_k   0.0
Fe_tri+3 = Fe_tri+3
        log_k   0.0
…
PHASES
Goethite
        Fe_triOOH + 3 H+ = Fe_tri+3 + 2 H2O
        log_k   -1.0
RATES
Fe_di_ox
  -start
  10  Fe_di = TOT("Fe_di")
  20  if (Fe_di <= 0) then goto 200
  30  p_o2 = SR("O2(g)")
  40  moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2 * p_o2) * Fe_di * TIME
  200 SAVE moles
  -end
END

SOLUTION 1
        pH  7.0
        pe 10.0
        O(0) 1 O2(g) -0.67
        Fe_di  0.1
        Na  10.
        Cl  10.  charge
EQUILIBRIUM_PHASES 1
        O2(g)           -0.67
INCREMENTAL_REACTIONS true
KINETICS 1
Fe_di_ox
        -formula  Fe_di  -1.0  Fe_tri  1.0
        -steps 100 400 3100 10800 21600 5.04e4
        -step_divide 1e-4
END

RUN_CELLS
-cell 1
END

EQUILIBRIUM_PHASES_MODIFY 1
   -component Al(OH)3(a)
   -moles 10
   -si    0
END

RUN_CELLS
-cell 1
-time_step 518400
END

   Below is the general idea of the PhreeqcRM version code .The `ex9_1.pqi` file is as I mentioned earlier. Additionally, I wrote the code based on the `SimpleAdvect.py` example. I had to delete the `status = phreeqc_rm.RunString(True, False, False, "DELETE; -all")` statement to successfully run `status = simulator.phreeqc_rm.RunString(True, True, True, modify_string)`. I'm not sure if deleting this line will have any adverse consequences. Currently, I don't see any difference in the calculation results.

Code: [Select]

# Initialization (many definitions omitted here)
    status = phreeqc_rm.RunFile(True, True, True, "ex9_1.pqi")
    # status = phreeqc_rm.RunString(True, False, False, "DELETE; -all")
    time = 0.0
    time_step = 0.0
    status = phreeqc_rm.SetTime(time)
    status = phreeqc_rm.SetTimeStep(time_step)
    status = phreeqc_rm.RunCells()

# Kinetic reactions
    time = 0.0
    time_steps = [100, 400, 3100, 10800, 21600, 50400]
    for time_step in time_steps:
        time += time_step
        status = phreeqc_rm.SetTime(time)
        status = phreeqc_rm.SetTimeStep(time_step)
        status = phreeqc_rm.SetConcentrations(concentrations)
        status = phreeqc_rm.RunCells()
        concentrations = phreeqc_rm.GetConcentrations()

# Modify solid phases
    input = "EQUILIBRIUM_PHASES_MODIFY 1
        -component Al(OH)3(a)
        -moles 10
        -si    0"
    status = phreeqc_rm.RunString(True, False, True, input)

# Continue kinetic reactions
    time = 86400
    time_steps = [8.64e4, 8.64e4,8.64e4,  8.64e4,8.64e4,  8.64e4]
    for time_step in time_steps:
        time += time_step
        status = phreeqc_rm.SetTime(time)
        status = phreeqc_rm.SetTimeStep(time_step)
        status = phreeqc_rm.SetConcentrations(concentrations)
        status = phreeqc_rm.RunCells()
        concentrations = phreeqc_rm.GetConcentrations()

# Clean up
    status = phreeqc_rm.CloseFiles()
    status = phreeqc_rm.MpiWorkerBreak()
  Thank you for your help very much .

[dlparkhurst]

I will just say that you can use EQUILIBRIUM_PHASES_MODIFY in Phreeqc or PhreeqcRM. You will have to decide if you are doing the calculations the way you want.

I will also say that changing the equilibrium-phase assemblage mid-calculation is fairly rare. PhreeqcRM is used most often for groundwater models, and changes in the minerals in the subsurface are unlikely, except for minerals that might precipitate because of changes in water chemistry. Usually, this situation is handled by including minerals in the phase assemblage (with zero moles) that could precipitate if they become saturated.

[Rensburg]

   Thank you for your reminder. My internet connection seemed unstable last night, so I couldn't reply to you in time.I will continue to learn more about this and conduct reasonability analysis.
   
   Finally, I would like to consult you on one matter: In PhreeqcRM, I can use RunString() and RunFile() to define the initial state the first time, and use RunCells() to assign them to respective cells. If I subsequently need to add, remove, or modify the initial definitions, is it correct that I only need to use RunString() or RunFile() again? Do I need to set time_step=0 and run RunCells() again to distribute them to the grid? Additionally, would this cause any of the original definitions to become invalid?

  "PhreeqcRM has three IPhreeqc instances: one for processing initial conditions (Initial Phreeqc), one for storing and running the contents of each cell (Worker), and one we don't need to worry about." When I first define the initial conditions, I should use `phreeqc_rm.RunFile(True, True, True, pqi_file)`. For subsequent additions of new definitions, what should be the correct combination of parameters? `phreeqc_rm.Runstring(False, True, True, pqi_file)`?

[dlparkhurst]

In the normal use of PhreeqcRM in reactive transport models, the initial conditions are set and you do not change the reactants mid calculation.

To get initial concentrations for the model cells, the strategy is to use the initial IPhreeqc instance to define a set of solutions and other reactants (EQUILIBRIUM_PHASES for example) that will be used as initial conditions and boundary conditions for your simulation. Usually, these solutions and reactants are defined by running a Phreeqc input file in the initial IPhreeqc instance.  (However, if you change the database definitions--SOLUTION_MASTER_SPECIES, SOLUTION_SPECIES, PHASES, EXCHANGE_MASTER_SPECIES, and others--you must run these definitions through the workers as well.) Next you distribute the initial solutions and other reactants to the model cells with InitialPhreeqc2Module. You then use RunCells (with no kinetics, 0 time step) to get the initial condition solutions and reactants in equilibrium. After that, GetConcentrations will give you the concentrations that you need to do the first transport step.

For boundary conditions, you would use one or more of the solutions in the initial IPhreeqc instance. You can obtain the concentrations needed for boundary conditions with the InitialPhreeqc2Concentrations method.

So, you transport to get new concentrations, submit those concentrations to PhreeqcRM with SetConcentrations, RunCells to calculate new concentrations as a result of reactions, and finally begin the process again with GetConcentrations and a new transpore calculation.

So, normally, you do not need to reset concentrations or reactants outside of this sequence. You can look at SimpleAdvect examples in the Tests subdirectory for the basic sequence of calculations.

Perhaps you are trying to do some other type of modeling. You will need to explain what you are trying to do. It is possible to use RunFile or RunString with the initial phreeqc instance (or simply use reactants that were previously defined in the initial phreeqc instance) and call InitialPhreeqc2Module at an intermediate point in the calculation, but it is not common.

[Rensburg]

  I have carefully studied the SimpleAdvect.py example, but it does not cover situations like suddenly wanting to change my EQUILIBRIUM_PHASES midway. How should I implement this? Can I simply write a relevant PHREEQC statement and use RunString?
   In the future, I want to simulate a surface water environment where, at time A, a certain solid phase (immobile, not transported with water flow) is introduced at a specific location and undergoes kinetic reactions for a period, and at time B, another solid phase is introduced at the same location for further kinetic reactions (also immobile). I think it is unreasonable to define all solid phases upfront, so I may need to add them midway using RunString(). Therefore, I would like to know if there is a standard workflow for such mid-simulation additions.

[dlparkhurst]

You can use the methods we have discussed.

I think you want to generate the equilibrium_phases_modify block(s), probably within your program that calls phreeqcrm, use RunString on the workers, and continue with the calculations.


To remove a phase, you may need to set the moles to zero, and set dissolve only.

[Rensburg]

 Bother you again.   I think I didn't express my doubt clearly. I mean that after the first definition of RunFile(), InitialPhreeqc2Module needs to be allocated.
  Then, when I need to supplement the definition for the next RunString, do I need to allocate InitialPhreeqc2Module again and balance it with Runcells() + time_step = 0?
  Also, for this supplementary definition, should I use RunString(True, False, False, input) or RunString(True, True, False ,input)? I try both, but it seems no difference in the result.
  Thank you.

[dlparkhurst]

If you define the EQUILIBRIUM_PHASES_MODIFY string and RunString(true, false, false, ...), then the equilibrium phases are modified and ready to use in the next RunCells calculation.

If you use RunString(false, true, false, ...), then you need first to have a definition in the InitialPhreeqc instance with the same cell number to modify, and then you would need to transfer the new EQUILIBRIUM_PHASES compositions to the workers with InitialPhreeqc2Module or InitialEquilibriumPhases2Module. You need to take care that you only transfer the equilibrium phases that you want to modify.

If you use RunString(true, true, false, ...), you are not accomplishing anything in the InitialPhreeqc instance that is related to your model calculation unless you also run InitialPhreeqc2Module (or InitialEquilibriumPhases2Module), in which case you will overwrite the changes made in common with RunString(true...).

Often, you only have one or two equilibrium phases definitions in the InitialPhreeqc instance, and these are then distributed to the cells of the model, which reside in the workers; the same definition from the InitialPhreeqc instance can be applied to multiple model cells. In this case, you would likely not have a definition in the InitialPhreeqc instance for each cell number., only one or two that apply to different zones of the model.

Whether you do the work in the InitialPhreeqc instance, or the workers, you do not need to make an extra call to RunCells. The equilibrium phases will be ready for the call to RunCells in the next time step.

[Rensburg]

  Thank you for your reply. This time I am very clear about what I should do. Once again, I am extremely grateful.