Coupling organic C degradation and redox processes in transport simulation

[Lukem]

Dear David and colleagues,

I am trying to simulate sea level rise effects in an initially oxidised soil column via using the TRANSPORT function. The main goal is to simulate reduction processes in the soil (i.e. reductive dissolution of nitrate, Fe and Mn oxides) as sea level rises, potentially enhanced by  the entry of terminal electron acceptors via the seawater (i.e. sulfate in particular). I have FeS included as a sulfide mineral phase that could potentially form (i.e. 0 concentration at start). To drive the reduction reactions in the model I am also trying to include microbial degradation of organic C to supply DOC/organic electron donors (via the "Organic_C" Monod rate expression included in the Phreeqc.dat database).

My code is below which runs but has some issues that I am having trouble solving. The main issue at present seems to be with coupling the KINETIC expression for microbial degradation of C with the TRANSPORT part, as the organic C input only seems to be being implemented at the end of my 30 year simulation (USER GRAPH 3 illustrates this), when I want this input occurring right from the start of the simulation. Would really appreciate your advice on this.

Code: [Select]

#sea level rise scenario
SOLUTION_SPECIES  #dummy species to help with redox stability
    H2O + 0.01e- = H2O-0.01
    -log_k -9

SOLUTION 1-30 #Oxidised fresher pore water (1/10 seawater)
        units   ppm
        pH      7.22
        pe      8.451
        density 1.023
        temp    25.0
        Ca              41.23
        Mg              129.18
        Na              1076.80
        K               39.91
        Fe              0.002
        Mn              0.0002
        Cl              1935.30 charge
        Alkalinity      141.682 as HCO3
        S(6)            271.20
        N(5)            0.32    gfw   62.0 #nitrate
        O(0)            1.0     O2(g) -0.7

EQUILIBRIUM_PHASES 1-30 #solid phases
 Calcite;  CO2(g) -3.4 #Equilibrium with atmospheric CO2
 Goethite  0 2.5e-3  #Fe oxide
 FeS(ppt)  0 0       #Fe sulfide
 Pyrolusite  0 4e-5  #Mn oxide

KINETICS 1-30 # Organic matter degradation, the RATE expression for Organic_C is in Phreeqc.dat
Organic_C
-formula C
-tol    1e-8
-m      5e-3   # SOC in mol
#-time 30 year in 30  #Turned off as time is specified in TRANSPORT block??
INCREMENTAL_REACTIONS true

SAVE solution 1-30
PRINT; -status true
END                                           

SOLUTION 0   #seawater composition entering the soil column
        units   ppm
        pH      8.22
        pe      8.451
        density 1.023
        temp    25.0
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Fe              0.002
        Mn              0.0002
        Cl              19353.0 charge
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
        N(-3)           0.32    as    NH4

TRANSPORT    # 30 yr flow, 3 cm/yr, 90cm
 -cells 30;               -length 0.03
 -time_step 3.15e7
 -flow_direction forward; -shifts 30
 -dispersivity 0.1;         -punch_frequency 30

SELECTED_OUTPUT
        -file            Column_model_1.csv
        -reset           false
        -solution
        -distance        true
        -totals          Cl Ca S(6) C Fe(2) S(-2)
        -equilibrium_phases  goethite FeS(ppt) pyrolusite

USER_GRAPH 1
 -heading dist Fe(2) S(-2) S(6) FeS
 -init false                                   # -init is shorthand for initial_solutions
 -plot_concentration_vs x
 -axis_titles "Distance / m" "mmol/L" "FeS, moles"
 -start
 10 graph_x dist
 20 graph_y tot("Fe(2)"), tot("S(-2)"), tot("S(6)")
 30 graph_sy equi("FeS(ppt)")             
 -end

USER_GRAPH 2
 -headings time Fe(2) S(-2) S(6) FeS
 -axis_titles "Time / years" "mmoL"
 -start
 10 graph_x total_time/3.1536e7 #converting seconds to years
 20 graph_y tot("Fe(2)"), tot("S(-2)"), tot("S(6)")
 30 graph_sy equi("FeS(ppt)")
 -end

USER_GRAPH 3
 -headings time Organic_C
 -axis_titles "Time / years" "Organic C/ mmoL"
 -start
 10 graph_x total_time/3.1536e7 #converting seconds to years
 20 graph_y KIN("Organic_C")
 -end

END
Kind regards,

Luke

[dlparkhurst]

First, I think one main source of confusion is that the boundary cell, cell 31, is plotted. I don't like to second guess Tony, but the boundary cell should not be plotted by default. If you set -punch_cells 1-30, then you are looking just at the column. The boundary cells are probably useful if constant boundary conditions are used, but not for flux or closed boundary conditions.

There are a lot of moving parts, but the organic decomposition reaction is being used in all cells, and the decomposition reaction consumes electron acceptors. The seawater has ammonia, which also consumes electron acceptors as it advects/diffuses into the column. The seawater also has a little bit of Fe(3) and Mn(2), which may cause some redox reactions.

I have made a few changes to the script to demonstrate the processes. First, just to simplify, I made the system diffusion only. I limited shifts to 8 to demonstrate the redox sequence. Finally, I rearranged the data block so that KINETICS was not used; the fresh water reacts with EQUILIBRIUM_PHASES to define solutions 1-30.

Three graphs are generated by the script below. The first graph shows, starting from the bottom of the column, a regions with dissolved constituents (1) both oxygen and nitrate, (2) nitrate, (3)  manganese, and (4) manganese and iron.

The second graph shows regions where pyrolusite is present and absent, and a decrease in goethite at the surface of the column. The amount of organic carbon remaining decreases monotonically down the column.

The third graph shows essentially rates of reaction for pyrolusite and goethite dissolution and precipitation. Also, the rate of organic decomposition is less in the upper part of the column, where only sulfate is present than deeper in the column where nitrate, or nitrate plus oxygen are present.

At no point is FeS(ppt) generated.

You can add back the effects of advection to TRANSPORT, although diffusion will still be almost as important as advection.

Code: [Select]

#sea level rise scenario
SOLUTION_SPECIES  #dummy species to help with redox stability
    H2O + 0.01e- = H2O-0.01
    -log_k -9
RATES
##########
#Organic_C
##########
#
# Example of KINETICS data block for SOC (sediment organic carbon):
#       KINETICS 1
#       Organic_C
#               -formula C
#               -tol    1e-8
#               -m      5e-3   # SOC in mol
#               -time 30 year in 15
Organic_C
 -start
1   REM      Additive Monod kinetics for SOC (sediment organic carbon)
2   REM      Electron acceptors: O2, NO3, and SO4

10  if (M <= 0) THEN GOTO 200
20  mO2   = MOL("O2")
30  mNO3  = TOT("N(5)")
40  mSO4  = TOT("S(6)")
50  k_O2  = 1.57e-9    # 1/sec
60  k_NO3 = 1.67e-11   # 1/sec
70  k_SO4 = 1.e-13     # 1/sec
80  rate  = k_O2 * mO2/(2.94e-4 + mO2)
90  rate  = rate + k_NO3 * mNO3/(1.55e-4 + mNO3)
100 rate  = rate + k_SO4 * mSO4/(1.e-4 + mSO4)
110 moles = rate * M * (M/M0) * TIME
200 SAVE moles
 -end
END
SOLUTION 1-30 #Oxidised fresher pore water (1/10 seawater)
        units   ppm
        pH      7.22
        pe      8.451
        density 1.023
        temp    25.0
        Ca              41.23
        Mg              129.18
        Na              1076.80
        K               39.91
        Fe              0.002
        Mn              0.0002
        Cl              1935.30 charge
        Alkalinity      141.682 as HCO3
        S(6)            271.20
        N(5)            0.32    gfw   62.0 #nitrate
        O(0)            1.0     O2(g) -0.7

EQUILIBRIUM_PHASES 1-30 #solid phases
 Calcite 0 10 
 CO2(g) -3.4 #Equilibrium with atmospheric CO2
 Goethite  0 2.5e-3  #Fe oxide
 FeS(ppt)  0 0       #Fe sulfide
 Pyrolusite  0 4e-5  #Mn oxide
SAVE solution 1-30
END

KINETICS 1-30 # Organic matter degradation, the RATE expression for Organic_C is in Phreeqc.dat
Organic_C
-formula C
-tol    1e-8
-m      5e-3   # SOC in mol
END                                           

SOLUTION 0   #seawater composition entering the soil column
        units   ppm
        pH      8.22
        pe      8.451
        density 1.023
        temp    25.0
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Fe              0.002
        Mn              0.0002
        Cl              19353.0 charge
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
        N(-3)           0.32    as    NH4
END

TRANSPORT    # 30 yr flow, 3 cm/yr, 90cm
-cells 30;               
-length 0.03
-time_step 3.15e7
#-flow_direction forward
#-bc flux flux
-dispersivity 0.1;
-shifts 8
-punch_frequency 8
-punch_cells 1-30
-flow diffusion_only
-bc constant closed

USER_GRAPH 1
 -heading dist O(0) N(5) Mn Fe
 -init false                                   # -init is shorthand for initial_solutions
 -plot_concentration_vs x
 -axis_titles "Distance, m" "log10 mol/kgw" ""
 -start
 10 graph_x dist
 20 graph_y LOG10[tot("O(0)")], LOG10[tot("N(5)")], LOG10[TOT("Mn")], LOG10[TOT("Fe")]
 -end

USER_GRAPH 2
 -heading dist Pyrolusite Goethite FeS(ppt) Organic_C
 -init false                                   # -init is shorthand for initial_solutions
 -plot_concentration_vs x
 -axis_titles "Distance, m" "moles" "Organic_C, moles"
 -start
 10 graph_x dist
 20 graph_y equi("Pyrolusite"), equi("Goethite"), equi("FeS(ppt)")
 30 graph_sy KIN("Organic_C")
 -end

USER_GRAPH 3
 -heading dist Delta_Pyrolusite Delta_Goethite Delta_FeS(ppt) Delta_organic_C
 -init false                                   # -init is shorthand for initial_solutions
 -plot_concentration_vs x
 -axis_titles "Distance, m" "moles" "Delta Organic_C, moles"
 -start
 10 graph_x dist
 20 graph_y equi_delta("Pyrolusite"), equi_delta("Goethite"), equi_delta("FeS(ppt)")
 30 graph_sy KIN_DELTA("Organic_C")
 -end
END



[Lukem]

Dear David

Thanks so much for your reply and help. I think the transport part is working well conceptually now thanks. In regards to the reduction reactions, I am not sure we are seeing those in the simulation they way I was expecting as I have plotted concentrations of Fe and Mn oxide minerals, Fe2+ and dissolved sulfide and there is little change. I think this is due to insufficient carbon addition in the model, so just oxygen being consumed. Looking at the plot of the organic carbon being consumed, the total amount of C consumed is quite low (approx. 2.5e-4 moles when I was expecting/wanting to supply the full 5e-3 moles of C added in the simulation input (m parameter). I tried to change some of the rate and MONOD parameters to see if that was the issue but still could not seem to get that full consumption of soil organic C, and supply of C to solution to occur. Appreciate any further insights you have.

Below is code I modified from your version:

Code: [Select]

#sea level rise scenario
SOLUTION_SPECIES  #dummy species to help with redox stability
    H2O + 0.01e- = H2O-0.01
    -log_k -9
RATES
Organic_C
 -start
1   REM      Additive Monod kinetics for SOC (sediment organic carbon)
2   REM      Electron acceptors: O2, NO3, and SO4

10  if (M <= 0) THEN GOTO 200
20  mO2   = MOL("O2")
30  mNO3  = TOT("N(5)")
40  mSO4  = TOT("S(6)")
50  k_O2  = 1.57e-9    # 1/sec   
60  k_NO3 = 1.67e-11   # 1/sec 
70  k_SO4 = 1.e-13     # 1/sec
80  rate  = k_O2 * mO2/(2.94e-4 + mO2)
90  rate  = rate + k_NO3 * mNO3/(1.55e-4 + mNO3)
100 rate  = rate + k_SO4 * mSO4/(1e-4 + mSO4)
110 moles = rate * M * (M/M0) * TIME
200 SAVE moles
 -end
END

SOLUTION 1-30 #Oxidised fresher pore water (1/10 seawater)
        units   ppm
        pH      7.22
        temp    25.0
        Ca              41.23
        Mg              129.18
        Na              1076.80
        K               39.91
        Cl              1935.30 charge
        Alkalinity      141.682 as HCO3
        S(6)            271.20
        N(5)            0.02    gfw   62.0 #nitrate
        O(0)            1.0     O2(g) -0.7

EQUILIBRIUM_PHASES 1-30 #solid phases
# Calcite 0 10
 CO2(g) -3.4 #Equilibrium with atmospheric CO2
 Goethite  0 2.5e-3  #Fe oxide
 FeS(ppt)  0 0       #Fe sulfide
 Pyrolusite  0 4e-5  #Mn oxide
SAVE solution 1-30
END

KINETICS 1-30 # Organic matter degradation, the RATE expression for Organic_C is in Phreeqc.dat
Organic_C
-formula C
-tol    1e-8
-m     5e-3   # Starting SOC in mol
END                                           

SOLUTION 0   #seawater composition entering the soil column
        units   ppm
        pH      8.22
        temp    25.0
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Cl              19353.0 charge
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
END

TRANSPORT    # 30 yr flow, 3 cm/yr, 90cm
-cells 30;               
-length 0.03
-time_step 3.15e7
#-flow_direction forward
#-bc flux flux
-dispersivity 0.1;
-shifts 8
-punch_frequency 8
-punch_cells 1-30
-flow diffusion_only
-bc constant closed

USER_GRAPH 1
 -heading dist Na Cl SO4 NO3 O2
 -init false                                   # -init is shorthand for initial_solutions
 -plot_concentration_vs x
 -axis_titles "Distance, m" "log10 mol/kgw" "O2"
 -start
 10 graph_x dist
 20 graph_y LOG10[tot("Na")], LOG10[tot("Cl")], LOG10[tot("S(6)")], LOG10[tot("N(5)")],
 30 graph_sy tot("O(0)")/2
 -end

USER_GRAPH 2
 -heading dist Fe(2) S(-2) S(6) FeS
 -init false                                   # -init is shorthand for initial_solutions
 -plot_concentration_vs x
 -axis_titles "Distance / m" "mmol/L" "FeS, moles"
 -start
 10 graph_x dist
 20 graph_y tot("Fe(2)"), tot("S(-2)"), tot("S(6)")           
 -end

USER_GRAPH 3
 -heading dist Pyrolusite Goethite FeS(ppt) Organic_C
 -init false                                   # -init is shorthand for initial_solutions
 -plot_concentration_vs x
 -axis_titles "Distance, m" "log10 moles" "Organic_C, moles"
 -start
 10 graph_x dist
 20 graph_y LOG10[equi("Pyrolusite")], LOG10[equi("Goethite")], LOG10[equi("FeS(ppt)")]
 30 graph_sy KIN("Organic_C")
 -end

END
Many thanks,

Luke

[dlparkhurst]

In your model, the consumption of S(6) depends on the rate constant (1e-13) and the amount of organic matter (5e-3). If you calculate the maximum rate per year it is

Code: [Select]

1e-13*5e-3*3.15e7 = 1.58e-8 mol/y in a 1 L vol

Not much relative to 0.03 mol/L sulfate. So consider first the amount of organic carbon available per liter of water to see if 5e-3 is appropriate. I think even if it is substantially larger, the above calculation will still be relatively small. That leaves the rate constant.

Many years ago when I worked with estuarine sediment pore waters, the rates of sulfate reduction we found were sufficient to reduce all of the available sulfate within the first few centimeters of sediment; so the rate was sufficient to maintain a diffusion gradient of say 30 mmol/5 cm or greater.

The rate constant depends on how juicy the organic matter. The gradient I gave is for pretty reactive organic matter (large rate constant) that is constantly accumulating in the sediment. In your model, the original water contains oxygen and you have a fairly small value for organic carbon. If the sediments have been exposed to an oxygenated environment, the organic matter may be fairly limited and refractory (small rate constant). With only 5  mmol of organic matter to work with, there will not be much effect on sulfate diffusion, no matter what rate constant you use. Even if it all reacts, diffusion will quickly replace all of the sulfate that does react.

[Lukem]

Thanks very much David for your help, yes I think that was the issue, I was looking at that initial organic carbon parameter the 'wrong way around' based on DOC concentration rather than Soil Organic Carbon (including particulate carbon) that then kinetically degrades and produces DOC. Plotting SOC on log scale also helped, there was a lot of consumption where sulfate was entering the column but was also some SOC consumption further up column. Best wishes, Luke