Kinetics and Rates with Transport Modelling of Tailings columns

[rfembilejr]

Hi,

I am working on a soil/tailings that I am leaching with pure water.
I have some models working (attached), however I just want to check if TRANSPORT, KINETICS and RATES block make sense.
I am using Pyrrhotite and Forsterites as the main phases in the tailings.
Looking at the resulting plots, if I use TRANSPORT only, the pH decrease by about 1 unit while if I use KINETICS, the pH drops to 4 initially and increased by half unit before steadily decreasing. If I add KINETICS-RATES and look at some species (Mg, Fe and SO4) through time, there is higher concentration of these species, I think this is probably due to the acidic pH generated from the dissolution of Pyrrhotite?
Attached are codes for TRANSPORT only and TRANSPORT with KINETICS that should generate the plots. I am using wateq4f database.

I also have some questions about the plots. I would like to model the pH and species concentration through time but I get pretty similar trends if I vary the time_step to 1, 10, 50 or 100 years. I was expecting I might be able to see some changes using different times.
When I adjust the time_step to 10 or 50 years (written in seconds), I don't see a change in the trend of pH and some species. Why is this happening? I am converting the X-axis units to the desired time (ie. 10 or 50 years) such that I don't see 10^6 for example but then again, the trend is pretty much the same. Is there something wrong with how I convert the total time?
I would appreciate if someone can also check my RATES BASIC code. I am not quite sure if I wrote the RATES for Pyrrhotite and Forsterite correctly.

I was also hoping to make a plot showing the changes in concentration of species and pH as you go down from the top of the column to the bottom but not really sure on how to proceed with this.

Thank you in advance for any help.

Rodrigo

[dlparkhurst]

I think you should skip all the output and TRANSPORT calculations and get a batch reaction working. You are not using PHASE correctly. First, the 1-10 in PHASES 1-10 is ignored. There is only one set of phases definitions for the stoichiometry and log K of a reaction. So your repeated definitions have no effect.

Defining something in PHASES does not produce a reaction. If you want pyrrhotite to react, then you must include it in REACTION, KINETICS, or EQUILIBRIUM_PHASES.

For example, the following would consider a system open to atmospheric oxygen. The input file simulates the oxidation of up to 10 mmol of pyrrhotite in this environment. I don't know your setting, but pyrite (or marcasite) may be more common than pyrrhotite. Also, you are not considering carbonates, which would be important if present, so I assume they are not present.

PHASES 1-10
Pyrrhotite
   Fe0.87S = +0.870Fe+2   -1.000H+   + 1.00HS-   -0.260e-
   log_k    -5.59   #dont matter in inverse
     
#Forsterite
#   Mg1.682Fe0.31Ni0.008SiO4 + 4H+ = H4SiO4 + 1.682Mg+2 + 0.31Fe+2 + 0.008Ni+2  #Balanced by averaging Olivines and manipulating :)
#      log_k            0
END
SOLUTION 1 Rana (column)   # very first rinse; initial solution in column (fized pO2)
    temp      20
    pH        5.29
    pe        4.32
    redox     pe
    units     mg/l
    density   1
    Alkalinity 4.14
    Ca        415.80
    Cl        40
    Cu(2)     0.026230
    Fe        0.0110
    Mg        731.97
    Na        57.464
    S(6)      4510.68
    Si        10
    -water    1 # kg
END
EQUILIBRIUM_PHASES 1-10
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
END
REACTION
Pyrrhotite 1
10 mmol in 10 steps
END
USER_GRAPH 2                #Concentration in cell
    -headings               pH
    -axis_titles            "Pyrrhotite reacted" "pH"
    -chart_title            "Concentration in 1 year_Transport only"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X RXN
20 GRAPH_Y -LA("H+")
1000 REM end
  -end
END
USE solution 1
USE equilibrium_phases 1
USE reaction 1
END

[rfembilejr]

Hi David,

First of all, thank you very much for your remarks.
1.  My columns contain tailings from magmatic Ni sulfide deposit which is mainly composed of Forsterite and Pyrrhotite. I rinse the column with pure water every week and collect a leachate. I also assume that the other minerals ie. pyroxene and amphiboles are very slow reacting that is why I am not including them in the PHASES or EQUILIBRIUM_PHASES. There is also very little chalcopyrite and almost negligible pyrite. I have defined the Pyrrhotite phases since it is not found in the wateq4f database.
The Forsterite stoichiometery is based on microprobe data so I just used the chemistry to come up with its chemical formula. Ideally this should not mess with the model as compared to just using Forsterite from the database, right? or should I just use Forsterite from the database which is just Mg2SiO4. I wanted to use my own chemical formula because it has Fe and Ni on it.
I think Pyrrhotite is oxidizing in the tailings while the Forsterite is also dissolving due to a decrease in pH.


2. If I get your point, I should first set up a batch reaction that is working. However, following from your example, how do I proceed with adding TRANSPORT and KINETICS? I guess, I am not quite sure how to set-up the code for this. I modified my previous code mainly by removing 1-10 from the PHASES block (which did not really make a difference). Skipping TRANSPORT could defeat my purpose of modelling the reactions within the column during rinsing. I also plan to relate this to actual field setting where the tailings are deposited and rainwater percolates through the tailings while changing different time horizons ie. 10 or 100 years. I have gone through several examples in this forum as well as phreeqc examples but setting up TRANSPORT together with KINETICS-RATES is very tricky. I would appreciate any suggestion or if there is a template on how to set up all these.


TITLE Transport with kinetics for RA-01


#### Infiltrating Solution ######

SOLUTION 0 Pure water
    temp      20
    pH        7.0
    pe        4
    redox     pe
    units     mg/l
    density   1
    Ca        0.001
    Cl        0.001
    Cu(2)     0.002
    Fe        0.001
    Mg        0.001
    Na        0.001
    S(6)      0.001
    Si        0.001
    -water    1 # kg

EQUILIBRIUM_PHASES 0       # Equilibrates with the atmosphere
    CO2(g)    -3.5 10
    O2(g)     -0.699 10

SAVE SOLUTION 0
END

##### Initial column conditions ############

PHASES
Pyrrhotite
   Fe0.87S = +0.870Fe+2   -1.000H+   + 1.00HS-   -0.260e-
   log_k    -5.59   #dont matter in inverse
     
Forsterite
   Mg1.682Fe0.31Ni0.008SiO4 + 4H+ = H4SiO4 + 1.682Mg+2 + 0.31Fe+2 + 0.008Ni+2  #Balanced by averaging Olivines and manipulating :)
      log_k            0
END


SOLUTION 1-10 Rana (column)   # very first rinse; initial solution in column (fized pO2)
    temp      20
    pH        5.29
    pe        4.32
    redox     pe
    units     mg/l
    density   1
    Alkalinity 4.14
    Ca        415.80
    Cl        40
    Cu(2)     0.026230
    Fe        0.0110
    Mg        731.97
    Na        57.464
    S(6)      4510.68
    Si        10
    -water    1 # kg
END
EQUILIBRIUM_PHASES 1-10
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
END

### RATES ###

RATES

Pyrrhotite
-start
1 rem PARM(1) = log10(A/V, 1/dm)
2 rem PARM(2) = exp for (m/mo)
10 if (m <= 0) then goto 200
20 if (SI("Pyrrhotite") >= 0) then goto 200
30 lograte = -10.04 + PARM(1) + PARM(2)*log10(m/m0) #+ PARM(3)*lm("Fe3+") + PARM(4)*lm("H+")    #from Palandri and Kharaka (2004) for monoclinic pyrrhotite
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 = -11.31 + PARM(1) + PARM(2)*log10(m/m0) + PARM(3)*lm("H+")    # from Palandri and Kharaka (2004) fo91 at 25C, pH=5.4, logK = -9.31 (Pokrovsky and Schott, 2000)
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 SAVE moles
-end
END

##### KINETICS #####

KINETICS 1-10
Pyrrhotite
   -formula Fe0.87S 1
      -m 1.49                 #from calculations similar to that of albite
      -m0 1.49
      -PARMS 2.27 0.1 #0.355 -0.597       #####first parm is from surface area calculated similar to albite.
-steps      1 10 100 1000 10000
-step_divide 1
-runge_kutta 6
-bad_step_max 500

Forsterite
   -formula Mg1.682Fe0.31Ni0.008SiO4
      -m 22.72                   #from calculations similar to that of albite
      -m0 22.72
      -parms 2.53 0.1 5.4        ######needs to be checked! #####not sure about the first and second parm here.
                  # Area (3.39)
      #-time 1.5 year in 40       #check this
-steps      1 10 100 1000 10000
-step_divide 1
-runge_kutta 6
-bad_step_max 500
END



##### TRANSPORT ####

TRANSPORT
    -cells                 10
    -shifts                10
    -time_step             31557600 # seconds , 1 yrs
    -lengths               0.06   #total column lenght is 0.3m
    -dispersivities        10*0.003 #0.0175L^1.46 (Flow velocity = cell length/time step; 0.3/0.125 = 2.4 m/day or 0.1 m/hr)
    -correct_disp          true
    -thermal_diffusion     10   3e-10
    -diffusion_coefficient 1e-9 #m2/s
    -print_cells           10
    -punch_cells           10
    -boundary_conditions   flux flux #at column ends



SELECTED_OUTPUT 1
    -file                 C:\Users\RodrigoJr\Desktop\Desktop 14Aug\PhD Project\PHREEQC Codes\RANA Experiments\TRANSPORT Codes\RA-01\RA-01_TEST_10 years no kinetics.xls
    -reset                true
    -pe                   false
    -totals               S(6)  Mg  Fe  Ca  Alkalinity  Ni
    -saturation_indices   Forsterite  Pyrrhotite  Fe(OH)3(a)  Goethite
    -kinetic_reactants    Forsterite  Pyrrhotite


### GRAPHS ###

USER_GRAPH 1                #Concentration in cell
    -headings               Concentration  Mg  Fe Ni pH S(6)
    -axis_titles            "Time (Days)" "Concentration(mg/L)" #"pH"
    -chart_title            "Concentration in 1 year 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/1000
20 GRAPH_Y TOT("S(6)")
30 GRAPH_Y TOT("Mg")
40 GRAPH_Y TOT("Fe")
50 GRAPH_Y TOT("Ni")
#60 GRAPH_SY -LA("H+")
1000 REM end
  -end


USER_GRAPH 2                #Concentration in cell
    -headings               pH
    -axis_titles            "Time (Days)" "pH"
    -chart_title            "Concentration in 1 year 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/1000
20 GRAPH_Y -LA("H+")
1000 REM end
  -end

END

[dlparkhurst]

I would say the next step is to replace the REACTION definition with your kinetic definition and let the reaction proceed for varying amounts of time. . I am not sure what TRANSPORT would give you that this batch reaction does not, considering that the reactions are the same in all cells.

[rfembilejr]

Thanks David.
I will work on this and see how it goes. Will post an update soon.

[rfembilejr]

Hi David,

I have now moved forward with the kinetics part of my models. I just want it checked for errors.
I have several questions while running the code and about the script and results in general. The script I am working on is attached.

####Just a quick overview of my project since it has been a while:
I am working with kinetic testing columns containing nickel sulfide ore. The tailings is 50% forsterite, 19% orthopyroxene (enstatite), 8.5% amphibole, plagioclase feldspar (11.5%), chlorite-clay minerals (8%), 2% pyrrhotite and some minor chalcopyrite. The columns are rinsed with DI water weekly for 70 weeks and the leachates analyzed for major ions (Ca, Mg, SO4) and metals (Ni, Fe, Cu). 
We think forsterite is the main dissolving mineral along with pyrrhotite. Also, based on reaction rates from the literature, the other silicates are dissolving by at least an order of magnitude less than forsterite. Although very low in composition, pyrrhotite is relatively fast reacting. However, this sulfide is usually enclosed within the forsterite, but once liberated, it will react quite fast. This is why in the script below, we are only using forsterite and pyrrhotite in the rates and kinetics. ####

1) When I was running the script, you can see the following warnings (also appears in the output).
-Negative moles in solution for Mg, (some numbers), Recovering
-Maximum iterations exceed
-Numerical method failed with this set of convergence parameters

The run is usually completed and I get the plots that seem to make sense. Should I worry about these warnings or it is just the program trying to fix it? Is there a way to avoid these especially the first one, if that is something of major concern?

2) For the RATES block, I am following Palandri and Kharaka's rate law format. For pyrrhotite there is no rate dependency to Fe3+ that I am incorporating in the rate law. I am using a "k" that was obtained from the plot of log(rate) vs log(H+) which we developed using INVERSE_MODELING. The reaction order with respect to pH was used as PARM(3). The same process is was used for forsterite.
The rate law equations we are using are as follows:

Rate (forsterite) = -8.2651 (H+)^0.413
Rate (pyrrhotite) = -7.785 (H+)^0.388

We think these values are more accurate to use for our modeling purposes instead of using the laboratory rates found in the literature which are typically much faster.

3) So our input solution is deionized water (we will also change this to rain water when applied to field setting) and for the solution in the column (or tailings itself), we are currently using the solution chemistry of the very first rinse to model subsequent trends in solute concentration and pH. Is this reasonable? Based on our actual kinetic testing data trend, the pH seem to reflect the trend, although not exactly the same values. However, the solute concentrations are kind of not really jiving with the ones we observe in the actual concentrations in time. For example, Mg started off with almost zero concentration and increases to about 580 mg/L and maintains that until week 70. For SO4, the values are kind of very high.

4) For the KINETICS block, what should be the PARM(2) or the exp for (m/m0)? I have seen values equal to 2/3 for uniformly dissolving spheres, which is what I am using here.

5) To calculate my values for m and m0, I used the approach in the calculation for KINETICS data block for albite in the database of llnl.dat. Is this reasonable?

6) I have had error regarding bad_step_max, and the suggestion was to increase the value, it was originally 500 and so I changed it to 10000 and it works, is this ok? I am not really familiar how runge_kutta and step_divide works, despite reading about it in the examples.

7) I have five solutes (SO4, Mg, Fe, Ni, Ca, Cu) indicated in my solution composition and I was trying to modify the plots in terms of symbols or color but it doesn't seem to work. I know I can also replot these data in excel but just wanted to know any tricks on how to solve this.

8) The concentration of solutes are automatically plotted as moles? And so in order for me to see values in mg/L, I am just multiplying them with MW and 1000 (pls. see USER_GRAPH block).

9) I would like to plot solute concentration changes in with depth but I could not really make it work using the USER_GRAPH block. Any examples would be very helpful.

Thank you in advance for your help.

Rodrigo

[dlparkhurst]

(1) These are just warnings. If something serious happens, PHREEQC will send an "ERROR" message and the calculations will stop.

(2) Seems justified to me.

(3) If your model does not match your data, then you will probably make some adjustments.

(4) Unless you dissolve a substantial fraction of the number of moles initially present, the factor is nearly 1, so the exponent is not very sensitive.

(5) If you have xray quantification, then you could use that to estimate the moles of mineral. Otherwise, you do what you can. M and M0 are usually equal in the KINETICS definition.

(6) If the program does not throw an ERROR, I think the integration was successful. The RK integration estimates the error over each integration time substep. If it is too large, then the step is repeated with smaller time substeps. 500 was the number of times it failed to meet the error tolerance. Step_divide enforces a smaller initial sub time step.

If RK is having a hard time (many bad steps), you may also want to try using the -cvode option. It is an implicit method that can be more efficient for some (stiff) sets of rates.

(7) Use PLOT_XY (rather than GRAPH_X and GRAPH_Y) to adjust colors and symbols. See the USER_GRAPH documentation.

(8) OK. You can use GFW("SO4") to get the gram formula weight for sulfate.

(9) I think you could do it with some work. You probably need to use the PUT and GET statements combined with CELL_NO to transfer information from one cell to the next. Might be easier to use Excel or R.

[rfembilejr]

Thanks David!
I really appreciate your quick response.
At least, now I know I am in the right direction.