Varying residence time using mix adn reactive transport modelling

[FlouretAlex]

Hello PREEQC user,
First of all, I wish to all of you a very happy new year full of PHREEQC adventure =).

I would like to know if it's possible to vary residence time during one experiment by using the MIX and TRANSPORT keywords.  Indeed, till now for my script, I'm using the same flow rate during my experiment for simplification. But now as everything is working great I would like to vary my transport parameter during the experiment but only that, thus keeping all the other data simulated till the change. Maybe there is an example close to what I want in the manual but which one?
Thank you for your help.

Here below you can find my script for better visualization.

Simulation has been done with PHREEPLOT for ease and safety with data manipulation. If you want a script that works with PHREEQC just let me know.

Code: [Select]

SPECIATION
  calculationType "custom"         
  calculationMethod -1                 
  database                C:\Users\International\Documents\PKDLM_BRGM_database_phreeqc_ThermoddemV1.10_06Jun2017.dat
 loopfile "1_RHIZO_H_millet.dat"     
debug 1


CHEMISTRY
SELECTED_OUTPUT
-reset false
-file C:\Users\International\Documents\3_MILLET_RHIZOplant-H_T<vs_i>.sel
-high_precision TRUE
USER_PUNCH
-headings q(nbre_volumes_de_cellulle) t   logCsaq   logKd   Csini   Csaq Kaq Naaq Caaq Mgaq   Kd   Kdsurf   Kdexch   Cssorb   Cssurf   Csexch   Csverif   illiKd   montKd   kaoKd Cs_plante Mg_plante Ca_plante Plt Plante_m
-start
10 Msoil = <m_s>/1000 # soil_masse(Kg)
20 Vwater = <v>/1000 # volume batch (L)
30 MVw = Msoil/Vwater #Kg/L
40 VwM = Vwater/Msoil #L/Kg
50 Csini = sys("Cs") # (mol/L)
60 Caaq = tot("Ca") #(mol/L)
70 Csaq = tot("Cs") #(mol/L)
80 Kaq = tot("K") #(mol/L)
90 Naaq = tot("Na") #(mol/L)
100 Mgaq = tot("Mg") #(mol/L)
110 surfCs = (surf("Cs","Illitex")+ surf("Cs","Montx")+ surf("Cs","Kaolinitex")) #(mol/L)
120 exchCs = (mol("CsIlli")+ mol("CsMont")+ mol ("CsKao")) #(mol/L)
130 illCs = (surf("Cs","Illitex")+ mol("CsIlli")) #(mol/L) #(mol/L)
140 montCs = (surf("Cs","Montx")+ mol("CsMont")) #(mol/L)
150 kaoCs = (surf("Cs","Kaolinitex")+ mol ("CsKao")) #(mol/L)
160 Csfix = surfCs + exchCs #(mol/L)
170 Kd = (Csfix/Csaq)*VwM #(L/Kg)
180 Kdsurf = (surfCs/Csaq)*VwM #(L/Kg)
190 Kdexch = (exchCs/Csaq)*VwM #(L/Kg)
200 Kdilli = (illCs/Csaq)*VwM #(L/Kg)
210 Kdmont = (montCs/Csaq)*VwM #(L/Kg)
220 Kdkao = (kaoCs/Csaq)*VwM #(L/Kg)
230 Cssorb = Csaq*Kd #(mol/Kg)
240 Cssurf = Csaq*Kdsurf #(mol/Kg)
250 Csexch = Csaq*Kdexch #(mol/Kg)
260 Csverif = Csaq + Cssorb*MVw #(mol/L)
270 dt= <t_steps>
280 T_R= <tr> # = cel renew time (s)
290 nb_cel= dt*(step_no)/ T_R #nb renewed cel
300 T_exp = dt*(step_no) #experiment time (s)
310 Cs_plante = KIN("Upcs") #mol
320 Mg_plante = KIN("Upmg") #mol
330 Ca_plante = KIN("Upca") #mol


400 punch nb_cel, (T_exp/86400)-3, log10(Csaq), log10(Kd), Csini, Csaq,Kaq,Naaq,Caaq,Mgaq, Kd, Kdsurf, Kdexch, Cssorb, Cssurf, Csexch, Csverif, Kdilli, Kdmont, Kdkao, Cs_plante, Mg_plante, Ca_plante, tot("Plt"), GET(1)
-end

SOLUTION_MASTER_SPECIES
Plt              Plt-    0.0 Plt 1.0

SOLUTION_SPECIES
Plt- = Plt-
        log_k   0.0

SURFACE 10
-equilibrate 10
-sites_units density 
Illitex <d_i> 97 <mv> # 0.5%
Montx <d_m> 800 <mv> # 1.1%
Kaolinitex <d_k> 10 <mv> # 0.7%
# ^site/nm^2 ^SSA (m^2/g) ^M/V (g/L)
EXCHANGE  10
-equilibrate 10
Illi <e_i>
Mont <e_m>
Kao <e_k>

# Equilibrage spike
SOLUTION 10
 units mol/L
pH <h>
Na <na>
K <k>
Ca <ca>
Mg <mg>
N(5) <n5>
Cl <cl>
S(6) <s6>
C(4) <c4>
EQUILIBRIUM_PHASES   10
Fix_H+ -<h> HCl 10
save surface 3
save exchange 3
END

#Contaminated soil ==> RIZO
SOLUTION 3 # ==> stagnant cell
units mol/L
pH <h>
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04
Cs <cs>
Plt 10
use surface 3 #number of surface and exchange need to be the same as the cell here solution 3 ==> cell 3
use Exchange 3
EQUILIBRIUM_PHASES   3
Fix_H+ -<h> HCl  10
END

SOLUTION 0-1 # inlet solution 0 # define solution 1 since PHREEQC wants at least 1 cell in the colum for transport ==> mobile cell.
 units mol/L
pH <h>
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04


Solution 100 #for evaporation/evapotranspiration

Rates
Upcs
-start
10 Z = TOTAL_TIME # s
20 if Z <= 259200 then goto 240 # soil/plant contact time => day 3 after soil/solution equilibration
100 km = 1E-7 # mol
110 Vmax = 1E-9 # mol/g-1.s-1
120 bm = GET(1) # biomass (g)

200 Vp = ((Vmax*mol("Cs+"))/(km + mol("Cs+"))) # Cs absorption rate (mol/g-1.s-1)
210 rate = Vp * bm # Cs absorption flux (mol/.s-1)
220 moles = -rate * TIME # Cs absorbed (mol)
230 if (m - moles) < 0 then moles = m
240 SAVE moles
-end

Upk
-start
10 Z = TOTAL_TIME # s
20 if Z <= 259200 then goto 230 # soil/plant contact time => day 3 after soil/solution equilibration
100 km = 9E-7 # mol
110 Vmax = 1E-10 # mol/g-1.s-1
120 bm = GET(1) # biomass (g)

200 Vp = ((Vmax*mol("K+"))/(km + mol("K+"))) # K absorption rate (mol/g-1.s-1))
210 rate = Vp * bm # K absorption flux (mol/.s-1)
220 moles = -rate * TIME # K absorbed (mol)
230 if (m - moles) < 0 then moles = m
240 SAVE moles
-end

Upna
-start
10 Z = TOTAL_TIME # s
20 if Z <= 259200 then goto 230 # soil/plant contact time => day 3 after soil/solution equilibration
100 km = 1E-2 # mol
110 Vmax = 1E-15 # mol/g-1.s-1
120 bm = GET(1) # biomass (g)

200 Vp = ((Vmax*mol("Na+"))/(km + mol("Na+"))) # Ca absorption rate (mol/g-1.s-1)
210 rate = Vp * bm # Ca absorption flux (mol/.s-1)
220 moles = -rate * TIME # Ca absorbed (mol)
230 if (m - moles) < 0 then moles = m
240 SAVE moles
-end

Upca
-start
10 Z = TOTAL_TIME # s
20 if Z <= 259200 then goto 230 # soil/plant contact time => day 3 after soil/solution equilibration
100 km = 1E-4 # mol
110 Vmax = 1E-11 # mol/g-1.s-1
120 bm = GET(1) # biomass (g)

200 Vp = ((Vmax*mol("Ca+2"))/(km + mol("Ca+2"))) # Ca absorption rate (mol/g-1.s-1)
210 rate = Vp * bm # Ca absorption flux (mol/.s-1)
220 moles = -rate * TIME # Ca absorbed (mol)
230 if (m - moles) < 0 then moles = m
240 SAVE moles
-end

Upmg
-start
10 Z = TOTAL_TIME # s
20 if Z <= 259200 then goto 240 # soil/plant contact time => day 3 after soil/solution equilibration
100 km = 5E-6 # mol
110 Vmax = 1E-8 # mol/g-1.s-1
120 bm = GET(1) # biomass (g)

200 Vp = ((Vmax*mol("Mg+2"))/(km + mol("Mg+2"))) # Mg absorption rate (mol/g-1.s-1)
210 rate = Vp * bm # Mg absorption flux (mol/.s-1)
220 moles = -rate * TIME # Mg absorbed (mol)
230 if (m - moles) < 0 then moles = m
240 SAVE moles
-end

Biomass
-start
10 Z = TOTAL_TIME
20 IF Z <= 259200 THEN GOTO 110
30 t = TOTAL_TIME- 259200 #s
40 ib = 1.35 #initial biomass (g)
50 r = 7.77514921513017E-07 # max grow rate (g/s)
60 k = 5.00337635961718 # max biomass (g)
70 y = k-EXP(-r*t)*(k-ib) # biomass (g) at time t
80 rate = y/t
90 moles = -rate * TIME
100 PUT(y,1)
110 SAVE moles
-end

KINETICS 3
Upcs
-formula PltCs
-m     0
-tol   1.e-16
Upk
-formula PltK
-m     0
-tol   1.e-16
Upna
-formula PltNa
-m     0
-tol   1.e-16
Upca
-formula Plt2Ca
-m     0
-tol   1.e-16
Upmg
-formula Plt2Mg
-m     0
-tol   1.e-16
Biomass
-formula  H  0
-m 0
-tol 1e-8

END


MIX 3 # simulate water evaporation/evapotranspiration
0 0.2
100 -0.2
3 1
END


TRANSPORT
-cells 1
-lengths 0.01 # 1 cm soil
-shifts <sh> # number of time_steps, total time is 24 days
-time_step <t_steps> # each time_step (s) the MIX is performed
-flow_direction diffusion_only
-boundary_conditions constant closed
-dispersivities  1
-diffusion_coefficient 2.3E-08
-stagnant 1
-punch_cells 3
-punch_frequency 1
END



[dlparkhurst]

Residence time depends only on the time step. I think if you want to change the residence time, you will have to change the time step. That means that you would need another TRANSPORT definition each time the residence time changes. The good thing is that you only need to define TRANSPORT; -time_step <t_steps>; END, and all of the other parameters remain unchanged. The bad thing is that requiring new definitions at each change may make it difficult to arrange in PhreePlot.

An alternative is to use IPhreeqc with a programming or scripting language. That would give you control over your variables, while allowing relatively easy manipulation of multiple TRANSPORT simulations.

[FlouretAlex]

Hi and thank you for your response.
 So if I understand correctly I just need to have a succession of transport one after another? Then his time will be incremented from the previous transport? As I have kinetics that is dependent on the total time of the simulation (biomass) I will need to have this to be continuous and not set to zero each time a new TRANSPORT data block is used. I'll try and see how it goes.

I would love to work with IPhreeqc coupled with R, but I didn't take much time to see how it works and I do not have much time for this at the moment. I will have a look after the end of my Ph.D. and I can't wait for it as It looks so much fun!

Thank you for your help!

[dlparkhurst]

I think you also will have to define -initial_time (set to the cumulative time to that point) at each TRANSPORT definition to get the value of the Basic function TOTAL_TIME that you want.

[FlouretAlex]

OK, I see.
Thank you for your help.