Stop Precipitation from Changing the Solution Isotope Ratio

[King]

The model: I want to create a calcite saturated solution with a known isotopic composition, react it incrementally with CO2 with a known isotopic solution, then have it incrementally precipitate calcite and see the isotopic solution of that calcite, and fluid during each addition of CO2.

The problem: The solution is so oversaturated with calcite, that when it precipitates calcite it also removes so much 13C from the solution in the form of Ca13CO3, that it modifies the isotopic composition of the solution and obviously changes the isotopic composition of the calcite added during the next addition of CO2.

The question: How can I get the solution to only precipitate a small/fixed amount of calcite during each step? Say some tiny amount 0.00001 moles of calcite.

What I've tried (that doesn't work):
- EQUILIBRIUM_PHASES: but since I'm dealing with a solid solution, I need to set the amount of calcite and Ca13CO3 that I want to precipitate, and thereby fix the isotope ratio, which isn't what I want.
- Using MIX to remove some water and react it separately: it only changes the mass of the fluid, but the calcite still precipitates in large amounts relative to the fluid.
- Using COPY: It's only called after everything has already run.

What I haven't tried, Pseudocode below (that could work but seems inelegant, although I'd face the same issue with so much calcite being removed from the system that it modifies the isotopic composition of the fluid):

Code: [Select]

USE SOLUTION 2
REACT with known amount of CO2
SAVE SOLUTION 3
USE SOLUTION 3
REACT with calcite solid solution
OUTPUT isotope ratio of calcite
END

USE SOLUTION 3
REACT with known amount of CO2
SAVE SOLUTION 4
USE SOLUTION 4
... and repeat for 50 steps
The code:

Create a calcite saturated fluid with a d13C of -14 permil (note the addition of 0.0002216 [13C], which is 0.011079879*0.02[moles of Na2CO3], and gets us to a d13C of -14).

*For those following along at home, in the iso.dat database, the standard value of d13C is 0.011802 and not the current VPDB value of 0.01123720. For my models, it has been changed to the latter VPDB value to calculate the correct d13C values.

Code: [Select]

SOLUTION 1
    temp      5
    pH        6.998
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
REACTION 1
CaCl2 0.02
Na2CO3 0.02
[13C] 0.0002216
NaOH 0.2
SAVE SOLUTION 2
END
Then to figure out how many moles of CO2 will equilibrate it with it at an atmospheric CO2 concentration of 420 ppm (again calculate the R value at -8 and multiply that by 420 ppm [~4ppm] and convert it to log value):

Code: [Select]

USE SOLUTION 2
EQUILIBRIUM_PHASE
CO2(g) -3.37675071
[13C]O2(g) -5.329580927
END
Then we run the following code, which reacts the known amount of CO2 with a fixed d13C in 50 steps, equilibrates it with a calcite solid solution, then exports the isotope ratios of everything to a file.

Code: [Select]

USE SOLUTION 2
REACTION 2
CO2(g) 0.989254
[13C]O2(g) 0.0111473024
0.14425 moles in 50 steps
SOLID_SOLUTION
Calcite
-comp Calcite 0
-comp Ca[13C]O3(s) 0
SELECTED_OUTPUT
-file   ex20F.sel
-high_precision true
-distance false
-time false
      -temperature
-totals [13C] C
-isotopes R(13C)
        -isotopes R(13C)_Calcite
END

[dlparkhurst]

This may take a few iterations.

Here is a file that is based on your descriptions. I used a GAS_PHASE with a very large volume to represent the atmosphere. It should be more or less invariant and using GAS_PHASE adds more isotope ratios in the output file.

The last calculation uses the GAS_PHASE and SOLID_SOLUTION, which should precipitate calcite while maintaining the solution in equilibrium with the atmosphere. The calcite will be in isotopic equilibrium with the solution, but at a different isotope ratio because of fractionation. I think that is what should happen in an ideal world.

That would be one end member. At the other extreme, perhaps it the calcite precipitates quickly and there is no exchange with the atmosphere, the calcite would precipitate with no fractionation relative to the solution?

Anyway, maybe you could comment on what you want for this solution after calcite precipitation and before the incremental CO2 reaction.

For the next incremental CO2 reaction, I am confused. If you are adding CO2 to the solution, you would be dissolving calcite rather than precipitating. Is that what you have in mind, or are you allowing the solution to degas CO2 resulting in calcite precipitation.

Code: [Select]

ISOTOPES
C       
     -isotope        [13C]   permil  0.01123720              # VPDB, Vienna Pee Dee Belemnite
                                                             # Chang and Li, 1990, Chinese Science Bulletin
     -isotope        [13C](4)        permil  0.01123720      # VPDB, Vienna Pee Dee Belemnite
                                                             # Chang and Li, 1990, Chinese Science Bulletin
     -isotope        [13C](-4)       permil  0.01123720      # VPDB, Vienna Pee Dee Belemnite
C(4)   
     -isotope        [13C](4)        permil  0.01123720      # VPDB, Vienna Pee Dee Belemnite
                                                             # Chang and Li, 1990, Chinese Science Bulletin
END

SOLUTION 1
    temp      5
    pH        6.998 charge
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
REACTION 1
CaCl2 0.02
Na2CO3 0.02
Na2[13C]O3 0.0002216
NaOH 0.2
SAVE SOLUTION 2
END

USE SOLUTION 2
#EQUILIBRIUM_PHASE
# CO2(g) -3.37675071  10
# [13C]O2(g) -5.329580927 10
GAS_PHASE 1
    -fixed_volume
    -pressure 1
    -volume 100000000
    -temperature 5
    CO2(g)     0.000419999999615196
    [13C]O2(g) 4.68186700724408e-06
SOLID_SOLUTIONS 1
    Calcite
        -comp Calcite      0
        -comp Ca[13C]O3(s) 0
END


[King]

Thanks for the suggestion of using a large GAS_PHASE to simulate the atmosphere, that actually comes in handy with another scenario I'm dealing with.

The code you provided results in a calcite that precipitated in isotopic equilibrium from a fluid isotopically equilibrated with CO2 from the atmosphere. This is what would happen if all the calcite instantaneously precipitated from the fluid in isotopic equilibrium (barring any kinetic effects, but I'm not interested in that right now). The thing is, this doesn't actually happen.

Instead, a small amount of calcite is initially formed from the solution but not enough to appreciably modify the [Ca] or [C] concentration. Then this calcite recrystallizes in the supersaturated fluid. What appears to happen though, is that this fluid slowly incorporated CO2, thus driving its pH and d13C downwards, and the resulting calcite (either precipitated or recrystallized) progressively had a lower d13C. So the isotopic composition of the calcite that formed, is a function of the amount of CO2 added to that super saturated fluid (obviously there's kinetic effects, and recrystallization is messy and complicated, but these factors are extraneous to what I want to model here).

The addition of a little CO2 doesn't lead to dissolution in this case, it just lowers the SI of the calcite in the fluid. The fluid remains oversaturated with respect to calcite.

I'm interested in knowing what would be the isotopic composition of a calcite after each tiny addition of CO2 (but without precipitating so much calcite that it changes the isotopic composition of the fluid). From our results, and the modelling, it appears that only a small amount of CO2 leads to a large change in the isotopic composition of the fluid and the resulting calcite.

So the ideal modelling steps could be:

Code: [Select]

Make SATURATED SOLUTION 1

ADD small amount of CO2 with [13C] to SOLUTION 1
PRECIPITATE small amount of Calcite with [13C]
    OUTPUT d13C of Calcite and SOLUTION 1
SAVE SOLUTION 2

ADD small amount of CO2 with [13C] to SOLUTION 2
PRECIPITATE small amount of Calcite with [13C]
    OUTPUT d13C of Calcite and SOLUTION 1
SAVE SOLUTION 3

ADD small amount of CO2 with [13C] to SOLUTION 3
PRECIPITATE small amount of Calcite with [13C]
    OUTPUT d13C of Calcite and SOLUTION 1
...
The difficulty is precipitating a tiny amount of calcite without changing the isotopic composition of the solution.

I just had a thought:
When we use:

Code: [Select]

SOLID_SOLUTIONS 1
    Calcite
        -comp Calcite      0
        -comp Ca[13C]O3(s) 0Is the isotopic composition of the solid solution in equilibrium with:
1.) the initial isotopic solution of the fluid (before the calcite precipitates)?
2.) the final isotopic solution of the fluid (after the calcite precipitates)?
or
3.) the integral of isotopic solution of that fluid (during the calcite precipitation)?

If the answer is (1) then my question would be solved, because then I could run the code twice. First, without precipitation, to get the isotopic composition of the fluid after each CO2 addition. Then second with precipitation, which gives me the isotopic composition of the calcite at each addition.

[dlparkhurst]

I think you want to kinetically precipitate calcite and Ca[13C]O3 at the same time you kinetically dissolve CO2(g) and [13C]O2(g).

Here is an attempt at half the reactions, where the solids precipitate kinetically. I used a simple RATES definition for the precipitation of Calcite (where calcite represents the 12C solid). I then calculated the alpha fractionation factor to determine the isotopic ratio of the solid that would be in equilibrium with the solution and precipitated the corresponding amount of Ca[13C]O3 that produced the the right isotopic ratio. So, there are two kinetic reactions, one that precipitates calcite, and one that precipitates Ca[13C]O3. The calcite rate was arbitrarily scaled to approach equilibrium in about 1 day.

Compared to the previous calculations with GAS_PHASE and SOLID_SOLUTION, the difference is that the solids do not continue to react to stay in equilibrium with the solution. Only the current incremental precipitate is in isotopic equilibrium with the solution.

I used the equations from the NETPATH manual (Plummer, Prestemon, and Parkhurst, 1994, ~p. 18) for the equations for alphas. Hopefully, I translated them correctly, but you can check.

The addition of the other two RATES expressions for CO2(g) and [13C]O2(g) is left as an exercise for the user. :-) The relative rates of precipitation and ingassing are probably unknown, so I'm not sure how you deal with that unknown. Equilibrium precipitation or ingassing may also be subject to question.

Code: [Select]

ISOTOPES
C       
     -isotope        [13C]   permil  0.01123720              # VPDB, Vienna Pee Dee Belemnite
                                                             # Chang and Li, 1990, Chinese Science Bulletin
     -isotope        [13C](4)        permil  0.01123720      # VPDB, Vienna Pee Dee Belemnite
                                                             # Chang and Li, 1990, Chinese Science Bulletin
     -isotope        [13C](-4)       permil  0.01123720      # VPDB, Vienna Pee Dee Belemnite
C(4)   
     -isotope        [13C](4)        permil  0.01123720      # VPDB, Vienna Pee Dee Belemnite
                                                             # Chang and Li, 1990, Chinese Science Bulletin
END
CALCULATE_VALUES
Alpha_13C_Calcite/Soln # equilibrium for 13C: R(calcite)/R(soln)
     -start
10 alpha = -9999.999
20 if (TOT("C(4)") <= 0) THEN GOTO 1000
50 tot_CO2 = SUM_SPECIES("*CO2*","C")
60 tot_HCO3 = SUM_SPECIES("*HCO3*","C")
70 tot_CO3 = SUM_SPECIES("*CO3","C") - sum_species("*HCO3*","C")
80 tot_C = tot_CO2 + tot_HCO3 + tot_CO3
90 n_CO2 = tot_CO2 / tot_C
100 n_HCO3 = tot_HCO3 / tot_C
110 n_CO3 = tot_CO3 / tot_C
120 alpha_CO2_HCO3 = EXP( [ -0.91 + 6300 / TK^2 ]/ 1000 )
130 alpha_CO3_HCO3 = EXP( [ -3.4 + 870000 / TK^2 ]/ 1000 )
140 alpha_cal_CO2g = EXP( [ -3.63 + 11940000 / TK^2 ]/ 1000 )
150 alpha_HCO3_CO2g = EXP( [ -4.54 + 10990000 / TK^2 ]/ 1000 )
160 alpha_cal_HCO3 = alpha_cal_CO2g / alpha_HCO3_CO2g
170 alpha = alpha_cal_HCO3 / (n_CO2 * alpha_CO2_HCO3 + n_HCO3 + n_CO3 * alpha_CO3_HCO3)
1000 SAVE alpha
     -end
END
RATES
Calcite
-start
10 k = parm(1)
20 rate = k * (1 - SR("Calcite"))
30 moles = rate * TIME
40 SAVE moles
50 PUT(moles, 1)
60 END
-end

Ca[13C]O3
-start
10 moles = GET(1)
200 alpha = CALC_VALUE("Alpha_13C_Calcite/Soln")
210 r_soln = CALC_VALUE("R(13C)")
220 r_calcite = r_soln * alpha
230 SAVE moles * r_calcite
END
SOLUTION 1
    temp      5
    pH        6.998 charge
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
REACTION 1
CaCl2 0.02
Na2CO3 0.02
Na2[13C]O3 0.0002216
NaOH 0.2
SAVE SOLUTION 2
END

INCREMENTAL_REACTIONS
USE solution 2
KINETICS
#-cvode
Calcite
-m 0
-parm 5e-9
Ca[13C]O3
-m 0
-formula Ca[13C]O3
-step 86400 in 100
USER_GRAPH 1
    -headings               time permil_soln permil_ppt_step permil_ppt_avg SI(Calcite)
    -axis_titles            "Time, seconds" "Permil" "SI Calcite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
30 permil_soln = (CALC_VALUE("R(13C)") / 0.01123720 - 1) * 1000
40 permil_calcite = ((KIN_DELTA("Ca[13C]O3")/KIN_DELTA("Calcite")) / 0.01123720 -1) * 1000
50 permil_calcite_avg = ((KIN("Ca[13C]O3")/KIN("Calcite")) / 0.01123720 - 1) * 1000
60 GRAPH_Y permil_soln, permil_calcite, permil_calcite_avg
70 GRAPH_SY SI("Calcite")
  -end
    -active                 true
END