Sulfate Reduction Biomass

[stfmtu]

Hello, I am trying to generate biomass through sulfate reduction and hydrogen consumption, as described in exercise 15. The solution already contains S(6) and the system is supposed to take hydrogen gas from the gas phase. However, it seems that the rate is not consuming the sulfate and hydrogen, even though I have inserted the desired formula in the kinetics.

Code: [Select]

SO_4^(2-) (aq)+5H_2 (aq)→H_2 S(aq)+H_2 O(l)Do I need to define a solution species for this, or am I doing something wrong? Should I ad a Equilibrium_Phases?

Thank you!

Code: [Select]

RATES
    Sulfate_formation
-start
10 qm = 1.21E-5
15 Y = 0.07
20 rate = -qm/Y * TOT("S(-2)")/(2.7515E-3 + TOT("S(-2)")) * kin("Biomass") * TOT("Hdg")/(2.9E-6 + TOT("Hdg"))
30 moles = rate * TIME
40 PUT(rate,1)            #save the rate for use in Biomass rate calculation
50 SAVE moles
-end

    Biomass
-start
   1 Y = 0.07          #SRB: Y=0.08, b=8.565e-07
   2 b = 1.21E-6
   10 rate = get(1)
   20 rate = -Y*rate - b*M  # negative rate increases M (biomass)
                              # positive rate (b*B) decreases M (biomass)
   30 moles = -rate * time
   40 if(M + moles) < 0 then moles = -M
   50 save moles
-end
End
 
SOLUTION 1
    temp      48          #optimal 48degC, max 113degC, min 10degC
    pH        7.01 charge   #optimal 7.01, max 11.5, min 1
    pe        4
    redox     pe
    units     umol/L
    S(6)       1000
    -water    1 # kg
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 200
    -volume 2
    -temperature 48      #optimal 48degC, max 113degC, min 10degC
    #CO2(g)      0
    H2S(g)      0
    Hdg(g)      200    #In Phreeqc Database is Hdg
   
END

INCREMENTAL_REACTIONS true

KINETICS 1
   
Sulfate_formation
    -formula  Hdg -5 SO4-2 -1 H2O 1 H2S 1
    -m0       1.409e+01        #all available moles of Hdg --> see phr.out file in the beginning when reaction starts
    -tol      1e-08
   
Biomass
   -formula       CH1.8O0.5N0.2 0.0
   -m              1.0e-5


   
-steps       17280000 in 100 steps  #20 Tage
#-steps       3456000  in 100 steps  #400 Tage
-cvode
#-step_divide 1
#-runge_kutta 3
#-bad_step_max 500
END

USE solution 1
USE kinetics 1
USE gas_phase 1

USER_GRAPH 1
-chart_title "Gaseous components during sulfate reduction"
    -headings               H2S(g) Hdg(g)
    -axis_titles            "Days"  "Pressure gas phase [atm]"
    -axis_scale y_axis      auto auto #log
    #-axis_scale x_axis      0 400
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, GAS("H2S(g)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, GAS("Hdg(g)"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
  -end
    -active                 true
USER_GRAPH 2
-chart_title      "Aqueous species during sulfate reduction"
    -headings               SO4(aq) H2S(aq) Hdg(aq)  Biomass(aq)
    -axis_titles            "Days" "Mole concentration of aqu. species" "Moles of Biomass in aqu. phase"
    -axis_scale y_axis      auto auto #log
    -axis_scale sy_axis      auto auto #log
    #-axis_scale x_axis      0 400
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOT("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("S(-2)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOT("Hdg"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
40 PLOT_XY TOTAL_TIME/86400, KIN("Biomass"), color=Green,symbol=Circle,symbol_size=6,y_axis=2,line_width=1
  -end
    -active                 true

USER_GRAPH 3
-chart_title "Total Moles in the System"
    -headings               SO4 H2S Hdg
    -axis_titles            "Days" "Moles" "Moles of H2S"
    -axis_scale y_axis      auto auto #log
    -axis_scale sy_axis      auto auto #log
    #-axis_scale x_axis      0 400
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(6)")+GAS("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("S(-2)")+GAS("H2S(g)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=2,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOT("Hdg")+GAS("Hdg(g)"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
  -end
    -active                 true

[dlparkhurst]

Here is a script you can work from. You will have to look carefully as to how it differs from yours. The main changes are
(1) in the signs of the reaction rates and the KINETICS formula, (2) the kinetic formula should simply remove Hdg and introduce H2, and (3) the inhibition factor for sulfate reduction should use S(6), not S(-2).

As for the stair step in Hdg, I am going to have to check with Tony Appelo. It looks like there is a tolerance somewhere, such that the gases are only recalculated if they exceed a tolerance.

Code: [Select]

RATES
    Sulfate_reduction
-start
10 qm = 1.21E-5
15 Y = 0.07
#20 rate = -qm/Y * TOT("S(-2)")/(2.7515E-3 + TOT("S(-2)")) * kin("Biomass") * TOT("Hdg")/(2.9E-6 + TOT("Hdg"))
20 rate = qm/Y * TOT("S(6)")/(2.7515E-3 + TOT("S(6)")) * KIN("Biomass") * TOT("Hdg")/(2.9E-6 + TOT("Hdg"))
30 moles = rate * TIME
40 PUT(rate,1)            #save the rate for use in Biomass rate calculation
50 SAVE moles
60 PRINT moles
-end

    Biomass
-start
   1 Y = 0.07          #SRB: Y=0.08, b=8.565e-07
   2 b = 1.21E-6
   10 rate = get(1)
   20 rate = -Y*rate + b*M  # negative rate increases M (biomass)
                              # positive rate (b*B) decreases M (biomass)
   30 moles = rate * time
   #40 if(M + moles) < 0 then moles = -M
   50 save moles
-end
End
SOLUTION 1
    temp      48          #optimal 48degC, max 113degC, min 10degC
    pH        7.01 #charge   #optimal 7.01, max 11.5, min 1
    pe        4
    redox     pe
    units     umol/L
    Na         1  charge
    S(6)       1000
    -water    1 # kg
END
GAS_PHASE 1
    -fixed_pressure
    -pressure 200
    -volume 2
    -temperature 48      #optimal 48degC, max 113degC, min 10degC
    #CO2(g)      0
    H2S(g)      0
    Hdg(g)      200    #In Phreeqc Database is Hdg
END
INCREMENTAL_REACTIONS true
KINETICS 1
Sulfate_reduction
    -formula Hdg -5 H2 2.5
    -m0       1.409e+01        #all available moles of Hdg --> see phr.out file in the beginning when reaction starts
Biomass
   -formula       CH1.8O0.5N0.2 0.0
   -m              1.0e-5
-steps     8640000 in 100 steps  #100 Tage
-cvode
END

USE solution 1
USE kinetics 1
USE gas_phase 1
USER_GRAPH 1
-chart_title "Gaseous components during sulfate reduction"
    -headings               H2S(g) Hdg(g)
    -axis_titles            "Days"  "Pressure gas phase [atm]"
    -axis_scale y_axis      auto auto auto auto log
    #-axis_scale x_axis      0 400
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, GAS("H2S(g)"), color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, GAS("Hdg(g)"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
  -end
    -active                 true
USER_GRAPH 2
-chart_title      "Aqueous species during sulfate reduction"
    -headings               SO4(aq) H2S(aq) Hdg(aq) Biomass(aq)
    -axis_titles            "Days" "Molality" "Moles of Biomass in aqu. phase"
    -axis_scale y_axis      auto auto auto auto log
    -axis_scale sy_axis      auto auto #log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOT("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("S(-2)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOT("Hdg"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
40 PLOT_XY TOTAL_TIME/86400, KIN("Biomass"), color=Green,symbol=Circle,symbol_size=6,y_axis=2,line_width=1
  -end
    -active                 true

USER_GRAPH 3
-chart_title "Total Moles in the System"
    -headings               S(6) S(-2) Hdg
    -axis_titles            "Days" "Molesf" "Moles of Hdg"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(-2)") + GAS("H2S(g)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOTMOLE("Hdg") + GAS("Hdg(g)"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=2,line_width=3
  -end
    -active                 true


[stfmtu]

Thank you, Mr. Parkhurst for your help.

The file sent is really helpful. I would greatly appreciate any advice you can give regarding the stair step in HDG.

My next task will involve iron reduction biomass, also generated through hydrogen consumption. Should I follow a similar approach or would it be better to use Example 9 as a starting point?

[stfmtu]

Hello,

Thank you for the file you shared. I have edited it to fit my needs. I have a question regarding the KINETICS formula. Why is Hdg 5 and H2 2.5?
It seems that H2S would be produced in the gas phase only if I add H2S 1 in the kinetics formula. Would this be correct, or is there another way?

Thank you!

Code: [Select]

RATES
    Sulfate_reduction
-start
10 qm = 1.21E-5
15 Y = 0.07
20 rate = qm/Y * TOT("S(6)")/(2.7515E-3 + TOT("S(6)")) * TOT("Hdg")/(2.9E-6 + TOT("Hdg")) * KIN("Biomass")
30 moles = rate * TIME
40 PUT(rate,1)           
50 SAVE moles
#60 PRINT moles
-end

    Biomass
-start
   1 Y = 0.07         
   2 b = 1.21E-7
   10 rate = get(1)
   20 rate = -Y*rate + b*M  # negative rate increases M (biomass)
                              # positive rate (b*B) decreases M (biomass)
   30 moles = rate * time
   40 save moles
-end
End

SOLUTION 1
    temp      25          #optimal 48degC, max 113degC, min 10degC
    pH        7.01        #optimal 7.01, max 11.5, min 1
    pe        4
    redox     pe
    units     mmol/L
    S(6)       1000
    -water    2 # kg
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 300
    -volume 5            #Gas Phase must be bigger than Solution Volume else H2S goes into Solution
    -temperature 25      #optimal 48degC, max 113degC, min 10degC
    H2S(g)      0
    Hdg(g)      300   
END


INCREMENTAL_REACTIONS true

KINETICS 1
Sulfate_reduction
    -formula      Hdg -1  H2 1  #H2S 1
    -m0           1.409e+10     
   
Biomass
   -formula       CH1.8O0.5N0.2 0.0
   -m              1.0e-5
   
-steps     3456000 in 100 steps  #400 Tage
-cvode
END

USE solution 1
USE kinetics 1
USE gas_phase 1

USER_GRAPH 1
-chart_title      "Gaseous components during Sulfate reduction"
    -headings                H2S(g)  Hdg(g)
    -axis_titles            "Days" "Pressure gas phase [atm]"
    -axis_scale y_axis       auto auto
   # -axis_scale x_axis       0 400
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, GAS("H2S(g)"), color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, GAS("Hdg(g)"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
  -end
    -active                 true
   
USER_GRAPH 2
-chart_title      "Aqueous species during Sulfate reduction"
    -headings                SO4(aq)  H2S(aq)  Hdg(aq)  Biomass(aq)
    -axis_titles            "Days" "Molality" "Moles of Biomass in aqu. phase"
    -axis_scale y_axis       auto auto
    -axis_scale sy_axis      auto auto
    #-axis_scale x_axis       0 400
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(-2)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOTMOLE("Hdg"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
40 PLOT_XY TOTAL_TIME/86400, KIN("Biomass"), color=Green,symbol=Circle,symbol_size=6,y_axis=2,line_width=1
  -end
    -active                 true

USER_GRAPH 3
-chart_title      "Total Moles in the System"
    -headings                SO4  H2S  Hdg
    -axis_titles            "Days" "Moles" "Moles of Hdg"
    -axis_scale y_axis       auto auto
    -axis_scale sy_axis      auto auto
    #-axis_scale x_axis       0 400
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(-2)") + GAS("H2S(g)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOTMOLE("Hdg") + GAS("Hdg(g)"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=2,line_width=3
  -end
    -active                 true

[dlparkhurst]

Actually, my mistake, Hdg and H2 should have the same coefficient; Hdg -1 H2 +1 as you have is correct.

Addition of H2 to the solution will reduce the most favorable (thermodynamically) electron acceptor. So, if O2(aq) is present, it will be reduced by the reaction 2H2+ O2 = 2H2O. Alternatively, NO3, NO2, Fe(3), or SO4 will be reduced. So if none of the more favorable electron acceptors are available, H2 will reduce SO4 and produce S(-2), for example, 4H2 + SO4-2 = S-2 + 4H2O. (Should be HS- in most cases with production of HCO3- or CO3-2, but I'm too lazy to balance that reaction.)

[stfmtu]

Thank you for your response, everything is working properly now.
I have a question regarding the time steps. My calculations finish after around 12 days, but I would like to keep the same volume of solution and gas in the system while extending the simulation to run for 10-20 months exactly with those parameters. Should I introduce TRANSPORT in the code? I have tried, but as I am a beginner, I still haven?t managed to succeed. What should I do to ensure my simulation runs for a couple of months instead of just days? Is that even possible?
Thank you!

Code: [Select]

RATES
    Sulfate_reduction
-start
10 qm = 1.042E-5                                #Specific maximum growth Rate
15 Y = 0.08                                     #Yield coefficient
20 rate = qm/Y * TOT("Hdg")/(5.018E-8 + TOT("Hdg")) * TOT("S(6)")/(1.800E-6 + TOT("S(6)")) * KIN("Biomass")  #Sulfate reduction using dual Monod Rate
#25 PRINT "Sulfate_reduction rate: ", rate
30 moles = rate * TIME
#35 PRINT "Sulfate_reduction moles: ", moles
40 PUT(rate,1)           
50 SAVE moles
-end

    Biomass
-start
   1 Y = 0.08                #Yield coefficient
   2 b = 1.042E-7            #Decay rate
   10 rate = get(1)
   #15 PRINT "Biomass rate_M: ", rate_M
   20 rate = -Y*rate + b*M    # negative rate increases Biomass
                              # positive rate (b) decreases Biomass
   #25 PRINT "Biomass rate: ", rate
   30 moles = rate * time
   #35 PRINT "Biomass moles: ", moles
   40 save moles
-end
End

SOLUTION 1
    temp      30               #optimal 30degC, max 113degC, min 10degC
    pH        7 charge          #optimal 7, max 11.5, min 0.8
    pe        4
    redox     pe
   #Fe_tri     2000
    S(6)       2000
    -water    3                 #Solution volume
END

GAS_PHASE 1
    -fixed_volume
    -pressure    300
    -volume      5              #Gas phase volume
    -temperature 30             #optimal 30degC, max 113degC, min 10degC
    Hdg(g)      200            #Inert Hydrogen
    CO2(g)      100
    H2S(g)      0
END


INCREMENTAL_REACTIONS true

KINETICS 1
Sulfate_reduction
    -formula      Hdg -5  H2 5          #Conversion of inert Hydrogen into reactive Hydrogen and conversion of HS- to H2S
    -m0           4.090e+01             #Initial available moles of Hdg
   
Biomass
    -formula       CH1.8O0.5N0.2 0.0      #Biomass generalized formula
    -m0              1e-04                 #Initial mole concentration of Biomass
   
-steps       1728000 in 100 steps         
#-steps       5184000 in 100 steps
-cvode
END

USE solution 1
USE kinetics 1
USE gas_phase 1

USER_GRAPH 1
-chart_title      "Gaseous components during Sulfate reduction"
    -headings                H2S(g)  Hdg(g)
    -axis_titles            "Time" "Mole concentration"
    -axis_scale y_axis       auto auto
    -axis_scale x_axis       0 20
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, GAS("H2S(g)")/5, color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, GAS("Hdg(g)")/5, color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
  -end
    -active                 true
   
USER_GRAPH 2
-chart_title      "Aqueous species during Sulfate reduction"
    -headings                SO4(aq)  H2S(aq)  Hdg(aq)  Biomass(aq)
    -axis_titles            "Time" "Mole concentration of aqu. species"
    -axis_scale y_axis       auto auto
    -axis_scale sy_axis      auto auto
    -axis_scale x_axis       0 20
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(-2)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOTMOLE("Hdg"), color=Blue,symbol=Diamond,symbol_size=0,y_axis=1,line_width=3
  -end
    -active                 true

USER_GRAPH 3
-chart_title      "Total Moles in the System"
    -headings                SO4 H2S Hdg
    -axis_titles            "Time" "Moles"
    -axis_scale y_axis       auto 9
    -axis_scale sy_axis      0 auto
    -axis_scale x_axis       0 20
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(6)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOTMOLE("S(-2)") + GAS("H2S(g)")/5,color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 PLOT_XY TOTAL_TIME/86400, TOTMOLE("Hdg") + GAS("Hdg(g)")/5, color=Blue,symbol=Diamond,symbol_size=0,y_axis=2,line_width=3
  -end
    -active                 true

[dlparkhurst]

There are a few ways to extend the simulation.

(1) You can simply add more time in the KINETICS definition, and the simulation will continue until the specified time.

(2) You can add SAVE solution 1; SAVE gas_phase 1; END at the end of the simulation. This will save the state of the solution and the gas phase (KINETiCS is saved automatically).  You can then define a new reaction like USE solution 1; USE kinetics 1; USE gas_phase 1, which would repeat the simulation for another 100 steps.

(3) You can also look at RUN_CELLS, which does USE and SAVE automatically.