Possible to add a legend to define time periods on porosity-distance plots?

[rworden]

I am trying to refine a USER_GRAPH of distance (x-axis) versus porosity (y-axis) during flow of CO2 through a porous limestone.  I have set the simulation to run for one year and I would like to show how the porosity varies  with distance and time.  I could take the .txt file and plot in R or excel but it would be easier and faster if I could produce working USER_GRAPHs, real-time in phreeqc, with a legend defining each time interval.  The USER_GRAPH that results from the code knows that each time interval is different as the code assigns a different line style to each time interval, so it seems not-impossible that each time interval could be referenced in a legend to the USER_GRAPH.  It would also be useful if the USER_GRAPH could show time step zero as a (flat) line instead of a single point.  Any advice or comments, especially on the first point, gratefully received.

Code: [Select]

TITLE Modeling Dissolution of Calcite with 300 atm (30 MPa) CO2 during advection through 100m of rock over 1 year

SOLUTION 0
    Na 1.0
    Ca 0.0
    Cl 1.0
    pH 7
    temp 25
    pressure 200
    pe 4
    redox pe
    units mmol/kgw
    density 1
    -water 1

GAS_PHASE 0
   -fixed_pressure       
   -volume 0.0375         
   -pressure 300         
   -temperature 25
   CO2(g)  300           

END

SOLUTION 1-10 Initial solution in column after CO2 equilibration   
    Na 1.0
    Ca 0.0
    Cl 1.0
    pH 7
    temp 25
    pressure 200
    pe 4
    redox pe
    units mmol/kgw
    density 1
    -water 1

GAS_PHASE 1-10
    -fixed_volume           
    -volume 0.01             
    -temperature 55
    CO2(g)  0.000001
    H2O(g)  0.01

RATES 1-10
Calcite
-start
  1   rem M = current number of moles of calcite
  2   rem M0 = number of moles of calcite initially present
  3   rem PARM(1) = A/V, cm^2/L
  4   rem PARM(2) = exponent for M/M0
  10  si_cc = SI("Calcite")
  20  if (M <= 0 and si_cc < 0) then goto 200
  30    k1 = 10^(0.198 - 444.0 / TK )
  40    k2 = 10^(2.84 - 2177.0 / TK)
  50    if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
  60    if TC > 25 then k3  = 10^(-1.1 - 1737.0 / TK )
  70    t = 1
  80    if M0 > 0 then t = M/M0
  90    if t = 0 then t = 1
  100   area = PARM(1) * (t)^PARM(2)
  110   rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
  120   rem 1e-3 converts mmol to mol
  130   rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
  140   moles = rate * TIME
  200 SAVE moles
-end

KINETICS 1-10
Calcite
    -formula      CaCO3  1
    -m0           81.235 # quantity equivalent to 25% porosity limestone
    -parms        5 0.3
    -tol          1.e-8

#END

TRANSPORT
    -cells 10                               
    -length 10                               
    -dispersivity 0.002                     
    -shifts 876.6                           
    -flow_direction forward
    -time_step 36000                         
    -boundary_conditions flux flux           
    -diffusion_coefficient 0
    -porosities 10*0.25
    -punch_cells 1-10
    -punch_frequency 87.6
    -print_cells 1-10

USER_PUNCH 1-10
    -headings water_density Pressure PP_CO2 si_cc TOTAL_TIME Cell Porosity
    -start
10 PUNCH RHO               
20 PUNCH PRESSURE           
30 PUNCH PR_P("CO2(g)")     
40 PUNCH SI("Calcite")     
45 PUNCH TOTAL_TIME
50 PUNCH CELL_NO
70 new_por = 1000 + 36.93 * (81.235 - KIN("Calcite"))
80 Porosity =  new_por / (new_por + 36.93 * KIN("Calcite"))
110 PUNCH Porosity   
   -end

SELECTED_OUTPUT 1
    -file Calcite_water_kinetics.txt
    -reset false
    -step true
    -pH true
    -reaction true
    -totals Ca C(4)
    -equilibrium_phases
    -saturation_indices Calcite
    -kinetic_reactants Calcite
    -distance true

USER_GRAPH 1
    -headings               
    -axis_titles            "Distance, metres" "Fractional Porosity" ""
    -axis_scale x_axis      auto auto auto auto   
    -chart_title            "CO2-rich water reacts with limestone - advection - with kinetics"
    -axis_scale y_axis      auto auto auto auto
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 new_por = 1000 + 36.93 * (81.235 - KIN("Calcite"))
20 graph_x (cell_no * 10) - 5 # "-5" to put the value in the middle of the cell
30 graph_y new_por / (new_por + 36.93 * KIN("Calcite"))
  -end
    -active                 true

END

[dlparkhurst]

I don't have a great answer for you. The following script is a little clunky, but it does what you wanted, I think.

There are other options, of using PUT and GET to save the data and ultimately plot it, or USER_PUNCH to write the code that does the TRANSPORT and USER_GRAPH data blocks, but those get even clunkier.

Code: [Select]

TITLE Modeling Dissolution of Calcite with 300 atm (30 MPa) CO2 during advection through 100m of rock over 1 year

SOLUTION 0
    Na 1.0
    Ca 0.0
    Cl 1.0
    pH 7
    temp 25
    pressure 200
    pe 4
    redox pe
    units mmol/kgw
    density 1
    -water 1

GAS_PHASE 0
   -fixed_pressure       
   -volume 0.0375         
   -pressure 300         
   -temperature 25
   CO2(g)  300           

END

SOLUTION 1-10 Initial solution in column after CO2 equilibration   
    Na 1.0
    Ca 0.0
    Cl 1.0
    pH 7
    temp 25
    pressure 200
    pe 4
    redox pe
    units mmol/kgw
    density 1
    -water 1

GAS_PHASE 1-10
    -fixed_volume           
    -volume 0.01             
    -temperature 55
    CO2(g)  0.000001
    H2O(g)  0.01

RATES 1-10
Calcite
-start
  1   rem M = current number of moles of calcite
  2   rem M0 = number of moles of calcite initially present
  3   rem PARM(1) = A/V, cm^2/L
  4   rem PARM(2) = exponent for M/M0
  10  si_cc = SI("Calcite")
  20  if (M <= 0 and si_cc < 0) then goto 200
  30    k1 = 10^(0.198 - 444.0 / TK )
  40    k2 = 10^(2.84 - 2177.0 / TK)
  50    if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
  60    if TC > 25 then k3  = 10^(-1.1 - 1737.0 / TK )
  70    t = 1
  80    if M0 > 0 then t = M/M0
  90    if t = 0 then t = 1
  100   area = PARM(1) * (t)^PARM(2)
  110   rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
  120   rem 1e-3 converts mmol to mol
  130   rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
  140   moles = rate * TIME
  200 SAVE moles
-end

KINETICS 1-10
Calcite
    -formula      CaCO3  1
    -m0           81.235 # quantity equivalent to 25% porosity limestone
    -parms        5 0.3
    -tol          1.e-8

END

TRANSPORT
    -cells 10                               
    -length 10                               
    -dispersivity 0.002                     
    -shifts  1 #876.6                           
    -flow_direction forward
    -time_step 0 #36000                         
    -boundary_conditions flux flux           
    -diffusion_coefficient 0
    -porosities 10*0.25
    -punch_cells 1-10
    -punch_frequency 1
    -print_cells 1-10
    -print_frequency 1

USER_PUNCH 1-10
    -headings water_density Pressure PP_CO2 si_cc TOTAL_TIME Cell Porosity
    -start
10 PUNCH RHO               
20 PUNCH PRESSURE           
30 PUNCH PR_P("CO2(g)")     
40 PUNCH SI("Calcite")     
45 PUNCH TOTAL_TIME
50 PUNCH CELL_NO
70 new_por = 1000 + 36.93 * (81.235 - KIN("Calcite"))
80 Porosity =  new_por / (new_por + 36.93 * KIN("Calcite"))
110 PUNCH Porosity   
   -end

SELECTED_OUTPUT 1
    -file Calcite_water_kinetics.txt
    -reset false
    -step true
    -pH true
    -reaction true
    -totals Ca C(4)
    -equilibrium_phases
    -saturation_indices Calcite
    -kinetic_reactants Calcite
    -distance true

USER_GRAPH 1
    -headings               dist 0_years
    -axis_titles            "Distance, metres" "Fractional Porosity" ""
    -chart_title            "CO2-rich water reacts with limestone - advection - with kinetics"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 new_por = 1000 + 36.93 * (81.235 - KIN("Calcite"))
20 graph_x DIST #(cell_no * 10) - 5   
30 graph_y new_por / (new_por + 36.93 * KIN("Calcite"))
  -end
    -active                 true
END
TRANSPORT
    -time_step 36000 
    -shifts 87
    -punch_frequency 87
    -print_frequency 87
USER_GRAPH 1
    -headings               dist 0.1_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.2_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.3_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.4_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.5_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.6_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.7_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.8_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 0.9_year
END
TRANSPORT
USER_GRAPH 1
    -headings               dist 1.0_year
END


[rworden]

Many thanks, David.  That certainly does work.  I appreciate your time and the effective code.