tracking 13C fractionation in a marl lake

[grundl]

Hello everyone,

I am trying to model the 13C content of the aqueous carbonate species and calcite in a groundwater fed lake that is actively precipitating calcite as the relatively acidic, CO2-rich groundwater enters the lake and degases CO2. This is analogous to the classic example given by Appelo in chapter 27 of Schulz and Hadeler book "Geochemical Processes in Soil and Groundwater".
I used Appelo's input file and changed the input solution from his solution (pure Ca2+ CO32-) to the actual makeup of the incoming groundwater but it does not work. The simulation runs and gives outputs that make sense (pH rises, pCO2 drops, calcite precipitates), but no matter how I change graph axes or the initial isotopic composition of the incoming groundwater, I cannot get a graph showing the 13C compositions.
I am not as adept as I should be with the KINETICS/RATE subroutines so maybe this is the problem. Can anyone help me out? I have attached my input file.

Cheers,

Tim Grundl

[dlparkhurst]

I have to admit that I did not try to debug the Appelo example. The permil values are about 1000, which you can see if you right click on a graph and select Set Scale to Default. Note that the example from the book should be run with phreeqc.dat, not iso.dat. Iso.dat uses the species H3O+, whereas phreeqc.dat uses H+. I don't have access to Tony's book right now, but I also would take a different approach that may or may not be simpler; depends on your point of view. 

My approach would be to use iso.dat. I have attached an example below. I'm not sure I have created the calculation that you want. I have assumed that the water eventually equilibrates with air, which may not be the case. However, the endpoint PCO2 and d13C can be adjusted in the parameters of KINETICS.

As background, iso.dat assumes that carbon-13 (denoted [13C]) is a separate component from carbon-12, so they react independently. Carbon-13 calcite has a slightly different equilibrium constant than carbon-12 calcite, same with Henry's law constants for CO2, and the kinetics of the two isotopes should be slightly different. I won't describe the input file in detail, but the basics are that both CO2 and [13C]O2 are lost from solution at rates that depend on the endpoint PCO2 (same as Appelo) and d13C (differs from Appelo, I think). The rate constant is 0.001 for CO2, but slightly smaller for [13C]O2, by the factor alpha--the fractionation factor. Calcite (carbon-12) precipitates according to the rate equation and Ca[13C]O3 precipitates according to equilibrium fractionation. The "apparent" d13C of gas is the d13C calculated from the partial pressures of the gases in solution; at the end, it is the actual d13C specified for air. I did not know the d13C of the groundwater, so I arbitrarily chose -13.0.

I'm not positive that I have done everything correctly. First, I may not have the right conceptual model of using an endpoint. Second, I'm not sure about the use of alpha in the rate equation; it does make a (relatively small) difference in the results. Still, the results seem plausible, and perhaps it is  close enough that you can modify the script to suit your needs.

Code: [Select]

RATES
CO2(g)
# d(CO2) / dt = -k * (P_aq - P_gas)
# parms 1= log(P_CO2)_gas. 2= reaction factor relative to calcite rate
 -start
  10 P_aq = SR("CO2(g)")
  20 dCO2 = parm(2) * (10^parm(1) - P_aq) * TIME
  30 SAVE dCO2
  50 PUT(parm(1), 15) # log10 PCO2 in atmosphere
  60 PUT(parm(2), 16) # rate constant
  70 PUT(dCO2, 5)
 -end

[13C]O2(g)
# d(13CO2) / dt = -k * (P_aq - P_gas)
# parm(1) is 13C of the atmosphere
-start
10 pCO2_atm = 10^GET(15)                                # PCO2 in atmosphere
20 p13CO2_atm = (parm(1)/1000 + 1)*0.0111802 * pco2_atm # P13CO2 in atmosphere
30 alpha = 10^(LK_NAMED("Log_alpha_13C_CO2(aq)/CO2(g)"))
40 d13CO2 = GET(16) * alpha * (p13CO2_atm - SR("[13C]O2(g)"]) * TIME 
50 SAVE d13CO2
-end


Calcite
# Calcite rate modified from Plummer et al., 1978
# parms 1= A/V, 1/dm. 2= exponent for m/m0.
 -start
  10 si_cc = si("Calcite")
  20 if (m <= 0  and si_cc < 0) then goto 200
  30  k1 = 10^(0.198 - 444.0 / (273.16 + tc) )
  40  k2 = 10^(2.84 - 2177.0 / (273.16 + tc) )
  50  if tc <= 25 then k3 = 10^(-5.86 - 317.0 / (273.16 + tc) )
  60  if tc > 25 then k3 = 10^(-1.1 - 1737.0 / (273.16 + tc) )
  70   t = 1
  80   if m0 > 0 then t = m/m0
  90   if t = 0 then t = 1
  100   moles = parm(1) * 0.1 * (t)^parm(2)
  110   moles = moles * (k1 * act("H+") + k2 * act("CO2") + k3 * act("H2O"))
  120   moles = moles * (1 - 10^(2/3*si_cc))
  130   moles = moles * time
  140  if (moles > m) then moles = m
  150 if (moles >= 0) then goto 200
  160  temp = tot("Ca")
  170  mc  = tot("C(4)")
  180  if mc < temp then temp = mc
  190  if -moles > temp then moles = -temp
  200 SAVE moles
  210 PUT(moles, 4)
 -end

Ca[13C]O3
# Calcite precipitates at equilibrium fractionation
-start
20 alpha = 10^(LK_NAMED("Log_alpha_13C_Calcite/CO2(aq)"))
30 dCa13CO3 = GET(4) * CALC_VALUE("R(13C)_CO2(aq)") * alpha
40 SAVE dCa13CO3
-end
 
SOLUTION 1 Xul Ha spring
    temp      28.7
    pH        6.62
    density   1
    C(4)      4.88   #244 ppm
    [13C](4)     -13
    Ca        13.62  #546 ppm
    Cl        5.55   #197 ppm
    K         .01    #0.3 ppm
    Mg        3.75   #91  ppm
    Na        1.74   #40  ppm
    S(6)      14.6   #1400 ppm

KINETICS 1 case A
 CO2(g)
   -m0 0
   -parm -3.4 0.001
 [13C]O2(g)
   -m0 0
   -parm -8
 Calcite
   -m0 0
   -parm 60 0.67
 Ca[13C]O3
   -m0 0
 -step 10*10 10*30 10*100 10*300
INCREMENTAL_REACTIONS
USER_GRAPH 1 Degassing CO2 and [13C]O2, calcite precip at eq fractionation
    -headings               time d13C_sol d13C_HCO3- d13C_CO2(aq) d13C_Calcite Apparent_d13C_gas pH
    -axis_titles            "Seconds" "d13C" "pH"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y ISO("[13C]")
30 GRAPH_Y ISO("R(13C)_HCO3-")
40 GRAPH_Y ISO("R(13C)_CO2(aq)")
50 GRAPH_Y (KIN_DELTA("Ca[13C]O3")/KIN_DELTA("Calcite")/0.0111802 - 1)*1e3
60 GRAPH_Y (SR("[13C]O2(g)")/SR("CO2(g)")/0.0111802 - 1) * 1000
70 GRAPH_SY -LA("H3O+")
  -end
    -active                 true
USER_GRAPH 2
    -headings               Time log10(PCO2) log10(P[13C]O2) SI("Calcite") Calcite
    -axis_titles            "Seconds" "SI" "Moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y SI("CO2(g)"), SI("[13C]O2(g)")
30 GRAPH_Y SI("Calcite")
40 GRAPH_SY KIN("Calcite")
  -end
    -active                 true
END