Transport - decaying concentration in infilling solution

[AntonioG]

#Hello,

# I'm new with phreeqc, hopefully someone could provide me bit of help on that one. I'm trying to setup a simple transport model of lithium based on a first order decay of concentrations in the infilling solution. The infilling solution attempts to simulate the release of lithium in groundwater from a known initial source ("Li_source" =19 mg/L) and a first order decay. To do so I created a rate "Li_release" and the associated kinetic block. It seems that Li concentration does not decrease with time in any of the cells.  I'm not sure whether this is the right approach and/or where the kinetic block must be defined. Does it have to be linked to the infilling solution? the first transport cell? or somewhere else in the Transport data block? Thanks

DATABASE LLNL.DAT

RATES
Li_release
  -start
   10  m_source = tot("Li_source")
   15  if (m_source <= 1e-20) then goto 200
   30  rate = m_source*parm(1)
   40  moles = rate * time
   50 if (moles > m) then moles = m
   200 SAVE moles
  -end
END

SOLUTION_MASTER_SPECIES
Li_source   Li_source      0.0   Li_source  7   7

SOLUTION_SPECIES
Li_source  = Li_source   
        log_k           0.0
   
SOLUTION 0 #Source composition
   -units   mg/L   
   Li_source   19
   #Li      19
   Li      0.0019
   -water   1.0
KINETICS 0
Li_release
   -parms 5.991E-2
   -formula Li_source -1 Li 1   
COPY solution 0 100
COPY kinetics 0 100
SAVE SOLUTION 0
END

# KINETICS 1
# Li_release
# -parms 5.991E-2
# -formula Li_source -1 Li 1.0
# COPY kinetics 1 100
# END

#Groundwater
SOLUTION 1-36
        -units     mg/L
   Li   0.019
   Li_source 0.00001
        -water    1

COPY solution 1 101
END

SOLUTION 37 #Source residual release
   units   mg/L   
   Li   0.95
   -water   1.0
COPY solution 37 111
END

#############################################
SELECTED_OUTPUT 1
        -file   output.xls
   -reset               true
   -pH                   false
   -pe                   false
   -temperature         false
   -alkalinity            false
   -ionic_strength      false
   -water                false
   -percent_error       false

USER_PUNCH 1
   -headings    Cum_Time_months Li
   -start
   11 PUNCH TOTAL_TIME/(86400*365.5)
   22 PUNCH TOT("Li")*7*1000

USER_GRAPH 2
        -headings Li Li_source
        -axis_titles "Years" "Concentration [mg/L]"
        -axis_scale x_axis auto #0 300
        -axis_scale y_axis auto #0 300
   -axis_scale secondary_y_axis auto
   -plot_concentration_vs t
   -connect_simulations    true
        -start
          10 x = TOTAL_TIME/(86400*365.5)
     20 if (step_no < 0) then goto 100 
     30 PLOT_XY x,TOT("Li")*7*1000, color = Blue, line_width = 1, symbol_size = 3
     #40 PLOT_XY x,TOT("Li_source")*7*1000, color = Green, line_width = 1, symbol_size = 3
     100 REM end
          -end
    
#############################################
TRANSPORT 1 TSF Seepage for first
    -cells                    36
    -shifts                   224
    -time_step                817 day
    -initial_time       0
    -lengths                  36*27.3
    -flow_direction      forward
    -porosities          36*0.0075
    -dispersivities          36*50
    -correct_disp           true
    -diffusion_coefficient    1.25e-9
    -boundary_conditions     flux flux
    -punch_cells          1# breakthrough in last cell only   
    -punch_frequency         1
   
COPY solution 111 0
#COPY kinetics 100 0 #
END

# TRANSPORT 1 Last XX days of Groundwater flow
        # -shifts               96
   # -time_step            60 day
# END

[dlparkhurst]

Here is my approach to what you have described. I set up the kinetics in cell 1 to create an exponentially decaying Li concentration. You can play with the constant "k" to suit your time sequence, but the script calculates a target concentration as a function of time. It then uses kinetics to approach that target concentration. The constant k2 determines how quickly it reaches the target concentration, but as long as k2 >> k, the actual concentration will be close to the target concentration.

If you just take out Li, unanticipated redox reactions will occur because you are removing Li metal. I've added Cl to the system; taking out LiCl does not cause any redox reactions.

Code: [Select]

RATES
Li_Function
  -start
   10 k = parm(1)/3.15e7 # 1/s
   20 target_Li = (0.019/GFW("Li")) * exp(-k*TOTAL_TIME)
   30 k2 = .00001       # 1/s but fast relative to k
   40 rate = (TOT("Li") - target_Li)*k2
   50 moles = rate * TIME
  200 SAVE moles
  -end
END
SOLUTION 0-1 #Source composition
-units   mg/L   
Li      19
Cl      1 charge
   -water   1.0

KINETICS 1
-cvode
Li_function
   -parms 5.991E-2 # 1/y ??
   -formula LiCl -1   
END

#Groundwater
SOLUTION 2-36
        -units     mg/L
   Li   0.019
   Cl   1 charge
  -water    1

END

TRANSPORT
    -cells                 36
    -shifts                72 #224
    -time_step             70588800 # seconds
    -lengths               27.3
    -dispersivities        50
    -correct_disp          true
    -diffusion_coefficient 1.25e-09
    -thermal_diffusion     2   1.25e-09
    -punch_cells           1 36
 
USER_GRAPH 1
    -headings               Li
    -axis_titles            "Years" "Concentration [mg/L]" ""
    -chart_title            "Cell 1"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 10 x = TOTAL_TIME/(86400*365.5)
 20 if (cell_no <> 1) then goto 100
 30 PLOT_XY x,TOT("Li")*7*1000, color = Blue, line_width = 1, symbol_size = 3
100 REM end
  -end 

USER_GRAPH 36
    -headings               Li
    -axis_titles            "Years" "Concentration [mg/L]" ""
    -chart_title            "Cell 36"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 10 x = TOTAL_TIME/(86400*365.5)
 20 if (cell_no <> 36) then goto 100
 30 PLOT_XY x,TOT("Li")*7*1000, color = Blue, line_width = 1, symbol_size = 3
100 REM end
  -end
END



[AntonioG]

Hi dlparkhurst,
Thanks very much for your reply, it will help me to move forward.

[AntonioG]

Hi dlparkhurst,
Thanks for your helpful reply, I would have a small additional question.
The decay of my source term is very quick (~650 days) vs the total duration of the simulation (~200 years).
Is it possible to "turn off" the kinetic block when concentrations become negligible and then increase the timestep and cell lengths to speed up the simulation? Thanks.

DATABASE LLNL.DAT

RATES
Li_Function
  -start
   10 k = parm(1)/3.15e7 # 1/s
   20 target_Li = (0.020/GFW("Li")) * exp(-k*TOTAL_TIME)
   30 k2 = 0.0001       # 1/s but fast relative to k
   40 rate = (TOT("Li") - target_Li)*k2
   50 moles = rate * TIME
  200 SAVE moles
  -end
END
SOLUTION 0-1 #Source composition
-units   mg/L   
Li      20
Cl      1 charge
   -water   1.0

KINETICS 1
-cvode
-bad_step_max 1000
Li_function
   -parms 0.593 #5.991E-2 # 1/y ??
   -formula LiCl -1   
END


#Groundwater
SOLUTION 2-950
        -units     mg/L
   Li   0.002
   Cl     1 charge
  -water    1

END

SOLUTION 951 #Background infilling concentration in cell 0
   -units   mg/L   
   Li   0.002
   -water   1.0
COPY solution 951 952
END


TRANSPORT 1 #initial transport block to capture source decay
    -cells                 950
    -shifts                10
    -time_step             37 days #28166400 # seconds
    -lengths               1
    -dispersivities       4
    -porosities       0.0075
    -correct_disp          true
    -diffusion_coefficient 1.25e-09
    -thermal_diffusion     2   1.25e-09
    -punch_cells           1 950
END

TRANSPORT 1
        -shifts               392
   -cells              95
   -lengths           6
   -time_step            373 day
    COPY solution 952 0 
END

USER_GRAPH 1
    -headings               Li
    -axis_titles            "Years" "Concentration [mg/L]" ""
    -chart_title            "Cell 1"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 10 x = TOTAL_TIME/(86400*365.5)
 20 if (cell_no <> 1) then goto 100
 30 PLOT_XY x,TOT("Li")*7*1000, color = Blue, line_width = 1, symbol_size = 3
100 REM end
  -end

USER_GRAPH 92
    -headings               Li
    -axis_titles            "Years" "Concentration [mg/L]" ""
    -chart_title            "Cell 950"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 10 x = TOTAL_TIME/(86400*365.5)
 20 if (cell_no <> 950) then goto 100
 30 PLOT_XY x,TOT("Li")*7*1000, color = Blue, line_width = 1, symbol_size = 3
100 REM end
  -end
END

[dlparkhurst]

I'll just say that I would never use 1000 cells in a simulation with PHREEQC.

If the decay is fast and you are looking at long times, you already now the answer.

I guess you could resample your cells, selecting 1, 11, 21, etc and changing transport to 95 cells, each with 10 times the length. (COPY  solution 11 2; COPY 21 3; ... and redefine TRANPORT). Seems like more work than you want to do. I'd just use the coarser grid from the start.

You can set some limits in RATES with a factor like c/(1e-10+c) on the rate, so the concentration changes slowly below a concentration of 1e-10, for example.