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

TL;DR

This study uses advanced quantum chemical methods to accurately predict the formation enthalpies of TBP and its degradation products, providing crucial thermodynamic data for nuclear safety assessments.

Contribution

It offers the first reliable, high-accuracy thermodynamic data for TBP and its degradation products using state-of-the-art ab initio calculations and solvation models.

Findings

Final uncertainties around 20 kJ/mol for formation enthalpies.

Predicted thermodynamic values for TBP and derivatives in gas and liquid phases.

Hydrolysis of TBP is nearly thermoneutral.

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-(H$_2$MBP}) 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-(H$_2$MBP) 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 kJ mol$^{-1}$. We thus propose with final uncertainty of about 20 kJ mol$^{-1}$ standard enthalpies of formation of TBP, HDBP, and H$_2$MBP which amounts to -1281.7$\pm$24.4, -1229.4$\pm$19.6 and -1176.7$\pm$14.8 kJ 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 kJ mol$^{-1}$, to which we may add about -22 kJ mol$^{-1}$ error from the COSMO-RS solvent model. From these data, we predict the complete hydrolysis of TBP to be nearly thermoneutral.

Figures

Tables

Table 1. Table 1: Known Gas/Liquid Phase Standard Enthalpies of Formation Δ f ​ H ⊖ subscript Δ 𝑓 superscript 𝐻 symmetric-difference \Delta_{f}H^{\ominus} ( 298.15   K times 298.15 kelvin 298.15\text{\,}\mathrm{K} ) in kJ   mol − 1 times kilojoule mole 1 \mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}

Molecule	${\Delta }_f{H}^{\ominus }$($298.15K$)	Reference
\ceHF(g)	$-273.3\pm 0.7$	9
\ceHCl(g)	$-92.31\pm 0.10$	9
\ceH2O(g)	$-241.826\pm 0.040$	9
\ceC4H9OH(g)	$-277\pm 5.$	10
\cePH3(g)	$5.47\pm 1.7$	11
\cePOF3(g)	$-1254.25\pm 8$	11
\cePOCl3(g)	$-559.82\pm 1.7$	11
\ceH3PO4(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

Table 2. Table 2: Reaction energies Δ r ​ E subscript Δ 𝑟 𝐸 \Delta_{r}E in   kJ   mol − 1 times absent times kilojoule mole 1 \text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1} for \ce 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

Table 3. Table 3: COSMO-RS vaporization enthalpies of all relevant species in   kJ   mol − 1 times absent times kilojoule mole 1 \text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1} , together with experimental values when available.

Molecule	liquid phase	${\Delta }_{vap}{H}^{\ominus }$($298.15K$)
		COSMO-RS	Exp.	Ref.
\ceC4H9OH	\ceC4H9OH	58.8	$52.34\pm 0.02$	40
\ceH2O	\ceH2O	49.2	$43.99\pm 0.07$	41
\ceTBP	\ceTBP	85.6	64.4	5
\ceHDBP	\ceTBP	118.5
	\ceHDBP	119.3
\ceH2MBP	\ceTBP	159.5
	\ceH2MBP	147.1
\ceH3PO4	\ceH3PO4	144.7

Table 4. Table 4: Δ f ​ H ⊖ subscript Δ 𝑓 superscript 𝐻 symmetric-difference \Delta_{f}H^{\ominus} ( 298.15   K times 298.15 kelvin 298.15\text{\,}\mathrm{K} ) of TBP, HDBP, and \ce 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 . Δ f ​ H ⊖ ¯ ¯ subscript Δ 𝑓 superscript 𝐻 symmetric-difference \overline{\Delta_{f}H^{\ominus}} ( 298.15   K times 298.15 kelvin 298.15\text{\,}\mathrm{K} ) represent the average between these two data sets. The standard deviations Δ Δ \mathrm{\Delta} Δ f ​ H ⊖ subscript Δ 𝑓 superscript 𝐻 symmetric-difference \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 Δ f ​ H ⊖ subscript Δ 𝑓 superscript 𝐻 symmetric-difference \Delta_{f}H^{\ominus} ( 298.15   K times 298.15 kelvin 298.15\text{\,}\mathrm{K} )= Δ f ​ H ⊖ ¯ ¯ subscript Δ 𝑓 superscript 𝐻 symmetric-difference \overline{\Delta_{f}H^{\ominus}} ( 298.15   K times 298.15 kelvin 298.15\text{\,}\mathrm{K} )- Δ v ​ a ​ p ​ H ⊖ subscript Δ 𝑣 𝑎 𝑝 superscript 𝐻 symmetric-difference \Delta_{vap}H^{\ominus} ( 298.15   K times 298.15 kelvin 298.15\text{\,}\mathrm{K} ). All results are in kJ   mol − 1 times kilojoule mole 1 \mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1} .

Species	Reaction	${\Delta }_f{H}^{\ominus }$($298.15K$)	$\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.15K$) =	-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.15K$) =	-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.15K$) =	-1176.7	14.8
H2MBP(l)		-1323.8

Table 5. Table 5: Ab initio LCCSD ​ ( T ) ​ ( CBS ∗ ) LCCSD T superscript CBS ∗ \mathrm{LCCSD(T)(CBS^{\ast})} Gibbs reaction energies for the hydrolysis reactions \ce H_3-nnBP(l) + H2O(l) -¿ H_4-n(n-1)BP(l) + C4H9OH(l) in the TBP and water solvents at 298.15   K times 298.15 kelvin 298.15\text{\,}\mathrm{K} and 363.15   K times 363.15 kelvin 363.15\text{\,}\mathrm{K} ( 90   °C times 90 celsius 90\text{\,}\mathrm{\SIUnitSymbolCelsius} ). All values are in kJ   mol − 1 times kilojoule mole 1 \mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1} .

Reaction	${\Delta }_r{H}^{\ominus }$($298.15K$)	${\Delta }_r{G}^{\ominus }$($298.15K$)	${\Delta }_r{G}^{\ominus }$($363.15K$)
	TBP liquid
\ceTBP(l) + H2O(sln) -¿ HDBP(sln) + C4H9OH(sln)	-27.4	0.4	5.7
\ceHDBP(sln) + H2O(sln) -¿ H2MBP(sln) + C4H9OH(sln)	-32.6	-2.6	3.1
\ceH2MBP(sln) + H2O(sln) -¿ H3PO4(sln) + C4H9OH(sln)	-35.8	-6.8	-1.2
\ceTBP(l) + 3 H2O(sln) -¿ H3PO4(sln) + 3 C4H9OH(sln)	-95.8	-9.0	7.7
	water
\ceTBP(aq) + H2O(l) -¿ HDBP(aq) + C4H9OH(aq)	-21.6	0.6	4.9
\ceHDBP(aq) + H2O(l) -¿ H2MBP(aq) + C4H9OH(aq)	-23.7	-1.8	2.5
\ceH2MBP(aq) + H2O(l) -¿ H3PO4(aq) + C4H9OH(aq)	-22.7	-3.8	-0.1
\ceTBP(aq) + 3 H2O(l) -¿ H3PO4(aq) + 3 C4H9OH(aq)	-68.0	-5.0	7.2

Equations

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

$$E_{\text{corr}}(n)=E_{\text{corr}}^{\text{CBS}}+A{n}^{-3}.$$

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

Full text

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 \ceUO2(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 (\ceH2MBP), and ultimately phosphoric acid (\ceH3PO4). The reaction by-products (butanol and butyl nitrate) may then undergo oxidation to form carboxylic acids and gases. In a similar way, \ceUO2(NO3)2(TBP)2 may be decomposed to \ceUO2(NO3)2TBP(HDBP) or \ceUO2(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 \ceH2MBP.

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 \ceH_3-nnBP

The standard enthalpies of formation in the gas and TBP liquid phases of the \ceH_3-nnBP molecules, \ceH2MBP (n=1), \ceHDBP (n=2), and \ceTBP (n=3) are derived from computed enthalpies of reactions.

[TABLE]

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

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 \ceH_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 \ceH2MBP, \ceHDBP and \ceTBP are displayed on Figure 1.

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 \ceH2MBP. 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)

[TABLE]

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

[TABLE]

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})}}$:

[TABLE]

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

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

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 \ceP=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 \ceH2MBP 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 \cePH3 (three single \ceP-H bonds) and in the product \ceH2MBP (one double \ceP=O bond, three single \ceP-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”.

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 \ceH3PO4, 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 \ceH3PO4 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 \ceH2MBP, \ceHDMP and \ceTBP in the TBP organic phase.

III Results and Discussion

First, the standard enthalpies of formation of \ceH2MBP, 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 \ceHDBP and \ceH2MBP (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, \ceHDBP(l) and \ceH2MBP(l).

III.2 Heat of formation of HDBP and \ceH2MBP and investigation of the hydrolysis reactions

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 \ceH2MBP 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 \ceP-O and \ceO-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 \ceH2MBP 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 \ceTBP and its degradation products \ceHDBP and \ceH2MBP 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 \ceTBP). 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 \ceTBP, 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 \ceHDBP and \ceH2MBP, 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 \ceTBP 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 \ceH3PO4(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.