HPx

#TEMPORARILY -> should be produced by GUI
GLOBAL_VARIABLE _TIMEUNIT  5    "Year"
GLOBAL_VARIABLE _SPACEUNIT 3    "Meter"
#INPUT PARAMETERS
GLOBAL_VARIABLE InertSolidPhase    0.35    ""
DATABLOCK
MineralProp
-origin 1
-headings Name MolarVolume density GFW InitialVolume rate
-rowid Name
-start
calcite 36.93336 2.71 100.09 30 5E-8
-end
#WRITE KINETIC BLOCK TO KINETICS.IN 
#using datablock MineralProp
CALCULATE_VALUES
definition
-state definition
-mybasic
-start
 tab$ = "RatesAq"
 mult$ = " * " 
 qt$ = chr(34)
 empty$ = "       "
 ds$="MineralProp"
 nrrec=get_defnrrec(ds$)
 def WriteKinetics()
   strkin = "KINETICS 1001" + EOL + "-material 1" + eol
   write("MINERAL REACTION NETWORK - Overview of kinetic parameters" + eol)
   for i=1 to nrrec
     mineral$=get_deflabel(ds$,i)
     volume=get_def(ds$,"InitialVolume",i)
     molarvolume=get_def(ds$,"MolarVolume",i)
     amount=volume*10/molarvolume
     rate=get_def(ds$,"rate",i)/CONV_T("s")
     strkin = strkin + mineral$ + eol
     strkin = strkin +  " -m0 " + str(amount) + eol
     strkin = strkin +  " -parms " + str(rate) + eol
     line$= mineral$ + "  Amount (mol/dm�) " + str(amount) + "   Rate(1/y)  " + str(rate)
     write(line$ + eol)
   next i
   strkin = strkin + "-cvode true" + EOL
   strkin = strkin + EOL+EOL
   openf("kinetics.in")
   writef(strkin)
   closef()
 enddef
 WriteKinetics()
-end
#DEFINE RATE EQUATIONS
RATES
calcite
-mybasic
-start
 REM parm(1) k_0 mol/dm�/s
 MolarVolume = Get_Def("MineralProp","MolarVolume","calcite")
 phi0 = m0()*MolarVolume # initial bulk volume
 phi = m()*MolarVolume
 k_red = (phi/phi0)^(2/3)
 k = k_red*parm(1)
 rate = k*(1 - SR("Calcite"))
 mole = rate*time()
 if (mole > m()) then mole = m()
 save(mole)
-end 
#UPDATE TRANSPORTPROPERTIES
Calculate_values
UpdateScaleFactors
-mybasic
-start
 solphase = get_var("InertSolidPhase") + sum_solid("kr","MineralProp","MolarVolume")/1000 
 porosity=1-solphase
 change_por(porosity)
 update_sfwater(porosity)
 if (total_time() < 1E-15) then
  a = ((1-porosity)^2 / porosity^3)
  put_def("IniFactor",a,cell_no())
 else
  scale = (porosity^3 / (1-porosity)^2) * get_def("IniFactor",cell_no())
  change_sfcond(scale)
 endif
 save(0)
-end
initial_por
-state initial
-mybasic
-start
 a=calc_value("UpdateScaleFactors")
 Change_por0(get_por())
 change_sfwater0(get_sfwater()) 
-end
cal_update
-state update
-mybasic
-start
 a = calc_value("UpdateScaleFactors")
-end

solution 1001 initial condition mobile phase
-temp 25
-units mol/kgw
-pH 9.375
Ca 1.566E-4
C 2.566E-4
S 1E-10
end
solution 3001 boundary condition
-temp 25
-units mol/kgw
-pH 3
Ca 1E-4
C 1E-2
S 6.457E-4

INCLUDE$ kinetics.in
reactive_transport
-update hp t
-openmp t

rt_output
 -profile gnuplot
user_punch
-headings Distance por HydrCond cal press waterflux
-plot_titles "Distance" "Porosity" "Hydraulic Conductivity" "Calcite vol fraction" "Hydraulic head" "Flux"
-plot_units "m" "-" "m s^{-1}" "m^3 m^{-3}" "m" "m^3 d^{-1}"
-mybasic
-start
  punch(-dist())
  punch(get_por())
  punch(transprop(5) * get_sfcond() * conv_t("s"))
  punch(kin("calcite")*Get_Def("MineralProp","MolarVolume","calcite") / 1000)
  punch(pressure_head() / conv_x("m"))
  punch(-water_flux() * conv_t("d"))
-end
global_variable MIN3PProfiles$ "MIN3P - Profiles (figure 2).txt"
global_variable HPxProfiles$ "nod_inf_chem.out"
global_variable MIN3PTime$ "MIN3P - Time (figure 2).txt"
global_variable HPxTime$ "obs_nod_chem81.out"
selected_graph 
#generate figure 2 from Xie et al. (2015) using HPx results and MIN3P
-start
 set linestyle 1 lc "red" lw 3 
 set linestyle 2 lc "blue" lw 3 
 set linestyle 3 lc "green" lw 3 
 set linestyle 4 lc "magenta" lw 3 
 set linestyle 51 lc "black" lw 2 pi 5 dt 3
 set linestyle 52 lc "black" lw 2 pi 50 dt 3
 tCreated = 'Created with @_HPVERSION$ on @_DATE$ - @_TIME$ '
 set label 1 tCreated at graph 1, screen 0.98 right font ",6" textcolor "gray40"
 tProj = '@_PROJECTNAME$ '
 set label 2 tProj at graph 1, screen 0.02 right font ",6" textcolor "gray40"
-end
-fig_nr 1
-fig_group "PPT-SS-DIS(B1)"
-fig_variabletype "Xie et al - 2015 - Fig 2"
-fig_name "Calcite"
-fig_x Distance
-fig_y VolCal10 <ds "@MIN3PProfiles$ " ls 1 k "MIN3P - 10a" > 
-fig_y VolCal100 <ds "@MIN3PProfiles$ " ls 2 k "MIN3P - 100a" > 
-fig_y VolCal120 <ds "@MIN3PProfiles$ " ls 3 k "MIN3P - 120a" > 
-fig_y cal <ds "@HPxProfiles$ " ls 51 i 1 k "HPx" style linespoints>
-fig_y cal <ds "@HPxProfiles$ " ls 51 i 2 vk false style linespoints>
-fig_y cal <ds "@HPxProfiles$ " ls 51 i 3 vk false style linespoints>
-fig_xt "Distance [m]"
-fig_yt "Volume fraction of caclite [m^3 m^{-3}]"
-fig_start
  set key at graph 0.45 , graph 0.25
  set yrange [0:0.4]
-fig_end
-fig_nr 2
-fig_group "PPT-SS-DIS(B1)"
-fig_variabletype "Xie et al - 2015 - Fig 2"
-fig_name "Porosity"
-fig_x Distance
-fig_y Por10 <ds "@MIN3PProfiles$ " ls 1 k "MIN3P - 10a" > 
-fig_y Por100 <ds "@MIN3PProfiles$ " ls 2 k "MIN3P - 100a" > 
-fig_y Por120 <ds "@MIN3PProfiles$ " ls 3 k "MIN3P - 120a" > 
-fig_y por <ds "@HPxProfiles$ " ls 51 i 1 k "HPx" style linespoints>
-fig_y por <ds "@HPxProfiles$ " ls 51 i 2 vk false style linespoints>
-fig_y por <ds "@HPxProfiles$ " ls 51 i 3 vk false style linespoints>
-fig_xt "Distance [m]"
-fig_yt "Porosity [-]"
-fig_start
  set key at graph 0.35, graph 0.75
  set yrange [0.3:0.7]
-fig_end
-fig_nr 3
-fig_group "PPT-SS-DIS(B1)"
-fig_variabletype "Xie et al - 2015 - Fig 2"
-fig_name "Hydraulic Head"
-fig_x Distance
-fig_y head10 <ds "@MIN3PProfiles$ " ls 1 k "MIN3P - 10a" > 
-fig_y head100 <ds "@MIN3PProfiles$ " ls 2 k "MIN3P - 100a" >
-fig_y head120 <ds "@MIN3PProfiles$ " ls 3 k "MIN3P - 120a" >
-fig_y press <ds "@HPxProfiles$ " ls 51 i 1 k "HPx" style linespoints>
-fig_y press <ds "@HPxProfiles$ " ls 51 i 2 vk false style linespoints>
-fig_y press <ds "@HPxProfiles$ " ls 51 i 3 vk false style linespoints>
-fig_xt "Distance [m]"
-fig_yt "Hydraulic Head [m]"
-fig_start
  set key inside
  set yrange [0:0.008]
-fig_end
-fig_nr 4
-fig_group "PPT-SS-DIS(B1)"
-fig_variabletype "Xie et al - 2015 - Fig 2"
-fig_name "Water Flux"
-fig_x time
-fig_y Outflow <ds "@MIN3PTime$ " ls 1 k "MIN3P" > 
-fig_y waterflux <ds "@HPxTime$ " ls 52 k "HPx" style linespoints>
-fig_xt "Time [years]"
-fig_yt "Flux of outflow [m^3 d^{-1}]"
-fig_start
  set key inside top left
  set xrange [0:150]
-fig_end