Theoretical study of deuterium isotope effects on acid–base equilibria under ambient and hydrothermal conditions

Abstract

Quantum electronic structure methods are applied for the first time to the study of deuterium isotope effects (DIE) on pK_a values under ambient (25 °C, 101.3 kPa) and hydrothermal (250 °C, 20.0 MPa) conditions. This work focuses on sixteen organic acids and explores several methodologies for calculating pK_a values and various pK_a differences in H₂O and D₂O under two sets of conditions. Two functionals are considered (B3LYP and BLYP) and solvent effects are accounted for by means of continuum solvation methods (PCM, CPCM, Onsager and SMD). Excellent agreement with experiment is obtained for the calculated DIE (ΔpK_a = pK_a(D₂O) − pK_a(H₂O)) at the B3LYP-PCM/6-311++G(d,p) level of theory for the two sets of conditions. These values, which are almost constant for a given set of temperature and pressure conditions, are determined by the difference between the Gibbs free energies of formation of the acid and its deuterated form in each solvent. However, accurate predictions under ambient conditions can also be made from zero-point energy differences. The average calculated ΔpK_a values under ambient (experimental average: 0.53) and hydrothermal conditions were 0.65 and 0.37, respectively. The mean absolute error between calculated and experimental ΔpK_a values under ambient conditions was 0.11. The methodology applied is a very important tool for accurately predicting DIE on pK_a values under both ambient and hydrothermal conditions, which can be used to make accurate pK_a predictions in D₂O.

---

1. Introduction

Heavy water (D₂O) has significant practical applications as it acts as a neutron moderator in the reactor core and as a heat-transfer fluid in the pressure tubes of the Canadian-designed CANDU® (Canadian Deuterium Uranium) nuclear reactors, which operate at temperatures from 250 to 300 °C (523.15 K to 573.15 K) and pressures of ∼10 MPa. Lithium hydroxide and hydrogen are added to the heavy water coolant as pH and redox control agents to minimize corrosion, fuel deposits, and corrosion product transport; however, these chemical treatments are largely based on high-temperature light-water experimental data, and 40 years of reactor operating experience. The operators of CANDU reactors have expressed interest in reducing the operating pH to mitigate the effects of flow-accelerated corrosion of outlet feeder pipes and have identified a need for a more thorough understanding of deuterium isotope effects on acid–base equilibria in high-temperature water. At room temperature, the dissociation constants (K_a, pK_a = −log K_a) of inorganic and organic acids are known to be greater in water than in heavy water, ΔpK_a = pK_a(D₂O) − pK_a(H₂O), 0.86 ∀ 0.23.¹ However, to date only a handful of quantitative studies to measure acid–base equilibrium constants above 100 °C have been reported in the literature.^1b,2–7

The determination of pK_a values at high temperatures and pressures in light and heavy water is a challenge experimentally because of the highly specialized equipment needed. As a result, there is much interest in developing computational predictive tools. Ab initio calculations for determining accurate pK_a values in water are also quite challenging, largely because of the difficulty of treating solvation effects;^8–10 the most widely applied and pragmatic approach is to use continuum solvation methods.¹¹ Sometimes, explicit solvent molecules are used in combination with these methods.^8a,12 Ab initio molecular dynamic simulations have also been used to determine pK_a values in solution.¹³ However, this alternative technique is also based on approximations and it is significantly more complex and expensive than continuum methods from a computational point of view.

This initial study focuses on sixteen organic acids: acetic acid, four thermally-stable colorimetric pH indicators (β-naphthol, protonated s-collidine, protonated acridine, and β-naphthoic acid), and eleven phenols. Their structures are shown in Fig. 1. The first five molecules are used in laboratory studies to examine deuterium isotope effects on reactor chemistry and corrosion product transport under reactor operating conditions. At 25 °C, the substituted phenols included in this study cover the range 3 < pK_a < 10, depending on the nature of the substituents. Many of these compounds are thermally stable, and have UV-visible spectra, so they are candidates for use as high-temperature pH indicators. Two of them, 2-nitrophenol and 4-nitrophenol, have been used in this application in light water.^6a

Fig. 1 Molecules studied.

This study aims to determine whether computational methods frequently used for modelling the dissociation of organic acids in light water at room temperature can be extended to heavy water and/or high temperatures and pressures. We would also like to investigate whether these methods can be used to predict the temperature-dependence of pK_a values in light and heavy water, and deuterium isotope effects (ΔpK_a = pK_a(D₂O) − pK_a(H₂O)) under different sets of temperature/pressure conditions. Furthermore, we are interested in investigating theoretically how heavy water affects acid–base equilibria under ambient and hydrothermal conditions, as a means of predicting the deuterium isotope effect (DIE) under nuclear reactor operating conditions. In the research reported below, we have examined the success of several computational methodologies in accurately reproducing the available experimental pK_a values in H₂O and D₂O for the acids mentioned above under both ambient [25 °C (298.15 K), 101.3 kPa], and hydrothermal conditions [250 °C (523.15 K), 20.0 MPa]. While DIE on pK_a values has been widely examined under ambient conditions, to the best of our knowledge, no previous theoretical studies of this topic have been carried out under ambient and elevated temperatures and pressures applying electronic structure methods.

2. Methodology

2.1. Computational details

Electronic structure calculations were carried out with the Gaussian03 software package.¹⁴ Initial preliminary calculations that focused on five acids (β-naphthol, protonated s-collidine, protonated acridine, β-naphthoic acid, and acetic acid) were performed using the B3LYP and BLYP functionals with the 6-311++G(d,p) basis set (see Appendix A in the ESI† section). All stationary points were characterized as minima by a vibrational frequency analysis using analytical second derivatives. Continuum solvent models (IEF-PCM,¹⁵ CPCM,¹⁶ Onsager¹⁷ and SMD¹⁸) were used to account for solvent effects through both single-point energy calculations and in geometry optimizations and frequency calculations. The IEF-PCM and CPCM calculations employ UAHF atomic radii when constructing the solvent cavity. Volumes for the Onsager calculations were determined from gas-phase geometries at the same level of theory that such geometries were optimized (functional/6-311++G(d,p)).

Water at 25 °C and 101.3 kPa (1 atm) is an explicitly defined solvent in Gaussian03. As displayed in Table 1, solvent parameters (dielectric constants, molar volumes, and numeral densities)¹⁹ for heavy water (and for both solvents at 250 °C and 20.0 MPa) were required to properly define the solvent as a continuum when applying the IEF-PCM and CPCM methods. Values for the density and dielectric constant of light water were calculated from the equations of state reported by Wagner and Pruss^19a and Fernandez et al.,^19b respectively. The density of heavy water was taken from Hill's equation of state^19c using software distributed by NIST.^19d,e The dielectric constant of heavy water was calculated from light water values, using the method reported by Trevani et al.^19f For Onsager and SMD calculations only the dielectric constants are needed. When calculations were done on the acids in heavy water, the acidic hydrogen atom in each compound was replaced by deuterium. Thermal corrections were calculated taking the desired isotopes, temperature and pressure conditions into account. Based on the preliminary results obtained, the eleven phenols were calculated at the B3LYP-IEF-PCM/6-311++G(d,p) level of theory.

Table 1 Solvent parameters used taken from ref. 19^a

		25 °C, 101.3 kPa		250 °C, 20.0 MPa
		H2O	D2O	H2O	D2O
a vmol = NA/rho, where NA is Avogadro's number; g03 reports vmol in units of Å^3 but the actual units are cm^3 mol^−1. Temperature is another parameter to be defined when different from 298.15 K.
Dielectric constant	eps	79.14	78.95	27.87	27.75
Molar volume^a (cm^3 mol^−1)	vmol	17.90	17.97	22.06	22.17
Numeral density (Å^−3)	rho	0.03365	0.03352	0.02730	0.02716

---

2.2. Equilibria and equations used for calculating pK_a values in solution

pK values can be calculated using eqn (1), where R is the ideal gas constant and T is the temperature in Kelvin. The standard Gibbs free energy change (ΔG°) is calculated as the difference between the sum of the standard Gibbs free energies of formation (Δ_fG°) of the products, and of the reactants (stoichiometric coefficients are assumed to be 1), as shown in eqn (2).

pK = ΔG°/RT ln(10) (1)

ΔG° = ∑Δ_fG°(products) − ∑Δ_fG°(reactants) (2)

The thermodynamic pK_a of an acid (HA) is the experimental value of the equilibrium quotient extrapolated at ionic strength zero. For the acid ionization reaction, shown in Scheme (1), the thermodynamic equilibrium constant is shown in eqn (3). For convenience in comparing equilibrium constants in H₂O to those in D₂O, experimentalists often follow the practice recommended by Laughton and Robertson,^1b and define molalities in both light and heavy water as moles of solute per 55.509 moles solvent, so that m° = 1 mol/55.509 mol solvent. This is the so-called “aquamolal” standard state.^1b

HA ⇆ A⁻ + H⁺ (S1)

From a theoretical point of view, the calculation of the dissociation constant for Scheme (1) requires the experimental Δ_fG° values for H⁺ in either the gas phase or aqueous solution.²⁰ This approach cannot be applied to determine the pK_a of a deuterated acid (DA) in heavy water because the required experimental data for D⁺ are not available. An acid is said to be deuterated in this study when its acidic H atom has been replaced by D; this is the case when the dissociation of the acid is considered in heavy water.

An alternative equilibrium (Schemes (2a) and (2b)) can be considered for determining pK_a values in H₂O and D₂O. If the pK of Schemes (2a) and (2b) is denoted pK₂, the thermodynamic pK_a,c value (relative to the hypothetical 1 molar standard state (1 mol L⁻¹)) of HA in H₂O and of DA in D₂O can be calculated using eqn (4a) and (4b), respectively, where ρ is the density (in kg L⁻¹) and M is the molar mass (in kg mol⁻¹).

HA_(H2O) + H₂O_(H2O) ⇆ A_(H2O)⁻ + H₃O_(H2O)⁺ (S2a)

DA_(D2O) + D₂O_(D2O) ⇆ A_(D2O)⁻ + D₃O_(D2O)⁺ (S2b)

pK_a,c(H2O) = pK_2(H2O) + pK_a(H3O⁺) = pK_2(H2O) − log[H₂O] = pK_2(H2O) − log(ρ_H2O/M_H2O) (4a)

pK_a,c(D2O) = pK_2(D2O) + pK_a(D3O⁺) = pK_2(D2O) − log[D₂O] = pK_2(D2O) − log(ρ_D2O/M_D2O) (4b)

In principle, experimental Δ_fG° values in aqueous solution for H₂O and H₃O⁺ could be used in Scheme (2a), but these values are not available for D₂O and D₃O⁺. Hence, the pK of these equilibria will be determined using the calculated Δ_fG° of the four species involved in each equilibrium.

Care must be taken when comparing solvation thermodynamics in different solvents. Although standard practice under ambient conditions¹ is to compare Henry's Law standard states defined at the same solvent volume (1 molar standard state) or same mole fraction, the thermodynamics of hydrothermal solutions is based on the hypothetical 1 molal standard state [1 mol kg⁻¹ = 1 mol/(55.509 mol H₂O)]. Our values from ab initio calculations are expressed in terms of the hypothetical 1 molar standard state, pK_a,c. Following the recommendations of Laughton and Robertson,^1b we have chosen to express these calculated values, and the experimental values from the literature, relative to the hypothetical 1 aquamolal standard state mentioned earlier. In H₂O, aquamolality = molality, thus K_a,aq = K_a,m = K_a,c(H2O)/ρ_H2O. However, in D₂O, aquamolality = molality × 1.1117 (number which is obtained from the molar mass of D₂O multiplied by 55.509 mol (1 kg) of solvent), thus K_a,aq = 1.1117K_a,m = 1.1117K_a,c(D2O)/ρ_D2O. From this, eqn (4a) and (4b) become eqn (5a) and (5b), respectively.

pK_a,aq = pK_a,c(H2O) + log ρ_H2O = pK_2(H2O) + log M_H2O (5a)

pK_a,aq = pK_a,c(D2O) + log(ρ_D2O/1.1117) = pK_2(D2O) + log(M_D2O/1.1117) (5b)

The final equation, used for calculating the thermodynamic pK_a,aq value of the acids (deuterated acids) under study relative to the hypothetical 1 aquamolal standard state in water (heavy water) using Schemes (2a) and (2b) at the two sets of conditions considered, is eqn (6a) and (6b). In this equation, ΔG^°_2a (ΔG^°_2b) refers to the standard Gibbs free energy change according to Schemes (2a) and (2b), determined with eqn (7a) and (7b) using calculated Δ_fG° values in water (heavy water).

pK_a,aq = ΔG^°_2a/RT ln(10) + log M_H2O (6a)

pK_a,aq = ΔG^°_2b/RT ln(10) + log(M_D2O/1.1117) (6b)

2.3. Calculated Δ_fG° values in solution

Depending on how solvent effects are accounted for, several Δ_fG° values (denoted G for simplicity) in solution can be considered when working with eqn (7a) and (7b). The simplest approach is to consider solvent effects on a single-point energy calculation using a gas-phase geometry at the functional-solvent method/6-311++G(d,p)//functional/6-311++G(d,p) level of theory. Three G values in solution were considered this way. One of them, labelled G_1w (see eqn (8)), combines the uncorrected energy in solution, E_w, with the gas-phase thermal correction to the Gibbs free energy, TCG_gas. G_2w (see eqn (9)) is calculated combining the gas-phase G value, G_gas, with the Gibbs free energy of solvation, ΔG_solv. When using the SMD solvation model, ΔG_solv is calculated by subtracting the energy in solution, E_w, and the gas-phase energy, E_gas, (as indicated in the User's Guide of the software used) which makes the calculation of G_1w equivalent to that of G_2w. Therefore, for calculations using the SMD (or Onsager) method, only G_1w is reported.

G_1w = E_w + TCG_gas (8)

G_2w = G_gas + ΔG_solv (9)

The recommended radii for calculations that request the determination of ΔG_solv values is UAHF and these radii were optimized at the HF/6-31G(d) level of theory. Hence, ΔG_solv values calculated in water at 25 °C at this level of theory are expected to be better than when calculated at the same level of theory at which the gas-phase geometries are obtained. G values (labelled G_3w, see eqn (10)) calculated by combining the previously mentioned G_gas values and ΔG_solv values calculated at the HF/6-31G(d)-solvent method//functional/6-311++G(d,p) level of theory, ΔG_solv(HF), are also considered. G_2w and G_3w values are calculated with IEF-PCM and CPCM. Solvent effects can also be considered in both geometry optimizations and frequency calculations. Δ_fG° values in solution obtained this way were labelled G_w.

G_3w = G_gas + ΔG_solv(HF) (10)

Since in the reaction scheme used (Schemes (2a) and (2b)) there are the same number of reactant and product species, there is no need to explicitly change the reference state of the calculated G values from 101.3 kPa (1 atm) to 1 mol L⁻¹. Eqn (6a) and (6b) can be directly applied to the Δ_fG° values of each species calculated in solution (G_1w, G_2w, G_3w and G_w), in water and heavy water, respectively.

The data needed for the calculations previously described for the sixteen acids studied at the B3LYP-IEF-PCM/6-311++G(d,p) level of theory appear in Table S1 of the ESI section.† The raw data used for the preliminary calculations that focused on five of the acids appear in Tables A1 to A3 of Appendix A in the ESI† section.

3. Results and discussion

Initial preliminary calculations that focused on five of the sixteen acids considered in this study were performed using the B3LYP and BLYP functionals with the 6-311++G(d,p) basis set. Four continuum solvation methods were applied (IEF-PCM, CPCM, Onsager and SMD). Tables with the results obtained and their discussion are included in Appendix A in the ESI section.† The best results overall were obtained at the B3LYP/6-311++G(d,p) level of theory applying the IEF-PCM solvation method. Hence, this was the level of theory chosen to study the dissociation of a larger set of acids in both light and heavy water under ambient and hydrothermal conditions.

3.1. Comparison of the calculated pK_a values with experiment

The available experimental pK_a values in H₂O and D₂O under ambient conditions for the molecules under study in various standard states are shown in Table 2. These values, which have been mainly taken from ref. 1b, have been converted to the hypothetical 1 aquamolal standard state. Experimental values without a clear indication of the reference state used to report them were not considered in this study. In the cases where more than one experimental value was found for a given acid, averaged values were calculated and reported in Table 3, together with the available experimental pK_a values under hydrothermal conditions.

Table 2 Experimental pK_a values under ambient conditions converted to the aquamolal standard state^a

Acid	Values as reported				Aquamolal^c
	pKa(H2O)	pKa(D2O)	Standard state^b	Ref.	pKa(H2O)	pKa(D2O)
a Values at hydrothermal conditions already follow the aquamolal standard state and have been excluded from this table for simplicity. Values reported (mostly from ref. 1b) without a clear indication of the reference state used were not taken into account. When more than one pKa value is listed for a given acid, average values were calculated and reported in Table 3.b Abbreviations used: M, molarity; m, molality; aq, aquamolality.c Equations used for the reference state conversions: in H2O, pKa,aq = pKa,m = pKa,c(H2O) + log ρH2O; in D2O, pKa,aq = pKa,m(D2O) − log 1.1117 = pKa,c(D2O) + log(ρD2O/1.1117).
Phenol	10.00	10.62	M	7a	10.00	10.62
3-Methoxyphenol	9.62	10.20	M	7a	9.62	10.20
4-Methoxyphenol	10.24	10.85	M	7a	10.24	10.85
4-Bromophenol	9.35	9.94	M	7a	9.35	9.94
2-Nitrophenol	7.25	7.82	m	7b	7.25	7.77
	7.23	7.81	M	7c	7.23	7.81
4-Nitrophenol	7.22	7.77	M	7a	7.22	7.77
	7.24	7.80	m	7b	7.24	7.76
2,4-Dinitrophenol	4.06	4.55	M	7a	4.06	4.55
	4.07	4.59	M	7d	4.07	4.59
2,5-Dinitrophenol	5.19	5.70	M	7a	5.19	5.70
	5.20	5.73	M	7d	5.20	5.73
	5.17	5.67	m	7b	5.17	5.62
2,6-Dinitrophenol	3.73	4.22	M	7d	3.73	4.22
3,5-Dinitrophenol	7.31	7.92	m	7e	7.31	7.87
	6.70	7.31	m	7f	6.70	7.27
4-Chloro-2,6-dinitrophenol	2.96	3.45	M	7d	2.96	3.45
	2.97	3.48	m	7f	2.97	3.44
β-Naphthol	9.47	10.06	m	7a	9.47	10.01
	9.63	10.17	aq	4	9.63	10.17
β-Naphthoic acid	4.21	4.68	m	7a	4.21	4.63
Acetic acid	4.76	5.31	m	7f	4.76	5.27
	4.74	5.23	aq	5	4.74	5.23
s-Collidine	7.43		aq	7g	7.43
Acridine	5.58		aq	21a	5.58

---

Table 3 Experimental^a and calculated^b pK_a values of the acids studied in H₂O and D₂O

Acid	25 °C, 101.3 kPa				250 °C, 20.0 MPa
	Exp. pKa(H2O)	Calc. pKa(H2O)	Exp. pKa(D2O)	Calc. pKa(D2O)	Exp. pKa(H2O)	Calc. pKa(H2O)	Exp. pKa(D2O)	Calc. pKa(D2O)
a The reference state is 1 aquamolal; some experimental values at ambient conditions are calculated averages when more than one value has been reported (detailed in Table 2) [values in brackets are estimates using calculated DIE values from Table 4: unknown pKa(D2O) = exp. pKa(H2O) + calc. DIE(T,p)].b Calculations at the B3LYP-IEF-PCM/6-311++G(d,p) level of theory (i.e., using Gw values).c Mean absolute error, excluding estimated values in brackets.d Ref. 6a.e Ref. 4.f Ref. 21b.g Ref. 5.h Ref. 7g.i Ref. 21a.
Phenol	10.00	12.89	10.62	13.52		8.31		8.68
3-Methoxyphenol	9.62	12.73	10.20	13.36		8.23		8.60
4-Methoxyphenol	10.24	13.84	10.85	14.46		8.89		9.26
4-Bromophenol	9.35	11.35	9.94	11.99		7.33		7.71
2-Nitrophenol	7.24	7.02	7.79	7.70	6.85^d	4.69	[7.23]	5.07
4-Nitrophenol	7.23	5.80	7.76	6.43	6.57^d	4.04	[6.93]	4.41
2,4-Dinitrophenol	4.07	1.12	4.57	1.80		1.27		1.65
2,5-Dinitrophenol	5.19	4.00	5.69	4.69		2.50		2.88
2,6-Dinitrophenol	3.73	2.65	4.22	3.32		1.13		1.51
3,5-Dinitrophenol	7.00	5.55	7.57	6.17		3.75		4.11
4-Chloro-2,6-dinitrophenol	2.97	1.24	3.44	1.91		2.30		2.70
β-Naphthol	9.55	12.26	10.09	12.90	8.97^e	7.90	9.36^e	8.27
β-Naphthoic acid	4.21	5.95	4.63	6.56	5.98^f	4.28	[6.34]	4.64
Acetic acid	4.75	6.13	5.25	6.75	6.00^g	4.30	6.44^g	4.66
s-Collidine	7.43	5.66	[8.14]	6.36	4.26^h	2.91	[4.64]	3.30
Acridine	5.58	3.67	[6.26]	4.35	3.41^i	1.86	[3.79]	2.24
MAE^c		1.95		1.93		1.72		1.44

---

The aqueous pK_a values of acridine and β-naphthoic acid at 250 °C have been obtained using the equations reported in ref. 21 which predict the pK_a values of these compounds at any temperature. The pK_a values of β-naphthol and acetic acid were determined using UV-visible spectroscopy with a high-pressure platinum flow cell⁴ and AC conductance techniques,⁵ respectively, up to 300 °C in both solvents. Aqueous pK_a values for s-collidine^7g and for 2- and 4-nitrophenol^6a at 250 °C have also been reported.

The calculated pK_a values (using the four types of Δ_fG° values indicated in Section 2.3: G_1w, G_2w, G_3w, and G_w) and their errors, expressed as mean absolute error (MAE) and average error (AE), in H₂O and D₂O at 25 °C, 101.3 kPa and at 250 °C, 20.0 MPa, are displayed in Tables S2–S5 and Fig. S1–S4 of the ESI† section.

Under ambient temperature and pressure conditions (see Tables S2 and S3, Fig. S1 and S2†), the best results are obtained using G_w values (MAE = 1.95 (H₂O) and 1.93 (D₂O)), i.e., when solvent effects are accounted for in geometry optimizations and frequency calculations (B3LYP-IEF-PCM/6-311++G(d,p)). These values are displayed in Table 3 and Fig. 2 together with the corresponding experimental values in H₂O and D₂O. Most of the calculated values have errors greater than 1 pK_a unit in both solvents. Calculations using G_3w values (MAE = 2.18 (H₂O) and 1.77 (D₂O)), with solvent effects considered only on the calculated energies, produce slightly similar errors. The other approaches (using G_1w and G_2w values) produce pK_a values with significantly larger errors. In any case, the smallest errors obtained are still too large to give the methods applied any useful predictive capability.

Fig. 2 Experimental and calculated (G_w) pK_a values under ambient conditions in H₂O and D₂O.

There are very few experimental values to compare with under hydrothermal conditions (see Tables 2, S4 and S5, Fig. S3 and S4†); hence, it is difficult to make relevant generalizations. The errors are slightly reduced relative to the pK_a calculations under ambient conditions using G_w values (MAE = 1.72 (H₂O) and 1.44 (D₂O), see Table 2). However, the calculations using G_1w (MAE = 1.05 (H₂O) and 0.68 (D₂O)) and, in particular, G_2w (MAE = 0.71 (H₂O) and 0.61 (D₂O)), give the smallest errors. Four of the seven calculated pK_a values in H₂O using G_2w have errors equal to or greater than 1 pK_a unit. Hence, based on this information, it seems that the methods applied are not adequate for directly predicting accurate pK_a values under hydrothermal conditions.

3.2. Comparison of calculated pK_a differences

It is of interest to explore the accuracy of the calculation of two types of pK_a differences. One of them reflects the temperature-dependence of pK_a values in a given solvent, while the other difference reflects the deuterium isotope effect (DIE) on pK_a values under ambient and hydrothermal conditions. Exploring the accuracy of both types of pK_a differences making use of continuum solvation methods are additional objectives of this study.

3.2.1. The temperature-dependence of pK_a (pK_a(25 °C) − pK_a(250 °C). The effects of temperature on the ionization constants of acids and bases have been reviewed by Tremaine et al.² and Mesmer et al.²² The thermodynamic contributions to the ionization process can be described by the Born–Haber cycle shown in Fig. S5.† The major temperature effects are associated with the solvation processes − Δ_solvG°(HA), Δ_solvG°(A⁻) and Δ_solvG°(H₃O⁺). Under ambient conditions, the major solvation effects are due to short-range solute–solvent interactions associated with hydrogen bonding in the primary and secondary hydration spheres. At elevated temperatures, the hydrogen-bonded “structure” of water breaks down, and long-range ion–solvent polarization effects become important.^23,24 At temperatures above about 300 °C, the high compressibility of water associated with classical near-critical effects causes long-range polarization by ion–solvent interactions to dominate. At temperatures above 350 °C, hydrophilic and hydrophobic hydration of neutral species both have a very strong effect on the thermodynamics of solvation.^2,22,23 Indeed, the decrease in the dissociation constants of most acids, bases and ion-pairs with increasing temperature at steam saturation pressures is caused by the negative entropy associated with the orientation of water molecules due to long-range polarization effects.^2,23

The calculated pK_a temperature-dependence (using the four types of Δ_fG° values) and their errors in H₂O and D₂O are displayed in Tables S6 (Fig. S6) and S7 (Fig. S7), respectively, of the ESI† section. The best results are obtained at the B3LYP-IEF-PCM/6-311++G(d,p) level of theory, using G_w values (MAE = 2.02 (H₂O) and 3.59 (D₂O)). Given that very few experimental values (seven in H₂O and two in D₂O) are available to judge the accuracy of these calculations, and that the errors are considerably large, we conclude that the methods applied are unable to describe the temperature-dependence (under ambient and hydrothermal conditions) of pK_a values in light and heavy water.

The continuum solvation models used here simulate all the solute–solvent interactions through the polarization of the surrounding dielectric continuum by multipoles associated with the local functional groups. Key parameters are the shape and size of the cavity occupied by the molecule, solvent properties previously mentioned (dielectric constant, molar volume), and the computational details of the treatment for calculating the reaction field.^7g This treatment is more suitable for non-polar solvents and non-hydrogen-bonded polar solvents. Our previous calculations^8b–d have shown that accurate results for the ionization of organic acids in H₂O at room temperature can sometimes be obtained by these methods, in part because the solute–solvent interactions of large organic groups are similar for the acid and its conjugate base, and in part because the cavity parameters in the software used (Gaussian)¹⁴ have been optimized to yield the best possible results for water at 25 °C. For non-electrolytes and ions with large hydrophobic groups, the challenge is that no continuum model so far has been able to reproduce the large increase in the partial molar volumes of hydrophobic solutes under near-critical conditions.

3.2.2. The deuterium isotope effect on pK_a (ΔpK_a = pK_a(D₂O) − pK_a(H₂O)). The DIE on pK_a values, which is the difference in pK_a between the two solvents at a given set of temperature and pressure conditions, is the quantity best calculated with the methods applied in this study. DIE calculations at the B3LYP-IEF-PCM/6-311++G(d,p) level of theory, i.e., using G_w values, under ambient and hydrothermal conditions, are shown in Table 4 (see the first column of calculated values for each set of conditions) together with the experimental data available, which is much more abundant under ambient conditions, as previously indicated. These values are also plotted in Fig. S8.† The agreement between experiment and theory is excellent, with MAE values of 0.11 (ambient conditions) and 0.06 (hydrothermal conditions), and this is, without doubt, a consequence of the cancellation of systematic errors present in the pK_a calculations in both solvents at a given set of temperature and pressure conditions. Similar excellent results have been obtained in a follow-up work on this topic, currently in preparation, dealing with more than seventy organic and inorganic acids with a variety of structures. The methodology applied is a very important tool for accurately predicting DIE on pK_a values under both ambient and hydrothermal conditions. Given that very few experimental pK_a values are available in water (light or heavy) under hydrothermal conditions, this finding becomes particularly important.

Table 4 Experimental and calculated DIE (ΔpK_a = pK_a(D₂O) − pK_a(H₂O)) for the acids studied

Acid	25 °C, 101.3 kPa				250 °C, 20.0 MPa
	Exp. ΔpKa^a	Calc. ΔpKa^a^,^b using eqn (12) (ΔΔG, ΔTCG)	Calc. ΔpKa^c using eqn (19) (ΔTCE)	Calc. ΔpKa^d using eqn (20) (ΔZPE)	Exp. ΔpKa^a	Calc. ΔpKa^a^,^b using eqn (12) (ΔΔG, ΔTCG)	Calc. ΔpKa^c using eqn (19) (ΔTCE)	Calc. ΔpKa^d using eqn (20) (ΔZPE)
a Calculated using the values shown in Table 3.b (eqn (12), (14) and (16) are equivalent).c .d .e Mean absolute error.
Phenol	0.62	0.63	0.45	0.53		0.37	0.12	0.22
3-Methoxyphenol	0.58	0.63	0.45	0.54		0.37	0.12	0.22
4-Methoxyphenol	0.61	0.62	0.45	0.54		0.37	0.12	0.21
4-Bromophenol	0.59	0.64	0.45	0.53		0.38	0.12	0.22
2-Nitrophenol	0.55	0.68	0.54	0.62		0.38	0.16	0.27
4-Nitrophenol	0.53	0.63	0.45	0.54		0.36	0.12	0.22
2,4-Dinitrophenol	0.50	0.67	0.53	0.61		0.38	0.15	0.26
2,5-Dinitrophenol	0.50	0.69	0.53	0.61		0.38	0.15	0.26
2,6-Dinitrophenol	0.49	0.67	0.53	0.61		0.37	0.15	0.26
3,5-Dinitrophenol	0.57	0.62	0.44	0.52		0.36	0.12	0.21
4-Chloro-2,6-dinitrophenol	0.47	0.67	0.52	0.61		0.40	0.15	0.26
β-Naphthol	0.54	0.63	0.46	0.54	0.39	0.37	0.12	0.22
β-Naphthoic acid	0.42	0.61	0.45	0.54		0.36	0.11	0.22
Acetic acid	0.50	0.62	0.46	0.54	0.44	0.35	0.12	0.22
s-Collidine		0.71	0.60	0.66		0.38	0.18	0.29
Acridine		0.68	0.57	0.64		0.38	0.17	0.28
Average	0.53	0.65	0.49	0.57	0.42	0.37	0.14	0.24
Range	0.42–0.62	0.61–0.71	0.44–0.60	0.52–0.66	0.39–0.44	0.36–0.40	0.10–0.22	0.21–0.29
MAE^e		0.11	0.08	0.07		0.06	0.29	0.19

---

DIE on pK_a values cannot be reproduced when solvent effects are accounted for only during single-point energy calculations (see Tables S8 and S9†), or when using the Onsager solvation method. The best preliminary results (see Appendix A in the ESI† section) were obtained at the level of theory chosen for this study. Calculations of accurate DIE values require working with the lowest Gibbs free energy conformation of both the acid and its conjugate base. Otherwise, significant variations on a given pK_a value are observed and the calculated DIE values would have greater errors when compared with experiment.

Given the success of the methodology applied in accurately reproducing the DIE on pK_a values under ambient and hydrothermal conditions, it is possible to additionally predict pK_a values in heavy water using the corresponding light water experimental pK_a value (see Table 3) and the calculated DIE value at the same temperature and pressure (see Table 4). Seven pK_a values in D₂O have been predicted for 2-nitrophenol, 4-nitrophenol, β-naphthoic acid, s-collidine, and acridine. These values are shown in Table 3 within brackets. As a test of this approach, predicted and experimental values in D₂O are shown in Table S10.† The MAE of these predictions coincides with the MAE previously reported for the DIE values under both sets of conditions (0.11 and 0.06). These are excellent results with important implications.

Another important aspect to observe in both the experimental and calculated DIE values under ambient and hydrothermal conditions is the fact that, in spite of the significant structural differences between the compounds considered, these values are very similar for a given set of temperature and pressure conditions. Inspecting Fig. 1, which focuses on ambient conditions, would lead us to the same realization: there is an almost constant and similar difference between the experimental and theoretical plots of pK_a values in light and heavy water. The experimental and theoretical DIE (ΔpK_a) range of values are 0.42–0.62 (0.39–0.44) and 0.61–0.71 (0.36–0.40), respectively, under ambient (hydrothermal) conditions, and the average values are 0.53 (0.42) and 0.65 (0.37) for both experiment and theory, respectively, under ambient (hydrothermal) conditions. The well-known increase in pK_a values for a given acid when going from H₂O to D₂O becomes smaller under hydrothermal conditions, i.e., smaller DIE values are calculated at high temperatures and pressures. The methodology chosen is able to reproduce this experimental trend.

The remarkably constant calculated ΔpK_a values, in excellent agreement with experiment, linked to the facts that substituted phenols can exhibit a wide range of pK_a values at 25 °C depending on the nature of the substituent (see Table 3), and that many of them are thermally stable and have UV-visible spectra, lead to a potentially useful application for these compounds. If the calculated ΔpK_a values under hydrothermal conditions could be verified experimentally by studying two or three representative systems, the substituted phenols and naphthols may well form a practical class of thermally-stable pH indicators for studying deuterium isotope effects at elevated temperatures.

3.3. Investigating the DIE on pK_a values

Bunton and Shiner, in their classical paper of 1961, derived an empirical equation for estimating DIE on pK_a values (see eqn (11)).²⁵ Neglecting long-range polarization effects, isotopic entropy effects, tunneling, anharmonicities, librations and bending modes, the authors took into account the differences in solute–solvent hydrogen bonding for the species involved in the acid–base equilibrium. In the absence of experimental data, they estimated O–H vibrational frequencies of these hydrogen bonds (the v_H values that appear in eqn (11)), making use of pK_a and pK_b values of the acids and bases involved in Scheme (2a). Their DIE estimate for acetic acid at ambient conditions was very good, but their model fails to make accurate predictions under hydrothermal conditions.

To the best of our knowledge, this paper reports the first calculations of DIE on acid dissociation constants under ambient and hydrothermal conditions applying electronic structure calculations. Hence, we felt curious to explore these calculations and previous observations from a different angle. The ab initio calculations presented here include all the vibrational (using the harmonic oscillator model), rotational (using the free-rotor model) and translational (using the free-particle-in-a-box model) degrees of freedom of a molecule. The continuum solvation methods applied replace the explicit presence of solvent molecules interacting with the solute by building a solvent cavity around the solute. Hydrogen bonding between the solute and solvent molecules is not explicitly considered. Long-range polarization effects are treated by multipole electrostatic interactions with a cavity that, except for the Onsager spherical cavity, conforms to the molecule's shape. The calculations ignore anharmonicities and treat the solvent as an incompressible medium.

Combining eqn (6a) and (6b) to calculate ΔpK_a values, we obtain eqn (12). The second term of this equation, , is a constant (Q₁) equal to 3.2396 × 10⁻⁴, which is independent of temperature and pressure. If eqn (7a) and (7b) are taken into account eqn (12) becomes eqn (13).

The term: is a temperature-dependent constant (see Fig. B1 in Appendix B of the ESI† section) that will be labelled Q₂. At 25 °C, Q₂ = −0.002004 au and at 250 °C, Q₂ = −0.002340 au.

As can be seen from Table B1,† the difference is basically zero both under ambient and hydrothermal conditions. Hence, eqn (13) can be further simplified to eqn (14). The remaining difference, , is almost constant for each acid, depending on temperature and pressure (see Table B1†). Under ambient conditions it is calculated in the range 0.00333–0.00354 au (with an average value of 0.00341 au and a standard deviation of 0.00006 au), while under hydrothermal conditions it is calculated in the range 0.00368–0.00380 au (with an average value of 0.00374 au and a standard deviation of 0.00003 au).

As shown, the main contributor to the DIE on pK_a values is the difference in Δ_fG° in solution for each acid and its deuterated form in the corresponding solvent. The Δ_fG° values in solution are calculated by adding the uncorrected energy in solution (the value at the bottom of the potential energy well obtained after the geometry optimization in solution has taken place, E_H2O or E_D2O) to the corresponding thermal correction to the Gibbs free energy (TCG_H2O or TCG_D2O). Taking this into account, eqn (14) becomes eqn (15).

As shown in Table B2,† the difference E_H2O(HA) − E_D2O(DA) is basically zero for the two sets of conditions; hence, eqn (15) can be further simplified, as shown in eqn (16). Next, the expression to calculate the TCG values can be further investigated. TCG values in general are calculated using eqn (17), where TCE is the thermal correction to the energy at a given temperature, k_B is the Boltzmann constant, T is the absolute temperature and S is the entropy of the system at this temperature. The TCE contains the ZPE and additional thermal corrections (ATC, which accounts for the fact that at temperatures greater than 0 K, additional vibrational states beyond v = 0 become available to the system) at the temperature of interest. Taking eqn (17) into account, eqn (16) becomes eqn (18). An alternative but equivalent derivation of these equations could make use of harmonic frequencies and molecular partition functions.

TCG = TCE + k_BT − TS = ZPE + ATC + k_BT − TS (17)

The values of ZPE, TCE, TS and TCG for each acid in H₂O and D₂O under ambient and hydrothermal conditions are displayed in Tables B3 and B4,† respectively. The differences between these quantities for a given set of conditions are displayed in Table B5.† It is of interest to note that the calculated differences in TCG, ΔTCG = TCG_H2O(HA) − TCG_D2O(DA), are equal to the differences between G values, , (see Table B1†) for each acid under both sets of conditions. This fact quantitatively confirms the validity of the simplification of eqn (14)–(16), previously derived. In other words, eqn (12), (14) and (16) are equivalent.

It can be seen that the difference TΔS (TS_H2O(HA) − TS_D2O(DA)) is much smaller than that of ΔTCE(TCE_H2O(HA) − TCE_D2O(DA)), so we could explore the effect of additionally simplifying eqn (18) to (19). Furthermore, given that the ZPE is the main contributor to the TCE (see Tables B3–B5†), we could also investigate the validity of eqn (20), in which the term ΔATC − TΔS = [ATC_H2O(HA) − ATC_D2O(DA)] − [TS_H2O(HA) − TS_D2O(DA)] is neglected.

ΔpK_a values calculated using eqn (19) and (20) are shown in Table 4 to facilitate their comparison with the experimental values available and those calculated using eqn (12) (same as eqn (14) and (16)). In simplifying eqn (16) into eqn (19), the error introduced (−TΔS; average = 0.00034 au) is small and the MAE is slightly better (0.08) than calculated using eqn (12) (MAE = 0.11) under ambient conditions. However, even though the experimental data for comparison under hydrothermal conditions are few, the predicted values using eqn (19) seem to be underestimated (MAE = 0.29) compared to those obtained using eqn (12) (MAE = 0.06). This is in line with a much greater error made in the simplification (−TΔS; average = 0.00088 au).

In simplifying eqn (16) into eqn (20), the error introduced (ΔTCG − ΔZPE = ΔATC − TΔS; average = 0.00017 au) is even smaller than in the previous simplification and the MAE is also slightly better (0.07) than calculated using eqn (12) (MAE = 0.11) under ambient conditions. However, the predicted values under hydrothermal conditions, even though better than when using eqn (19), are still a bit lower (MAE = 0.19) than when using eqn (12) (MAE = 0.06). The error introduced (ΔTCG − ΔZPE = ΔATC − TΔS; average = 0.00049 au) is also larger than for ambient conditions but smaller than when using eqn (19). In all cases, the error made in simplifying eqn (12) into eqn (20) is smaller than when making the simplification into eqn (19).

The description above indicates that the difference in ZPE of the acid and its deuterated form in H₂O and D₂O seems to account quite well for the DIE on pK_a values under ambient conditions, which is in agreement with previous empirical work done on this topic.^25–27 However, this approximation would be insufficient under hydrothermal conditions where the effects of vibrational excitation and entropy changes (ΔATC − TΔS) seem to play a more important role. TΔS (a negative quantity) decreases when going from ambient (average: −0.00034 au) to hydrothermal (average: −0.00088 au) conditions (see Table B5†). ΔATC (also a negative quantity) gets reduced as well (ambient average: −0.00017 au; hydrothermal average: −0.00039 au), but to a lesser degree, which causes an overall increase in the error introduced (see above) when attempting to simplify eqn (16) into eqn (20). That is, the approximation ΔZPE ≈ ΔTCG works quite well under ambient conditions, but it does not work under hydrothermal ones. It should also be pointed out that within the framework of the harmonic approximation, ZPE values are calculated by means of eqn (21) using the 3N − 6 (or 3N − 5 for linear systems of N atoms) vibrational frequencies (ν) of a molecular system, which are not temperature-dependent. Hence, regardless of temperature and pressure, the calculated ΔZPE of any system will always be a fixed quantity (see Table B5†).

The observation that pK_a values in D₂O are greater than in H₂O for a given acid (i.e., acids in D₂O are weaker than in H₂O) has been partially explained by several authors in terms of zero-point energy (ZPE) differences in the O–H and O–D bonds.^1b,25–27 One would assume that such an explanation could also be extended to acids in which the acidic atom is attached to atoms other than oxygen (e.g., nitrogen, carbon, etc.). The simplest way to explain this fact resembles the way kinetic isotope effects are explained. Replacing the acidic H atom with D increases the reduced mass (μ) of the system. The reduced-mass increment is much greater when we focus on the bond between this acidic atom and the rest of the acid species. This leads to a reduction of the vibrational frequencies (, k is the force constant), particularly related to this bond. In turn, this causes an overall decrease of the ZPE, which leads to an increase in the energy required to break (dissociate) the acidic bond. Hence, the acid becomes weaker when dissociating in D₂O. The calculations reported in this paper clearly show (see Tables B3–B5†) that the ZPE of an acid in H₂O decreases when in D₂O. Other secondary factors involving solvation differences in H₂O and D₂O for the acid and the ions formed after dissociation (the conjugate base of the acid and either H⁺ or D⁺) could also be discussed, but the explanation provided agrees with the experimental fact (see Table 3) and the derivations previously made.

A mathematical explanation for why the DIE on pK_a values gets reduced as temperature increases from ambient to hydrothermal conditions (see Table 4) can be found by inspecting the previously derived equations. If we focus on the non-simplified eqn (16), it can be seen that while the numerator, ΔTCG + Q₂ = ΔZPE + ΔATC − TΔS + Q₂, slightly increases with temperature, it does so at a much lower rate than the denominator (RT ln 10). Other explanations are based on the entropic effects of long-range solvent polarization (i.e., solvent compressibility effects), as discussed by Mesmer et al.²²

From the work done so far it can be concluded that for any set of conditions, the DIE on pK_a values can be accurately calculated from differences between thermal corrections to the Gibbs free energy (or differences between the standard Gibbs free energies of formation) of the acid and its deuterated form in H₂O and D₂O. The set of quantum-mechanical calculations to be performed at a given temperature and pressure conditions does not need to involve the conjugate base of the acid under study. Accurate predictions under ambient conditions can also be made from ZPE differences between of the acid and its deuterated form in H₂O and D₂O. The application of continuum solvation methods (e.g., PCM and related ones) on both geometry optimizations and frequency calculations of the acids capture the subtle but almost-constant difference in pK_a values for a given acid in light and heavy water under ambient and hydrothermal conditions. The results obtained seem to indicate that the difference in the number and strength of solute–solvent hydrogen bonds is not a determining factor when quantifying DIE on pK_a values.

4. Conclusions

Quantum electronic structure methods are applied for the first time to the study of deuterium isotope effects on pK_a values under ambient (25 °C, 101.3 kPa) and hydrothermal (250 °C, 20.0 MPa) conditions. Several methodologies for calculating pK_a values and various pK_a differences in H₂O and D₂O are explored. After preliminary calculations, the B3LYP-PCM/6-311++G(d,p) level of theory was applied to the study of sixteen organic acids (acetic acid, β-naphthol, s-collidine, acridine, β-naphthoic acid and eleven phenols). None of the methods applied are adequate for directly predicting accurate pK_a values for these compounds.

When solvent effects are accounted for in geometry optimizations and frequency calculations, excellent agreement with experiment is obtained with the calculated DIE on pK_a values (ΔpK_a = pK_a(D₂O) − pK_a(H₂O)) under both ambient and hydrothermal conditions. Using the calculated DIE values and the experimental pK_a values in H₂O, excellent predictions of pK_a values in D₂O can be made for both sets of conditions studied. Following this approach, pK_a values in D₂O for 2-nitrophenol, 4-nitrophenol, β-naphthoic acid, s-collidine and acridine are predicted to be 7.23, 6.93, 6.34, 4.64 (8.14), 3.79 (6.26), respectively, at 250 °C and 20.0 MPa (25 °C and 101.3 kPa). The experimental and calculated DIE values are almost constant for a given set of temperature and pressure conditions. The average calculated ΔpK_a values under ambient (experimental average: 0.53) and hydrothermal conditions were 0.65 and 0.37, respectively. The mean absolute error between calculated and experimental ΔpK_a values under ambient conditions was 0.11.

It has been shown that continuum solvation methods, frequently used to account for solvent effects under ambient conditions in light water, can be successfully applied to predict DIE on pK_a values under ambient and hydrothermal conditions. Furthermore, it has been demonstrated that DIE on pK_a values are determined by the difference between the standard Gibbs free energies of formation (or the difference between the thermal corrections to the Gibbs free energy) of the acid and its deuterated analogue in each solvent. However, accurate predictions at room temperature can also be made from zero-point energy differences.

Acknowledgements

We gratefully acknowledge the Natural Sciences and Engineering Research Council of Canada (NSERC), the Basque Foundation for Science (Ikerbasque), Thompson Rivers University (CUEF U-REAP program) and the University Network of Excellence in Nuclear Engineering (UNENE) for financial support. The preliminary initial calculations reported in the Appendix of the ESI† were performed by Dr Elena Formoso. Earlier work by Tiffany S. Fransbergen, a former TRU undergraduate student, is also acknowledged.

References

1. (a) R. Bates, in Solute–Solvent Interactions, ed. J. F. Coetzee and C. D. Ritchie, Marcel Dekker, New York, 1969, ch. 4; (b) P. M. Laughton and R. E. Robertson, in Solute–Solvent Interactions, ed. J. F. Coetzee and C. D. Ritchie, Marcel Dekker, New York, 1969, ch. 7, and references there in.
2. P. R. Tremaine, K. Zhang, P. Benezeth and C. Xiao, in Aqueous Systems at Elevated Temperatures and Pressures: Physical Chemistry in Water, Steam and Aqueous Solutions, ed. D. A. Palmer, R. Fernandez-Prini and A. H. Harvey, Elsevier Academic Press, Amsterdam, 2004, ch. 13.
3. (a) M. H. Lietzke and R. W. Stoughton, J. Phys. Chem., 1963, 67, 652; (b) D. W. Shoesmith and L. Woon, Can. J. Chem., 1976, 54, 3553; (c) R. E. Mesmer and D. L. Herting, J. Solution Chem., 1978, 7, 901.
4. E. Bulemela and P. R. Tremaine, J. Solution Chem., 2009, 38, 805.
5. K. M. Erickson, H. Arcis, D. Raffa, G. H. Zimmerman and P. R. Tremaine, J. Phys. Chem. B, 2011, 115, 3038.
6. (a) J. Ehlerova, L. N. Trevani, J. Sedlbauer and P. Tremaine, J. Solution Chem., 2008, 37, 857; (b) B. Balodis, M. Madekufamba, L. N. Trevani and P. R. Tremaine, Geochim. Cosmochim. Acta, 2012, 93, 182.
7. (a) E. L. Wehry and L. B. Rogers, J. Am. Chem. Soc., 1966, 88, 351; (b) D. C. Martin and J. A. V. Butler, J. Chem. Soc., 1939, 1366; (c) L. Pentz and E. R. Thornton, J. Am. Chem. Soc., 1967, 89, 6931; (d) R. P. Bell and A. T. Kuhn, Trans. Faraday Soc., 1963, 59, 1789; (e) R. A. Robinson, J. Chem. Eng. Data, 1969, 14, 247; (f) R. A. Robinson, M. Paabo and R. G. Bates, J. Res. Natl. Bur. Stand., Sect. A, 1969, 73, 299; (g) W. T. Wofford, Spectroscopic studies of acid–base behavior and pH measurement in supercritical water, Ph.D. Dissertation, University of Texas, 1997.
8. See for example: (a) G. C. Shields and P. G. Seybold, QSAR in Environmental and Health Sciences: Computational Approaches for the Prediction of pK_a Values, CRC Press, London, GBR, 2013; (b) T. N. Brown and N. Mora-Diez, J. Phys. Chem. B, 2006, 110, 9270; (c) T. N. Brown and N. Mora-Diez, J. Phys. Chem. B, 2006, 110, 20546; (d) I. E. Charif, S. M. Mekelleche, D. Villemin and N. Mora-Diez, J. Mol. Struct., 2007, 818, 1.
9. (a) J. R. Pliego and J. M. Riveros, J. Phys. Chem. A, 2002, 106, 7434; (b) C. O. Silva, E. C. da Silva and M. A. C. Nascimento, J. Phys. Chem. A, 1999, 103, 11194; (c) C. O. Silva, E. C. da Silva and M. A. C. Nascimento, J. Phys. Chem. A, 2000, 104, 2402; (d) E. Soriano, S. Cerdán and P. Vallesteros, J. Mol. Struct., 2004, 684, 121.
10. (a) M. J. Citra, Chemosphere, 1999, 38, 191; (b) K. R. Adam, J. Phys. Chem. A, 2002, 5, 187; (c) K. Murłowska and N. Sadlej-Sosnowska, J. Phys. Chem. A, 2005, 109, 5590.
11. See for example: (a) J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999; (b) F. J. Luque, C. Curutchet, J. Muñoz-Muriedas, A. Bidon-Chanal, I. Soteras, A. Morreale, J. L. Gelpı and M. Orozco, Phys. Chem. Chem. Phys., 2003, 5, 3827; (c) J. Ho and M. L. Coote, Theor. Chem. Acc., 2010, 125, 3; (d) J. Ho, M. L. Coote, M. Franco-Perez and R. Gomez-Balderas, J. Phys. Chem. A, 2010, 114, 11992; (e) J. Ho and M. L. Coote, J. Chem. Theory Comput., 2009, 5, 295; (f) K. S. Alongi and G. C. Shields, in Annual Reports in Computational Chemistry, Elsevier B. V., 2010, vol. 6, ch. 8.
12. J. R. Pliego and J. M. Riveros, J. Phys. Chem. A, 2001, 105, 7241.
13. (a) I. Ivanov and M. L. Klein, J. Am. Chem. Soc., 2002, 124, 13380; (b) I. Ivanov, B. Chen, S. Raugei and M. L. Klein, J. Phys. Chem., 2006, 110, 6365.
14. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery Jr, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez, and J. A. Pople, Gaussian 03, revision C.02, Gaussian, Inc., Wallingford CT, 2004.
15. (a) M. T. Cancès, B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032; (b) M. Cossi, V. Barone, B. Mennucci and J. Tomasi, Chem. Phys. Lett., 1998, 286, 253; (c) B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 106, 5151.
16. (a) V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995; (b) M. Cossi, N. Rega, G. Scalmani and V. Barone, J. Comput. Chem., 2003, 24, 669.
17. (a) M. W. Wong, K. B. Wiberg and M. J. Frisch, J. Am. Chem. Soc., 1992, 114, 523; (b) J. G. Kirkwood, J. Chem. Phys., 1934, 2, 351; (c) L. Onsager, J. Am. Chem. Soc., 1936, 58, 1486.
18. (a) A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378; (b) C. J. Cramer and D. G. Truhlar, SMx Continuum Models for Condensed Phases, in Trends and Perspectives in Modern Computational Science; Lecture Series on Computer and Computational Sciences, ed. G. Maroulis and T. E. Simos, Brill/VSP, Leiden, 2006, pp. 112–140.
19. (a) W. Wagner and A. Pruss, J. Phys. Chem. Ref. Data, 2002, 31, 387; (b) D. P. Fernandez, A. R. H. Goodwin, E. W. Lemmon, J. M. H. Levelt-Sengers and R. C. Williams, J. Phys. Chem. Ref. Data, 1997, 26, 1125; (c) P. G. Hill, R. D. C. MacMillan and V. A. Lee, J. Phys. Chem. Ref. Data, 1982, 11, 1; (d) “ASME and IAPWS Formulation for Water and Steam”, NIST Standard Ref. Database 10, 2.2; (e) “REFPROP: Equations of State for Pure and Binary Fluids” NIST Standard Ref. Database 22, 8.0; (f) L. N. Trevani, E. Balodis and P. R. Tremaine, J. Phys. Chem. B, 2007, 111, 2015.
20. M. Palascak and G. C. Shields, J. Phys. Chem. A, 2004, 108, 3692.
21. (a) E. T. Ryan, T. Xiang, K. P. Johnston and M. A. Fox, J. Phys. Chem. A, 1997, 101, 1827; (b) T. Xiang and K. P. Johnston, J. Solution Chem., 1997, 26, 13.
22. R. E. Mesmer, W. L. Marshall, D. A. Palmer, J. M. Simonson and H. F. Holmes, J. Solution Chem., 1988, 17, 699.
23. T. M. Seward and T. Driesner, in Aqueous Systems at Elevated Temperatures and Pressures: Physical Chemistry in Water, Steam and Aqueous Solutions, ed. D. A. Palmer, R. Fernandez-Prini and A. H. Harvey, Elsevier Academic Press, Amsterdam, 2004, ch. 5.
24. R. Fernández Prini, M. L. Japas and E. J. Marceca, J. Supercrit. Fluids, 2010, 55, 472.
25. C. A. Bunton and V. J. Shiner, J. Am. Chem. Soc., 1961, 83, 42.
26. G. N. Lewis and P. W. Schutz, J. Am. Chem. Soc., 1934, 56, 1913.
27. O. Halpern, J. Chem. Phys., 1935, 3, 456.

---

Footnote

† Electronic supplementary information (ESI) available: Calculated raw data with the B3LYP functional and the IEF-PCM (UAHF) solvation method (Table S1). Calculated pK_a values and their errors in H₂O and D₂O at 25 °C/101.3 kPa and 250 °C/20 MPa using the B3LYP functional and the IEF-PCM (UAHF) solvation method (Tables S2–S5, Fig. S1–S4). A Born–Haber cycle for the ionization of an acid in light and heavy water (Fig. S5). Calculated temperature-dependence of pK_a and errors in H₂O and D₂O (pK_a(25 °C) − pK_a(250 °C)) using the B3LYP functional and the IEF-PCM solvation method (Tables S6 and S7, Fig. S6 and S7). Calculated DIE on pK_a values (ΔpK_a = pK_a(D₂O) − pK_a(H₂O)) and errors under ambient and hydrothermal conditions using the B3LYP functional and the IEF-PCM solvation method (Tables S8 and S9, Fig. S8). Predicted pK_a values in D₂O using experimental pK_a values in H₂O and calculated DIE values under ambient and hydrothermal conditions (Table S10). Preliminary initial calculations working with a subset of acids at different levels of theory (Appendix A, 42 pages). Data obtained while investigating the numeric source of the DIE on pK_a values (Appendix B, 7 pages). Cartesian coordinates of the acids and conjugate bases studied calculated at the B3LYP-IEF-PCM/6-311++G(d,p) level of theory. See DOI: 10.1039/c4ra14087g

---