Methanogenic bacteria and Methane formation

[stfmtu]

Dear all,

I am trying to generate biomass that should be linked with methane formation, but the biomass doesn't appear in the plot. When I attempt to change the number of days from 10 to 100, it seems that after the 10th day, I only get constant values for all the plotted data. Could someone explain what I am doing wrong?
As my model was partially adapted from this forum, I would also like to ask why, in the kinetics of methane formation, only hydrogen is mentioned and there's no mention of other elements such as oxygen, water, or methane itself.

SOLUTION_MASTER_SPECIES
   Biomass       Biomass          0     Biomass         1
SOLUTION_SPECIES

Biomass = Biomass
   -gamma 3
   log_k 0
   

RATES
    Methane_formation
-start
10 if (m < 0) then goto 70
20 rate = 2.3e-09*(TOT("C(4)")/1e-3 + TOT("C(4)"))
30 moles = rate * TIME
40 if (moles > m) then moles = m
50 SAVE moles
60 put(rate,1)
70 END
-end

    Biomass
-start
   1 Y = parm(1)
   2 d = parm(2)
   10 rate_M = get(1)
   20 B = kin("Biomass")
   30 rate = -Y*rate_M - d*B
   40 dB = rate * time
   50 save -dB
-end
End
 
SOLUTION 1
    temp      37
    pH        7 charge
    pe        4
    redox     pe
    units     ppm
    density   1
    -water    1 # kg
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 37
    CH4(g)    0
    CO2(g)    40
    Hdg(g)    160 #In Phreeqc Database is Hdg
   
END

INCREMENTAL_REACTIONS true

KINETICS 1
   
Methane_formation
    -formula  H2  1 Hdg  -1
    -m0       1
    -tol      1e-08
   
Biomass
   -formula       H 0.0
   -m              1.36e-4
   -parms          0.89 5.78E-7   
   
   
-steps       864000 in 100 steps # seconds; 100 Tage
-step_divide 1
-runge_kutta 3
-bad_step_max 500
END

USE solution 1
USE kinetics 1
USE gas_phase 1

USER_GRAPH 1
    -headings               time TDIC CH4(aq) Hdg(aq)  Biomass CO2(g) Hdg(g) CH4(g)
    -axis_titles            "Days" "Mole concentration" "Pressure gas phase [atm]"
    -axis_scale y_axis      1e-6 2 auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 86400
30 GRAPH_Y TOT("C(4)"), TOT("C(-4)"), TOT("Hdg"), TOT("Biomass")
40 GRAPH_SY GAS("CO2(g)"), GAS("Hdg(g)"), GAS("CH4(g)")
  -end
    -active                 true
   
 Thank you!

[dlparkhurst]

The -formula for biomass is H 0, which is fine. There will be no change in the solution composition from this reaction. However, it implies that the moles of kinetic reactant Biomass is the quantity representing the biomass. As such, a positive value for the SAVE quantity will decrease the biomass, and a negative value for the SAVE quantity will

Code: [Select]

increase biomass. Therefore the growth rate should be negative and the death rate should be positive. Sorry, maybe it should be the other way, but that is the way it is; it makes more sense if the reaction were dissolving a mineral, where concentrations in solution increase with positive rates.

By converting Hdg to H2, you are adding reductant to the solution. Electron acceptors will be used in the sequence determined by thermodynamics. This is generally the way you expect redox reactions to proceed. In simple terms, O2 will be consumed then NO3- then SO4-2, and finally C(4) will be reduced. If you don't have any other electron acceptors, then you will only see methanogenesis.

Code: [Select]

RATES
    Methane_formation
-start
10 if (m < 0) then goto 70
20 rate = 2.3e-09*(TOT("C(4)")/1e-3 + TOT("C(4)"))
30 moles = rate * TIME
40 if (moles > m) then moles = m
50 SAVE moles
60 put(rate,1)
70 END
-end

    Biomass
-start
   1 Y = parm(1)
   2 d = parm(2)
   10 rate_M = get(1)
   20 B = kin("Biomass")
   30 rate = -Y*rate_M + d*B  # negative rate increases M (biomass)
                              # positive rate (d*B) decreases M (biomass)
   40 dB = rate * time
   50 save dB
-end
End
 
SOLUTION 1
    temp      37
    pH        7 charge
    pe        4
    redox     pe
    units     ppm
    density   1
    -water    1 # kg
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 37
    CH4(g)    0
    CO2(g)    40
    Hdg(g)    160 #In Phreeqc Database is Hdg
   
END

INCREMENTAL_REACTIONS true

KINETICS 1
   
Methane_formation
    -formula  H2  1 Hdg  -1
    -m0       1
    -tol      1e-08
   
Biomass
   -formula       H 0.0
   -m              1.36e-4
   -parms          0.89 5.78E-7   
   
   
-steps       864000 in 100 steps # seconds; 100 Tage
-step_divide 1
-runge_kutta 3
-bad_step_max 500
END

USE solution 1
USE kinetics 1
USE gas_phase 1

USER_GRAPH 1
    -headings               time TDIC CH4(aq) Hdg(aq)  Biomass CO2(g) Hdg(g) CH4(g)
    -axis_titles            "Days" "Mole concentration" "Pressure gas phase [atm]"
    -axis_scale y_axis      1e-6 2 auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 86400
30 GRAPH_Y TOT("C(4)"), TOT("C(-4)"), TOT("Hdg"), KIN("Biomass")
40 GRAPH_SY GAS("CO2(g)"), GAS("Hdg(g)"), GAS("CH4(g)")
  -end
    -active                 true


[stfmtu]

The following code is based on the recommendations given. Would you be so kind as to take a look at why the HDG has a negative value and results in an error?

RATES
    Methane_formation
-start
10 if (m < 0) then goto 70
20 rate = 1.643e-5*TOT("C(4)")/(1.094e-7 + TOT("C(4)"))
30 moles = rate * TIME
40 if (moles > m) then moles = m
50 SAVE moles
60 put(rate,1)
70 END
-end

    Biomass
-start
   1 Y = parm(1)
   2 d = parm(2)
   10 rate_M = get(1)
   20 B = kin("Biomass")
   30 rate = -Y*rate_M + d*B  # negative rate increases M (biomass)
                              # positive rate (d*B) decreases M (biomass)
   40 dB = rate * time
   50 save dB
-end
End
 
SOLUTION 1
    temp      45          #optimal 45degC, max 122degC, min 10degC
    pH        7.7 charge   #optimal 7.7, max 10.2, min 4.1
    pe        4
    redox     pe
    units     ppm
    density   1
    -water    1 # kg
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 45      #optimal 45degC, max 122degC, min 10degC
    CH4(g)    0
    CO2(g)    40
    Hdg(g)    160 #In Phreeqc Database is Hdg
   
END

INCREMENTAL_REACTIONS true

KINETICS 1
   
Methane_formation
    -formula  H2  1 Hdg  -1
    -m0       6.118        #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.36e-4
   -parms          0.03 8.565e-7   # Methanogeic bacteria: Y=0.03, d=8.565e-07
   
   
#-steps       8640000 in 100 steps # seconds; 100 Tage
#-steps       86400000 in 100 steps # seconds; 1000 Tage
#-steps       1728000 in 100 steps # seconds; 20 Tage
-steps       34560000 in 100 steps # seconds; 400 Tage
-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 methanogenesis"
    -headings               CO2(g) CH4(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("CO2(g)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, GAS("CH4(g)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 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 methanogenesis"
    -headings               CO2(aq) CH4(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("C(4)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("C(-4)"),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               CO2 CH4 Hdg
    -axis_titles            "Days" "Moles" "Moles"
    -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("C(4)")+GAS("CO2(g)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("C(-4)")+GAS("CH4(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]

The sign of moles in SAVE moles in the methane_formation rate is positive. The coefficient of Hdg in KINETICS is negative. The combination is trying to remove Hdg from solution, but there is no Hdg to be removed. The rate calculation cannot procede because no matter what time step is used, it will make a negative concentration of Hdg.

[stfmtu]

Hello,

I have tried changing the sign of moles or the one for Hdg, but neither way works. I have also tried modifying the KINETICS formula to: H2 -4 CO2 -1 CH4 1 H2O 2, but still no results. My code should consume the available hydrogen and produce methane and a biomass that will decay the moment hydrogen or carbon dioxide is not available.

Any help on identifying my mistakes would be greatly appreciated.

[dlparkhurst]

Sorry, I was wrong in my previous comment. There is a source of Hdg in the GAS_PHASE. So, the reaction can proceed. However, it only proceeds until the point that Hdg is completely consumed, at which time, the reaction can go no farther without generating negative Hdg concentrations.

I added an inhibitor factor for Hdg in the methane rate equation. Also, the Hdg is consumed in about 4 days, after which the biomass decays for the remaining time of the simulation. I ran for 40 days.

You may need to consider your GAS_PHASE definition. You have defined 1 L of volume at 45 C. PHREEQC calculates the amount of each gas component roughly with the ideal gas law (although Peng-Robinson for phreeqc.dat), for example P(Hdg)*1 = 160 = n*R*318. When the reaction occurs the gases dissolve and reequilibrate with a gas phase with volume such that the volume is adjusted until the sum of the partial pressures of the gas components is 100 atm.As the Hdg is consumed, the gas phase shrinks. Alternatives include using a larger volume of gas by which the moles of Hdg and CO2 in the initial gas are calculated, or specifying a partial pressure of Hdg (and or CO2) with EQUILIBRIUM_PHASES.

Code: [Select]

RATES
    Methane_formation
-start
10 if (m < 0) then goto 70
20 rate = 1.643e-5*TOT("C(4)")/(1.094e-7 + TOT("C(4)"))*TOT("Hdg")/(1e-10 + TOT("Hdg"))
30 moles = rate * TIME
40 if (moles > m) then moles = m
50 SAVE moles
60 put(rate,1)
70 END
-end

    Biomass
-start
   1 Y = parm(1)
   2 d = parm(2)
   10 rate_M = get(1)
   20 B = kin("Biomass")
   30 rate = -Y*rate_M + d*B  # negative rate increases M (biomass)
                              # positive rate (d*B) decreases M (biomass)
   40 dB = rate * time
   50 save dB
-end
End
 
SOLUTION 1
    temp      45          #optimal 45degC, max 122degC, min 10degC
    pH        7.7 charge   #optimal 7.7, max 10.2, min 4.1
    pe        4
    redox     pe
    units     ppm
    density   1
    -water    1 # kg
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 45      #optimal 45degC, max 122degC, min 10degC
    CH4(g)    0
    CO2(g)    40
    Hdg(g)    160 #In Phreeqc Database is Hdg
   
END

INCREMENTAL_REACTIONS true

KINETICS 1
   
Methane_formation
    -formula  H2  1 Hdg  -1
    -m0       6.118        #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.36e-4
   -parms          0.03 8.565e-7   # Methanogeic bacteria: Y=0.03, d=8.565e-07
   
   
#-steps       8640000 in 100 steps # seconds; 100 Tage
#-steps       86400000 in 100 steps # seconds; 1000 Tage
#-steps       1728000 in 100 steps # seconds; 20 Tage
-steps       3456000 in 100 #34560000 in 100 steps # seconds; 400 Tage
-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 methanogenesis"
    -headings               CO2(g) CH4(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("CO2(g)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, GAS("CH4(g)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 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 methanogenesis"
    -headings               CO2(aq) CH4(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("C(4)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("C(-4)"),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               CO2 CH4 Hdg
    -axis_titles            "Days" "Moles" "Moles"
    -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("C(4)")+GAS("CO2(g)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("C(-4)")+GAS("CH4(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


[stfmtu]

Dear dlparkhurst,

Thank you for the advice. It seems that when I change the volume of the gas phase, the simulation runs until the desired number of days. However, the only issue is that the biomass does not decrease once the hdg and CO2 are consumed. It stagnates at the same value. What can I do in this situation?

I have also tried to add Equilibrium_Phases, but it seems it did not do anything for me. I am most likely doing something wrong. Could you possibly show the right way as reference for my code to consume CO2 and Hdg (or H2) and produce CH4?

Thank you.

[dlparkhurst]

Biomass decreases in the code I provided. You must have done something differently.  You will have to provide your script, or better yet, look at the output file and figure out what you have done.

[stfmtu]

Dear Mr. David,

Thank you for your advice. Yes, I have readjusted the data, and now it is working properly. I also managed to achieve good results without the Hdg inhibitor. However, the issue is that even though the biomass has a decreasing rate, it does not decrease once the hydrogen and CO2 are no longer consumed. Could you please have a look?

Thank you!

RATES
    Methane_formation
-start
10 if (m < 0) then goto 70
20 rate = 1.643e-5*TOT("C(4)")/(1.094e-7 + TOT("C(4)"))
30 moles = rate * TIME
40 if (moles > m) then moles = m
50 SAVE moles
60 put(rate,1)
70 END
-end

    Biomass
-start
   1 Y = parm(1)
   2 d = parm(2)
   10 rate_M = get(1)
   20 B = kin("Biomass")
   30 rate = -Y*rate_M + d*B  # negative rate increases M (biomass)
                              # positive rate (d*B) decreases M (biomass)
   40 dB = rate * time
   50 save dB
-end
End
 
SOLUTION 1
    temp      45          #optimal 45degC, max 122degC, min 10degC
    pH        7.7 charge   #optimal 7.7, max 10.2, min 4.1
    pe        4
    redox     pe
   # units     ppm
    density   1
    -water    1 # kg
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 250
    -volume 8.5
    -temperature 45      #optimal 45degC, max 122degC, min 10degC
    CH4(g)    0
    CO2(g)    62.5
    Hdg(g)    187.5 #In Phreeqc Database is Hdg
   
END

INCREMENTAL_REACTIONS true

KINETICS 1
   
Methane_formation
    -formula  H2  1 Hdg  -1
    -m0       6.000e+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.36e-4
   -parms          0.03 8.565e-7   # Methanogeic bacteria: Y=0.03, d=8.565e-07
   
   
#-steps       17280000 in 100 steps # seconds; 200 Tage
#-steps       86400000 in 100 steps # seconds; 1000 Tage
#-steps       1728000 in 100 steps # seconds; 20 Tage
-steps       4320000 in 100 steps # seconds; 50 Tage
-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 methanogenesis"
    -headings               CO2(g) CH4(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("CO2(g)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, GAS("CH4(g)"),color=Orange,symbol=Triangle,symbol_size=0,y_axis=1,line_width=3
30 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 methanogenesis"
    -headings               CO2(aq) CH4(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("C(4)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("C(-4)"),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               CO2 CH4 Hdg
    -axis_titles            "Days" "Moles" "Moles"
    -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("C(4)")+GAS("CO2(g)"), color=Black,symbol=Square,symbol_size=0,y_axis=1,line_width=3
20 PLOT_XY TOTAL_TIME/86400, TOT("C(-4)")+GAS("CH4(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]

You can use PRINT statements in the RATES definitions to debug your rate expressions.

In this case, you are using PUT(rate), where rate is positive, even though the SAVEd value is zero. In the biomass rate, the rates reach steady state where growth (based on the PUT value) is balanced by death.