Code for sulfure oxidation (DMA) simulation

[SCH]

HI All,
See below a code I tried, but my goal is not reach. I want to get 0.0027 mol/l of SO4 and 0.0016 of Ca after 13 days.
Can you check for me and let me know if I need to modify something.
Thanks

############
TITLE Transport with kinetics


######################################   INFILTRATING SOLUTION   #####################################################

SOLUTION 0 # Distilled water
     temp     25
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1
    Al        1e-20
    Ca        2.2e-05
    Cl(-1)    0.0001
    Cu(2)     1e-10
    Fe(2)     1e-10
    Fe(3)     2e-05
    K         2.2e-05
    Mg        8.2e-07
    Na        2e-05
    Ni        1e-10
    S(-2)     1e-12
    S(6)      1e-05
    -water    1 # kg

EQUILIBRIUM_PHASES 0
    CO2(g)    -3.5 0.001
    O2(g)     -0.699 0.000
SAVE SOLUTION 0
END



SOLUTION 1-50 Rana (column)
     temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1
    Al        1e-20
    Ca        2.2e-05
    Cl(-1)    0.0001
    Cu(2)     1e-10
    Fe(2)     1e-10
    Fe(3)     2e-05
    K         2.2e-05
    Mg        8.2e-07
    Na        2e-05
    Ni        1e-10
    S(-2)     1e-12
    S(6)      1e-05
    -water    1 # kg
   
EQUILIBRIUM_PHASES 1-50
    CO2(g)    -3.5 0.001
    Gibbsite  0 0.1
    Goethite  0 0.1
    Magnesite 0 0
    O2(g)     -0.699 0.01
    Tenorite  0 0
SAVE SOLUTION 1-2


###############################################  RATES  ###############################################################
### Rmin = r_min *(A0/V)*(m/m0)^n*(1-SR "mineral")
### For kinetically controlled reactions, the solution is far from equilibrium thus, the term (1-SR) is approx 1.

RATES

Pyrite
-start
1 rem   parm(1) = log10(A/V, 1/dm)   parm(2) = exp for (m/m0)
2 rem   parm(3) = exp for O2      parm(4) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Pyrite") >= 0) then goto 200
30 lograte = -11.48*parm(1)*parm(3)*lm("O2") * parm(4)*lm("H+") * parm(2)*log10(m/m0)
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
-end
#
Ilmenite
-start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Ilmenite") >= 0) then goto 200
30 lograte = -8.35*PARM(1)*PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
-end
#
        Hematite
-start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Hematite") >= 0) then goto 200
30 lograte = -9.39 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
-end
#
        Labradorite
        -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Labradorite") >= 0) then goto 200
30 lograte = -11 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end
#
         Enstatite
          -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Enstatite") >= 0) then goto 200
30 lograte = -11 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end
#
         Orthoclase
         -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Orthoclase") >= 0) then goto 200
30 lograte = -11 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end
#
          Chlorite
          -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Chlorite") >= 0) then goto 200
30 lograte = 13 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end
#
         Spinel
        -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Spinel") >= 0) then goto 200
30 rate = -36.3 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
         -end
#
         Rutile
        -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Rutile") >= 0) then goto 200
30 lograte = 9.6 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
        -end
#
        muscovite
        -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("muscovite") >= 0) then goto 200
30 lograte = -13 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
        -end
#
        Biotite ###########
        -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("biotite") >= 0) then goto 200
30 lograte = 43.3 * PARM(1) * PARM(2)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
        -end
#
       Chalcopyrite ###########
        -start
1 rem   parm(1) = log10(A/V, 1/dm)   parm(2) = exp for (m/m0)
2 rem   parm(3) = exp for O2      parm(4) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("chalcopyrite") >= 0) then goto 200
30 lograte = -32.6 + PARM(1) * PARM(2)*log10(m/m0) * PARM(3)*lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
         -end

##############################################  KINETICS  ###########################################################
### Calculations similar to that of Albite in Database of llnl.dat
# Assume tailings major components that are reacting to generate the leachate chemistry are Forsterite and Pyrrhotite
# Assume tailings is 51.19% Forterite by mass in 1 mm spheres (radius 0.5 mm)
# Assume tailings is 1% Pyrrhotite, which accounts for those that are liberated and exposed to react.
# Assume density of tailings is 2600 kg/m^3 = 2.6 kg/L. This is primarily due to the high concentration (about 50%) of forsterite which makes bulk
# of the tailings
## second parm is assumed to be 2/3 for uniformly dissolving spheres and

KINETICS 0-50
Pyrite
    -formula  FeS2  1
    -m        0.005
    -m0       0.005
    -parms    1.03 0.67 0.5 -0.11
-steps       86400 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End


Biotite
    -formula  KMg2FeAlSi3O10(OH)2  1
    -m        0.005
    -m0       0.005
    -parms    1.58 0.25
-steps       3600 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End

Chalcopyrite
    -formula  CuFeS2  1
    -m        0.0009
    -m0       0.0009
    -parms    0.5 0.66 0.5 -0.11
   -steps       3600 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End

muscovite
    -formula  KAl2(AlSi3O10)(OH)2  1
    -m        0.0027
    -m0       0.0027
    -parms    1 0.37
 -steps       3600 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End

Chlorite
    -formula  (Mg2Fe3Al2)Si3O10(OH)8  1
    -m        0.01
    -m0       0.01
    -parms    1.89 0.25
-steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End
Labradorite
    -formula  (CaNa)Al3Si5O16  1
    -m        0.21
    -m0       0.21
    -parms    3.2 0.5
    -steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End
Enstatite
    -formula  Mg1.686Fe0.322SiO3  1
    -m        0.044
    -m0       0.044
    -parms    2.21 0.6
-steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End

Hematite
    -formula  Fe2O3  1
    -m        0.056
    -m0       0.056
    -parms    2 1
-steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
End
Ilmenite
    -formula  FeTiO3  1
    -m        0.18
    -m0       0.18
    -parms    2.7 0.42
-steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
Orthoclase
    -formula  (KAl)Si3O8  1
    -m
       0.03
    -m0       0.03
    -parms    2.12 -0.2
-steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000

Spinel
    -formula  Al2MgO4  1
    -m        0.003
    -m0       0.003
    -parms    0.83 0.3
-steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000

Rutile
    -formula  TiO2  1
    -m        0.001
    -m0       0.001
    -parms    0.26 0.3
    -tol      1e-08
-steps       600000 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 5000000
END

#################################################   SURFACE  ###########################################################

SURFACE 1-50
        Hfo_w 100  600   1-50     
        -equilibrate with solution 1-5
END

 #First number is site concentration (mol/L)
 #Second number is surface area (m2/g)             is originally 600
 #Third number is mass of sorbent (g/L)

##############################################  TRANSPORT  ###########################################################
## Column height is 0.12 m

TRANSPORT
    -cells                 50
    -shifts                1000
    -time_step             604800 # seconds = 7 days
    -boundary_conditions   flux constant
    -lengths               50*0.5
    -correct_disp          true
    -diffusion_coefficient 1e-09
    -thermal_diffusion     2   1e-09
    -print_cells           50
    -punch_cells           50



SELECTED_OUTPUT 1
    -file                 :)
    -reset                true
    -pe                   false
    -totals               S(6)  Mg  Fe  Ca  Ni  Alkalinity
    -saturation_indices   Fe(OH)3  Forsterite  Pyrrhotite  Goethite
    -kinetic_reactants    Forsterite  Pyrrhotite  Goethite  Labradorite
                          muscovite  Orthoclase  Pyrite

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

##############################################  GRAPHS  ##############################################################


# pH #############################################  pH  ##############################################################

USER_GRAPH 1
    -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/10000, -LA("H+"), color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 2
    -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/10000, TOT("S(6)")*96.02*1000, color = Red, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true


USER_GRAPH 3
    -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/10000, TOT("Mg")*24.305*1000, color = Green, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 4
    -headings               Fe
    -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/10000, TOT("Fe")*55.845*1000, color = Orange, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 5
    -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/10000, TOT("Ni")*58.6934*1000, color = Yellow, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 6
    -headings               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/10000, TOT("Ca")*40.078*1000, color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 7
    -headings               Cu
    -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/10000, TOT("Cu")*63.546*1000, color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

[dlparkhurst]

My guess is that the major reactions are pyrite oxidation, iron oxyhydroxide precipitation, and calcite dissolution. Oxidation of 1.35 mmol of pyrite, precipitation of all the iron, and dissolution of 1.6 mmol of calcite would produce the concentrations that you give. I doubt many of the other minerals you have included in equilibrium phases (magnesite) and kinetics will react much over a few weeks.

Source of the trace metals will be hard to determine between sorption reactions and mineral dissolution unless their concentrations are significant.

[SCH]

Hi Dr. Thanks for your feedback.

Since 2 days I tried something but I think that It have not reaction in my cell beacause the concentrations on the graph of SO4-2 and Ca are the same with the solution concentration.

How I can do to get a reaction in my column. Thanks for all





############
TITLE Transport with kinetics for RA-01_ACTUAL COLUMNS_not field setting
#### All mineral components with 5th rinse
#### Using BET method for WEATHERED grains


######################################   INFILTRATING SOLUTION   #####################################################

SOLUTION 0 # Distilled water
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1
    Al        1e-20
    Ca        2.2e-05
    Cl(-1)    0.0001
    Cu(2)     1e-10
    Fe(2)     1e-10
    Fe(3)     2e-05
    K         2.2e-05
    Mg        8.2e-07
    Na        2e-05
    Ni        1e-10
    S(-2)     1e-12
    S(6)      1e-05
    -water    1 # kg

EQUILIBRIUM_PHASES 0       ### Equilibrates with the atmosphere
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
 
SOLUTION 1-100 Rana (column)
     temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1
    Al        1e-20
    Ca        2.2e-05
    Cl(-1)    0.0001
    Cu(2)     1e-10
    Fe(2)     1e-10
    Fe(3)     2e-05
    K         2.2e-05
    Mg        8.2e-07
    Na        2e-05
    Ni        1e-10
    S(-2)     1e-12
    S(6)      1e-05
   
EQUILIBRIUM_PHASES 1-100       ### Equilibrates with the atmosphere
    CO2(g)    -3.5 10
    O2(g)     -0.699 10
   
COPY solution 0 100 # for use later on, and in
COPY solution 1 101 # 20 cells model

END

INCREMENTAL_REACTIONS true
KINETICS 1-10
Pyrite
    -formula  FeS2  1
    -m        0.00135
    -m0       0.00135
    -parms    -5.65 0.67 0.5 -0.11

 Labradorite
    -formula  (CaNa)Al3Si5O16  1
    -m         0.0016
    -m0        0.0016
    -parms    -3.38 6.626

Chlorite(14A)
    -formula  (Mg2Fe3Al2)Si3O10(OH)8  1
    -m        0.0002
    -m0       0.0002
    -parms    -5.93 0.5

COPY kinetics 1 101 # to use with 20 cells               
END

###############################################  RATES  ###############################################################
### Rmin = r_min *(A0/V)*(m/m0)^n*(1-SR "mineral")
### For kinetically controlled reactions, the solution is far from equilibrium thus, the term (1-SR) is approx 1.

RATES

       Pyrite
       -start
1 rem   parm(1) = log10(A/V, 1/dm)   parm(2) = exp for (m/m0)
2 rem   parm(3) = exp for O2      parm(4) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Pyrite") >= 0) then goto 200
30 lograte = 24.65 * parm(1) * parm(3)+lm("O2") * parm(4)+lm("H+") * parm(2)+log10(m/m0)
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
        -end
#
      Labradorite
        -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Labradorite") >= 0) then goto 200
30 lograte = -7.87 * PARM(1) * PARM(2)+lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end
 Chlorite(14A)
          -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Chlorite(14A)") >= 0) then goto 200
30 lograte = -11.11 * PARM(1) * PARM(2)+lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end

##############################################  KINETICS  ###########################################################
### Calculations similar to that of Albite in Database of llnl.dat
# Assume tailings major components that are reacting to generate the leachate chemistry are Forsterite and Pyrrhotite
# Assume tailings is 51.19% Forterite by mass in 1 mm spheres (radius 0.5 mm)
# Assume tailings is 1% Pyrrhotite, which accounts for those that are liberated and exposed to react.
# Assume density of tailings is 2600 kg/m^3 = 2.6 kg/L. This is primarily due to the high concentration (about 50%) of forsterite which makes bulk
# of the tailings
## second parm is assumed to be 2/3 for uniformly dissolving spheres and

#################################################   SURFACE  ###########################################################

SURFACE 1-3
        Hfo_w 100  600   1-3     
        -equilibrate with solution 1-10
END


#time step is 1/3 of a week since we have 3 cells and one full cell equals one week of movement. So since three cells, each cell has to have a third
#of the total time step of 1 week
#shift is 3 times 70 since it moves the water 3 times in a column and we have 70 columns
#flow velocity will become 1.98E-7 m/s
#the actual flow is not equivalent to the one that we are modelling, bc there are other reactions occuring in the microflow that contributes to
#the chemistry

## Total simulation time = shifts*time_step :
## Flow velocity = cell length/time step;
## Dispersivity = 0.0175L^1.46 where L is column height, = 1.32e-3
## Shifts = Total simulation time/time step
## Pore Volume = shifts/cells:
## Time step: Pore volume/Irrigation rate

TRANSPORT
    -cells                 20
    -shifts                13
    -lengths               20*0.125
    -dispersivities        20*0.05
    -print_cells           20
    -print_frequency       10
    -punch_cells           20
    -punch_frequency       2
COPY cell 101 0
END
TRANSPORT
    -shifts                55
    -time_step             604800 1 # seconds
    -lengths               20*0.125
    -dispersivities        20*0.05
    -print_cells           20
    -punch_cells           20
    -multi_d               true 2.39e-09 0.35 0 1

##############################################  OUTPUT  ##############################################################

SELECTED_OUTPUT 1
    -file                  :)
    -reset                true
    -pe                   false
    -solution             true
    -distance             true
    -totals               S(6)  Mg  Fe  Ca  Ni  Alkalinity
    -saturation_indices   Fe(OH)3  Forsterite  Pyrrhotite  Goethite
    -kinetic_reactants    Forsterite  Pyrrhotite  Goethite   


##############################################  USER_PUNCH  ##############################################################

USER_PUNCH 1
    -headings Sulfate_mole/l Mole_mg/l Fe_mole/L Ni_mole/L Ca_mole/L Cu_mole/L
    -start
10 PUNCH TOT("S(6)")
20 PUNCH TOT("Mg")
30 PUNCH TOT("Fe")
40 PUNCH TOT("Ni")
50 PUNCH TOT("Ca")
60 PUNCH TOT("Cu")
-end

##############################################  GRAPHS  ##############################################################


# pH #############################################  pH  ##############################################################

USER_GRAPH 1
    -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/1000, -LA("H+"), color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 2
    -headings               SO4 Mg Fe Ni Ca
    -axis_titles            "Time (Weeks)" "Concentration (mole/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/10000, TOT("S(6)"), color = Red, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true


USER_GRAPH 3
    -headings               Mg
    -axis_titles            "Time (Weeks)" "Concentration (mole/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/10000, TOT("Mg"), color = Green, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 4
    -headings               Fe
    -axis_titles            "Time (Weeks)" "Concentration (mole/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/10000, TOT("Fe"), color = Orange, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 5
    -headings               Ni
    -axis_titles            "Time (Weeks)" "Concentration (mole/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/10000, TOT("Ni"), color = Yellow, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 6
    -headings               Ca
    -axis_titles            "Time (Weeks)" "Concentration (mole/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/10000, TOT("Ca"), color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

USER_GRAPH 7
    -headings               Cu
    -axis_titles            "Time (Weeks)" "Concentration (mole/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/10000, TOT("Cu"), color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end
    -active                 true

End

[dlparkhurst]

I suggest you try to debug a simpler setup.

I don't know your experimental setup, but I assume here that the column is closed to oxygen; oxygen only enters with solution 0. It doesn't matter because your pyrite rate is approximately 0 (see results of PRINT statement added to the RATES definition), so I think you have not coded the rate equation correctly.

If pyrite is oxidized, I think you need to include a ferric oxyhydroxide phase.

You will have to decide whether you have picked the correct reactants for your system.

Code: [Select]


############
TITLE Transport with kinetics for RA-01_ACTUAL COLUMNS_not field setting
#### All mineral components with 5th rinse
#### Using BET method for WEATHERED grains

###############################################  RATES  ###############################################################
### Rmin = r_min *(A0/V)*(m/m0)^n*(1-SR "mineral")
### For kinetically controlled reactions, the solution is far from equilibrium thus, the term (1-SR) is approx 1.

RATES

       Pyrite
       -start
1 rem   parm(1) = log10(A/V, 1/dm)   parm(2) = exp for (m/m0)
2 rem   parm(3) = exp for O2      parm(4) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Pyrite") >= 0) then goto 200
30 lograte = 24.65 * parm(1) * parm(3)+lm("O2") * parm(4)+lm("H+") * parm(2)+log10(m/m0)
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
210 PRINT "Pyrite ", moles, 24.65 * parm(1) * parm(3), lm("O2") * parm(4), lm("H+") * parm(2), log10(m/m0)
        -end
#
      Labradorite
        -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
#20 if (si("Labradorite") >= 0) then goto 200
30 lograte = -7.87 * PARM(1) * PARM(2)+lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end
 Chlorite(14A)
          -start
1 rem   parm(1) = log10(A/V, 1/dm)
2 rem   parm(2) = exp for H+
10 if (m <= 0) then goto 200
20 if (si("Chlorite(14A)") >= 0) then goto 200
30 lograte = -11.11 * PARM(1) * PARM(2)+lm("H+")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 save moles
          -end
######################################   INFILTRATING SOLUTION   #####################################################

SOLUTION 0 # Distilled water

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

KINETICS 1-20
-cvode
Pyrite
    -formula  FeS2  1
    -m        0.00135
    -m0       0.00135
    -parms    -5.65 0.67 0.5 -0.11

 Labradorite
    -formula  (CaNa)Al3Si5O16  1
    -m         0.0016
    -m0        0.0016
    -parms    -3.38 6.626

Chlorite(14A)
    -formula  (Mg2Fe3Al2)Si3O10(OH)8  1
    -m        0.0002
    -m0       0.0002
    -parms    -5.93 0.5
       
END
EQUILIBRIUM_PHASES 1-20
Goethite 0 0

##############################################  KINETICS  ###########################################################
### Calculations similar to that of Albite in Database of llnl.dat
# Assume tailings major components that are reacting to generate the leachate chemistry are Forsterite and Pyrrhotite
# Assume tailings is 51.19% Forterite by mass in 1 mm spheres (radius 0.5 mm)
# Assume tailings is 1% Pyrrhotite, which accounts for those that are liberated and exposed to react.
# Assume density of tailings is 2600 kg/m^3 = 2.6 kg/L. This is primarily due to the high concentration (about 50%) of forsterite which makes bulk
# of the tailings
## second parm is assumed to be 2/3 for uniformly dissolving spheres and

#################################################   SURFACE  ###########################################################



#time step is 1/3 of a week since we have 3 cells and one full cell equals one week of movement. So since three cells, each cell has to have a third
#of the total time step of 1 week
#shift is 3 times 70 since it moves the water 3 times in a column and we have 70 columns
#flow velocity will become 1.98E-7 m/s
#the actual flow is not equivalent to the one that we are modelling, bc there are other reactions occuring in the microflow that contributes to
#the chemistry

## Total simulation time = shifts*time_step :
## Flow velocity = cell length/time step;
## Dispersivity = 0.0175L^1.46 where L is column height, = 1.32e-3
## Shifts = Total simulation time/time step
## Pore Volume = shifts/cells:
## Time step: Pore volume/Irrigation rate

TRANSPORT
    -time_step             604800
    -cells                 20
    -shifts                1
    -lengths               20*0.125
    -dispersivities        20*0.05
    -print_cells           1-20
    -print_frequency       1
    -punch_cells           1-20
    -punch_frequency       1
USER_GRAPH 1
    -headings               time SO4 Ca Fe3 Fe2
    -axis_titles            "Days" "Molality" ""
  -start
10 GRAPH_X dist #total_time / 86400
20 GRAPH_Y TOT("S(6)"), TOT("Ca"), TOT("Fe(3)"), TOT("Fe(2)")
  -end
USER_GRAPH 2
    -headings               time Pyrite Labradorite Chlorite
    -axis_titles            "Days" "DELTA MOLES" ""
  -start
10 GRAPH_X dist #total_time / 86400
20 GRAPH_Y KIN_DELTA("Pyrite"), KIN_DELTA("Labradorite"), KIN_DELTA("Chlorite(14A)")
  -end


[SCH]

Hi Dr.
Thank you for your feedback and help
It is is to simulate a deposit of mining waste rock. It is is a field test which contains 30m3 of waste rock with a density of 4.2 g / cm3 (the height of the cell is 20.5m), for a porisity of 0.35. The pyrite (density 4.27) is 1.6% of the Mass fraction  and labratorite (density 2.69) 16.1%. It is open to oxygen. We recovered the effleunt at the 13th day with 0.0027 mole /l of SO4-2 and 0.0016 mole/l of Ca. I am trying to calibrate my model according to this information. I do not know if I miscalculate the percentage of mole and A0 / V at the start. I am confused because I spend whole days there without results and I turn in circles

[dlparkhurst]

OK, you are right that O2 and CO2 should be included in each cell.

Here is your log rate expression

Code: [Select]

30 lograte = 24.65 * parm(1) * parm(3)+lm("O2") * parm(4)+lm("H+") * parm(2)+log10(m/m0)

And here are your parameters

Code: [Select]

-parms    -5.65 0.67 0.5 -0.11

Calculate the rate; I think there is an error somewhere.

When in doubt, simplify and try alternatives. Scale back to a single batch reaction rather than TRANSPORT. Try rate equation for pyrite from phreeqc.dat. Ultimately, you will probably adjust the surface area to get a rate that you like.