ChemistryLab.jl

Appearance

Calculation of thermodynamic properties of calcite dissolution ​

Before this page

Getting started, whose calcite example this page rebuilds without a database.

This example is equivalent to quickstart except that no database is needed.

Let us recall the example which consisted of studying the dissolution of calcite in water.

It is possible to calculate the thermodynamic properties of the reaction, in particular the solubility constant of the reaction () which is related to the Gibbs free energy of the reaction (). This solubility constant is a function of temperature and the calculation is performed at a reference temperature of 298 K and at a pressure of 1 Atm using the following equation:

where is deduced from the Gibbs energies of formation () of the other chemical species involved in the reaction:

First step: species declaration ​

The first step, therefore, is to construct each of the species present in the reaction. This can be done with the Species function. Its simplest use is as follows:

julia
using ChemistryLab

calcite = Species("CaCO3", aggregate_state=AS_CRYSTAL, class=SC_COMPONENT)
Ca²⁺ = Species("Ca+2", aggregate_state=AS_AQUEOUS, class=SC_COMPONENT)
CO₃²⁻ = Species("CO3-2", aggregate_state=AS_AQUEOUS, class=SC_COMPONENT)
CO₃²⁻
Species{Int64}
           name: CO3-2
         symbol: CO3-2
        formula: CO3-2 ◆ CO₃²⁻
          atoms: C => 1, O => 3
         charge: -2
aggregate_state: AS_AQUEOUS
          class: SC_COMPONENT
     properties: M = 0.06000799997926252 kg mol⁻¹

The created object contains a certain amount of information whose properties can be entered a posteriori (or during the construction of the Species).

Calculation of Molar Mass

It can be noted that during the construction of the species, a calculation of the molar mass is systematically performed.

Second step: calculation of the thermodynamic properties of each species ​

For each species, it is possible to assign thermodynamic properties, such as the Gibbs energy of formation or the heat capacity. These data come from a database, here Thermoddem (Blanc et al., 2012), whose values for the three species at 25 °C are transcribed in data/literature/Blanc2012.json:

julia
using DynamicQuantities, Printf

thermoddem(sp) = literature_row("Blanc2012", "individual_properties", sp)
@printf("%-8s %12s %12s %12s %12s %12s\n", "", "ΔfG° kJ/mol", "ΔfH° kJ/mol", "S° J/K/mol", "Cp° J/K/mol", "V° cm³/mol")
for sp in ("Calcite", "Ca+2", "CO3-2")
    r = thermoddem(sp)
    @printf("%-8s %12.3f %12.3f %12.2f %12.2f %12.2f\n", sp, ustrip(us"kJ/mol", r.dfG), ustrip(us"kJ/mol", r.dfH),
            ustrip(us"J/(mol*K)", r.S), ustrip(us"J/(mol*K)", r.Cp), ustrip(us"cm^3/mol", r.V))
end
          ΔfG° kJ/mol  ΔfH° kJ/mol   S° J/K/mol  Cp° J/K/mol   V° cm³/mol
Calcite     -1129.109    -1207.605        91.78        83.47        36.93
Ca+2         -552.806     -543.000       -56.20       -26.38       -18.15
CO3-2        -527.900     -675.230       -50.00      -276.88        -5.28

Thermodynamic properties of formation ​

The first step involves associating the values of the thermodynamic properties of formation for each species. The values are read from data/literature/Blanc2012.json rather than typed, so that the page cannot drift from its source. Typing them by hand, as 83.47u"J/K/mol" for the heat capacity, gives the same entry. For calcite:

julia
cal = thermoddem("Calcite")
th_prop_0_calcite = Dict(:Cp⁰ => cal.Cp, :ΔₐH⁰ => cal.dfH, :S⁰ => cal.S, :ΔₐG⁰ => cal.dfG, :V⁰ => cal.V)
Dict{Symbol, Quantity{Float64, Dimensions{FRInt32}}} with 5 entries:
  :V⁰   => 3.6934e-5 m³ mol⁻¹
  :S⁰   => 91.78 m² kg s⁻² K⁻¹ mol⁻¹
  :ΔₐH⁰ => -1.2076e6 m² kg s⁻² mol⁻¹
  :ΔₐG⁰ => -1.12911e6 m² kg s⁻² mol⁻¹
  :Cp⁰  => 83.47 m² kg s⁻² K⁻¹ mol⁻¹

Heat capacity, enthalpy and free energy as a function of temperature ​

The second step is to describe the evolution of heat capacity as a function of temperature for each species. For calcite, the database expresses the heat capacity as a function of temperature,    , whose coefficients are transcribed in the same file. A reference temperature can then be defined in order to construct thermodynamic functions, such as heat capacity, entropy, enthalpy and free enthalpy. For calcite, this can be done as follows:

julia
mk = literature_row("Blanc2012", "maier_kelley", "Calcite")
params_Cp_calcite = Dict(:a₀ => mk.a, :a₁ => mk.b, :a₂ => mk.c)
T_ref = Dict(:T => 298.15u"K")
params_calcite = merge(th_prop_0_calcite, params_Cp_calcite, T_ref)
dtf_calcite = build_thermo_functions(:cp_ft_equation, params_calcite)
OrderedCollections.OrderedDict{Symbol, SymbolicFunc{1, @NamedTuple{T::Quantity{Float64, Dimensions{FRInt32}}}, Float64, Dimensions{FRInt32}}} with 4 entries:
  :Cp⁰  => 99.72 + 0.02692T + -2.158e6 / (T^2) [m² kg s⁻² K⁻¹ mol⁻¹] ◆ vars=(T) ◆ T=298.15 K
  :ΔₐH⁰ => -1.24577e6 + 99.72T + 2.158e6 / T + 0.01346(T^2) [m² kg s⁻² mol⁻¹] ◆ vars=(T) ◆ T=298.15 K
  :S⁰   => -496.549 + 0.02692T + 99.72log(T) + 1.079e6 / (T^2) [m² kg s⁻² K⁻¹ mol⁻¹] ◆ vars=(T) ◆ T=298.15 K
  :ΔₐG⁰ => -1.13991e6 + 596.269T + 1.079e6 / T - 0.01346(T^2) - 99.72T*log(T) [m² kg s⁻² mol⁻¹] ◆ vars=(T) ◆ T=298.15 K

Expression of Cp as a function of temperature

In the Thermoddem database, the expression for heat capacity as a function of temperature is written as:

However, the expression given in Thermoddem for this species is a simplification of the following more complete function:

This expression is implemented in ChemistryLab and can be called by passing :cp_ft_equation as an argument.

Calling the thermo_functions (here build_thermo_functions) allows the calculation of the expressions for the different thermodynamic properties as a function of temperature, according to the following expressions:

where and are the apparent enthalpy and free energy (Gibbs) at T.

The expressions for the thermodynamic properties of calcite can be added to the species calcite as follows. The symbolic expression is computed for each property, and the last line shows the one of the Gibbs energy of calcite:

julia
calcite.Cp⁰ = dtf_calcite[:Cp⁰]
calcite.ΔₐH⁰ = dtf_calcite[:ΔₐH⁰]
calcite.S⁰ = dtf_calcite[:S⁰]
calcite.ΔₐG⁰ = dtf_calcite[:ΔₐG⁰]
SymbolicFunc:
  Expression: -1.1399107839328907e6 + 596.2686733939438T + 1.079e6 / T - 0.01346(T^2) - 99.72T*log(T) [m² kg s⁻² mol⁻¹]
  Variables: T
  References: T=298.15 K
julia
using Plots

p1 = plot(xlabel="Temperature [°C]", ylabel="ΔₐG⁰ [J.mol⁻¹]", title="Gibbs energy of calcite \nas a function of temperature")
plot!(p1, θ -> calcite.ΔₐG⁰(T = 273.15+θ), 0:0.1:100, label="ΔₐG⁰ of calcite")

Similarly, we can provide the properties of the species and , printed above. Unlike that of calcite, the heat capacity of an aqueous ion follows the Helgeson-Kirkham-Flowers equations of state (HKF (Helgeson et al., 1981)). ChemistryLab has them (:solute_hkf88_reaktoro), but they take the HKF coefficients of each ion, which are not among the values transcribed here; we therefore hold the heat capacity at its value at 25 °C, -26.38 and -276.88 J mol⁻¹ K⁻¹ respectively.

julia
ca = thermoddem("Ca+2")
th_prop_0_Ca²⁺ = Dict(:Cp⁰ => ca.Cp, :ΔₐH⁰ => ca.dfH, :S⁰ => ca.S, :ΔₐG⁰ => ca.dfG, :V⁰ => ca.V)
params_Cp_Ca²⁺ = Dict(:a₀ => ca.Cp)   # the heat capacity at 25 °C, held constant
params_Ca²⁺ = merge(th_prop_0_Ca²⁺, params_Cp_Ca²⁺, T_ref)
dtf_Ca²⁺ = build_thermo_functions(:cp_ft_equation, params_Ca²⁺)
Ca²⁺.ΔₐG⁰ = dtf_Ca²⁺[:ΔₐG⁰]

co3 = thermoddem("CO3-2")
th_prop_0_CO₃²⁻ = Dict(:Cp⁰ => co3.Cp, :ΔₐH⁰ => co3.dfH, :S⁰ => co3.S, :ΔₐG⁰ => co3.dfG, :V⁰ => co3.V)
params_Cp_CO₃²⁻ = Dict(:a₀ => co3.Cp)
params_CO₃²⁻ = merge(th_prop_0_CO₃²⁻, params_Cp_CO₃²⁻, T_ref)
dtf_CO₃²⁻ = build_thermo_functions(:cp_ft_equation, params_CO₃²⁻)
CO₃²⁻.ΔₐG⁰ = dtf_CO₃²⁻[:ΔₐG⁰]
SymbolicFunc:
  Expression: -460255.728 - 1804.4305786067766T + 276.88T*log(T) [m² kg s⁻² mol⁻¹]
  Variables: T
  References: T=298.15 K

A constant heat capacity for the ions

The expressions for enthalpies and free energies remain temperature-dependent thanks to the integration performed on Cp. Within the temperature and pressure ranges typically considered in ChemistryLab (0–100 °C, 1 atm), assuming a constant Cp has little impact on the solubility product.

Third step: writing the reaction ​

Now we can write the dissolution/precipitation reaction of calcite.

julia
r = Reaction([calcite, Ca²⁺, CO₃²⁻]; equal_sign='↔')
  equation: CaCO₃ ↔ Ca²⁺ + CO₃²⁻
 reactants: CaCO₃ => 1
  products: Ca²⁺ => 1, CO₃²⁻ => 1
    charge: 0
properties: ΔᵣG⁰ = 117958.22293289064 - 2521.1818533574337T + -1.079e6 / T + 0.01346(T^2) + 402.98T*log(T) [m² kg s⁻² mol⁻¹] ◆ vars=(T) ◆ T=298.15 K
            logK⁰ = (1.079e6 - 117958.22293289064T + 2521.1818533574337(T^2) - 402.98(T^2)*log(T) - 0.01346T*(T^2)) / (19.144757680815896(T^2)) [m² kg s⁻² mol⁻¹] ◆ vars=(T) ◆ T=298.15 K

During the construction of this reaction, the thermodynamic properties of the reaction are calculated.

julia
using Plots

p1 = plot(xlabel="Temperature [°C]", ylabel="pKs", title="Solubility product (pKs) of calcite \nas a function of temperature")
plot!(p1, θ -> r.ΔᵣG⁰(T = 273.15+θ) / (R_GAS * (273.15+θ)) / log(10), 0:0.1:100, label="pKs")

Notes and next steps ​

This example demonstrates the manual workflow: create species, attach thermodynamic data from an external source, build a reaction and evaluate its temperature-dependent properties.

In practice, loading species from a database (see Importing thermodynamic databases) is faster and less error-prone:

julia
all_species = build_species(datapath("cemdata18-thermofun.json"))
species = speciation(all_species, split("Cal H2O@");
              aggregate_state=[AS_AQUEOUS], exclude_species=split("H2@ O2@ CH4@"))

Where to go next ​

The equilibrium state of the same system, with its species amounts, its pH and the saturation index of calcite, follows from a ChemicalSystem built on these species and passed to equilibrate; the tutorial Chemical Equilibrium takes that workflow from end to end.