Modelling Ca-carbonate precipitation

[Bulbasaur]

Hi there, I am trying to make a code that will be able to change the amount of phosphorus in my solutions. I want to graph the results as concentration of P in solution (X axis) and the calcite SI (y-axis).

Below is the code I have so far:

Code: [Select]

TITLE Modelling calcite SI with increasing [P]

SOLUTION 1
    temp      7.6
    pH        12.36
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.14
    C         190
    Ca        0.029
    Cl        190
    Fe        4.8e-05
    K         3
    Li        0.32
    Mg        0.027
    Na        520
    S(6)      1.1e-05
    Si        57
    Sr        0.00099
    Zn        4.9e-05
    P         1E-5
    -water    1 # kg
END

KNOBS
-iter 100000
Use solution 1
REACTION 1
P 1
0.0001 mol in 100 steps
USER_GRAPH 1
-headings P calcite SI
-axis titles "P moles" "calcite SI"
-initial_solutions false
-connect_simulations true
-plot_concentration_vs x
-start
10 GRAPH_X P
20 GRAPH_Y SI("calcite")
-end
-active true
END

Thanks in advance,

B

[dlparkhurst]

Maybe you want to add pure phosphorus to your solution, but I expect something like Na2HPO4 is a more reasonable reactant. That is what is added in the following script.

Not sure exactly why the numerical method gets hung up, but I think it is because there is nothing to set the redox state. The SOLUTION_SPECIES definition adds a little bit of a redox species that creates redox buffer and allows the numerical method to converge.

Code: [Select]

TITLE Modelling calcite SI with increasing [P]
SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01; -log_k -9 # aids convergence
SOLUTION 1
    temp      7.6
    pH        12.36
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.14
    C         190
    Ca        0.029
    Cl        190
    Fe        4.8e-05
    K         3
    Li        0.32
    Mg        0.027
    Na        520
    S(6)      1.1e-05
    Si        57
    Sr        0.00099
    Zn        4.9e-05
    P         1E-5
    -water    1 # kg
END

Use solution 1
REACTION 1
Na2HPO4 1
0.0001 mol in 10 steps
USER_GRAPH 1
    -headings               P calcite SI
    -axis_titles            "P moles" "Calcite SI" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOT("P")
20 GRAPH_Y SI("calcite")
  -end
    -active                 true
END


[Bulbasaur]

Hi Dr. Parkhurst, thank you for the reply this is very helpful.

I have a couple more questions.

1) Is it possible to add the concentration of another ion (e.g. calcium to the secondary x-axis)? It is not working for me in this code below.

Code: [Select]

TITLE Modelling calcite SI with increasing [P]
SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01; -log_k -9 # aids convergence
SOLUTION 1
    temp      7.6
    pH        12.36
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.14
    C         190
    Ca        0.029
    Cl        190
    Fe        4.8e-05
    K         3
    Li        0.32
    Mg        0.027
    Na        520
    S(6)      1.1e-05
    Si        57
    Sr        0.00099
    Zn        4.9e-05
    P         1E-5
    -water    1 # kg
END

Use solution 1
REACTION 1
NaHPO4 1
0.01 mol in 10 steps
USER_GRAPH 1
    -headings               P    calcite SI   Ca
    -axis_titles            "P moles" "Calcite SI" "Ca"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOT("P"),
20 GRAPH_Y SI("calcite")
30 GRAPH_SX TOT("Ca")
  -end
    -active                 true
END
2) Would it be possible to simulate the mixing of this solution (solution 1) with another solution that contains NaHPO4 but add that solution in steps (possibly by volume)?

Thanks,

B

[dlparkhurst]

The secondary axis is "GRAPH_SY".

The volume correction is often pretty small, so adding Na2HPO4 with REACTION may be adequate. (Note Na2HPO4, not NaHPO4, which would have some odd redox reactions.)

It is possible to add a solution in steps, but it is a little arcane to do in general. The script below converts a SOLUTION into a REACTION by writing the mol/L of each element into a REACTION definition in a file. Once the file is included with INCLUDE$. The amount of reaction is then added in the script following the statement that includes the file. Adding 1 mol of the REACTION adds 1 L of the 1 mol/L Na2HPO4 solution (although the script will work with any solution). In this example, 0.01 moles of reaction is equivalent to adding 10 mL in 1 mL steps.

Code: [Select]

TITLE Modelling calcite SI with increasing [P]
SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01; -log_k -9 # aids convergence
SOLUTION 1
    temp      7.6
    pH        12.36
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.14
    C         190
    Ca        0.029
    Cl        190
    Fe        4.8e-05
    K         3
    Li        0.32
    Mg        0.027
    Na        520
    S(6)      1.1e-05
    Si        57
    Sr        0.00099
    Zn        4.9e-05
    P         1E-5
    -water    1 # kg
END
SOLUTION 2
-pH  7 charge
-units mol/L
-density 1 calc
Na 2
P  1
SELECTED_OUTPUT 1
-reset false
-file solution_as_reaction.pqi
USER_PUNCH 1
10 s$ = "REACTION 1" + EOL$
20 t = SYS("elements", count, name$, type$, moles )
30 FOR i = 1 TO count
40  IF (type$(i) <> "dis") OR (INSTR(name$(i), "(") > 0) THEN GOTO 100
50  s$ = s$ + PAD(name$(i), 10) + STR_E$(moles(i)/SOLN_VOL, 20, 10) + EOL$
100 NEXT i
110 PUNCH s$
END
SELECTED_OUTPUT 1
-active false
END
Use solution 1
INCLUDE$ solution_as_reaction.pqi
0.010 in 10 steps # 10 mL in 1 mL increments

USER_GRAPH 1
    -headings               P    SI(Calcite)   Ca
    -axis_titles            "mL of 1 M Na2HPO4" "Calcite SI" "Ca, mol/kgw"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN * 1000
20 GRAPH_Y SI("calcite")
30 GRAPH_SY TOT("Ca")
  -end
    -active                 true
END


[MichaelZ20]

Hi David!
It seems that there is something strange with element concentrations produced by the USER_PUNCH

[dlparkhurst]

Ah, the conversion from /L to /kgw is a little tricky. The density is used to get mass of solution in a liter, and the mass of dissolved solutes is subtracted to get the mass of water. Using the default value of the gfw for P (31) does not give the right mass for solutes for this relatively concentrated solution. Instead we should use the gfw of HPO4 in the SOLUTION definition, as follows:

Code: [Select]

P  1 as HPO4

The result is then

Code: [Select]

REACTION 1
H             1.0196716354e+02
O             5.4483581752e+01
Na            1.9999999902e+00
P             9.9999999512e-01


[Bulbasaur]

Thank you, Dr. Parkhurst.

This is very helpful to me!