Unusual peaks in transport plots

[rfembilejr]

Hi,

I was able to model transport of solutes through my columns. However, I am seeing these unusual peaks (see attached) in the plots. My other columns do not have the same peaks so I am really puzzled how and why these peaks appear.
Any thoughts?

Also, I am a bit surprised why SO4 drops immediately. I tried checking the output for any possible complexes that may consume SO4 but did not really found something that is saturated enough to control SO4. I would appreciate any thoughts on this as well. Below is my code if someone wants to take a look.

TITLE RA-04
SOLUTION 0 Pure water
    temp      25
    pH        7.0
EQUILIBRIUM_PHASES 0       # Equilibrates with the atmosphere
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
SAVE SOLUTION 0
END

PHASES
Pyrrhotite                              
      Fe0.903Ni0.011S = +0.903Fe+2   +0.011Ni+2  -1.000H+   + 1.00HS-   -0.172e-
   log_k        0
Forsterite                              
   Mg1.675Fe0.32Ni0.005SiO4 + 4H+ = 2H2O + SiO2 + 1.675Mg+2 + 0.32Fe+2 + 0.005Ni+2 
   log_k     0
Enstatite                              
   Mg1.686Fe0.322SiO3 + 2H+ = 1H2O + SiO2 + 1.686Mg+2 + 0.322Fe+2  + 2.016e-
   log_k   0
END

SOLUTION 1 Rana (column)   # fifth rinse; initial solution in column (fixed pO2)
    temp      20
    pH        4.02      
    pe        9.58         
    redox     pe
    units     mg/l
    density   1
    Alkalinity 0.1
    Ca        8.87             
    Cu        0.1515          
    Fe        0.0121       
    Mg        40.899
    S(6)      228.195       
    Si        10
    Ni        2.222
    -water    0.9
EQUILIBRIUM_PHASES 1
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
    Fe(OH)3 0 10
    #Gypsum 0 0
    #FeSO4 0 0.01

SAVE SOLUTION 1
END

RATES
Pyrrhotite
-start
1 rem PARM(1) = log10(A/V, 1/dm)
2 rem PARM(2) = exp for (m/mo)
3 rem PARM(3) = exp for H+
4 rem PARM(4) = exp for O2
10 if (m <= 0) then goto 200
20 if (SI("Pyrrhotite") >= 0) then goto 200
30 lograte = -5.79 + PARM(1) + PARM(2)*log10(m/m0) + PARM(3)*lm("H+") + PARM(4)*lm("O2")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 SAVE moles
-end
Forsterite
-start
1 rem PARM(1) = log10(A/V, 1/dm)
2 rem PARM(2) = exp for (m/mo)     
3 rem PARM(3) = exp for H+
10 if (m <= 0) then goto 200
20 if (SI("Forsterite") >= 0) then goto 200
30 lograte = -5.96 + PARM(1) + PARM(2)*log10(m/m0) + PARM(3)*lm("H+")    ## k from Inverse modelling
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 SAVE moles
-end
END

KINETICS 1
#PARM(2) is of 2/3 for uniformly dissolving spheres/cubes(Parkhurst and Appelo, 1999)
Pyrrhotite
   -formula Fe0.903Ni0.011S
      -m 0.25                 
      -m0 0.25
      -PARMS -0.71 0.67 0.25 0.1 # 0.1 is the optimum O2 reaction order
   -tol 1e-9               
Forsterite
   -formula Mg1.675Fe0.32Ni0.005SiO4
      -m 10.71                         
      -m0 10.71
      -parms 3.09 0.67 0.28             
   -tol 1e-9
-step_divide 1
-runge_kutta 6
-bad_step_max 1000000000
END

EXCHANGE 1
-equilibrate 1
X 0.1   #tailings impoundments CEC can be 1-10 meq/100g
END

SURFACE 1
       Hfo_wOH Fe(OH)3 equilibrium_phase 0.2 54300 # 0.3 800 18.5 
       Hfo_sOH Fe(OH)3 equilibrium_phase 0.005
      -equilibrate with solution 1
END

REACTION 1
H2O -1.0 # -1.0 is the relative moles of the reaction.
2.8 moles #this is the number of moles to be removed, but  must be done every week; assume only half is gone
END

TRANSPORT
    -cells                 1
    -shifts                65      
    -time_step             604800      
    -lengths               0.12      
    -dispersivities        7.92e-4
    -correct_disp          true
    -diffusion_coefficient 1e-9 #m2/s
    -print_cells           1
    -punch_cells           1   
    -flow_direction        forward
    -boundary_conditions   flux flux

USER_PUNCH 1
    -headings Sulfate_mg/l Mg_mg/l Fe_mg/L Ni_mg/L Ca_mg/L Cu_mg/L
    -start
10 PUNCH TOT("S(6)")*96.06*1000
20 PUNCH TOT("Mg")*24.305*1000
30 PUNCH TOT("Fe")*55.845*1000
40 PUNCH TOT("Ni")*58.693*1000
50 PUNCH TOT("Ca")*40.078*1000
60 PUNCH TOT("Cu")*63.546*1000
-end

USER_GRAPH 1                #Concentration in cell
    -headings               pH
    -axis_titles            "Time (Weeks)" "pH"
    -chart_title            "RA-01_pH in 70 weeks WITH KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, -LA("H+"), color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
USER_GRAPH 2                #Concentration in cell
    -headings               pH
    -axis_titles            "Time (Weeks)" "pH"
    -chart_title            "RA-01_pH in 70 weeks with kinetics"
    -initial_solutions      true
    #-axis_scale y_axis      0 10         #concentration   
    #-axis_scale x_axis       0 30 5 1    #cell number     
    #-axis_scale sy_axis     6.6 7.1      #pH scale
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
#5 if cell_no < 30 then goto 1000
10 GRAPH_X total_time/604800
20 GRAPH_Y -LA("H+")
1000 REM end
  -end

USER_GRAPH 2                #Solute concentration in column
    -headings               SO4 Mg  Fe Ni Ca
    -axis_titles            "Time (Weeks)" "Concentration (mg/L)"
    -chart_title            "RA-01_Solute Concentration in 70 weeks_with KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, TOT("S(6)")*96.02*1000, color = Red, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end

USER_GRAPH 3               #Solute concentration in column
    -headings               Mg
    -axis_titles            "Time (Weeks)" "Concentration (mg/L)"
    -chart_title            "RA-01_Solute Concentration in 70 weeks_with KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, TOT("Mg")*24.305*1000, color = Green, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
-end

USER_GRAPH 5               #Solute concentration in column
    -headings               Ni
    -axis_titles            "Time (Weeks)" "Concentration (mg/L)"
    -chart_title            "RA-01_Solute Concentration in 70 weeks_with KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, TOT("Ni")*58.6934*1000, color = Yellow, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end

[dlparkhurst]

First, are you sure TRANSPORT is working the way you expect? At each shift, solution 0 is advected into cell 1, and then 2.8 mol of water are removed. So at all times, the volume of water is about 55.5 - 2.8 = 52.7 mol water or about 0.95 kgw. So, the loss of water is not cumulative.

Then, I think you need to be more careful with your PHASES definition. You need to define log K, if you use SI in the RATES definitions.

Also, the definition for Pyrrhotite has e- in the definition. I think this causes some strange redox and pH reactions. I would normalize so that the cations add up to 1.0. I have simply increased the iron stoichiometry. I am not sure, but I suspect that the stoichiometry of Pyrrhotite is the reason for the odd blips in your calculation. So consider what you want for PHASES definitions, the pyrrhotite stoichiometry will affect the pH and the rates of reactions.

Lastly, your rate for Forsterite is fast. I have added factors 1 - SR so that rates slow down and approach equilibrium Even so, the entire pool of forsterite will dissolve with the changes that I made. At each shift a lot dissolves and is then advected out of the system at the next shift.

TITLE RA-04
SOLUTION 0 Pure water
    temp      25
    pH        7.0
EQUILIBRIUM_PHASES 0       # Equilibrates with the atmosphere
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
SAVE SOLUTION 0
END

PHASES
Pyrrhotite                             
        #FeS + 1.0000 H+  =  + 1.0000 Fe++ + 1.0000 HS-
        Fe0.989Ni0.011S = +0.989Fe+2   +0.011Ni+2  -1.000H+   + 1.00HS- 
        log_k           -3.7193
#Pyrrhotite                             
#      Fe0.903Ni0.011S = +0.903Fe+2   +0.011Ni+2  -1.000H+   + 1.00HS-   -0.172e-
#   log_k        0
Forsterite                             
   Mg1.675Fe0.32Ni0.005SiO4 + 4H+ = 2H2O + SiO2 + 1.675Mg+2 + 0.32Fe+2 + 0.005Ni+2 
   log_k           27.8626
#   log_k     0
Enstatite                             
   Mg1.686Fe0.322SiO3 + 2H+ = 1H2O + SiO2 + 1.686Mg+2 + 0.322Fe+2  + 2.016e-
   log_k   0
END

SOLUTION 1 Rana (column)   # fifth rinse; initial solution in column (fixed pO2)
    temp      20
    pH        4.02     
    pe        9.58         
    redox     pe
    units     mg/l
    density   1
    Alkalinity 0.1
    Ca        8.87             
    Cu        0.1515         
    Fe        0.0121       
    Mg        40.899
    S(6)      228.195       
    Si        10
    Ni        2.222
    -water    0.9
EQUILIBRIUM_PHASES 1
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
    Fe(OH)3 0 10
    #Gypsum 0 0
    #FeSO4 0 0.01

SAVE SOLUTION 1
END

RATES
Pyrrhotite
-start
1 rem PARM(1) = log10(A/V, 1/dm)
2 rem PARM(2) = exp for (m/mo)
3 rem PARM(3) = exp for H+
4 rem PARM(4) = exp for O2
10 if (m <= 0) then goto 200
#20 if (SI("Pyrrhotite") >= 0) then goto 200
30 lograte = -5.79 + PARM(1) + PARM(2)*log10(m/m0) + PARM(3)*lm("H+") + PARM(4)*lm("O2")
40 moles = (10^lograte)*TIME * (1 - SR("Pyrrhotite"))
#50 if (moles > m) then moles = m
200 SAVE moles
-end
Forsterite
-start
1 rem PARM(1) = log10(A/V, 1/dm)
2 rem PARM(2) = exp for (m/mo)     
3 rem PARM(3) = exp for H+
10 if (m <= 0) then goto 200
#20 if (SI("Forsterite") >= 0) then goto 200
30 lograte = -5.96 + PARM(1) + PARM(2)*log10(m/m0) + PARM(3)*lm("H+")    ## k from Inverse modelling
40 moles = (10^lograte)*TIME * (1 - SR("Forsterite"))
#50 if (moles > m) then moles = m
200 SAVE moles
-end
END

KINETICS 1
#PARM(2) is of 2/3 for uniformly dissolving spheres/cubes(Parkhurst and Appelo, 1999)
Pyrrhotite
   #-formula Fe0.903Ni0.011S
      -m 0.25                 
      -m0 0.25
      -PARMS -0.71 0.67 0.25 0.1 # 0.1 is the optimum O2 reaction order
   #-tol 1e-9               
Forsterite
   #-formula Mg1.675Fe0.32Ni0.005SiO4
      -m 10.71                         
      -m0 10.71
      -parms 3.09 0.67 0.28             
   #-tol 1e-9
-cvode
-step_divide 1
-runge_kutta 6
-bad_step_max 1000000000
END

EXCHANGE 1
-equilibrate 1
X 0.1   #tailings impoundments CEC can be 1-10 meq/100g
END

SURFACE 1
       Hfo_wOH Fe(OH)3 equilibrium_phase 0.2 54300 # 0.3 800 18.5 
       Hfo_sOH Fe(OH)3 equilibrium_phase 0.005
      -equilibrate with solution 1
END

REACTION 1
H2O -1.0 # -1.0 is the relative moles of the reaction.
2.8 moles #this is the number of moles to be removed, but  must be done every week; assume only half is gone
END

TRANSPORT
    -cells                 1
    -shifts                65     
    -time_step             604800     
    -lengths               0.12     
    -dispersivities        7.92e-4
    -correct_disp          true
    -diffusion_coefficient 1e-9 #m2/s
    -print_cells           1
    -punch_cells           1   
    -flow_direction        forward
    -boundary_conditions   flux flux

USER_PUNCH 1
    -headings Sulfate_mg/l Mg_mg/l Fe_mg/L Ni_mg/L Ca_mg/L Cu_mg/L
    -start
10 PUNCH TOT("S(6)")*96.06*1000
20 PUNCH TOT("Mg")*24.305*1000
30 PUNCH TOT("Fe")*55.845*1000
40 PUNCH TOT("Ni")*58.693*1000
50 PUNCH TOT("Ca")*40.078*1000
60 PUNCH TOT("Cu")*63.546*1000
-end

USER_GRAPH 1                #Concentration in cell
    -headings               pH
    -axis_titles            "Time (Weeks)" "pH"
    -chart_title            "RA-01_pH in 70 weeks WITH KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, -LA("H+"), color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
USER_GRAPH 2                #Concentration in cell
    -headings               pH
    -axis_titles            "Time (Weeks)" "pH"
    -chart_title            "RA-01_pH in 70 weeks with kinetics"
    -initial_solutions      true
    #-axis_scale y_axis      0 10         #concentration   
    #-axis_scale x_axis       0 30 5 1    #cell number     
    #-axis_scale sy_axis     6.6 7.1      #pH scale
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
#5 if cell_no < 30 then goto 1000
10 GRAPH_X total_time/604800
20 GRAPH_Y -LA("H+")
1000 REM end
  -end

USER_GRAPH 2                #Solute concentration in column
    -headings               SO4 Mg  Fe Ni Ca
    -axis_titles            "Time (Weeks)" "Concentration (mg/L)"
    -chart_title            "RA-01_Solute Concentration in 70 weeks_with KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, TOT("S(6)")*96.02*1000, color = Red, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end

USER_GRAPH 3               #Solute concentration in column
    -headings               Mg
    -axis_titles            "Time (Weeks)" "Concentration (mg/L)"
    -chart_title            "RA-01_Solute Concentration in 70 weeks_with KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, TOT("Mg")*24.305*1000, color = Green, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
-end

USER_GRAPH 5               #Solute concentration in column
    -headings               Ni
    -axis_titles            "Time (Weeks)" "Concentration (mg/L)"
    -chart_title            "RA-01_Solute Concentration in 70 weeks_with KINETICS"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 PLOT_XY total_time/604800, TOT("Ni")*58.6934*1000, color = Yellow, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end

[rfembilejr]

I believe so, but now, I am thinking twice. This simulation is for a kinetic testing column that is rinsed with DI water once a week. So the 2.8 moles of H2O is what I infer to be evaporating every week. The 0.95 kgw is simply the mass of water collected. We collect roughly 0.90 kgw leachate and the other 0.28 moles remain in the tailings pore water while the other 0.28 moles would evaporate.

I did not know I had to define log K. Since we are trying to use our own calculated rate law for forsterite and pyrrhotite dissolution, is it ok to use the log K from the database? I was thinking this log K may not be representative for the forsterite and pyrrhotite present in our tailings. The dissolution rates we are using are from our kinetic test experiment. This is probably why it seems faster than typical. If I revert back to my original RATES definition, would it still be ok? But then, I would not need the log k, right?

After running the script that you modified, it becomes more difficult to fit our actual kinetic test data. pH has become very high from the actual pH of only 4.5. Similarly, SO4, Mg and Ni are elevated in the first week and then drops to almost nil afterwards. I am attaching some plots of our actual kinetic test data and the prior model (with unusual peaks). As  you can see, they are still very off. I decided not to plot the data from your modified script because it is too high and will mask the kinetic test and Geo2 model.

Apologies for these issues. I have been working on these for too long now and thought I got it to work mostly, until now.
I highly appreciate any feedback.

[dlparkhurst]

I do not think you should use SI in the rate expression if you do not know the log K. I believe the rate for fosterite is very fast, so it will bump into SI > 0. In that case the rate is limited by your arbitrary log K and depends very little on the parameters of the rate expression.

I suspect you are right that the surplus of S in pyrrhotite allows more acid generating and a lower pH.

I tried quite a few things to try to understand the blips in your model. What seems to work is to add -donnan to the SURFACE definition. I think the sequence of transport reactions could cause charge to accumulate on the surface. Pure speculation, but there may be two numerical solutions to the problem when this charge accumulates and occasionally, a high-pH numerical solution is found. -donnan will require the surface plus donnan layer to be charge balanced throughout the calculation. See if that helps. (May not be necessary, but I would probably charge balance SOLUTION 1, which is used to determine the composition of the surface.)

[rfembilejr]

Sorry, I am a bit confused when you said forsterite will bump into SI>0 knowing that its dissolution is very fast. Should't it be the other way around? SI>0 would mean precipitation, at least that's what I thought. But now, I agree that if I use SI without log K value, my rate law parameters will be bypassed, which is not something I would like.

I am not quite familiar on how to implement -donnan on the SURFACE block other than the one in the SURFACE example from the usgs website.
Should it be something like this:

SURFACE 1
       Hfo_wOH Fe(OH)3 equilibrium_phase 0.2 54300
       Hfo_sOH Fe(OH)3 equilibrium_phase 0.005
      -equilibrate with solution 1
      -Donnan   #default is 10E-8
END

I also charge balanced SOLUTION one using SO4 since this was the issue with my inverse model, the solution has excess negative charge. With this set-up, the blips (unwanted peaks) disappear, but still, the concentrations of SO4, Ni and Mg are initially too high but tapers down halfway. The pH is still also too high.

HOWEVER, if I change the surface area of forsterite to almost the same as of pyrrhotite, keeping all else the same (with the donnan), I got almost the model that I want. Quite interesting! which would suggest the reactive surface area of forsterite may actually be smaller than what I thought/calculated.

I guess, I would need to play with the surface area/initial moles more, unless you have other suggestions.
Thanks a lot!

Rodrigo

[dlparkhurst]

Your original RATES caused the rate to be zero if SI >=0, so the reaction stops at SI=0. My point was that if the log K is artificially small, the SI will be artificially large and the rate will stop. If you remove the IF statement, then I think you will need a much smaller rate for forsterite.

I would simply add -donnan.

You will have to decide on the parameters, perhaps you want to use parameter estimation. My only concern was the odd blips.

[rfembilejr]

I see. Yes, I have performed parameter estimation as well. So I will just use that in conjunction with this.
Thanks again.