Using PhreeqPy to interrupt a model and calculate new values for subsequent simulations

[Tom Kibblewhite]

Hi Mike, David, et al.

In the PhreeqC input script attached (after some kinetic mineral reactions to form my initial solution), I remove 5% of the water volume by evaporation (“REACTION”) and refill with the same volume (“MIX”) of my initial solution, before allowing for the precipitation of calcite and sepiolite. I repeat these steps 100s of times by automating the input with Python: 

For step in range(total_steps): 

              evap = str(evaporation_moles[step])      #moles of water to be removed 

              mix = str(mix_volume[step])                     #volume of initial solution to be mixed in 

              leak = str(leak_volume[step])                   #volume of solution to be removed (currently zero) 

              co2 = str(atmospheric_co2[step])            #saturation index of atmospheric co2 

              o2 = str(atmospheric_o2[step])               #saturation index of atmospheric o2

              temp = str(temperature[step])        #temperature of solution

              line =     ('\n #' + str(int(step)+1 

                             '\n USE solution 1; REACTION; H2O -1; ' + evap + ' moles; SAVE solution 1;    #Evaporating previous solution'

                             '\n MIX; 1 1; ' + mix + 'SAVE solution 1;                              #Adding input solution(s)'

                             '\n EQUILIBRIUM_PHASES 1; CO2(g) ' + co2 + ' 1000;'

                             '\n\t\t\t\t O2(g) ' + o2 + ' 1000; SAVE solution 1; END                   #Equilibrating with atmospheric CO2'

                             '\n SOLUTION_MODIFY 1; -temp ' + temp + '; SAVE solution 1;            #Equilibrating to atmospheric temperature'

                             '\n MIX; 1 ' + leak + '; SAVE solution 1;                               #Loss of water by leakage/spillover'

                             '\n USE solution 1; USE kinetics 1; SAVE solution 1; END                 #Calculating mineral precipitation' 

                             ) 

              phreeqc_input.write(line + '\n')

I had initially expected the mass of water to remain at 1kg throughout the simulation, but quickly discovered that sepiolite precipitation was removing water from the solution, such that the % of water removed at each step steadily increases. 

I realise that many people, on PhreeqcUsers and indeed in the few posts on this PhreeqPy forum, for a way of maintaining a constant 1kg of water. It appears that the only/best answer is to calculate the exact “MIX” volume required outside of Phreeqc after each step, which brings me to PhreeqPy. 

I have read the PhreeqcPy examples and located the relevant lines for running the attached script (as the string “phreeqc_input”) in one go: 

              import phreeqpy.iphreeqc.phreeqc_dll as phreeqc_mod

              phreeqc = phreeqc_mod.IPhreeqc() 

              phreeqc.set_selected_output_file_on() 

              phreeqc.load_database(database_path) 

              phreeqc.run_string(phreeqc_input)

This produces the same selected output as in PhreeqC, and is surprisingly much faster (I suspect because the output file is not being written?).

Partly because I usually write my code in order without defining any functions, I am not sure which parts of the examples I should use to stop and start the model? To run just the above for loop part of my input one iteration at a time? Getting the current “mass_H2O” from IPhreeqc to calculate a new value for my “mix_volume” between each iteration.

I notice that including the line "phreeqc.set_selected_output_file_on()" in the "initialize" function of the PhreeqPy example produces only a file for the final shift. Presumably appending this to an array/dataframe in Python would produce the full selected output file as it would appear if run directly in PhreeqC? Pehaps I should then include "set_selected_output_file_on()", the code for reading the "mass_H2O" from it, and "run_string()" within my for loop? Would IPhreeqc remember all of the solutions, kinetics etc. until I start a new instance of IPhreeqc?

N.B.: Your 2011 paper (Mueller_etal_MODFLOWandMORE2011_Proceedings.pdf (phreeqpy.com)) directs you to phreeqpy.com for the code? Is it the same code as in the phreeqpy.com example, or is the code no longer available?

Thanks all,

Tom