# Facing the challenge of predicting the standard formation enthalpies of
n-butyl-phosphate species with ab initio methods

Mohamad Saab  Laboratoire PhLAM, CNRS UMR 8523, Université de Lille, 59655
Villeneuve d’Ascq Cedex, France  Florent Réal  Laboratoire PhLAM, CNRS UMR
8523, Université de Lille, 59655 Villeneuve d’Ascq Cedex, France  Martin Šulka
Slovak University of Technology in Bratislava, ATRI, Faculty of Materials
Science and Technology in Trnava, Bottova 25, 917 24 Trnava, Slovak Republic
Laurent Cantrel  Institut de Radioprotection et de Sûreté Nucléaire (IRSN),
PSN-RES, Cadarache, St Paul Lez Durance, 13115, France  François Virot
Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-RES, Cadarache,
St Paul Lez Durance, 13115, France  Valérie Vallet  valerie.vallet@univ-
lille1.fr Laboratoire PhLAM, CNRS UMR 8523, Université de Lille, 59655
Villeneuve d’Ascq Cedex, France

###### Abstract

Tributyl-phosphate (TBP), a ligand used in the PUREX liquid-liquid separation
process of spent nuclear fuel, can form explosive mixture in contact with
nitric acid, that might lead to violent explosive thermal runaway. In the
context of safety of a nuclear reprocessing plant facility, it is crucial to
predict the stability of TBP at elevated temperatures. So far, only the
enthalpies of formation of TBP is available in the literature with a rather
large uncertainties, while those of its degradation products, di-(HDBP) and
mono-(H2MBP) are unknown. In this goal, we have used state-of-the art quantum
chemical methods to compute the formation enthalpies and entropies of TBP and
its degradation products di-(HDBP), mono-(H2MBP) in gas and liquid phases.
Comparisons of levels of quantum chemical theory revealed that there are
significant effects of correlation on their electronic structures, pushing for
the need of not only high level of electronic correlation treatment, namely
local coupled cluster with single and double excitation operators and
perturbative treatment of triple excitations [LCCSD(T)], but also
extrapolations to the complete basis to produce reliable and accurate
thermodynamics data. Solvation enthalpies were computed with the conductor
like screening model for real solvents [COSMO-RS], for which we observe errors
not exceeding $22\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. We thus
propose with final uncertainty of about
$20\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ standard enthalpies of
formation of TBP, HDBP, and H2MBP which amounts to $-1281.7\pm 24.4$,
$-1229.4\pm 19.6$ and $-1176.7\pm
14.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, respectively, in the
gas phase. In the liquid phase, the predicted values are $-1367.3\pm 24.4$,
$-1348.7\pm 19.6$ and $-1323.8\pm
14.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, to which we may add
about $-22\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ error from the
COSMO-RS solvent model. From these data, the complete hydrolysis of TBP is
predicted as an exothermic phenomena but showing a slightly endergonic
process.

##  I Introduction

All spent nuclear fuel reprocessing plants use the PUREX process (Plutonium
Uranium Refining by Extraction), which is a liquid-liquid extraction method.
The organic extracting solvent is a mixture of tri-n-butyl phosphate (TBP) and
hydrocarbon solvent such as hydrogenated tetra-propylene (TPH). By chemical
complexation, uranium and plutonium (from spent fuel dissolved in nitric acid
solution), are separated from fission products and minor actinides. During a
normal extraction operation, uranium is extracted in the organic phase as the
UO2(NO3)2(TBP)2 complex.

The TBP extractant can form an explosive mixture called red oil when it comes
in contact with nitric acid. The formation of this unstable organic phase
originates from the reaction between TBP and its degradation products on the
one hand, and nitric acid, its derivatives and heavy metal nitrate complexes
on the other hand. The decomposition of the red oil can lead to violent
explosive thermal runaway. These hazards are at the origin of several
accidents such as, the two in the United States in 1953 and 1975 (Savannah
River) Colven et al. (1953); Durant (1983) and, more recently, the one in
Russia in 1993 (Tomsk). IAE (1998) The accidental consequences were the
rupture of the equipment, followed by a significant radioactive release to the
environment. This raises the question of the exothermicity of reactions that
involve TBP and all other degradation products, and calls for a better
knowledge of the underlying chemical phenomena.

TBP is relatively stable, but at elevated temperatures, its degradation may
occur, as a result of the following phenomena: radiolysis, pyrolysis, acid
catalyzed hydrolysis, as well as the de-alkylation (see the therein reference
for a detailed reaction scheme.Wright and Paviet-Hartmann (2010)) The two last
lead to the stepwise formation of degraded organic products such as dibutyl
phosphate (HDBP), mono-butyl phosphate (H2MBP), and ultimately phosphoric acid
(H3PO4). The reaction by-products (butanol and butyl nitrate) may then undergo
oxidation to form carboxylic acids and gases. In a similar way,
UO2(NO3)2(TBP)2 may be decomposed to UO2(NO3)2TBP(HDBP) or UO2(NO3)2(HDBP)2.
To understand the chemistry of these systems, the knowledge of some physical
properties of such species is necessary, such as thermodynamic functions or
kinetics rate law of degradation. However, such data are rather scarce in the
available literature. Even if there are a plethora of studies on TBP-nitric
acid systems (see reviews Wright and Paviet-Hartmann, 2010; Burger, 1990),
some of the fundamental properties, like standard enthalpies of formation or
kinetics of by-product degradations, are poorly known.

Only three articles have reported experimental investigations on standard
enthalpy of formation of TBP and there is no data for degraded organic
products. These published data were obtained by combustion method in a bomb
calorimeter. Starostin et al.Starostin et al. (1966) first studied the TBP
combustion with static isothermal calorimeter, the value, $-1458\pm
12\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ was obtained. Later,
KindleKindle (1974) reported measurements carried on adiabatic bomb
calorimeter targeting the amount of heat of combustion in a nuclear facility
for safety analysis. Even if initially the data are not intended for the
calculation of enthalpy of formation, an evaluation is possible by supposing a
complete combustion: $-1250\pm
69\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. The most recent study,
carried out by Erastov and TarasovErastov and Tarasov (1992) derived from
rotating bomb calorimeter measurements, the standard enthalpy of formation of
TBP; the value $-1382\pm 12\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$
was retained. Despite similar approaches, experimental data exhibit a large
uncertainty over a range of about
$200\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. Therefore the aim of
our safety strategy of red oil issues is split into several steps where the
first investigation is the calculation of standard enthalpy of formation of
TBP in the TBP solvent as well as its degraded products, HDBP and H2MBP.

In the first part of this article, the quantum chemistry calculations are
introduced. The choice of electronic correlation treatment and a benchmark of
basis sets are discussed as the methodology to evaluate the solvent
contributions. In the second part, the standard enthalpy of formation of
target species in both gas and liquid phases, as well as the hydrolysis
reaction free energies, are presented and discussed regarding the experimental
data.

##  II Computational Details

###  II.1 Chemical reactions used to derive the standard enthalpies of
formation of H_3-nnBP

The standard enthalpies of formation in the gas and TBP liquid phases of the
H_3-nnBP molecules, H2MBP (n=1), HDBP (n=2), and TBP (n=3) are derived from
computed enthalpies of reactions.

(1) (2) (3)

together with the standard enthalpies of formation reported in Table 1.

Table 1: Known Gas/Liquid Phase Standard Enthalpies of Formation $\Delta_{f}H^{\ominus}$($298.15\text{\,}\mathrm{K}$) in $\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ Molecule |  $\Delta_{f}H^{\ominus}$($298.15\text{\,}\mathrm{K}$) | Reference
---|---|---
HF(g) | $-273.3\pm 0.7$ | 9
HCl(g) | $-92.31\pm 0.10$ | 9
H2O(g) | $-241.826\pm 0.040$ | 9
C4H9OH(g) | $-277\pm 5.$ | 10
PH3(g) | $5.47\pm 1.7$ | 11
POF3(g) | $-1254.25\pm 8$ | 11
POCl3(g) | $-559.82\pm 1.7$ | 11
H3PO4(l) |  $-1271.66\pm 2.9$111Details on how the reported uncertainty was estimated are given in the supplementary material. Further comparisons with computed gas-phase formation enthalpies are also included in supplementary material | 9

Note that all these reactions involve closed-shell molecules. As discussed by
WattsWatts (2016) and illustrated by HehreHehre (2003) it tends to be easier
to obtain accurate reaction enthalpies of isodesmic reactions than for an
isogyric reaction. However, only one isodesmic chemical reaction, reaction 1,
was identified leading to the formation of TBP and its derivatives. To verify
the accuracy of the derived enthalpies of formation, also considered three
isogyric reactions 2, 3 and LABEL:rxn:H2, that include chemical compounds with
known enthalpies of formation.

###  II.2 Choice of the atomic basis sets and the treatment of electron
correlation

The optimal structures of all species were optimized in the gas phase without
symmetry constraints using the B3LYP density functional.Becke (1993); Lee et
al. (1988); Vosko et al. (1980); Stephens et al. (1994) Triple-$\zeta$ valence
basis sets def2-TZVP basis setsWeigend et al. (1998); Weigend and Ahlrichs
(2005) were used for all atoms. To describe butyl-phosphate H_3-nnBP(l)
molecules in the liquid phase, solvent effects have to be considered. However,
we noted that the TBP solvent, treated by the COSMO continuum model,Klamt and
Schüürmann (1993) induce only minor changes on the optimal structures of the
molecules. Hence, gas-phase optimized geometries will be used throughout to
compute the harmonic frequencies that enter, without being scaled, the
vibrational partition functions at $298.15\text{\,}\mathrm{K}$ and 0.1 MPa.
The latter are necessary to calculate the enthalpic and entropic contributions
to the gas-phase energies. Perspective views of the optimal geometries of
H2MBP, HDBP and TBP are displayed on Figure 1.

(a) H2MBP

(b) HDBP

(c) TBP

Figure 1: Perspective views of the geometries of H2MBP, HDBP and TBP optimized
at the B3LYP-def2-TZVP level of theory. Hydrogen atoms are in white, carbon
atoms in black, oxygen atoms in red, and phosphorous in blue.

The derivation of unknown formation enthalpies requires to compute reaction
enthalpies, which are composed of reaction energies and enthalpic corrections.
The accuracy is mostly determined by the errors in the computed electronic
reaction energies with respect to both the basis set (basis-set error) and the
treatment of electron correlation (intrinsic error with respect to the exact
energy), by considering the reactions 1, 3 and LABEL:rxn:H2 leading to the
formation of the smallest butyl phosphate H2MBP. The other di and tributyl
phosphate are difficult in this respect due to their size, which makes it
impossible to consider huge basis sets very close to complete basis set (CBS)
limit, at least for the correlated WFT methods, although large basis sets are
much less prone to basis set superposition errors (BSSE). On the other hand,
smaller basis sets have to be checked for the remaining basis set
incompleteness error (BSIE) and BSSE. While the BSIE can be estimated with
extrapolation schemes, the intramolecular BSSE is not directly accessible for
the investigated formation reactions listed above.

We have used either the triple and quadruple-$\zeta$ valence basis sets of
def2-type (def2-TZVP and def2-QZVP)Weigend et al. (1998); Weigend and Ahlrichs
(2005) to extrapolate to the CBS limit,Neese and Valeev (2011) or the sequence
of augmented correlation consistent basis sets, aug-cc-pV$n$Z, from triple
($n$=3) to quintuple-$\zeta$ ($n$=5) quality with subsequent CBS
extrapolation. Correlation effects are treated either with the B3LYP
functional or with correlated wave-function theory (WFT) methods of increasing
accuracy, namely MP2 and the single and double coupled cluster theory with
inclusion of a perturbative estimation for triple excitation [CCSD(T)], the
latter representing the “gold standard”. All WFT calculations are performed
with either the parallel resolution of the identity approximation,Hättig and
Weigend (2000); Hättig et al. (2006) or density fitting correlated
methods,Werner and Schütz (2011) with the appropriate atomic auxiliary basis
functions,Weigend et al. (1998); Weigend and Ahlrichs (2005) and the frozen-
core approximation (chemical core, that is, only the valence electrons were
correlated). For the aug-cc-pVnZ basis sets, extrapolation to the CBS limit
was carried out according to the three-point exponential formula for HF
energies $E_{\text{HF}}$:Feller (1992, 1993)

$E_{\text{HF}}(n)=E^{\text{CBS}}_{\text{HF}}+A\exp{(-Bn)},$ (13)

and a two-point extrapolation for the total DFT energies or the WFT
correlation energies $E_{\text{corr}}$:Helgaker et al. (1997)

$E_{\text{corr}}(n)=E^{\text{CBS}}_{\text{corr}}+An^{-3}.$ (14)

If calculations on the monobutyl-phosphate molecule (20 atoms) can be
performed with ”standard” (using canonical molecular orbitals) post-Hartree-
Fock MP2 and CCSD(T) methods, their prohibitive scaling, $O(N^{5})$ and
$O(N^{7})$ with $N$ the measure of the molecular size or the number of
correlated electrons, renders calculations on the dibutyl and tributyl
phosphate ligands unfeasible. However, the efficient local CCSD(T) (LCCSD(T))
method with density fitting approximation of the integrals offers a favorable
$N$ scaling.Werner and Schütz (2011) To reach the highest accuracy in local
coupled cluster treatments, Schwilk et al.Schwilk et al. (2015) advocate to
include all close pairs amplitudes in the LCCSD (keepcls=2 option in
MOLPRO.Werner et al. (2015, 2012)) As discussed by Werner and Schütz,Werner
and Schütz (2011) for reactions involving large molecules, one might observe
significant differences between the results of local and canonical
calculations, originating either from the local approximations or the larger
BSSE effects in the canonical calculations. Therefore, it is recommended to
add a correction $\Delta
E_{\mathrm{MP2(CBS)}}=E_{\mathrm{MP2(CBS)}}-E_{\mathrm{LMP2}}$, to the
LCCSD(T) energies, where $E_{\mathrm{MP2(CBS)}}$ is the CBS extrapolated
canonical MP2 energy using aug-cc-pV$n$Z basis sets ($n$ = 3–5), and
extrapolation formulas reaction 13 and reaction 14. The LMP2 and LCCSD(T)
calculations are computed with triple-$\zeta$ basis sets. We will note this
energy $\Delta E_{\mathrm{LCCSD(T)(CBS^{\ast})}}$:

$\Delta
E_{\mathrm{LCCSD(T)(CBS^{\ast})}}=E_{\mathrm{LCCSD(T)}}+E_{\mathrm{MP2(CBS)}}-E_{\mathrm{LMP2}}$
(15)

All calculations were performed with Turbomolepro  and MOLPRO.Werner et al.
(2015, 2012)

###  II.3 Benchmarking basis sets and correlation methods for the H2MBP
formation reactions 1, 3 and LABEL:rxn:H2

Table 2:  Reaction energies $\Delta_{r}E$ in $\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ for H2MBP formation reactions 1, 3 and LABEL:rxn:H2, using different methods and basis sets. Method | Basis set |  $\Delta_{r}E$ (1) |  $\Delta_{r}E$ (3) |  $\Delta_{r}E$(LABEL:rxn:H2)
---|---|---|---|---
B3LYP | def2-TZVP | 0.0 | -131.1 | -105.5
B3LYP | def2-QZVP | -2.3 | -124.2 | -103.3
B3LYP | CBS | -4.0 | -119.3 | -101.8
B3LYP | def2-QZVPP | -2.3 | -124.2 | -103.3
B3LYP | aug-cc-pVTZ | -3.9 | -108.9 | -51.7
B3LYP | aug-cc-pVQZ | -3.7 | -116.3 | -78.1
B3LYP | aug-cc-pV5Z | -3.2 | -121.1 | -97.1
B3LYP | CBS | -2.8 | -126.0 | -117.0
MP2 | def2-TZVP | -8.0 | -135.7 | -202.5
MP2 | def2-QZVP | -10.0 | -119.0 | -155.9
MP2 | CBS | -11.4 | -106.7 | -121.9
MP2 | aug-cc-pVTZ | -13.8 | -112.0 | -110.5
MP2 | aug-cc-pVQZ | -12.3 | -112.8 | -141.5
MP2 | aug-cc-pV5Z | -11.6 | -113.3 | -161.6
MP2 | CBS | -11.1 | -111.1 | -177.9
LCCSD(T) | CBS∗ | -10.3 | -124.8 | -136.0

Reaction 1 is isodesmic, implying that the type of chemical bonds broken in
the reactants are the same as the type of bonds formed in the reaction
products. Such type of reactions are used preferentially in thermochemistry to
calculate accurate enthalpies of formation. As we mentioned previously, as
only one isodesmic reaction was found, we considered as well three non-
isodesmic reactions that are isogyric: both reactions 3 and 2 retain the P=O
bond, while LABEL:rxn:H2 does not. We can thus expect better error
cancelations in the former two than in the latter. Table 2 gathers the
computed energies for reactions 1, 3 and LABEL:rxn:H2, leading to the
formation of H2MBP with different basis sets and DFT and WFT methods.

Let us begin by scrutinizing the results for reaction 1. Both def2-basis and
aug-cc-pV$n$Z basis set families yield very similar CBS values with the B3LYP
functional. This is also true at the MP2 level, though the CBS reaction
energies are about $-7.5\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ more
exothermic that at the B3LYP level. The LCCSD(T) agrees within
$0.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ with the MP2 one. This
suggests that electron correlation effects essentially cancel out in this
isodesmic reaction, and confirms that the fact that error compensations in
such reactions are optimal.

With respect to the isogyric reactions, we first focus on reaction 3. With the
B3LYP functional, the reaction energy exhibits an unexpected strong basis set
dependence with the def2 basis sets. The reaction energies are not converged
with the largest basis set QZVP. Extending the def2-QZVP basis set with extra
polarization functions (def2-QZVPP) has no effect on the computed values. The
other surprising observation is that the aug-cc-pV$n$Z basis sets also
converge slowly to the basis set limit, together with the fact that both def2-
and cc-type basis sets lead to CBS values that differ by
$6.7\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. Whenever large basis
set dependences are observed at the DFT level, one may suspect large BSSE even
with very extended atomic orbital basis that results from numerical
instabilities of the DFT functionals, as discussed recently by Hansen et
al.Hansen et al. (2014) This makes us conclude that BSIE and BSSE effects are
crucial with the B3LYP functional for the various reactants and products of
reactions 3 and LABEL:rxn:H2. The electron correlation effects appear to be
significantly different in the two sides of the two reaction equations, thus
making error compensations not favorable, and forcing us to push the treatment
of dynamic correlation up to the CBS limit.

At the MP2 level, while the def2 basis sets reaction energy values converge
very slowly to its CBS limit, the aug-cc-pV$n$Z basis sets are far less prone
to BSIE for reaction 3, being only
$0.9\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ apart from the CBS level
with triple-$\zeta$ basis set. MP2/aug-cc-pV$n$Z energies for LABEL:rxn:H2
converge more slowly to its CBS limit. Because the def2-CBS limit
significantly differs from the aug-cc-pV$n$Z one, we cross-checked these
results with relativistic (Douglas Kroll Hess, DK) correlation consistent
basis sets (aug-cc-pV$n$Z-DK) and confirmed the non-relativistic aug-cc-pV$n$Z
values. We thus infer that the very different bonding around the phosphorous
atom in the reactant PH3 (three single P-H bonds) and in the product H2MBP
(one double P=O bond, three single P-O bonds) cannot be accurately captured by
def2-basis sets.

To push the treatment of electron correlation further, we performed LCCSD(T)
calculations extrapolated to the CBS limit as explained earlier in this
section. The resulting reaction energies estimated to the CBS limit are
$-124.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ and
$-136.0\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ for reaction 3 and
LABEL:rxn:H2, respectively. Though we cannot compare this value to
experimental values, we trust LCCSD(T)/CBS results as the best estimate for
CCSD(T)/CBS accuracy, and we will see in the next section that this level of
theory helps us deriving accurate standard enthalpies of formation.

###  II.4 COSMO-RS method to estimate solvent contributions to thermodynamic
quantities

Whenever the molecules considered are in the liquid phase, we have to estimate
their thermodynamic properties in the corresponding phase. While most
continuum solvent models only estimate Gibbs Free energy of solvation, the
Conductor-like screening model-Real Solvents (COSMO-RS) modelEckert and Klamt
(2002) is the only one that allows us to access enthalpies of solvation, which
in a pure solvent are strictly equal to the opposite of the enthalpies of
vaporization enthalpy, or for a solvated molecule can be extracted as the
derivative of the solvation free energy by finite differences with two
different temperatures ($25\text{\,}\mathrm{\SIUnitSymbolCelsius}$,
$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$).Diedenhofen et al. (2007) In order
to obtain reliable thermophysical property data with COSMO-RS predictions,
accurately optimized structures with corresponding gas-phase and COSMO state
(ideal dielectric continuum by using COSMO with $\epsilon_{r}=\infty$)
energies are necessary.Diedenhofen et al. (2007) These are performed with the
BP functional of the densityBecke (1988); Perdew (1986a, b) and def-TZVP basis
sets.Schäfer et al. (1994) All statistical COSMO-RS thermodynamic
calculationsEckert and Klamt (2002) were performed with COSMOtherm,Eckert and
Klamt  using the parameter file BP TZVP“C30_1701”.

Table 3:  COSMO-RS vaporization enthalpies of all relevant species in $\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, together with experimental values when available. Molecule | liquid phase |  $\Delta_{vap}H^{\ominus}$($298.15\text{\,}\mathrm{K}$) |
---|---|---|---
|  | COSMO-RS | Exp. | Ref.
C4H9OH |  C4H9OH | 58.8 | $52.34\pm 0.02$ | 40
H2O |  H2O | 49.2 | $43.99\pm 0.07$ | 41
TBP |  TBP | 85.6 | 64.4 | 5
HDBP |  TBP | 118.5 |  |
|  HDBP | 119.3 |  |
H2MBP |  TBP | 159.5 |  |
|  H2MBP | 147.1 |  |
H3PO4 |  H3PO4 | 144.7 |  |

To evaluate the accuracy of COSMO-RS out-coming thermodynamic quantities, we
first start by comparing the calculated enthalpies of vaporization
$\Delta_{vap}H^{\ominus}$ to available experimental values for butanol, water
and TBP. The results reported in Table 3 shows that the COSMO-RS
$\Delta_{vap}H^{\ominus}$ values overestimate by 6.5 and
$5.2\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ the experimental ones
for butanol and water, respectively. For TBP, it overshoots by
$21\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. Such deviations are
within the range of the reported range of accuracy.Diedenhofen et al. (2007);
Schröder and Coutinho (2014) For H3PO4, there does not exist an experiment
value for the enthalpy of formation in the gaseous phase. Thus, we have
calculated the enthalpy of vaporization $\Delta_{vap}H^{\ominus}$ for H3PO4
using COSMO-RS and deduced the enthalpy of formation in the gaseous phase from
that in the liquid phase and the calculated enthalpy of vaporization. This
leads us to a gas-phase $\Delta_{f}H^{\ominus}$($298.15\text{\,}\mathrm{K}$)
value of $-1126.96\pm 2.9\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, in
very good agreement with the value of
$-1141.4\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ calculated at the
CCSD(T)-CBS level by Alexeev et al.Alexeev et al. (2005)

We therefore conclude that COSMO-RS can predict solvation contributions to a
reaction energy with good confidence and help us to determine the unknown
formation enthalpies of H2MBP, HDMP and TBP in the TBP organic phase.

##  III Results and Discussion

Table 4: $\Delta_{f}H^{\ominus}$($298.15\text{\,}\mathrm{K}$) of TBP, HDBP, and H2MBP in the gaseous phase calculated from the isodesmic reaction 1 and the average of the three isogyric reactions 2, 3 and LABEL:rxn:H2 listed in LABEL:tab:Hf:isogyric. $\overline{\Delta_{f}H^{\ominus}}$($298.15\text{\,}\mathrm{K}$) represent the average between these two data sets. The standard deviations $\mathrm{\Delta}$$\Delta_{f}H^{\ominus}$ includes the uncertainty on the average value as well as the experimental error bars. Values in the liquid state are given by $\Delta_{f}H^{\ominus}$($298.15\text{\,}\mathrm{K}$)= $\overline{\Delta_{f}H^{\ominus}}$($298.15\text{\,}\mathrm{K}$)-$\Delta_{vap}H^{\ominus}$($298.15\text{\,}\mathrm{K}$). All results are in $\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. Species | Reaction |  $\Delta_{f}H^{\ominus}$($298.15\text{\,}\mathrm{K}$) |  $\mathrm{\Delta}$$\Delta_{f}H^{\ominus}$
---|---|---|---
TBP(g) | reaction 1 | -1278.4 | 18.0
reaction 2–LABEL:rxn:H2 | -1285.0 | 23.8
$\overline{\Delta_{f}H^{\ominus}}$($298.15\text{\,}\mathrm{K}$) = | -1281.7 | 24.4
TBP(l) |  | -1367.3
HDBP(g) | reaction 1 | -1226.1 | 13.0
reaction 2–LABEL:rxn:H2 | -1232.7 | 19.0
$\overline{\Delta_{f}H^{\ominus}}$($298.15\text{\,}\mathrm{K}$) = | -1229.4 | 19.6
HDBP(l) |  | -1348.7
H2MBP(g) | reaction 1 | -1173.4 | 7.9
reaction 2–LABEL:rxn:H2 | -1180.0 | 14.3
$\overline{\Delta_{f}H^{\ominus}}$($298.15\text{\,}\mathrm{K}$) = | -1176.7 | 14.8
H2MBP(l) |  | -1323.8

First, the standard enthalpies of formation of H2MBP, HDBP and TBP species are
determined in gaseous state from reactions 1, 2, 3 and LABEL:rxn:H2.
Electronic energies of involved species are computed with the LCCSD(T)/CBS
approach whereas the vibrational contributions are calculated at the
B3LYP/def2-TZVP level (like geometry optimization) to ultimately derive the
enthalpies of reaction. Note that we are confident in the quality of the
partition functions and computed vibrational spectra, because of the estimated
values of enthalpy changes, entropies and heat capacities reported in Tables
S3 and S5 of the supplementary material agree excellently with experimental
data. To transform the standard enthalpy of reaction towards the standard
enthalpy of formation, the standard heats of formation at
$298.15\text{\,}\mathrm{K}$ of key species (See Table 1) are taken into
account. Secondly the computation of solvation enthalpies of target species is
performed with the COSMO-RS method, as introduced in Section II.4. Hence the
combination of two last computed values leads the standard enthalpy of
formation in the TBP solvent.

###  III.1 Heat of formation of TBP

To validate our approach in computing enthalpies of formation, the type of the
chosen chemical reaction might be critical. Inspection of Table 4 reveals that
the three isogyric reactions 2, 3 and LABEL:rxn:H2 yields values of enthalpies
of formation that are very close to each other. For TBP in the gaseous state,
the average over these three reactions (See Table S2 in the supplementary
material) is $-1285.0\pm
23.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, which compares
extremely well with the value derived with the isodesmic reaction 1,
$-1278.4\pm 18.0\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. The same is
true for the formation enthalpies of HDBP and H2MBP (Comparison of the
isogyric averages in Table S2 of the supplementary material and the isodesmic
values in Table 4). We thus decided to define our best estimates from the
average of reaction 1 and the isogyric reactions 2, 3 and LABEL:rxn:H2,
$\overline{\Delta_{f}H^{\ominus}}$($298.15\text{\,}\mathrm{K}$), with a
standard deviation $\mathrm{\Delta}$$\Delta_{f}H^{\ominus}$ that includes the
uncertainty on the average value as well as the experimental error bars listed
in Table 1.

For gaseous TBP, we predict a value of $-1281.7\pm
24.4\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, which is in line with
the value obtained from the group parameter correlation method,
$-1310\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$.Kertes (1972) As there
is no other data available for $\Delta_{f}H^{\ominus}$ of TBP, our data could
be considered as the most accurate one.

The $\Delta_{f}H^{\ominus}$ of TBP in the liquid state are obtained from data
in the gaseous state and its corresponding heat of solvation
($-85.6\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, see Table 3 in
Section II.4), the value $-1367.3\pm
24.4\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ is obtained. The
experimental data are scattered, thereby the heats of combustion of TBP
obtained by bomb calorimeter were reported in only three experimental studies.
Starostin et al. (1966); Kindle (1974); Erastov and Tarasov (1992) The derived
standard enthalpies of formation exhibit a large uncertainty over a range of
about $200\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. These
disagreements shall be assessed from the details given by the authors because
the calorimetry requires a very accurate analysis regarding correction factors
of each experimental step, therefore, a study lacking all details has to be
cautiously considered.

Kindle Kindle (1974) investigated the heat release of the TBP combustion for a
safety analysis of flammability. Initially, these data were not intended for
the evaluation of thermodynamic quantities. Nevertheless, an estimation is
possible by considering a complete combustion, as performed in a review of TBP
properties:Burger (1990) $-1250\pm
69\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. But this experiment is
tainted by unburnt TBP as well as soot formation. Then author does not mention
any analysis of the purity of the TBP sample and the post-combustion products
are not fully investigated. From these statements, the value extracted from
Kindle’s report is rejected even if the calorimeter calibration, as well as
explanation of each correction factor, have been perfectly performed.

Two other investigations Starostin et al. (1966); Erastov and Tarasov (1992)
used a bomb calorimeter with an isothermal shell. However, as opposed to
Starostin,Starostin et al. (1966) Erastov and TarasovErastov and Tarasov
(1992) has employed a rotating bomb. The main benefit of this facility is to
ensure a complete combustion process and a good thermal transfer between
sample zone and shell. To go further, the analyses carried out by Starostin
Starostin et al. (1966) do not show the hydrolysis of pyrophosphoric acid
certainly due to the use of static shell. In addition, as the both experiments
are similar, the two reported net heat released can be compared; Starostin has
measured the lowest heat release highlighting also the consequence of static
bomb. These findings suggest that Erastov and Tarasov’s measurements Erastov
and Tarasov (1992) are the most relevant. Otherwise, Erastov and Tarasov have
purified the TBP followed by a check of impurities content (moisture $<$
0.1%mol and no impurity). Then each sample is sealed in a container made of
polyester to prevent any reaction with the environment before combustion run.
Furthermore, a complete analysis at the end of runs have been performed in
order to fully examine the composition of by-products. To conclude, the
standard heat of formation extracted from Erastov and Tarasov’s investigations
is considered as the most accurate and equal to $-1382\pm
12\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$.

Finally, our theoretical standard enthalpy of formation is in agreement with
the selected experimental data. The discrepancy between them lies below
$14.7\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, therefore, validating
our quantum chemistry approach and its ability in predicting thermochemical
data. Thus, it gives confidence for applications to degraded products, HDBP(l)
and H2MBP(l).

###  III.2 Heat of formation of HDBP and H2MBP and investigation of the
hydrolysis reactions

Table 5:  Ab initio $\mathrm{LCCSD(T)(CBS^{\ast})}$ Gibbs reaction energies for the hydrolysis reactions H_3-nnBP(l) + H2O(l) -¿ H_4-n(n-1)BP(l) + C4H9OH(l) in the TBP and water solvents at $298.15\text{\,}\mathrm{K}$ and $363.15\text{\,}\mathrm{K}$ ($90\text{\,}\mathrm{\SIUnitSymbolCelsius}$). All values are in $\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$. Reaction |  $\Delta_{r}H^{\ominus}$($298.15\text{\,}\mathrm{K}$) |  $\Delta_{r}G^{\ominus}$($298.15\text{\,}\mathrm{K}$) |  $\Delta_{r}G^{\ominus}$($363.15\text{\,}\mathrm{K}$)
---|---|---|---
| TBP liquid
TBP(l) + H2O(sln) -¿ HDBP(sln) + C4H9OH(sln) | -27.4 | 0.4 | 5.7
HDBP(sln) + H2O(sln) -¿ H2MBP(sln) + C4H9OH(sln) | -32.6 | -2.6 | 3.1
H2MBP(sln) + H2O(sln) -¿ H3PO4(sln) + C4H9OH(sln) | -35.8 | -6.8 | -1.2
TBP(l) + 3 H2O(sln) -¿ H3PO4(sln) + 3 C4H9OH(sln) | -95.8 | -9.0 | 7.7
| water
TBP(aq) + H2O(l) -¿ HDBP(aq) + C4H9OH(aq) | -21.6 | 0.6 | 4.9
HDBP(aq) + H2O(l) -¿ H2MBP(aq) + C4H9OH(aq) | -23.7 | -1.8 | 2.5
H2MBP(aq) + H2O(l) -¿ H3PO4(aq) + C4H9OH(aq) | -22.7 | -3.8 | -0.1
TBP(aq) + 3 H2O(l) -¿ H3PO4(aq) + 3 C4H9OH(aq) | -68.0 | -5.0 | 7.2

The calculated values of standard enthalpy of formation for both degraded
products are reported in Table 4. In the liquid state, $\Delta_{f}H^{\ominus}$
of HDBP and H2MBP are $-1348.7\pm 19.6\text{\,}$ and $-1323.8\pm
14.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, respectively. Again,
there is no available data to compare these results.

To go further from these standard heats of formation, the hydrolysis phenomena
can be analyzed in both TBP and water solvents Table 5 summarizes the heat and
Gibbs energies of reaction for different temperatures. The hydrolysis process
(see reactions in Table 5) corresponds to the cleavage of P-O and O-H bonds
together with the formation of the same bonds regarding each of the three
intermediate steps. Consequently, the step-wise heats of reaction and the
overall one are expected to be nearly thermoneutral. Obtained heats of
reaction are weakly exothermic in both solvents which is consistent with the
later assumption. The trend is that the larger butyl phosphate, the lower
exothermic its degradation. Nevertheless, it is commonly admitted that TBP
slowly degrades by an acid catalyzed hydrolysis even at room temperature
Burger (1958) leading to a negative Gibbs energy of reaction. However, our
data sometimes highlight slightly positive value of Gibbs energy of reaction
meaning an unfavorable process in contradiction with the TBP knowledge. At
room temperature, whatever the solvent, there is only the first hydrolysis
step that is weakly endergonic but the value is smaller than the modeling
accuracy of COSMO-RS. The overall hydrolysis reaction remains a favorable
process. At temperature just below boiling point, the Gibbs reaction energies
become more endergonic except for hydrolysis of mono-butyl phosphate. These
disagreements can be explained by the fact that our theoretical model reaches
its accuracy limits to predict a weak thermal effect. In addition, the entropy
of such species is very sensitive data and requires a detailed vibrational
analysis of molecular motions in solution. Therefore its contribution to Gibbs
reaction energy may be responsible too. Thereby, the TBP, HDBP and H2MBP
hydrolysis should moderately contribute to the heat release during a thermal
runaway.

##  IV Conclusions

In this paper, we have established a methodology to predict the thermodynamic
parameters of TBP and its degradation products HDBP and H2MBP in both gas and
liquid phases. Basis set and solvent effects are considered and incorporated
into a composite thermochemistry scheme based on B3LYP/def2-TZVP geometries
and partition functions, and LCCSD(T) extrapolated to the complete basis set
limit for calculating accurate reaction enthalpies that are used to determine
heats of formations for the butyl-phosphate molecules. The excellent
similarities of the derived standard enthalpies of formation irrespective to
the chosen formation reaction, makes us propose that the uncertainty in the
herewith reported $\Delta_{f}H^{\ominus}$ is in the range of
$20\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ in the gas phase. In the
solvent, the hereafter quoted standard deviations have to be complemented by
the error induced by the COSMO-RS solvent model (about
$22\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ for TBP). This work makes
us propose the value $-1367.3\pm
24.4\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ for the heat of
formation ($298.15\text{\,}\mathrm{K}$) of TBP, confirming the value of
$-1382\pm 12\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ proposed by
Erastov and TarasovErastov and Tarasov (1992) In the absence of experimental
data for HDBP and H2MBP, this work predicts their heats of formation in the
liquid phase to be $-1348.7\pm 19.6$ and $-1323.8\pm
14.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, respectively. These
data allow us to predict that the complete hydrolysis of TBP is clearly
exothermic.

##  V Supplementary material

See supplementary material for 1) details on the estimate of the uncertainty
for the $\Delta_{f}H^{\ominus}$($298.15\text{\,}\mathrm{K}$) of H3PO4(l); 2)
Energies and enthalpic contributions for reactions 1, 2, 3 and LABEL:rxn:H2;
3) Enthalpies of formation derived from the isogyric reactions; 4) Molecular
thermodynamics data for the studied systems.

##  VI Acknowledgments

The members of the PhLAM laboratory acknowledge support from the CaPPA project
(Chemical and Physical Properties of the Atmosphere) that is funded by the
French National Research Agency (ANR) through the PIA (Programme
d’Investissement d’Avenir) under contract“ANR-11-LABX-0005-0” and by the
Regional Council “Hauts de France” and the ”European Funds for Regional
Economic Development” (FEDER) through the Contrat de Projets Etat-Région
(CPER) CLIMIBIO (Changement climatique, dynamique de l’atmosphère, impacts sur
la biodiversité et la santé humaine). We also acknowledge financial support
from Slovak grants APVV-15-0105 and VEGA 1/0465/15, the computational
resources of HPC Cluster of Slovak University of Technology (projects ITMS
26230120002 and 26210120002) and funding by the ERDF under the project
”University Scientific Park Campus MTF STU - CAMBO” ITMS: 26220220179.

## References

  * Colven et al. (1953) T. J. Colven, G. M. Nichols, and T. H. Siddall, Tech. Rep. DP-MS-83-142; CONF-840392-1; ON: DE84008185, Atomic energy division-Savannah River laboratory, Augusta, Ga (1953).
  * Durant (1983) W. S. Durant, Tech. Rep. DP-MS-83-142; CONF-840392-1; ON: DE84008185, Savannah River laboratory, Aiken, SC (USA) (1983).
  * IAE (1998) Tech. Rep., International Atomic Energy Agency, Vienna (1998).
  * Wright and Paviet-Hartmann (2010) A. Wright and P. Paviet-Hartmann, Sep. Sci. Technol. 45, 1753 (2010).
  * Burger (1990) L. L. Burger, in _Science and Technology of Tributyl Phosphate_ , edited by W. W. Schulz, L. L. Burger, and J. D. Navratil (CRC Press, London, 1990), vol. 3, chap. 3, pp. 25–67.
  * Starostin et al. (1966) A. D. Starostin, A. V. Nikolaev, and Y. A. Afanas’ev, Russ. Chem. Bull. 15, 1255 (1966).
  * Kindle (1974) C. H. Kindle, Tech. Rep. ARH–2862, Atlantic Richfield Hanford Co., Richland, Wash. (USA) (1974).
  * Erastov and Tarasov (1992) P. A. Erastov and A. V. Tarasov, Russ. J. Phys. Chem. 66, 2591 (1992).
  * Cox et al. (1984) J. Cox, D. Wagman, and V. Medvedev, _CODATA Key Values for Thermodynamics_ (Hemisphere Publishing Corp., 1984), p. 1.
  * Afeefy et al. (http://webbook.nist.gov) H. Afeefy, J. Liebman, and S. Stein, in _NIST Chemistry WebBook, NIST Standard Reference Database Number 69_ , edited by P. J. Linstrom and W. G. Mallard (National Institute of Standards and Technology, Gaithersburg MD, 20899, http://webbook.nist.gov).
  * Chase (1998) M. W. Chase, Jr., J. Phys. Chem. Ref. Data 25, 1 (1998).
  * Watts (2016) J. D. Watts, _Chemical Reactions: Thermochemical Calculations_ (Springer Netherlands, Dordrecht, 2016), pp. 1–18, ISBN 978-94-007-6169-8.
  * Hehre (2003) W. J. Hehre, _A guide to molecular mechanics and quantum chemical calculations._ (Wavefunction, Inc., Irvine, 2003).
  * Becke (1993) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
  * Lee et al. (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  * Vosko et al. (1980) S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
  * Stephens et al. (1994) P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Phys. Chem. 98, 11623 (1994).
  * Weigend et al. (1998) F. Weigend, M. Häser, H. Patzelt, and R. Ahlrichs, Chem. Phys. Lett. 294, 143 (1998), ISSN 0009-2614.
  * Weigend and Ahlrichs (2005) F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005).
  * Klamt and Schüürmann (1993) A. Klamt and G. Schüürmann, J. Chem. Soc. Perkin Trans. 2 pp. 799–805 (1993).
  * Neese and Valeev (2011) F. Neese and E. F. Valeev, J. Chem. Theory Comput. 7, 33 (2011).
  * Hättig and Weigend (2000) C. Hättig and F. Weigend, J. Chem. Phys. 113, 5154 (2000).
  * Hättig et al. (2006) C. Hättig, A. Hellweg, and A. Köhn, Phys. Chem. Chem. Phys. 8, 1159 (2006).
  * Werner and Schütz (2011) H.-J. Werner and M. Schütz, J. Chem. Phys. 135, 144116 (2011).
  * Feller (1992) D. Feller, J. Chem. Phys. 96, 6104 (1992).
  * Feller (1993) D. Feller, J. Chem. Phys. 98, 7059 (1993).
  * Helgaker et al. (1997) T. Helgaker, W. Klopper, H. Koch, and J. Noga, J. Chem. Phys. 106, 9639 (1997).
  * Schwilk et al. (2015) M. Schwilk, D. Usvyat, and H.-J. Werner, J. Chem. Phys. 142, 121102 (2015).
  * Werner et al. (2015) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona, R. Lindh, et al., _Molpro, version 2015.1, a package of ab initio programs_ (2015), see http://www.molpro.net.
  * Werner et al. (2012) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, WIREs Comput. Mol. Sci. 2, 242 (2012), ISSN 1759-0884.
  * (31) _TURBOMOLE V7.1 2016, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com._
  * Hansen et al. (2014) A. Hansen, C. Bannwarth, S. Grimme, P. Petrović, C. Werlé, and J.-P. Djukic, ChemistryOpen 3, 177 (2014), ISSN 2191-1363.
  * Eckert and Klamt (2002) F. Eckert and A. Klamt, Angew. Chem. Int. Ed 48, 369 (2002), ISSN 1547-5905.
  * Diedenhofen et al. (2007) M. Diedenhofen, A. Klamt, K. Marsh, and A. Schäfer, Phys. Chem. Chem. Phys. 9, 4653 (2007).
  * Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  * Perdew (1986a) J. P. Perdew, Phys. Rev. B 33, 8822 (1986a).
  * Perdew (1986b) J. P. Perdew, Phys. Rev. B 34, 7406 (1986b).
  * Schäfer et al. (1994) A. Schäfer, C. Huber, and R. Ahlrichs, J. Chem. Phys. 100, 5829 (1994).
  * (39) F. Eckert and A. Klamt, _COSMOtherm, Version C3.0, Release 17.01, COSMOlogic GmbH & Co. KG, Leverkusen, Germany, 2015._
  * Chickos et al. (1995) J. S. Chickos, S. Hosseini, and D. G. Hesse, Thermochim. Acta 249, 41 (1995).
  * Marsh (1987) K. N. Marsh, ed., _Recommended Reference Materials for the Realization of Physicochemical Properties_ (Boston : Blackwell Scientific Publications, Oxford, 1987).
  * Schröder and Coutinho (2014) B. Schröder and J. A. P. Coutinho, Fluid Phase Equilibr. 370, 24 (2014), ISSN 0378-3812.
  * Alexeev et al. (2005) Y. Alexeev, T. L. Windus, C.-G. Zhan, and D. A. Dixon, Int. J. Quantum Chem. 102, 775 (2005).
  * Kertes (1972) A. Kertes, J. Inorg. Nucl. Chem. 34, 796 (1972).
  * Burger (1958) L. L. Burger, Tech. Rep. HW-40910, General Electric Co. Hanford Atomic Products Operation, Richland, Wash. (USA) (1958).