13C isotopes and TDIC on addition of soil CO2 to an infiltrating rainwater

[catafsi]

Hello everyone,
I want to model the effect on the 13C isotopic signature and TDIC off adding progressively soil CO2 to an infiltrating rainwater (according to Frondini et al., 2009). Also on each step soil CO2 is added to the water calcite dissolution at equilibirum conditions is allowed. I tried with adding 0.1 L of CO2 100 times. I made the code modifying a code that modeled precipitation of calcite due to degassification from Perez-Moreno et al.,(2021) and with information I found here on the posts of the forum.
I could run the code, but I'm not getting something that makes sense.

This is the code:
SOLUTION 1 #infiltrating rain water in equilibrium with atm CO2 (gas phase 1)
    temp      70
    pH        5 charge
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    C         0.000025
    -water    1 # kg
END
USE solution 1
GAS_PHASE 1 #atm CO2
    -fixed_volume
    -pressure 0.0001 #partial pressure of atm CO2
    CO2(g)         1
    [13C]O2(g)     0.0112416911 # -5.5 permil,  delta = (1 - Rsmp/Rstd)*1000, Rstd = 0.0111802
SAVE solution 1
END

USE solution 1 
GAS_PHASE 2 #soil CO2 addition
    -fixed_volume
    -pressure 0.01 #partial pressure of soil CO2
    -volume 0.1
    CO2(g)         1
    [13C]O2(g)     0.011459705 # -25 permil,  delta = (1 - Rsmp/Rstd)*1000, Rstd = 0.0111802
SOLID_SOLUTION 1
Calcite
-Comp Calcite 9.8283e-01
-Comp Ca[13C]O3(s) 1.1007e-02
END

#Generate graph
USER_GRAPH 1
-headings TC 13C(aq) 13C(CO2(g)) pH
-axis_titles "TDIC, mmoles" "PERMIL" "pH"
-initial_solutions false
-connect_simulations true
-plot_concentration_vs x
-start
10 GRAPH_X TOT("C")*1e3
20 GRAPH_Y ISO("[13C]"), ISO("R(13C)_CO2(g)")
30 GRAPH_SY -LA("H3O+")
-end
-active true
END
#add progressively soil CO2 until reaching equilibrium, copy and paste this part until reaching equilibrium
USE gas_phase 2
USE solution 1
SAVE solution 1
END
(...)
USER_PRINT
10 PRINT "13C units:            ", ISO_UNIT("[13C]")
20 PRINT "R(13C)_CO2(aq) units: ", ISO_UNIT("R(13C)_CO2(aq)")
END(code)

I would appreciate very much if you could help me,
Than you.

[dlparkhurst]

Your delta formula is incorrect. It should be delC = (Rsample/Rstd - 1)*1000.

Here is a script that generates a solution (~1L) in equilibrium with a basically infinite atmosphere at PCO2 -3.4, 25C, and isotopic composition -5.5 permil. To make the calculation in one run, the script calculates what the GAS_PHASE composition should be and writes it to a file, and then it imports that file as part of the calculation. Hopefully, this will get you started in the right direction.

Code: [Select]

SOLUTION 1 #infiltrating rain water in equilibrium with atm CO2 (gas phase 1)
    temp      25
    pH        5 charge
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
SELECTED_OUTPUT 101
-file gas_phase.pqi
USER_PUNCH 101
10 REM delC = (Rsample/Rstd - 1)*1000
20 Rsample = [-5.5 / 1000 + 1]*0.0111802       
30 P12CO2 = 10^-3.4
40 P13CO2 = P12CO2*Rsample
50 s$ = "GAS_PHASE 1" + EOL$
60 s$ = s$ + "  -fixed_volume" + EOL$
70 s$ = s$ + "  -volume 1e10" + EOL$
80 s$ = s$ + "CO2(g)     " + STR_E$(P12CO2, 15, 10) + EOL$ 
90 s$ = s$ + "[13C]O2(g) " + STR_E$(P13CO2, 15, 10) + EOL$     
100 PUNCH s$
END
SELECTED_OUTPUT 101
-active false
END
USE solution 1
INCLUDE$ gas_phase.pqi
SAVE solution 1
END


[catafsi]

Thank you very much! I'll try that

[dlparkhurst]

Note that you get the same solution composition if you use EQUILIBRIUM_PHASES to specify the partial pressures of CO2(g) and [13C]O2(g).

Code: [Select]

SOLUTION 1
EQUILIBRIUM_PHASES 1
CO2(g)     -3.4
[13C]O2(g) -5.3539456399253
END